模拟建模

1. 模拟建模

2. 模拟的概念及作用 • 现实系统的数学或逻辑模型可能十分复杂，例如大多数具有随机因素的复杂系统，其中的一些随机性因素很难用准确的数学公式表述，从而也无法对整个系统采用解析法求解。模拟是处理这类实际问题的有力工具。

3. 应用领域： • 运输系统模拟 • 摩天大楼安全疏散系统模拟 • 国民经济发展模拟 • 人口增长系统模拟 • 供水系统模拟

4. 蒙特卡洛法 蒙特卡洛(MonteCarlo)法，或称统计试验法、计算机随机模拟方法 是一种基于“随机数”，采用统计抽样方法，近似求解数学问题或物理问题的过程。 • 当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。

5. 例. 蒲丰氏问题 • 为了求得圆周率π值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（l＜a）的平行线相交的频率代替概率P，再利用准确的关系式： • 求出π值 • 其中Ｎ为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。

6. 为了得到具有一定精确度的近似解，所需试验的次数是很多的，通过人工方法作大量的试验相当困难，甚至是不可能的。因此，蒙特卡罗方法的基本思想虽然早已被人们提出，却很少被使用。本世纪四十年代以来，由于电子计算机的出现，使得人们可以通过电子计算机来模拟随机试验过程，把巨大数目的随机试验交由计算机完成，使得蒙特卡罗方法得以广泛地应用，在现代化的科学技术中发挥应有的作用。为了得到具有一定精确度的近似解，所需试验的次数是很多的，通过人工方法作大量的试验相当困难，甚至是不可能的。因此，蒙特卡罗方法的基本思想虽然早已被人们提出，却很少被使用。本世纪四十年代以来，由于电子计算机的出现，使得人们可以通过电子计算机来模拟随机试验过程，把巨大数目的随机试验交由计算机完成，使得蒙特卡罗方法得以广泛地应用，在现代化的科学技术中发挥应有的作用。

7. 计算机模拟试验过程 • 计算机模拟试验过程，就是将试验过程（如投针，射击）化为数学问题，在计算机上实现。 • 由上面例题看出，蒙特卡罗方法常以一个“概率模型”为基础，按照它所描述的过程，使用由已知分布抽样的方法，得到部分试验结果的观察值，求得问题的近似解。 • 从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。

8. 具有同时计算多个方案与多个未知量的能力 • 对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。 • 另外，使用蒙特卡罗方法还可以同时得到若干个所求量，而不像常规方法那样，需要逐一计算所求量。

9. 计算机模拟：在已经建立的数学、逻辑模型的基础之上，通过计算机试验，对一个系统按照一定的决策原则或作业规则，由一个状态变换为另一个状态的行为进行描述和分析。计算机模拟：在已经建立的数学、逻辑模型的基础之上，通过计算机试验，对一个系统按照一定的决策原则或作业规则，由一个状态变换为另一个状态的行为进行描述和分析。 • 系统状态随时间而变化的动态写照

10. 模拟的作用 • 对于很难用解析方法加以处理的问题，模拟是一种有效的技术； • 对建模过程中的假设进行鉴定，对理论研究的结论加以检验； • 对不同的实现方案进行多次模拟，按照既定的目标函数对不同方案进行比较，从中选择最优方案。

11. 模拟的分类 系统状态随时间而变化的动态写照 • 通常，模拟时间是模拟的主要自变量。其他变量作为因变量来看待。 • 推进模拟时间的基本方法： • 下次事件法：将模拟时间由一个事件发生的时间点推进到紧接着的下一次事件发生的时间点。---系统的状态仅在 • 固定时间步长法：模拟时间每次均以相等的固定步长向前推进，每到达一个新的模拟时间点需检查相应时间段内是否发生了事件。需根据实际问题合理设置模拟时间发生改变的步长。

12. 模拟的分类 • 根据模拟过程中因变量的变化情况进行分类： • 第一类：离散型模拟：因变量在与事件时间有关的具体模拟时间点呈离散性变化。大多数系统（如排队服务系统）可采用离散型模拟。 • 排队系统通常采用离散型模拟模型。 任意时刻的排队的顾客数 只有新顾客到来或有顾客完成服务后离去，函数值才发生变化

13. 其中，发生系统状态变化的事件有两个： 一是有顾客到达；二是服务员完成服务。 将最近发生上述两种事件之一的时刻设置为下次事件发生点，就可将服务过程描述为下图所示的模拟模型。

