第3讲 离散事件系统仿真原理及程序

1、第3讲 离散事件系统仿真原理及程序主要内容（以排队系统为例）u3.1问题的提出u3.2 排队系统的组成u3.3 事件调度法u3.4系统的统计性能指标u3.5 实体到达模式与泊松过程3.1 问题的提出n 排队排队常常是件很令人恼火的事情尤其是在我们这样的人口大国。电话亭1978年在北京15%的电话要在1小时后才能接通。在电报大楼打电话的人还要带着午饭去排队。银行窗口，ATM医院火车售票交通理发游乐场的游乐项目3.1 问题的提出3.1 问题的提出n怎样缩短排队的等待时间？银行的排队叫号机。只是有序的组织了顾客，并没有减少等待时间。如果能实现知道轮到自己需要等待多少时间，再选择合适的时间来，岂不很好

2、？3.1 问题的提出例，单人理发馆系统, 设上午9:00开门, 下午5:00关门。特点：顾客到达时间一般是随机的, 为每个顾客服务的时间长度也是随机的。系统的状态：服务台的状态(忙或闲)、顾客排队等待的队长也是随机的。 状态量的变化只能在离散的随机时间点上发生。3.2排队系统的组成n排队系统的三个基本组成部分到达模式到达模式指临时实体临时实体按怎样的规律到达。服务模式服务模式指同一时刻有多少服务台服务台可以接纳临时实体，它们的服务要多少时间。排队规则排队规则服务台完成当前服务后，从队列中选择下一个实体服务的原则。n排队系统的基本结构到达模式排队服务机构到达离去3.2排队系统的组成nM/M/1的

3、含义到达时间间隔服从指数分布, 用M表示服务时间的分布服从指数分布, 用M表示单队单服务台3.2排队系统的组成n稳态平均延误时间n实体通过系统的稳态平均滞留时间wn稳态平均队长Qn系统中稳态平均实体数LL(t)=Q(t)+S(t)nDdniin/lim1nSDnWwniiinniin/ )(lim/lim11TdttQQTT/)(lim0LL t dtTQ tS t dtTTTTTlim( )/lim( ( )( )/00External definitions for single-server#include #include #include lcgrand.h /*Header fil

4、e for random-number generator.*/#define Q_LIMIT 100 /* Limit on queue length. */#define BUSY 1 /* Mnemonics for servers being busy */#define IDLE 0 /* and idle. */External definitions for single-server (continued)int next_event_type, num_custs_delayed, num_delays_required, num_events,num_in_q, serve

5、r_status;float area_num_in_q, area_server_status, mean_interarrival, mean_service,sim_time, time_arrivalQ_LIMIT + 1, time_last_event, time_next_event3, total_of_delays;FILE *infile, *outfile;void initialize(void);void timing(void);void arrive(void);void depart(void);void report(void);void update_tim

6、e_avg_stats(void);float expon(float mean);3.3事件调度法n事件（Event）引起系统状态发生变化的行为，系统是由事件来驱动的。 “顾客到达顾客到达”为一类事件顾客到达引起系统状态：服务员的“状态”从闲变到忙（如果无人排队)；或者另一系统状态排队的人数发生变化（队列人数增加)。 “顾客离去顾客离去”为一类事件服务完毕后离开系统服务台“状态”由忙变成闲。3.3事件调度法n定义系统事件类型: 类型1 顾客到达事件 类型2 顾客接受服务事件（本程序未考虑）类型3 顾客服务完毕并离去事件n定义程序事件:仿真运行到150个时间单位(例如分钟)结束。n假定已经得到

7、到达时间间隔随机变量的样本值为: n系统初始状态: 思考思考：请列出该排队：请列出该排队系统的事件表。系统的事件表。,22,40,24,32,1554321AAAAASSSS123443363428,qZ0000, 3.3事件调度法n事件表main() /* Main function. */ /* Open input and output files. */ infile = fopen(mm1.in, r); outfile = fopen(mm1.out, w); /* Specify the number of events for the timing function. */ n

8、um_events = 2; /* Read input parameters. */ fscanf(infile, %f %f %d, &mean_interarrival, &mean_service, &num_delays_required); /* Write report heading and input parameters. */ fprintf(outfile, Single-server queueing systemnn); fprintf(outfile, Mean interarrival time%11.3f minutesnn, mean

9、_interarrival); fprintf(outfile, Mean service time%16.3f minutesnn, mean_service); fprintf(outfile, Number of customers%14dnn, num_delays_required);Main function/* Initialize the simulation. */ initialize();/* Run the simulation while more delays are still needed. */ while (num_custs_delayed num_del

10、ays_required) /* Determine the next event. */ timing(); /* Update time-average statistical accumulators. */ update_time_avg_stats(); /* Invoke the appropriate event function. */ switch (next_event_type) case 1: arrive(); break; case 2: depart(); break; Main function (continued)Main function (continu

11、ed) /* Invoke the report generator and end the simulation. */ report(); fclose(infile); fclose(outfile); return 0;void initialize(void) /* Initialization function. */ /* Initialize the simulation clock. */ sim_time = 0.0; /* Initialize the state variables. */ server_status = IDLE; num_in_q = 0; time

