【分毫析厘】确定性仿真结果检验和因子分析方法

4月前浏览5239

作者：傅惠民, 邬文娟北京航空航天大学

目前仿真技术已在航空、航天、船舶、电子、通信、材料、机械、建筑、生物工程、经济、金融等各个领域发挥着重要作用, 例如, 核仿真、结构强度仿真、振动仿真、可靠性仿真、经济仿真等.但是仿真结果的正确性与可信性问题一直困扰着人们。

为此,文献建立了随机仿真结果统计检验和因子分析方法, 解决了随机仿真结果小样本检验问题。鉴于确定性仿真(不考虑随机因素影响的仿真)在工程中的广泛应用, 本文进一步建立一种确定性仿真结果检验方法。虽然确定性仿真未考虑随机因素的影响, 仿真结果并不存在随机性, 但由于试验结果存在分散性, 不可能与仿真结果完全吻合, 那么仿真结果与试验结果之间的误差是由于仿真不正确所引起的系统误差, 还是由于试验结果的随机性引起的偶然误差? 目前将试验结果与仿真结果进行简单比对的方法不能判断该误差是系统误差还是偶然误差, 也就是说, 不能科学合理地判断仿真结果正确与否。若分别在每一状态对仿真结果进行假设检验, 则每一状态都要做大量试验, 使总的试验量很大, 而且每一状态下都存在弃真概率, 这将导致所有状态的总的显著性水平无法控制。

本文方法则能够将各个状态的仿真结果与试验数据作为一个整体进行统计推断, 对其误差进行小样本检验, 确定该误差是系统误差还是偶然误差, 从而解决确定性仿真结果的正确性与可信性问题。

在仿真过程中需要对众多复杂的影响因素(如温度、应力等)进行模拟。当仿真结果不正确时, 如何判断哪些影响因素已被正确仿真, 哪些影响因素尚未被正确模拟, 这是当前亟待解决的难题。为此, 本文建立一种确定性仿真结果多因子分析方法, 可以通过少量试验, 分析各个因素对仿真结果的影响, 确定哪些影响因素已被正确仿真, 哪些尚未被正确模拟, 从而指导仿真软件的编制。

多状态确定性仿真结果检验方法

1 .1 　基本模型

　设μ′i 为第i 个状态下真值μi 的仿真结果, i=1 , 2 , … ,m。当该仿真结果正确时, 有

　为了检验仿真结果μ′_i是否满足式(1), 在第i个状态下进行了n_i 个试样的独立试验, 试验值记为X _{i j} , j =1 ,2 , …, n_i.通常试验结果存在分散性, 并遵循正态分布, 即

则当仿真结果μ′_i正确, 即式(1)成立时, 有

因此, 通过检验Y _{i j} 是否来自于均值为零的正态母体, 可以判断仿真得到的μ′_i正确与否。

1 .2 　仿真结果检验方法

　　对于给定的显著性水平α, 若

则可以断定仿真结果不正确。式中

1 .3 　方差估计

　　方差σ²的无偏估计量 σ²由下式给出

1 .4 　仿真检验算例

　　设某产品性能在第i 个状态下的真值为μ_i ,其仿真结果为μ′_i , i =1 , 2 , …, 5 。由于在工程实际中真值μ_i是未知的, 因此无法将μ′_i与μ_i进行比较, 只能将μ′_i与试验结果比较。本文方法能够很好地判断仿真结果μ′_i和试验结果之间的误差是由于仿真不正确引起的系统误差, 还是由于试验结果的随机性引起的偶然误差.下面给出一个算例。

