排队系统仿真python 排队系统建模仿真

排队模型和排队系统仿真

Gary哥哥 2021.1.31

排队论又称随机服务系统，是研究系统随机聚散现象和随机服务系统工作过程的数学理论和方法，是运筹学的一个分支。排队论的基本思想是 1909 年丹麦数学家 A.K. 埃尔朗在解决自动电话设计问题时开始形成的，当时称为话务理论。现实生活中如排队买票、病人排队就诊、轮船进港、高速路上汽车排队通过收费站、机器等待修理等都属于排队论问题。

定义

• 通过对服务对象到来及服务时间的统计研究，得出这些数量指标（等待时间、排队长度、忙期长短等）的统计规律
• 然后根据这些规律来改进服务系统的结构或重新组织被服务对象，使得服务系统既能满足服务对象的需要，又能使机构的费用最经济或某些指标最优。

数模应用：

• CUMCM 2009B 的眼科病床的合理安排问题
• MCM 2005B 收费站最佳配置问题
• ICM 2017D 机场安检问题

…

基本构成和指标

基本构成

• 输入过程：描述顾客按照怎样的规律到达排队系统。顾客总体（有限/无限）、到达的类型（单个/成批）、到达时间间隔。
• 排队规则：指顾客按怎样的规定次序接受服务。常见的有等待制、损失制、混合制、闭合制。
• 服务机构：服务台的数量; 服务时间服从的分布。

数量指标

• 队长：系统中的平均顾客数（包括正在接受服务的顾客）。
• 等待队长：系统中处于等待的顾客的数量。
• 等待时间：等待时间包括顾客的平均逗留时间。
• 忙期：连续保持服务的时长。

等待机制模型

A/B/C/n

A 输入过程，B 服务时间，C 服务台数，n 系统容量。

example：

• M/M/S/∞

• 输入过程是 Poisson 流
• 服务时间服从负指数分布
• 系统有 S 个服务台平行服务（并行）
• 系统容量为无穷大的等待制排队系统（容量无穷大）

案例

服务台不再是一个之后

S==1时：

来访人员按照 Poisson 流到达，到达速率为 µ = 20 人/小时。接待人员的服务速率间服 λ = 9 人/小时的负指数分布。为使来访问者等待不超过半小时，最少应配置几名接待员？

lambda = 20; mu = 9; s = 3;
rho = lambda/(s*mu);
k=0:s-1
p0 = 1./( sum((s*rho).^k./factorial(k)) + ... 
(s*rho)^s/(factorial(s)*(1-rho)) );
Ls = s*rho + (s*rho)^s*rho/(factorial(s)*(1-rho)^2)*p0;
Ws = Ls/lambda;
Wq = Ws - 1/mu

• 应用，收费站

其他模型

系统容量K

• 损失制模型 M/M/S/S

• 顾客到达服从泊松分布，服务台服务时间服从负指数分布，当 S 个服务台被占用后，顾客自动离开，不再等待。

• 混合制模型 M/M/S/K

• 顾客到达服从泊松分布，服务台服务时间服从负指数分布，系统容量为 K，当 K 个位置被占用时，顾客自动离开。

• 闭合制模型 M/M/S/K/K

• 顾客到达服从泊松分布，服务台服务时间服从负指数分布，系统容量和潜在的顾客数都为 K。

模拟

单服务台

• 服务时刻(i) = max { 到达时刻(i), 离开时刻(i -1)}
• 离开时刻(i) = 服务时刻(i) -服务时长(i)
• 等待时长(i) = 离开时刻(i) -到达时刻(i)

多服务台

尽快完成服务的目的

• 服务时刻(i) = max { 到达时刻(i), min*{服务台空闲时刻}*}
• 所使用服务台(i) = k,其中 k 使 服务台空闲时刻(k) = min (第k个服务台空闲)

• 离开时刻(i) = 服务时刻(i) + 服务时长(i)
• 服务台空闲时刻(k) = 离开时刻(i)
• 等待时长(i) = 离开时刻(i) -到达时刻(i)

自动取款问题

• 银行计划安置取款机, A 机价格和平均服务率都是 B 机的 2 倍. 应购置 1 台 A 机还是 2 台 B 机?

