建模与仿真

第八章 系统仿真结果分析

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

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

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

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

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

§8.1终态仿真的结果分析

8.1.1 重复运行法

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

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

n

n

j 1

Xj tn 1,1

S(n)/n

2

(8.1)

2

建模与仿真

其中， S2 n

n

X(n) X / n 1

j

j 1

n

2

(8.2)

X

1n

j 1

Xj

(8.3)

α为置信水平。

根据中心极限定理，若产生的样本点Xj越多，即仿真运行的次数越多，则Xj越接近于正态分布，因此在终态仿真中使用仿真方法运行的重复次数n不能选取得太小。

8.1.2序贯程序法

在终态仿真结果分析得重复运行法中，通过规定次数得仿真 可以得到随机变量取值的置信区间，置信区间的长度与仿真次数的平方根成反比。显然，若要缩小置信区间的长度就必然增加仿真次数n。这样就产生了另一个方面的问题，即在一定的精度要求下，规定仿真结果的置信区间，无法确定能够达到精度要求的仿真次数。这样就可以对置信区间的长度进行控制，避免得出不适用的结论。

一般说来，在同样精度要求下，采用序贯程序法得出的仿真重复运行次数比利用解析法得到的次数要少。

由式(8.1)可知，样本X的100(1-α)%置信区间的半长为：

tn 1, /2

X

(8.4) (8.5)

式中 X S/n

S为样本的标准差，n为重复运行次数。设给定一准确的临界值ε，即限定置信区间的长度为 X ,X ，并给定精度(1-α)。为了达到此精度要求，需要取足够大的仿真运行次数n，使之满足：

建模与仿真

P(X X ) 1

(8.6)

假设仿真已经重复运行了n0次(n0 2)，为了满足置信区间半长的临界值，必须选择重复运行次数n，使得：

n n0

且

(8.7) (8.8)

tn 1, /2 S0

n

初始运行仿真运行的次数应当至少大于2，最好取4或5。由式8.8可以推出n应当满足

n

tn 1, /2 S0

2

(8.9)

显然n的解就是满足式8.9的最小整数。

n min i:i

tn 1, /2 S0

2

(8.10)

注意这里假定n次独立重复运行结果总体方差σ2的估计值S2(n)随着增加n次运行没有显著的变化，因此可以用n0的总体方差代替。

实际上，利用n0次仿真运行的方差S2(n0)来替代n次仿真运行的方差，会使得计算得出的n值偏大。为了消除这种影响，一般采用序贯程序法，其步骤为：

1）预定独立仿真运行的初始次数n0

2，置

n=n0，独立运行n次；

2）计算该n次运行的样本X1,X2, ,Xn以及相应的S2(n)； 3）利用下式计算 值

tn 1, /2

S(n)n

2

如果

，则得到置信度为

1

的满足精度要求的置信区间

X(n) ,X(n) ，从而确定了相应的仿真次数n；

建模与仿真

4）否则令n=n+1，进行仿真得到样本值Xn 1； 5）返回步骤2）。

8.2稳态仿真的结果分析

研究系统的稳态性能，需要研究一次运行时间很长的仿真。在仿真运行过程中，每隔一段时间即可获得一个观测值Yi，从而可以得到一组自相关时间序列的采样值Y1,Y2, ,Yn，其稳态平均值定义为：

lim

n

1n

n

Yi

i 1

(8.11)

如果ν的极值存在，则ν与仿真的初始条件无关。

8.2.1批均值法

批均值法的基本思想是：设仿真运行时间足够长，可以得到足够多的观测值Y1,Y2, ,Ym，将Yi i 1,2, ,m 分为n批，每一批中有l个观测值，则每批观测数据如下：

第一批：Y1,Y2, ,Yl 第二批：Yl 1,Yl 2, ,Y2l

第n批：Y(n 1)l 1,Y(n 1)l 2, ,Ynl

首先对每批数据进行处理，分别得出每批数据的均值

l

Y

j

1

l Y(j 1)l k

k 1

(8.13)

建模与仿真

由此可得总得样本均值为：

n

Y

1n Y

j 1

j

1m

Yi

i 1

m

(8.14)

此即ν的点估计。

为了构造ν的置信区间，需要假定Yj是独立的且服从正态分布的随机变量，并具有相同的均值和方差。此时ν的近似置信区间的计算公式为：

Y tn 1,1 /2

Sj(n)n

2

2

(8.15) (8.16)

式中 S(n)

2j

n1n 1

(Yj

j 1

Y)

n为观测值的批数。

8.2.2稳态序贯法

在利用批均值法进行计算时，假定每批观测值的均值是独立的，但实际上

Y1,Y2, ,Yn

是相关的。为了得到不相关的Yj，直观的做法是：保持批数n不变，

不断增大l，直到满足不相关的条件为止。

但是如果n选择过小，则Yj的方差加大，结果得到的置信区间就会偏大，为此n也必须足够大。这样为了达到精度要求就必须选择足够大的n和l，使得样本总量m素。

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

