clear

clc

%*****************************************

%初始化顾客源

%*****************************************

%总仿真时间

Total_time = 10

%队列最大长度

N = 10000000000

%到达率与服务率

lambda = 10

mu = 6

%平均到达时间与平均服务时间

arr_mean = 1/lambda

ser_mean = 1/mu

arr_num = round(Total_time*lambda*2)

events = []

%按负指数分布产生各顾客达到时间间隔

events(1,:) = exprnd(arr_mean,1,arr_num)

%各顾客的到达时刻等于时间间隔的累积和

events(1,:) = cumsum(events(1,:))

%按负指数分布产生各顾客服务时间

events(2,:) = exprnd(ser_mean,1,arr_num)

%计算仿真顾客个数，即到达时刻在仿真时间内的顾客数

len_sim = sum(events(1,:)<= Total_time)

%*****************************************

%计算第 1个顾客的信息

%*****************************************

%第 1个顾客进入系统后直接接受服务，无需等待

events(3,1) = 0

%其离开时刻等于其到达时刻与服务时间之和

events(4,1) = events(1,1)+events(2,1)

%其肯定被系统接纳，此时系统内共有

%1个顾客，故标志位置1

events(5,1) = 1

%其进入系统后，系统内已有成员序号为 1

member = [1]

for i = 2:arr_num

%如果第 i个顾客的到达时间超过了仿真时间，则跳出循环

if events(1,i)>Total_time

break

else

number = sum(events(4,member) >events(1,i))

%如果系统已满，则系统拒绝第 i个顾客，其标志位置 0

if number >= N+1

events(5,i) = 0

%如果系统为空，则第 i个顾客直接接受服务

else

if number == 0

%其等待时间为 0

2009.1516

%PROGRAMLANGUAGEPROGRAMLANGUAGE

events(3,i) = 0

%其离开时刻等于到达时刻与服务时间之和

events(4,i) = events(1,i)+events(2,i)

%其标志位置 1

events(5,i) = 1

member = [member,i]

%如果系统有顾客正在接受服务，且系统等待队列未满，则 第 i个顾客进入系统

else len_mem = length(member)

%其等待时间等于队列中前一个顾客的离开时刻减去其到 达时刻

events(3,i)=events(4,member(len_mem))-events(1,i)

%其离开时刻等于队列中前一个顾客的离开时刻加上其服

%务时间

events(4,i)=events(4,member(len_mem))+events(2,i)

%标识位表示其进入系统后，系统内共有的顾客数

events(5,i) = number+1

member = [member,i]

end

end

end

end

%仿真结束时，进入系统的总顾客数

len_mem = length(member)

%*****************************************

%输出结果

%*****************************************

%绘制在仿真时间内，进入系统的所有顾客的到达时刻和离

%开时刻曲线图（stairs：绘制二维阶梯图）

stairs([0 events(1,member)],0:len_mem)

hold on

stairs([0 events(4,member)],0:len_mem,'.-r')

legend('到达时间 ','离开时间 ')

hold off

grid on

%绘制在仿真时间内，进入系统的所有顾客的停留时间和等

%待时间曲线图（plot：绘制二维线性图）

figure

plot(1:len_mem,events(3,member),'r-*',1: len_mem,events(2,member)+events(3,member),'k-')

legend('等待时间 ','停留时间 ')

grid on

function r=randexp(lambda)r = -lambda*log(rand)

说明同上

%mmm simulation in matlab

clcclear

ST_Idle=0

ST_Busy=1

EV_NULL=0

EV_Arrive=1

EV_Depart=2

EV_LEN=3

Q_LIMIT=1000

% next_event_type=[]

% next_depart=[]

% num_custs_delayed=[]

% num_delays_required=[]

% num_events=[]

% num_in_q=[]

% server_status=[]

% area_num_in_q=[]

% area_server_status=[]

% mean_interarrival=[]

% mean_service=[]

% sim_time=[]

% time_last_event=[]

% total_of_delays=[]

%

time_arrival=[]%到达时刻

time_next_event=zeros(1,EV_LEN)

%仿真参数

num_events=EV_LEN-1

num_server=3%M/M/m m=2

mean_interarrival=1

mean_service=.5

num_delays_required=2000 %

outfile=fopen('MM m.txt','w')

fprintf(outfile, 'Multiple-server queueing system\n\n')

fprintf(outfile, 'Mean interarrival time%11.3f minutes\n\n',mean_interarrival)

fprintf(outfile, 'Mean service time%16.3f minutes\n\n', mean_service)

fprintf(outfile, 'Number of servers%20d\n\n', num_server)

fprintf(outfile, 'Number of customers%14d\n\n', num_delays_required)

%part1 initialize

sim_time=0.0

%/* Initialize the state variables. */

server_status =zeros(1,num_server) %idle

num_in_q= 0

time_last_event = 0.0

%/* Initialize the statistical counters. */

num_custs_delayed = 0

total_of_delays= 0.0

total_of_time= 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_event(EV_Arrive) = sim_time + randexp(mean_interarrival)

time_next_event(EV_Depart) = 1.0e+230

time_depart=zeros(1,num_server)

next_depart=0

%part2

while (num_custs_delayed <num_delays_required)

%Run the simulation while more delays are still needed.

%/* Determine the next event. */

min_time_next_event = 1.0e+290

next_event_type = 0