表1第2列和第3列分别给出了真值μ_i和仿真结果μ′_i , i =1 , 2 , …, 5 .第4 列给出了各个状态的模拟试验数据(标准差σ=5)。现采用本文方法对表1 第3列和第4列数据进行分析处理, 结果列于表2 。从中可以看到, 在显著性水平α=0.05下, 经本文方法检验可以得出表1 第3列仿真结果是正确的结论。由于该算例采用了MonteCarlo 模拟试验数据, 真值μ_i已知, 所以通过表1第2 列真值μ_i和第3列仿真结果μ′_i的直接对比, 可知本文方法的检验结论是正确的。从表1 还可知, 第5个状态的仿真结果μ′₅ =75 与试验结果86.2 之间的相对误差为14 .9 %, 如按照简单比对的方法,就有可能得出该仿真结果不正确的错误结论。

若表1 第2 列真值μ_i的仿真结果为μ′_i =1 .1μ_i , i =1 , 2 , … , 5 , 见表3。显然, 此时的仿真结果μ′_i是不正确的.现采用本文方法对表1 第4 列模拟试验数据和表3 仿真结果进行分析处理, 结果列于表4 。从中可以看到, 在显著性水平α=0 .05下, 经本文方法检验, 可以得出表3 仿真结果不正确的结论.从表1 第4 列和表3 第6 列还可知, 此时第5 个状态的仿真结果μ′₅ =82 .5 与试验结果86.2 之间的相对误差仅为4 .5 %, 如按照简单比对的方法,就有可能得出该仿真结果正确的错误结论。

确定性仿真结果多因子分析方法

仿真结果因子分析旨在对仿真中遇到的众多影响因素(如温度、应力等)进行分析,以确定哪些影响因素已被正确仿真, 而哪些尚未被正确模拟,从而指导仿真软件的编制。对于单一因子的情况,只要把上节多状态仿真结果检验方法中的状态换成因子, 即可得到仿真结果的单因子分析方法。下面将给出多因子分析方法。

2 .1 　基本模型

　　对于多因子试验, 设试验总次数为n , 试验值为X_i , 且有X_i ～ N (μ_i, σ²), i =1 , 2 , …, n .设μ′_i为真值μ_i的仿真结果。当μ_i的各个影响因素被正确仿真时, 有

则有

，因此, 通过检验Y_i 是否来自于均值为零的正态母体, 可以判断真值μ_i的各个影响因素是否被正确仿真。

2 .2 　多因子分析

2 .2 .1 　因子及其交互效应分析

　　对于给定的显著性水平α, 若

则可以断定仿真结果μ′_i(i =1 , 2 , …, n)未正确考虑因子A 对真值μ_i的影响。

同样, 在显著性水平α下, 当

时, 可以断定仿真结果μ′_i未正确考虑交互作用A×B 对真值μ_i的影响(i =1 , 2, …, n).式中SE为误差平方和, SA 为因子A 的效应平方和, S_A_×B为因子A 和因子B 的交互作用效应平方和, f_E 、f _A 和f_A_×B分别为其自由度.下面给出它们的计算方法。

2 .2.2 　总离差平方和

　　设总离差平方和为

式中QT 为所有试验数据的平方和, T 为所有试验数据的总和, Y 为总平均值, 即

可以证明, 当μ′_i(i =1 , 2 , …, n)满足式(10)时，有

2 .2 .3 　因子效应平方和

设n_{A i} 为因子A 的第i(i =1 , 2 , … ,mA)个水平下的试验次数, X_{A i k}为其第k 个试验结果,Y _{A i k}为其经式(11)变换得到的数据.则因子A 的效应平方和由下式计算

式中T_Ai 为因子A 第i 个水平下的所有试验数据的总和, Y_Ai为其平均值, 即

可以证明, 当μ′i(i =1 , 2 , …,n)满足式(10)时, 有

2 .2 .4 　两因子交互效应平方和

　　因子A 和因子B 的交互作用效应平方和S_A_×B由下式计算

