计算机仿真-2014数学建模

1、计算机仿真的基本概念随机数的产生时间步长法事件步长法Monte Carlo 方法建模实例一、计算机仿真的基本概念一、计算机仿真的基本概念仿真仿真仿真就是将所研究的对象用其它手段加以模仿的一种活动。实物仿真非实物仿真如曹冲称象、军事演习、飞行器风洞试验、核爆炸试验等，属于实物仿真的例子。1.1.计算机仿真的基本概念计算机仿真的基本概念计算机仿真是一种非实物仿真方法计算机仿真是一种非实物仿真方法 计算机仿真通过建立数学模型、编制计算机程序计算机仿真通过建立数学模型、编制计算机程序实现对真实系统的模拟，对系统的结构和行为进行动实现对真实系统的模拟，对系统的结构和行为进行动态演示态演示, ,从而了解系

2、统随时间变化的行为或特性，以从而了解系统随时间变化的行为或特性，以评价或预测系统的行为效果。它是解决较复杂的实际评价或预测系统的行为效果。它是解决较复杂的实际问题的一条有效途径问题的一条有效途径。计算机仿真反映出新的科学技术的时代特征计算机仿真反映出新的科学技术的时代特征, ,它的应用为各个领域带来新气象和成果。它的应用为各个领域带来新气象和成果。应用的领域有：应用的领域有： 航空管理, 公交车的调度, 飞机设计, 动画设计, 三峡的安全、生态, 道路的修建, 医疗保险, 国债的发行, 家居装修, 炼钢的温度估计, 发电厂的操作训练, 飞行员训练， 鼠疫的检测和预报。 黑死病菌寄生于老鼠身上,

3、 是由跳蚤传染给人类，又叫鼠疫。病菌随着跳蚤叮咬进入人体，约2-5天的潜伏期之后，患者鼠蹊部及其他淋巴结开始红肿、 疼痛， 随之开始发高烧、疲倦、 皮肤变黑，故被称为黑死病。死亡率高达60-90% 鼠疫的传播鼠疫的传播最原始的生化武器最原始的生化武器 1346年，蒙古大将去攻打黑海边富庶的卡法城，久攻不下，这时蒙古军中发生鼠疫，士兵死亡无数，眼看就要无法而退了，这时蒙古将军想出一个方法，把死亡士兵的尸体用弹弩投入城中，迫使城中流行鼠疫，城门自然不攻而破。 在城破时，一位意大利热内亚的富商，帶着妻小和金银珠宝乘船逃了出來。他在地中海各国漂流很久，沒有国家敢收留他们，大家都害怕鼠疫的传染。最后回到

4、家乡热內亚，他把所有的财富全部推在甲板上，对着守城的人说：“我离开卡法城已经六個月了，我若感染鼠疫早就死了，但我并沒有死，可见我并沒有瘟疫。假如你让我进城，甲板上的珠宝就是你们的。” 我们现在知道鼠疫是由老鼠身上的跳蚤传染的，通常老鼠躲在船底污秽处，人们不易察觉。 热内亚人打开城门让这艘船进來后，鼠疫就从热内亚传播开来，传遍整个欧洲，包括北方的斯堪尼亚半岛都无法幸免。每天黃昏时有人推著独轮车，手里摇着铃喊着：“bring out the dead , bring out the dead.”（把尸体拿出來，把尸体拿出來）家家户户就把尸体搬出来丟到城外焚烧，说死尸如山是一点都不为过。鼠疫肆虐欧洲

5、一百多年，使得三分之一的人口死亡。欧洲鼠疫流行时死亡无数鼠疫期间贵族纷纷弃城逃亡 疫如此让人恐怖，那么有没有什么好的预测方式呢？ 计算机仿真就是一个很好的预测方法。 研究发现，鼠疫是由老鼠身上一种特殊的跳蚤传播的。跳蚤的多少决定是否发生鼠疫。跳蚤跳蚤 老鼠老鼠 水草水草 我们可以用计算机根据一个地区的气候模拟出当年此地水草的情况就可以预测出是否有鼠疫要发生。 三峡水库总库容393 亿立方米，总装机容量1820万千瓦，将是世界上最大的水电站。 但是三峡的安全问题是一个很重要的问题，我们不可能等到建好后再看它的安全性，用计算机仿真就可以很好的解决这一问题。 长江三峡工程长江三峡工程飞机设计中有一个

6、重要环节：风洞试验。实际的风洞试验费用巨大。使用计算机仿真进行模拟风洞试验，使费用大大降低。飞机设计飞机设计2.2.为什么要进行计算机仿真为什么要进行计算机仿真 便于重复进行试验，便于控制参数，时间短，代价小。 可以在真实系统建立起来之前，预测其行为效果，从而可以从不同结构或不同参数的模型的结果比较之中，选择最佳模型。 对于缺少解析表示的系统，或虽有解析表示但无法精确求解的系统，可以通过仿真获得系统运行的数值结果。 对于随机性系统，可以通过大量的重复试验，获得其平均意义上的特性指标。3.3.适用计算机仿真解决的问题适用计算机仿真解决的问题 难以用数学公式表示的系统，或者没有建立和求解数学模型的

