CN105512404A - 基于混沌多项式展开的时变可靠性全局灵敏度分析方法 - Google Patents

Abstract

一种基于混沌多项式展开的时变可靠性全局灵敏度分析方法，包括以下步骤：一，在基于故障机理的可靠性与性能一体化模型的基础上，计算产品退化过程离散时刻的性能混沌多项式展开；二，基于正交实验设计方案，计算产品退化过程离散时刻的可靠性混沌多项式展开；三，根据移动最小二乘原理，计算产品退化过程时变可靠性混沌多项式展开；四，将上述时变可靠性混沌多项式展开的系数按照索博尔(Sobol’)分解重组，计算全局灵敏度Sobol’指标。该方法能够高效地计算复杂工程模型在退化过程中可靠性全局灵敏度分析问题，具有精度较高、适用范围广等特点。

Description

基于混沌多项式展开的时变可靠性全局灵敏度分析方法

(一)技术领域：

本发明提供一种基于混沌多项式展开的可靠性全局灵敏度分析方法，注重于解决复杂工程模型的可靠全局灵敏度分析问题，属于系统优化设计与可靠性设计的交叉技术领域。

(二)背景技术：

系统可靠性与性能一体化技术是一项在产品设计过程中考虑故障和环境扰动，采用可靠性优化和不确定性分析等方法来实现可靠性与性能综合设计分析的新技术。由于产品单元退化机理在时间轴上的发展及其在不同单元间的传播，产品性能输出随着时间的增长表现为一个逐渐退化的随机过程，时变可靠性作为可实时反映系统可靠性及质量特性的指标，逐渐得到一体化分析设计人员的重视。曾声奎、陈云霞等人对时变可靠性建模仿真技术开展了大量的工作。

灵敏度分析是在一项产品设计阶段研究系统设计变量对产品特性影响的方法。它包括局部灵敏度分析和全局灵敏度分析。局部灵敏度分析主要考虑由单个设计变量在标称值处的线性梯度所引起产品性能变化的比率；而全局灵敏度分析则能够衡量所有设计变量的不确定性对产品性能的综合作用，为设计人员优化设计方案提供一种有效手段。1990年，Sobol’提出一种基于方差的全局灵敏度分析指标(Sobol’指标)，它能够快速简便地计算出设计变量高阶交叉影响项，在工程界得到广泛应用。但是该方法需要利用蒙特卡洛积分，面对具有复杂耦合关系和随机过程的复杂工程模型，计算负担很大。后来Saltelli、吕震宙等人将Sobol’指标向可靠性领域进行扩展，但是为获得某些小失效概率的结果，往往也需要大量蒙特卡洛仿真作为支撑，进一步增加了计算成本。混沌多项式展开(Polynomialchaosexpanse,PCE)采用相互正交的多项式作为基底，它将系统性能输出投影到概率空间，用标准随机变量来表示设计变量不确定性，从而建立起产品性能与设计变量之间关系的代理模型，具有较高的精度且较小的计算量。2007年，BrunoSudret首次提出将PCE按照Sobol’分解、重组，然后利用重组后的PCE系数直接求解Sobol’指标，引起学术界和工程界的广泛关注。但是这些研究的对象没有涉及产品的退化过程，不能反映产品在整个寿命周期内的真实情况。

如果能提供一种针对时变可靠性灵敏度分析方法，即使对于复杂工程模型，也能就其退化过程进行高效、精确的分析，将显著地提高灵敏度分析的工程实用性，扩大其应用范围。

(三)发明内容：

(1)目的：

针对上述问题，本发明提出了一种基于混沌多项式展开的时变可靠性全局灵敏度分析方法，为可靠性与性能一体化设计提供一种客观、合理的分析技术手段。

(2)技术方案：

本发明是一种基于混沌多项式展开的时变可靠性全局灵敏度分析方法，该方法包括如下四个步骤：

步骤一：根据PCE与基于故障机理的可靠性与性能一体化仿真模型之间的关系，计算产品退化过程中各离散时刻的性能PCE，用来表示性能输出与关键设计变量的关系；

步骤二：根据蒙特卡洛抽样原理，计算各离散时刻的可靠性PCE，用来描述可靠性与关键设计变量的关系；

步骤三：根据移动最小二乘原理，计算退化过程的时变可靠性PCE，其中，所述的时变是指PCE系数随时间变化；

步骤四：在上述时变可靠性PCE的基础上，计算时变全局灵敏度的Sobol’s指标。

其中，在步骤一中所述的“计算产品退化过程中各离散时刻的性能PCE”，其计算步骤可细分为两步。

步骤1)确定性能PCE的基底和阶数

首先根据关键设计变量的分布类型选择相应的标准正交基底，如下列表1所示，其中，标准正交基底是相应标准随机变量的多项式函数。注意，当模型设计变量存在多种分布类型时，只能根据实际情况选择一种主要的正交基底形式。然后根据工程经验确定阶数p，并分别计算p＝k和p＝k+1两种情况下PCE及其相应的误差估计值，如两者误差差别不大，则可将阶数最终确定为k+1，否则再计算k+2，直至相邻两阶PCE的误差估计值基本一致，取最高阶PCE作为最终结果。这里k为大于1的正整数。

表1混沌多项式类型及对应的随机变量

表中，Hermite表示厄米多项式正交基底，N(0,1)表示均值为0、方差为1的正态分布；Legendre表示勒让德多项式正交基底，U[-1,1]表示上下界分别为1和-1的均匀分布；Laguerre表示拉盖尔多项式正交基底；GeneralizedLaguerre表示广义拉盖尔多项式正交基底，Γ(α+1,1)表示分布参数为α+1和1的伽马分布。上述正交基底的具体形式可以通过公开文献获得，这里不再赘述。

由此可以得到时刻t性能输出y(ξ；t)的PCE的一般形式为：

$y(\mathrm{\ξ};t)={\mathrm{PCE}}_t=\underset{j=0}{\overset{N-1}{\Σ}}{c}_j\left(t\right){\mathrm{\ψ}}_j\left(\mathrm{\ζ}\right),for\mathrm{\ζ}=({\mathrm{\ξ}}_1,...,{\mathrm{\ξ}}_n)---\left(1\right)$

式(1)中n是关键设计变量个数；p是PCE展开阶数；N为PCE所包含系数的总个数，它由n和p来确定，即：

$N=\left(\begin{array}{c}n+p\\ p\end{array}\right)=\frac{(n+p)!}{n!p!}---\left(2\right)$

${\left\{c_j\left(t\right)\right\}}_{j=0}^{N-1}$ 是PCE在各离散时刻t的系数， $\begin{array}{cc}{\mathrm{\ψ}}_j\left(\mathrm{\ξ}\right)=\underset{k=1}{\overset{n}{\Π}}{\mathrm{\φ}}_{{\mathrm{\α}}_k^j}\left({\mathrm{\ξ}}_k\right)& \left(\right|{\mathrm{\α}}^j|=\underset{k=1}{\overset{j}{\Σ}}{\mathrm{\α}}_k^j\≤p)\end{array},$ 其中是标准正交多项式基底，ξ_i～N(0,1)，i＝1,…,n。

步骤2)计算性能PCE的系数

根据待定系数个数N－1，在比PCE阶数度高一阶的标准正交多项式基底的根中选择合适数量的配点其中，ξ^k是PCE基底对应的标准随机变量的一组样本点；然后将配点转化变为仿真模型的变量输入，并求解系统响应y(ξ^k；t)。以Hermit正交基底为例，由于Hermit正交基底对已的标准随机变量ξ_h服从标准正态分布N(0,1)，则配点与仿真模型设计变量输入之间的转换关系如下列表2所示。

