DM1_Picaud_Sammier
Nicolas
2024-03-31
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, with a few tweaks. It is included to compare against the G/G/k implementation.
Question 1 : The G/G/K code
# This function plots the average response time of a M/G/k queue for various scheduling disciplines.
# As input, it takes the job service times (Service), their mean value (ES),
# a range of arrival rates for the Poisson process (lambdas), and the number of servers (k),
# as well as optional parameters to toggle logging and plotting.
plot_response_times_GGk <- function(Service, ES, lambdas, k, loglevel=1, plot=TRUE) {
N = length(Service); # number of jobs to simulate;
debug=0;
scheduling_disciplines = c(1:4); # 'FCFS'=1,'LCFS'=2,'RO-S'=3,'SRPT'=4
# 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-PR","LCFS-PR","RO-S-PR","SRPT-PR");
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
if (loglevel>=2){cat("Arrival=", Arrival, "\nService=", Service, "\n")}
for (scheduling_discipline in scheduling_disciplines)
{
if(ONLY_FCFS && scheduling_discipline!= 1) {break;}
if (loglevel>=2){cat("\tSimulating discipline ", rownames(R_time)[scheduling_discipline], "\n")}
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;
CurrentTasks = rep(k, x=NA);
NextArrival = 1;
while (TRUE) {
if (loglevel>=3){
print("*********************");
print(cat("Arrival = ", Arrival));
print(cat("Service = ", Service));
print(cat("Remaining = ", Remaining));
}
dtA = NA; # Time to next arrival
dtC = NA; # Time to next completion
if(length(Arrival[Arrival>t])>0) { # if an arrival exists after t
NextArrivalTime = head(Arrival[Arrival>t],n=1) ; # time to next arrival
dtA = NextArrivalTime - t;
}
if(length(CurrentTasks[!is.na(CurrentTasks)]) > 0) { # if a task is running
# time to next completion: min of all remaining times for non-NA current tasks
dtC = min(Remaining[CurrentTasks[!is.na(CurrentTasks)]]);
}
if(is.na(dtA) & is.na(dtC)) {
break; # No event, exit
}
dt = min(dtA,dtC,na.rm=T); # Time to next "event" (arrival or completion)
# Move to next event, update system variables
# t = t + dt; # WARNING: this is not the same as the line below
# Because we use a "==" condition to determine if the simulation time
# has reached the NextArrival time, we need to be extremely careful with
# floating-point precision errors. If we detect that dt is dtA, we avoid
# the addition altogether and set the time to exactly the NextArrival time.
t = ifelse(!is.na(dtA) & dt==dtA, NextArrivalTime, t+dt);
if (is.na(t)) browser()
AvgJobs=AvgJobs+dt*sum(!is.na(Remaining));
# If the event is a new task arrival
if((NextArrival <=N) & (Arrival[NextArrival] == t)) {
# This task's Remaining is its length, because it only just started
Remaining[NextArrival] = Service[NextArrival];
NextArrival = NextArrival + 1; # The next arrival is the next task
}
# Job departure
# for each server that has a running task:
for (i in which(!is.na(CurrentTasks))) {
task = CurrentTasks[i]
# Decrease its Remaining
Remaining[task] = Remaining[task] - dt;
if (Remaining[task] <= 0) {
# If its Remaining hits zero, save completion time and clear task
Completion[task] = t;
Remaining[task] = NA;
CurrentTasks[i] = NA;
}
}
if (loglevel>=2){
cat("Sim. time:", t, "| Jobs arrived:", NextArrival-1, "| Jobs completed: ", length(Completion[!is.na(Completion)]), "\n");
}
# Scheduling
# Our various scheduling disciplines each have 2 versions.
# - preemptive: running tasks are removed from the servers and treated
# just like new tasks, which means they can be stopped or moved
# - non-preemptive: running tasks are excluded from scheduling, they
# run on the same server until completion
# We implement the preemptive version here
CurrentTasks = rep(k, x=NA);
# for each server:
for (i in 1:k) {
# Therefore the WaitingList = all tasks up to NextArrival that exist in Remaining,
# excluding the ones that exist in CurrentTasks
WaitingList=(1:NextArrival)[!is.na(Remaining)];
WaitingList=WaitingList[!(WaitingList %in% CurrentTasks)];
if (length(WaitingList)>0) {
# FCFS
if (scheduling_discipline==1){
CurrentTasks[i] = head(WaitingList,n=1);
}
# LCFS
if (scheduling_discipline==2){
CurrentTasks[i] = tail(WaitingList,n=1);
}
# ROS: Random Order of Service
if (scheduling_discipline==3){
if (length(WaitingList)>1) {
CurrentTasks[i] = sample(WaitingList,size=1);
} else{
CurrentTasks[i] = WaitingList;
}
}
# SRPT: Shortest Remaining Processing Time
if (scheduling_discipline==4){
# Find minimum in Remaining, excluding already scheduled tasks
rem=Remaining;
for (task in CurrentTasks) { rem[task]=NA; }
CurrentTasks[i] = which.min(rem);
}
}
if (loglevel>=2){cat("Current Tasks = ", CurrentTasks, "\n")}
} # end of scheduling
} # end of main loop
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 (loglevel>=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
title = paste("G/G/k,", N, "jobs,", k, "servers")
if (loglevel>=1){
cat("\n", title, "\n")
cat("Response time matrix:\n")
print(R_time)
}
if(plot){
matplot(t(R_time), type = "l", xlab="Load (lambda*E[S])", ylab="Avg Response Time", main=title, 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));
}
return(R_time);
} # end functionQuestion 2 : Proving that the code is correct
We take \(k = 3\), the number of servers. \(N = 10\) the number of jobs simulated. \(\mu = 1\) and \(ES = 1/\mu\). The number of arrival rates simulated are \(5\)
Simulate the dynamics induced by 10 arrivals :
## Sim. FCFS-PR (rho= 0.17 ) completed. Little law verification: 0.1668 =N=lambda*R= 0.2391
## Sim. LCFS-NP (rho= 0.17 ) completed. Little law verification: 0.1668 =N=lambda*R= 0.2391
## Sim. RO-S-NP (rho= 0.17 ) completed. Little law verification: 0.1668 =N=lambda*R= 0.2391
## Sim. SRPT-PR (rho= 0.17 ) completed. Little law verification: 0.1668 =N=lambda*R= 0.2391
## Sim. FCFS-PR (rho= 0.33 ) completed. Little law verification: 1.5473 =N=lambda*R= 0.9704
## Sim. LCFS-NP (rho= 0.33 ) completed. Little law verification: 1.5754 =N=lambda*R= 0.988
## Sim. RO-S-NP (rho= 0.33 ) completed. Little law verification: 1.7435 =N=lambda*R= 1.0934
## Sim. SRPT-PR (rho= 0.33 ) completed. Little law verification: 1.2008 =N=lambda*R= 0.7531
## Sim. FCFS-PR (rho= 0.5 ) completed. Little law verification: 2.2738 =N=lambda*R= 1.7888
## Sim. LCFS-NP (rho= 0.5 ) completed. Little law verification: 1.9961 =N=lambda*R= 1.5704
## Sim. RO-S-NP (rho= 0.5 ) completed. Little law verification: 1.6393 =N=lambda*R= 1.2897
## Sim. SRPT-PR (rho= 0.5 ) completed. Little law verification: 1.4584 =N=lambda*R= 1.1474
## Sim. FCFS-PR (rho= 0.67 ) completed. Little law verification: 2.5598 =N=lambda*R= 2.4875
## Sim. LCFS-NP (rho= 0.67 ) completed. Little law verification: 3.3434 =N=lambda*R= 2.5775
## Sim. RO-S-NP (rho= 0.67 ) completed. Little law verification: 2.4744 =N=lambda*R= 2.4045
## Sim. SRPT-PR (rho= 0.67 ) completed. Little law verification: 2.2413 =N=lambda*R= 1.6195
## Sim. FCFS-PR (rho= 0.83 ) completed. Little law verification: 2.09 =N=lambda*R= 2.9305
## Sim. LCFS-NP (rho= 0.83 ) completed. Little law verification: 1.6836 =N=lambda*R= 2.2235
## Sim. RO-S-NP (rho= 0.83 ) completed. Little law verification: 1.5535 =N=lambda*R= 2.1782
## Sim. SRPT-PR (rho= 0.83 ) completed. Little law verification: 1.6339 =N=lambda*R= 1.8656
##
## G/G/1, 10 jobs
## Response time matrix:
## rho=0.167 rho=0.333 rho=0.5 rho=0.667 rho=0.833
## FCFS-PR 1.434482 2.911081 3.577672 3.731247 3.516596
## LCFS-NP 1.434482 2.963926 3.140731 3.866218 2.668167
## RO-S-NP 1.434482 3.280152 2.579390 3.606817 2.613832
## SRPT-PR 1.434482 2.259158 2.294784 2.429253 2.238756
## Sim. FCFS-PR (rho= 0.17 ) completed. Little law verification: 0.1668 =N=lambda*R= 0.2391
## Sim. LCFS-PR (rho= 0.17 ) completed. Little law verification: 0.1668 =N=lambda*R= 0.2391
## Sim. RO-S-PR (rho= 0.17 ) completed. Little law verification: 0.1668 =N=lambda*R= 0.2391
## Sim. SRPT-PR (rho= 0.17 ) completed. Little law verification: 0.1668 =N=lambda*R= 0.2391
## Sim. FCFS-PR (rho= 0.33 ) completed. Little law verification: 0.7997 =N=lambda*R= 0.4782
## Sim. LCFS-PR (rho= 0.33 ) completed. Little law verification: 0.7997 =N=lambda*R= 0.4782
## Sim. RO-S-PR (rho= 0.33 ) completed. Little law verification: 0.7997 =N=lambda*R= 0.4782
## Sim. SRPT-PR (rho= 0.33 ) completed. Little law verification: 0.7997 =N=lambda*R= 0.4782
## Sim. FCFS-PR (rho= 0.5 ) completed. Little law verification: 1.0786 =N=lambda*R= 0.7172
## Sim. LCFS-PR (rho= 0.5 ) completed. Little law verification: 1.0786 =N=lambda*R= 0.7172
## Sim. RO-S-PR (rho= 0.5 ) completed. Little law verification: 1.0786 =N=lambda*R= 0.7172
## Sim. SRPT-PR (rho= 0.5 ) completed. Little law verification: 1.0786 =N=lambda*R= 0.7172
## Sim. FCFS-PR (rho= 0.67 ) completed. Little law verification: 1.1514 =N=lambda*R= 0.9563
## Sim. LCFS-PR (rho= 0.67 ) completed. Little law verification: 1.1514 =N=lambda*R= 0.9563
## Sim. RO-S-PR (rho= 0.67 ) completed. Little law verification: 1.1514 =N=lambda*R= 0.9563
## Sim. SRPT-PR (rho= 0.67 ) completed. Little law verification: 1.1514 =N=lambda*R= 0.9563
## Sim. FCFS-PR (rho= 0.83 ) completed. Little law verification: 1.6653 =N=lambda*R= 1.2312
## Sim. LCFS-PR (rho= 0.83 ) completed. Little law verification: 1.7076 =N=lambda*R= 1.2312
## Sim. RO-S-PR (rho= 0.83 ) completed. Little law verification: 1.7318 =N=lambda*R= 1.2487
## Sim. SRPT-PR (rho= 0.83 ) completed. Little law verification: 1.7076 =N=lambda*R= 1.2312
##
## G/G/k, 10 jobs, 3 servers
## Response time matrix:
## rho=0.167 rho=0.333 rho=0.5 rho=0.667 rho=0.833
## FCFS-PR 1.434482 1.434482 1.434482 1.434482 1.477436
## LCFS-PR 1.434482 1.434482 1.434482 1.434482 1.477436
## RO-S-PR 1.434482 1.434482 1.434482 1.434482 1.498382
## SRPT-PR 1.434482 1.434482 1.434482 1.434482 1.477436
## rho=0.167 rho=0.333 rho=0.5 rho=0.667 rho=0.833## GG1 FCFS-PR 1.4344816 2.9110814 3.5776721 3.7312468 3.5165955
## GGk FCFS-PR 1.4344816 1.4344816 1.4344816 1.4344816 1.4774364
## Ratio GGk/GG1 1.0000000 0.4927659 0.4009539 0.3844510 0.4201326
##
## GG1 LCFS-NP 1.4344816 2.9639255 3.1407308 3.8662179 2.6681668
## GGk LCFS-PR 1.4344816 1.4344816 1.4344816 1.4344816 1.4774364
## Ratio GGk/GG1 1.0000000 0.4839803 0.4567350 0.3710297 0.5537272
##
## GG1 RO-S-NP 1.4344816 3.2801524 2.5793904 3.6068167 2.6138321
## GGk RO-S-PR 1.4344816 1.4344816 1.4344816 1.4344816 1.4983817
## Ratio GGk/GG1 1.0000000 0.4373216 0.5561320 0.3977140 0.5732509
##
## GG1 SRPT-PR 1.4344816 2.2591580 2.2947841 2.4292533 2.2387560
## GGk SRPT-PR 1.4344816 1.4344816 1.4344816 1.4344816 1.4774364
## Ratio GGk/GG1 1.0000000 0.6349629 0.6251052 0.5905031 0.6599363Little’s law problem
As we can see with the plots, the results are not consistent and the Little’s law is not well calibrated, this is due to the low amount of jobs provided to the system. With \(N=1000\) we can verify that the code is indeed correct.
## Sim. FCFS-PR (rho= 0.17 ) completed. Little law verification: 0.1918 =N=lambda*R= 0.207
## Sim. LCFS-NP (rho= 0.17 ) completed. Little law verification: 0.1919 =N=lambda*R= 0.2072
## Sim. RO-S-NP (rho= 0.17 ) completed. Little law verification: 0.1918 =N=lambda*R= 0.2071
## Sim. SRPT-PR (rho= 0.17 ) completed. Little law verification: 0.1704 =N=lambda*R= 0.1839
## Sim. FCFS-PR (rho= 0.33 ) completed. Little law verification: 0.5491 =N=lambda*R= 0.5473
## Sim. LCFS-NP (rho= 0.33 ) completed. Little law verification: 0.533 =N=lambda*R= 0.5312
## Sim. RO-S-NP (rho= 0.33 ) completed. Little law verification: 0.5443 =N=lambda*R= 0.5425
## Sim. SRPT-PR (rho= 0.33 ) completed. Little law verification: 0.446 =N=lambda*R= 0.4446
## Sim. FCFS-PR (rho= 0.5 ) completed. Little law verification: 1.3485 =N=lambda*R= 1.3864
## Sim. LCFS-NP (rho= 0.5 ) completed. Little law verification: 1.4171 =N=lambda*R= 1.457
## Sim. RO-S-NP (rho= 0.5 ) completed. Little law verification: 1.3201 =N=lambda*R= 1.3572
## Sim. SRPT-PR (rho= 0.5 ) completed. Little law verification: 0.8411 =N=lambda*R= 0.8647
## Sim. FCFS-PR (rho= 0.67 ) completed. Little law verification: 2.1744 =N=lambda*R= 2.2312
## Sim. LCFS-NP (rho= 0.67 ) completed. Little law verification: 2.1818 =N=lambda*R= 2.2388
## Sim. RO-S-NP (rho= 0.67 ) completed. Little law verification: 2.1905 =N=lambda*R= 2.2478
## Sim. SRPT-PR (rho= 0.67 ) completed. Little law verification: 1.2651 =N=lambda*R= 1.2981
## Sim. FCFS-PR (rho= 0.83 ) completed. Little law verification: 13.6781 =N=lambda*R= 13.0284
## Sim. LCFS-NP (rho= 0.83 ) completed. Little law verification: 13.7431 =N=lambda*R= 13.0787
## Sim. RO-S-NP (rho= 0.83 ) completed. Little law verification: 13.0694 =N=lambda*R= 12.4393
## Sim. SRPT-PR (rho= 0.83 ) completed. Little law verification: 4.4586 =N=lambda*R= 4.2415
##
## G/G/1, 1000 jobs
## Response time matrix:
## rho=0.167 rho=0.333 rho=0.5 rho=0.667 rho=0.833
## FCFS-PR 1.242160 1.641848 2.772808 3.346814 15.634131
## LCFS-NP 1.243027 1.593652 2.913908 3.358245 15.694438
## RO-S-NP 1.242306 1.627561 2.714479 3.371626 14.927181
## SRPT-PR 1.103618 1.333719 1.729468 1.947204 5.089855## Sim. FCFS-PR (rho= 0.17 ) completed. Little law verification: 0.1667 =N=lambda*R= 0.1799
## Sim. LCFS-PR (rho= 0.17 ) completed. Little law verification: 0.1667 =N=lambda*R= 0.1799
## Sim. RO-S-PR (rho= 0.17 ) completed. Little law verification: 0.1667 =N=lambda*R= 0.1799
## Sim. SRPT-PR (rho= 0.17 ) completed. Little law verification: 0.1667 =N=lambda*R= 0.1799
## Sim. FCFS-PR (rho= 0.33 ) completed. Little law verification: 0.3655 =N=lambda*R= 0.3611
## Sim. LCFS-PR (rho= 0.33 ) completed. Little law verification: 0.3649 =N=lambda*R= 0.3604
## Sim. RO-S-PR (rho= 0.33 ) completed. Little law verification: 0.3647 =N=lambda*R= 0.3603
## Sim. SRPT-PR (rho= 0.33 ) completed. Little law verification: 0.3646 =N=lambda*R= 0.3602
## Sim. FCFS-PR (rho= 0.5 ) completed. Little law verification: 0.5319 =N=lambda*R= 0.5465
## Sim. LCFS-PR (rho= 0.5 ) completed. Little law verification: 0.5313 =N=lambda*R= 0.546
## Sim. RO-S-PR (rho= 0.5 ) completed. Little law verification: 0.5319 =N=lambda*R= 0.5466
## Sim. SRPT-PR (rho= 0.5 ) completed. Little law verification: 0.5295 =N=lambda*R= 0.5441
## Sim. FCFS-PR (rho= 0.67 ) completed. Little law verification: 0.7529 =N=lambda*R= 0.7311
## Sim. LCFS-PR (rho= 0.67 ) completed. Little law verification: 0.7505 =N=lambda*R= 0.7287
## Sim. RO-S-PR (rho= 0.67 ) completed. Little law verification: 0.7485 =N=lambda*R= 0.7268
## Sim. SRPT-PR (rho= 0.67 ) completed. Little law verification: 0.7467 =N=lambda*R= 0.7251
## Sim. FCFS-PR (rho= 0.83 ) completed. Little law verification: 0.9327 =N=lambda*R= 0.923
## Sim. LCFS-PR (rho= 0.83 ) completed. Little law verification: 0.9349 =N=lambda*R= 0.9252
## Sim. RO-S-PR (rho= 0.83 ) completed. Little law verification: 0.9318 =N=lambda*R= 0.9221
## Sim. SRPT-PR (rho= 0.83 ) completed. Little law verification: 0.9231 =N=lambda*R= 0.9135
##
## G/G/k, 1000 jobs, 3 servers
## Response time matrix:
## rho=0.167 rho=0.333 rho=0.5 rho=0.667 rho=0.833
## FCFS-PR 1.079463 1.083207 1.093083 1.096673 1.107595
## LCFS-PR 1.079463 1.081345 1.092029 1.093104 1.110234
## RO-S-PR 1.079463 1.080848 1.093158 1.090180 1.106474
## SRPT-PR 1.079463 1.080460 1.088168 1.087632 1.096231## rho=0.167 rho=0.333 rho=0.5 rho=0.667 rho=0.833## GG1 FCFS-PR 1.2421603 1.6418475 2.7728080 3.3468142 15.6341313
## GGk FCFS-PR 1.0794629 1.0832069 1.0930828 1.0966731 1.1075955
## Ratio GGk/GG1 0.8690206 0.6597487 0.3942151 0.3276767 0.0708447
##
## GG1 LCFS-NP 1.2430271 1.5936524 2.9139081 3.3582448 15.6944382
## GGk LCFS-PR 1.0794629 1.0813451 1.0920289 1.0931036 1.1102345
## Ratio GGk/GG1 0.8684146 0.6785326 0.3747644 0.3254985 0.0707406
##
## GG1 RO-S-NP 1.2423062 1.6275614 2.7144791 3.3716260 14.9271812
## GGk RO-S-PR 1.0794629 1.0808482 1.0931580 1.0901798 1.1064735
## Ratio GGk/GG1 0.8689185 0.6640906 0.4027137 0.3233395 0.0741247
##
## GG1 SRPT-PR 1.1036184 1.3337190 1.7294677 1.9472043 5.0898548
## GGk SRPT-PR 1.0794629 1.0804596 1.0881678 1.0876324 1.0962309
## Ratio GGk/GG1 0.9781124 0.8101104 0.6291923 0.5585610 0.2153757We want to test the theory and the simulation, to do this we make these asumptions
- Arrivals occur at rate \(\lambda\) according to a Poisson process and move the process from state \(i\) to \(i+1\).
- Service times have an exponential distribution with parameter \(\mu\). If there are fewer than \(k\) jobs, some of the servers will be idle. If there are more than \(k\) jobs, the jobs will queue up in a buffer.
- The buffer is of infinite size, so there is no limit on the number of jobs it can hold
If the traffic is greater than 1 then the queue will grow without bound but if server utilization \(\rho = \frac{\lambda} {k * \mu} < 1\) then the system has a stationary distribution with a probability mass function.
For our theoretical calculations we will use the Erlang’s C formula, namely an implementation found on Github
\[ {\displaystyle {\text{ C}}(c,\lambda /\mu )={\frac {1}{1+\left(1-\rho \right)\left({\frac {c!}{(c\rho )^{c}}}\right)\sum _{k=0}^{c-1}{\frac {(c\rho )^{k}}{k!}}}}} \] The response time of our system can now be calculated with the Response time formula
\[ {\displaystyle {\frac {{\text{ C}}(c,\lambda /\mu )}{c\mu -\lambda }}+{\frac {1}{\mu }}.} \]
# This function calculates the Erlang B formula, which is used to find the probability of all servers being busy in a given time.
erlang_b <- function(number_of_servers, intensity) {
inverse_b <- 1.0
# Loop through each server
for(i in 1:number_of_servers) {
inverse_b <- 1.0 + inverse_b * i / intensity
}
# Return the inverse of the final calculation, which gives us the Erlang B formula
return(1.0 / inverse_b)
}
# This function calculates the Erlang C formula, which is used to find the number of calls that are blocked or delayed.
erlang_c <- function(number_of_agents, intensity) {
# Calculate Erlang C using the Erlang B formula
ec <- number_of_agents * erlang_b(number_of_agents, intensity) / (number_of_agents - intensity * (1 - erlang_b(number_of_agents, intensity)))
return(ec)
}
# Define the number of servers
K_values = seq(1, 3, by = 1);
# Define the number of jobs
N = 10;
# Set ONLY_FCFS to TRUE to only use First Come First Serve discipline
ONLY_FCFS = TRUE;
# Define the number of arrival rates
num_arrival_rates=1;
# Initialize a 2D matrix to store f_rtm_ggk, Wr, and error ratio
results_matrix <- matrix(0, num_arrival_rates * length(K_values), 3);
# Initialize index
index <- 0;
# Loop through each number of servers
for(j in seq_along(K_values)) {
# Set the number of servers
K=K_values[j];
# Set the service rate
mu = K;
# Calculate the expected service time
ES = 1/mu;
# Initialize matrices to store f_rtm_ggk and Wr
f_rtm_ggk <- matrix(0, num_arrival_rates, length(K_values));
Wr_matrix <- matrix(0, num_arrival_rates, length(K_values))
# Calculate the arrival rates
lambdas_a = ((1/ES) * c(1:num_arrival_rates)/(num_arrival_rates+1));
# Set the seed for reproducibility
set.seed(global_seed);
# Generate service times
Service = rexp(N,rate=mu);
# Calculate the response times for the G/G/k queue
set.seed(global_seed);
rtm_ggk = plot_response_times_GGk(Service, ES, lambdas_a, K, loglevel=0, plot=FALSE);
f_rtm_ggk[,j] <- rtm_ggk[1, , drop = TRUE];
# Loop through each arrival rate
for(i in 1:num_arrival_rates){
# Set the arrival rate
lambda = lambdas_a[i];
# Calculate the traffic intensity
p <- lambda / (K * mu);
# Check if system is stable
if (p >= 1) {
cat("System is not stable for arrival rate", lambda, "\n")
Wr = Inf
} else {
# Calculate the probability of having K customers in the system using Erlang-C formula
C_K = erlang_c(K, p);
# Calculate average number of customers in the queue using Little's law
Lq = p * C_K / (1 - p)
# Calculate average time spent in the queue using Little's law
Wq = Lq / lambda
# Calculate average service time
Ws = 1 / mu
# Calculate average response time as sum of average time spent in the queue and average service time
Wr = Wq + Ws
Wr_matrix[i, j] <- Wr
}
# Increment index
index <- index + 1;
# Store the results
results_matrix[index, 1] = f_rtm_ggk[i,j];
results_matrix[index, 2] = Wr;
results_matrix[index, 3] = abs(f_rtm_ggk[i,j] - Wr) /Wr;
}
}
# Create a vector of indices for the x-axis
indices <- 1:nrow(results_matrix)
# Create a matrix with f_rtm_ggk and Wr columns
plot_matrix <- cbind(results_matrix[, 1], results_matrix[, 2])
# Create a vector of colors for each row
colors <- c("red", "blue")
# Plot the results matrix
title = paste("Response time comparison for variating \n number of servers and ", N, "jobs")
matplot(indices, plot_matrix, type = "l", lty = 1, col = colors, xlab = "Number of servers", ylab = "Response Times", main = title)
# Add a legend to the plot
legend("topright", legend = c("Sim", "Formulas"), col = colors, lty = 1, title = "Response Time Values")
For 10 jobs and K ∈ {1,2,3} we can see that the theory and the real simulation are still a bit far, but with more jobs, we directly see that our calculations are correct.
Comparing the G/G/k queue with a G/G/1 queue having one server that is k times faster.
For increasing values of the load (i.e., arrival rate divided by service rate) and k = 5, we plot the ratio between the average response time obtained by a G/G/k/FCFS queue and the average response time obtained by a G/G/1/FCFS queue with a server that operates at speed k.
N = 1e3;
K = 5;
mu = 1;
ES_GG1 = 1/mu;
ES_GGK = 1/mu;
num_arrival_rates=5;
set.seed(global_seed);
ServiceGGK = rexp(N,rate=mu);
ServiceGG1 = ServiceGGK/K # All jobs are just K times shorter
Trials = 6;
Factor = 2;
# Initialize 3D arrays
f_rtm_gg1 <- matrix(0, Trials, num_arrival_rates);
f_rtm_ggk <- matrix(0, Trials, num_arrival_rates);
ONLY_FCFS = TRUE;
for(i in 1:Trials){
# On simule pour plusieurs espérances différentes de taux d'arrivée
lambdas_gg1 = ((1/ES_GG1) * c(1:num_arrival_rates)/(num_arrival_rates+1)) * (i*Factor);
lambdas_ggk = ((1/ES_GGK) * c(1:num_arrival_rates)/(num_arrival_rates+1)) * (i*Factor);
set.seed(global_seed)
rtm_gg1 <- plot_response_times_GG1(ServiceGG1, ES_GG1, lambdas_gg1, loglevel=0, plot=FALSE)
set.seed(global_seed)
rtm_ggk <- plot_response_times_GGk(ServiceGGK, ES_GGK, lambdas_ggk, K, loglevel=0, plot=FALSE)
# Store the first row of rtm_gg1 and rtm_ggk in the ith layer of the 3D array
f_rtm_gg1[i,] <- rtm_gg1[1, , drop = FALSE]
f_rtm_ggk[i,] <- rtm_ggk[1, , drop = FALSE]
# Create a matrix with f_rtm_gg1 and f_rtm_ggk stacked vertically for this trial
f_rtm <- rbind(f_rtm_gg1[i,], f_rtm_ggk[i,])
# Create a vector of colors for this trial
colors <- c("red", "blue")
# Create a vector of line types for this trial
lty <- c(1, 2)
# Plot f_rtm using matplot for this trial
matplot(t(f_rtm), type = "l", lty = lty, col = colors, xlab = "Arrival rate", ylab = "Response time", xaxt = "n")
# Add a legend
legend("topright", c("GG1", "GGk"), col = c("red", "blue"), lty = c(1, 2))
axis(side=1,at=1:ncol(f_rtm),labels=round(lambdas_gg1*ES_GG1, digits=2));
# Add a title to the plot with the trial number
title(main = paste("Trial", i))
}
The result show us that the 2 policies are coming closer together as the load increases. This shows that even though the G/G/1 queue is running K times faster than the G/G/K one, when the load increases there is a limit and response times for both are converging to \(1\).
The GG1 model is a single-server queueing model, while the GGK model is a multi-server queueing model with k servers. The main difference between the two models is the number of servers. In the GG1 model, there is only one server, so customers must wait in line if the server is busy. In the GGK model, there are k servers, so customers can be served by any available server. This means that the GGK model can handle a higher load of customers than the GG1 model, as there are more servers available to serve customers.
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.2670567 0.3114101 0.3467650 0.3715889 0.4461399
## [2,] 0.2876802 0.3905408 0.4867681 0.6145179 0.8075711
## [3,] 0.3115170 0.5279784 0.6924491 0.9799191 0.9921905
## [4,] 0.3355186 0.6606590 0.9833863 0.9929205 0.9979455
## [5,] 0.3601067 0.8220261 0.9951935 0.9988562 0.9988027
## [6,] 0.3836634 0.9809875 0.9953320 0.9993364 0.9990856
We can see thanks to the graph that when the Arrival rate augments, the difference between the G/G/1 queue and the G/G/C queue is lowering untill reaching a ratio of \(1\).
This part seems a bit odd at first but the “overall performance” of each system are the same. Although we will find the system empty more often in the G/G/1 queue but the wait time and expected queue length are less for the G/G/K.
This can be shown with the following formulas. Take \(G/G/1\) and \(G/G/2\) with \(\mu_1 = 2\mu_2\)
\[ \rho_2 = \frac{\lambda}{2\mu_2} = \rho_1 = \frac {\lambda}{\mu_1} \] Now for the G/G/1 queue we find \[ L_q = \frac {\rho²}{1-\rho} ; W_q = \frac {\rho²}{\lambda(1-\rho)} ; \pi_0 = 1 - \rho \] G/G/2 \[ L_q = \frac {\rho³}{1-\rho²} ; W_q = \frac {\rho³}{\lambda(1-\rho²)} ; \pi_0 = \frac {1 - \rho}{1+\rho} \] For a understanding of the utilization, we have twice as many servers, each working half as fast as the alone server, so the total amount being done is the same.
For a general M/M/n queue the utiliztion is \[ \rho_k=\frac{\lambda}{k\mu} \]
Finding the value of the ratio whilethe load approaches the maximum.
We plotted the theoretical values for the G/G/1 and G/G/K queue and spotted the moment when the G/G/1 server response time is higher than the G/G/K queue.
As we can see on the plot, we can find the point where the two curves touch each other
# Define parameters
K <- 5
mu <- K
mu_1 <- K * mu
# Define functions for response times
Wr_5 <- function(lambda) {
p <- lambda / (K * mu)
if (p >= 1) return(Inf)
C_K <- erlang_c(K, p)
Lq <- p * C_K / (1 - p)
Wq <- Lq / lambda
Ws <- 1 / mu
Wr <- Wq + Ws
return(Wr)
}
Wr_1 <- function(lambda) {
p_1 <- lambda / mu_1
RespTime <- p_1 / (mu_1 - lambda)
return(RespTime)
}
# Define range of lambda values
# The range should be up to the capacity of the system, not beyond it
lambda_range <- seq(0, min(K * mu, mu_1), length.out = 1000)
# Calculate response times for each lambda value
Wr_5_values <- sapply(lambda_range, Wr_5)
Wr_1_values <- sapply(lambda_range, Wr_1)
# Remove infinite values and keep only the common indices
valid_indices <- !is.infinite(Wr_5_values) & !is.infinite(Wr_1_values)
Wr_5_values <- Wr_5_values[valid_indices]
Wr_1_values <- Wr_1_values[valid_indices]
lambda_range <- lambda_range[valid_indices]
filter=1;
# Filter response times higher than 5
Wr_5_values_filtered <- Wr_5_values[Wr_5_values < filter]
Wr_1_values_filtered <- Wr_1_values[Wr_1_values < filter]
# Get corresponding lambda values
lambda_range_filtered_5 <- lambda_range[Wr_5_values < filter]
lambda_range_filtered_1 <- lambda_range[Wr_1_values < filter]
# Plot filtered response times
plot(lambda_range_filtered_1, Wr_1_values_filtered, type = "l", xlab = "Arrival rate (lambda)", ylab = "Response time", col = "blue")
lines(lambda_range_filtered_5, Wr_5_values_filtered, type = "l", col = "red")
legend("topright", c("K = 5", "K = 1 (faster server)", "GGK sim", "GG1 sim"), col = c("red", "blue", "purple", "green"), lty = 1)
K=5;
set.seed(global_seed);
ServiceGGK = rexp(N,rate=mu);
ServiceGG1 = ServiceGGK/K;
ES = 1/mu;
f_rtm_ggk <- matrix(0, 1, 1);
f_rtm_gg1 <- matrix(0, 1, 1);
set.seed(global_seed);
rtm_ggk = plot_response_times_GGk(ServiceGGK, ES, 21, K, loglevel=0, plot=FALSE);
f_rtm_ggk[1,] <- rtm_ggk[1, , drop = FALSE]
set.seed(global_seed);
rtm_gg1 = plot_response_times_GG1(ServiceGG1, ES, 24, loglevel=0, plot=FALSE);
f_rtm_gg1[1,] <- rtm_gg1[1, , drop = FALSE]
abline(h = f_rtm_ggk[1,], col="purple", ylab="GGK");
abline(h = f_rtm_gg1[1,], col="green", ylab="GG1");
abline(v = 21, col="black");
axis(side=1,at=21, labels="λ*");
Sources : - Github Erlang’s formula - M/M/K queues formulas