14. 对任意一个模拟，首先要作的是: • 找出能完全描述任意时刻的系统的状态变量的集合： 这个问题中有两类状态变量： 在等待的顾客的人数（离散的非负整数）； 服务员Ai是否正在工作（是或否）； • 给出能从时刻 t 的状态变量算出时刻 t+1 的新的状态变量的程序。 • 从而得到平均队伍的长度，最长的队伍，顾客等待的平均时间以及两个理发员的忙闲程度等，

15. 模拟的分类 • 第二类：连续型模拟。因变量随时间的改变呈连续性变化。在大多数计算机模拟过程中，按固定的步长推进模拟时间。 通常需建立一系列的由系统状态变量组成的状态方程组，以描述状态变量与模拟时间的关系。 • 第三类：混合型模拟。因变量随时间的推移而作连续性的变化并具有离散性的突变，如库存控制系统。

16. 模拟的方式 • 终态模拟：在规定的时间T内进行模拟运行，时间达到T时，模拟终止。其性能指标明显取决于系统的初始状态。 • 稳态模拟：随着模拟时间的推移，系统的性能逐渐趋于平稳。其目的是研究非终态系统长期运行条件下的稳态性能，模拟时间的长短取决于能否获得系统性能的优良估计（可由模拟输出的精度确定）。

17. 模拟的一般步骤 • 明确问题，建立模型。 明确模拟目的，确定模拟输出结果的目标函数 分析各状态变量之间关系，建立系统模型 • 收集和整理数据资料。特别是随机性资料。 分析收集的随机数据，确定系统中随机性因素的概率分布特征，以此为依据产生抽样数据 • 编制程序，模拟运行。 编程、设定初始状态，模拟运行时间、随机样本量、模拟运行次数 • 分析模拟输出结果： 模拟结果的统计特性（样本均值、方差、置信区间等） 灵敏性分析（输入参数做微小变化是否导致输出结果巨大改变） 选择最优方案。

18. 系统 否 是否用仿真 用其他方法 是 构造模型 修改模型 数据准备 编程 修改模型策略或重复运行 程序运行 否 是 分析仿真结果 否 模型是否有效 结果是否充分 是 停止

19. 系统模拟注意事项： • 一次模拟结果毫无意义！ • 模拟是试验性的,是思维结果的验证。 • 必须进行足够多次的模拟，并对结果进行统计分析。

20. 随机数的产生 • 数学方法常见产生均匀分布随机数的几种方法的有平方取中法、倍积取中法、乘同余法、二阶与三阶线性同余法 • 由均匀分布产生各种分布的随机数、反函数法、取舍法、Box-Muller方法和极方法。 • 。

22. randperm(n)：产生一个1到n的随机顺序。 • randint(m,n,[1 N])：生成m×n的在1到N之间的随机整数矩阵， • rand(m,n)：生成0到1之间的m×n的随机数矩阵 • randn是均值为0方差为1的正态分布

24. 例在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏．为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点．例在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏．为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点． 经过长期观察发现，我方指挥所对敌方目标的指示有50％是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部消灭敌人． 现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值。 分析：这是一个概率问题，可以通过理论计算得到相应的概率和期望值.但这样只能给出作战行动的最终静态结果，而显示不出作战行动的动态过程. 为了能显示我方20次射击的过程，现采用模拟的方式。

25. 1. 问题分析 需要模拟出以下两件事： [1] 观察所对目标的指示正确与否 模拟试验有两种结果，每一种结果出现的概率都是1/2． 因此，可用投掷一枚硬币的方式予以确定，当硬币出现正面时为指示正确，反之为不正确． [2] 当指示正确时，我方火力单位的射击结果情况 模拟试验有三种结果：毁伤一门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6)． 这时可用投掷骰子的方法来确定： 如果出现的是１、２、３三个点：则认为没能击中敌人； 如果出现的是４、５点：则认为毁伤敌人一门火炮； 若出现的是６点：则认为毁伤敌人两门火炮．

26. 2. 符号假设 初始化:i=0,k1=0,k2=0,k3=0 i=i+1 Y N 硬币正面? 1，2，3 骰子点数? 4，5 6 k1=k1+1 k1=k1+1 k2=k2+1 k3=k3+1 Y i＜20? N E=(k2+k3)/20 E1=0*k1/20+1*k2/20+2*k3/20 停止 i：要模拟的打击次数； k1：没击中敌人火炮的射击总数； k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数． E：有效射击比率； E1：20次射击平均每次毁伤敌人的火炮数． 3. 模拟框图

27. R2=unidrnd(6)