表2常见分布与标准正态分布关系

其中，是高斯误差函数，y＝exp(x)＝e^x是指数函数。

再将配点代入到式(1)中多项式部分，就能得到多组计算样本

${\{{\mathrm{\ψ}}_0\left({\mathrm{\ξ}}^k\right),{\mathrm{\ψ}}_1\left({\mathrm{\ξ}}^k\right),...,{\mathrm{\ψ}}_{N-1}\left({\mathrm{\ξ}}^k\right);y({\mathrm{\ξ}}^k;t)\}}_{k=1}^M,$ 利用多元线性回归求取系数 ${\left\{c_j\left(t\right)\right\}}_{j=0}^{N-1},$

即：

式中M为配点数目，为保证系数矩阵的条件数，样本数量要求不小于未知系数个数的两倍(M≥2N)，同时还要增加一种配点为零的选择方案，并且在布置配点吋，应尽量关于原点对称。

按上述步骤不断重复，获得各离散时刻的性能PCE。

$\begin{array}{ccccccc}t:& 0& t_1& t_2& t_3& \mathrm{...}& t_{M_t}\\ \↓& \↓& \↓& \↓& \↓& & \↓\\ PCE:& {\mathrm{PCE}}_0& {\mathrm{PCE}}_1& {\mathrm{PCE}}_2& {\mathrm{PCE}}_2& \mathrm{...}& {\mathrm{PCE}}_{M_t}\end{array}---\left(4\right)$

式中，M_t为离散时刻数目。

其中，在步骤二中所述的“计算各离散时刻的可靠性PCE”，其目的是获得关键设计变量的分布参数与产品可靠度之间的关系。其中，可以认为关键设计变量的分布参数在取值范围内的任何取值都是等可能的，即关键设计变量分布参数在取值范围内服从均匀分布。根据表1选择可靠性PCE的基底为Legendre基底；而可靠性PCE的阶数确定方法与性能PCE的阶数确定方法相同，但根据实际经验，一般和性能PCE同阶数。各离散时刻可靠性PCE的形式为：

$R(\mathrm{\ζ};t)={\mathrm{RPCE}}_t=\underset{j=0}{\overset{N_r-1}{\Σ}}{c}_{rj}\left(t\right){\mathrm{\ψ}}_{Lj}\left(\mathrm{\ζ}\right),for\mathrm{\ζ}={({\mathrm{\ζ}}_1,...,{\mathrm{\ζ}}_n)}^T---\left(5\right)$

式中， $\begin{array}{cc}{\mathrm{\ψ}}_{Lj}\left(\mathrm{\ξ}\right)=\underset{k=1}{\overset{n}{\Π}}{\mathrm{\φ}}_{{\mathrm{L\α}}_k^j}\left({\mathrm{\ξ}}_k\right)& \left(\right|{\mathrm{\α}}^j|=\underset{k=1}{\overset{j}{\Σ}}{\mathrm{\α}}_k^j\≤p)\end{array},$ 是Legendre正交多项式， ${\left\{c_{ri}\left(t\right)\right\}}_{i=0}^{N_r-1}$ 是可靠性PCE系数，N_r是可靠性PCE系数的总个数。其具体计算步骤如下：

步骤1)离散时刻可靠性PCE配点设计

根据工程经验，确定关键设计变量的分布参数的取值范围。然后将各关键设计变量的分布参数水平均划分为m等，进行随机搭配，得到可靠性PCE的试验方案。其中，总试验方案数目M′＝mⁿ，且M′不小于2N。然后将配点映射到区间[-1,1]上，得到均匀分布标准随机变量与试验方案的关系。一般地，令关键设计变量的分布参数μ_i服从区间[a,b]的均匀分布，即，μ_i～U[a,b]；均匀分布标准随机变量ζ_i～U[-1,1]。则有

${\mathrm{\ζ}}_i=\frac{2{\mathrm{\μ}}_i-a-b}{b-a}---\left(6\right)$

步骤2)各配点对应可靠度的计算

依据上述试验方案进行蒙特卡洛抽样，以获得相应的可靠度。该步骤需要利用到各离散时刻性能PCE，即y(ξ；t)＝PCE_t。该计算过程包括以下两个步骤：

I)建立试验方案中各个配点的关键设计变量的分布参数与步骤一中配点的关键设计变量的初始分布参数以及均匀分布标准随机变量之间的关系，得到反映试验方案的性能PCE修正随机变量和均匀分布标准随机变量为后续计算提供抽样分布。

例如，在以Hermit正交多项式为基底的PCE中，设某设计变量的初始分布参数为即均值为μ_i、方差为的正态分布，对应离散时刻性能PCE的标准正态变量ξ_i＝(X-μ)/σ；在第k个试验方案中则该试验方案对应的性能PCE修正随机变量为：

${\mathrm{\ξ}}_i^{\′\left(k\right)}=\frac{X_i^{\left(k\right)}-{\mathrm{\μ}}_i-{\mathrm{\Δ\μ}}_i^{\left(k\right)}}{{\mathrm{\σ}}_i}={\mathrm{\ξ}}_i-\frac{{\mathrm{\Δ\μ}}_i^{\left(k\right)}}{{\mathrm{\σ}}_i}---\left(7\right)$

又如，在以Hermit正交多项式为基底的PCE中，设某设计变量的初始分布参数为X_i～U[a_i,b_i]，即上下限为a_i和b_i的均匀分布；在第k个试验方案中则该试验方案对应的性能PCE修正随机变量为：

${\mathrm{\ξ}}^{\′}=\sqrt{2}{\mathrm{erf}}^{-1}\left(erf\right(\mathrm{\ξ}/\sqrt{2})-\frac{2{\mathrm{\Δ\μ}}_i}{b-a})---\left(8\right)$

式中， ${\mathrm{erf}}^{-1}\left(x\right)\≈\frac{2x}{\sqrt{\mathrm{\π}}}+\frac{2(4-\mathrm{\π})x^3}{2{\mathrm{\π}}^{1.5}}+\frac{(3{\mathrm{\π}}^2-40\mathrm{\π}+96)x^5}{15{\mathrm{\π}}^{2.5}}.$

此时，根据式(6)，第k个试验方案中对应的均匀分布标准随机变量值为：

${\mathrm{\ζ}}_i^k=\frac{2({\mathrm{\μ}}_i+{\mathrm{\Δ\μ}}_i^{\left(k\right)})-a-b}{b-a}---\left(9\right)$

II)利用各离散时刻的性能PCE和各试验方案的修正随机变量，进行蒙特卡洛抽样。修正随机变量反映了设计变量分布参数变化对失效概率的影响。例如，根据式(7)中对配点进行抽样时，PCE中随机变量的分布不再服从标准正态分布，而是服从参数为的正态分布。

然后根据抽样的结果计算可靠度：

$R_k\left(t\right)\≈\frac{num\left(y\right({\mathrm{\ξ}}^{\′\left(k\right)};t|MC)\∈{\mathrm{\Ω}}_f)}{N_{sim}}---\left(10\right)$

式(10)中，y(ξ^′(k)；t|MC)表示在t时刻对第k个试验方案的性能PCE修正随机变量进行蒙特卡洛抽样，N_sim是抽样次数，num(y(ξ′；t)∈Ω_f)表示仿真结果落入失效域的数目，R_k(t)表示对应于第k个试验方案的可靠度。注意，与均匀分布标准随机变量对应。

