M/M/1 Simulation with Matlab

A simple simulation of M/M/1 queue with Matlab. The distribution of the number of “packets” in the system is computed and compared with the theoretical result.

delta=0.1; % simulation step in sec
lambda=0.1; % arrival rate in packets per second
mu=0.2; % departure rate in packets per second
rho=lambda/mu;
M=50*3600/delta; % number of simulation step for 50 hours
a=zeros(M,1); % arrival log
d=zeros(M,1); % departure log
n=zeros(M,1); % number of packets at each simulation time step
n_cur=0; % current number of packets in the system
for i=1:M
 n_cur; % # packets before the process
 n(i)=n_cur; % record the # of packet before the process
 if rand < lambda*delta % probability of an arrival occurred = lambda * delta
  a(i)=1;
  n_cur=n_cur+1; 
 end
 if rand < mu*delta && n_cur > 0 % probability of departure occurred = mu * delta provided that there is a packet in the system
  d(i)=1;
  n_cur=n_cur-1;
 end
 if mod(i,1000)==0 % show progress for every 1000 simulation steps
  i 
 end
end
figure;
edges=0:10;
count=histc(n,edges);
bar(0:length(count)-1,count/sum(count));
hold; 
nlist=0:length(count)-1;
plot(nlist, (1-rho)*rho.^(nlist),'r-+');
legend('experiment','theory');
title('Distribution of number of packets in the system');

Now, let say if we want to find the delay distribution for the packet arrived when there are 2 packets in the system.

figure;
ainds=find(a==1); % the indices when we have an arrival
atimes=ainds*delta; % arrival times in seconds
dtimes=find(d==1)*delta; % departure times in second
delays=dtimes-atimes(1:length(dtimes)); % list of delays in second 
inds=find(n(ainds)==2); % find indices of packets arrived when # packets in system is 2
h=histc(delays(inds),0:60);
bar(0:60,h/sum(h));
hold;
ezplot(mu^3/2*t^2*exp(-mu*t),[0,60]);
legend('experiment','theory');
title('Delay distribution for packets arrived when # packets in system is 2');
xlabel('Time in second');
ylabel('Probability');