8第八章-系统仿真结果分析

8第八章-系统仿真结果分析
第八章 系统仿真结果分析
采用统计方法来估计系统的性能，利用统计分析方法要求样本数据具有统计独立性，但实际上在很多情况下这个条件并不能满足。

解决这一难题的途径无非两条：一是对样本序列进行处理，使之尽量满足统计独立性条件；二是在经典统计方法的基础上进行修正使之适合于处理相关的样本序列。

终态仿真是指仿真实验在某个持续事件段上运行。

稳态仿真则是通过系统的仿真实验，希望的得到一些系统性能测度指标在系统达到稳态时的估计值。

有必要采用方差减小技术，即在相同的仿真运行次数下获得较小方差的仿真输出结果。

§8.1终态仿真的结果分析
8.1.1 重复运行法
所谓重复运行方法是指选用不同的独立随机数序列，采用相同的参数、初始条件以及用相同的采样次数n 对系统重复进行仿真运行。

对于一终态仿真的系统，由于每次运行是相互独立的，因此可以认为每次仿真运行结果()n i X i ,,2,1⋅⋅⋅=是独立同分布的随机变量，是服从正态分布的随机变量。

随机变X 量的期望值E (X )地估计值μ为：
S 为样本的标准差，n 为重复运行次数。

设给定一准确的临界值ε，即限定
置信区间的长度为[]
εε+-X X
,，并给定精度(1-α)。