MM1无穷的问题：

单服务台

n = 100000;  %模拟顾客数量
mu = 1; muA = 0.9; %到达率和服务率
tarr = cumsum(exprnd(mu,1,n));% 到达时刻（指数分布）
tsrv = exprnd(muA,1,n);%服务时间长度
tsta = zeros(1,n);%初始化服务时刻
tlea = zeros(1,n);%初始化离开的时刻
twat = zeros(1,n);%初始化等待时长
tsta(1) = tarr(1);%首位顾客服务时刻=到达时刻
tlea(1) = tsta(1) + tsrv(1);%首位顾客离开时刻
twtime(1) = tlea(1) - tarr(1);%首位顾客等待时长
for i = 2:n
    %服务时刻=max{到达时刻，上一个顾客离开时刻}
    tsta(i) = max(tarr(i),tlea(i-1));
    tlea(i) = tsta(i) + tsrv(i);% 离开时刻=服务时刻+服务时长
    twat(i) = tlea(i) - tarr(i);% 等待时长=离开时刻-到达时刻
end
hist(twat)
sum(twat)/n

多服务台

• 顾客平均每分钟到达 1 位，A 型机的平均服务时间为 0.9, B型机为 1.8 分钟, 顾客到达间隔和服务时间都服从指数分布.

n = 100000; 
mu = 1; 
muB = 1.8;
tarr = cumsum(exprnd(mu,1,n));
tsrv = exprnd(muB,1,n);
tsta = zeros(1,n);%服务时刻
tlea = zeros(1,n);%离开时刻
twat = zeros(1,n);%等待时刻
last = [0 0];%服务台结束服务时刻
for i = 2:n
    [minemp, k] = min(last); %找出最快结束服务的服务台时刻
    tsta(i) = max(tarr(i),minemp);
    tlea(i) = tsta(i) + tsrv(i);
    last(k) = tlea(i);
    twat(i) = tlea(i) - tarr(i);
end
hist(twat)
sum(twat)/n

多服务台就在last中加就好了

数学建模真题

银行服务问题

银行经理正试图通过提供更好的服务来提高顾客满意度. 管理层期待顾客平均等待时间小于 2 分钟, 平均队列长度 (等待队列的长度) 是 2 人或者更少. 银行估算, 每天大约为 150 名顾客提供服务. 现有的到达和服务时间如下表所示:

根据经理的指引, 确定当前顾客是否对服务满意. 如果不满意, 通过模型对服务进行微小的改变, 以达到经理的目标.

除了比赛格式的论文以外, 给经理写一封简短的, 非技术性的, 大约 1-2页的信, 给出你的最终建议.

取时间间隔的方式比较巧妙

n = 150;
ta = [5 4 3 2 1 0]; pa = [0.05 0.25 0.35 0.10 0.15 0.10];
ts = [ 4 3 2 1 ]; ps = [ 0.15 0.40 0.20 0.25 ];
pacum = cumsum(pa);
pscum = cumsum(ps);
Tarrival = rand(1,n);
for i = 1:length(pa)
	Tarrival(Tarrival<pacum(i)) = ta(i);%一次判断，上面的处理
end
%服务的时间长度
Tarrival = cumsum(Tarrival);
Tservice = rand(1,n);
for i = 1:length(ps)
	Tservice(Tservice<pscum(i)) = ts(i);
end

Tstart = zeros(1,n); Tleave = zeros(1,n);
Twait = zeros(1,n); line = zeros(1,n);%队长

%第一个人
Tstart(1) = Tarrival(1);
Tleave(1) = Tstart(1) + Tservice(1);
Twait(1) = Tleave(1) - Tarrival(1) - Tservice(1);
line(1) = 0;

for i = 2:n
	Tstart(i) = max(Tleave(i-1), Tarrival(i));
	Tleave(i) = Tstart(i) + Tservice(i);
	Twait(i) = Tleave(i) - Tarrival(i) - Tservice(i);
	%计算队长
	k = i-1;
	%找到前面为离开的那个人，直到不在队列的人
	while ( k>0 )&&( Tarrival(i)<Tleave(k) )
		line(i) = line(i) + 1;
		k = k - 1;
	end