31. 用蒙特卡洛法解非线性规划问题

32. 基本假设 试验点的第j个分量xj服从[aj ,bj]内的均匀分布． 符号假设 求解过程 先产生一个随机数作为初始试验点,以后则将上一个试验点的第j个分量随机产生，其它分量不变而产生一新的试验点．这样，每产生一个新试验点只需一个新的随机数分量．当K>MAXK或P>MAXP时停止迭代．

33. 初始化:给定MAXK,MAXP;k=0,p=0,Q:大整数 xj=aj+R(bj-aj) j=1,2,…,n j=0 框 图 j=j+1,p=p+1 Y N f(X)≥Q? Y N P>MAXP? xj=aj+R(bj-aj) Y gi(X)≥0? i=1,2…n N Y N j<n? X*=X,Q=f(X) k=k+1 Y N k>MAXK? 输出X,Q,停止

34. % lpconst.m function lpc=lpconst(x) % 约束条件 ± 0 . 5 if 3*x(1)+x(2) - 10<=0.5 & 3*x(1)+x(2) - 10>= - 0.5 % 约束条件的误差为 lpc=1; else lpc=0; end 在Matlab软件包中编程，共需三个Ｍ－文件:randlp.m, mylp.m,lpconst.m.主程序为randlp.m. % mylp.m function z=mylp(x) 　　　　　　　%目标函数 z=2*x(1)^2+x(2)^2-x(1)*x(2)-8*x(1)-3*x(2);　 %转化为求最小值问题

35. % randlp.m function [sol,r1,r2]=randlp(a,b,n)　　 %随机模拟解非线性规划 debug=1; a=0;　　　　　　　　　　　　　 %试验点下界 b=10;　　　　　　　　　　　　 　%试验点上界 n=1000;　　　　　　　　　　　　 %试验点个数 r1=unifrnd(a,b,n,1);　　　　　　%n×1阶的［a,b］均匀分布随机数矩阵 r2=unifrnd(a,b,n,1); sol=[r1(1) r2(1)]; z0=inf; for i=1:n x1=r1(i); x2=r2(i); lpc=lpconst([x1 x2])； if lpc==1 z=mylp([x1 x2])； if z<z0 z0=z； sol=[x1 x2]； end end end

36. 随机数 • 由具有已知分布的总体中抽取简单子样，在蒙特卡罗方法中占有非常重要的地位。总体和子样的关系，属于一般和个别的关系，或者说属于共性和个性的关系。由具有已知分布的总体中产生简单子样，就是由简单子样中若干个性近似地反映总体的共性。 • 随机数是实现由已知分布抽样的基本量，在由已知分布的抽样过程中，将随机数作为已知量，用适当的数学方法可以由它产生具有任意已知分布的简单子样。

37. 由于随机数在蒙特卡罗方法中占有极其重要的位置，我们用专门的符号ξ表示。由随机数序列的定义可知，ξ1，ξ2，…是相互独立且具有相同单位均匀分布的随机数序列。也就是说，独立性、均匀性是随机数必备的两个特点。由于随机数在蒙特卡罗方法中占有极其重要的位置，我们用专门的符号ξ表示。由随机数序列的定义可知，ξ1，ξ2，…是相互独立且具有相同单位均匀分布的随机数序列。也就是说，独立性、均匀性是随机数必备的两个特点。 由于随机数在蒙特卡罗方法中所处的特殊地位，它们虽然也属于由具有已知分布的总体中产生简单子样的问题，但就产生方法而言，却有着本质上的差别。

38. 伪随机数 • 在计算机上产生随机数最实用、最常见的方法是数学方法，即用如下递推公式： • 产生随机数序列。对于给定的初始值ξ1,ξ2…，ξk，确定ξn+k，ｎ=1，2，…。

39. 伪随机数存在的两个问题 • 用数学方法产生的随机数，存在两个问题： • 递推公式和初始值ξ1,ξ2…，ξk确定后，整个随机数序列便被唯一确定。不满足随机数相互独立的要求。 • 由于随机数序列是由递推公式确定的，而在计算机上所能表示的[0，1]上的数又是有限的，因此，这种方法产生的随机数序列就不可能不出现无限重复。 • 随机数序列便出现了周期性的循环现象。

