数学建模协会第二讲(模拟)

基于matlab的计算机模拟舒兴明TelQ：117562750一、模拟是什么模拟(也称仿真)：就是用计算机程序在计算机上模仿各种实际系统的运行过程，并通过计算了解系统随时间变化的行为或特性。计算机模拟：是在已经建立的数学、逻辑模型之上，通过计算机实验，对一个系统按照一定的决策原则或作业规则，由一个状态变换为另外一个状态的行为进行描述和分析。实际问题数学、逻辑模型计算机模型数学、计算机解实际解二、为何要模拟（1）实际系统建立之前，要对系统的行为或结果进行分析研究；（2）有些真实系统做实验会影响系统运行，例如，在生产中任意改变工艺系数可能导致废品，在经济活动中随意将一个决策付诸行动会导致经济混乱；（3）在系统上做多次试验，很难保证每次试验的操作条件相同，因而对实验结果好坏很难作出正确的判断；（4）当人是系统的一部分时，他的行为往往实验结果有所影响，这时，最好进行模拟研究；（5）实验时间太长，费用太大，或者有危险，使得试验不容易进行；（6）有些系统一旦建立起来后无法复原，例如，建立大型企业，要分析社会和经济效益，不能用建立起来试试看的办法。三、哪些问题适合模拟研究（1）难以用数学公式表示的系统，或者没有求解数学模型的有效方法；（2）虽然可以用解析的方法解决问题，但是问题的分析与计算过于复杂，这时计算机仿真可能提供简单可行的求解方法；（3）希望在较短时间内观察到系统的发展全过程，以估计某参数对系统行为的影响；（4）难以在实际环境中进行试验和观察，计算机仿真是唯一的方法；（5）需要对系统或过程进行长期运行的比较，从大量方案中寻找最优。四、模拟的分类模拟是系统状态随时间而变化的动态写照，因此，通常时间是模拟的主要自变量，其它的变量为因变量。（1）按照模拟过程中因变量的变化情况，可以将模拟分为离散、连续、混合3种类型；（2）如果采用模拟计算机、采用数字计算机以及联合使用则分为模拟仿真、数字仿真以及混合仿真；（3）根据仿真变量的特征分为随机模拟仿真和模糊模拟仿真。五、随机模拟的基础基础知识：线性代数部分向量矩阵的表述；微积分的基本积分和微分表述；概率论的基本分布；基本统计知识等等。模拟（蒙特卡罗方法）基本原理：概率论的大数定律独立随机序列的和按概率收敛于这些随机序列的数学期望的和（常数），即六、matlab统计工具箱中产生随机数的函数1、概率密度函数pdf分布类型调用格式备注二项分布binoy=binopdf(x,n,p)返回参数为p，n的二项分布在x处的分布律值指数分布expy=exppdf(x,mu)返回参数为mu的指数分布在x处的密度函数值Χ2分布chi2y=chi2pdf(x,v)返回参数为v的卡方分布在x处的密度函数值F分布fy=fpdf(x,v1,v2)返回参数为v1,v2的F分布在x处密度函数值几何分布geoy=geopdf(x,p)返回参数为p的几何分布在x处的分布律值超几何分布hygey=hygepdf(x,M,k,n)返回参数为M,k,n的超几何分在x处的分布律值常见随机变量的密度函数表分布类型调用格式备注正态分布normnormpdf(x,mu,sig)返回参数为mu和sig的正态分布在x处的密度值泊松分布poisspoisspdf(x,lamda)返回参数为lamda的泊松分布在x处的分布律值T分布ttpdf(x,v)返回参数为v的t分布在x处的密度函数值连续均匀分布unifunifpdf(x,a,b)返回在[a,b]上服从均匀分布的密度函数值离散均匀分布unidunidpdf(x,N)返回1~N的均匀分布概率分布律值密度函数续注意：二项分布(bino)，几何分布(geo)，超几何分布(hyge),泊松分布(poiss)，指数分布(exp)，正态分布(norm)，卡方分布(chi2)，t分布(t)，F分布(f)，均匀连续分布(unif)，均匀离散分布(unid)。2、常见的分布函数cdf分布类型调用格式备注二项分布y=binocdf(x,n,p)返回参数为p，n的二项分布在x处的分布值指数分布y=expcdf(x,mu)返回参数为mu的指数分布在x处的分布函数值Χ2分布y=chi2cdf(x,v)返回参数为v的卡方分布在x处的分布函数值F分布y=fcdf(x,v1,v2)返回参数为v1,v2的F分布在x处分布函数值几何分布y=geocdf(x,p)返回参数为p的几何分布在x处的分布值超几何分布y=hygecdf(x,M,k,n)返回参数为M,k,n的超几何分在x处的分布值常见随机变量的分布函数数表分布类型调用格式备注正态分布normcdf(x,mu,sig)返回参数为mu和sig的正态分布在x处的分布值泊松分布poisscdf(x,lamda)返回参数为lamda的泊松分布在x处的分布值T分布tcdf(x,v)返回参数为v的t分布在x处的分布函数值连续均匀分布unifcdf(x,a,b)返回在[a,b]上服从均匀分布的分布函数值离散均匀分布unidcdf(x,N)返回1~N的均匀分布概率分布值分布函数续将分布函数的概率调用p=cdf(‘name’,x,a1,a2)调过来用，就是逆分布函数x=icdf(‘name’,p,a1,a2)。>>x1=icdf('norm',0.95,0,1)x1=1.644853626951473>>x2=icdf('t',0.95,10)x2=1.812461122811677>>x3=icdf('exp',0.95,3)x3=8.987196820661970>>x4=icdf('chi2',0.95,10)x4=18.3070380532751433、模拟随机数的产生格式1：x=random(‘name’,A1,A2,A3,m,n)产生参数为A1，A2，A3的name分布的m行n列的矩阵。>>x=random('norm',0,1,3,4)x=0.53770.8622-0.43362.76941.83390.31880.3426-1.3499-2.2588-1.30773.57843.0349>>y=random('unif',1.2,2.2,3,4)y=2.15721.34191.99221.23571.68541.62182.15952.04912.00032.11571.85572.1340产生3行4列的参数为0和1的正态分布随机数产生在[1.2,2.2]上服从均匀分布的随机数格式2：x=namernd(A1,A2,m,n)产生m行，n列的服从参数为A1，A2的name分布的随机数。>>x=trnd(3,4)x=0.6880-1.9018-0.6416-1.14291.0967-1.64313.3100-2.16210.9627-1.1476-0.7321-1.18490.07281.10551.62900.6073>>y=trnd(3,3,4)y=-4.3874-1.88741.0409-1.5729-1.08610.8758-2.57105.65671.75710.1385-0.25441.2256产生自由度为3的服从t分布的随机数（4阶方阵）。产生3行4列自由度为3的服从t分布的随机数。七、随机样本的简单描述性统计量1、中心趋势（位置）度量样本的几何平均m=geomean(X)样本的调和平均m=harmmean(X)样本的算术平均m=mean(X)样本的中位数m=median(X)样本的剔出极端数据的平均m=trimmean(X,percent)>>X=random('unif',1,8,10,1);>>[geomean(X),harmmean(X),mean(X),median(X),trimmean(X,0.2)]ans=4.59033.74005.24016.06785.24012、分散度量样本四分位数间距m=iqr(X)样本数据的平均绝对偏差m=mad(X)即mean(abs(X-mean(X)))样本极差m=range(X)样本方差m=var(X)样本标准差m=std(X)>>X=random('unid',20,10,1);>>m=[iqr(X),mad(X),range(X),var(X),std(X)]m=9.00005.160018.000035.73335.9777八、随机模拟计算积分例1如下图所示，在正方形内有1/4单位圆。向正方形内投小石头，假设每次都能够投进正方形内且等可能落在正方形内任何一点。问，小石头落在1/4单位圆内(包含边界)的概率多大？01x1y分析：假设投入正方形内的石头有n块，有k块落入了1/4单位圆内。P为小石头落入1/4单位圆内的概率。那么根据Bernoull(伯努利)大数定律，有即，当实验次数n充分大时，频率和概率之差小于任意数ε的概率趋于1。而另外一方面由几何概型有在实际操作中，实验次数n不可能趋于无穷大，所以有(数学模型)对于估计只有不断重复做实验，这种试验可以具体去操作，(均匀投石块，然后数数，这样需要较高成本)。也可以让计算机去重复试验，但是需要将数学模型转化为计算机模拟模型(让计算机完成均匀投石块，自动计数，也需要成本，但成本比起前者，几乎可以忽略。)。用计算机模拟投石块过程和步骤如下：1、自动生成随机点[0,1]x[0,1]，模拟石块在正方形内的任意位置，用(xi,yi)表示，共n个点；2、判别(xi,yi)是否满足xi2+yi2≤1，即判别石块是否落在1/4单位圆内，共k个点满足；3、整理、统计模拟结果，用4k/n估计π。符号设置：符号说明n试验次数（投石块次数）k石块落入1/4圆内次数(xi,yi)第i块石头落点坐标pai圆周率计算机模拟的流程图如下初始化：k=0,i=1i>n是Pai=4k/n否在[0,1]x[0,1]上产生随机点(xi,yi)r=xi^2+yi^2r<=1是k=k+1i=i+1否设计一个m文件，专门针对给定n，模拟计算pai的估计值functionpai=fpai(n)ifn<=0error('nhastobeapositiveinteger');endk=0;i=1;whilei<=nx=random('unif',0,1,1,2);r=sum(x.^2);ifr<=1k=k+1;endi=i+1;endpai=4*k/n;模拟的稳定性描述>>clear>>tic>>n=50:10:5000;>>fork=1:length(n)pai(k)=fpai(n(k));end>>plot(n,pai,'*')TocElapsedtimeis223.847533seconds.>>mean(pai)ans=3.1405！！值得注意的是，每次模拟即或程序相同，不同人不同计算时间，结果一般不同(因为随机)。functionpai=fpai1(n)ifn<=0error('nhastobeapositiveinteger');endk=0;x=random('unif',0,1,n,2);fork1=1:nx1=x(k1,:);r=sum(x1.^2).^0.5;ifr<=1k=k+1;endendpai=4*k/n;把while语句改为for语句，同时不用每次都调用random函数，m文件如下>>clearticn=50:10:5000;fork=1:length(n)pai(k)=fpai1(n(k));endplot(n,pai,'*')tocElapsedtimeis7.372741seconds.时间大大缩短了！！计算定积分分析：若随机变量X的概率分布密度是P(x)(a≤x≤b)，则随机变量Y=f(X)的数学期望为当X是[a,b]上的均匀分布时，有将p(x)代入上面的期望算子，有例2而另外一方面，根据大数定理，设y1,y2,…,yn为来自总体Y的一组样本，有其中所以有设x1,x2,…,xn为来自服从[a,b]上均匀分布总体X的一组样本。且yi=f(xi)，i=1,2,…,n。那么计算的步骤如下：1、在[2，3]上抽样x1,x2,…,xn；2、计算3、整理估计定积分；为了对比模拟结果，先把这个积分求出来：>>int('(1+x^2)^0.5',2,3)ans=2.6947540047794075420843777595616然后根据上面的模拟步骤设计一个m文件functionI=jifen1(n)x=random('unif',2,3,1,n);y=(1+x.^2).^0.5;I=(3-2)/n*sum(y);>>clear>>tic;n=50:10:5000;fork=1:length(n)I(k)=jifen1(n(k));endplot(n,I,'*');tocElapsedtimeis0.334083seconds.例3求在空间坐标系下，两个圆柱体相交部分的体积。两个圆柱体的柱面方程分别为根据二重积分的知识，相交部分的体积为xyz=16/3计算的模拟步骤如下（1）在空间[0,1]x[0,1]x[0,1]上产生均匀随机点(xi,yi,zi)，i=1,2,…,n；（2）判断是否成立，若有k个；（3）用8k/n来估计V的值。把以上步骤编写成m文件functionV=erchjfen(n)x=random('unif',0,1,n,3);k=0;fork1=1:nr1=x(k1,1)^2+x(k1,2)^2;r2=x(k1,1)^2+x(k1,3)^2;if(r1<=1)&(r2<=1)k=k+1;endendV=8*k/n;>>clear>>tic>>n=20:20:4000;>>fork=1:length(n)v(k)=erchjfen(n(k));end;plot(n,v,'*');tocElapsedtimeis72.712497seconds.九、评估事件发生的概率有些随机变量，由于不容易找到（或者不存在）准确的概率分布律或者密度函数，不能算出此类随机变量表达的随机的事件的概率，需要用到模拟。其中X1服从均匀分布U[2,5]，X2服从指数分布exp(3)，X3和X4分别服从正态分布N(3,2)和N(1,1)。例4求下面事件表达的概率1、在总体X=(X1,X2,X3,X4)中抽样，得到样本(x1k,x2k,x3k,x4k)，k=1,2,…,N；2、判断随机样本是否满足不等式组，假设满足不等式组的有n个；3、用n/N估计p。（大数定律作保障，或者点估计原理）步骤：按照上述步骤编写m文件functionp=gailv1(n)x1=random('unif',2,5,n,1);x2=random('exp',3,n,1);x3=random('norm',3,2,n,1);x4=random('norm',1,1,n,1);k1=0;fork=1:nr1=x1(k)+x2(k)^2;r2=x3(k)+x4(k)^2;if(r1>=3)&(r2<=9)k1=k1+1;endendp=k1/n;>>n=20:20:4000;>>tic>>fork=1:length(n)p(k)=gailv1(n(k));end;plot(n,p,'*');tocElapsedtimeis38.580145seconds.>>p(length(n))ans=0.8375十、求随机变量的临界值（上下分为数）求使下式成立的最大f值步骤：1、产生N个独立随机序列(X1k,X2k,X3k),k=1,2,…,N；2、计算fk=f(X1k,X2k,X3k),k=1,2,…,N;3、从从大到小的序列f1,f2,…,fN中取f0.8N作为f的估计值，或从小到大的f0.2N作为f的估计值。例5其中，X1服从均匀分布U[1,3]，X2服从指数分布exp(1)，X3服从正态分布N(2,1)。根据上面的流程，编写m文件如下functionf=linjie(n)x1=random('unif',1,3,n,1);x2=random('exp',1,n,1);x3=random('norm',2,1,n,1);fork=1:ny(k)=x1(k)+x2(k)^2+x3(k)^3;endz=sort(y);k1=round(0.2*n);f=z(k1);>>clearn=50:50:10000;ticfork=1:length(n)f(k)=linjie(n(k));end;plot(n,f,'*');tocElapsedtimeis11.335281seconds.十一、报童问题有个报童每天从报社批发报纸，每份报纸1元，然后再以零售价1.5元每份售出，如果卖不出，则退回报社，退回价格为每份1.2元。某个社区某个月记录在案的每天报纸的需求量如下表需求量50454038353028天数2367732频率0.06670.10.20.23330.23330.10.0667问：报童每天该批发多少分报纸才是最佳策略？合理假设（1）报童兜售报纸的社区比较稳定，即每个月的报纸的需求量的分布抽样差不多；（2）把上表给的结果看成社区每天报纸的需求分布律。例6符号假设a订购价（元/份）b零售价（元/份）c退回价（元/份）r每天报纸的需求量P(r)每天报纸的需求量的抽样分布律x报童报纸的批发量（份）y报童每天的获利（元）建立模型每天报纸的需求量r是个随机变量，报童每天的盈利由r和批发量x决定，也是随机变量，即目标1每天盈利期往最大maxE(y)数学模型1（数学期望模型）r~P(r)模型1求解将报童每天盈利函数编写成m文件functionEy=baotong(x)globalabcPrifb<a|c<aerror('bandcmustbelargerthana');endn=size(Pr,2);r=Pr(1,:);p=Pr(2,:);fork=1:nifx<=r(k)y(k)=(b-a)*x;elsey(k)=(b-a)*r(k)-(c-a)*(x-r(k));endendEy=sum(y.*p);在命令窗口调用上述函数，计算得>>clear>>globalabcPr>>a=1;b=1.5;c=1.2;>>Pr=[50454038353028;0.06670.10.20.23330.23330.10.0667];>>x=1:1:60;>>fork=1:length(x)Ey(k)=baotong(x(k));end>>plot(x,Ey,'*')>>Ey1=max(Ey);>>m=find(Ey==Ey1)m=40最佳批发为40份报纸。数学模型2（最低盈利机会最大模型）目标2r~P(r)每天至少赚16元的概率最大（机会最大）将每天赚16元的机会计算程序编制为m文件：functionp=baotongjh(x,ybar)globalabcPrifb<a|c<aerror('bandcmustbelargerthana');endn=size(Pr,2);r=Pr(1,:);pr=Pr(2,:);fork=1:nifx<=r(k)Ey(k)=(b-a)*x;elseEy(k)=(b-a)*r(k)-(c-a)*(x-r(k));endendpy=Ey>ybar;p=sum(py.*pr);>>clear>>globalabcPr>>a=1;b=1.5;c=1.2;>>Pr=[50454038353028;0.06670.10.20.23330.23330.10.0667];>>ybar=15;>>x=1:60;>>fork=1:60p(k)=baotongjh(x(k),ybar);end>>plot(x,p,'*')由此可见，当报纸批发量在30~45之间时，报童每天可以赚15元钱的概率最大有83%。十二、利用模拟法求解非线性规划例7用随机模拟方法求解给定的规划问题求解此问题的流程图如下初始化n=0,A=[]n>1000YN在[0,10]上产生均匀随机数x1，x2=(10-x1)/3，y=x1+2x2y>2NYn=n+1,A=[A;[x1,x2]]z=0,k=1:1000,xx=[0,0]x1=A(k,1),x2=A(k,2)z1=f(x1,x2)z1>zYz=z1Nk=k+1xx=[x1,x2]依照上述流程编写一个求此非线性规划最大值的m文件：function[z,xx]=nlpf(N)n=0;A=[];whilen<=Nx1=random('unif',0,10);x2=(10-x1)/3;y=x1+2*x2;ify>2n=n+1;A=[A;[x1,x2]];endendz=0;xx=[0,0];fork=1:Nx11=A(k,1);x22=A(k,2);z1=-2*x11^2-x22^2+x11*x22+8*x11+3*x22;ifz1>zz=z1;xx=[x11,x22];endend>>clear>>A=[];>>N=20:5:400;>>fork=1:length(N)[z,xx]=nlpf(N(k));A=[A;z,xx];end>>plot(N,A(:,1),'*')十三、一个排队系统的模拟例7一个排队系统的仿真某店只有一个收款台，顾客到达收款台的时间间隔服从均值为4.5的负指数分布，每个顾客的服务时间服从均值为3.2，标准差为0.6的正态分布，对100位顾客去收款台缴款的排队过程进行仿真。分析：排队系统涉及的参考变量为：顾客的等待时间Wq，顾客的逗留时间W,排队长Lq,队长L，服务台的空闲时间I，服务台的繁忙时间B，排队系统的状态概率P。最后需要统计出有关指标的平均值均方差。合理假设1、进入系统的顾客没有人因为不愿意等而离开系统；2、先到先服务；模拟步骤：1、产生初始事件表，并从初始事件表中找出最靠前的事件；2、启动仿真时钟；3、判别当前事件是哪一个事件，如果是顾客到达事件，则启动顾客到达子程序；如果是服务事件，则启动服务结束子程序；4、重复2-3，直到仿真完毕；5、整理、分析并输出结果。顾客到达的子程序1、产生下一个顾客到来的时刻，记入事件表；2、判别服务台是否忙s=1?，若是，则顾客则排队长L=L+1;否则，服务台空闲,L=0;3、产生一个服务结束时刻，记入事件表；4、统计所需的数据；服务结束子程序1、结束服务，产生结束时间，记入事件表；2、判别是否有顾客排队，L=0?，若L=0，则服务台空闲；否则，排队长L=L-1;3、产生服务结束时刻，记入事件表；4、统计所需的数据；产生初始事件表从事件表中找出最靠前的事件仿真时钟步进是哪一类事件顾客到达服务结束产生下一个顾客到达纪录，记入事件表结束服务，产生结束时间，记入事件表出纳是否忙？出纳空闲，s=0顾客排队，Lq=Lq+1否是产生服务结束时刻，记入事件表是否顾客排队队长减1，Lq-1出纳空闲，s=0否产生服务结束时刻，记入事件表是统计所需数据仿真是否完毕是结束程序否按照流程，编写模拟程序function[ndaoda,nlikai,t,LLq,LL,BI,tdaoda,tfuwu]=paiduisimu(n,a,b,c)t=0；LL=[];LLq=[];BI=[];tdaoda=[];tfuwu=[];ndaoda=0;nlikai=0;L=0;Lq=0;t1=exprnd(a);t2=+inf;s=0;whilenlikai<n;ift1<t2t=t1;r1=exprnd(a);t1=t+r1;ifs==1L=L+1;Lq=max(L-1,0);elseL=1;s=1;Lq=max(0,L-1);r2=normrnd(b,c);t2=t+r2;endLL=[LL,L];LLq=[LLq,Lq];tdaoda=[tdaoda,t];

