This code simulates the behavior of a M/M/1 queue with the FIFO discipline. However, mutatis mutandis, it can be easily changed to simulate the dynamics of a G/G/1 queue with an arbitrary scheduling discipline.

knitr::opts_chunk$set(echo = TRUE)

set.seed(17); # Set seed for reproducibility
N = 1e2; # number of jobs to simulate;

lambda = 0.5; # arrival rate
mu = 1; # service rate

Arrival = cumsum(rexp(n=N, rate = lambda)); # Arrival times
Service = rexp(N,rate=1/mu); # Service times
Remaining = rep(N, x=NA);    # Remaining service times of each job
Completion = rep(N, x=NA);   # Completion time of each job

t = 0; # simulation time

CurrentTask = NA;
NextArrival = 1;
while (TRUE) {
    # print(t);
    # 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;
    if((NextArrival <=N) & (Arrival[NextArrival] == t)) {
        Remaining[NextArrival] = Service[NextArrival];
        NextArrival = NextArrival + 1;
    }
    if(!is.na(CurrentTask)) {
        Remaining[CurrentTask] = Remaining[CurrentTask] - dt ;
        if(Remaining[CurrentTask] <= 0) {
            Completion[CurrentTask] = t;
            Remaining[CurrentTask] = NA;
        }
        CurrentTask = NA;
    }
    # FIFO scheduling discipline. Change here to implement other scheduling disciplines
    WaitingList=(1:NextArrival)[!is.na(Remaining)]; # Optimize this line!
    if(length(WaitingList)>0) {
        CurrentTask = head(WaitingList,n=1);
    }
}

ResponseTime = mean(Completion-Arrival); #average response time

print(ResponseTime)
## [1] 1.656654
# Make a plot of the ResponseTime as a function of lambda.