系统工程第四章系统仿真

67第四章系统仿真4.1基本概念系统仿真也叫系统模拟，就是通过建立和运行系统的计算机仿真模型，来模仿实际系统的运行状态及其随时间变化的规律，以实现在计算机上对实际系统的结构和行为进行动态实验的全过程。通过对仿真过程运行过程的观察和统计，得到实际系统的仿真输出参数和基本特性，可以此来估计和推断实际系统的真实参数和真实性能，以便掌握实际系统运动变化的规律，找到最优的或满意的解决实际问题的办法。特别是比较复杂的系统很难建立数学表达式时，往往可以通过计算机仿真技术来进行研究。例如，在图4-1中，已知采区来煤量q1是一个随机变量，根据统计结果，它服从正态分布，已知其均值μ=1.45t/min，其标准差σ=0.095t/min。大巷是胶带输送机，其运输量是常数：q2=1.33t/min。求8h后，采区最大煤仓容量。象这种问题，就可通过计算机仿真方法来寻求解答。下面再举一例进一步说明计算机仿真的基本概念。早在18世纪就有人作了随机模拟实验，投掷硬币4040次，出现正面2043次，从而求得出现正面的频率为P=2043/4040=0.5080。后来有人投了24000次，出现正面教育2012次，求得出现正面的频率为0.5005，更接近0.5。上述投掷硬币实验是人工作的，很费时间和精力，实际上完全可由计算机来进行，这就要建立一个计算机能够运算的模拟模型。建立这个模型的关键是用随机数去模拟投掷硬币的实验。所谓随机数就是一组随机出现的数列，通常用的是[0，1]区间的随机数。这时要求[0，1]区间的随机数必须具备下列性质：一是它的均匀性，即这些随机数落在[0，1]区间内任一位置的概率是相同的；另一个是它的随机性，即这些数在[0，1]区间内任一位置出现与否是随机的、独立的。如以R表示随机数，则出现0≤R≤0.5和0.5≤R≤1的概率相等，都为0.5，这样就可以模拟硬币出现的情况。如以随机数R≤0.5模拟出现正面，以R≥0.5模拟出现反面，则模拟模型就建立起来了。至于用什么办法获得[0，1]区间的均匀随机数，后面将专门论述。设I是投掷次数，N为正面出现次数，R为随机数，P是正面出现频率，下面画出模拟模型框图，如图4-2所示。根据研究的系统对象的性质，系统仿真可分为q2→q1→煤仓图4-1采区煤仓系统示意图I24000?R0.50.5实验开始I=I+1产生随机数R记录出现正面次数N=N+1计算并输出结果N，P模拟结束初始状态I=0，N=0YNYN图4-2投掷硬币实验模拟框图68连续系统仿真和离散事件系统仿真。连续系统是指系统状态随时间连续变化的系统，系统行为通常是一些连续变化的过程。连续系统模型通常是用一组方程式描述，如微分方程、差分方程等。差分方程形式上是时间离散的，但变量的变化过程本质上是时间连续的，如人口的变化过程、导弹运动、化工过程等。因此，连续系统仿真的主要任务就是如何求解系统运动方程组。离散时间系统中，表征系统性能的状态只能在随机的时间点上发生跃变，且这种变化是由随机事件驱动的，在两个时间点之间，系统状态不发生任何变化。4.2蒙特卡罗仿真蒙特卡罗(Monte-Carlo)仿真是一种特殊的数值计算仿真方法,它是充分利用计算机计算能力的随机实验方法。下面通过一个单变量积分问题作简要介绍。设)(xf是区间（a,b）内的连续函数，且已知在该区间内0)(xf，要求计算badxxf)(的值。当)(xf的形式比较复杂且是非可积函数时，用解析法往往十分困难。但却可以构成一个十分简单的离散、随机性仿真模型。因0)(xf，令c为)(xf之上限，则函数)(xf所表示的曲线将包含在边长为c和(b-a)的矩形之内，如图4-3所示。若在（a，b）和（0，c）区间中分别产生两个均匀分布的随机数xi和yi，由此构成的坐标点Pi(xi，yi)必定在矩形abde内。为了计算该积分的值，可规定以下准则：若yi≤f(xi)，则接受该点并作计数统计；若yif(xi)，则舍弃该点不作计数统计。设共产生N个随机仿真点，如果共接受个m点，（即恰在曲线上或在曲线以下的点），当N充分大时有badxxfcabNm)()(曲线下面积（4-1）由上式可见，通过产生随机数的方法来统计落在曲线下的坐标点数目，即可计算出定积分的值，计算的精度将随N的增大而提高。这就是典型的蒙特卡罗方法。将蒙特卡罗方法划分为一种仿真方法，主要是因为许多仿真问题都包含了随机数应用的缘故。但系统仿真和蒙特卡罗方法是两种数值计算技术，前者主要是用于动态模拟计算，而后者主要用于静态模型的计算。求定积分badxxf)(本身是一个确定型的连续模型，但由于采用了随机实验的蒙特卡罗方法，使其变成了一个离散的随机性问题。由此可以看出，在仿真领域，仿真模型的类别有时并不是绝对的。连续系统往往可以通过离散技术建模，而确定型问题也可以建立成随机性模型。4.3均匀随机数的产生方法4.3.1均匀随机数所有随机事件模拟都离不开均匀随机数。所谓均匀随机数是指服从均匀分布的随机变量拒绝接受dxxbxjxia0cef(xi)f(xj)图4-3求定积分的仿真模型69的样本值。经常用到的[0，1]区间上的均匀随机数，它的密度函数f(x)和F(x)分别为：其他0101)(xxf（4-2）111000)(xxxxxF（4-3）均匀随机数的均匀性是指在长度相同的每一子区间出现的概率值是一个常数；它的随机性则表现在[0，1]区间内，出现在某一子区间是没有规律的、任意的。4.3.2均匀随机数的产生方法产生随机数的方法很多，其中应用最广的是用公式产生随机均匀随机数。用数学公式产生的均匀随机数，不是客观随机现象的真实的随机数，所以人们称之为伪随机数。但产生这种随机数非常简单，并具有随机数的统计性质，因此得到了广泛应用。最常用的公式方法是线性同余法。先解释一下什么是同余数。如有两个数13和23，13÷10余3，23÷10也余3，则称13和23是以10为模的同余数。一般地，若有a、b两个整数被某个整数M除后，它们的余数相同，则称a、b两数对模M同余，记作a=b(ModM)。进一步，由于a=cM+t，b=dM+ta-b=(c-d)M则有整数dcMba所以也可以换一种说法：若整数a与b之差是自然数M的倍数，则称a、b对模数M同余。求余数的公式为MMyyc（4-4）这里My是对方括号中的数值取最大整数，例如1510505151010515515c产生随机数的线性同余法是一个递推公式：niMcaZZii,2,1))(mod(1（4-5）按照求余数的公式，上式可写为：MMcaZcaZZiii11（4-6）式中乘数a、增量c、模数M及随机数种子Z0都是非负整数，并且M0，Ma，Mc，MZ0。根据式（4-6），给出Z0可求出Z1，把Z1代入该式可求出Z2，再把Z2代入又可求出Z3…这样就可得到一个随机数列。更进一步，可通过下式得到落在[0，1]区间内的随机数γ：70Mxrn1（4-7）理论研究和实践表明，要使通过式（4-6）计算获得的数列具有较好的均匀随机性，M应尽量大，并选好a、c和Z0。两个较好的计算公式如下：)2)(mod15(35115iiZZ（4-8）)2)(mod453806246314159269(311iiZZ（4-9）4.3.3均匀随机数检验方法用公式产生的随机数还应经过检验来判断是否具备良好的均匀性和随机性（独立性）。检验随机数均匀性常用x2检验法。如果有N个随机数，把它们从小到大排列起来并分成互不重叠的等区间的k段，从理论上讲，这些随机数应该均匀分布在每一段上，即在每一段上的理论频数应该相等。理论频数为kikNmi,2,1,如果统计出这些随机数实际落在每段的频数为ni，则令kiiiimmnx122)(（4-10）若x2=0，说明理论频数与实际频数完全吻合，该随机数的均匀性很好。但实际情况并不为零，而是大于零的常数。x2值越大表明随机数列越不均匀。究竟多大才符合要求？有一个判定标准称α为显著度。可以根据自由度r=k-1，从x2表中查得值，如果计算所得的x2，就认为该随机数列在显著度α的情况下，其均匀性符合要求。例如，假设由公式（4-8）求得了100个随机数，将其分为5段（k=5）。找出该随机数列的最大值为31252，最小值为1488，并求出每段区间数值为5952.8，便可找出各段的范围。已知各段的理论频数mi=100/50=20，统计出各段的实际频数ni，最后求出x2=2。查x2表，根据自由度r=5-1=4，显著度α=0.05查分布临界值表，得49.9205.0x。由于205.02xx则认为在显著度0.05情况下，该随机数列的均匀性符合要求。随机数的独立性是指一个随机数的产生与它前一个或前几个随机数没有关联。检验方法之一是求出它们的相关系数。前后距离为j个数的相关系数为：202/)(1SZZZNjNijiij（4-11）式中，Z—随机数的均值，若为[0，1]区间的均匀随机数，其理论值为21；S2—随机数的方差，对[0，1]区间均匀分布随机数，有212)21(11NiiZNS（4-12）N—被检验的随机数目。2x2x71当N充分大时（N-j50），若取零假设H0：ρ=0，并取统计量jNuj(近似服从N（0，1）标准正态分布)，同时选定显著度α=0.05，则根据概率论，当|u|1.96时为差异显著，拒绝H0假设；当|u|≤1.96时接受假设，即检验的随机数符合独立性要求。4.3.4计算机函数通常，计算机高级语言有产生随机数的标准函数，可直接调用。如在有的BASIC语言版本中，就有标准函数RND（x）用于产生随机数。4.4随机变量的产生方法在生产实际中遇到的随机变量往往不服从均匀分布，而是其他各种形式的分布，如负指数分布、正态分布等。这些非均匀的随机变量，可以利用[0，1]区间的均匀数，通过一定的数学公式把它们表现出来。下面介绍随机变量几种常用产生方法。4.4.1逆变换法概率论中有重要结论：设随机变量Y的分布函数F（y）是连续函数，而X是在（0，1）上均匀分布的随机变量，令Z=F-1（X），则Z与Y有相同分布。事实上有),())(()}({})({)(1zFzFGzFXPzXFPZFz式中：)(xG是随机变量X的分布函数：xxxxxG1,110,0,0)(若已知Y的概率密度为，由Y=F-1（X）可得YdyyfYFX)()(对给出的（0，1）上均匀分布随机数ir，则具有给定分布的随机数iy，可由方程iyidyyfr)(中解出。已知某随机变量的分布函数F（x）（如图4-4），并知其定义域为[0，1]区间，这时可用[0，1]区间的均匀分布随机数R来表示F（x），即F（x）=R，由此可得到随机变量的样本值x=F-1(R)。F-1为逆变换的表达式。例如，已知负指数分布的密度函数为：000)(xxexfx其分布函数为01)(0xedtexFxxtxx0R1F(x)图4-4随机变量逆变换72令uReexFRxx11)(因R为[0，1]区间均匀分布随机数，1-R=u同样也是一随机数。上式两边取对数后有：则uxuxlnln4.4.2近似法很多随机变量都服从正态分布，但正态分布的函数形式比较复杂，用逆变换法很困难，一般常用近似法。正态分布的密度函数为xeaxfx,021)(222)(（4-12）分布函数为xtdtexF222)(21)(（4-13）其中μ为均值，σ为方差。当μ=0，σ=1时称标准正态分布，其密度函数和分布函数为：2221)(xexf（4-14）xtdtexF2221)(（4-15）从概率论可知，任何非标准的正态分布随机变量Y，通过下式都可以化为标