7、有效方法 虽然可以用解析的方法解决问题，但数学的分析与计算过于复杂，这时计算机仿真可能提供简单可行的求解方法 希望能在较短的时间内观察到系统发展的全过程，以估计某些参数对系统行为的影响 难以在实际环境中进行实验和观察时，计算机仿真是唯一可行的方法，例如太空飞行的研究 需要对系统或过程进行长期运行比较，从大量方案中寻找最优方案4.4.仿真常用术语仿真常用术语系统系统一些具有特定功能相互之间以一定的规律联系着的物体所组成的总体为了限制所研究问题涉及的范围, 用系统边界把所研究的系统与影响系统的环境区分开来.系统的对象、系统的组成元素都可以称为实体系统系统边界边界实体实体属性反映实体的某些性质，它可

8、以是文字型、数字型或逻辑型系统的状态是指在某一时刻实体及其属性值的集合导致系统状态变化的一个过程为活动活动反映了系统变化的规律属性属性状态状态活动活动活动是指一段过程，即在一段时间内发生的情况事件是指一个时间的情况，系统发生变化的瞬间就发生了事件为了使仿真程序能如实地模拟实际系统的变化，在某些离散事件的仿真中，采用事件表的形式进行调度事件表一般是一个有序的记录列，每个记录包括事件发生时间、事件类型等一些内容事件事件事事件件表表研究系统一般是为认识其状态随时间变化的规律，所以需要一个仿真时间变量对连续系统仿真时，常在均匀时间点上展现其状态值，这样，仿真时钟的步进是一个常数对离散系统仿真时，只有在

9、事件发生时，系统的状态才会发生变化，才有必要展现出系统的状态，此时仿真时钟的步进根据事件发生的时刻变化。仿真时钟5.5.仿真研究的步骤仿真研究的步骤明确问题和提出总体方案明确问题和提出总体方案 把被仿真系统的内容表达清楚； 弄清仿真的目的、系统的边界； 确定问题的目标函数和可控变量； 找出系统的实体、属性和活动等。系统分析建立模型建立模型 选择合适的仿真方法（如时间步长法、事件表法等）；确定系统的初始状态；设计整个系统的仿真流程图。收集数据收集数据编写程序、程序验证编写程序、程序验证模型确认模型确认模型构造仿真研究的步骤仿真研究的步骤运行运行 确定具体的运行方案,如初始条件、参数、步长、重复次

10、数等,然后输入数据,运行程序。改进改进 将得出的仿真结果与实际系统比较,进一步分析和改进模型,直到符合实际系统的要求及精度为止。模型的运行与改进结果输出结果输出 设计出结构清晰的仿真结果输出。包括提供文件的清单,记录重要的中间结果等。 输出格式要有利于用户了解整个仿真过程 ,分析和使用仿真结果.设计格式输出仿真结果仿仿真真研研究究步步骤骤问题的阐述设置目标及完整的项目研究计划建立模型收集数据编程序程序验证模型确认试验设计运行与分析进一步运行仿真结束输出结果是是是是否否否否系统分析模型构造模型运行输出结果计算机仿真举例计算机仿真举例例例1.1 1.1 混凝土搅拌中心的位置 有K1(x1,y1),

11、 K2(x2,y2), ., Kn(xn,yn)共n个工地，各需混凝土Q1, Q2,.,Qn吨，混凝土每吨公里的运费为C元。如何确定混凝土搅拌中心的位置K0(x0,y0) 使得费用最少？计算机仿真举例计算机仿真举例方法1：数学计算记第K个工地的位置为(xk,yk),中心的位置为(x0,y0)， 则目标函数为2200001()()min (,)niiiicQxxyyf x y000221000(,)()0()()niiiiif xyQ xxcxxxyy计算机仿真举例计算机仿真举例方法2：计算机的模拟。计算机仿真举例计算机仿真举例方法3：工程模拟 物理、工程仿真的特点：形象、直观、便于类比；但由于

