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');