CN102033994A - 基于马尔可夫链蒙特卡罗的舵机可靠性仿真抽样方法 - Google Patents

Abstract

一种基于马尔可夫链蒙特卡罗的舵机可靠性仿真抽样方法，包括四个阶段，第一为马尔可夫过程模拟，有4个步骤：选取马尔可夫链的初始状态，确定随机转移抽样概率密度函数，确定马尔可夫链的下一个状态，不断重复，产生极限分布为渐进最优的重要抽样密度函数的随机样本点；第二为核密度估计，有3个步骤：选取核密度函数，确定窗口宽度参数和局部带宽因子，依据马尔可夫状态点，采用自适应宽核密度估计法，产生混合重要抽样概率密度函数；第三为重要抽样，根据第二阶段产生的混合重要抽样函数进行重要抽样；第四为统计计算，根据第三阶段产生的重要样本点，进行失效概率估计，并计算系统的失效概率。本方法有效解决了仿真效率低、精度低以及混合系统问题。

Description

基于马尔可夫链蒙特卡罗的舵机可靠性仿真抽样方法

技术领域

本发明提供一种基于马尔可夫链蒙特卡罗的舵机可靠性仿真抽样方法，它属于系统可靠性仿真分析领域的一种高效、高精度仿真方法，注重于解决含有离散变量的混合系统问题，如四余度舵机系统等。

背景技术

系统可靠性与性能一体化设计是一项在系统设计阶段利用故障及扰动注入、系统可靠性仿真分析和优化设计等方法来实现可靠性与性能综合分析与设计的新技术。实施可靠性与性能综合分析与设计可以实现在设计阶段进行性能设计的同时，得到相关的可靠性指标，为设计人员提供可靠性分析数据，为尽早发现系统的设计缺陷和优化设计方案提供一种有效手段。系统可靠性仿真分析作为一体化设计的关键环节之一，如何提高仿真效率、精度以及解决混合系统问题日益受到分析设计人员的重视。

在系统可靠性仿真分析中，蒙特卡罗法因为普适性强和简单易行被广泛应用，但由于其结果的准确性和收敛性均由大数定理保证，致使要获得某些小失效概率的系统可靠性结果时，需要庞大的仿真量作为支持，计算成本很高。为此发展了很多高效仿真算法，如对偶抽样法、条件期望抽样法、重要抽样法、分层抽样法、控制变数法和相关抽样法等。其中应用最多也最为有效的是重要抽样法，但传统的重要抽样法需要先根据极限状态方程确定抽样中心，对于很多复杂系统无法明确写出极限状态方程，从而很难求解抽样中心。1999年，Au S.K.首次提出了基于马尔可夫链蒙特卡罗法(MCMC)的重要抽样法，它将马尔可夫随机过程引入到蒙特卡罗模拟中，实现抽样分布随模拟的进行而改变，具有自适应的特点，不需要求解出极限状态方程。后来吕震宙等人就多失效模式、灵敏度求解等方向进行了算法的拓展改进。但其所研究的对象都是连续变量系统，并没有涉及存在故障不确定性的混合变量系统。在混合变量系统中，由于单元故障的发生，可能会导致系统构型的变化，从而影响重要抽样分布的稳定性。因此，针对此问题提出本发明方法。

发明内容

(1)目的：本发明的目的是提供一种基于马尔可夫链蒙特卡罗的可靠性仿真抽样方法，以提高仿真效率、精度和解决含有离散变量的混合变量系统仿真问题。

(2)技术方案

本发明一种基于马尔可夫链蒙特卡罗的舵机可靠性仿真抽样方法，该方法首先通过初始样本点在失效空间(因为关键故障而产生的不同失效域组成的集合)中游走来构造马尔可夫链模拟样本；然后综合考虑离散变量与连续变量，利用核密度估计构建混合核抽样密度函数；其次根据该密度函数进行重要抽样仿真；最后计算系统的失效概率和可靠度。该方法极大地提高了仿真效率和精度，并有效解决了离散变量与连续变量共存的混合变量系统的可靠性仿真分析问题，从而其在可靠性与性能一体化设计中具有更广泛的应用。