设仿真运行观测值的批长度为l，已经有观测值 n批（ 为i的两批观测值批均值的相关系数

i(l) Cov[Yj,Yj 1],(j 1,2, ,n 1)

2

n l

特别大，而仿真过程中时间的消耗也是必须考虑的重要因

），考察相隔

建模与仿真

i(l)随

l的变化规律大致有三种情况：

1） i(l)为递减函数（见图8.1）；

2） i(l)的值一次或多次改变方向，然后严格地减少到0（见图8.2）； 3） i(l)<0或者随着l变化没有一定的规律。

图8.1 i(l)为单调递减函数

l

l

图8.2 i(l)多次改变方向然后递减

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

2） 进行长度为mi的仿真运行，获得mi个观测值Y1,Y2, ,Ym。

i

3） 令l mi4） 如果

/(n f)

，计算Yk(k

1,2, ,nf)、 j(nf,l)(可以取j 1)

。

2mi 1，

j(nf,l) u

，则说明mi太小，需加大，可以令i=i+1，且mi

返回第2步获取其余mi 1个观测值。

5） 如果 第8步。

6） 如果0

j(l)

j(nf,l) 0

，则表明增长仿真运行长度无助于 j(l)的判断，执行

j(nf,l) u，计算Yk(2l),(k 1,2, ,nf/2)、 j(nf/2,2l)(j 1)

，判断

是否具有第2类特征；如果

则说明该

j(l)j(nf/2,2l) j(nf,l)，

确实具有

建模与仿真

第2类特征，需要进一步加大mi，令i=i+1，且mi余mi 1个观测值。

7） 如果

j(nf/2,2l) j(nf,l)，则说明 j(l)

2mi 1，返回第

2步获取其

已经具有第1类特征，而且达的值满足独立性要求，此时

到 j(l)判断值n的l已经得到，可以相信

j(n,fl)

用批均值法计算该n批长度为fl的置信区间。

8） 计算Yk(n,

fl),Y(n,fl)

以及置信区间的半长

Y(n,fl)

tn 1,1 /2

S

2

n

，最后得

2mi 1，返回第

9） 如果 ，说明精度不满足要求，令i=i+1，且mi获取其余mi 1个观测值。

10） 止。

如果 ，则精度满足要求，可以令估计值

2步

Y(n,fl)

，仿真停

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

8.2.3再生法

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

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

在仿真过程中，随着仿真时钟的推进，系统的状态变量在不断地发生变化。

建模与仿真

如果在某一时刻观测到了系统一组状态变量的数值，而在其后的若干时间之后又重新观测到系统的完全相同的一组状态变量的数值，则称所观测到的系统为再生系统。也就是说，在稳态仿真中，系统从某一初始状态开始运行，若干时间后重新达到该状态；这时可以认为系统重新达到该状态后的过程相对于以前的过程是独立的，这就相当于系统在此时重新运行。显然在若干时间后这种情况将重新发生，因此这个重复的过程称为系统的再生周期，而系统初始状态重复出现的时刻点称为系统的再生点。。

再生法的缺点在于系统再生点的数量要求足够多，而且每个再生周期应该是独立的。而实际系统的仿真运行中可能不存在再生点或者再生周期过长，这样就要求仿真运行的总长度要足够大。

假设在M/M/1系统的观测中有p个完整的再生周期，令Yj为第j个再生周期中各个实体等待时间的总和：

nj

Yj

n

kj

k 1

j

为第j个再生周期中受到服务的实体个数。{Yj}和{nj}都是独立同分布

的随机序列，然而Yj和nj并不相互独立，因为较大的Yj值可指望有较大的nj值伴随产生。

假设总观测次数为N，各个实体的等待时间分别为 1, 2, , N，则实体的平均等待时间的估计值由下式给出：

N

W

1

N

i

i 1

如果将各个实体等待时间根据再生周期进行分组，则上式又可以写为：

W

Y1 Y2 YNn1 n2 nN

Yn

建模与仿真

式中：

p

Y

p

Y

j 1p

j

j

n

1p

n

j 1

Y

是一个再生周期中实体等待时间综合的估计值，n是一个再生周期中受到

服务的实体个数的估计值。当p足够大时，W是渐近无偏的，即：

lim

P

E(W) E(W)

而实际上，W对W的估计值是有偏的，因而需要估计统计值W的方差，以确定平均等待时间的置信区间，由于Yj和nj皆为随机变量，为了避免直接处理随机变量之比，引入变量Vj：

Vj Yj E(W)nj

这是一个独立同分布的随机变量序列，同时我们可以得到：

E(Vj) E(Yj) E(W)E(nj) 0

设 2为随机变量Vj的方差，根据中心极限定理，当p 时，下列随机变量：

Z

V

p

为收敛于标准正态分布的随机变量。式中

p

p

j

p

j

V

1

p

V

j 1

1p

Y

j 1

p

E(W)

1p

n

j 1

j

Y E(W)n

(8.17) (8.18)

从而有 P( Z1 /2

V

Z1 /2) 1

建模与仿真