end

优化机场安检

ICM2017-D

• 串并联问题

• 建立一个或多个模型，研究旅客通过安检口的流量，确定瓶颈，明确判断当前流程问题区域位置。
• 设计两个或更多对现有系统德潜在改进，提高旅客通信，减少等待时间。模拟这些变化展示改进如何影响流程。

排队系统

并联后串联部分

function [tlea, twat, qlen] = mms(tarr, type, mus)
% MMS Stochastic simulation for M/M/c queue
%
% [tlea, twat, qlen] = mms(tarr, type, mus)
%     tarr = arrival time of customers
%     type = customer type parameters
%     mus  = serere rate of servers
%     tlea = leaving time of servers
%     twat = waiting time of servers
%     qlen = length of the queue (length of the waiting line) for customers
% 
% Zhou Lvwen: [e-mail]
% January 21, 2017

narr = length(tarr);        % number of customers
nsvr = length(mus);         % number of servers

% last time at which a customer left a particular server，标志服务台
last = zeros(nsvr,1);

[tsta, tlea, twat, qlen] = deal(zeros(narr,1));%deal全部赋值为0

rndm = zeros(nsvr,narr);    % rndm(k,i) = service time for i-th customer，生成每个人你的服务长度
for k = 1:nsvr; rndm(k,:) = exprnd(mus(k)*type); end
%套路操作
for i = 1:narr
    % find booth service was/will be emptied soonest and record
    [minemp, ksvr(i)] = min(last); 
    
    % start time = max{arrival time, minemp}
    tsta(i) = max(tarr(i), minemp); 
    
    % severe time = exponential random number with mean parameter mu
    tsvr(i) = rndm(ksvr(i),i);
    
    % leaving time = start time + service time
    tlea(i) = tsta(i) + tsvr(i);
    
    % last time of k-th server = leaving time of i-th customer 
    last(ksvr(i)) = tlea(i);
    
    % waiting time = leaving time - arrival time
    twat(i) = tlea(i) - tarr(i);
    
    % queue length for i customer
    j = i - 1;
    while j>0 && tarr(i)<tlea(j)
        if ksvr(j)==ksvr(i); qlen(i) = qlen(i) + 1; end
        j = j - 1;
    end
end

main.m

n1 = 2;  n2 = 3; n3 = 3;
mu1 = 12; mu2 = 9; mu3 = 16;
muR = 10; muB = 13;

nR = ceil(24*3600/muR); nB = ceil(24*3600/muB);
tArrR = cumsum(exprnd(muR,nR,1));
tArrB = cumsum(exprnd(muB,nB,1)); 
tArr = [tArrR; tArrB];
type = [0.8*ones(nR,1); 1.2*ones(nB,1)];

[tLeaR, tWatR, qLenR] = mms(tArrR, ones(nR,1), mu1*ones(n1,1));
[tLeaB, tWatB, qLenB] = mms(tArrB, ones(nB,1), mu2*ones(n2,1));

[tArrG, order] = sort([tLeaR; tLeaB]);
[tLeaG, tWatG, qLenG] = mms(tArrG, type(order), mu3*ones(n3,1));
tLeaG(order) = tLeaG;
tWatG(order) = tWatG;
qLenG(order) = qLenG;


figure('position',[50,50,1200,600])
subplot(2,3,1); hist(qLenR); ylabel('Frequency'); 
xlabel('length of the waiting line'); title('Red')
subplot(2,3,4); hist(tWatR); ylabel('Frequency'); 
xlabel('waiting time'); title('Red')


subplot(2,3,2); hist(qLenB); ylabel('Frequency');
xlabel('length of the waiting line'); title('Blue')
subplot(2,3,5); hist(tWatB); ylabel('Frequency'); 
xlabel('waiting time'); title('Blue')

subplot(2,3,3); hist(qLenG); ylabel('Frequency');
xlabel('length of the waiting line'); title('Green')
subplot(2,3,6); hist(tWatG); ylabel('Frequency'); 
xlabel('waiting time'); title('Green')

完结散花