步骤3)计算离散时刻可靠性PCE的系数

在各离散时刻，根据上述步骤1)和步骤2)，对各试验方案对应的修正随机变量进行蒙特卡罗抽样计算可靠度的同时，将该试验方案对应的均匀分布标准随机变量代入(5)中的多项式部分，可以得到多组样本点然后类似步骤二中求解各离散时刻性能PCE系数的方法，利用多元线性回归求取可靠性PCE的系数。

式(11)中，为各试验方案对应的均匀分布标准随机变量值。

按上述步骤不断重复，获得各离散时刻的可靠性

$\begin{array}{ccccccc}t:& 0& t_1& t_2& t_3& \mathrm{...}& t_{M_t}\\ \↓& \↓& \↓& \↓& \↓& & \↓\\ \mathrm{RPC}E:& {\mathrm{RPCE}}_0& {\mathrm{RPCE}}_1& {\mathrm{RPCE}}_2& R{\mathrm{PCE}}_2& \mathrm{...}& {\mathrm{RPCE}}_{M_t}\end{array}---\left(12\right)$

式中，M_t为离散时间点数目。

其中，在步骤三中所述的“退化过程的时变可靠性PCE”，其计算过程如下。

这里的“时变”是指可靠性PCE的所有系数分别随时间变化，即展开系数分别是时间的函数。该函数利用已知的各离散时刻的值，通过移动最小二乘拟合，即

${\hat{c}}_{rj}\left(t\right)=p^T\left(t\right)b_j,j=0,1,...,N-1---\left(13\right)$

式(13)中，p(t)是移动最小二乘拟合基函数，一般可取p(t)＝(1,t,t²)^T，b_j是移动最小二乘拟合系数。其中移动最小二乘的计算原理可从公开资料中获得，其计算结果为：

$b_j={\[P^{T}w\left(t_{1:M_t}\right)P\]}^{-1}{p}^{T}w\left(t_{1:M_t}\right)c_{rj}\left(t_{1:M_t}\right)$

$P={\[p\left(t_1\right),p\left(t_2\right),...,p\left(t_{M_t}\right)\]}^T$ (14)

$w\left(t_{1:M_t}\right)=diag\left(w\right(t_1),w(t_2),...,w(t_{M_t}\left)\right)$

$c_{rj}\left(t_{1:M_t}\right)={\[c_{rj}\left(t_1\right),c_{rj}\left(t_2\right),...,c_{rj}\left(t_{M_t}\right)\]}^T$

式(14)中，M_t为离散时间点数目；w(t)＝w(||t-t_i||)是权函数，比较常用的有样条函数，径向基函数，高斯函数等。代入式(13)得到时变可靠性PCE：

$R(\mathrm{\ζ},t)=\underset{j=0}{\overset{N_r-1}{\Σ}}{\hat{c}}_{rj}\left(t\right){\mathrm{\ψ}}_{Lj}\left(\mathrm{\ζ}\right)---\left(15\right)$

式中，为通过MLS估计得到的时变可靠性PCE系数。

其中，在步骤四中所述的“计算时变全局灵敏度的Sobol’指标”，其具体计算步骤如下：

首先根据时变可靠性PCE计算时变可靠性方差D(R(ζ,t))：

$D\left(R\right(\mathrm{\ζ},t\left)\right)=\underset{j=1}{\overset{N-1}{\Σ}}{\hat{c}}_{rj}^2\left(t\right)E\left({\mathrm{\ψ}}_{Lj}^2\right(\mathrm{\ζ}\left)\right)---\left(16\right)$

式中， $\begin{array}{cc}E\left({\mathrm{\ψ}}_{Lj}^2\right)=\underset{k=1}{\overset{n}{\Π}}E\left({\mathrm{\φ}}_{{\mathrm{L\α}}_k^j}^2\right)=\underset{k=1}{\overset{n}{\Π}}\frac{2}{2{\mathrm{\α}}_k^j+1}& \left(\right|{\mathrm{\α}}^j|=\underset{k=1}{\overset{j}{\Σ}}{\mathrm{\α}}_k^j\≤p)\end{array}.$

然后将式(15)的系数按照Sobol’分解的形式进行重组，通过重组后的PCE的系数直接计算，具体步骤如下：

步骤1)PCE系数重组

系数重组按照如下规则进行：

(1)被加项不展开，将单独的各个变量的一次项、二次项……放在一堆。

(2)被加项不展开，将具有两个变量(交互作用)的被加项的一次项、二次项……放在一堆。

(3)被加项不展开，将具有三个变量(交互作用)的被加项的一次项、二次项……放在一堆。

重组后数学表达式为：

式(17)中，α＝(α₁,…,α_n)是一个整数序列序列，是的具体值，满足

步骤2)Sobol’s指标

然后利用PCE基底的正交性等优良性质，可以直接得到Sobol’s指标：

(3)本发明的优点及功效：

该方法能够高效地计算复杂工程模型在退化过程中可靠性全局灵敏度分析问题，具有精度较高、适用范围广等特点。

(四)附图说明：

图1是连杆机构结构示意图。

图2本发明所述方法流程图。

图3时变可靠性PCE系数曲线。

图4各关键设计变量一阶时变可靠性Sobol’指标曲线。

图5各关键设计变量总的时变可靠性Sobol’指标曲线。

图中符号说明如下：

(X_i，Y_i)第i点的横纵坐标

PCE_tt时刻性能混沌多项式展开

RPCE_tt时刻可靠性混沌多项式展开

R_tkt时刻第k个试验方案对应的可靠度值

时变可靠性混沌多项式展开系数

S_i第i个关键设计变量一阶的可靠性Sobol’指标

ST_i第i个关键设计变量总的可靠性Sobol’指标

(五)具体实施方式：

见图1-5，下面将结合某连杆装置对本发明作进一步的详细说明。该收放装置包括滑动杆、后主推臂、后支撑臂、铰链、前主推臂、前支撑臂、液压装置及负载构成，如下列表3所示，主要设计变量数据如表3所示。故障判据为最大摩擦阻力Max-Force不超过了液压装置所能输出力的最大值。研究目标是分析2000次收放过程中，连杆装置的时变可靠性对设计变量均值的灵敏度。

表3主要设计变量参数值

符号	含义	名义值	符号	含义	名义值
						XA	A点的横坐标	-910.49mm	Gd	负载重量	2600N
XB	B点的横坐标	-150mm	W1	铰链上连杆宽	50
						YC/YD	C点的纵坐标	-264.67mm	D1	铰链上连杆深	20
Rlug	铰链耳片半径	15mm	R1	作动筒半径	50
						L1	铰链上连杆长	70mm	R2	作动筒滑动杆半径	25
L2	作动筒长	500mm	W2/3	前/后支撑臂宽度	50
						L3	作动筒滑动杆	1580.16mm	D2	前/后支撑臂深度	20
L4/5	前/后支撑臂	50mm	W3/4	前/后主推臂宽度	50
						L6/7	前/后主推臂长	260mm	D3/4	前/后主推臂深度	20

本发明一种基于混沌多项式展开的时变可靠性全局灵敏度分析方法，如下列表4所示，其中已经建立基于故障机理的可靠性与性能一体化仿真模型。它包括如下步骤：

步骤一，计算产品退化过程中各离散时刻的性能PCE