ndaoda=ndaoda+1;nlikai=nlikai;BI=[BI,s];elseift1>t2t=t2;tfuwu=[tfuwu,t];ifLq>0s=1;L=L-1;Lq=Lq-1;r2=normrnd(b,c);t2=t+r2;elses=0;L=0;Ls=max(L-1,0);t2=+inf;endLL=[LL,L];LLq=[LLq,Lq];BI=[BI,s];nlikai=nlikai+1;ndaoda=ndaoda;

elseift1==t2t=t1;tdaoda=[tdaoda,t];tfuwu=[tfuwu,t];ndaoda=ndaoda+1;nlikai=nlikai+1;L=L;Lq=Lq;s=s;LL=[LL,L];LLq=[LLq,Lq];BI=[BI,s];t1=t+exprnd(a);t2=t+normrnd(b,c);endend>>[ndaoda,nlikai,t,LLq,LL,BI,tdaoda,tfuwu]=paiduisimu(n,a,b,c);>>[mean(LLq),mean(LL),mean(BI)]ans=0.96021.76620.8060由此可见，这个服务窗口一直处以繁忙状态（空闲率)的概率为1-0.8060;平均排队长（mean(L)=0.9602人，不用再增加一个窗口来缓解压力。十四、航班延误问题某航空公司经营A，B，C三个城市的航线，这些航线每天班次起飞与到达时间如下表所示。设飞机在机场停留的损失费大致与停留时间成正比，又每架飞机从降落到下班起飞至少需2小时准备时间，试决定一个使停留费用损失为最小的分派飞行方案。另外，根据抽样资料，三个城市A，B，C之间的航班延误时间分别服从负指数分布：城市A到城市B的航班延误时间服从均值为1的指数分布；城市A到城市C的航班延误时间服从均值为1.5的指数分布；城市B到城市C的航班延误时间服从均值为0.8的指数分布；城市B到城市A的航班延误时间服从均值为1.2的指数分布；城市C到城市A的航班延误时间服从均值为0.8的指数分布；城市C到城市B的航班延误时间服从均值为0.5的指数分布；航班号

