计算机仿真-2014数学建模

计算机仿真,计算机仿真的基本概念随机数的产生时间步长法事件步长法MonteCarlo方法建模实例,一、计算机仿真的基本概念,仿真,仿真就是将所研究的对象用其它手段加以模仿的一种活动。,实物仿真,非实物仿真,如曹冲称象、军事演习、飞行器风洞试验、核爆炸试验等，属于实物仿真的例子。,1.计算机仿真的基本概念,计算机仿真是一种非实物仿真方法计算机仿真通过建立数学模型、编制计算机程序实现对真实系统的模拟，对系统的结构和行为进行动态演示,从而了解系统随时间变化的行为或特性，以评价或预测系统的行为效果。它是解决较复杂的实际问题的一条有效途径。,计算机仿真反映出新的科学技术的时代特征,它的应用为各个领域带来新气象和成果。应用的领域有：航空管理,公交车的调度,飞机设计,动画设计,三峡的安全、生态,道路的修建,医疗保险,国债的发行,家居装修,炼钢的温度估计,发电厂的操作训练,飞行员训练，鼠疫的检测和预报。,黑死病菌寄生于老鼠身上,是由跳蚤传染给人类，又叫鼠疫。病菌随着跳蚤叮咬进入人体，约2-5天的潜伏期之后，患者鼠蹊部及其他淋巴结开始红肿、疼痛，随之开始发高烧、疲倦、皮肤变黑，故被称为黑死病。死亡率高达60-90%,鼠疫的传播,最原始的生化武器,1346年，蒙古大将去攻打黑海边富庶的卡法城，久攻不下，这时蒙古军中发生鼠疫，士兵死亡无数，眼看就要无法而退了，这时蒙古将军想出一个方法，把死亡士兵的尸体用弹弩投入城中，迫使城中流行鼠疫，城门自然不攻而破。在城破时，一位意大利热内亚的富商，帶着妻小和金银珠宝乘船逃了出來。他在地中海各国漂流很久，沒有国家敢收留他们，大家都害怕鼠疫的传染。最后回到家乡热內亚，他把所有的财富全部推在甲板上，对着守城的人说：“我离开卡法城已经六個月了，我若感染鼠疫早就死了，但我并沒有死，可见我并沒有瘟疫。假如你让我进城，甲板上的珠宝就是你们的。”我们现在知道鼠疫是由老鼠身上的跳蚤传染的，通常老鼠躲在船底污秽处，人们不易察觉。热内亚人打开城门让这艘船进來后，鼠疫就从热内亚传播开来，传遍整个欧洲，包括北方的斯堪尼亚半岛都无法幸免。每天黃昏时有人推著独轮车，手里摇着铃喊着：“bringoutthedead,bringoutthedead.”（把尸体拿出來，把尸体拿出來）家家户户就把尸体搬出来丟到城外焚烧，说死尸如山是一点都不为过。鼠疫肆虐欧洲一百多年，使得三分之一的人口死亡。,欧洲鼠疫流行时死亡无数,鼠疫期间贵族纷纷弃城逃亡,疫如此让人恐怖，那么有没有什么好的预测方式呢？计算机仿真就是一个很好的预测方法。研究发现，鼠疫是由老鼠身上一种特殊的跳蚤传播的。跳蚤的多少决定是否发生鼠疫。跳蚤老鼠水草我们可以用计算机根据一个地区的气候模拟出当年此地水草的情况就可以预测出是否有鼠疫要发生。,三峡水库总库容393亿立方米，总装机容量1820万千瓦，将是世界上最大的水电站。但是三峡的安全问题是一个很重要的问题，我们不可能等到建好后再看它的安全性，用计算机仿真就可以很好的解决这一问题。,长江三峡工程,飞机设计中有一个重要环节：风洞试验。实际的风洞试验费用巨大。使用计算机仿真进行模拟风洞试验，使费用大大降低。,飞机设计,2.为什么要进行计算机仿真,便于重复进行试验，便于控制参数，时间短，代价小。可以在真实系统建立起来之前，预测其行为效果，从而可以从不同结构或不同参数的模型的结果比较之中，选择最佳模型。对于缺少解析表示的系统，或虽有解析表示但无法精确求解的系统，可以通过仿真获得系统运行的数值结果。对于随机性系统，可以通过大量的重复试验，获得其平均意义上的特性指标。,3.适用计算机仿真解决的问题,难以用数学公式表示的系统，或者没有建立和求解数学模型的有效方法虽然可以用解析的方法解决问题，但数学的分析与计算过于复杂，这时计算机仿真可能提供简单可行的求解方法希望能在较短的时间内观察到系统发展的全过程，以估计某些参数对系统行为的影响难以在实际环境中进行实验和观察时，计算机仿真是唯一可行的方法，例如太空飞行的研究需要对系统或过程进行长期运行比较，从大量方案中寻找最优方案,4.仿真常用术语,系统,一些具有特定功能相互之间以一定的规律联系着的物体所组成的总体,为了限制所研究问题涉及的范围,用系统边界把所研究的系统与影响系统的环境区分开来.,系统的对象、系统的组成元素都可以称为实体,系统边界,实体,属性反映实体的某些性质，它可以是文字型、数字型或逻辑型,系统的状态是指在某一时刻实体及其属性值的集合,导致系统状态变化的一个过程为活动活动反映了系统变化的规律,属性,状态,活动,活动是指一段过程，即在一段时间内发生的情况事件是指一个时间的情况，系统发生变化的瞬间就发生了事件,为了使仿真程序能如实地模拟实际系统的变化，在某些离散事件的仿真中，采用事件表的形式进行调度事件表一般是一个有序的记录列，每个记录包括事件发生时间、事件类型等一些内容,事件,事件表,研究系统一般是为认识其状态随时间变化的规律，所以需要一个仿真时间变量对连续系统仿真时，常在均匀时间点上展现其状态值，这样，仿真时钟的步进是一个常数对离散系统仿真时，只有在事件发生时，系统的状态才会发生变化，才有必要展现出系统的状态，此时仿真时钟的步进根据事件发生的时刻变化。,仿真时钟,5.仿真研究的步骤,明确问题和提出总体方案把被仿真系统的内容表达清楚；弄清仿真的目的、系统的边界；确定问题的目标函数和可控变量；找出系统的实体、属性和活动等。,系统分析,建立模型选择合适的仿真方法（如时间步长法、事件表法等）；确定系统的初始状态；设计整个系统的仿真流程图。收集数据编写程序、程序验证模型确认,模型构造,仿真研究的步骤,运行确定具体的运行方案,如初始条件、参数、步长、重复次数等,然后输入数据,运行程序。改进将得出的仿真结果与实际系统比较,进一步分析和改进模型,直到符合实际系统的要求及精度为止。,模型的运行与改进,结果输出设计出结构清晰的仿真结果输出。包括提供文件的清单,记录重要的中间结果等。输出格式要有利于用户了解整个仿真过程,分析和使用仿真结果.,设计格式输出仿真结果,仿真研究步骤,系统分析,模型构造,模型运行,输出结果,计算机仿真举例,例1.1混凝土搅拌中心的位置有K1(x1,y1),K2(x2,y2),.,Kn(xn,yn)共n个工地，各需混凝土Q1,Q2,.,Qn吨，混凝土每吨公里的运费为C元。如何确定混凝土搅拌中心的位置K0(x0,y0)使得费用最少？,计算机仿真举例,方法1：数学计算记第K个工地的位置为(xk,yk),中心的位置为(x0,y0)，则目标函数为,计算机仿真举例,方法2：计算机的模拟。,计算机仿真举例,方法3：工程模拟物理、工程仿真的特点：形象、直观、便于类比；但由于是实物仿真,所以费用大,速度慢。,计算机仿真举例,例1.2：追击问题我辑私舰雷达发现距c里处有一艘走私船正以匀速度a沿直线行驶,辑私舰立即以最大的速度b追赶,若用雷达进行跟踪,保持船的瞬时速度方向始终指向走私船,试求辑私舰追逐路线和追上的时间,并且给出缉私舰和走私船的路线图。,选取坐标系如图:,假设缉私艇初始位置在点(c,0),走私船初始位置在点(0,0),走私船的行驶方向为y方向.,设缉私艇为动点D,走私船为动点R.在时刻t,缉私艇的位置是D(x,y),走私船的位置是R(0,at).,取时间间隔(步长)为,则在时刻t+,D的位置是,(*),计算机仿真举例,算法:,赋初值:初始时刻t0,时间步长,速度a,b,初始位置c,循环:,由公式(*)计算动点D在各时刻的坐标D,计算动点R在各时刻的坐标R.,终止:当D,R的距离小于给定值时终止,计算机仿真举例,将上述算法中求出的点D用直线段连接起来，就得到追赶曲线。循环终止时的即为追赶时间。可以看出，步长选取得越小，所求的曲线越精确。,计算机仿真举例,例1.3理发店的服务过程仿真一个理发店有两位服务员A和B，顾客随机地到达该理发店,每分钟有一个顾客到达和没有顾客到达的概率均是1/2,其中60%的顾客理发仅用5分钟,另外40%的顾客用8分钟.试对前10分钟的情况进行仿真。,解：假设开始无顾客，顾客到达、服务开始和结束都在每分钟开始时进行，产生顾客：抛硬币，正面-有顾客到来；反面-无顾客到来时间长短：摸球，3白2黑，白球，5分钟，黑球，8分钟。仿真过程:.,计算机仿真举例,计算机仿真举例,例1.7射击命中率在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点经过长期观察发现，我方指挥所对敌方目标的指示有50是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部消灭敌人现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值。,计算机仿真举例,分析：这是一个概率问题，可以通过理论计算得到相应的概率和期望值.但这样只能给出作战行动的最终静态结果，而显示不出作战行动的动态过程.为了能显示我方20次射击的过程，现采用模拟的方式。,计算机仿真举例,问题分析:需要模拟出以下两件事：(1)指挥所对目标的指示正确与否对目标的指示有两种结果，出现的概率都是1/2因此，用投掷一枚硬币的方式予以确定，当硬币出现正面时为指示正确，反之为不正确(2)当指示正确时，我方火力的射击结果模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6),计算机仿真举例,毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6)这时可用投掷骰子的方法来确定：出现、点：则认为没击中敌人；出现、点：则认为击毁敌人一门火炮；出现点：则认为击毁敌人两门火炮,计算机仿真举例,2.符号假设,i：要模拟的打击次数；k1：没击中敌人火炮的射击总数；k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数E：有效射击比率；E1：20次射击平均每次毁伤敌人的火炮数,计算机仿真举例,3.模拟框图,计算机仿真举例,4.模拟结果,计算机仿真举例,计算机仿真举例,5.理论计算,计算机仿真举例,计算机仿真举例,6.结果比较,模拟结果与理论计算虽然结果不完全一致，可这反映了事件发生的随机性。只要多次实验求平均值，模拟值就会很接近理论值。,二、随机数的产生,现实世界充满不确定性，我们所研究的现实对象往往难以摆脱随机因素的影响要使我们的数学模型能够较真实地刻画实际对象，必须面对这个现实概率论是用数学的思想和方法处理和研究随机现象的一个有效的工具但是它还难以用来处理复杂系统中的随机性这里介绍使用计算机来模拟随机现象的方法，它经常应用于复杂系统的动态仿真的研究当中是处理复杂系统中随机性的计算机模型也是使用计算机研究和解决实际问题的一条重要途径,1均匀分布的随机数及其产生,对随机现象进行模拟，实质上是要给出随机变量的模拟也就是说利用计算机随机地产生一系列数值，它们的出现服从一定的概率分布,则称这些数值为随机数。最常用的是在(0,1)区间内均匀分布的随机数，也就是我们得到的这组数值可以看作是(0,l)区间内均匀分布的随机变量的一组独立的样本值以后我们将指出其它分布的随机数可利用均匀分布的随机数产生,记XU(0,1),0,1)区间上均匀分布的随机变量X的概率密度f(x)和概率分布函数F(x)分别为:,数学期望:E(X)=1/2;方差:=1/12,下图为均匀分布的密度函数和分布函数示意图：,均匀随机数是产生其它随机数的基础。例如，抛硬币、抽签、统计经验分布都可以由它产生。,产生随机数的方法：,（1）随机数表1927年,4万随机数表,以后有100万随机数表（可以输入内存，随时调用）；,（2）硬件设备从真实物理现象的随机因素中产生随机数，放射性粒子的放射源，电子晶体管的固有噪音等，单位时间内放射出的粒子数是随机的。优点：真正的随机数；缺点：外部设备，无法重复,（3）数学公式用数学公式或位移寄存器的移位操作来产生的随机数,称为伪随机数.因为真实的随机数,只能从客观真实的随机现象本身产生出来,所以产生理想的伪随机数列不是一件容易的事.,一般对于产生伪随机数的方法,有如下几点要求:1)要求伪随机数列有较理想的随机性和均匀性,就是对其随机性和均匀性进行统计检验时,有合乎要求的精度;2)产生伪随机数的程序应当简短、运算速度快,占用计算机的内存单元少;3)伪随机数列的循环周期应当尽可能地大,以满足模拟的需要;4)伪随机数列中,前后之间和各子列之间,要求相互是独立无关的。,当我们需要一系列均匀分布的随机数时,可以按照一定的算法由计算机计算产生的伪随机数产生随机数的方法很多，其中以乘同余法使用较广用以产生均匀分布随机数的乘同余法的递推公式为：,xn+1=xn(modM),rn=xn/M,其中，是乘因子，M是模数前面式子的右端称为以M为模数（modulus）的同余式。给定了一个初值x0(称为种子)后，计算出的rn即为(0,l)上均匀分布的随机数,C语言中与随机数有关的函数：randomize(void)初始化x=random(intM)产生0M之间的随机数x=rand(void)产生0215-1之间的随机数Matlab中与随机数有关的函数：Rand(n)n个0,1之间均匀分布随机数Rand(m,n)m*n个0,1之间均匀分布随机数randn(n)n个N(0,1)标准正态分布随机数randn(m,n)m*n个N(0,1)标准正态分布随机数,无论用哪一种方法产生的随机数都存在这样的问题，必须对它进行统计检验看看它们是否具有较好的独立性和均匀性一般在计算机(或计算器)及其使用的算法语言中都有随机数生成的命令，它们所生成的随机数都是经过检验并且可用的.这里就不再详细介绍检验的方法了,2随机变量的模拟,利用均匀分布的随机数可以产生具有任意分布的随机变量的样本，从而可以对随机变量的取值情况进行模拟,P（p(n-1)Rp(n)）=p(n)-p(n-1)=pi，n=1,2,（l）离散型随机变量的模拟,由此可知，事件(p(n-1)Rp(n)和事件(X=xn)有相同的发生的概率因此我们可以用随机变量R落在小区间内的情况来模拟离散的随机变量X的取值情况具体执行的过程是：产生一个(0,1)上均匀分布的随机数r(简称随机数)，若p(n-1)rp(n)则理解为发生事件(X=xn)于是就可以模拟随机变量的取值情况,（2）连续型随机变量的模拟,一般说来，具有给定分布的连续型随机变量可以利用在区间(0,1)上均匀分布的随机数来模拟最常用的方法是反函数法,由概率论的理论可以证明，若随机变量Y有连续的分布函数F(y)，而X是区间(0,1)上均匀分布的随机变量，令，则Z与Y有相同的分布由此，若已知Y的概率密度为，由可得,反函数法,是区间(0,l)上均匀分布的随机变量如果给定区间(0,1)上均匀分布的随机数ri，则具有给定分布的随机数yi可由方程解出,例当需要模拟服从参数为的指数分布时,由,可得.因为(1-ri)和ri同为(0,1)区间上的均匀分布的随机数，故上式可以简化为:,反函数法是一种普通的方法，但当反函数不存在或难以求出时，就不宜于使用了舍选法是另一种方法，其实质是从许多均匀随机数中选出一部分，使之成为具有给定分布的随机数步骤如下：,设随机变量X有密度f(x)，又存在实数ab，使Pr(aXb)=1.,1).选取常数，使f(x)l，x(a,b);,2).产生均匀的随机数r1和r2，令y=a(b-a)r1；,舍选法,3).若r2f(y)，则令x=y，否则剔除r1和r2，重返步骤2，如此重复循环，产生的随机数x1,x1,.,xn,服从的概率分布由密度函数f(x)确定,若不存在上述的有限区间，可以选择一个有限区间(a1,b1)使得Pr(a11-其中是充分小的正数重复上面的步骤，所产生的随机数仅会出现较小的误差,正态随机数的模拟,除了可用反函数和舍选法模拟正态随机变量外,还有两种常用的方法,坐标变换法:设r1,r2是相互独立的均匀随机数,令x1(-2lnr1)1/2cos(2r2)x2=(-2lnr1)1/2sin(2r2)则x1,x2是相互独立的标准正态的随机数,中心极限定理：设是n个相互独立的(0,1)上的均匀的随机变量，有,由中心极限定理知,将渐近服从正态分布N(0,l)因此，取n个均匀的随机数ri,则有,y将可看成是标准正态分布N(0,1)的随机数这里的n应取得足够大，一般取n=10即可若取n=12,则上式简化为再由z=y+得到正态分布N(,2)的随机数z.,计算机仿真举例,例1.4赶火车过程仿真一列火车从A站经过B站开往C站，某人每天赶往B站乘这趟火车。已知火车从A站到B站运行时间为均值30分钟、标准差为2分钟的正态随机变量.火车大约在下午1点离开A站。离开时刻的频率分布为,用计算机仿真火车开出、火车到达B站、这个人到达B站情况，并给出能赶上火车的仿真结果。,这个人到达B站时的频率分布为：,引入以下变量：T1火车从A站开出的时刻；T2此人到达B站的时刻；T3火车从A站运行到B站所需要的时间；,s1=0;s2=0;s3=0;x=rand(10000,1);fori=1:10000ifx(i)0.9s3=s3+1;endends1/10000,1-s1/10000-s3/10000,s3/10000,火车从A站出发时刻T1的仿真测试,s1=0;s2=0;s3=0;s4=0;x=rand(10000,1);fori=1:10000ifx(i)0.3s1=s1+1;elseifx(i)0.7s2=s2+1;elseifx(i)0.9s3=s3+1;elses4=s4+1;endendends1/10000,s2/10000,s3/10000,s4/10000,人到达时刻T2仿真测试,x=randn(10000,1);fori=1:10000y(i)=30+2*x(i);end,火车运行时间T3的仿真测试,赶上火车的仿真结果s=0;x1=rand(10000,1);x2=rand(10000,1);x3=randn(10000,1);fori=1:10000ifx1(i)0.7T1=0;elseifx1(i)0.9T1=5;elseT1=10;endT3=30+2*x3(i);,ifx2(i)0.3T2=28;elseifx2(i)0.7T2=30;elseifx2(i)0?,Q=Q-1,W=2+(4-2)*RAN(1),T2T?,结束,忙,否,xT1)if(Q0)Q-=1;x=rand()%1000;x/=1000;W=(2+(4-2)*x);T3+=W;T1=T2+W;elseT4+=1;while(!(T2=T);P=T3/T;printf(P=%10.5fn,P);,时间步长法例子,运行结果：第1次：0.50627第2次：0.60152第3次：0.59153第4次：0.62537第5次：0.66309第6次：0.57581第7次：0.57048第8次：0.65604第9次：0.53166第10次：0.56407平均：0.58858,应用举例-池水含盐量问题,例1:池水含盐量问题某水池有2000m3水，其中含盐2kg，以每分钟6m3的速率向水池内注入含盐率为0.5kg/m3的盐水，同时又以每分钟4m3的速率从水池流出搅拌均匀的盐水试用计算机仿真该水池内盐水的变化过程，并每隔10min计算水池中水的体积、含盐量和含盐率欲使池中盐水的含盐率达到0.2kg/m3，需经过多少时间？,应用举例-池水含盐量问题,池水含盐量随时间变化的示意图,应用举例-池水含盐量问题,（1）理论分析当系统开始运行后，池中的水和含盐量均随时间变化。记x(t)为t时刻池水中含盐量，考察t,t+t时间区间内盐的变化，这些变化是由于流入的盐水和流出的盐水所引起的，于是有平衡方程:x(t+t)-x(t)=60.54x(t)/(2000+2t)t两边同时除以t,再令t趋于0得下面微分方程的初始值问题。dx(t)/dt=3-4x(t)/(2000+2t)x(0)=2.,利用Mathematica求解此问题。指令为DSolvext=3-4*xt/(2000+2*t),x0=2,xt,t运行后得到的结果为x(t)=(2000000+3000000t+3000t2+t3)/(1000+t)2,应用举例-池水含盐量问题,（2）计算机仿真分析:系统中，实体是水，属性是水的体积、含盐量和含盐率，活动是水的注入与流出，由于注入和流出活动的作用，池中水的体积与含盐量、含盐率均随时间变化，初始时刻含盐率为0.001kg/m3,以后每分钟注入含盐率为0.5kg/m3的水6m3，流出混合均匀的盐水4m3,当池中水的含盐率达到0.2kg/m3时,结束仿真过程.,应用举例-池水含盐量问题,为了能定量地考察系统在任一时刻的属性，引入下列记号：注水速度：WI6m3/min排水速度：WO4m3/min注入水的含盐率：SI0.5kg/m3最终含盐率：SF0.2kg/m3T时刻水的体积：VTm3T时刻水的含盐量：STkgT时刻水的含盐率：SR=ST/VTkg/m3,应用举例-池水含盐量问题,对于这样一个连续系统仿真时，必须在一系列离散时间点t0t1t2=SF)break;if(T%10=0)printf(Time=%d,VT=%.f,ST=%.4f,SR=%.5fn,T,VT,ST,SR);while(1);printf(Time=%d,VT=%.f,ST=%.4f,SR=%.5fn,T,VT,ST,SR);return;,池水含盐量仿真程序,%Matlab程序为h=1;s0=2;w0=2000;r0=s0/w0;s(1)=s0+0.5*6*h-4*h*r0;w(1)=w0+2*h;r(1)=s(1)/w(1);t(1)=h;y(1)=(2000000+3000000*h+3000*h2+h3)/(1000+h)2;fori=2:200t(i)=i*h;s(i)=s(i-1)+0.5*6*h-4*h*r(i-1);w(i)=w(i-1)+2*h;r(i)=s(i)/w(i);y(i)=(2000000+3000000*t(i)+3000*t(i)2+t(i)3)/(1000+t(i)2;m=floor(i/10);ifi/10-m0.2t02=i*h;r02=r(i);breakendendt02,r0210*tm,sm,rmplot(t,s,blue,t,y,red),Matlab程序,池水含盐量仿真结果,应用举例-池水含盐量问题,表给出了不同时刻池中水的体积、含盐量及含盐率，最后一行的数据给出了池水中含盐率达到0.2kg/m3时的时间和含盐量,在这个问题中，池水的体积和含盐量、含盐率都是随时间连续变化的故这是一个连续系统的仿真问题，我们是把这个连续过程离散化后用时间步长法进行仿真的我们还可以用微分方程建立这个问题的数学模型，并求出它的解析解，这个解析解就是问题的精确解有兴趣可以按照这一思路求出该问题的精确解，考察相应时刻精确解与仿真解的差异也可进一步调整前面仿真过程的时间步长，通过与精确解的比较来研究时间步长的大小对仿真精度的影响,应用举例-面积计算,例2:面积计算求由曲线x2+y2=16,x2/36+y2=1,以及(x-2)2+(y+1)2=9所围成图形的面积。,用Matlab语句作出图形得区域图t=0:0.01:2*pi;x=sin(t);y=cos(t);plot(4*x,4*y,6*x,y,2+3*x,3*y-1)axis(-6,6,-6,6),应用举例-面积计算,应用举例-面积计算,此区域形状复杂，理论分析困难，可以用计算机仿真实现。将可能的区域等分，考察每个小区域是否在此区域中，将在此区域中的小面积相加即可其仿真图如下：,应用举例-面积计算,Matlab程序为x=-2:0.01:6;y=-2:0.01:2;s=0;h=0.01;fori=1:800forj=1:400 xx=-2+i*h;yy=-2+j*h;ifxx2+yy2=16ifxx2/36+yy2=1if(xx-2)2+(yy+1)2=9s=s+h2;endendendendends运行后给出面积的值8.8310。Matlab程序为,应用举例-库存问题,例:库存问题在物资的供应过程中，由于到货与销售不可能做到同步、同量，故总要保持一定的库存储备如果库存过多，就会造成积压浪费以及保管费的上升；如果库存过少，会造成缺货如何选择库存和订货策略，就是一个需要研究的问题库存问题有多种类型，一般比较复杂，下面讨论一种简单的情况,某自行车商店的仓库管理人员采取一种简单的订货策略，当库存量降低到P辆自行车时就向厂家订货，每次订货Q辆，如果某一天的需求量越过了库存量，商店就有销售损失和信誉损失，但如果库存量过多，将会导致资金积压和保管费增加若现在已有如下表所示的五种库存策略,试比较选择一种策略以使总费用最少,应用举例-库存问题,这个问题的已知条件是：1）从发出订货到收到货物需隔3天2）每辆自行车保管费为0.75元/天，每辆自行车的缺货损失为1.80元/天，每次的订货费为75元3）每天自行车需求量是0到99之间均匀分布的随机数.4）原始库存为115辆，并假设第一天没有发出订货.,应用举例-库存问题,这一问题用解析法讨论比较麻烦，但用计算机按天仿真仓库货物的变动情况却很方便我们以150天为例，依次对这五种方案进行仿真，最后比较各方案的总费用，从而就可以作出决策,输入一些常数和初始数据后，以一天为时间步长进行仿真.1)检查这一天是否为预定到货日期，如果是，则原有库存量加Q，并把预定到货量清为零；如果不是，则库存量不变;2)接着仿真随机需求量，这可用计算机语言中的随机函数得到若库存量大于需求量，则新的库存量减去需求量；反之，则新库存量变为零，并且要在总费用上加一缺货损失3)检查实际库存量加上预定到货量是否小于重新订货点P，如果是，则需要重新订货，这时就加一次订货费如此重复运行150天，即可得所需费用总值,库存问题仿真流程,应用举例-库存问题,按照图的流程编写程序，上机运行后得这五种库存策略的总费如表所示,比较这五种方案的总费用，可以看出方案2最好，即库存管理人员应取最低订发点为125辆、每次订150辆自行车的方案这时在150天中总的费用是17147.85元,应用举例库存问题,#include#includevoidmain(void)/重新订货点floatP5=125,125,150,175,175;/重新订货量floatQ5=100,150,250,250,300;/库存量floatS5=115,115,115,115,115;floatZ5=0,0,0,0,0;/总费用intT=0;/仿真时间int,day5=0,0,0,0,0;/预定日期intx;/需求量inti;/临时变量intYD5=0,0,0,0,0;/预定量randomize();/初始化随机数,doT+=1;for(i=0;i150);for(i=0;i5;i+)printf(%d%.2fn,i+1,Zi);,应用举例-报童的策略,例4报童每天清晨从报社购进报纸零售,晚上将没有卖掉的报纸退回.每份报纸的购进价为1.3元,零售价为2元,退回价为0.2元.报童售出一份报纸赚0.7元,退回一份报纸赔1.1元.报童每天如果购进的报纸太少,不够卖时会少赚钱,如果购得太多卖不完时要赔钱.试为报童筹划每天应如何确定购进的报纸数使得收益最大.报纸每捆10张，只能整捆购买，报纸可以分为3种类型的新闻日：好、一般、差，它们的概率分别为0.35,0.45和0.2,在这些新闻日中每天对报纸的需求分布的统计结果下图：,应用举例-报童的策略,试确定每天报童应该订购的报纸数量,应用举例-报童的策略,解：我们通过计算机仿真来解决此问题。最优策略应该是每天的利润最大。利润=销售收入-报纸成本-损失+残值这是一个随机现象的计算机仿真问题，故先确定各种情况的随机数的对应关系。新闻日和需求量对应的随机数分别如下面两个表格所示,应用举例-报童的策略,计算机仿真的流程：1）令每天的报纸订购数变化，40-100；2）让时间从1开始变化（循环）到360；3）产生新闻种类的随机数，确定当天的新闻类型；4）产生需求量随机数，确定当天的报纸需求量；5）计算当天的收入，计算累积利润，6）比较得出最优定货量。,应用举例-报童的策略,具体的计算由Matlab编程计算实现。x1=rand(365,1);x2=rand(365,1);forn=4:10paper=n*10;sb(n)=0;fori=1:365ifx1(i)0.35ifx2(i)0.03news=40;elseifx2(i)0.08news=50;elseifx2(i)0.23news=60;elseifx2(i)0.43news=70;,elseifx2(i)0.78news=80;elseifx2(i)0.93news=90;elsenews=100;endelseifx1(i)0.8ifx2(i)0.10news=40;elseifx2(i)0.28news=50;elseifx2(i)0.68news=60;elseifx2(i)0）。每个工人技术相同,且修理时间相互独立并服从同一负指数分布,平均修理时间为1/u(u0).修复后的机器其寿命分布不变。机器停止运转每单位时间的损失费为c1元，工人单位时间的产值为c2元。若机器的等待时间为E，工人总的空闲时间为F，则系统总的损失费为S=c1E+c2F试求当机器数m固定时，为使系统的总损失费最小，应配备多少工人为最优？,事件步长法-例子,解:采用最短时间事件步长法进行仿真，首先要确定一个初始状态。不妨假定开始时所有的机器都正常运转，工人处于空闲状态。设T为仿真终止时间。依次仿真计算修理工人数C分别取1、2、3、时的系统损失费。最后根据系统损失费S的极小值来确定最优工人数。根据以上思想并采用最短时间事件步长法仿真。,事件步长法-例子,事件步长法-例子,框图中各标识符号含义如下：T:仿真终止时间m:机器数c:工人数B(i):第i台机器所处状态A(i):第i台机器现在状态的当前时刻D(j):第j个工人修复机器的释放时刻c1:一个工人单位时间的产值c2:一台机器单位故障时间的损失费,B(i)=,1机器正常运转,0机器发生故障,事件步长法-例子,U:一个工人单位时间平均能修复的机器数V:一台机器在单位运转时间内发生故障的平均次数S:总损失费用E:机器故障时间累计F:工人空闲时间累计K:最短时间机器的当前时刻H:最短时间机器对应的机器号L:工人修完机器的释放时刻x:最短释放时间工人的当前时刻y:最短释放时间对应的工人序号,事件步长法-例子,下表列出了当m=86，1/V=500小时，1/U=34小时，c1=3.46元/小时，c23.2元/小时时的仿真结果，其中每次预定仿真时间为一万个小时，连续仿真五次取其平均值作为仿真结果。由仿真结果可知当工人数c10时，最优工人数为7,应用举例-可靠性问题,例：一设备上三个相同的轴承，每个轴承正常工作寿命为随机变量，其概率分布如表所示,在任何一个轴承损坏都可以使设备停止工作，从有轴承损坏，设备停止工作，到检修工到达开始更换部件为止，称为一个延迟时间延迟时间也是随机变量，其概率分布如下表所示,应用举例-可靠性问题,设备停工时每分钟损失5元,检修工每小时工时费12元，轴承每个成本16元更换一个轴承需要20min，同时更换两个轴承需要30min，同时更换三个轴承需要40min,现在有两种方案：方案一是损坏一个更换一个；方案二是一旦有轴承损坏就全部更换试通过计算机仿真对这两种方案做出评价,在这一问题中，轴承寿命在1000到1900h之间，而延迟在5到15min之间，故若用时间步长法时，步长选取有些困难步长小浪费很大，步长大又不精确，所以采用事件步长法在事件发生时再考虑系统状态的变化情况，这就比较合理,应用举例-可靠性问题,为了进行仿真，首先对轴承寿命和延迟时间与随机数对应，对应规则分别如下两个表,轴承寿命随机数,延迟时间随机数,应用举例-可靠性问题,由于在这一问题中各个轴承的寿命完全决定了系统的运行状态，也即决定了两个方案的费用大小，故我们选择轴承发生故障作为事件，这三个轴承发生故障的事件分别记为A、B、C,（1）方案一的仿真,产生初始事件表,应用举例-可靠性问题,仿真时钟步进，计算费用，产生下一个事件由表上看出，最早发生的事件是A，所以,t1400hcost(5+20)54116145元，,下一个A事件发生的时刻为第2400小时25分钟（随机产生轴承寿命为1000h），刷新事件表，即删去老的A事件，产生新的A事件刷新后的事件表如下表,应用举例-可靠性问题,寻找事件表中的最早事件进行处理,由上表看出，B、C事件同时发生在第1500h，故同时处理时钟步进为t1500h，再根据费用的计算方法得：,cost=145(15+30)56216408元最后利用随机数产生新事件B和C，刷新事件表，得到的新事件表:,应用举例-可靠性问题,重复，t=2400h25min，cost=408(5+20)54116=553元,重复这一过程，一直需要的时间结束即可得到方案的费用,应用举例-可靠性问题,（2）方案二的仿真,方案二与方案一的区别就是一旦故障发生，就更换3个轴承设初始事件表仍为方案一初始表表中最早的事件是A，处理事件A时要考虑延时，更换3个轴承的时间和费用,cost=(5+40)58316281元,应用举例-可靠性问题,根据下一次三个轴承发生故障的时刻刷新后的事件表:,再重复上述过程，累加费用，即可得到方案二的总费用,最后比较两种方案的费用大小即可确定选取那一种,应用举例-可靠性问题,程序运行结果：T=100000h方案一：cost=32705元方案二：cost=24429元从而得出方案二较方案一优。,应用举例-可靠性问题,/可靠性问题#include#include#include#defineTIMES100000structtablefloatbegin;intdelay;intmake;structtabletab3;floatcost=0;intsort3;intmake_event1(void)intfirst=0,sum;sort0=sort1=sort2=-1;floatb=tab0.begin;if(btab1.begin)first=1,b=tab1.begin;if(btab2.begin)first=2,b=tab2.begin;inti=0;,应用举例-可靠性问题,if(first=0)sorti+=0;if(tab1.begin=tab0.begin)sorti+=1;if(tab2.begin=tab0.begin)sorti+=2;elseif(first=1)sorti+=1;if(tab2.begin=tab1.begin)sorti+=2;elseif(first=2)sorti+=2;for(sum=0,i=0;i=0)sum+;switch(sum)case1:tabsort0.make=20;cost+=(20+tabsort0.delay)*5+4+16;break;case2:tabsort0.make=30;cost+=(30+tabsort0.delay)*5+6+32;break;case3:tabsort0.make=40;cost+=(40+tabsort0.delay)*5+8+48;break;returnsum;,应用举例-可靠性问题,intmake_event2(void)intfirst=0,sum;sort0=sort1=sort2=-1;floatb=tab0.begin;if(btab1.begin)first=1,b=tab1.begin;if(btab2.begin)first=2,b=tab2.begin;inti=0;if(first=0)sorti+=0;if(tab1.begin=tab0.begin)sorti+=1;if(tab2.begin=tab0.begin)sorti+=2;elseif(first=1)sorti+=1;if(tab2.begin=tab1.begin)sorti+=2;elseif(first=2)sorti+=2;tabsort0.make=40;cost+=(40+tabsort0.delay)*5+8+48;returnsum;,应用举例-可靠性问题,voidcreate_event1(structtable*t)floatx,y;structtable*t1;x=random(1000);x/=1000;t-begin+=(t-delay+t-make)/60.0;if(xbegin+=1000;elseif(xbegin+=1200;elseif(xbegin+=1300;elseif(xbegin+=1400;elseif(xbegin+=1500;elseif(xbegin+=1600;elseif(xbegin+=1700;elseif(xbegin+=1800;elseif(xbegin+=1900;for(t1=,y=random(1000);y/=1000;if(ydelay=5;elseif(ydelay=10;elseif(ydelay=15;return;,应用举例-可靠性问题,voidcreate_event2(structtable*t,structtable*t0)floatx,y;structtable*t1;x=random(1000);x/=1000;if(t0!=0)t-begin=t0-begin+(t0-delay+t0-make)/60.0;elset-begin+=(t-delay+t-make)/60.0;if(xbegin+=1000;elseif(xbegin+=1200;elseif(xbegin+=1300;elseif(xbegin+=1400;elseif(xbegin+=1500;elseif(xbegin+=1600;elseif(xbegin+=1700;elseif(xbegin+=1800;elseif(xbegin+=1900;for(t1=,y=random(1000);y/=1000;if(ydelay=5;elseif(ydelay=10;elseif(ydelay=15;return;,应用举例-可