根据工程经验，在本案例中关键设计变量为A点横坐标X_A、C点纵坐标Y_C、负载重量G_d以及铰链半径R_lug四项，它们的不确定性输入参数的概率分布如表4所示：

表4各输入变量的概率分布

符号	描述	设计均值/ai	标准差/bi	分布类型
					XA	A点横坐标	-910.49	3	正态分布
YC	C点纵坐标	-264.67	1	正态分布
					Gd	负载重量	2600	10	正态分布
Rlug	铰链半径	14.9	15.1	均匀分布

根据表4，选择Hermit多项式作为PCE的基底，其第零、一、二、三、四阶Hermit多项式的形式分别为1，ξ，ξ²-1，ξ³-3ξ和ξ⁴-6ξ²+3。相应的标准正交随机变量为标准正态随机变量。然后选择阶数p＝2和p＝3时分别构造PCE，其对应的系数个数N分别为15和35。

然后设计配点以2阶PCE为例，其配点的坐标从三阶Hermit多项式ξ³-3ξ的根中选，每个配点是这三个根的组合，配点方案的总数目M至少为30。由于根中已经包含0点，无需再额外增加配点为零的选择方案。考虑到配点尽量关于原点对称，其设计方案如下列表5所示。

表5二阶性能PCE配点方案

根据设计变量与标准正态随机变量之间的转化关系如下列表6所示，其中ξ_i～N(0,1)，i＝1,2,3,4，将然后将配点转化变为仿真模型的变量输入。以配点方案10为例，即仿真模型输入为(905.2938,-262.9379,2600,14.9)^T。

表6设计变量与标准变量之间转换关系

根据各配点所对应的仿真模型输入，可以计算得到各离散时刻系统的响应输出其中M_t为离散时刻的数目，对本案例而言，将一次收放结束定义为离散时间间隔，即M_t＝2000。再将配点代入到式(1)中多项式部分得到回归样本其中， $\begin{array}{cc}{\mathrm{\ψ}}_j\left(\mathrm{\ξ}\right)=\underset{k=1}{\overset{n}{\Π}}{\mathrm{\φ}}_{{\mathrm{\α}}_k^j}\left({\mathrm{\ξ}}_k\right)& \left(\right|{\mathrm{\α}}^j|=\underset{k=1}{\overset{j}{\Σ}}{\mathrm{\α}}_k^j\≤p)\end{array}.$ 以二阶PCE的配点方案10在初始时刻的回归样本为为例，具体形式为

即

$\{1,-\sqrt{3},\sqrt{3},0,0,2,2,-1,-1,-3,0,0,0,0,0;8067.17\}$

利用式(3)，可以求出各离散时刻性能PCE。以在初始时刻三阶PCE为例，其回归方程组为

计算结果为：

$\begin{array}{c}y(t=0)=8065.707-97.314{\mathrm{\ξ}}_1-132.979{\mathrm{\ξ}}_2+26.375{\mathrm{\ξ}}_3+13.718{\mathrm{\ξ}}_4\\ -10.074({\mathrm{\ξ}}_1^2-1)-3.623({\mathrm{\ξ}}_2^2-1)+2.37({\mathrm{\ξ}}_3^2-1)-2.316({\mathrm{\ξ}}_4^2-1)\\ +10.718{\mathrm{\ξ}}_1{\mathrm{\ξ}}_2-10.762{\mathrm{\ξ}}_1{\mathrm{\ξ}}_3+1.316{\mathrm{\ξ}}_1{\mathrm{\ξ}}_4-7.982{\mathrm{\ξ}}_2{\mathrm{\ξ}}_3-1.235{\mathrm{\ξ}}_2{\mathrm{\ξ}}_4\\ +15.073{\mathrm{\ξ}}_3{\mathrm{\ξ}}_4-4.713({\mathrm{\ξ}}_1^3-3{\mathrm{\ξ}}_1)+4.615({\mathrm{\ξ}}_2^3-3{\mathrm{\ξ}}_2)-2.656({\mathrm{\ξ}}_3^3-3{\mathrm{\ξ}}_3)\\ +2.443({\mathrm{\ξ}}_4^3-3{\mathrm{\ξ}}_4)+2.497({\mathrm{\ξ}}_1{\mathrm{\ξ}}_2^2-{\mathrm{\ξ}}_1)-6.03({\mathrm{\ξ}}_1{\mathrm{\ξ}}_3^2-{\mathrm{\ξ}}_1)-0.27({\mathrm{\ξ}}_1{\mathrm{\ξ}}_4^2-{\mathrm{\ξ}}_1)\\ +6.046({\mathrm{\ξ}}_2{\mathrm{\ξ}}_1^2-{\mathrm{\ξ}}_2)-1.863({\mathrm{\ξ}}_2{\mathrm{\ξ}}_3^2-{\mathrm{\ξ}}_2)+3.91({\mathrm{\ξ}}_2{\mathrm{\ξ}}_4^2-{\mathrm{\ξ}}_2)+5.946({\mathrm{\ξ}}_3{\mathrm{\ξ}}_1^2-{\mathrm{\ξ}}_3)\\ -1.516({\mathrm{\ξ}}_3{\mathrm{\ξ}}_2^2-{\mathrm{\ξ}}_3)+3.229({\mathrm{\ξ}}_3{\mathrm{\ξ}}_4^2-{\mathrm{\ξ}}_3)+3.243({\mathrm{\ξ}}_4{\mathrm{\ξ}}_1^2-{\mathrm{\ξ}}_4)-2.058({\mathrm{\ξ}}_4{\mathrm{\ξ}}_2^2-{\mathrm{\ξ}}_4)\\ -4.044({\mathrm{\ξ}}_4{\mathrm{\ξ}}_3^2-{\mathrm{\ξ}}_4)+5.847{\mathrm{\ξ}}_1{\mathrm{\ξ}}_2{\mathrm{\ξ}}_3-2.327{\mathrm{\ξ}}_1{\mathrm{\ξ}}_2{\mathrm{\ξ}}_4+2.378{\mathrm{\ξ}}_1{\mathrm{\ξ}}_3{\mathrm{\ξ}}_4+1.581{\mathrm{\ξ}}_2{\mathrm{\ξ}}_3{\mathrm{\ξ}}_4\end{array}---\left(22\right)$

同时，由2阶PCE、3阶PCE以及MCS得到结果的误差在0.1％左右，如下列表7所示，满足精度要求，故最终选择PCE的阶数为p＝3，基底个数即PCE系数总个数为35。

表7各阶PCE对最大力估计值的对比

方法	调用模型次数	均值	方差(×105)
				MCS	10000	8074.602	6.3366
2阶PCE	30	8038.331	5.9674
				相对误差	---	0.0045	0.0583
3阶PCE	70	8056.497	6.2845
				相对误差	---	0.0022	0.0082

步骤二，计算产品退化过程中各离散时刻的可靠性PCE。

可靠性PCE的基底为Legendre正交多项式，其第零、一、二、三、四阶Hermit多项式的形式分别为1，ζ，和相应的标准正交随机变量为均匀分布标准随机变量ζ～U[-1,1]，可靠性PCE的阶数取3阶。

首先构造试验方案。根据工程经验，设计变量的设计均值变化范围为设计均值的±1％。将各实验因素(关键设计变量均值)在其变化范围内等分为3个水平，进行随机搭配，得到共计3⁴＝81个试验方案，达到了PCE基底个数的2倍，满足要求；其次，将各试验方案映射到[-1,1]，获得其对应的均匀分布标准随机变量值；同时计算各试验方案对应性能PCE修正随机变量值，如下列表8所示。