本发明一种基于马尔可夫链蒙特卡罗的舵机可靠性仿真抽样方法，该方法依次按照下述四个阶段进行；

1、马尔可夫过程模拟：通过马尔可夫随机过程模拟得到一些样本点(失效样本)；

2、核密度估计：应用上一步得到的样本点，对失效域进行预估计，拟和出混合核抽样密度函数；

3、重要抽样：应用上一步得到的结果作为重要抽样密度函数，进行重要抽样，得到高效样本点；

4、统计计算：利用上一步得到的样本点统计计算失效率和可靠度。

其详细步骤如下：

第一阶段：马尔可夫过程模拟

马尔可夫过程模拟主要包括以下四个步骤：

①选取马尔可夫链的初始状态X₀：

一般可依据工程经验或简单的数值方法确定失效域中的一点作为X₀。

②确定随机转移抽样概率密度函数：

定义混合型随机转移抽样概率密度函数为P^*(X_D ^(j))p^*(ξ|X_C ^(j))，混合概率密度函数P^*(X_D ^(j))p^*(ξ|X_C ^(j))与当前样本点X^(j)有关，用于定义X^(j)→X^(j+1)的产生过程；一般P^*(X_D ^(j))p^*(ξ|X_C ^(j))应具有对称性，在此可以选择较为简单的均匀分布，则：

P^*(X_D ^(j))p^*(ξ|X_C ^(j))