12、是实物仿真,所以费用大,速度慢。计算机仿真举例计算机仿真举例例例1.2 1.2 ：追击问题 我辑私舰雷达发现距c里处有一艘走私船正以匀速度a沿直线行驶,辑私舰立即以最大的速度 b追赶 , 若用雷达进行跟踪,保持船的瞬时速度方向始终指向走私船 , 试求辑私舰追逐路线和追上的时间 ,并且给出缉私舰和走私船的路线图。 选取坐标系如图:假设缉私艇初始位置在点(c,0),走私船初始位置在点(0,0),走私船的行驶方向为y方向.Oxy(c,0)D(x,y).R(0,at)设缉私艇为动点D,走私船为动点R.在时刻t,缉私艇的位置是D(x,y),走私船的位置是R(0,at).)()0(/ )0(22yatxx

13、Cos)()0(/ )(22yatxyatSin取时间间隔(步长)为 ,则在时刻 t+ ,D的位置是 ,tt),(yyxxtSinbytCosbx追赶方向(D的运动方向)为DR.用方向余弦表示为:(*)计算机仿真举例计算机仿真举例算法算法:t赋初值赋初值:初始时刻 t0,时间步长 ,速度a,b,初始位置c循循 环环:),(kkyx由公式(*)计算动点D在各时刻的坐标 D 22122100)(,)(0,kkkkkkkkkkyatxyattbyyyatxxtbxxycx计算动点R在各时刻的坐标R. ),(kkyxtakyxkk0终终 止止: 当D,R的距离小于给定值 时终止计算机仿真举例计算机仿真

14、举例将上述算法中求出的点D用直线段连接起来，就得到追赶曲线。循环终止时的即为追赶时间。可以看出，步长选取得越小，所求的曲线越精确。计算机仿真举例计算机仿真举例例例1.3 1.3 理发店的服务过程仿真 一个理发店有两位服务员A和B，顾客随机地到达该理发店,每分钟有一个顾客到达和没有顾客到达的概率均是1/2 , 其中60%的顾客理发仅用5分钟,另外40%的顾客用8分钟 . 试对前10分钟的情况进行仿真。解：解：假设开始无顾客，顾客到达、服务开始和结束都在每分钟开始时进行， 产生顾客产生顾客：抛硬币，正面-有顾客到来； 反面-无顾客到来 时间长短时间长短：摸球， 3白2黑，白球，5分钟， 黑球，8分

15、钟。 仿真过程仿真过程: .计算机仿真举例计算机仿真举例 时间排队等待人数服务员A服务员B T=0 T=1 T=2 T=3 T=4 T=5 T=6 T=7 T=8 T=9 T=10计算机仿真举例计算机仿真举例例例1.7 1.7 射击命中率 在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点经过长期观察发现，我方指挥所对敌方目标的指示有50是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部消灭敌人 现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来，

16、确定有效射击的比率及毁伤敌方火炮的平均值。计算机仿真举例计算机仿真举例分析分析： 这是一个概率问题，可以通过理论计算得到相应的概率和期望值.但这样只能给出作战行动的最终静态结果，而显示不出作战行动的动态过程. 为了能显示我方20次射击的过程，现采用模拟的方式。计算机仿真举例计算机仿真举例问题分析问题分析: 需要模拟出以下两件事：(1) (1) 指挥所对目标的指示正确与否指挥所对目标的指示正确与否对目标的指示有两种结果，出现的概率都是1/2因此，用投掷一枚硬币的方式予以确定，当硬币出现硬币出现正面时为指示正确，反之为不正确正面时为指示正确，反之为不正确(2) (2) 当指示正确时，我方火力的射击

17、结果当指示正确时，我方火力的射击结果模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6)计算机仿真举例计算机仿真举例毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6)这时可用投掷骰子的方法来确定可用投掷骰子的方法来确定：出现、点：则认为没击中敌人；出现、点：则认为没击中敌人；出现、点：出现、点： 则认为击毁敌人一门火炮；则认为击毁敌人一门火炮；出现点：出现点： 则认为击毁敌人两门火炮则认为击毁敌人两门火炮计算机仿真举例计算机仿真举例2. 符号假设符号假设i：

18、要模拟的打击次数；k1：没击中敌人火炮的射击总数； k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数E：有效射击比率； E1：20次射击平均每次毁伤敌人的火炮数计算机仿真举例计算机仿真举例初始化初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子骰子点数点数?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i20?E=(k2+k3)/20 E1=0*k1/20+1*k2/20+2*k3/20停止停止硬币正面硬币正面?YNNY1，2，34，563. 模拟框图模拟框图计算机仿真举例计算机仿真举例4. 模拟结果模拟结果消消灭灭敌敌人人火火炮炮数数 试试验验 序序号号

19、投投硬硬币币 结结 果果 指指示示 正正确确 指指 示示 不不正正确确 掷掷骰骰子子 结结 果果 正 正 反 正 正 反 正 正 反 反 计算机仿真举例计算机仿真举例消消灭灭敌敌人人火火炮炮数数 试试验验 序序号号 投投硬硬币币 结结 果果 指指示示 正正确确 指指 示示 不不正正确确 掷掷骰骰子子 结结 果果 正 反 正 反 正 正 正 正 反 正 从以上模拟结果可计算出： E=7/20=0.35 20322041201301E=0.5 计算机仿真举例计算机仿真举例5. 理论计算理论计算计算机仿真举例计算机仿真举例则由全概率公式： E = P(A0) = P(j=0)P(A0j=0) + P

20、(j=1)P(A0j=1) = 25. 02121021 P(A1) = P(j=0)P(A1j=0) + P(j=1)P(A1j=1) = 613121021 P(A2) = P(j=0)P(A2j=0) + P(j=1)P(A2j=1) = 1216121021 E1 = 33. 01212611 计算机仿真举例计算机仿真举例6. 结果比较结果比较 理论计算和模拟结果的比较 分类 项目 无效射击 有效射击 平均值 模 拟 理 论 模拟结果与理论计算虽然结果不完全一致，模拟结果与理论计算虽然结果不完全一致，可这反映了事件发生的随机性。只要多次实验求可这反映了事件发生的随机性。只要多次实验求平