12、_last_event = 0.0; /* Initialize the statistical counters. */ num_custs_delayed = 0; total_of_delays = 0.0; area_num_in_q = 0.0; area_server_status = 0.0; /* Initialize event list. Since no customers are present, the departure (service completion) event is eliminated from consideration. */ time_next

13、_event1 = sim_time + expon(mean_interarrival); time_next_event2 = 1.0e+30;Function Initializevoid timing(void) /* Timing function. */ int i; float min_time_next_event = 1.0e+29; next_event_type = 0; /* Determine the event type of the next event to occur. */ for (i = 1; i = num_events; +i) if (time_n

14、ext_eventi Q_LIMIT) /* The queue has overflowed, so stop the simulation. */ fprintf(outfile, nOverflow of the array time_arrival at); fprintf(outfile, time %f, sim_time); exit(2); /* There is still room in the queue, so store the time of arrival of the arriving customer at the (new) end of time_arri

15、val. */ time_arrivalnum_in_q = sim_time; Function Arriveelse /* Server is idle, so arriving customer has a delay of zero. (The following two statements are for program clarity and do not affect the results of the simulation.) */ delay= 0.0; total_of_delays += delay; /* Increment the number of custom

16、ers delayed, and make server busy. */ +num_custs_delayed; server_status = BUSY; /* Schedule a departure (service completion). */ time_next_event2 = sim_time + expon(mean_service); Function Arrive (continued)Flowchart for departure routinevoid depart(void) /* Departure event function. */ int i; float

17、 delay; /* Check to see whether the queue is empty. */ if (num_in_q = 0) /* The queue is empty so make the server idle and eliminate the departure (service completion) event from consideration. */ server_status = IDLE; time_next_event2 = 1.0e+30; Function Departelse /* The queue is nonempty, so decr

18、ement the number of customers in queue. */ -num_in_q; /* Compute the delay of the customer who is beginning service and update the total delay accumulator. */ delay = sim_time - time_arrival1; total_of_delays += delay; /* Increment the number of customers delayed, and schedule departure. */ +num_cus

19、ts_delayed; time_next_event2 = sim_time + expon(mean_service); /* Move each customer in queue (if any) up one place. */ for (i = 1; i = num_in_q; +i) time_arrivali = time_arrivali + 1; Function Depart (continued)Function Reportvoid report(void) /* Report generator function. */ /* Compute and write e

20、stimates of desired measures of performance. */ fprintf(outfile, nnAverage delay in queue%11.3f minutesnn, total_of_delays / num_custs_delayed); fprintf(outfile, Average number in queue%10.3fnn, area_num_in_q / sim_time); fprintf(outfile, Server utilization%15.3fnn, area_server_status / sim_time); f

21、printf(outfile, Time simulation ended%12.3f minutes, sim_time);Function update_time_avg_statsvoid update_time_avg_stats(void) /* Update area accumulators for time- average statistics. */ float time_since_last_event; /* Compute time since last event, and update last-event-time marker. */ time_since_l

22、ast_event = sim_time - time_last_event; time_last_event= sim_time; /* Update area under number-in-queue function. */ area_num_in_q+= num_in_q * time_since_last_event; /* Update area under server-busy indicator function. */ area_server_status += server_status * time_since_last_event;3.5实体到达模式与泊松过程n平稳