式(1)中S为离散变量的状态维数。l_i是n-1维以X_C ^(j)为中心的超立方体在第i维的边长，它决定了X_C ^(j+1)偏离X_C ^(j)的最大允许距离。在给定步长数目的情况下，l_i影响着马尔可夫链样本覆盖区域的大小。l_i越大，样本覆盖的区域也越大，但是过大的l_i将增加无效重复样本的数量。按照X_C ^(j+1)偏离X_C ^(j)的最大允许距离为三倍p^*(ξ|X_C ^(j)的标准差来确定l_i，则：

l_i＝6σ_iM^-1/(n+3)    (2)

式(2)中n为连续变量数目，M为马氏链步长即模拟样本数目，σ_i为h_opt(X_C)相对于X_C ^(j)的标准差的近似，可依经验确定。

③确定马尔可夫链的第j+1个状态

马尔可夫链的第j+1个状态X^(j+1)是在前一个状态X^(j)的基础上，由分布P^*(X_D ^(j))p^*(ξ|X_C ^(j))和Metropolis准则来确定的。基于X^(j)，依分布P^*(X_D ^(j))p^*(ξ|X_C ^(j))产生备选状态ξ，计算备选状态ξ和X^(j)状态的条件概率密度函数的比值r，即

r＝q(ξ)/q(X^(j))                 (3)

式(3)中q(X)＝I[X]P(X_D)f(X_C)。    (4)

然后根据Metropolis接受准则，确定马尔可夫转移的下一状态点：

$X^{(j+1)}=\left\{\begin{array}{ccc}\mathrm{\&xi;}& \mathrm{if}& r\&GreaterEqual;1\\ \mathrm{\&xi;}& \mathrm{if}& r<1,\mathrm{and\; r}\&GreaterEqual;\mathrm{random}\left[\mathrm{0,1}\right]\\ X^{\left(j\right)}& \mathrm{if}& r<1,\mathrm{and\; r}<\mathrm{random}\left[\mathrm{0,1}\right]\end{array}\right.---\left(5\right)$

④按上述步骤不断重复，产生极限分布为渐进最优的重要抽样密度函数的M个随机样本点{X⁽¹⁾，X⁽²⁾…，X^(M)}。

第二阶段：核密度估计

核密度估计主要包括以下三个步骤：

①选取核密度函数：

通常核函数属于对称的密度函数族P，从减小积分均方误差的角度来看，Silverman及Pracase Rao等指出P族中不同的核函数无明显差别，所以在此为计算方便，选择Gaussian密度函数，具体形式为：

$K\left(X_C\right)=\frac{1}{\sqrt{{\left(2\mathrm{\&pi;}\right)}^{n-1}\left|S\right|}}\exp (-\frac{1}{2}{X_C}^T{S}^{-1}{X}_C)---\left(6\right)$

式(6)中S为样本点集{X⁽¹⁾，X⁽²⁾…，X^(M)}的协方差阵，主要描述每个样本点在不同方向和范围上的数据分散性

$S=\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}({X_C}^{\left(j\right)}-{\stackrel{\&OverBar;}{X}}_C){({X_C}^{\left(j\right)}-{\stackrel{\&OverBar;}{X}}_C)}^T---\left(7\right)$

②确定窗口宽度参数和局部带宽因子：

窗口宽度过大，可能导致核函数对最优重要抽样函数的逼近过于光滑；窗口宽度过小，可能导致核函数对最优重要抽样函数的逼近产生不必要的噪声。此外若在整个样本区域内窗口宽度都是固定的，则当最优抽样函数的尾部较长时，核抽样函数可能会出现伪波动，影响整个计算的精度和效率。所以考虑在低(高)概率密度区域选用较大(小)的窗口宽度，即自适应宽核密度法的基本思想。具体分为：

局部带宽因子λ_j的求解：

${\mathrm{\&lambda;}}_j={\{{\left[\underset{k=1}{\overset{M}{\mathrm{\&Pi;}}}f\left(x^{\left(k\right)}\right)\right]}^{1/M}/f\left(x^{\left(j\right)}\right)\}}^{\mathrm{\&alpha;}}---\left(8\right)$

式(8)中0≤α≤1为灵敏因子，通常取α＝0.5。

窗口宽度参数w的求解：

$w={M_d}^{-\frac{1}{n+3}}---\left(9\right)$

式(9)中n为连续参数的个数，M_d为不同样本的个数(M_d≤M)。

③依据{X⁽¹⁾，X⁽²⁾…，X^(M)}，采用自适应宽核密度估计法，产生混合重要抽样概率密度函数k(X)：

$k\left(X\right)=P^h\left(X_D\right)\left[\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{1}{{\left(w{\mathrm{\&lambda;}}_j\right)}^{n-1}}K\left(\frac{X_C-{X_C}^{\left(j\right)}}{w{\mathrm{\&lambda;}}_j}\right)\right]---\left(10\right)$

在上式(10)中

$P^h\left(X_D\right)=\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{I}_{X_D}\left({X_D}^{\left(j\right)}\right),$ X_D＝1，2，…，S

(11)

在上式(11)中

第三阶段：重要抽样

从{1，2…，M}中均匀产生一个离散随机整数u，如果u＝j，则X_D＝X_D ^(j)，选取第j个分量的核抽样概率密度函数k_j(X_C)，来产生样本X_C，重复上述过程，直到得到N个样本点{X₍₁₎，X₍₂₎，…，X_(i)，…，X_(N)}。

$k_j\left(X_C\right)=\frac{1}{{\left(w{\mathrm{\&lambda;}}_j\right)}^{n-1}}K\left(\frac{X_C-{X_C}^{\left(j\right)}}{w{\mathrm{\&lambda;}}_j}\right)$

(13)

$=\frac{1}{{\left(\sqrt{2\mathrm{\&pi;}}w{\mathrm{\&lambda;}}_j\right)}^{n-2}}\frac{1}{\sqrt{2\mathrm{\&pi;}}w{\mathrm{\&lambda;}}_j\sqrt{S}}\exp (-\frac{(X_C-{X_C}^{\left(j\right)})^2}{2{\left(w{\mathrm{\&lambda;}}_j\right)}^{2}S})$

第四阶段：统计计算

根据重要抽样仿真样本点，进行失效概率估计：

${\hat{P}}_f=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}I\left[X_{\left(i\right)}\right]\frac{P\left(X_{D\left(i\right)}\right)f\left(X_{C\left(i\right)}\right)}{k\left(X_{\left(i\right)}\right)}---\left(14\right)$

由上式(14)得到系统可靠度：