21、均值，模拟值就会很接近理论值。平均值，模拟值就会很接近理论值。 二、随机数的产生二、随机数的产生 现实世界充满不确定性，我们所研究的现实对象往往难以摆脱随机因素的影响要使我们的数学模型能够较真实地刻画实际对象，必须面对这个现实概率论是用数学的思想和方法处理和研究随机现象的一个有效的工具但是它还难以用来处理复杂系统中的随机性这里介绍使用计算机来模拟随机现象的方法，它经常应用于复杂系统的动态仿真的研究当中是处理复杂系统中随机性的计算机模型也是使用计算机研究和解决实际问题的一条重要途径1 1均匀分布的随机数及其产生均匀分布的随机数及其产生 对随机现象进行模拟，实质上是要给出随机变量的模拟也就是说利用

22、计算机随机地产生一系列数值，它们的出现服从一定的概率分布,则称这些数值为随机数。最常用的是在(0,1)区间内均匀分布的随机数，也就是我们得到的这组数值可以看作是(0,l)区间内均匀分布的随机变量的一组独立的样本值以后我们将指出其它分布的随机数可利用均匀分布的随机数产生记 XU(0,1),0,1)区间上均匀分布的随机变量X的概率密度f(x)和概率分布函数F(x)分别为:111000)(0101)(xxxxxFxxf其它数学期望: E(X)=1/2; 方 差: = 1/12 下图为均匀分布的密度函数和分布函数示意图： 均匀随机数是产生其它随机数的基础。例如，抛硬币、抽签、统计经验分布都可以由它产生

23、。 产生随机数的方法：（1）随机数表随机数表 1927年,4万随机数表,以后有100万随机数表（可以输入内存，随时调用）；（2）硬件设备硬件设备 从真实物理现象的随机因素中产生随机数，放射性粒子的放射源，电子晶体管的固有噪音等，单位时间内放射出的粒子数是随机的。 优点：真正的随机数； 缺点：外部设备，无法重复（3）数学公式数学公式 用数学公式或位移寄存器的移位操作来产生的随机数,称为伪随机数.因为真实的随机数, 只能从客观真实的随机现象本身产生出来,所以产生理想的伪随机数列不是一件容易的事. 一般对于产生伪随机数的方法,有如下几点要求:1) 要求伪随机数列有较理想的随机性和均匀性,就是对其随机

24、性和均匀性进行统计检验时, 有合乎要求的精度;2) 产生伪随机数的程序应当简短、运算速度快,占用计算机的内存单元少;3) 伪随机数列的循环周期应当尽可能地大,以满足模拟的需要;4) 伪随机数列中,前后之间和各子列之间,要求相互是独立无关的。 当我们需要一系列均匀分布的随机数时,可以按照一定的算法由计算机计算产生的伪随机数产生随机数的方法很多，其中以乘同余法使用较广用以产生均匀分布随机数的乘同余法的递推公式为：xn+1= xn(mod M), rn=xn/M其中，是乘因子，M是模数前面式子的右端称为以M为模数（modulus）的同余式。给定了一个初值x0(称为种子)后，计算出的rn即为(0,l)

25、上均匀分布的随机数C语言中与随机数有关的函数：语言中与随机数有关的函数：randomize(void) 初始化初始化x=random(int M) 产生产生0 0M M 之间的随机数之间的随机数x=rand(void) 产生产生0 02 21515-1-1之间的随机数之间的随机数MatlabMatlab中与随机数有关的函数：中与随机数有关的函数：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

26、)标准正态分布随机数标准正态分布随机数 无论用哪一种方法产生的随机数都存在这样的问题，必须对它进行统计检验看看它们是否具有较好的独立性和均匀性一般在计算机(或计算器)及其使用的算法语言中都有随机数生成的命令，它们所生成的随机数都是经过检验并且可用的.这里就不再详细介绍检验的方法了2随机变量的模拟随机变量的模拟 利用均匀分布的随机数可以产生具有任意分布的随机变量的样本，从而可以对随机变量的取值情况进行模拟P（p(n-1)R p(n)）= p(n) - p(n-1)= pi ，n=1,2,设随机变量X的分布律为 P(X=xi)=pi,i=1,2,令 p(0)=0，p(n)= pi，n=1,2,，将