为了达到此精度要求，需
要取足够大的仿真运行次数n ，使之满足：
αε-≥<-1)(X X P
(8.6)
假设仿真已经重复运行了n 0次(n 0≥2)，为了满足置信区间半长的临界值，必须选择重复运行次数n ，使得：
n ≥n 0
(8.7)
且 εβα≤⋅=
-n
S t n 0
2/,1
(8.8)
初始运行仿真运行的次数应当至少大于2，最好取4或5。

由式8.8可以推出n 应当满足
(
)2
02/,1ε
αS t n n ⋅-≥
(8.9)
显然n 的解就是满足式8.9的最小整数。

(
){
}2
02/,1: min ε
αS t n i i n ⋅-≥
=
(8.10)
注意这里假定n 次独立重复运行结果总体方差σ2的估计值S 2(n )随着增加n 次运行没有显著的变化，因此可以用n 0的总体方差代替。

实际上，利用0n 次仿真运行的方差)0(2n S 来替代n 次仿真运行的方差，会使得计算得出的n 值偏大。

为了消除这种影响，一般采用序贯程序法，其步骤为：
1） 预定独立仿真运行的初始次数20≥n ，置n=0n ，独立运行n 次； 2） 计算该n 次运行的样本n X X X ,,2,1 以及相应的)(2n S ；
3） 利用下式计算β值
n
n S n t )
(2/,12•
-=αβ
如果εβ≤，则得到置信度为α-1的满足精度要求的置信区间
[]ββ+-)(,)(n X n X ，从而确定了相应的仿真次数n ；
4） 否则令n=n+1，进行仿真得到样本值1+n X ； 5） 返回步骤2）。

8.2稳态仿真的结果分析
研究系统的稳态性能，需要研究一次运行时间很长的仿真。

在仿真运行过程中，每隔一段时间即可获得一个观测值i Y ，从而可以得到一组自相关时间序列的采样值n Y Y Y ,,,21 ，其稳态平均值定义为：
∑=∞→=n
i i Y n n 1
1lim
ν (8.11)
如果ν的极值存在，则ν与仿真的初始条件无关。

8.2.1批均值法
批均值法的基本思想是：设仿真运行时间足够长，可以得到足够多的观测值m Y Y Y ,,,21 ，将()m i Y i ,,2,1 =分为n 批，每一批中有l 个观测值，则每批观测数据如下：
第一批：l Y Y Y ,,,21
第二批：l l l Y Y Y 221,,, ++
第n 批：nl l n l n Y Y Y ,,,2)1(1)1( +-+-
首先对每批数据进行处理，分别得出每批数据的均值
∑=+-=
l
k k l j l
j Y Y 1
)1(1
(8.13)
由此可得总得样本均值为：
∑∑===
=
n
j m
i i m
j
n
Y Y Y 1
1
11
(8.14)
此即ν的点估计。

为了构造ν的置信区间，需要假定j Y 是独立的且服从正态分布的随机变量，并具有相同的均值和方差。

此时ν的近似置信区间的计算公式为：
n
n S n j t Y )(2
/1,12αν--±=
(8.15) 式中 2
1
1
12)
()(∑=--=n
j j n j
Y Y n S
(8.16)
n 为观测值的批数。

8.2.2稳态序贯法
在利用批均值法进行计算时，假定每批观测值的均值是独立的，但实际上n Y Y Y ,,2,1 是相关的。

为了得到不相关的j
Y ，直观的做法是：保持批数n 不变，
不断增大l ，直到满足不相关的条件为止。

但是如果n 选择过小，则j Y 的方差加大，结果得到的置信区间就会偏大，为此n 也必须足够大。

这样为了达到精度要求就必须选择足够大的n 和l ，使
得样本总量l n m ⨯=特别大，而仿真过程中时间的消耗也是必须考虑的重要因素。

稳态序贯法是一种尽可能减少m 的方法，较好地解决了批长度的确定以及仿真运行总长度的确定问题，并能满足规定的置信区间精度的要求。

设仿真运行观测值的批长度为l ，已经有观测值n •λ批（2≥λ），考察相隔为i 的两批观测值批均值的相关系数
)1,,2,1(],1,[)(-=+=n j j Y j Y Cov l i ρ
)(l i ρ随l 的变化规律大致有三种情况：
1） )(l i ρ为递减函数（见图8.1）；
2） )(l i ρ的值一次或多次改变方向，然后严格地减少到0（见图8.2）；
3） )(l i ρ<0或者随着l 变化没有一定的规律。

根据)(l i ρ的以上3种特性，基于批均值法的稳态序贯法原理如下：
1） 给定批数因子n 、f 以及仿真长度1m （1m 是f n •的整数倍），
)(l i ρ的判断值为u ，置信区间的相对精度γ，置信水平α。

令i =1。

2） 进行长度为i m 的仿真运行，获得i m 个观测值i m Y Y Y ,,2,1 。

l
l
)
(l i ρ
)
(l i ρ
图8.1 )(l i ρ为单调
图8.2 )(l i ρ多次改
3） 令)/(f n i m l •=，计算)1)(,(),,2,1(==j l nf j nf k k Y 可以取、ρ 。

4） 如果u l nf j ≥),(ρ，则说明i m 太小，需加大，可以令i =i +1，且12-=i m i m ，返回第2步获取其余1-i m 个观测值。

5） 如果0),(≤l nf j ρ，则表明增长仿真运行长度无助于)(l j ρ的判断，执行第8步。

6） 如
果
u
l nf j <<),(0ρ，计算)1)(2,2/()2/,,2,1(),2(==j l nf j nf k l k Y ρ、 ，判断)(l j ρ是否具有第
2类
特征；如果),()2,2/(l nf j l nf j ρρ≥，则说明该)(l j ρ确实具有第2类特征，需要进一步加大i m ，令i =i +1，且12-=i m i m ，返回第2步获取其余1-i m 个观测值。

7） 如果),()2,2/(l nf j l nf j ρρ<，则说明)(l j ρ已经具有第1类特征，而且达到)(l j ρ判断值n 的l 已经得到，可以相信),(fl n j ρ的值满足独立性要求，此时用批均值法计算该n 批长度为fl 的置信区间。

8） 计算),(),,(fl n Y fl n k Y 以及置信区间的半长
n S n t 2
2/1,1αδ--=，最后得
)
,(ˆfl n Y δ
γ=
9） 如果γγ>ˆ，说明精度不满足要求，令i =i +1，且12-=i m i m ，返回第2步获取其余1-i m 个观测值。

10） 如果γγ≤ˆ，则精度满足要求，可以令估计值δυ±=),(fl n Y ，
仿真停止。

稳态序贯法较好地解决了批长度的确定以及仿真运行总长度的确定问题，并能满足规定的置信区间精度的要求。

8.2.3再生法
在批均值法中，选取批长度的原则尚未完全确定，因此有必要考虑其它有效的方法。

再生法的思想就是要找出稳态仿真过程中系统的再生点，由每个再生点开始的再生周期中所获得的统计样本都是独立同分布的，可以采用经典统计方法对参数进行评估并构造参数值的置信区间
在仿真过程中，随着仿真时钟的推进，系统的状态变量在不断地发生变化。

如果在某一时刻观测到了系统一组状态变量的数值，而在其后的若干时间之后又重新观测到系统的完全相同的一组状态变量的数值，则称所观测到的系统为再生系统。

也就是说，在稳态仿真中，系统从某一初始状态开始运行，若干时间后重新达到该状态；这时可以认为系统重新达到该状态后的过程相对于以前的过程是独立的，这就相当于系统在此时重新运行。

显然在若干时间后这种情况将重新发生，因此这个重复的过程称为系统的再生周期，而系统初始状态重复出现的时刻点称为系统的再生点。

再生法的缺点在于系统再生点的数量要求足够多，而且每个再生周期应该是独立的。

而实际系统的仿真运行中可能不存在再生点或者再生周期过长，这样就要求仿真运行的总长度要足够大。

Y为第j个再生周假设在M/M/1系统的观测中有p个完整的再生周期，令
j
期中各个实体等待时间的总和：
∑==j
n k kj j Y 1
ϖ
j n 为第j 个再生周期中受到服务的实体个数。

{j Y }和{j n }都是独立同分布
的随机序列，然而j Y 和j n 并不相互独立，因为较大的j Y 值可指望有较大的j n 值伴随产生。

假设总观测次数为N ，各个实体的等待时间分别为N ϖϖϖ,,,21 ，则实体的平均等待时间的估计值由下式给出：
∑==
N
i i
N
W 1
1
ϖ
如果将各个实体等待时间根据再生周期进行分组，则上式又可以写为：
n
Y
n n n Y Y Y W N N =++++++=
2121
式中：
∑∑===
=p
j j
p
p
j j
p
n
n Y
Y 1
111
Y 是一个再生周期中实体等待时间综合的估计值，n 是一个再生周期中受
到服务的实体个数的估计值。

当p 足够大时，W 是渐近无偏的，即：
)()(lim W E W E P =∞
→
而实际上，W 对W 的估计值是有偏的，因而需要估计统计值W 的方差，以确定平均等待时间的置信区间，由于j Y 和j n 皆为随机变量，为了避免直接处理随机变量之比，引入变量j V ：
j j j n W E Y V )(-=
这是一个独立同分布的随机变量序列，同时我们可以得到：
0)()()()(=-=j j j n E W E Y E V E
设2σ为随机变量j V 的方差，根据中心极限定理，当∞→p 时，下列随机变量：
p V
Z σ=
为收敛于标准正态分布的随机变量。

式中
n W E Y n W E Y V
V p j p j j p j p p j j p )()(1
11111-=•-==∑∑∑=== (8.17) 从而有 ααασ-=≤≤---1)(2/12/1Z Z P p V (8.18)
式中2/1α-Z 为对应显著水平为α的标准正态分布的临界限。

将式（8.17）代入式（8.18），可以得出：
ααασσ-=≤-≤--
-
1))((2/12/1p Z p Z n W E Y P
即 ααασσ-=+≤≤---1))((2/12/1p n Z p n Z W W E W P
从而得到平均等待时间的近似)%1(100α-置信区间为： p n Z W 2/1ασ-±
§8.3 方差减小技术
8.3.1公用随机数法(CRN)
公用随机数法是应用于对两个或者几个不同的系统模型进行比较的情况。

采用公用随机数法的目的就是在其它环境条件完全相同的情况下，尽量消除因
为选取随机数造成的仿真运行结果的差异，而使得所观测到的差异仅仅只是来源于系统模型本身的差异。

公用随机数法的思想为：在不同模型的仿真运行过程中，采用相同的单位均匀分布种子随机数流。

考虑两个模型，设j X 1和j X 2分别是从第1个模型和第2个模型的仿真运行中得到的第j 个独立再生周期中的数据，对)2()1(j X E j X E -=α进行估计。

如果对每个模型产生了n 个再生周期，并且设定
n j j X j X j Z ,,2,1,21 =-=
则α=)(j Z E ，而∑==n j j Z n n Z 1
1)( 是α的一个无偏估计。

由于j Z 是独立同分布的随机变量，因此，我们可以得出：
n j X j X Cov j X D j X D n j Z D n Z D )2,1(2)2()1()()]([-+==
如果两个模型得运行是独立得，则j X 1和j X 2是独立的，即0)2,1(=j X j X Cov ；而如果能够使得j X 1和j X 2是正相关的，即使0)2,1(>j X j X Cov ，这样得到的估计)(n Z 的方差就减小了。

为了实施公用随机数法，需要使各个模型中的随机数同步，即在一个模型中使用于一个具体目的的随机数，在所有其它模型中也应该使用于同一目的，在仿真中达到这种同步的一般原则为：
1） 如果能够有几个可以同时工作的不同随机数发生器，则可以用
一个发生器专门为一个指定的随机变量产生种子。

不同的随机变量用不同的随机数发生器；
2） 实现产生出所需要的随机数并存储起来，在对各个模型仿真运
行时按照需要取用这些随机数；
3） 使用逆变换法产生随机变量，因为这种方法每产生一个随机变
量仅仅只需要一个单位均匀分布的随机数。

8.4.2对偶变量法(A V)
对偶变量法是一种应用于单个系统模型仿真运行时的方差减小技术。

对于同一个系统模型，每一次仿真运行中得到的观测数据时存在差异的，同样这种差异可能由随机数的选取而引起，采用对偶变量法的目的就是尽量消除这种差异。

对偶变量法的中心思想就是在系统模型的两次仿真运行过程中，设法使得第1次运行中的小观测值能够被第2次仿真运行中的大观测值所补偿，或者是反过来。

这就相当于采用两次运行中观测值的平均值作为分析的基准数据点，而这个平均值与所估计的观测值的期望更加接近。

一般情况下，对偶变量法使用互补的随机数驱动系统模型的两次运行。

也就是说，如果k U 是用于第1次运行中某一具体目的的单位均匀分布随机数，则在第2次运行中将（1-k U ）用于同一目的。

考察系统模型所进行的两次仿真运行，设定每次运行产生n 个再生周期，
这样可以构成一系列观测值对：))2(,)1((,),)2(1
,)1(1(n X n X X X 。

各观测值对相互独立，令
n j j X j X j X ,,2,1),)2(,)1((2
1 == 而 ∑==n j j X n n X 1
1)( 为)(j X E =μ的点估计，由于j X 是独立同分布的随机变量，因此有
n
j X j X Cov j X D j X D n j X D n X D ))2(,)1((2))2(())1(()()]([++== 如果两次运行是相互独立的，则0))2(,)1((=j X j X Cov 。

如果能设法使得
)2()1(j X j X 和之间形成负相关，也就是使0))2(,)1((<j X j X Cov ，方差便会减小。

8.4.3控制变量法
控制变量法是利用随机变量之间的相关性来实现方差衰减的目的。

设Y 是某一个输出随机变量，我们要估计)(Y E =μ，而X 是一个均值)(X E =υ已知，且与Y 相关的随机变量。

令：
)(υα--=X Y c Y
其中：α为任意实数，易见c Y 也是μ的无偏估计。

由于
[]),(2)()()()(22Y X Cov X Var Y Var X Y E Y Var c ααυα-+=--= 如果0),(2)(2<-Y X Cov X Var αα，则c Y 的方差 将小于Y 的方差，称X 为Y 的控制变量，因为X 部分地控制了Y 。

从上面地讨论可以看出，控制变量法地两个关键问题是寻找适当的控制变量和确定系数α，使)(c Y Var 尽可能小。

优良的控制变量X 应当与Y 是强相关的，这取决于模型本身的结构。

例如，在排队系统中取输出变量Y 为顾客的平均等待时间，则控制变量可以选用已知分布的输入变量。

当选用到达时间间隔时间为控制变量，则它与Y 之间构成负相关关系。

如果选用服务时间为控制变量，则它与Y 之间构成正相关关系。

利用这种控制变量的优点是它们在仿真中必须产生，因此不必另外花时间去专门产生这种随机变量。

1、对于杂货店仿真模型，在相同初始条件下，经过4次独立运行，给出各
次顾客平均逗留和平均服务时间，并求它们的95%置信区间。

（重复运行法）
2、对于杂货店仿真模型，从上午9：00开始到下午4：00进行仿真运行，
以两个小时为分批区间，给出各批顾客平均逗留和平均服务时间，并求它们的95%置信区间。

（批均值法）
3、从基于Agent的建模与仿真、分布交互仿真、基于Petri网的建模与仿
真中选择一个较小的主题，写一篇3000字的文献综述，列出参考文献。