$\hat{R}=1-{\hat{P}}_f---\left(15\right)$

综上所述，本发明所述方法共分为四大阶段，每一阶段又细分为几个步骤，共有9个步骤。马尔可夫过程模拟得到的样本点为下一阶段核密度估计进行失效域预估计之用；核密度估计得到的混合核抽样密度函数又是下一阶段重要抽样的抽样函数；通过重要抽样得到的重要区域样本点又是下一阶段统计计算所必须的。依次进行的四大阶段环环相扣，缺一不可。

(3)功效、优点

本发明所述的方法进一步完善了系统可靠性与性能一体化设计仿真技术。其功效主要在于以下三方面：

1.在相同的计算精度要求下，与传统蒙特卡罗法相比减少近1/15的计算量，计算效率明显提高。

2.在同样的计算量的要求下，与传统方法蒙特卡罗法相比，计算精度明显提高。

3.该方法可以直接处理离散变量与连续变量共存的混合变量系统，打破了传统基于马尔可夫蒙特卡罗的自适应重要抽样法的局限性，具有更广泛的适用性。

附图说明

图1本发明所述方法的流程框图

图2本发明所述方法的仿真结果图，即本发明方法与蒙特卡罗法失效概率变化对比曲线图

图3本发明所述方法的仿真结果图，即本发明方法与蒙特卡罗法方差缩减速度对比曲线

图中符号说明如下：

M马尔可夫链长度

N重要抽样样本点数

具体实施方式

本发明一种基于马尔可夫链蒙特卡罗的舵机可靠性仿真抽样方法，该方法依次按照下述四个阶段进行；其方法流程框图见图1所示，具体实施方式详述如下：

第一阶段：马尔可夫过程模拟

马尔可夫过程模拟主要包括以下四个步骤：

①选取马尔可夫链的初始状态X₀：

在四余度舵机系统中中，初始状态选择的是下表中各量的均值。

表1四余度舵机系统可靠性设计随机参数

电机和泵的转动惯量	～N(1.6E-3，2E-4)
		电机电枢电阻	～N(0.5，0.06)
电机电枢电感	～N(1E-2，9E-4)
		电机阻尼系数	～N(3E-4，2E-5)
泵的容积效率	～N(0.85，0.05)
		负载阻尼系数	～N(1000，50)
负载弹性系数	～N(5E6，1E5)
		泵排量	～N(1，0.08)

②确定随机转移抽样概率密度函数：

定义混合型随机转移抽样概率密度函数为P^*(X_D ^(j))p^*(ξ|X_C ^(j))。混合概率密度函数P^*(X_D ^(j))p^*(ξ|X_C ^(j))与当前样本点X^(j)有关，用于定义X^(j)→X^(j+1)的产生过程。一般P^*(X_D ^(j))p^*(ξ|X_C ^(j))应具有对称性，在此可以选择较为简单的均匀分布，则：

