排隊論的matlab仿真(包括仿真代碼)

1、Wireless NetworkExperiment Three:Queuing TheoryABSTRACTTins expenment is designed to leani the fiuidamentals of the queuing theory. Mainly about the M/NI/S and M/M/ii/n queuing models.KEY WORDS: queuing theory,M/M/n/n, Erlang B, Eilang C.INTRODUCTIONA queue is a waiting line and queueing thcoiy is t

2、lie mathematical tlicory of waiting lines. More generally, queueing tlicory is conccmed witli tlic niatlicinatical modeUng and analysis of systems that provide service to random demands, hi conumuiication networks, queues aic encountered everywhere. For example, the incoming data packets arc randoml

3、y arrived and buffered, waiting for the router to deliver. Such situation is considered as a queue. A queueing model is an absti act desciiption of such a system. Typically, a queueing model represents (1) the system's physical configuration by specifying the munber and arrangement of the server

4、s, and (2) tlie stochastic nature of the demands, by specifying the vaiiability in the amval process and in the service processTlie essence of queueing thcoiy is tliat it takes into accoiuit the raiidonuiess of tlie anival process and the randoiimcss of the service process The most conunon assmnptio

5、n about tlie anival process is tliat the customer arrivals follow a Poisson process, whci c the times between arrivals are exponentially distributed. Tlie probability of tlie exponential distribution functionis f(t) = Ae-Xt. Erlsuig B modelOne of the most impoitant queueing models is tlic Erlang B m

6、odel (i.c., M/M/11/n). It assumes that the airivals follow a Poisson process and have a finite n servers hi Erlaiig B model, it assumes that the anival customers arc blocked and clcaied when all tlie scivers arc busy. Tlic blocked probability of a Erlang B model is given by the famous Erlang B fonmi

7、la,AZ(1.1)"(%)=畀空,where n is the inunber of scivers and A=A/p is the offered load in Erlangs, A is the anival rate and 1/p is the average service time. Fonniila (1.1) is hard to calculate directly from its right side when n and A arc large. However, it is easy to calculate it using the followin

8、g iterative scheme:(m = 1,2, ，；P(0. A) = 1) Erlang C modelThe Erlang delay model (M/M/11) is similar to Erlang B model, except tliat now it assumes tliat tlic arrival customers aic waiting in a queue for a sciver to become available witliout coiisideling tlie lengtli of tlie queue. Tlic probability

9、ofblocking (all tlic seivers arc busy) is given by the Erlang C Formula,a(1-3)Pc(n,A)=n!(p)寸”一 1屮心厶人二o TT 十 n!(i-p)Wlicrc p = 1 if A > n and p = if A < n. The quantity p indicates the server utilization.(14)The Erlang C fomnila (1.3) cau be easily calculated by tlic following itciativc schemen