表8试验方案及对应均匀分布标准随机变量值

然后，计算试验方案对应的性能PCE修正随机变量，以试验方案1为例，

关键设计变量均值的取值为(μ₁,μ₂,μ₃,μ₄)^T＝(905.94,-263.35,2587.00,14.83)^T，均值增量分别为 ${\mathrm{\Δ\μ}}_1^{\left(1\right)}=-4.55,{\mathrm{\Δ\μ}}_2^{\left(1\right)}=1.32,{\mathrm{\Δ\μ}}_3^{\left(1\right)}=-13.00$ 以及 ${\mathrm{\Δ\μ}}_4^{\left(1\right)}=-\mathrm{0.07.}$ 则性能PCE修正随机变量分别为：

${\mathrm{\ξ}}_1^{\′}={\mathrm{\ξ}}_1-\frac{{\mathrm{\Δ\μ}}_1^{\left(1\right)}}{{\mathrm{\σ}}_1}={\mathrm{\ξ}}_1+1.52$

${\mathrm{\ξ}}_2^{\′}={\mathrm{\ξ}}_2-\frac{{\mathrm{\Δ\μ}}_2^{\left(1\right)}}{{\mathrm{\σ}}_2}={\mathrm{\ξ}}_2-1.32$ (23)

${\mathrm{\ξ}}_3^{\′}={\mathrm{\ξ}}_3-\frac{{\mathrm{\Δ\μ}}_3^{\left(1\right)}}{{\mathrm{\σ}}_3}={\mathrm{\ξ}}_3+1.3$

${\mathrm{\ξ}}_4^{\′}=\sqrt{2}{\mathrm{erf}}^{-1}\left(erf\right({\mathrm{\ξ}}_4/\sqrt{2})+0.47)$

通过对第k个试验方案对应的性能PCE修正随机变量进行蒙特卡洛抽样，结合故障判据，按照式(10)计算一系列的可靠度数值其中N_sim取10000。再将对应的均匀分布标准随机变量值代入到式(5)中多项式部分得到回归样本

${\{{\mathrm{\ψ}}_{L0}\left({\mathrm{\ζ}}^k\right),{\mathrm{\ψ}}_{L1}\left({\mathrm{\ζ}}^k\right),...,{\mathrm{\ψ}}_{L,N-1}\left({\mathrm{\ζ}}^k\right);R_k\left(t\right)\}}_{k=1}^{M^{\′}},$

其中， $\begin{array}{cc}{\mathrm{\ψ}}_{Lj}\left(\mathrm{\ξ}\right)=\underset{k=1}{\overset{n}{\Π}}{\mathrm{\φ}}_{{\mathrm{L\α}}_k^j}\left({\mathrm{\ξ}}_k\right)& \left(\right|{\mathrm{\α}}^j|=\underset{k=1}{\overset{j}{\Σ}}{\mathrm{\α}}_k^j\≤p)\end{array}.$ 以试验方案1在初始时刻的回归样本为为例，回归样本的具体形式为：

即

$\left\{\begin{array}{c}1,-0.5,0.5,-0.5,-0.5,-0.125,-0.125,-0.125,-0.125,-0.25,\\ 0.25,0.25,-0.25,-0.25,0.25,0.44,-0.44,0.44,\\ 0.44,0.06,0.06,-0.125,-0.06,\\ -0.06,0.125,0.06,0.06,-0.125,0.06,\\ 0.06,0.06,0.125,0.125,-0.125,0.125;0.9997\end{array}\right\}$

最后，根据回归方程组计算可靠性PCE。以初始时刻为例，其余不再赘述。

得可靠性PCE：

$\begin{array}{c}R(t=0)=0.9929-0.143{\mathrm{\ζ}}_1-0.076{\mathrm{\ζ}}_2-0.086{\mathrm{\ζ}}_3+0.085{\mathrm{\ζ}}_4\\ -0.099(3{\mathrm{\ζ}}_1^2-1)-0.076(3{\mathrm{\ζ}}_2^2-1)+0.623(3{\mathrm{\ζ}}_3^2-1)-0.028(3{\mathrm{\ζ}}_4^2-1)\\ +0.028{\mathrm{\ζ}}_1{\mathrm{\ζ}}_2-0.134{\mathrm{\ζ}}_1{\mathrm{\ζ}}_3^{\′}+0.038{\mathrm{\ζ}}_1{\mathrm{\ζ}}_4+0.37{\mathrm{\ζ}}_2{\mathrm{\ζ}}_3+0.087{\mathrm{\ζ}}_2{\mathrm{\ζ}}_4-0.014{\mathrm{\ζ}}_3{\mathrm{\ζ}}_4\\ -0.037(5{\mathrm{\ζ}}_1^3-3{\mathrm{\ζ}}_1)-0.021(5{\mathrm{\ζ}}_2^3-3{\mathrm{\ζ}}_2)+0.062(5{\mathrm{\ζ}}_3^3-3{\mathrm{\ζ}}_3)+0.031(5{\mathrm{\ζ}}_4^3-3{\mathrm{\ζ}}_4)\\ -0.025(3{\mathrm{\ζ}}_1{\mathrm{\ζ}}_2^2-{\mathrm{\ζ}}_1)-0.04(3{\mathrm{\ζ}}_1{\mathrm{\ζ}}_3^2-{\mathrm{\ζ}}_1)+0.07(3{\mathrm{\ζ}}_1{\mathrm{\ζ}}_4^2-{\mathrm{\ζ}}_1)+0.036(3{\mathrm{\ζ}}_2{\mathrm{\ζ}}_1^2-{\mathrm{\ζ}}_2)\\ -0.047(3{\mathrm{\ζ}}_2{\mathrm{\ζ}}_3^2-{\mathrm{\ζ}}_2)+0.032(3{\mathrm{\ζ}}_2{\mathrm{\ζ}}_4^2-{\mathrm{\ζ}}_2)-0.04(3{\mathrm{\ζ}}_3{\mathrm{\ζ}}_1^2-{\mathrm{\ζ}}_3)-0.077(3{\mathrm{\ζ}}_3{\mathrm{\ζ}}_2^2-{\mathrm{\ζ}}_3)\\ +0.001(3{\mathrm{\ζ}}_3{\mathrm{\ζ}}_4^2-{\mathrm{\ζ}}_3)+0.07(3{\mathrm{\ζ}}_4{\mathrm{\ζ}}_1^2-{\mathrm{\ζ}}_4)+0.014(3{\mathrm{\ζ}}_4{\mathrm{\ζ}}_2^2-{\mathrm{\ζ}}_4)+0.102(3{\mathrm{\ζ}}_4{\mathrm{\ζ}}_3^2-{\mathrm{\ζ}}_4)\\ -0.098{\mathrm{\ζ}}_1{\mathrm{\ζ}}_2{\mathrm{\ζ}}_3-0.046{\mathrm{\ζ}}_1{\mathrm{\ζ}}_2{\mathrm{\ζ}}_4+0.018{\mathrm{\ζ}}_1{\mathrm{\ζ}}_3{\mathrm{\ζ}}_4-0.016{\mathrm{\ζ}}_2{\mathrm{\ζ}}_3{\mathrm{\ζ}}_4\end{array}---\left(25\right)$

步骤三，计算产品退化过程时变可靠性PCE。

选择三次样条函数作为移动最小二乘的权重函数，即：