式中Z1 /2为对应显著水平为 的标准正态分布的临界限。 将式（8.17）代入式（8.18），可以得出：

P(

Z1

p

/2

Y E(W)

Z1

n

/2

Z1

p

/2

) 1

即 P(W

Z1

n

/2

p

E(W) W

p

) 1

Z1

n

/2

从而得到平均等待时间的近似100(1 )%置信区间为： W

p

§8.3 方差减小技术

8.3.1公用随机数法(CRN)

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

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

考虑两个模型，设X1j和X2j分别是从第1个模型和第2个模型的仿真运行中得到的第j个独立再生周期中的数据，对 对每个模型产生了n个再生周期，并且设定

Zj X1j X2j,j 1,2, ,nn

E(X1j) E(X2j)

进行估计。如果

则E(Z

j)

，而Z(n)

1n

Zj

j 1

是 的一个无偏估计。由于Zj是独立同分布的随机变量，因此，我们可以得出：

建模与仿真

D[Z(n)]

D(Zj)

n

D(X1j) D(X2j) 2Cov(X1j,X2j)

n

如果两个模型得运行是独立得，则X1j和X2j是独立的，即Cov(X1j,X2j) 0；而如果能够使得X1j和X2j是正相关的，即使Cov(X1j,X2j) 0，这样得到的估计Z(n)的方差就减小了。

为了实施公用随机数法，需要使各个模型中的随机数同步，即在一个模型中使用于一个具体目的的随机数，在所有其它模型中也应该使用于同一目的，在仿真中达到这种同步的一般原则为：

1） 如果能够有几个可以同时工作的不同随机数发生器，则可以用一个发生器专门为一个指定的随机变量产生种子。不同的随机变量用不同的随机数发生器；

2） 实现产生出所需要的随机数并存储起来，在对各个模型仿真运行时按照需要取用这些随机数；

3） 使用逆变换法产生随机变量，因为这种方法每产生一个随机变量仅仅只需要一个单位均匀分布的随机数。

8.4.2对偶变量法(AV)

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

对偶变量法的中心思想就是在系统模型的两次仿真运行过程中，设法使得第1次运行中的小观测值能够被第2次仿真运行中的大观测值所补偿，或者是反过来。这就相当于采用两次运行中观测值的平均值作为分析的基准数据点，而这个平均值与所估计的观测值的期望更加接近。

建模与仿真

一般情况下，对偶变量法使用互补的随机数驱动系统模型的两次运行。也就是说，如果Uk是用于第1次运行中某一具体目的的单位均匀分布随机数，则在第2次运行中将（1-Uk）用于同一目的。

考察系统模型所进行的两次仿真运行，设定每次运行产生n个再生周期，

(1)(2)(1)(2)

这样可以构成一系列观测值对：(X1,X), ,(Xn,Xn)。各观测值对相互独1

立，令

Xj

12(X

(1)(2)

,X),j 1,2, ,njj1n

而 为

X(n)

Xj

j 1

n

E(Xj)的点估计，由于Xj

是独立同分布的随机变量，因此有

D(X

(1)(2)(1)(2)

) D(X) 2Cov(X,X)jjjj

n

D[X(n)]

D(Xj)n

如果两次运行是相互独立的，则Cov(X(j1),X(j2)) 0。如果能设法使得

X

(1)(2)(1)(2)

和XCov(X,X) 0之间形成负相关，也就是使jjjj

，方差便会减小。

8.4.3控制变量法

控制变量法是利用随机变量之间的相关性来实现方差衰减的目的。 设Y是某一个输出随机变量，我们要估计 已知，且与Y相关的随机变量。令：

Yc Y (X )

E(Y)

，而X是一个均值

E(X)

其中： 为任意实数，易见Yc也是 的无偏估计。由于

Var(Yc) E Y (X ) Var(Y) Var(X) 2 Cov(X,Y)

2

2

建模与仿真

如果 2Var(X) 2 Cov(X,Y) 0，则Yc的方差 将小于Y的方差，称X为Y的控制变量，因为X部分地控制了Y。

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

优良的控制变量X应当与Y是强相关的，这取决于模型本身的结构。例如，在排队系统中取输出变量Y为顾客的平均等待时间，则控制变量可以选用已知分布的输入变量。当选用到达时间间隔时间为控制变量，则它与Y之间构成负相关关系。如果选用服务时间为控制变量，则它与Y之间构成正相关关系。利用这种控制变量的优点是它们在仿真中必须产生，因此不必另外花时间去专门产生这种随机变量。

建模与仿真

1、 对于杂货店仿真模型，在相同初始条件下，经过4次独立运行，给出各

次顾客平均逗留和平均服务时间，并求它们的95%置信区间。（重复运行法）

2、 对于杂货店仿真模型，从上午9：00开始到下午4：00进行仿真运行，

以两个小时为分批区间，给出各批顾客平均逗留和平均服务时间，并求它们的95%置信区间。（批均值法）

3、 从基于Agent的建模与仿真、分布交互仿真、基于Petri网的建模与仿真

中选择一个较小的主题，写一篇3000字的文献综述，列出参考文献。