10、 A(1 PD(n, A)wlicrc Pq (11, A) is defined in Eq.(l.l).DESCRIPTION OF THE EXPERIMENTS1. Using the formulsi (1.2), calculate the blocking probability of the Erlang B model. Draw the relationsliip of the blocking probability PB(nA) and offered traffic A with n = 1厶 1°9 20, 30, 40, 50, 60, 70, 80,

11、90, 100. Compare it with tlie table in the text book (P.281, table 103).From tlie introduction, we know tliat when the n and A aic large、it is easy to calculate the blocking probability using the fonmila 1.2 as follows.PB(6 A)=APb(h-IA)m+APB(n-l,A)it use the tlicoiy of rcciusion for tlie calculation

12、. But the denominator and the inunerator of the fonmila both need to recins( PB (n 1,A) when doing the matlab calculation, it waste time and reduce tlie matlab calculation efficient. So wc change tlie fonmila to be :Pb®A) = "APb；" = "(I + aPb(；7A)APb(h-U)Tlicn the calculation onl

13、y need recurs once time and is more efficientHie liiatlab code for tlie for mil a is:erlsuigban% File: erlanb b.m% A = offered traffic in Erlangs% n = number of truncked channels% Pb is the result blocking probabilityfunction Pb = erlang_b( Az n ) if n=0Pb=l;% P(0, A)=lelsePb=l/ (1+n/(A*erlang_b(A,

14、n-1);% use recursion nerlang (Ar n-1)endendAs wc can see from the table on the text books, it uses die logaiitlun coordinate, so wc also use tlic logarithm coordinate to plot tlic result. Wc divide tlic inunber of scivers(11) into three parts, for each part wc can define a intcival of the traffic in

15、tcnsity(A) based on the figure on the text books :1. when 0<ii<10. 0.1<A<10.2. when 10<n<20,3<A<20.3. when 30<n<100,13<A<120.For each part, use tlic "erlang b” function to calculate and then use "loglog" function to figure tlic logantlun coordinateT

16、he niatlab code is :% for the three parts % n is the number servers% A is the traffic indensity% P is the blocking pr*obability.n_l = 1:2;A_1 = linspace (0. lr 10r 50) ; % 50 points bmtv/een 0l and 10n_2 = 10:10:20;A_2 = linspace(3/ 20f 50);n_3 = 30:10:100;A_3 = linspace(13/120r 50);% for each part,

17、 call the erlang_b () function.for i = 1:length(n_l)for j = 1: length (A_l)P_1 (jr i) = erlang_b (A_l (j), n_l (i);end endfor 丄=1: length (n_2)for j = 1: length (A_2)P_- (jr i) = erlang_b (A_2 (j), n_2 (i); endendfor i = 1: length (n_3)for j = 1: length (A_3)p_3 (j, i) = erlang_b (A_3 (j), n_3 (i);

18、endend% use loglog to figure the result within logantlun coordinate.loglog(A_lf p_lr fk-f, A_2, p_2r1k-1, A_3rp_3r );xlabel(Traffic indensity in Erlangs (A) 1) ylabel(f Probability of Blocking (P)T)axis(0.1 120 0.001 0.1)text(. 115, . 115r n=：L')text( 6, 115f 1n=2f)text(7, .115,110f)text(17,text

19、 (二7,text(45,.115, f20f).115f f30f).115f f50f)text(100/ 115, 1100f)J10rr2103)204060 70 03 90 1CDio1io1imc r cfeftily n briar 影血vr10£rJ£*BA 二 _qeTdTlic figure oil the text books is as follow:01oxno0.001Qof Trunked Chanatle (C)G051.0 100 Ioteoaicy h Brkop100.0Figure 36 The Erlang Bcmk sbOMvi

20、ngthc probability of blocking ac functions of th。ntmbor of chanrwls and traffc intercity " Erlangs.Hie madab code for tlic formula is :erlang canHie madab code for tlic formula is :erlang canWc can sec fiom tlic two pictures tliat, tlicy arc exactly tlic same witli each otlicr except tliat the

21、result of the experiincnt have not considered tlic situation with n=3,4,5,.,12,14,16,18 2. Using the formula (1.4), calculate the blocking probability of the Erlang C modeL Draw the relatioiisliip of tlie blockiiig probability PC(nA) and offered traiTic A with n = 1,2,10,20,30,40,50,60, 70,80,90,100

22、From tlic introduction, wc know tliat tlic foniuila 1.4 is :Pq (n. A)=nPB(n>A)n-A(l-PB(n,A)Hie madab code for tlic formula is :erlang canSince each time wc calculate tlie PB (n. A), wc need to recurs n times, so the fonmila is not efficient. Wc change it to be:Pc(幾 A)=nPB(nA)nA(lPB(n,A)=1/h_A(1_P

23、b (mA)nPg(n,A)"/InAnPB(nA)Hie madab code for tlic formula is :erlang canTlicn we only need recurs once. Pg(11, A) is calculated by the "crlang b" fiuictioii as step 1.Hie madab code for tlic formula is :erlang can% File: erianb_b.m % A = offered traffic in Erlangs% n = number of trunc

24、ked channels.% Pb is the resuJLt blocking probability% erlang_b (Az n) is the function of step 1.function Pc = erlang_c( Az n )Pc=l/(A/n) + (n-A)/(n*erlang_b(Ar n);endTlicn to figure out the tabic in the logarithm coordinate as what shown in the step 1.Hie matlab code is :% for the three parts.% n i

25、s the number servers% A is the traffic indensity.% P_c is the blocking probability of erlangC modeJL50 points between 0l and 10n_l = 1:2；A_1 = linspace(0lr10r SO);n_2 = 10:10:20;A_2 = linspace(3/20f 50);n_3 = 30:10:100;A_3 = linspace(13/120r 50);% for each partr call the erlang_c() function.for i =

26、1:length(n_l)for j = 1: length (A_l)P_1_C (j, i) = erlang_c (A_l (j) , n_l (i); %p-OA°Eyeriang_cendendfor i = 1: length (n_2)for j = 1: length (A_2)p_2_c (j, i) = erlang_c (A_2 (j) , n_2 (i);endendfor i = 1: length (n_3)for j = 1: length (A_3)p_3_c (jr i) = erlang_c (A_3 (j) r n_3 (i);endend% u

27、se loglog to figure the resuJLt within logantlun coordinate.*費仃六loglog (A_l, p_l_cz f g* - f f A_2, p_2_c； f g*-f / A_3, p_3_cr 1 g*-1);xlabel(fTraffic indensity in Erlangs (A)1) ylabel(Probability of Blocking (P)1) axis(0.1 120 0.001 0.1)text(.115, .115,fn=lf)text ( 6rtext(G, 115r f n=2 f).115, 110

28、f)text(14/text (20f.115, f20f).115, f30f)text (30ftext (39/.115, f40f).115f f50f)text(47,115f f60f)text (55/.115,f70f)text(65rtext (75/.115,f80f).115,190f)text(85/.115f f100f)Tlic result of blocking probability table of erlang C model.n=2102D 刃4053 60 TO &O 100d 二匚 WEB 2£虧36Th»lk »

29、;de -irty tiEdet 盧)lerTlicn wc put the table of erlang B and erlang C in tlic one figure, to compare their characteristic io1w1Irrfic imdensil/ in Erbns p*)10°Tlic line witli4 * is the erlang C model, the line witliout4 is tlie erlaiig B model. Wc can see from the pictiue that, for a constant t

30、raffic intensity (A), the erlaiig C model has a liiglicr blocking probability than erlang B model. Tlic blocking probability is increasing with traffic intensity. The system perfonns better when has a larger n.ADDITIONAL BONUSWrite a program to simulate a M/M/k queue system with input parameters of

31、lainda, mg k.hi tins part, wc will firstly simulate the M/M/k queue system use matlab to get the figure of the perfonnance of the system such as tlie leave time of each customer and tlie queue length of tlic system.About tlic simulation, wc fiistly calculate tlic anivc time and the leave time for ea

32、ch customer Tlicn analysis out the queue length and the wait time for each customer use “for" loops.Tlicn we let the input to be lainda = 3, nm = 1 and S = 3、and analysis perfonnance of the system for the first 10 customers in detail.Finally, wc will do two test to compared the perfonnance of t

33、lic system witli input lamda = 1, inn = 1 and S = 3 and the input lamda = 4, nm = 1 and S = 3The niatlab code is:mms fiinction9mfiuictionblock_iatcjise_rate=MNIS_fimction(mcaii_aiT4ncaii_scn;pco_inuii,scivcr_inmi) % %first part: compute the arriving time interval, service time%intcival,waiting time,

34、 leaving time during tlic whole service intcival% %statc=zeros(5 ,peo_mun);%rcprcscnt the state of each customer by%using a 5*pco_mun matrix%tlic meaning of each line is: arriving time interval, service time%intcival, waiting time, queue length wiicn NO.ncustomcr%anivc, leaving timestatc(l.: )=cxpm

35、d( 1 /me an_ air, 1 .peoniun);%aniving time intcival between each%customcr follows expouctial distiibutionstatc(2,:)=expmd(l/mcan _scrv;l,pco_inun);%scnice time of each customer follows exp one ti al distributionfor i=l:sciver_nmnstatc(3 J: servcr_umn)=O;endaiT_time=ciunsmn(state(l、:)；%acciuniilate

36、arriving time inteival to compute%aniving time of each customerstatc(l ,:)=an_timc;state(5,1: sciver_inun)=siun(state(: ,1: servxrinun);%computc living time of first NO.scrvxr niun%customer by using fbimilar aniviiig time + service time scrv_dcsk=state(5,1:server num);%crcatc a vector to store leavi

37、ng time of customers wliich is in servicefor i=(s crvxr_niini+1): p c oinunif an time (i)>uiin(seiv desk)statc(3,i)=0;elsestatc(3 ,i)=uiin(sciv_dcsk) an_timc(i);%wiicn customer NO.i anives and die%servcr is all busy, the waiting time can be compute by%niimis amving time from the niiniinum leaving

38、 timeendstatc(5 j)=siun(statc(: j);for j=l:scivcr_inunif scrvr_clcsk(j)=niiii(sciv_dcsk)s erv _dc sk(j)=st at c (5 4)； breakend%rcplacc tlie nuiumiun leaving time by the first waiting customer's leaving timeendend%sccond part: compute tlie queue length during the whole service intcival% %zero_ti

39、me=0;%zero_timc is used to identify wiiich scivcr is emptyscrv_dcsk(l:seivcr_niun)=zcro time;block_nimi=0;block_line=0;for i=l:peo_inunif block_line=0find_max=0;for j=l: servernumif seiv_dcsk(j)=zero_tiinc find_max=l; %mcans tlicrc is empty sciver breakelse continueendendif find_max=l%updatc serv de

40、skseiv_desk(j)=state(5 j);for k=l: servernumif scrv dcsk(k)<an time(i) %bcforc the NO.i customer actually arrives tlicrc maybe some customer leaves crv_ d c sk(k)=zcr otimc;else continueendendelseif an_timc(i)>niin(scnrclcsk)%if a customer will leave before the NO.i%customcr anivefor k=l:seive

41、r_inunif i«T_timc(i)>scrv_dcsk(k) scrvr_dcsk(k)=state(5 a)； breakelse continueendendfor k=l:seiver_niunif an time(i)>sciv desk(k) scnr_desk(k)=zcro_tiinc;else continueendendelse%if no customer leave before die NO.i customer anivebl o ck_inuw=bl o ck_mmi+1; block-liiie=block_line+1;endende

42、lse %tlic situation tliat the queue length is not zeron=0;%computc the inunber of leaing customer before the NO.i customer arrives for k=l:scrvrcr_niunif an time(i)>sciv desk(k)11=11+;seiv_dcsk(k)=zcro_tiinc;else continueendendfor k=l:block_lincif an_ti me (i)> st ate (5 4 -k)n=n+l;else contin

43、ueendendif n<block_liiie+l% n<block_linc+l means the queue lcngtli is still not zerobio ck_iuun=blo ck_nuin+1;for k=0:n-lif statc(5 j -bl o ck_line+k)> aiitiinc(i)for m=l:seiver_niunif scnr_dcsk(m)=zcro_timcscrv_dcsk(m)=statc(5 j-blockJinc+k)breakelse continue endendelsecontinueendendbl o c

44、k_liiic=bl o cklinc -11+1;else %ii>=block_liiie+1 means the queue lciigtli is zero%update seiv desk and queue lengtlifor k=0:block_lineif an_timc(i)<statc(5 j-k)for m=l:sciver_niunif sciv_dcsk(m)=zcro_tinicserv_desk(m)=statc(5 j k)breakelse continueendendelsecontinueendendblock_linc=0;endendst

45、atc(4j)=block_linc;endplot(statc(l,:),*-)；figureplot(statc(2,:),g);figureplot(statc(3,:)尸 *);figureplot(statc(4,:),gfigureplot(statc(5,:),*-)；Since the system is M/M/S wliich means the aniving rate and service rate follows Poisson distribution wliilc tlie number of scivcr is S and die buffer length

46、is infinite、we can computeall the arriving time, service time、waiting time and leaving time of each customer.We can test tlic code with uiput lainda = 3, mu = 1 and S = 3 Figures ai c below.cutonrri nurrr er°oID2D60 customer ruroe7090txArriving time of each customer86532Service time of each cus

47、tomerWaiting time of each customercustomer numtet10.16 5 4 3 tauJfalndnb2l oQueue length when each customer ai riveLeaving time of eadi customer1010As lamda = mu*scrvrcr_niun, the load of tlic system could be very liigh. Then we will zoom iii die result pictuies to analysis the performance of the sy

48、stem for the firstly 10 customer.45 -35The first customer enter the system at about Is.3 5 2 5 2 1 £?£<' J 刖 OEL nii-ibcr10Arriving time of first 10 customer土 er 3希The queue length is 1 for the 7th customer.910ctslomer nurrbe-Queue length of first 10 customerLeaving time of first 10

49、 customer1. As wc have 3 seivcr in tliis test, tlie first 3 customer will be served witliout any delay.2. Tlie aniving time of customer 4 is about 1.4 and the niiiiinnun leaving time of customer in service is about 1.2. So customer 4 will be served immediately and the queue length is still 0.3. Cust

50、omer 1,4,3 is in service4. Tlic amving time of custouier 5 is about 1.8 and tlic niiiiiiinuii leaving time of customer in service is about 1.6. So customer 5 will be served inunediately and tlic queue length is still 0.5. Customer 1, 5 is in service6. Tlic amving time of customer 6 is about 2.1 and

51、tlierc is a empty seiver. So customer 6 will be saved inunediately and the queue lcngtli is still 07. Customer 1, 5,6 is in service8. Tlic amving time of customer 7 is about 2.2 and tlie niiiiiimun leaving time of customer in sendee is about 2.5. So customer 7 caiuiot be served immediately and tlic

52、queue length will be 1.9 Customer 1, 5,6 is in service and customer 7 is waiting10. Tlie arriving time of customer 8 is about 2.4 and tlic niiiiiiniiin leaving tiine of customer in sendee is about 2.5. So customer 8 caiuiot be served ininiediatcly and tlic queue length will be 2.11. Customer 1, 5,6 is in