% Determine m depart event

min_time_depart=1e290

next_depart=-1

for i=1:num_server

if(server_status(i)==1 &time_depart(i)<min_time_depart)

min_time_depart=time_depart(i)

next_depart=i

end

end

time_next_event(2)=min_time_depart

%

%/* Determine the event type of the next event to occur. */

for i = 1: num_events

if (time_next_event(i) <min_time_next_event)

min_time_next_event = time_next_event(i)

next_event_type = i

end

end

%/* Check to see whether the event list is empty. */

if (next_event_type == 0)

%/* The event list is empty, so stop the simulation. */

fprintf(outfile, '\nEvent list empty at time %f', sim_time)

exit(1)

end

%/* The event list is not empty, so advance the simulation clock. */

sim_time = min_time_next_event

%/* Update time-average statistical accumulators. */

double time_since_last_event

%/* Compute time since last event, and update last-event-time marker. */

time_since_last_event = sim_time - time_last_event

time_last_event = sim_time

%/* Update area under number-in-queue function. */

area_num_in_q=area_num_in_q + num_in_q * time_since_last_event

%/* Update area under server-busy indicator function. */

for i=1:num_server

area_server_status =area_server_status + server_status(i) * time_since_last_event

end

%/* Invoke the appropriate event function. */

%arrival

if(next_event_type==EV_Arrive)

double delay

%/* Schedule next arrival. */

time_next_event(1) = sim_time + randexp(mean_interarrival)

%/* Check to see whether server is busy. */

s_idle=-1

for i=1:num_server

if (server_status(i) == ST_Idle)

s_idle=i

break

end

end

%/* all Server is busy, so increment number of customers in queue. */

if(s_idle== -1 )

num_in_q=1+num_in_q

%/* Check to see whether an overflow condition exists. */

if (num_in_q >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)

end

%/* There is still room in the queue, so store the time of arrival of the arriving customer at the (new) end of time_arrival. */

time_arrival(length(time_arrival)+1)=sim_time

else

%/* 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 =total_of_delays + delay

%/* Increment the number of customers delayed, and make server busy. */

num_custs_delayed = 1 + num_custs_delayed

server_status(s_idle) = ST_Busy

%/* Schedule a departure (service completion). */

time_depart(s_idle) = sim_time + randexp(mean_service)

end %if (server_status == ST_Busy)

%depart

else

double 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(next_depart) = ST_Idle

time_depart(next_depart) = 1.0e+230

%check_depart()

min_time_dapart=1e290

next_depart=-1

for i=1:num_server

if(server_status(i)==1 &time_depart(i)<min_time_depart)

min_time_depart=time_depart(i)

next_depart=i

end

end

time_next_event(2)=min_time_depart

else

%/* The queue is nonempty, so decrement the number of customers in queue. */

num_in_q=num_in_q-1

%/* Compute the delay of the customer who is beginning service and update the total delay accumulator. */

delay = sim_time - time_arrival(1)

total_of_delays =total_of_delays + delay

%/* Increment the number of customers delayed, and schedule departure. */

num_custs_delayed = 1 + num_custs_delayed

serv_time=randexp(mean_service)

time_depart(next_depart) = sim_time + serv_time

total_of_time = total_of_time + delay + serv_time

%check_depart()

min_time_dapart=1e290

next_depart=-1

for i=1:num_server

if(server_status(i)==1 &time_depart(i)<min_time_depart)

min_time_depart=time_depart(i)

next_depart=i

end

end

time_next_event(2)=min_time_depart

%/* Move each customer in queue (if any) up one place. */

tempForPop=time_arrival(2:length(time_arrival))

time_arrival=tempForPop

end %if (num_in_q == 0)

end %if(next_event_type==EV_Arrive)

end %while

%%%%%%%%%% part 3

%/* Invoke the report generator and end the simulation. */

fprintf(outfile, '\n\nAverage delay in queue%11.3f minutes\n\n',total_of_delays / num_custs_delayed)

fprintf(outfile, '\n\nAverage delay in system%11.3f minutes\n\n',total_of_time / num_custs_delayed)

fprintf(outfile, 'Average number in queue%10.3f\n\n',area_num_in_q / sim_time)

fprintf(outfile, 'Server utilization%15.3f\n\n',area_server_status / sim_time)

fprintf(outfile, 'Time simulation ended%12.3f minutes', sim_time)

fclose(outfile)

B

排队 不排队

A 排队 pr ,(1-p)r p(r-l) , (1-p)r+pl

不排队 pr+(1-p)l , (1-p)(r-l) pr-l, (1-p)r-l

从表上可以看出，A与B之间的博弈不存在纯粹均衡策略，因为(Y表示排队，N表示不排队，（Y，N）表示A选择Y，B选择N时的收益）

对A来说，（N，Y)>Y,Y>Y,N>N,N

对B来说，Y,N>Y,Y>N,Y>N,N

所以此博弈不存在纯粹均衡。

但是这里是一次博弈的情形，说明在只有1次博弈机会的情况下，每个人都有选择不排队的动力，而最终对所有人自己造成伤害。分析到这里并不能说明我们都应该排队或者插队，只能告诉我们在别人排队的时候插队是有利可图的。

如果需要推导出我们应该排队的结论，你可以从多次重复博弈的总收益来思考。

本人能力有限，希望能对你有所帮助。