水房拥挤问题.doc

目录一、问题的提出 2二、数据和建模设 2三、模型的建立与解 5四、模型分析 74.1 灵敏度分析 74.2 系统的最优化 84.3 两种排队方式的比较 10五、本模型小结11六、参考文献 11七、附录（MATLAB程序） 12一：问题的提出北京师范大学在校学生约六千余人,只有一个开水房,供水时间为早晨6:30到8:00,中午110到12:30,下午17:00到130水房内共有22个水龙头供人们使用,可排队的空地面积约10平方米.烧开水的锅炉容量较小,使得送水管道较细,水流量受到一定的限制,且水管易被水垢堵塞,使水流减小甚至状如细线。水房内常有排队现象，人们也抱怨水房的拥挤，这就给我们提出了问题：水房设计是否合理，为什么拥挤，拥挤程度如何？怎样进行改进？显然水房就是一个随即服务系统，可以应用排队论的方法对系统运行状态做定量的描述。但首要的是收集数据，提出合理假设，选择排队模型及估计参数；然后再对模型进一步讨论以描述拥挤现象，分析拥挤原因，研究改进措施。为此本文将综合运用多种工具建立一个完整的数学模型，用以讨论着一水房问题。二：数据和建模假设经过一周（星期一到星期五）在 11：40 到 12：20 之间的观测，得到了顾客到达的情况的数据；在管道通畅及堵塞严重时分别观测得到服务时间的数据。而且发现管道通畅是几无人排队等待，堵塞时水房内十分拥挤，许多人提着空壶离去！表1 每10秒到达人数及频数到达人数（个/10秒）012345678相应的频数6613213111050221043表2 管道通畅时顾客打水时间及频数打水时间（秒）30354045505560657075808595105125155245相应的频数12322139453542211表3 管道阻塞使顾客打水时间及频数打水时间（秒）3040455565707580859095100105110125130相应的频数2133234111412213打水时间（秒）135140145155160175185190200205215240255265300相应的频数222211211211111基于上述观察提出以下假设：1） 由于水房开放时间有限，可以认为顾客源无限，顾客单个到来且相互独立。并假设顾客流在观测期间平稳，不考虑出现高峰期与空闲期的可能性。2） 排队方式为单一队列的等候制，先到先服务。虽然水房内有多个服务台，每服务台都有自己的队列，但同时顾客总是自由转移到最短的队列上，不可能出现有顾客排队而服务器空闲的情况。这样就队列长度变化而言，这种排队方式与单一队列并没有区别。本文最后对两种排队方式的比较也表明这一假设是合理的。3） 顾客流满足参数为A的f泊松分布。其中人为单位时间到达的顾客平均数。4） 原水房共有22个并联的服务台（水龙头），舍每个服务台的服务时间服从某个相同的分布，T和分别是服务时间的均值和均方差，=/ t为偏离系数。由于锅炉及输水管容量的限制，使t依赖于正在进行服务的水龙头个数m，设此时平均服务时间t(m)。且存在一临界值当mm0 时，t(m)为常数t0;mm0时，管道中的水便分给 m 个龙头流出，从而 t(m) t0,且 t(m)是 m 的单增函数。 5） 并估单位时间为 10 秒。上面所提出的假设 1），2），4）都是合理的，下面对假设 3）进行拟合优度检验，观测数据（表 1表 3）得出有关参数的最大拟然估计；顾客到达率l =2.2 然后做分布拟合的c 2- 检验，取c 2统计量为：c 2=i=0k(ni-Npi)2Npi其中k为组数，N为样本数，n i与Npi分别为第 i组实测频数与拟合的理论频数。由20.052(k-1)只能说明我们没有足够（95%）的把握拒绝原假设，为进一步检验实际数据与理论分布的符合程度，定义拟合优度：P(k,a)=P(2 (k -1) a)(其中a为2 统计量的具体值），表明假设成立时，偏差2取值不小于a的可能性。显然P(k,a)越大，表明实际数据与理论分布拟合得越好，但因2-分布是中间隆起的单峰分布，典型的符合良好得情况是P(k,a)在0.30.7左右的范围内。 由表一知 N=528, k=8, l=2.2,于是： 2=6.1190.052(7)=14.067，0.5p(8，6.119)m0则相当于服务台之间可以相互帮助的服务系统，平均服务时间 t 为正在服务的服务台数 m 的函数。考虑一简单情形：当 m m0时，t(m)=t0;当 m0 m 22=c,因此现有的水房系统服从M/G/22的排队模型。这种非马氏过程的排队模型，目前还不存在任何简单的解。Tijms等给出了Lq及P0的近似解，公式如下： P0=i=0c-1(t)2i!+(t)2c!(1-)-1;P0=1-1c!1-(t)cp0; Lq=1-1c!1-(t)cp01+22;另外有公式L=Lq+c, Wg=Lq/，W=L/.对M/Er/c及M/D/c的近似解与精确解进行比较，即r不是很小时，相对误差小于5%，而且公式十分简单，我们只需要知道服务时间的均值与方差就行了，而不需要了解具体的分布，这在实际应用中十分方便。(注意公式中要求r1,否则系统永远不能到稳定状态，队长将趋于无穷大)对水房系统，l=2.2,c=22, 当管道通畅时，t1=7.58,1=3.60, r=0.7581,水房爆满，许多人提着分析将试我们进一步了解拥挤情况，拥挤原因以及缓解的办法。四：模型分析：4.1 灵敏度分析由近似公式可知直接影响系统个运行指标的参数为c,r=lt及g。下面将分析参数变动对参数结果的影响。从中可以得到对拥挤程度的描述，并分析缓解拥挤的方法。 显然，Lq与g的关系最为简单，为一元二次函数关系。当g=0时，平均队长最小。因此，若顾客流服从均匀分布，便可使一定的r 及c下的队长达到最小，g 对P0没有影响。若对原有的服务台数c=22,在假设g =1t1=0.455保持不变，便可得到r 与Lq，P0的关系曲线（图一）。若再令l=2.2则图1反映了平均服务时间t与Lq及P0的关系；而t是由管道堵塞程序及锅炉内水位影响的，因此，图1实际上反映了服务台（水龙头），在不同状态下水房的拥挤程度。P011616.51717.51818.51919.52020.52100.5图一 1616.51717.51818.51919.52020.521050Lq100C=22X=*t1617181920212223242500.5P0 116171819202122232425051015Lq20=2.2t=7.58m0=23图二22由图 1，2 可知，当 r/ c 趋于 1 时，队长将趋于无穷，增加服务台或者 r 必可缓解拥挤 状态。但增加服务台不仅会占用空间，还会使空闲概率增大造成浪费；由于顾客到达情况是 不容易控制的，为降低 r 须通过经常维修以保证管道畅通或增锅炉的容量及高度以提高服务速度，这就导致了设备管理费用的增加。怎样以最少的支出获得最大的收益呢?这就是下面的优化问题4.2 系统的最优化服务系统的最优化分为系统设计的最优化和系统控制的最优化。水房的最优化设计实际上就是要确定最佳临界服务台数m0,并使系统的实际服务台为m0个，因为若服务台少于m0，则管道内水没有被充分利用，而多于m0服务台对系统运营是没有影响的。如果从费用结构考虑，通常服务系统中有两种费用，一是每位顾客在系统中逗留单位时间的等候费；一种是每个服务台单位时间的服务台。若能给出平均总费用的计算公式，就可确定使之极小化的最优服务台数。但这里如何估算等待费方面缺乏必要的数据和材料，不可能对顾客的等待时间用固定的货币值度量，因此应用费用的模型是不合适的。现在的另外一种方法是:在两种互相矛盾的度量（平均等待时间Wg和服务台空闲概率P0)之间取折中值。即对Wg和P0规定上限值a和b，然后再决定满足。这两个规定上限的服务台数，以之作为最佳临界服务台m0进一步也决定了锅炉及输水管的规模。由不同的服务台数c给出Wg和P0,得到了图3 ，称之为设计曲线。若规定a=6,b=0.4，则最优服务台数在t1=7.58与t2=9.902两种情况下分别为：c1=18, c2=23。但顾客到达率在不同时段内是不一样的，因此不存在一致最优的服务台数。在这种情况下，应该通过控制服务速率来达到系统最优化。对水房而言，在保证管道通畅的条件下，通过调节锅炉内水位高度便可实现对服务速度的控制，在缓解拥挤的同时，可以避免因水位过高而造成的浪费。为此，在图3中将Wg换成Lq,对不同的r 做设计曲线得到一族r1 =16.676及r 2=21.637相应的控制曲线，即为图3中的t1 =7.58，t2=9.902时的设计曲线。c固定，给出Lq和P0的上限a，b，对已知的某个l，找到适合a，b的r，从而由t=r/ c求出平均服务时间t,在据此调节服务设备即可。例如c=22,取 a=10，b=0.4，得到r=21，若l=2，则t=10.5.另外控制曲线族还可用于具有一定r值得多个服务员的服务系统的最优化设计。 4.3 两种排队方式的比较如果假定系统中 c 个队列间没有顾客转移，则每个队列平均到达率为l / c ，从而成为 c 个 M/G/1 系统，平均服务时间 t 不变，l=0.1，仍利用原公式计算得到多队时系统的运行指标。见下表 4： 表 4，单队与多队的比较 运行指标 服务时间 排队方式L=7.58单队0.29416.970.1340.945多对1.432248.414.30.998=9.837单队32.7354.3715.30.089多对35.35（每个子系统）799.0（整个子系统）353.50.307显然单队时等待队长，等待时间都比多队时低，而服务台利用率比多队的高。因此，具有明显的优越性。同时，多队时系统内平均人数在 t=7.58 时为 48.4 人，与实际情况不符，这说明所做的单队假设是合理的。如果将现行排队方式改为严格的单队，就可避免队列间拥挤碰撞，有助于改进水房的服务状况。同时可得出一个结论：建一个水房优于建二个，即联合服务优于独立通道服务。五：本模型小结本问题的提出完全出去平时生活中的观察和思考，自己在观察数据方面花了很多的时间，然后再建立模型的时候，还是遇到了不少的麻烦！例如在这个排队论中，自己并不是很了解，查阅很多排队论的资料，然后结合该问题运用！然后再图形分析上，遇到了很大的问题，怎么写程序，怎么比较有意义？以及最优化的提出，都花了相当的多的时间，但是经过近二个星期的工作，终于有所结果！最后分析的结果，我自己观察的结果是很符合的！很有成就感！最重要的是：体会了一项工作的全部过程！很有感觉！谢谢批阅！六:参考文献：作者名 苏兆龙 编著 题名 排队论基础 = Fundamentals of queueing theory / 苏兆龙编著 eng排队论基础 = Fundamentals of queueing theory / 苏兆龙编著 eng 出版发行 成都 : 成都科技大学出版社, 1998 作者名 孙荣恒 Sun Rong Heng 著 题名 排队论基础 Pai Dui Lun Ji Chu / 孙荣恒, 李建平著排队论基础 Pai Dui Lun Ji Chu / 孙荣恒, 李建平著 出版发行 北京 : 科学出版社, 七：附录计算程序如下：（MATLAB编程）1: 顾客到达率 l的计算程序a=zeros(1,66); b=ones(1,132); c=2.*ones(1,131); d=3.*ones(1,110); e=4.*ones(1,50); f=5.*ones(1,22); g=6.*ones(1,10); h=7.*ones(1,4); i=8.*ones(1,3); x=a,b,c,d,e,f,g,h,i; x1=mean(x) 我们存为 thesis0.m 结果为 ：x1 =2.1705 即：l =2.22:服务时间的均值和样本方差在管道畅通时，为：t1=7.58, s1=3.60（以十秒为单位）。计算程序如下： x=30 35 35 40 40 40 45 45 50 50 55 60 60 60 . 65 65 65 65 65 65 65 65 65 70 70 70 70 70 . 80 80 80 85 85 85 85 85 95 95 95 95 105 105 . 125 125 155 245 x1=mean(x) x2=std(x) 我们存为 thesis1.m 结果为：x1 =75.7609x2 =35.9933即： t1=7.58,s 1= 3.60(以10秒为单位）3:服务时间的均值和样本方差在管道堵塞时，t1=7.58,s1=3.60（以十秒为单位）。计算程序如下：x=30 30 40 45 45 45 55 55 55 65 65 70 70 70 75 75 . 75 75 80 85 90 95 95 95 95 100 105 105 110 110 . 125 130 130 130 135 135 140 140 145 145 155 155 . 160 175 185 185 190 195 200 205 205 215 240 255 . 265 300; x1=mean(x) x2=std(x) 我们存为 thesis2.m 即：t1= 12.2,s 1=6.40(以10秒为单位）4:模型求解的计算程序如下：c=22;s1=1;s2=0;p=0.758;a=2.2;t1=7.58 d=1 for m=1:21 d=m*d;s1=(a*t1)m/d+s1;endc1=22*d;s2=(a*t1)22/c1/(1-p); s=s1+s2;p0=1/sP0=1-1/c1/(1-p)*(a*t1)22*p0 结果如下： P0 = 0.8445 Lq =0.2940 L = 16.9700 Wq = 0.1336 W = 7.7136 x1=3.448; b=x1/t1; Lq=p/(1-p)/c1/(1-p)*(a*t1)22*p0*(1+b2)/2 L=Lq+c*p Wq=Lq/a W=L/a 我们存为 thesis4.m5:1）灵敏度分析的图形一的程序：x=16:0.1:21;c=22; s1=1; s2=0; b=0.445; p=x./c ;d=1; for m=1:21 d=m*d; s1=x.2/d+s1; end c1=22*d; s2=x.2/c1/(1-p); s=s1+s2; p0=1./s x1=x.22 x2=p0.*x1 n=1/c1./(1-p) P0=(1-n.*x2) j1=1./(1-p); j2=p.*j1; b1=(1+b2)/2; Lq=j2.*n.*x2*b1 subplot(2,1,1),plot(x,P0),axis(16,21,0,1),title （图 1-1) subplot(2,1,2),plot(x,Lq),axis(16,21,0,100),title(图 1-2) 我们存为 thesis5.m 2)灵敏度分析的图形二的程序x=16.676; c=16:1:25;s1=1; s2=0; b=0.455; for i=1:10 p(i)=x/c(i); endfor i=1:10 s1(i)=1; d(i)=1; for m=1:c(i)-1 d(i)=m*d(i); s1(i)=xm/d(i)+s1(i)end end for i=1:10 c1(i)=c(i)*d(i); f(i)=xc1(i); s2(i)=f(i)/c1(i)/(1-p(i); s(i)=s1(i)+s2(i); p0(i)=1/s(i); x1(i)=xc(i); x2(i)=p0(i)*x1(i); n(i)=1/c1(i)/(1-p(i); P0(i)=(1-n(i)*x2(i); j1(i)=1/(1-p(i); j2(i)=p(i)*j1(i); b1=(1+b2)/2;Lq(i)=j2(i)*n(i)*x2(i)*b1;end subplot(2,1,1),plot(c(i),P0),axis(16,25,0,1),title(图 1-1) subplot(2,1,2),plot(c(i),Lq),axis(16,25,0,50),title(图 1-2)3）灵敏度分析的图形三的程序： x=16.676; a=2.2; c=16:1:25 s1=1; s2=0; b=0.455; for i=1:10 p(i)=x/c(i); end for i=1:10 s1(i)=1; d(i)=1; for m=1:c(i)-1 d(i)=m.*d(i); s1(i)=x.m/d(i)+s1(i) end end for i=1:10 c1(i)=c(i)*d(i); f(i)=xc(i); s2(i)=f(i)/c1(i)/(1-p(i); s(i)=s1(i)+s2(i); p0(i)=1/s(i); x1(i)=xc(i); x2(i)=p0(i)*x1(i); n(i)=1/c1(i)/(1-p(i) P0(i)=(1-n(i)