0

I wrote a Monte Carlo simulation for M/M/1 queuing system since it's deterministic. Now I tried modifying the same code for a G/G/1 queueing system with interarrival times and service times having a triangular distribution. Instead of drawing samples from Poisson and exponential distribution, I am drawing samples from the triangular distribution (for G/G/1). However, I am getting results that are hard to make sense of. The analytical approximation results say that the mean waiting time in the queue is 0.9708 which is kinda off from the result (~0.75) I am getting with the simulation.

clc
clear
tic
Nc = 1008;          
Tsys = zeros(1,Nc);    
Tia = zeros(1,Nc);     
Ta = zeros(1,Nc);      
Tsrv = zeros(1,Nc);    
Tnsrv = zeros(1,Nc);   
Txsys = zeros(1,Nc);    
WTqmc = zeros(1,Nc);   
MTqmc = zeros(1,Nc);   
MTsys = zeros(1,Nc);   
Tri1_min = 2; 
Tri1_mid = 3; 
Tri1_max = 5; 
Tri2_min = 1.8;
Tri2_mid = 2.8; 
Tri2_max = 4.5; 
pdat  = makedist('Triangular','a',Tri1_min,'b',Tri1_mid,'c',Tri1_max);
pdser = makedist('Triangular','a',Tri2_min,'b',Tri2_mid,'c',Tri2_max); 
Tia(1) = random(pdat);       
Ta(1) = Tia(1);              
Tsys(1) = Tsrv(1);           
Tnsrv(1) = Ta(1);            
Tsrv(1) = random(pdser);     
Txsys(1) = Ta(1) + Tsrv(1);  
WTqmc(1) = 0;                
MTqmc(1) = WTqmc(1);         
for i = 2:1:Nc

    Tia(i) = random(pdat);   
    Ta(i) = Ta(i-1) + Tia(i);
    Tsys(i) = Ta(i);
    if Ta(i) < Txsys(i-1)       
        Tnsrv(i) = Tnsrv(i-1) + Tsrv(i-1);
    else                        
        Tnsrv(i) = Ta(i);
    end
    Tsrv(i) = random(pdser);            
    Txsys(i) = Tnsrv(i) + Tsrv(i);      
    WTqmc(i) = Tnsrv(i)-Ta(i);          
    Tsys(i) = Txsys(i) - Ta(i);         
    MTqmc(i) = ((i-1)* MTqmc(i-1) + WTqmc(i))/i;    
    % MTsys(i) = ((i-1)* MTsys(i-1) + Tsys(i))/i;   
end
Tend = Txsys(Nc);
Nts = floor(Tend);          
dt = 1;
Time = zeros(1,Nts);        
Lq = zeros(1,Nts);
Time(1) = 0;
Lq(1) = 0;                  
for j = 2:1:Nts
    Time(j) = Time(j-1)+dt;
    Lq(j) = 0;
    for i = 1:1:Nc
        if (Ta(i) < Time(j) && Txsys(i) > Time(j) && Tnsrv(i) > Time(j))
            Lq(j) = Lq(j)+1;
        end
    end
end
MLqmc = mean(Lq)
MTqmc1 = MTqmc(Nc)
MTqmc_mean = mean(WTqmc)
SDmtamc = std(WTqmc)
SEmtamc = SDmtamc/sqrt(Nc)
toc

mean1= (2+3+5)/3;
mean2=(1.8+2.8+4.5)/3;
var1=(4+9+25-2*3-2*5-3*5)/18;
var2=(1.8*1.8+2.8*2.8+4.5*4.5-1.8*2.8-1.8*4.5-2.8*4.5)/18;
rho= (1/mean1)/(1/mean2);
ca2= var1/mean1^2;
cs2=var2/mean2^2;
Lq=((rho^2)*(1+cs2)*(ca2+(rho^2)*cs2))/(2*(1-rho)*(1+(rho^2)*cs2));
wq=Lq/(1/mean1)

%Plot Outcomes
figure(1)
plot(Ta,MTqmc)
xlabel('Arrivals')
ylabel('Mean Time in Queue')
title('Queue Length vs Time')
sevebebe
  • 11
  • 2
  • Yes, _but_: you lose the memoryless property of the exponential distribution, so it is important you pair up the arrival and service times properly for the system to work well. But undoubtedly your code already has this set-up well. In general, yes, you can change the distributions easily, but various theoretical bounds might no longer apply/be tight. – Nelewout May 12 '20 at 20:08
  • What language is that? It looks Matlab-y to me, but I'm not too familiar with Matlab to be sure. The idea looks solid, and nothing changes to your discrete-time recursions other than the distribution of the inter-arrival and service times. You should be OK. – Nelewout May 12 '20 at 20:18
  • Yea that is MATLAB. Thank you. I will report back after the changes – sevebebe May 12 '20 at 20:19
  • @N.Wouda I made some changes to my M/M/2 queue to convert it to G/G/1 queue but I am getting weird results. Could you please help me out – sevebebe May 14 '20 at 00:32
  • @N.Wouda The analytical approximation results say that the mean waiting time in the queue is 0.9708 which is nowhere the results I am getting with the simulation. – sevebebe May 14 '20 at 01:41
  • I do not immediately recognise your approximation, where did you find that? For the G/G/c queue, you can use Sakasegawa's formula. E.g. something like `((C_a^2 + C_s^2) / 2) * ((rho^(sqrt(2(c + 1)) - 1) / c(1 - rho)) * E[S]`. (I hope I got all the brackets right!). There is no closed-form solution for the G/G/c queue, however, so this too is but an approximation. – Nelewout May 14 '20 at 08:31
  • @N.Wouda I used the formula found here https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=7&cad=rja&uact=8&ved=2ahUKEwj0-ZeTurbpAhVumXIEHXToBtgQFjAGegQIBRAB&url=http%3A%2F%2Fweb.mst.edu%2F~gosavia%2Fqueuing_formulas.pdf&usg=AOvVaw1zi-EUzpqx4sbCJ5a96KRA – sevebebe May 15 '20 at 18:04
  • I think that might just be a poor approximation formula. You could use the general Sakasegawa formula, but again, these are approximations - they might not be all that close to what you observe empirically. – Nelewout May 16 '20 at 17:04

0 Answers0