27、p(n)作为分点，把区间(0 , 1)分为一系列小区间(p(n-1), p(n)对于均匀的随机变量RU(0,l)，则有ni=1（l l）离散型随机变量的模拟）离散型随机变量的模拟 由此可知，事件(p(n-1)R p(n)和事件(X=xn) 有相同的发生的概率因此我们可以用随机变量R落在小区间内的情况来模拟离散的随机变量X的取值情况具体执行的过程是： 产生一个(0,1)上均匀分布的随机数 r(简称随机数)，若p(n-1)r p(n)则理解为发生事件(X=xn)于是就可以模拟随机变量的取值情况（2 2）连续型随机变量的模拟）连续型随机变量的模拟 一般说来，具有给定分布的连续型随机变量可以利用在区间

28、(0,1)上均匀分布的随机数来模拟最常用的方法是反函数法反函数法 YdyyfYFX)()(由概率论的理论可以证明，若随机变量Y有连续的分布函数F(y)，而X是区间(0,1)上均匀分布的随机变量，令 ，则Z与Y有相同的分布由此，若已知Y的概率密度为 ，由 可得)(XFY1 )(yf)(XFZ1 反函数法反函数法 是区间(0,l)上均匀分布的随机变量 如果给定区间(0,1)上均匀分布的随机数ri，则具有给定分布的随机数yi可由方程解出 iyidyyfYFr0)()( YdyyfYFX)()(例当需要模拟服从参数为的指数分布时,由可得 .因为(1-ri) 和ri 同为(0,1)区间上的均匀分布的随机

29、数，故上式可以简化为:iiyyyiedyer 10)ln(iiry 11 iiryln 1 反函数法是一种普通的方法，但当反函数不存在或难以求出时，就不宜于使用了舍选法舍选法是另一种方法，其实质是从许多均匀随机数中选出一部分，使之成为具有给定分布的随机数步骤如下：设随机变量X有密度f(x)，又存在实数ab，使 Pr(aXb)=1.1).选取常数，使f(x) l，x(a,b);2).产生均匀的随机数 r1和 r2，令 y=a(b-a)r1；舍选法舍选法3).若r2 f(y)，则令x=y，否则剔除r1和r2，重返步骤2，如此重复循环，产生的随机数x1, x1,., xn,服从的概率分布由密度函数f