23、泊松过程在(t, t+s) 内到达的实体数k的概率为: N(t)表示在(0, t)区间内到达实体的个数 为到达速率!)()()(ksektNstNPks0, 0,2 , 1 , 0:stk其中3.5实体到达模式与泊松过程n定理如果 N(t) , t0是具有速率的泊松过程，那么其相应的到达时间间隔A1, A2 是具有均值1/ 的IID指数随机变量。n平稳泊松过程到达时间间隔服从指数分布, 其密度函数为：其中 为到达时间间隔均值。f teett( )/1()t 0/1Function Exponfloat expon(float mean) /* Exponential variate gener

24、ation function. */ /* Return an exponential random variate with mean mean. */ return -mean * log(lcgrand(1);Random-number generateThe default seeds for all 100 streams)1)31, 2()(mod(1630360016(powiZiZstatic long zrng = 1, 1973272912, 281629770, 20006270,1280689831,2096730329,1933576050, 913566091, 2

25、46780520,1363774876, 604901985,1511192140,1259851944, 824064364, 150493284, 242708531, 75253171,1964472944,1202299975, 233217322,1911216000, 726370533, 403498145, 993232223,1103205531, 762430696,1922803170,1385516923, 76271663, 413682397, 726466604, 336157058,1432650381,1120463904, 595778810, 877722

26、890,1046574445, 68911991,2088367019, 748545416, 622401386,2122378830, 640690903, 1774806513,2132545692,2079249579, 78130110, 852776735,1187867272, 1351423507,1645973084,1997049139, 922510944,2045512870, 898585771, 243649545,1004818771, 773686062, 403188473, 372279877,1901633463, 498067494,2087759558

27、, 493157915, 597104727,1530940798,1814496276, 536444882,1663153658, 855503735, 67784357,1432404475, 619691088, 119025595, 880802310, 176192644,1116780070, 277854671,1366580350, 1142483975,2026948561,1053920743, 786262391,1792203830,1494667770, 1923011392,1433700034,1244184613,1147297105, 539712780,1

28、545929719, 190641742,1645390429, 264907697, 620389253,1502074852, 927711160, 364849192,2049576050, 638580085, 547070247 ;Random-number generateFor example: lcgrand(1)float lcgrand(int stream).zi = zrngstream;.Simulation Output and DiscussionSingle-server queueing systemMean interarrival time 1.000 m

29、inutesMean service time 0.500 minutesNumber of customers 1000Average delay in queue 0.430 minutesAverage number in queue 0.418Server utilization 0.460Time simulation ended 1027.915 minutesSimulation Output and DiscussionuThe results are functions of the input parameters:The mean interarrival timeThe

30、 mean service timeThe n=1000 stop rule and also determined by the numbers the random-number generated(with another “seed” or “stream”).Alternative Stopping RulesuTwo types of terminationWhen the number of customers delayed became equal to 1000.The final value of the simulation clock is a random vari

31、able.After some fixed amount of time.The number of customers delayed is a random variable.The idea can be implemented in the programs by making changs to the following routines:The main programThe initialization routineThe report generatorAlternative Stopping Rules/* External definitions for single-

32、server queueing system, fixed run length. */#include #include #include lcgrand.h /* Header file for random-number generator. */#define Q_LIMIT 100 /* Limit on queue length. */#define BUSY 1 /* Mnemonics for servers being busy */#define IDLE 0 /* and idle. */int next_event_type, num_custs_delayed, nu

33、m_events, num_in_q, server_status; /* num_delays_required */ float area_num_in_q, area_server_status, mean_interarrival, mean_service, sim_time, time_arrivalQ_LIMIT + 1, time_end, time_last_event, time_next_event4 , total_of_delays; FILE *infile, *outfile;void initialize(void);void timing(void);void

34、 arrive(void);void depart(void);void report(void);void update_time_avg_stats(void);float expon(float mean);Alternative Stopping Rulestime_next_eventnext_event_type: next_event_type=1,2,31arrive();2depart();3report(); /the loop keeps repeating itself as long as the type of event just executed is not

35、3; after a type 3 event is chosen for execution , the loop ends and the simulation stops.main() /* Main function. */ /* Open input and output files. */ infile = fopen(mm1alt.in, r); outfile = fopen(mm1alt.out, w); /* Specify the number of events for the timing function. */ num_events = 3; /* Read in

36、put parameters. */ fscanf(infile, %f %f %f, &mean_interarrival, &mean_service, &time_end); /* Write report heading and input parameters. */ fprintf(outfile, Single-server queueing system with fixed run); fprintf(outfile, lengthnn); fprintf(outfile, Mean interarrival time%11.3f minutesnn,

37、 mean_interarrival); fprintf(outfile, Mean service time%16.3f minutesnn, mean_service); fprintf(outfile, Length of the simulation%9.3f minutesnn, time_end); /* Initialize the simulation. */ initialize();Main function /* Run the simulation until it terminates after an end-simulation event (type 3) oc

38、curs. */ do timing(); /* Determine the next event. */ update_time_avg_stats(); /* Update time-average statistical accumulators. */ switch (next_event_type) /* Invoke the appropriate event function. */ case 1: arrive(); break; case 2: depart(); break; case 3: report(); break; /* If the event just exe

39、cuted was not the end-simulation event (type 3), continue simulating. Otherwise, end the simulation. */ while (next_event_type != 3); fclose(infile); fclose(outfile); return 0;Main function (continued)void initialize(void) /* Initialization function. */ /* Initialize the simulation clock. */ sim_tim

40、e = 0.0; /* Initialize the state variables. */ server_status = IDLE; num_in_q = 0; time_last_event = 0.0; /* Initialize the statistical counters. */ num_custs_delayed = 0; total_of_delays = 0.0; area_num_in_q = 0.0; area_server_status = 0.0; /* Initialize event list. Since no customers are present, the departure (service completion) event is eliminated from consideration. The