起飞城市

起飞时间

到达城市

到达时间101A9：00B

12：00102

A10：00

B

13：00103

A

22：00

C

2：00（次日)104

B4：00A7：00

105

B11：00A

14：00

106

C7：00

A

11：00

107

B13：00

C

18：00

108C

15：00B

20：00各个航班往返城市及其时间如下表：合理假设（1）城市A，B，C的到达与起飞航班数目一样，每个城市都用到达航班的飞机作为起飞航班的飞机;（2）没有多余的飞机可用；（3）只算一天内的延误时间；（3）飞机到达后2个小时任何时刻都可以起飞；（4）先到先服务原则；（5）概率小于0.0001为小概率事件，不考虑。符号设置xijxij=1,城市A的第i到达作为第j离开，否则xij=0yijyij=1,城市B的第i到达作为第j离开，否则yij=0zijzij=1,城市C的第i到达作为第j离开，否则zij=0t1ijA城市第i到达作为第j离开的等待时间t2ijB城市第i到达作为第j离开的等待时间t3ijC城市第i到达作为第j离开的等待时间TEF第i城市到第j城市的延误时间E=A,B,C,F=A,B,C,E≠F同时根据假设（5），所有航班延误超过某个时刻后是小概率事件，认为不发生。1、可行解的产生方法k<=nNstop初始化：y=[],k=1Y产生随机数a1~U(0,1)a1<=1/3x=[100]YNa1>2/3Yx=[001]Nx=[010]产生随机数a2~U(0,1)a2<=0.5确定x的0所在列第二行1的位置确定剩1行1列1的位置y=[y;x]k=k+1functiony=hbzhipai3(n)y=[];k=1;whilek<=na1=random('unif',0,1);ifa1<=1/3x=[1,0,0];a2=random('unif',0,1);ifa2<=0.5x=[x;[0,1,0]];x=[x;[001]];elsex=[x;[001]];x=[x;[010]];endelseifa1>2/3x=[001];a2=random('unif',0,1);ifa2<=0.5x=[x;[100]];x=[x;[010]];elsex=[x;[010]];x=[x;[100]];endelsex=[010];a2=random('unif',0,1);ifa2<=0.5x=[x;[100]];x=[x;[001]];elsex=[x;[001]];x=[x;[100]];endendy=[y;x];k=k+1;end但是，根据这个问题的特殊性，城市A的航班安排只有厦门6个可行解：x1=[100;010;001];x2=[100;001;010];x3=[010;001;100];x4=[010;100;001];x5=[001;010;100];x6=[001;100;010];2、模拟在给定延误情况下的等待时间航班号\航班号101(9)102(10)103(22)B：104(7)2315B：105(14)19208C：106(11)222311城市A的航班间的不延误等待时间t1ij2.1航班延误时A城市的各个等待时间表由于B和C到A城市的航班延误时间服从的分布分别为当TBA和TCA取定时，航班飞机的等待时间就跟着确定，如下流程初始化：tba1~e(1.2)tba1<=1Yt1(104,102)=3-tba1t1(104,103)=15-tba1t1(104,101)=25-tba1Ntba1>1,tba1<=13Yt1(104,103)=15-tba1t1(104,101)=24-tba1t1(104,102)=25-tba1Nt1(104,101)=26-tba1t1(104,102)=17-tba1t1(104,103)=39-tba1初始化：tba2~e(1.2)tba2<=6Yt1(105,102)=20-tba2t1(105,103)=8-tba2t1(105,101)=19-tba2Ntba2>6,tba2<=17Yt1(105,103)=24-tba2t1(105,101)=19-tba2t1(105,102)=20-tba2Nt1(105,101)=43-tba2t1(105,102)=20-tba2t1(105,103)=32-tba2tba2>17,tba2<=18YNt1(105,101)=43-tba2t1(105,102)=44-tba2t1(105,103)=32-tba2初始化：tba3~e(0.8)tba3<=9Yt1(106,102)=23-tba3t1(106,103)=11-tba3t1(106,101)=22-tba3Ntba3>9,tba3<=20Yt1(106,103)=33-tba3t1(106,101)=22-tba3t1(106,102)=23-tba3Nt1(106,101)=46-tba3t1(106,102)=23-tba3t1(106,103)=35-tba3tba3>20,tba3<=21YNt1(106,103)=35-tba3t1(106,101)=46-tba3t1(106,102)=47-tba3按照这三个流程，编写M文件elseif(tba2>17)&(tba2<=18)t12=[432032]-tba2;elset12=[434432]-tba2;endtba3=random('exp',0.8);iftba3<=9t13=[222311]-tba3;elseif(tba3>9)&(tba3<=20)t13=[222333]-tba3;elseif(tba3>20)&(tba3<=21)t13=[462335]-tba3;elset13=[464735]-tba3;endtr1=[t11;t12;t13];t1=[t1;tr1];k=k+1;endfunctiont1=hbwaita(m)t1=[];k=1;whilek<=mtba1=random('exp',1.2);iftba1<=1t11=[25315]-tba1;elseif(tba1>1)&(tba1<=13)t11=[242515]-tba1;elset11=[262739]-tba1;endtba2=random('exp',1.2);iftba2<=6t12=[19208]-tba2;elseif(tba2>6)&(tba2<=17)t12=[192024]-tba2;2.2寻找城市A的最小费用航班安排模型按照普通指派问题，数学模型为步骤如下：（1）针对每一个可行解x；（2）产生m个航班延误时间t1；（3）计算可行解x在各种延误时间t1的总的安排时间，并把这m个时间平均，记录其为x对应的目标；（4）比较每个解所对应的目标，选最小的最为参考最优解。2.3寻找最优解按照上述流程编写m文件：function[tt]=zuiyoua1(x,m)tt1=[];fork1=1:mt1=hbwaita(1);yt=x.*t1;tt11=sum(sum(yt));tt1=[tt1;tt11];endtt=mean(tt1);在命令窗口计算，结果显示如下>clearm=10000;x1=[100;010;001];x2=[100;001;010];x3=[010;001;100];x4=[010;100;001];x5=[001;010;100];x6=[001;100;010];tt1=zuiyoua1(x1,m);tt2=zuiyoua1(x2,m);tt3=zuiyoua1(x3,m);tt4=zuiyoua1(x4,m);tt5=zuiyoua1(x5,m);tt6=zuiyoua1(x6,m);[tt1,tt2,tt3,tt4,tt5,tt6]ans=52.366652.480439.511739.058953.812553.8071计算结果显示，可行解x4比较好！X3其实也不错。对于城市B和城市C，可以用相同的方法模拟计算结果。DVD在线租赁