P^*(X_D ^(j)p^*(ξ|X_C ^(j))

式(16)中S选择为离散变量个数5。

l_i＝6σ_iM^-1/(n+3)    (17)

式(17)中n为连续变量的数目20，M为马氏链步长即模拟样本数目，在仿真中依次选择为100、200、300、400、500、600。

σ_i为h_opt(X_C)相对于X_C ^(j)的标准差的近似，选择为上述服从正态分布的连续变量的标准差(见表1)。

③确定马尔可夫链的第j+1个状态

马尔可夫链的第j+1个状态X^(j+1)是在前一个状态X^(j)的基础上，由分布P^*(X_D ^(j))p^*(ξ|X_C ^(j))和Metropolis准则来确定的。基于X^(j)，依分布P^*(X_D ^(j))p^*(ξ|X_C ^(j))产生备选状态ξ，计算备选状态ξ和X^(j)状态的条件概率密度函数的比值r，即

r＝q(ξ)/q(X^(j))                  (18)

式(18)中q(X)＝I[X]P(X_D)f(X_C)。    (19)

然后根据Metropolis接受准则，确定Markov转移的下一状态点：

$X^{(j+1)}=\left\{\begin{array}{ccc}\mathrm{\&xi;}& \mathrm{if}& r\&GreaterEqual;1\\ \mathrm{\&xi;}& \mathrm{if}& r<1,\mathrm{and\; r}\&GreaterEqual;\mathrm{random}\left[\mathrm{0,1}\right]\\ X^{\left(j\right)}& \mathrm{if}& r<1,\mathrm{and\; r}<\mathrm{random}\left[\mathrm{0,1}\right]\end{array}\right.---\left(20\right)$

上述各参数都是计算得来的。

④按上述步骤不断重复，产生极限分布为渐进最优的重要抽样密度函数的M个随机样本点{X⁽¹⁾，X⁽²⁾…，X^(M)}。

第二阶段：核密度估计

核密度估计主要包括以下三个步骤：

①选取核密度函数：

通常核函数属于对称的密度函数族P，从减小积分均方误差的角度来看，Silverman及Pracase Rao等指出P族中不同的核函数无明显差别，所以在此为计算方便，选择Gaussian密度函数，具体形式为：

$K\left(X_C\right)=\frac{1}{\sqrt{{\left(2\mathrm{\&pi;}\right)}^{n-1}\left|S\right|}}\exp (-\frac{1}{2}{X_C}^T{S}^{-1}{X}_C)---\left(21\right)$

式(21)中S为样本点集{X⁽¹⁾，X⁽²⁾…，X^(M)}的协方差阵，主要描述每个样本点在不同方向和范围上的数据分散性。

$S=\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}({X_C}^{\left(j\right)}-{\stackrel{\&OverBar;}{X}}_C){({X_C}^{\left(j\right)}-{\stackrel{\&OverBar;}{X}}_C)}^T---\left(22\right)$

上式的值都是计算得来的。

②确定窗口宽度参数和局部带宽因子：

局部带宽因子λ_j的求解：

${\mathrm{\&lambda;}}_j={\{{\left[\underset{k=1}{\overset{M}{\mathrm{\&Pi;}}}f\left(x^{\left(k\right)}\right)\right]}^{1/M}/f\left(x^{\left(j\right)}\right)\}}^{\mathrm{\&alpha;}}---\left(23\right)$

λ_j直接取为1。

窗口宽度参数w的求解：

$w={M_d}^{-\frac{1}{n+3}}---\left(24\right)$

式(24)中n为连续参数的个数20，M_d为不同样本的个数8。

③依据{X⁽¹⁾，X⁽²⁾…，X^(M)}，采用自适应宽核密度估计法，产生混合重要抽样概率密度函数k(X)：

$k\left(X\right)=P^h\left(X_D\right)\left[\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{1}{{\left(w{\mathrm{\&lambda;}}_j\right)}^{n-1}}K\left(\frac{X_C-{X_C}^{\left(j\right)}}{w{\mathrm{\&lambda;}}_j}\right)\right]---\left(25\right)$

在上式(25)中

$P^h\left(X_D\right)=\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{I}_{X_D}\left({X_D}^{\left(j\right)}\right),$ X_D＝1，2，…，S

(26)

在上式(26)中

第三阶段：重要抽样

从{1，2…，M}中均匀产生一个离散随机整数u，如果u＝j，则X_D＝X_D ^(j)，选取第j个分量的核抽样概率密度函数k_j(X_C)，来产生样本X_C，重复上述过程，直到得到N个样本点{X₍₁₎，X₍₂₎，…，X_(i)，…，X_(N)}，N分别取为1000和10000。

$k_j\left(X_C\right)=\frac{1}{{\left(w{\mathrm{\&lambda;}}_j\right)}^{n-1}}K\left(\frac{X_C-{X_C}^{\left(j\right)}}{w{\mathrm{\&lambda;}}_j}\right)$

(28)

$=\frac{1}{{\left(\sqrt{2\mathrm{\&pi;}}w{\mathrm{\&lambda;}}_j\right)}^{n-2}}\frac{1}{\sqrt{2\mathrm{\&pi;}}w{\mathrm{\&lambda;}}_j\sqrt{S}}\exp (-\frac{(X_C-{X_C}^{\left(j\right)})^2}{2{\left(w{\mathrm{\&lambda;}}_j\right)}^{2}S})$

第四阶段：统计计算

根据重要抽样仿真样本点，进行失效概率估计：

${\hat{P}}_f=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}I\left[X_{\left(i\right)}\right]\frac{P\left(X_{D\left(i\right)}\right)f\left(X_{C\left(i\right)}\right)}{k\left(X_{\left(i\right)}\right)}---\left(29\right)$