30、(x)确定若不存在上述的有限区间，可以选择一个有限区间(a1,b1)使得 Pr(a1X1-其中 是充分小的正数重复上面的步骤，所产生的随机数仅会出现较小的误差正态随机数的模拟,除了可用反函数和舍选法模拟正态随机变量外,还有两种常用的方法坐标变换法坐标变换法:设r1,r2是相互独立的均匀随机数,令 x1(-2lnr1)1/2cos(2r2) x2 = (-2lnr1)1/2sin(2r2)则x1, x2是相互独立的标准正态的随机数中心极限定理中心极限定理：设 是n个相互独立的(0,1)上的均匀的随机变量，有nixi,.,21 nixDxEii,.,)(,)(2112121 )(21211nxny

31、nii )(21211nrnynii 由中心极限定理知将渐近服从正态分布N(0,l)因此，取n个均匀的随机数ri,则有 y将可看成是标准正态分布N(0,1)的随机数这里的n应取得足够大，一般取n=10即可若取n=12,则上式简化为 再由z=y+得到正态分布N(, 2)的随机数z.61 niiry计算机仿真举例计算机仿真举例出发时刻（T） 1：00 1：051：10 频率 0.7 0.2 0.1例1.4 赶火车过程仿真一列火车从一列火车从A A站经过站经过B B站开往站开往C C站，某人每天赶往站，某人每天赶往B B站站乘这趟火车。已知火车从乘这趟火车。已知火车从A A站到站到B B站运行时间为

32、均值站运行时间为均值3030分钟、标准差为分钟、标准差为2 2分钟的正态随机变量分钟的正态随机变量. .火车大约火车大约在下午在下午1 1点离开点离开A A站。离开时刻的频率分布为站。离开时刻的频率分布为 0.1 0.2 0.4 0.3 频率1：34 1:32 1:30 1:28到达时刻（T）用计算机仿真火车开出、火车到达B站、这个人到达B站情况，并给出能赶上火车的仿真结果。这个人到达B站时的频率分布为： 0.1 0.2 0.4 0.3 P(频率)34323028T2(分)引入以下变量：T1 火车从A站开出的时刻；T2 此人到达B站的时刻；T3火车从A站运行到B站所需要的时间； 0.1 0.2

33、 0.7 P(频率) 10 5 0 T1(分),(22303NT s1=0; s2=0; s3=0; x=rand(10000,1); for i=1:10000 if x(i)0.9 s3=s3+1; end ends1/10000, 1-s1/10000-s3/10000,s3/10000 0.1 0.2 0.7 P(频率) 10 5 0 T1(分)火车从火车从A A站出发时刻站出发时刻T1T1的仿真测试的仿真测试 T2（分） 28 30 32 34 P( 频率) 0.3 0.4 0.2 0.1 s1=0; s2=0; s3=0;s4=0; x=rand(10000,1); for i=1

34、:10000 if x(i)0.3 s1=s1+1; elseif x(i)0.7 s2=s2+1; else if x(i)0.9 s3=s3+1; else s4=s4+1; end endends1/10000, s2/10000,s3/10000,s4/10000人到达时刻人到达时刻T2T2仿真测试仿真测试x=randn(10000,1);for i=1:10000 y(i)=30+2*x(i);end火车运行时间火车运行时间T3的仿真测试的仿真测试),(22303NT赶上火车的仿真结果赶上火车的仿真结果s=0;x1=rand(10000,1);x2=rand(10000,1);x3=

35、randn(10000,1);for i=1:10000if x1(i)0.7 T1=0;elseif x1(i)0.9 T1=5;else T1=10;endT3=30+2*x3(i);if x2(i)0.3 T2=28;elseif x2(i)0.7 T2=30; else if x2(i)0.9 T2=32; else T2=34; endendif T2 T1?Q0?Q=Q-1W=2+(4-2)*RAN(1)T2T?结束忙否x0.2T3= T3+WT4= T4+1 是否否是是闲T1= T2+WP= T3/T输出PT1=0,T2=0,T3=0,T4=0,Q=0 x=RND(1)T:预定仿

36、真时间T1:保管员的释放时刻T2:仿真时钟时刻T3:保管员忙时累计T4:保管员空闲时间累计Q:当前时刻要求领料的工人数W:领料时间P:保管员忙的概率时时间间步步长长法法例例子子/仿真程序，C语言编写#include #include void main(void) float T,T1,T2,T3,T4; int Q; double x,P,W,S; scanf(“%f”,&T);/ T=1000; T1=T2=T3=T4=Q=W=P=0; do T2+=1; x=rand()%1000; x/=1000; if(xT1) if(Q0) Q-=1; x=rand()%1000; x/=

37、1000; W=(2+(4-2)*x); T3+=W; T1=T2+W; else T4+=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:池水含盐量问题某水池有 2 000m3水，其中含盐 2 kg，以每分钟 6

38、m3的速率向水池内注入含盐率为0.5 kg/m3的盐水，同时又以每分钟4m3的速率从水池流出搅拌均匀的盐水试用计算机仿真该水池内盐水的变化过程，并每隔10min计算水池中水的体积、含盐量和含盐率欲使池中盐水的含盐率达到0.2kg/m3，需经过多少时间？应用举例应用举例- -池水含盐量问题池水含盐量随时间变化的示意图应用举例应用举例- -池水含盐量问题（1）理论分析）理论分析 当系统开始运行后，池中的水和含盐量均随时间变化。记x(t)为t时刻池水中含盐量，考察 t, t+t时间区间内盐的变化，这些变化是由于流入的盐水和流出的盐水所引起的，于是有平衡方程: x(t+ t)-x(t)= 6 0.5

39、4 x(t)/(2000+2t) t两边同时除以t,再令t趋于0得下面微分方程的初始值问题。 dx(t)/dt=3- 4x(t)/(2000+2t) x(0)=2.利用Mathematica 求解此问题。 指令为DSolve xt=3-4*xt/(2000+2*t), x0=2, xt,t运行后得到的结果为x(t)=(2000000 + 3000000 t + 3000 t2 + t3)/(1000 + t)2,应用举例应用举例- -池水含盐量问题（2） 计算机仿真计算机仿真分析分析:系统中，实体是水，属性是水的体积、含盐量和含盐率，活动是水的注入与流出，由于注入和流出活动的作用，池中水的体积

40、与含盐量、含盐率均随时间变化，初始时刻含盐率为0.001kg/m3,以后每分钟注入含盐率为0.5 kg/m3的水6 m3，流出混合均匀的盐水4m3,当池中水的含盐率达到0.2kg/m3时,结束仿真过程.应用举例应用举例- -池水含盐量问题为了能定量地考察系统在任一时刻的属性，引入下列记号： 注水速度： WI6 m3/min 排水速度： WO 4 m3/min 注入水的含盐率：SI0.5 kg/m3 最终含盐率：SF0.2 kg/m3 T时刻水的体积：VTm3 T时刻水的含盐量：ST kg T时刻水的含盐率：SR=ST/VT kg/m3应用举例应用举例- -池水含盐量问题对于这样一个连续系统仿真

41、时，必须在一系列离散时间点t0 t1 t2 tn上来进行考察，这些离散时间点之间的间隔T= ti- ti-1(i=1,2,n)就是时间步长若取步长为1min时，就要每隔1min考察一次系统的状态特性，并作出相应的记录与分析在本题中每隔1min池水的动态变化过程是： 每分钟水的体积增加 6-42 m3； 每分钟向池内注入盐 60.5 3 kg； 每分钟向池外流出盐 4SR kg； 每分钟池内增加盐 3-4SR kg取T0作为系统仿真的初始时刻，取T=0时系统的属性为系统的初始状态取池中盐水的含盐率达到 0.2kg/m3的时刻作为仿真的终止时刻池池水水含含盐盐量量仿仿真真流流程程初始化按时间步长前

42、进1min计算池水体积VT,含盐量ST,含盐率SR是否到达10min？结束输出结果是否含盐率达到SF?打印次数加1输出时间T,水体积VT含盐量ST,含盐率SR下一个10min开始，计时单元清零是否池池水水含含盐盐量量仿仿真真程程序序/*C语言片编写的仿真程序 #include #include void main(void) int T=0; float WI,WO,SI,SF,VT,ST,SR; SF=0.2; VT=2000.0;ST=2.0;SR=2.0/2000; do T+=1; VT=VT+(6-4); ST=ST+(3-4*SR); SR=ST/VT; if(SR=SF)brea

43、k; 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

44、)2;for i=2:200 t(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); if i/10-m0.2 t02=i*h; r02=r(i); break endendt02,r0210*tm,sm,rmplot(t,s,blue,t,y,red)Matlab程序池池水水含含盐盐量量仿仿真真结结果果时间T/min池水体积VT/m3含盐量ST/kg含盐率

45、SR/kg/m310202031.693650.0156920204060.810080.0290930206089.371620.04338402080117.39960.05644502100144.91410.06901602120171.93460.08110702140198.47950.09275802160224.56630.10397902180250.21180.114781002200275.43190.125191102220300.2420.135241202240324.65660.144941302260348.68970.154291402280372.35460

46、.163311502300395.66420.172031602320418.63060.180441702340441.26550.188571802360463.58020.196431852370474.62080.20026应用举例应用举例- -池水含盐量问题 表给出了不同时刻池中水的体积、含盐量及含盐率，最后一行的数据给出了池水中含盐率达到0.2kg/m3时的时间和含盐量在这个问题中，池水的体积和含盐量、含盐率都是随时间连续变化的故这是一个连续系统的仿真问题，我们是把这个连续过程离散化后用时间步长法进行仿真的我们还可以用微分方程建立这个问题的数学模型，并求出它的解析解，这个解析解就是

47、问题的精确解有兴趣可以按照这一思路求出该问题的精确解，考察相应时刻精确解与仿真解的差异也可进一步调整前面仿真过程的时间步长，通过与精确解的比较来研究时间步长的大小对仿真精度的影响应用举例应用举例- -面积计算例例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)应用举例应用举例- -面积计算应用举例应用举例- -面积计算此区域形状复杂，理

48、论分析困难，可以用计算机仿真实现。将可能的区域等分， 考察每个小区域是否在此区域中 ，将在此区域中的小面积相加即可其仿真图如下：应用举例应用举例- -面积计算Matlab程序为x=-2:0.01:6;y=-2:0.01:2;s=0;h=0.01;for i=1:800 for j=1:400 xx=-2+i*h; yy=-2+j*h; if xx2+yy2=16 if xx2/36+yy2=1 if (xx-2)2+(yy+1)2=9 s=s+h2; end end end endends运行后给出面积的值 8.8310。 Matlab程序为应用举例应用举例- -库存问题例例: :库存问题在物

49、资的供应过程中，由于到货与销售不可能做到同步、同量，故总要保持一定的库存储备如果库存过多，就会造成积压浪费以及保管费的上升；如果库存过少，会造成缺货如何选择库存和订货策略，就是一个需要研究的问题库存问题有多种类型，一般比较复杂，下面讨论一种简单的情况 某自行车商店的仓库管理人员采取一种简单的订货策略，当库存量降低到P辆自行车时就向厂家订货，每次订货Q辆，如果某一天的需求量越过了库存量，商店就有销售损失和信誉损失，但如果库存量过多，将会导致资金积压和保管费增加若现在已有如下表所示的五种库存策略,试比较选择一种策略以使总费用最少应用举例应用举例- -库存问题这个问题的已知条件是：1）从发出订货到收

50、到货物需隔3天2）每辆自行车保管费为0.75元/天，每辆自行车的缺货损失为1.80元/天，每次的订货费为75元3）每天自行车需求量是0到99之间均匀分布的随机数.4）原始库存为115辆，并假设第一天没有发出订货.方案编号重新订货点P/辆重新订量Q/辆方案1125100方案2125150方案3150250方案4175250方案5175300应用举例应用举例- -库存问题 这一问题用解析法讨论比较麻烦，但用计算机按天仿真仓库货物的变动情况却很方便我们以150天为例，依次对这五种方案进行仿真，最后比较各方案的总费用，从而就可以作出决策 输入一些常数和初始数据后，以一天为时间步长进行仿真.1)检查这一

51、天是否为预定到货日期，如果是，则原有库存量加Q，并把预定到货量清为零；如果不是，则库存量不变;2)接着仿真随机需求量，这可用计算机语言中的随机函数得到若库存量大于需求量，则新的库存量减去需求量；反之，则新库存量变为零，并且要在总费用上加一缺货损失3)检查实际库存量加上预定到货量是否小于重新订货点P，如果是，则需要重新订货，这时就加一次订货费 如此重复运行150天，即可得所需费用总值库库存存问问题题仿仿真真流流程程初始化库存+Q预定到货量置零需求量库存量？结束输出结果是今天是到货期？当天日期+1是否产生随机需求量计算新库存量，累加保管费实际库存+预定到货量p？满150d?总费用+短缺费预定到货量

52、=Q预定到货日期=今天+3是否否是应用举例应用举例- -库存问题按照图的流程编写程序，上机运行后得这五种库存策略的总费如表所示比较这五种方案的总费用，可以看出方案2最好，即库存管理人员应取最低订发点为125辆、每次订150辆自行车的方案这时在150天中总的费用是17147.85元方案编号总费用/元方案124476.70方案217147.85方案321824.40方案424708.75方案527275.85应应用用举举例例库存问题 #include #include void main(void) /重新订货点 float P5=125,125,150,175,175; /重新订货量 float

53、 Q5=100,150,250,250,300; /库存量 float S5=115,115,115,115,115; float Z5=0,0,0,0,0;/总费用 int T=0;/仿真时间 int ,day5=0,0,0,0,0;/预定日期 int x;/需求量 int i;/临时变量 int YD5=0,0,0,0,0;/预定量 randomize();/初始化随机数do T+=1; for(i=0;i5;i+) if(dayi=T)Si+=Qi,YDi=0; x=random(100); for(i=0;i5;i+) Zi+=0.75*Si;/保管费 if(!(x=Si) Zi+=(

54、x-Si)*1.8;/缺货费 Si=0; else Si-=x; if(Si+YDi150); for(i=0;i5;i+) printf(%d %.2fn,i+1,Zi);应用举例应用举例- -报童的策略例例 4 4 报童每天清晨从报社购进报纸零售,晚上将没有卖掉的报纸退回.每份报纸的购进价为1.3元,零售价为2元,退回价为0.2元.报童售出一份报纸赚0.7元 ,退回一份报纸赔1.1元.报童每天如果购进的报纸太少,不够卖时会少赚钱,如果购得太多卖不完时要赔钱. 试为报童筹划每天应如何确定购进的报纸数使得收益最大.报纸每捆10张，只能整捆购买，报纸可以分为3种类型的新闻日：好、一般、差，它们的

55、概率分别为0.35,0.45和0.2,在这些新闻日中每天对报纸的需求分布的统计结果下图：应用举例应用举例- -报童的策略需求量好新闻的需求概率一般新闻的需求概率 差新闻的需求概率 40 0.03 0.10 0.44 50 0.05 0.18 0.22 60 0.15 0.40 0.16 70 0.20 0.20 0.12 80 0.35 0.08 0.06 90 0.15 0.04 0.00 100 0.07 0.00 0.00试确定每天报童应该订购的报纸数量应用举例应用举例- -报童的策略解：我们通过计算机仿真来解决此问题。最优策略应该是每天的利润最大。利润利润= =销售收入销售收入- -报

56、纸成本报纸成本- -损失损失+ +残值残值这是一个随机现象的计算机仿真问题，故先确定各种情况的随机数的对应关系。新闻日和需求量对应的随机数分别如下面两个表格所示新闻种类 出现概率对应的随机数区间好新闻 0.35 (0.00,0.35)一般新闻 0.45 0.35,0.80)差新闻 0.20 0.80,1.00)应用举例应用举例- -报童的策略计算机仿真的流程：1）令每天的报纸订购数变化，40-100；2）让时间从1开始变化（循环）到360；3）产生新闻种类的随机数，确定当天的新闻类型；4）产生需求量随机数，确定当天的报纸需求量；5）计算当天的收入，计算累积利润，6）比较得出最优定货量。需求量好

57、新闻的随机数区间一般新闻的随机数区间差新闻的随机数区间 40 （0.00,0.03 （0.00,0.10 （0.00,0.44 50 0.03,0.08) 0.10,0.28) 0.44,0.66) 60 0.08,0.23) 0.28,0.68) 0.80,1.00) 70 0.23,0.43) 0.68,0.88) 0.80,1.00) 80 0.43,0.78) 0.880.96) 0.80,1.00) 90 0.78,0.93) 0.96,1.00) 100 , 102 100 0.93,1.00 8 , 100 104 , 125应用举例应用举例- -报童的策略具体的计算由具体的计算

58、由Matlab编程计编程计算实现。算实现。x1=rand(365,1);x2=rand(365,1);for n=4:10 paper=n*10; sb(n)=0; for i=1:365 if x1(i)0.35 if x2(i)0.03 news=40; elseif x2(i)0.08 news=50; elseif x2(i)0.23 news=60; elseif x2(i)0.43 news=70;elseif x2(i)0.78 news=80; elseif x2(i)0.93 news=90; else news=100; endelseif x1(i)0.8 if x2(i

59、)0.10 news=40; elseif x2(i)0.28 news=50; elseif x2(i)0.68 news=60; elseif x2(i)0.88 news=70; elseif x2(i)0.96 news=80;应用举例应用举例- -报童的策略else news=90; end else if x2(i)0.44 news=40; elseif x2(i)0.66 news=50; elseif x2(i)0.82 news=60;elseif x2(i)=news sale=news; remand=paper-news; 应用举例应用举例- -报童的策略else s

60、ale=paper; remand=0; end sb(n)=sb(n)+2*sale-1.3*paper+0.2*remand; endendoptnews=40;optmoney=sb(4);40,sb(4)/365for n=5:10 if sb(n)=optmoney optnews=n*10; optmoney=sb(n); end n,sb(n)/365endoptnews,optmoney,optmoney/365 Matlab程序 经过计算机仿真后得到最优购货量是每天60份，平均每天利润34.4元。四、事件步长法四、事件步长法事件步长法是以事件发生的时间为增量，按照事件发生的时间顺序，