$w\left(s\right)=\left\{\begin{array}{cc}\frac{2}{3}-4s^2+4s^3& s\&le;\frac{1}{2},\\ \frac{4}{3}-4s+4s^2-\frac{4}{3}{s}^3& \frac{1}{2}<s\&le;1,\\ 0& s>1\end{array}\right.---\left(26\right)$

式(26)中，s＝||t-t_i||，t_i是离散时刻。

然后结合上述离散时刻可靠性PCE的系数，根据式(14)，计算时变系数表达式。以式(15)中系数为例，其计算结果如图3所示。图中，离散的点是经过各离散时刻PCE的系数，曲线是利用移动最小二乘原理估计得到的曲线。

步骤四，计算时变全局灵敏度的Sobol’指标。

首先，根据式(16)求出时变可靠性PCE计算时变可靠性方差D(R(ζ,t))。

然后，将上述时变可靠性PCE按照式(17)重新展开。以初始时刻为例时变可靠性PCE为例，其系数重组后形式如式(27)所示。

$\begin{array}{c}R(t=0)=0.9929\\ -0.143{\mathrm{\&zeta;}}_1-0.076{\mathrm{\&zeta;}}_2-0.086{\mathrm{\&zeta;}}_3+0.085{\mathrm{\&zeta;}}_4\\ -0.099(3{\mathrm{\&zeta;}}_1^2-1)-0.076(3{\mathrm{\&zeta;}}_2^2-1)+0.623(3{\mathrm{\&zeta;}}_3^2-1)-0.028(3{\mathrm{\&zeta;}}_4^2-1)\\ -0.037(5{\mathrm{\&zeta;}}_1^3-3{\mathrm{\&zeta;}}_1)-0.021(5{\mathrm{\&zeta;}}_2^3-3{\mathrm{\&zeta;}}_2)+0.062(5{\mathrm{\&zeta;}}_3^3-3{\mathrm{\&zeta;}}_3)+0.031(5{\mathrm{\&zeta;}}_4^3-3{\mathrm{\&zeta;}}_4)\\ +0.028{\mathrm{\&zeta;}}_1{\mathrm{\&zeta;}}_2-0.134{\mathrm{\&zeta;}}_1{\mathrm{\&zeta;}}_3^{\&prime;}+0.038{\mathrm{\&zeta;}}_1{\mathrm{\&zeta;}}_4+0.37{\mathrm{\&zeta;}}_2{\mathrm{\&zeta;}}_3+0.087{\mathrm{\&zeta;}}_2{\mathrm{\&zeta;}}_4-0.014{\mathrm{\&zeta;}}_3{\mathrm{\&zeta;}}_4\\ -0.025(3{\mathrm{\&zeta;}}_1{\mathrm{\&zeta;}}_2^2-{\mathrm{\&zeta;}}_1)-0.04(3{\mathrm{\&zeta;}}_1{\mathrm{\&zeta;}}_3^2-{\mathrm{\&zeta;}}_1)+0.07(3{\mathrm{\&zeta;}}_1{\mathrm{\&zeta;}}_4^2-{\mathrm{\&zeta;}}_1)\\ +0.036(3{\mathrm{\&zeta;}}_2{\mathrm{\&zeta;}}_1^2-{\mathrm{\&zeta;}}_2)-0.047(3{\mathrm{\&zeta;}}_2{\mathrm{\&zeta;}}_3^2-{\mathrm{\&zeta;}}_2)+0.032(3{\mathrm{\&zeta;}}_2{\mathrm{\&zeta;}}_4^2-{\mathrm{\&zeta;}}_2)\\ -0.04(3{\mathrm{\&zeta;}}_3{\mathrm{\&zeta;}}_1^2-{\mathrm{\&zeta;}}_3)-0.077(3{\mathrm{\&zeta;}}_3{\mathrm{\&zeta;}}_2^2-{\mathrm{\&zeta;}}_3)+0.001(3{\mathrm{\&zeta;}}_3{\mathrm{\&zeta;}}_4^2-{\mathrm{\&zeta;}}_3)\\ +0.07(3{\mathrm{\&zeta;}}_4{\mathrm{\&zeta;}}_1^2-{\mathrm{\&zeta;}}_4)+0.014(3{\mathrm{\&zeta;}}_4{\mathrm{\&zeta;}}_2^2-{\mathrm{\&zeta;}}_4)+0.102(3{\mathrm{\&zeta;}}_4{\mathrm{\&zeta;}}_3^2-{\mathrm{\&zeta;}}_4)\\ -0.098{\mathrm{\&zeta;}}_1{\mathrm{\&zeta;}}_2{\mathrm{\&zeta;}}_3-0.046{\mathrm{\&zeta;}}_1{\mathrm{\&zeta;}}_2{\mathrm{\&zeta;}}_4+0.018{\mathrm{\&zeta;}}_1{\mathrm{\&zeta;}}_3{\mathrm{\&zeta;}}_4-0.016{\mathrm{\&zeta;}}_2{\mathrm{\&zeta;}}_3{\mathrm{\&zeta;}}_4\end{array}---\left(27\right)$

最后，根据式(19)和式(20)可直接求出Sobol’指标，如图4和图5所示。

Claims (5)

1.一种基于混沌多项式展开的时变可靠性全局灵敏度分析方法，其特征在于：该方法包括如下四个步骤：

步骤一：根据PCE与基于故障机理的可靠性与性能一体化仿真模型之间的关系，计算产品退化过程中各离散时刻的性能PCE，用来表示性能输出与关键设计变量的关系；

步骤二：根据蒙特卡洛抽样原理，计算各离散时刻的可靠性PCE，用来描述可靠性与关键设计变量的关系；

步骤三：根据移动最小二乘原理，计算退化过程的时变可靠性PCE，其中，所述的时变是指PCE系数随时间变化；

步骤四：在上述时变可靠性PCE的基础上，计算时变全局灵敏度的Sobol’s指标。

2.根据权利要求1所述的一种基于混沌多项式展开的时变可靠性全局灵敏度分析方法，其特征在于：在步骤一中所述的“计算产品退化过程中各离散时刻的性能PCE”，计算步骤细分为两步：

步骤1.1确定性能PCE的基底和阶数

首先，根据关键设计变量的分布类型选择相应的标准正交基底，如下所列；

表1所示，标准正交基底是相应标准随机变量的多项式函数，然后确定阶数p，并分别计算p＝k和p＝k+1两种情况下PCE及其相应的误差估计值，如两者误差差别不大，则将阶数最终确定为k+1，否则再计算k+2，直至相邻两阶PCE的误差估计值一致，取最高阶PCE作为最终结果，这里k为大于1的正整数；

表1混沌多项式类型及对应的随机变量

设计变量分布类型 混沌多项式基底形式 标准随机变量 支持区间 正态分布 Hermite ξ～N(0,1) [-∞,∞] 均匀分布 Legendre ξ～U[-1,1] [-1,1] 对数正态分布 Hermite ξ～N(0,1) [-∞,∞] 伽马分布 Generalized Laguerre ξ～Γ(α+1,1) [0,∞]

表中，Hermite表示厄米多项式正交基底，N(0,1)表示均值为0、方差为1的正态分布；Legendre表示勒让德多项式正交基底，U[-1,1]表示上下界分别为1和-1的均匀分布；Laguerre表示拉盖尔多项式正交基底；GeneralizedLaguerre表示广义拉盖尔多项式正交基底，Γ(α+1,1)表示分布参数为α+1和1的伽马分布；上述正交基底的具体形式能通过公开文献获得；

由此能得到时刻t性能输出y(ξ；t)的PCE的形式为：

$\begin{array}{ccc}y(\mathrm{\&xi;};t)={\mathrm{PCE}}_t=\underset{j=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{c}_j\left(t\right){\mathrm{\&psi;}}_j\left(\mathrm{\&xi;}\right),& for& \mathrm{\&xi;}=({\mathrm{\&xi;}}_1,...,{\mathrm{\&xi;}}_n)\end{array}---\left(1\right)$

式(1)中n是关键设计变量个数；p是PCE展开阶数；N为PCE所包含系数的总个数，它由n和p来确定，即：

$N=\left(\begin{array}{c}n+p\\ p\end{array}\right)=\frac{(n+p)!}{n!p!}---\left(2\right)$

是PCE在各离散时刻t的系数， ${\mathrm{\&psi;}}_j\left(\mathrm{\&xi;}\right)=\underset{k=1}{\overset{n}{\&Pi;}}{\mathrm{\&phi;}}_{{\mathrm{\&alpha;}}_k^j}\left({\mathrm{\&xi;}}_k\right)\left(\right|{\mathrm{\&alpha;}}^j|=\underset{k=1}{\overset{j}{\&Sigma;}}{\mathrm{\&alpha;}}_k^j\&le;p),$ 其中，是标准正交多项式基底，ξ_i～N(0,1)，i＝1,…,n；

步骤1.2计算性能PCE的系数

根据待定系数个数N－1，在比PCE阶数度高一阶的标准正交多项式基底的根中选择合适数量的配点其中，ξ^k是PCE基底对应的标准随机变量的一组样本点；然后将配点转化变为仿真模型的变量输入，并求解系统响应y(ξ^k；t)；由于Hermit正交基底对已的标准随机变量ξ_h服从标准正态分布N(0,1)，则配点与仿真模型设计变量输入之间的转换关系如下列表2所示

表2常见分布与标准正态分布关系

其中，是高斯误差函数，y＝exp(x)＝e^x是指数函数；

再将配点代入到式(1)中多项式部分，

就能得到多组计算样本 ${\{{\mathrm{\&psi;}}_0\left({\mathrm{\&xi;}}^k\right),{\mathrm{\&psi;}}_1\left({\mathrm{\&xi;}}^k\right),...,{\mathrm{\&psi;}}_{N-1}\left({\mathrm{\&xi;}}^k\right);y({\mathrm{\&xi;}}^k;t)\}}_{k=1}^M,$ 利用多元线性回归求取系数即：

式中M为配点数目，为保证系数矩阵的条件数，样本数量要求不小于未知系数个数的两倍，即M≥2N，同时还要增加一种配点为零的选择方案，并且在布置配点吋，应关于原点对称；

按上述步骤不断重复，获得各离散时刻的性能PCE；

$\begin{array}{ccccccc}t:& 0& t_1& t_2& t_3& ...& t_{M_t}\\ \&DownArrow;& \&DownArrow;& \&DownArrow;& \&DownArrow;& \&DownArrow;& & \&DownArrow;\\ PCE:& {\mathrm{PCE}}_0& {\mathrm{PCE}}_1& {\mathrm{PCE}}_2& {\mathrm{PCE}}_2& ...& {\mathrm{PCE}}_{M_t}\end{array}---\left(4\right)$

式中，M_t为离散时刻数目。

3.根据权利要求1所述的一种基于混沌多项式展开的时变可靠性全局灵敏度分析方法，其特征在于：在步骤二中所述的“计算各离散时刻的可靠性PCE”，其目的是获得关键设计变量的分布参数与产品可靠度之间的关系；其中，认为关键设计变量的分布参数在取值范围内的任何取值都是等可能的，即关键设计变量分布参数在取值范围内服从均匀分布；

根据表1选择可靠性PCE的基底为Legendre基底；而可靠性PCE的阶数确定方法与性能PCE的阶数确定方法相同，各离散时刻可靠性PCE的形式为：

$\begin{array}{ccc}R(\mathrm{\&xi;};t)={\mathrm{RPCE}}_t=\underset{j=0}{\overset{N_r-1}{\mathrm{\&Sigma;}}}{c}_{rj}\left(t\right){\mathrm{\&psi;}}_{Lj}\left(\mathrm{\&zeta;}\right),& for& \mathrm{\&zeta;}={({\mathrm{\&zeta;}}_1,...,{\mathrm{\&zeta;}}_n)}^T\end{array}---\left(5\right)$

式中， ${\mathrm{\&psi;}}_{Lj}\left(\mathrm{\&zeta;}\right)=\underset{k=1}{\overset{n}{\&Pi;}}{\mathrm{\&phi;}}_{{\mathrm{L\&alpha;}}_k^j}\left({\mathrm{\&zeta;}}_k\right)\left(\right|{\mathrm{\&alpha;}}^j|=\underset{k=1}{\overset{j}{\&Sigma;}}{\mathrm{\&alpha;}}_k^j\&le;p),$ 是Legendre正交多项式，是可靠性PCE系数，N_r是可靠性PCE系数的总个数；具体计算步骤如下：

步骤2.1离散时刻可靠性PCE配点设计

确定关键设计变量的分布参数的取值范围，然后将各关键设计变量的分布参数水平均划分为m等，进行随机搭配，得到可靠性PCE的试验方案；其中，总试验方案数目M′＝mⁿ，且M′不小于2N；然后将配点映射到区间[-1,1]上，得到均匀分布标准随机变量与试验方案的关系；令关键设计变量的分布参数μ_i服从区间[a,b]的均匀分布，即，μ_i～U[a,b]；均匀分布标准随机变量ζ_i～U[-1,1]，则有

${\mathrm{\&zeta;}}_i=\frac{2{\mathrm{\&mu;}}_i-a-b}{b-a}---\left(6\right)$

步骤2.2各配点对应可靠度的计算

该步骤需要利用到各离散时刻性能PCE，即y(ξ；t)＝PCE_t；该计算过程包括以下两个步骤：

步骤2.21建立试验方案中各个配点的关键设计变量的分布参数与步骤一中配点的关键设计变量的初始分布参数以及均匀分布标准随机变量之间的关系，得到反映试验方案的性能PCE修正随机变量和均匀分布标准随机变量为后续计算提供抽样分布；

在以Hermit正交多项式为基底的PCE中，设某设计变量的初始分布参数为即均值为μ_i、方差为的正态分布，对应离散时刻性能PCE的标准正态变量ξ_i＝(X-μ)/σ；在第k个试验方案中则该试验方案对应的性能PCE修正随机变量为：

${\mathrm{\&xi;}}_i^{\&prime;\left(k\right)}=\frac{X_i^{\left(k\right)}-{\mathrm{\&mu;}}_i-{\mathrm{\&Delta;\&mu;}}_i^{\left(k\right)}}{{\mathrm{\&sigma;}}_i}={\mathrm{\&xi;}}_i-\frac{{\mathrm{\&Delta;\&mu;}}_i^{\left(k\right)}}{{\mathrm{\&sigma;}}_i}---\left(7\right)$

在以Hermit正交多项式为基底的PCE中，设某设计变量的初始分布参数为X_i～U[a_i,b_i]，即上下限为a_i和b_i的均匀分布；在第k个试验方案中则该试验方案对应的性能PCE修正随机变量为：

${\mathrm{\&xi;}}^{\&prime;}=\sqrt{2}{\mathrm{erf}}^{-1}\left(erf\right(\mathrm{\&xi;}/\sqrt{2})-\frac{2{\mathrm{\&Delta;\&mu;}}_i}{b-a})---\left(8\right)$

式中， ${\mathrm{erf}}^{-1}\left(x\right)\&ap;\frac{2x}{\sqrt{\mathrm{\&pi;}}}+\frac{2(4-\mathrm{\&pi;})x^3}{2{\mathrm{\&pi;}}^{1.5}}+\frac{(3{\mathrm{\&pi;}}^2-40\mathrm{\&pi;}+96)x^5}{15{\mathrm{\&pi;}}^{2.5}};$

根据式(6)，第k个试验方案中对应的均匀分布标准随机变量值为：

${\mathrm{\&zeta;}}_i^k=\frac{2({\mathrm{\&mu;}}_i+{\mathrm{\&Delta;\&mu;}}_i^{\left(k\right)})-a-b}{b-a}---\left(9\right)$

步骤2.22利用各离散时刻的性能PCE和各试验方案的修正随机变量，进行蒙特卡洛抽样；修正随机变量反映了设计变量分布参数变化对失效概率的影响，根据式(7)中对配点进行抽样时，PCE中随机变量的分布不再服从标准正态分布，而是服从参数为的正态分布；

然后计算可靠度：

$R_k\left(t\right)\&ap;\frac{num\left(y\right({\mathrm{\&xi;}}^{\&prime;\left(k\right)};t\left|MC\right)\&Element;{\mathrm{\&Omega;}}_f)}{N_{sim}}---\left(10\right)$

式(10)中，y(ξ′^(k)；t|MC)表示在t时刻对第k个试验方案的性能PCE修正随机变量进行蒙特卡洛抽样，N_sim是抽样次数，num(y(ξ′；t)∈Ω_f)表示仿真结果落入失效域的数目，R_k(t)表示对应于第k个试验方案的可靠度；注意，与均匀分布标准随机变量对应；

步骤2.3计算离散时刻可靠性PCE的系数

在各离散时刻，根据上述步骤2.1和步骤2.2，对各试验方案对应的修正随机变量进行蒙特卡罗抽样计算可靠度的同时，将该试验方案对应的均匀分布标准随机变量代入(5)中的多项式部分，得到多组样本点 ${\{{\mathrm{\&psi;}}_{L0}\left({\mathrm{\&zeta;}}^k\right),{\mathrm{\&psi;}}_{L1}\left({\mathrm{\&zeta;}}^k\right),...,{\mathrm{\&psi;}}_{L,N-1}\left({\mathrm{\&zeta;}}^k\right);R_k\left(t\right)\}}_{k=1}^{M^{\&prime;}},$ 然后根据步骤二中求解各离散时刻性能PCE系数的方法，利用多元线性回归求取可靠性PCE的系数；

式(11)中，为各试验方案对应的均匀分布标准随机变量值；

获得各离散时刻的可靠性

$\begin{array}{ccccccc}t:& 0& t_1& t_2& t_3& ...& t_{M_t}\\ \&DownArrow;& \&DownArrow;& \&DownArrow;& \&DownArrow;& \&DownArrow;& & \&DownArrow;\\ RPCE:& {\mathrm{PPCE}}_0& {\mathrm{RPCE}}_1& {\mathrm{RPCE}}_2& {\mathrm{RPCE}}_2& ...& {\mathrm{RPCE}}_{M_t}\end{array}---\left(12\right)$

式中，M_t为离散时间点数目。

4.根据权利要求1所述的一种基于混沌多项式展开的时变可靠性全局灵敏度分析方法，其特征在于：在步骤三中所述的“退化过程的时变可靠性PCE”，计算过程如下：

这里的“时变”是指可靠性PCE的所有系数分别随时间变化，即展开系数分别是时间的函数:，该函数利用已知的各离散时刻的值，通过移动最小二乘拟合，即

${\hat{c}}_{rj}\left(t\right)=p^T\left(t\right)b_j,j=0,1,...,N-1---\left(13\right)$

式(13)中，p(t)是移动最小二乘拟合基函数，取p(t)＝(1,t,t²)^T，b_j是移动最小二乘拟合系数；其中，移动最小二乘的计算原理能从公开资料中获得，计算结果为：

$b_j={\&lsqb;P^{T}w\left(t_{1:M_t}\right)P\&rsqb;}^{-1}{P}^{T}w\left(t_{1:M_t}\right)c_{rj}\left(t_{1:M_t}\right)$

$P={\&lsqb;p\left(t_1\right),p\left(t_2\right),...,p\left(t_{M_t}\right)\&rsqb;}^T$

(14)

$w\left(t_{1:M_t}\right)=diag\left(w\right(t_1),w(t_2)...,w(t_{M_t}\left)\right)$

$c_{rj}\left(t_{1:M_t}\right)={\&lsqb;c_{rj}\left(t_1\right),c_{rj}\left(t_2\right),...,c_{rj}\left(t_{M_t}\right)\&rsqb;}^T$

式(14)中，M_t为离散时间点数目；w(t)＝w(||t-t_i||)是权函数，代入式(13)得到时变可靠性PCE：

$R(\mathrm{\&zeta;},t)=\underset{j=0}{\overset{N_r-1}{\&Sigma;}}{\hat{c}}_{rj}\left(t\right){\mathrm{\&psi;}}_{Lj}\left(\mathrm{\&zeta;}\right)---\left(15\right)$

式中，为通过MLS估计得到的时变可靠性PCE系数。

5.根据权利要求1所述的一种基于混沌多项式展开的时变可靠性全局灵敏度分析方法，其特征在于：在步骤四中所述的“计算时变全局灵敏度的Sobol’指标”，具体计算步骤如下：

首先，根据时变可靠性PCE计算时变可靠性方差D(R(ζ,t))：

$D\left(R\right(\mathrm{\&zeta;},t\left)\right)=\underset{j=1}{\overset{N-1}{\&Sigma;}}{\hat{c}}_{rj}^2\left(t\right)E\left({\mathrm{\&psi;}}_{Lj}^2\right(\mathrm{\&zeta;}\left)\right)---\left(16\right)$

式中， $E\left({\mathrm{\&psi;}}_{Lj}^2\right)=\underset{k=1}{\overset{n}{\&Pi;}}E\left({\mathrm{\&phi;}}_{{\mathrm{L\&alpha;}}_k^j}^2\right)=\underset{k=1}{\overset{n}{\&Pi;}}\frac{2}{2{\mathrm{\&alpha;}}_k^j+1},\left(\right|{\mathrm{\&alpha;}}^j|=\underset{k=1}{\overset{j}{\&Sigma;}}{\mathrm{\&alpha;}}_k^j\&le;p);$

然后将式(15)的系数按照Sobol’分解的形式进行重组，通过重组后的PCE的系数直接计算，具体步骤如下：

步骤4.1PCE系数重组

系数重组按照如下规则进行：

(1)被加项不展开，将单独的各个变量的一次项、二次项……放在一堆；

(2)被加项不展开，将具有两个变量交互作用的被加项的一次项、二次项……放在一堆；

(3)被加项不展开，将具有三个变量交互作用的被加项的一次项、二次项……放在一堆；

重组后数学表达式为：

式(17)中，α＝(α₁,…,α_n)是一个整数序列序列，是的具体值，满足

步骤4.2Sobol’s指标

然后利用PCE基底的正交性，直接得到Sobol’s指标：