系统可靠度为：

$\hat{R}=1-{\hat{P}}_f---\left(30\right)$

案例仿真结果：

蒙特卡罗法仿真10000次，其失效概率为0.0749，失效概率对应的方差为5.2312e-5，失效概率达到稳定精确值用新方法的仿真次数见表2。

表2不同的MCMC次数下重要抽样法所需的模拟次数

取不同的M值时的失效概率变化和方差变化与蒙特卡罗法进行对比，结果如图2、3所示。

Claims (1)

1.一种基于马尔可夫链蒙特卡罗的舵机可靠性仿真抽样方法，其特征在于：

该方法依次按照下述四个阶段进行；

第一阶段：马尔可夫过程模拟

马尔可夫过程模拟主要包括以下四个步骤：

①选取马尔可夫链的初始状态X₀：

依据工程经验及简单的数值方法确定失效域中的一点作为X₀；

②确定随机转移抽样概率密度函数：

定义混合型随机转移抽样概率密度函数为P^*(X_D ^(j))p^*(ξ|X_C ^(j))，在此可以选择简单的均匀分布，则：

P^*(X_D ^(j))p^*(ξ|X_C ^(j))

式(1)中S为离散变量的状态维数；l_i是n-1维以X_C ^(j)为中心的超立方体在第i维的边长，它决定了X_C ^(j+1)偏离X_C ^(j)的最大允许距离，按照X_C ^(j+1)偏离X_C ^(j)的最大允许距离为三倍p^*(ξ|X_C ^(j))的标准差来确定l_i，则：

l_i＝6σ_iM^-1/(n+3)    (2)

式(2)中n为连续变量数目，M为马氏链步长即模拟样本数目，σ_i为h_opt(X_C)相对于X_C ^(j)的标准差的近似，依经验确定；

③确定马尔可夫链的第j+1个状态

马尔可夫链的第j+1个状态X^(j+1)是在前一个状态X^(j)的基础上，由分布P^*(X_D ^(j))p_*(ξ|X_C ^(j))和Metropolis准则来确定的；基于X^(j)，依分布P^*(X_D ^(j))p^*(ξ|X_C ^(j))产生备选状态ξ，计算备选状态ξ和X^(j)状态的条件概率密度函数的比值r，即

r＝q(ξ)/q(X^(j))             (3)

式(3)中q(X)＝I[X]P(X_D)f(X_C)。(4)

然后根据Metropolis接受准则，确定Markov转移的下一状态点：

$X^{(j+1)}=\left\{\begin{array}{ccc}\mathrm{\&xi;}& \mathrm{if}& r\&GreaterEqual;1\\ \mathrm{\&xi;}& \mathrm{if}& r<1,\mathrm{and\; r}\&GreaterEqual;\mathrm{random}\left[\mathrm{0,1}\right]\\ X^{\left(j\right)}& \mathrm{if}& r<1,\mathrm{and\; r}<\mathrm{random}\left[\mathrm{0,1}\right]\end{array}\right.---\left(5\right)$

④按上述步骤不断重复，产生极限分布为渐进最优的重要抽样密度函数的M个随机样本点{X⁽¹⁾，X⁽²⁾…，X^(M)}；

第二阶段：核密度估计

核密度估计，包括以下三个步骤：

①选取核密度函数：

通常核函数属于对称的密度函数族P，为计算方便，选择Gaussian密度函数，具体形式为：

$K\left(X_C\right)=\frac{1}{\sqrt{{\left(2\mathrm{\&pi;}\right)}^{n-1}\left|S\right|}}\exp (-\frac{1}{2}{X_C}^T{S}^{-1}{X}_C)---\left(6\right)$

式(6)中S为样本点集{X⁽¹⁾，X⁽²⁾…，X^(M)}的协方差阵，描述每个样本点在不同方向和范围上的数据分散性

$S=\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}({X_C}^{\left(j\right)}-{\stackrel{\&OverBar;}{X}}_C){({X_C}^{\left(j\right)}-{\stackrel{\&OverBar;}{X}}_C)}^T---\left(7\right)$

