新建WinRARZIP压缩文件

多个服务台的排队模型的蒙特卡罗模拟 已有 1174 次阅读 2009-11-1 21:48 | 个人分类 :未分类 | 系统分类 :教学心得 % 多个服务台的排队模型的蒙特卡罗模拟（以 4 个为例） clear TjMm=[]; % 顾客到达的时间间隔序列 sTj=0; %顾客到达时间序列 while (sTj<=480) Tjp=exprnd(10);% 参数为 0.1 的指数分布 % Tfp1=unifrnd(4,15); %Tfp2=unifrnd(4,15); TjMm=[TjMm;Tjp]; sTj=sTj+Tjp; end n=length(TjMm); n=n-1;% 一天中顾客的个数 for i=1:n TjM(i)=TjMm(i); end TfM=unifrnd(50,80,[4,n]); %4 个服务台的服务时间序列， 50 到 80 的均匀分布 dengtime=zeros(1,n); % 顾客等待的时间序列 likaitime= zeros(1,n); % 顾客离开的时间序列 fuwutime= zeros(1,n); % 顾客服务的时间序列 daodatime(1)=TjM(1);% 第 1 个顾客到达的时间 fuwutime(1)=TfM(1,1); % 第 1 个顾客服务的时间 likaitime(1)=daodatime(1)+fuwutime(1); % 第 1 个顾客离开的时间 deskkongxiantime(1,1)=likaitime(1); % 第 1 个服务台空闲的时间 deskkongxiantime(2,1)=daodatime(1); % 第 2 个服务台空闲的时间 deskkongxiantime(3,1)=daodatime(1); % 第 3 个服务台空闲的时间 deskkongxiantime(4,1)=daodatime(1); % 第 4 个服务台空闲的时间 for i=2:n daodatime(i)=sum(TjM (1:i)); % 第 i 个人到达的时间 kongxianmin=min(deskkongxiantime(:,i-1));% 服务台最小的空间时间 for i1=1:4 if deskkongxiantime(i1,i-1)==kongxianmin s=i1;break % 空闲时间最长服务台 end end deskkongxiantime(:,i)=deskkongxiantime(:,i-1);% 第 i 个人时服务台空闲的时间（ s 例外） fuwutime(i)=TfM(s,i); % 第 i 个顾客服务的时间 if deskkongxiantime(s,i-1)<=daodatime(i) dengtime(i)=0; else dengtime(i)=deskkongxiantime(s,i-1)-daodatime(i);% 第 i 个顾客等待的时间 end likaitime(i)=daodatime(i)+fuwutime(i)+dengtime(i); % 第 i 个顾客离开的时间 deskkongxiantime(s,i)=likaitime(i);% 第 i 个人时第 s 个服务台空闲的时间 end pingjundengtime=sum(dengtime)/n; % 平均等待时间 n=2000; % 随机点数 (可增加点数 ) x=12+(62-12)*rand(1,n); % 产生 2000 个 12 到 62 的随机数 xx=12:2:62; % 画概率密度图的区间 nx=histc(x,xx); %计算 x 在 xx 每个小区间内的点数。 px=nx/n; sumpx=cumsum(px); subplot(1,2,1) bar(xx(1:end-1),px(1:end-1)); title('概率密度 ') subplot(1,2,2) plot(xx(1:end-1),sumpx(1:end-1)); title('累积概率密度 ') andn('state', sum(100*clock)); %利用时钟设置随机种子，这样每次产生的随机数就不同了 N = 1000000; % 设置样本个数 m = 0.43 + sqrt(0.008) * randn(N,1); n = 0.19 + sqrt(0.008) * randn(N,1); L = 0.5 + sqrt(0.008) * randn(N,1); %随机生成 m, n, L ，根据他们各自的分布 g = 16.4445 .* (m.*m.*L./(m.*m+n.*n))-8.2; %计算 g fprintf('the probability of g>0 is %.2f\n', sum(g>0)/length(g)); %统计 g>0 的概率，结果为 0。18 已经测试通过，可以运行，结果为 the probability of g>0 is 0.18 现在我通过纯数学的方式给出排队论的若干结论 （假定卷积符号是 @）（鉴于 M/M/n 理论推倒太变态暂不论，等下蒙卡罗特再说） （到达时间负指数分布 /服务时间负指数分布 /单台 FCFS M/M/1 ） （其实后面统一化为 X/Y/Z/A/B/C 模型才是正确说法） 到达率 namda 服务率 mu (t 时间间隔 n 人到达的概率 t>0) p(n,t)=(namda*t)^n*exp(-namda*t)/n! 到达间隔服从负指数分布，可以由顾客到达服从泊松分布推出 下面分析排队情况，鉴于排队过程稳定，人数不变，所以 p 与 t 无关，有 namda*p(0,t)=mu*p(1,t) namda*p(n-1,t)+mu*p(n+1,t)=namda*p(n,t)+mu*p(n,t) p(n,t) 的累加和为 1 解得 p(n,t)=(1-(namda/mu))*-(namda/mu)^n 平均队长就是 np(n,t) 的累加和，求和得 namda/(mu-namda) 下面就介绍怎么使用蒙卡罗特法拟合 先是 M/M1 首先先把时间分为片段，比如 1 秒钟，仿真时间为 1 天 你用 rand() 产生随机数与 p(n,t)=(namda*t)^n*exp(-namda*t)/n! 比较大小 这是分段函数，因为 p(n,t) 对 n 的累加和为 1，而 rand() 产生 0-1 的随机数 这就是假定来了的人得个数，而在 1 秒内被服务走的人数就是 mu*t 这样就得到了一个排队人数， 这个人数将变化 60*60*24 次 因为一天有这么多秒 你可以把它存在数组里，然后取中间或偏后的一些点的平均值 那就是你仿真得到的排队长 对于 M/M/n 其实也是一样 你排队的就是 n 个数组 来的人则是判定加入队最短的那个服务台 然后一样就可以了 添 加 函 数 newMM2andMM1， 参 数 的 解 释 也 在 源 程 序 里 说 明 了 。 则 按 意 输 入 meant=newMM2andMM1(1,1.5,2,0.2,1,300) ，得到结果为时间为 meant =2.2418 ； 各个人得到的结果是不太一样的， 是正确的， 因为程序模拟的随机分布的数， 是符合实际的。 源程序如下： function [meant]=newMM2andMM1(mean_arr,mean_serv,mean_serv1,mean_serv21,mean_se rv22,peo_num) %mean_arr, 到达的时间参数 %mean_serv,服务台 a 的服务时间参数 %mean_serv1,服务台 b 的服务时间参数 %mean_serv21,mean_serv22, 二级服务台的均匀分布的参数 %peo_num,总服务人数 nt=exprnd(mean_arr,1,peo_num); %各顾客到达时间间隔服从指数分布 state_a=zeros(3,peo_num); %用一个三行矩阵示 a 台每个顾客的状态 %三行依次为：到达时间间隔，服务时间，等待时间 state_b=zeros(3,peo_num); %用一个三行矩阵表示 b 台每个顾客的状态 %三行依次为：到达时间间隔，服务时间，等待时间 state_a(2,:)=exprnd(mean_serv,1,peo_num); %生成 a 台各顾客服务时间的矩阵 state_b(2,:)=exprnd(mean_serv1,1,peo_num); %生成 b 台各顾客服务时间的矩阵 state_a(3,1)=0; state_b(3,1)=0; a=1;%a 台服务的人数 b=1;%b 台服务的人数 arr_time=cumsum(nt); %到达时间由时间间隔变成连续时间 state_b(1,1)=arr_time(1); state_a(1,1)=arr_time(2); %state(1,:)=arr_time; lea_time_a(1)=sum(state_a(:,1)); %先计算前 1 名顾客的离开时间 lea_time_b(1)=sum(state_b(:,1)); %先计算第 2 名顾客的离开时间 for i=3:peo_num if lea_time_a(a)