其中T_ABij为因子A 第i 个水平和因子B第j 个水平组合(A_i , B _j )下的所有试验数据的总和,Y_ABij 为其平均值, n_ABij为两因子水平组合(A_i ,B_j )下的试验次数, Y _ABijk为其第k 个试验结果X_ABijk经式(11)变换得到的数据i =1 ,2, … ,m_A;j =1,2 , …,m_B ;k =1 ,2 , … , n_ABij .S_A_×B的自由度等于因子A 和因子B 的自由度f _A和f_B 乘积

可以证明, 当μ′_i(i =1 , 2 , …,n)满足式(10)时, 有

2 .2 .5 　采用正交表时的两因子交互效应平方和

若采用正交表安排多因子有交互作用的试验,其中因子A 和因子B的交互作用在正交表上占了若干列,那么该交互作用效应平方和S_A_×B等于正交表上这若干列的效应平方和之和,其中每列的效应平方和可将每列看成单一因子按式(20)计算。此外也可按式(25)直接计算S_A_×B。

2 .2.6 　误差平方和　

　误差平方和为

式中S_F为所有因子的效应平方和之和, S _I为所有交互作用的效应平方和之和.S _E的自由度为

式中f_F 为所有因子平方和的自由度之和, f_I为所有交互作用效应平方和的自由度之和.可以证明

2 .3 　均值检验

当对所有影响因子以及它们两两之间的交互作用,式(12)和(13)都不成立时, 还需要E(Y)=0 , 仿真结果μ′_i (i =1 , 2 , … , n)才可能满足式(10)。若

则在显著性水平α下, 可以断定仿真结果μ′_i(i =1 , 2 , …,n)不正确。此时,所有仿真结果μ′_i都比其真值μ_i(i=1 , 2 , … , n)增大了E(Y)(工程中可用 Y 近似).式中 σ由式(35)计算。

2 .4 　方差估计

方差σ²的无偏估计量 σ²由下式给出

2 .5 　因子分析算例

设某产品性能μ_ij受到因子A (如温度)和因子B(如载荷)的影响, 且μ_ij与因子水平A _i和B _j之间存在如下关系

式中a0 =100 , a1 =2 , a2 =5.现对产品性能μ_ij进行仿真, 得到仿真结果为

式中a′0 =100 , a′₁ =2 , a′₂=5 .由于在工程实际中式(36)是未知的, 所以不能通过式(36)与式(37)之间的比较来确定仿真结果正确与否。下面采用本文多因子分析方法对仿真结果式(37)进行检验。考虑到本例中因子A 和B各有三个水平, 见表5 ,为合理安排试验, 可采用正交表L9(3⁴)进行正交试验设计, 将A ,B 这两个因子安排在表的前两列上, 按照L9(3⁴)的试验号进行试验。为了便于比较, 采用Monte Carlo 模拟方法,在水平组合(A_i ,B _j )下以真值μ_ij为均值, 并取标准差σ=15,得到模拟试验数据X_ij 列于表6,i=1,2 , 3 ;j=1 , 2 , 3 .由式(36)和式(37)给出的真值μ_ij和仿真结果μ′_ij也一并列于表6。

现采用本文方法对表6 第5 列和第6 列数据进行分析处理,结果列于表7 。从中可以看到, 在显著性水平α=0 .05 下, 经本文方法因子分析可以得出式(37)已正确考虑因子A和B的影响, 并且表6 第5 列仿真结果是正确的结论。通过式(36)与式(37)直接对比, 可知本文方法的分析结论是正确的。

若式(37)中仿真结果a′₀ =100, a′₁ =2 , a′₂ =5 。5 ,则仿真结果μ′_ij由表8给出,显然,此时式(37)未正确考虑因子B的影响,仿真结果μ′_ij是不正确的。现采用本文方法对表6第6列模拟试验数据和表8仿真结果进行分析处理, 结果列于表9 。从中可以看到, 在显著性水平α=0。05下,经本文方法因子分析,可知式(37)对因子A 的仿真正确, 但对因子B 的仿真不正确, 从而导致表8仿真结果也不正确。