②确定窗口宽度参数和局部带宽因子：

在低/高概率密度区域选用较大/小的窗口宽度，具体分为：

局部带宽因子λ_j的求解：

${\mathrm{\&lambda;}}_j={\{{\left[\underset{k=1}{\overset{M}{\mathrm{\&Pi;}}}f\left(x^{\left(k\right)}\right)\right]}^{1/M}/f\left(x^{\left(j\right)}\right)\}}^{\mathrm{\&alpha;}}---\left(8\right)$

式(8)中0≤α≤1为灵敏因子，取α＝0.5；

窗口宽度参数w的求解：

$w={M_d}^{-\frac{1}{n+3}}---\left(9\right)$

式(9)中n为连续参数的个数，M_d为不同样本的个数(M_d≤M)；

③依据{X⁽¹⁾，X⁽²⁾…，X^(M)}，采用自适应宽核密度估计法，产生混合重要抽样概率密度函数k(X)：

$k\left(X\right)=P^h\left(X_D\right)\left[\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{1}{{\left(w{\mathrm{\&lambda;}}_j\right)}^{n-1}}K\left(\frac{X_C-{X_C}^{\left(j\right)}}{w{\mathrm{\&lambda;}}_j}\right)\right]---\left(10\right)$

在上式(10)中

$P^h\left(X_D\right)=\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{I}_{X_D}\left({X_D}^{\left(j\right)}\right),X_D=\mathrm{1,2},...,S---\left(11\right)$

在上式(11)中

第三阶段：重要抽样

从{1，2…，M}中均匀产生一个离散随机整数u，如果u＝j，则X_D＝X_D ^(j)，选取第j个分量的核抽样概率密度函数k_j(X_C)，来产生样本X_C，重复上述过程，直到得到N个样本点{X₍₁₎，X₍₂₎，…，X_(i)，…，X_(N)}；

$k_j\left(X_C\right)=\frac{1}{{\left(w{\mathrm{\&lambda;}}_j\right)}^{n-1}}K\left(\frac{X_C-{X_C}^{\left(j\right)}}{w{\mathrm{\&lambda;}}_j}\right)$

(13)

$=\frac{1}{{\left(\sqrt{2\mathrm{\&pi;}}w{\mathrm{\&lambda;}}_j\right)}^{n-2}}\frac{1}{\sqrt{2\mathrm{\&pi;}}w{\mathrm{\&lambda;}}_j\sqrt{S}}\exp (-\frac{(X_C-{X_C}^{\left(j\right)})^2}{2{\left(w{\mathrm{\&lambda;}}_j\right)}^{2}S})$

第四阶段：统计计算

根据重要抽样仿真样本点，进行失效概率估计：

${\hat{P}}_f=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}I\left[X_{\left(i\right)}\right]\frac{P\left(X_{D\left(i\right)}\right)f\left(X_{C\left(i\right)}\right)}{k\left(X_{\left(i\right)}\right)}---\left(14\right)$

系统可靠度为：

$\hat{R}=1-{\hat{P}}_f---\left(15\right)$

上述四大阶段，每一阶段又细分为几个步骤，共有9个步骤；马尔可夫过程模拟得到的样本点为下一阶段核密度估计进行失效域预估计之用；核密度估计得到的混合核抽样密度函数又是下一阶段的重要抽样函数；通过重要抽样得到的重要区域样本点又是下一阶段统计计算所必须的，依次进行的四大阶段环环相扣，缺一不可。

Images