2005年国赛B题舒兴明TeQ117562750原问题考虑如下的在线DVD租赁问题。顾客缴纳一定数量的月费成为会员，订购DVD租赁服务。会员对哪些DVD有兴趣，只要在线提交订单，网站就会通过快递的方式尽可能满足要求。会员提交的订单包括多张DVD，这些DVD是基于其偏爱程度排序的。网站会根据手头现有的DVD数量和会员的订单进行分发。每个会员每个月租赁次数不得超过2次，每次获得3张DVD。会员看完3张DVD之后，只需要将DVD放进网站提供的信封里寄回（邮费由网站承担），就可以继续下次租赁。请考虑以下问题：（1）网站正准备购买一些新的DVD，通过问卷调查1000个会员，得到了愿意观看这些DVD的人数（表1给出了其中5种DVD的数据）。此外，历史数据显示，60%的会员每月租赁DVD两次，而另外的40%只租一次。假设网站现有10万个会员，对表1中的每种DVD来说，应该至少准备多少张，才能保证希望看到该DVD的会员中至少50%在一个月内能够看到该DVD？如果要求保证在三个月内至少95%的会员能够看到该DVD呢？（2）表2中列出了网站手上100种DVD的现有张数和当前需要处理的1000位会员的在线订单（表2的数据格式示例如下表2，具体数据请从/mcm05/problems2005c.asp下载），如何对这些DVD进行分配，才能使会员获得最大的满意度？请具体列出前30位会员（即C0001~C0030）分别获得哪些DVD。（3）继续考虑表2，并假设表2中DVD的现有数量全部为0。如果你是网站经营管理人员，你如何决定每种DVD的购买量，以及如何对这些DVD进行分配，才能使一个月内95%的会员得到他想看的DVD，并且满意度最大？（4）如果你是网站经营管理人员，你觉得在DVD的需求预测、购买和分配中还有哪些重要问题值得研究？请明确提出你的问题，并尝试建立相应的数学模型。表1对1000个会员调查的部分结果DVD名称DVD1DVD2DVD3DVD4DVD5愿意观看的人数200100502510表2现有DVD张数和当前需要处理的会员的在线订单DVD编号D001D002D003D004…DVD现有数量10401520…会员在线订单C00016000…C00020000…C00030003…C00040000…………………注：D001~D100表示100种DVD,C0001~C1000表示1000个会员,会员的在线订单用数字1,2,…表示，数字越小表示会员的偏爱程度越高，数字0表示对应的DVD当前不在会员的在线订单中。合理假设（1）抽样调查的1000人具有代表性，即频率可以代替概率；（2）DVD租赁业务按规则进行，客户在规定时间内归还DVD；（3）每次租赁都是3张DVD；（4）每个会员每月的租赁次数与对DVD的偏好是独立的；（5）各个会员之间的选择和偏好也是相互独立的；符号设置Xi会员i每月租赁的次数，Xi=1,2,i=1,2,…,100000;Yij会员i对第j类DVD的选择，Yij=3,0,j=1,2,3,4,5;i=1,2,…,100000Rj第j类DVD的最少定购张数，j=1,2,3,4,5;问题（1）的模型根据条件和假设（5）会员i的每月定购次数X1,X2,…,X100000是服从同分布的随机变量，且分布规律为Xi21P(Xi)0.60.4且X1,…,X100000相互独立。在100000个会员中抽取1000个会员，样本容量n=1000足够大，会员对DVD的选择频率可以直接作为选择概率的点估计