40. 由于这两个问题的存在，常称用数学方法产生的随机数为伪随机数。对于以上存在的两个问题，作如下具体分析。由于这两个问题的存在，常称用数学方法产生的随机数为伪随机数。对于以上存在的两个问题，作如下具体分析。 • 关于第一个问题，不能从本质上加以改变，但只要递推公式选得比较好，随机数间的相互独立性是可以近似满足的。至于第二个问题，则不是本质的。因为用蒙特卡罗方法解任何具体问题时，所使用的随机数的个数总是有限的，只要所用随机数的个数不超过伪随机数序列出现循环现象时的长度就可以了。 • 用数学方法产生的伪随机数容易在计算机上得到，可以进行复算，而且不受计算机型号的限制。因此，这种方法虽然存在着一些问题，但仍然被广泛地在计算机上使用，是在计算机上产生伪随机数的主要方法。

41. 1999年A题 自动化车床管理 一道工序用自动化车床连续加工某种零件，由于刀具损坏等原因该工序会出现故障，其中刀具损坏故障占95%，其它故障仅占5%。工序出现故障是完全随机的，假定在生产任一零件时出现故障的机会均相同。工作人员通过检查零件来确定工序是否出现故障。 现积累有100次刀具故障记录，故障出现时该刀具完成的零件数如附表。现计划在刀具加工一定件数后定期更换新刀具。 已知生产工序的费用参数如下： 故障时产出的零件损失费用 f = 200元/件； 进行检查的费用 x = 10元/次； 发现故障进行调节使恢复正常的平均费用 d = 3000元/次 未发现故障时更换一把新刀具的费用 k = 1000元/次。 1）假定工序故障时产出的零件均为不合格品，正常时产出的零件均为合格品，试对该工序设计效益最好的检查间隔和刀具更换策略。 2）如果该工序正常时产出的零件不全是合格品，有2%为不合格品：而工序故障时产出的零件有40%为合格品，60%为不合格品。工序正常而误认为有故障停机产生的损失费用为1500元/次。对该工序设计效益最好的检查间隔和刀具更换策略。

42. 较复杂的非线性优化问题 1、验证模型的正确性，采用计算机模拟 • 以零件的序列为模拟序列 • 确定随机事件，如主要随机事件--刀具的损害作为，可分析证明刀具寿命服从正态分布 • 编制程序模拟事件发生过程 • 对不同的给定检查间隔和换刀间隔进行模拟，得到单位零件平均工作周期内的资金损耗 • 比较模型结果和模拟结果

43. 2、检验假设 • 如第二问中，加入了随机因子，问题更为复杂，建模过程中不免加入一些假设使问题更简单化。 • 如由于非刀具故障占总故障的百分比较小，为5%，故假设其对结果影响不大，从而不考虑非刀具故障。 • 可通过计算机模拟对此假设进行验证。

44. 3、对多个策略进行仿真比较 • 如对问题二通过分析提出如下几个检查策略（h为检查间隔） 策略1、检查第h个零件，如为正品，则继续生产；如为废品，停止生产， 换刀。 策略2、检查第h个零件，如为正品，则继续生产；如为废品，则检查第 h+1个，若正品，继续生产，如废品，换刀。 策略3、检查第h个零件，如为废品，则换刀；如为正品，则检查第h+1 个，若正品，继续生产，如废品，换刀。 策略4、检查第h和h+1个零件，如两件都为正品，则继续生产；如两件都 为废品，换刀；如一正一费，则检查第h+2个，若正品，继续生 产，如废品，换刀。 • 至多到T个产品一定换刀。 通过计算机模拟比较不同策略下的单位正品最小费用，对不同策略进行选择。且可利用模拟进行参数灵敏度分析。

45. 结论：策略IV较优

46. 倒媒台的操作方案 • 某煤矿公司有一个大型煤台，用于向运媒列车装煤。倒煤台的容量是1.5列标准列车。装满一个空的倒煤台需要一个小组 6 个小时的时间，费用是 9000元/小时。为提高装煤速度可以以 12000元/小时的代价动用第二个小组。铁道部门每天向这个倒煤台发三列空的标准车。这些列车可在上午5 点到下午8点之间的任何时刻到达。给一列标准车装满煤需要3小时，向倒煤台装煤和从倒煤台向列车装煤不能同时进行。如果列车到达后因等待装煤二停滞，铁道部门将征收每车1500元/小时的滞费。 • 此外，每星期四上午5 点到下午8点之间还有一列大容量列车到达，其容量为标准的列车的2倍，滞期费为 25000元/小时。 • 问：（1）安标准规则操作可使装煤费最低？费用使多少？ （2）如果标准列车能在指定的时间到达，什么样的调度安排最经济 ？ 可在方案的优化程度和简明性之间做一个折中，加入一些规则获得接近最优的解