CN111177974B - 基于双层嵌套寻优和子集模拟的结构小失效概率计算方法 - Google Patents

Abstract

本发明提供一种基于双层嵌套优化支持向量机和子集模拟结构小失效概率计算方法，步骤如下：一：确定研究对象；二：确定影响关键件的随机变量，建立有限元模型；三：根据有限元模型求出其对应的响应；四：根据当前的试验设计，构建初始模型；五：采用有效集进行内层优化构建初始模型；六：通过多路径粒子群寻优，对支持向量机中惩罚因子和核函数两参数进行外层优化，得到优化后的支持向量机参数；七：根据优化的参数构建最后的支持向量机回归模型，并得到最终的极限状态方程；八：利用子集模拟法对最终的极限状态方程进行失效概率求解，得到最终的失效概率；本发明所述方法科学，工艺性好，具有广阔推广应用价值。

Description

基于双层嵌套寻优和子集模拟的结构小失效概率计算方法

技术领域

本发明名称为“基于双层嵌套寻优和子集模拟的结构小失效概率计算方法”，它是提供一种基于双层嵌套优化支持向量机和子集模拟结构小失效概率计算方法，是涉及一种结构小失效概率计算方法，涉及到失效模型建模以及小失效概率求解，属于机械结构小失效概率计算领域。

背景技术

随着全球先进制造技术的迅猛提升，以航空航天、医疗器械等机电产品为代表的的高端装备创新领域对产品的失效概率提出了很高的要求。与此同时，单纯依靠现场试验已经无法验证产品的失效概率。数据收集困难、样本量小、对机械结构失效的影响因素多和小失效概率成为现在结构安全设计分析的难点问题。

一次二阶矩、高阶矩等方法是目前分析结构小失效概率的经典方法，但是这些方法均建立在极限状态函数是显式方程的情况下，无法进行直接求解，且对非线性较高的问题计算精度不高。随着概率的降低，采用蒙特卡罗法对模型进行小失效概率求解抽样次数往往需要达到10⁷这一量级，造成了很大的计算成本，很难被高端装备所接受。

作为在小失效概率计算中广泛应用的数值模拟方法，重要抽样法将抽样密度函数的抽样中心移到设计验算点，可以使更多的样本点集中在失效域，提高了抽样的效率。但在未知验算点的情况下，只能通过蒙特卡罗或解析近似法求解，这极大地限制了它们的使用范围。在建模方面，目前响应面法在高维非线性程度较高的情况下具有一定的局限性。而支持向量机通常在小样本和高维非线性的回归问题上表现较好，通过少量的训练样本就能得到较小的误差。

综上所述，针对目前高端装备领域的结构样本量小的问题，本发明提出了一种双层嵌套寻优支持向量机进行关键件的失效建模，接着用子集模拟方法计算失效概率。最后用蒙特卡罗算法验证了该方法的精度、计算效率以及可行性，对于提高产品的性能具有一定的工程实用价值。

发明内容

(一)本发明的目的

本发明方法利用支持向量机法建立可靠性模型，拟合隐式极限状态函数，通过多路径寻优和有效集法对支持向量机进行双层嵌套寻优，扩大本地搜索区域，在局部搜索和全局搜索设定不同路径，提高建模的精度和效率。在得到极限状态函数后，利用子集模拟通过引入合理的中间失效事件，将小失效概率表达为一系列较大的条件失效概率的乘积，而较大的条件失效概率可利用马尔可夫链模拟的条件样本点来高效估计，从而提高了小失效概率的概率的效率，从而有效的解决了高维、小失效概率事件以及隐式、非线性极限状态函数的概率求解问题，可以有效求解10^-4～10^-9下的小失效概率求解问题。

(二)技术方案

为了实现上述目的，本发明采用的技术方案是：

本发明一种基于双层嵌套寻优和子集模拟的结构小失效概率计算方法，即一种基于双层嵌套优化支持向量机和子集模拟结构小失效概率计算方法，该方法具体步骤为：

步骤一：确定研究对象和其关键失效模式；

步骤二：确定影响关键件失效模式的随机变量及其分布，并建立有限元模型：在此给出了典型关键件小失效概率计算中几种常见的特征和载荷/环境变量的概率表征方法；随机变量一般可分为随机变量(简单量)、随机场、随机过程以及随机场和随机过程的结合；随机因素很多，如材料性能不均匀、加工公差、表面完整性的差异、载荷的变化等，可以根据随机因素来源的不同，将这些随机变量分成不同的类型；根据关键件设计分析的特点，选取结构几何类随机变量、加工工艺类随机变量、载荷类随机变量和材料参数类随机变量该4类随机变量进行研究；在确定好相关变量的随机统计特性(变量的分布类型和均值)后，建立关键件的有限元模型；

步骤三：选择合适的试验设计抽取样本点，并根据有限元模型求出其对应的响应：

采用拉丁超立方抽样法抽取样本点，建立有限元模型，得到各样本点响应，组成样本集；

步骤四：根据当前的试验设计，利用支持向量机构建初始模型：

1.寻找一个核函数K(s,t)使得K(x_i,x_j)＝<φ(x_i),φ(x_j)>，高维情况下多采用高斯核函数；

式中，K(s,t)是核函数，φ(x_i)和φ(x_j)是输入变量的参数；

2.求优化问题

的解α_i,

式中：C是惩罚因子，C＞0，核函数K(x_i·x_j)＝Φ(x_i)·Φ(x_j)，当α_i不为0时，其对应的样本点即支持向量，ε表示精度；

步骤五：针对上述最优化问题采用有效集进行内层优化从而构建初始模型：

假设求解

其中

为核函数矩阵；

为精度控制；

为拉格朗日因子；

1.选取初值：给定初始可行点x₀∈Rⁿ，令k＝0；

2.解子问题：确定相应的有效集S_k＝E∪I(x_k)，求解子问题

得极小点d_k和拉格朗日乘子向量λ_k；若d_k≠0转步骤五中的4；否则，转步骤五中的3；

3.检验终止准则：计算拉格朗日乘子

λ_k＝B_kg_k

其中

g_k＝Hx_k+c,

令

若(λ_k)_t≥0，则x_k是全局极小点，停算；否则，若(λ_k)_t＜0，则令S_k:＝S_k\{t}，转步骤五中的2；

式中，λ_k为拉格朗日乘子向量，d_k为所得的极小点，x_k是全局极小点；

4.确定步长α_k.令

其中

令x_k+1＝x_k+α_kd_k.

式中，α_k为确定好的步长；

5.若α_k＝1，则令S_k+1:＝S_k；否则，若α_k＜1，则令S_k+1＝S_k∪{j_k}，其中j_k满足

6.令k＝k+1，转步骤五中的1；

7.最后输出优化后的模型；

步骤六：通过多路径粒子群寻优，对支持向量机中惩罚因子C和核函数σ两参数进行外层优化，得到优化后的支持向量机参数；

对支持向量机回归结果影响较大的参数有惩罚因子C和核函数σ，惩罚因子C是对错分的样本的惩罚程度的控制，越大表示惩罚越重，但其泛化能力也会同时降低；核函数σ是核函数的宽度参数，表示对径向范围的控制；合适的惩罚因子C和核函数σ对支持向量机回归性能有决定性影响；因此采用多路径粒子群寻优对这两参数进行寻优；

1.生成和初始化PSO及其各个参数(如粒子个数，惯性权重，学习因子，最大迭代次数等)，并求出适应度值，确定起始个体极值p_i和全局极值p_g；

2.将粒子群随机分成多组(本专利采用四组粒子群)，自定义每一组的c₁,c₂迭代路径(c₁控制了“自我认知”部分，即粒子自身值钱的飞行经验对之后飞行方向的影响，c₂控制了“社会认知”部分，即种群中所有粒子的飞行经验对每个粒子之后飞行方向的影响；自定义 c₁,c₂的迭代路径能够有效的平衡局部搜索和全局搜索的能力)；

3.通过多路径的循环迭代更新x、p_i和p_g：

①确定惯性权重w的取值(当w不是常数时)；

②根据

更新粒子的速度

若速度中的某一维超过了V_max，则取为V_max；

③根据

更新自变量x，若x的取值超过其定义域，则在其定义域内重新初始化；

④每组粒子群根据c₁,c₂自定义的迭代路径按照不同路径更新c₁,c₂；

⑤求得适应度值，通过比较更新个体极值p_i和全局极值p_g；

4.判断是否满足终止条件，若不满足，则转入步骤六中的3，若满足，则输出优化后的支持向量机参数，此步骤结束；

式中，

表示粒子当前的速度，

表示粒子当前位置，p_i表示个体极值，p_g表示群体极值，V_max为最大速度，c₁控制了“自我认知”部分，即粒子自身值钱的飞行经验对之后飞行方向的影响，c₂控制了“社会认知”部分，即种群中所有粒子的飞行经验对每个粒子之后飞行方向的影响；

步骤七：根据优化的惩罚因子C和核函数σ，构建最后的支持向量机回归模型，并得到最终的极限状态方程；

1.计算出最终的

2.构造非线性极限状态函数即最终的极限状态方程

式中，b为极限状态方程的参数变量，ε表示精度，K(x_i·x_j)表示核函数，当α_i不为0时，其对应的样本点即支持向量；

步骤八：利用子集模拟法对最终的极限状态方程进行失效概率求解，得到最终的失效概率；

1.用直接Monte Carlo模拟法产生N个服从联合概率密度函数为f_X(x)的相互独立的样本

2.通过功能函数g(x)得到这N个样本点对应的响应值

把这N个响应值从小到大排序，记为

取第Np₀个响应值作为中间事件 F₁＝{x:g(x)≤b₁}的临界值b₁，即

同时可知F₁区域失效概率的估计值 P₁＝P(F₁)＝p₀；

3.将落在F_i-1(i＝2,3,…,m)域内的Np₀个样本作为初始马尔可夫样本点，利用MCMC模拟N个服从密度函数q(x|F_i-1)的条件样本点

4.通过功能函数g(x)得到这N个条件样本对应的响应值

并对响应值进行从小到大排序，记为

取第Np₀个响应值作为中间事件F_i＝{x:g(x)≤b_i} 的临界值b_i，即

同时得到F_i-1发生的条件下F_i的条件失效概率的估计值 P_i＝P(F_i|F_i-1)＝p₀和F_i区域的失效概率估计值

5.重复步骤八中的3和4的过程，直到某一层(记为m层)的功能函数值从小到大排序后的第Np₀个响应值

值小于0，则令b_m＝0，自动分层结束；统计服从密度函数q(ε|F_m-1) 的条件样本点中落入失效域F中的个数N_f，则条件失效概率的估计值 P_m＝P(F_m|F_m-1)＝N_f/N；

6.分层结束后可得到失效概率估计值P_f为

式中，f_X(x)为联合概率密度函数，

为N个样本点对应的响应值，p₀为条件概率，N_f为落入失效域F中的个数，q(ε|F_m-1)为条件密度函数。

(三)本发明方法的优点和积极效果在于：

1)本发明通过双层嵌套优化来减少支持向量机建模的波动性；

2)多路径寻优可以扩大本地搜索空间的范围，平衡局部搜索和全局搜索，增加搜索效率和精度；

3)基于支持向量机的架构可以有效处理小子样情况下的建模问题，预测效率高效。

4)本发明给出了一种基于小样本情况下对发动机关键件建模和小失效概率求解的新方法。

5)本方法可以有效求解10^-4～10^-9下的小失效概率求解问题。

6)本方法同样适用于极限状态函数高度非线性的情况。

7)本发明所述方法科学，工艺性好，具有广阔推广应用价值。

附图说明

图1为支持向量机升维示意图。

图2为多路径粒子群寻优与与常规粒子群方法(PSO)、变异粒子群方法(MPSO)、协同进化粒子群方法(TACPSO)、约束粒子群方法(IPSO)进行比较图。

图3(a)、(b)、(c)为子集模拟示意图。

图4为蒙特卡罗法与子集模拟失效概率对比图。

图5为本发明所述方法流程图。

图中序号、符号、代号说明如下：

图2中，PSO为常规粒子群方法，MPSO为变异粒子群方法，TACPSO为协同进化粒子群方法，IPSO为约束粒子群方法。

图3中，F₁是自动抽样的第一层，b₁是第一层对应的临界值。F₂是自动抽样的第二层， b₂是第二层对应的临界值。

图4中，P_f是失效概率，b是阈值设定，MCS是Monte Carlo模拟法，SS是本文中所提的方法。

具体实施方式

见图1-图4，本发明一种基于双层嵌套寻优和子集模拟的结构小失效概率计算方法，即一种基于双层嵌套优化支持向量机和子集模拟结构小失效概率计算方法，本发明方法利用支持向量机法建立可靠性模型，拟合隐式极限状态函数，通过多路径寻优和有效集法对支持向量机进行双层嵌套寻优，扩大本地搜索区域，在局部搜索和全局搜索设定不同路径，提高建模的精度和效率；在得到极限状态函数后，利用子集模拟通过引入合理的中间失效事件，将小失效概率表达为一系列较大的条件失效概率的乘积，而较大的条件失效概率可利用马尔可夫链模拟的条件样本点来高效估计，从而提高了小失效概率的概率的效率，从而有效的解决了高维、小失效概率事件以及隐式、非线性极限状态函数的概率求解问题。

本发明一种基于双层嵌套寻优和子集模拟的结构小失效概率计算方法，即一种基于双层嵌套优化支持向量机和子集模拟结构小失效概率计算方法，见图5所示，该方法具体步骤为：

步骤一：确定研究对象和其关键失效模式：

发动机是复杂产品的典型代表，在此选用某型发动机低压压气机三级轮盘为关键件；其主要失效模式为低循环疲劳断裂；

步骤二：确定影响关键件失效模式的随机变量及其分布，并建立有限元模型：在此给出了典型关键件小失效概率计算中几种常见的特征和载荷/环境变量的概率表征方法；随机变量一般可分为随机变量(简单量)、随机场、随机过程以及随机场和随机过程的结合。随机因素很多，如材料性能不均匀、加工公差、表面完整性的差异、载荷的变化等，可以根据随机因素来源的不同，将这些随机变量分成不同的类型。根据关键件设计分析的特点，选取结构几何类随机变量、加工工艺类随机变量、载荷类随机变量和材料参数类随机变量4类随机变量进行研究。在确定好相关变量的随机统计特性(变量的分布类型和均值)后，建立关键件的有限元模型；

步骤三：选择合适的试验设计抽取样本点，并根据有限元模型求出其对应的响应：

采用拉丁超立方抽样法抽取样本点，建立有限元模型，得到各样本点响应，组成样本集；

步骤四：根据当前的试验设计，利用支持向量机构建初始模型：

1.寻找一个核函数K(s,t)使得K(x_i,x_j)＝<φ(x_i),φ(x_j)>，高维情况下多采用高斯核函数；

式中，K(s,t)是核函数，φ(x_i)和φ(x_j)是输入变量的参数；

2.求优化问题

的解α_i,

式中：C是惩罚因子，C＞0，核函数K(x_i·x_j)＝Φ(x_i)·Φ(x_j)，当α_i不为0时，其对应的样本点即支持向量，ε表示精度；

步骤五：针对上述最优化问题采用有效集进行内层优化从而构建初始模型：

假设求解

其中

为核函数矩阵；

为精度控制；

为拉格朗日因子；

1.选取初值：给定初始可行点x₀∈Rⁿ，令k＝0；

2.解子问题：确定相应的有效集S_k＝E∪I(x_k)，求解子问题

得极小点d_k和拉格朗日乘子向量λ_k；若d_k≠0转步骤五中的4；否则，转步骤五中的3；

3.检验终止准则：计算拉格朗日乘子

λ_k＝B_kg_k

其中

g_k＝Hx_k+c,

令

若(λ_k)_t≥0，则x_k是全局极小点，停算；否则，若(λ_k)_t＜0，则令S_k:＝S_k\{t}，转步骤五中的2；

式中，λ_k为拉格朗日乘子向量，d_k为所得的极小点，x_k是全局极小点；

4.确定步长α_k.令

其中

令x_k+1＝x_k+α_kd_k.

式中，α_k为确定好的步长；

5.若α_k＝1，则令S_k+1:＝S_k；否则，若α_k＜1，则令S_k+1＝S_k∪{j_k}，其中j_k满足

6.令k＝k+1，转步骤五中的1；

7.最后输出优化后的模型；

步骤六：通过多路径粒子群寻优，对支持向量机中惩罚因子C和核函数σ两参数进行外层优化，得到优化后的支持向量机参数；

对支持向量机回归结果影响较大的参数有惩罚因子C和核函数σ，惩罚因子C是对错分的样本的惩罚程度的控制，越大表示惩罚越重，但其泛化能力也会同时降低；核函数σ是核函数的宽度参数，表示对径向范围的控制；合适的惩罚因子C和核函数σ对支持向量机回归性能有决定性影响；因此采用多路径粒子群寻优对这两参数进行寻优；

1.生成和初始化PSO及其各个参数(如粒子个数，惯性权重，学习因子，最大迭代次数等)，并求出适应度值，确定起始个体极值p_i和全局极值p_g；

2.将粒子群随机分成多组(本专利采用四组粒子群)，自定义每一组的c₁,c₂迭代路径(c₁控制了“自我认知”部分，即粒子自身值钱的飞行经验对之后飞行方向的影响，c₂控制了“社会认知”部分，即种群中所有粒子的飞行经验对每个粒子之后飞行方向的影响；自定义 c₁,c₂的迭代路径能够有效的平衡局部搜索和全局搜索的能力)；

3.通过多路径的循环迭代更新x、p_i和p_g：

①确定惯性权重w的取值(当w不是常数时)；

②根据

更新粒子的速度

若速度中的某一维超过了V_max，则取为V_max；

③根据

更新自变量x，若x的取值超过其定义域，则在其定义域内重新初始化；

④每组粒子群根据c₁,c₂自定义的迭代路径按照不同路径更新c₁,c₂；

⑤求得适应度值，通过比较更新个体极值p_i和全局极值p_g；

4.判断是否满足终止条件，若不满足，则转入步骤六中的3，若满足，则输出优化后的支持向量机参数，此步骤结束；

式中，

表示粒子当前的速度，

表示粒子当前位置，p_i表示个体极值，p_g表示群体极值，V_max为最大速度，c₁控制了“自我认知”部分，即粒子自身值钱的飞行经验对之后飞行方向的影响，c₂控制了“社会认知”部分，即种群中所有粒子的飞行经验对每个粒子之后飞行方向的影响；

步骤七：根据优化的惩罚因子C和核函数σ，构建最后的支持向量机回归模型，并得到最终的极限状态方程；

1.计算出最终的

2.构造非线性极限状态函数即最终的极限状态方程

式中，b为极限状态方程的参数变量，ε表示精度，K(x_i·x_j)表示核函数，当α_i不为0时，其对应的样本点即支持向量；

步骤八：利用子集模拟法对最终的极限状态方程进行失效概率求解，得到最终的失效概率；

1.用直接Monte Carlo模拟法产生N个服从联合概率密度函数为f_X(x)的相互独立的样本

2.通过功能函数g(x)得到这N个样本点对应的响应值

把这N个响应值从小到大排序，记为

取第Np₀个响应值作为中间事件 F₁＝{x:g(x)≤b₁}的临界值b₁，即

同时可知F₁区域失效概率的估计值 P₁＝P(F₁)＝p₀；

3.将落在F_i-1(i＝2,3,…,m)域内的Np₀个样本作为初始马尔可夫样本点，利用MCMC模拟N个服从密度函数q(x|F_i-1)的条件样本点

4.通过功能函数g(x)得到这N个条件样本对应的响应值

并对响应值进行从小到大排序，记为

取第Np₀个响应值作为中间事件F_i＝{x:g(x)≤b_i} 的临界值b_i，即

同时得到F_i-1发生的条件下F_i的条件失效概率的估计值 P_i＝P(F_i|F_i-1)＝p₀和F_i区域的失效概率估计值

5.重复步骤八中的3和4的过程，直到某一层(记为m层)的功能函数值从小到大排序后的第Np₀个响应值

值小于0，则令b_m＝0，自动分层结束；统计服从密度函数q(ε|F_m-1) 的条件样本点中落入失效域F中的个数N_f，则条件失效概率的估计值 P_m＝P(F_m|F_m-1)＝N_f/N；

6.分层结束后可得到失效概率估计值P_f为

式中，f_X(x)为联合概率密度函数，

为N个样本点对应的响应值，p₀为条件概率，N_f为落入失效域F中的个数，q(ε|F_m-1)为条件密度函数。

下面结合附图和实施例子对本发明做进一步说明。

任务描述：

轮盘是航空发动机的主要组成部分，是发动机的关键转动部件之一。常见的导致轮盘失效的因素主要是低循环疲劳和蠕变等。低压压气机轮盘三维实体模型，该模型共有37个榫槽，而且在榫槽处存在局部应力集中。有限元建模时轮盘材料参数设置为TC11，半空心销为 3Cr13。以某型低压压气机三级轮盘作为研究对象进行案例验证。

由式

可知压气机轮盘在0～29189.6N的脉动载荷下的疲劳强度与寿命曲线方程为

lgN_life＝lgK-mlg(S_max-703.84)

式中：N_life为疲劳载荷循环次数，σ_a为循环疲劳应力幅值，S_max为应力场强，TC11的材料参数K＝6.2784×10¹⁵，m＝4.736。

构建极限状态函数，有

式中：Z＝[K m S_max]，

为压气机轮盘性能强度，

为性能退化量，n_max为在S_max下的循环次数。本文中只考虑单个疲劳应力水平下的概率风险分析，即

压气机轮盘各参数详见表。

轮盘有限元模型。

采用拉丁超立方抽样法抽取样本点，建立有限元模型，得到各样本点响应，组成样本集。

根据当前的试验设计，利用支持向量机构建初始模型：

为了能够实现非线性映射，同时考虑参数的数量较多，模型十分复杂，在此选用高斯核函数(RBF核函数)，如图1所示。

针对上述最优化问题采用有效集进行内层优化从而构建初始模型：

设定了四组粒子群，每组粒子群c₁,c₂具有不同的迭代路径，具体路径的函数方程式如下所示。

第一组：c₁＝(-2*(t³)/M³)+2.6

c₂＝(1/M)*t+1.2

第二组：c₁＝(-2*(t³)/M³)+2.6

c₂＝(2*(t³)/M³)+0.6

第三组：c₁＝(-2.13/M)*t+2.57

c₂＝(1/M)*t+1.2

第四组：c₁＝(-2.13/M)*t+2.57

c₂＝(2*(t³)/M³)+0.6

最终结果如下所示

pg＝2.935827145680287 1.245925317150316

Fv_last＝ 0.54478769385023825901960549453491

根据优化的惩罚因子C和核函数σ，构建最后的支持向量机回归模型，并得到最终的极限状态方程；

F(x)

＝2.2856e-42*exp(-1.16*x1^2)*exp(-1.16*x2^2)*exp(-1.16*x3^2)*exp(-1.16*x4^2)*ex p(-1.16*x5^2)*exp(-1.16*x6^2)*exp(-1.16*x7^2)*...

*exp(-1.16*x4^2)*exp(-1.16*x5^2)*exp(-1.16*x6^2)

*exp(2.2387*x7)*exp(9.4349*x3)*exp(4.9082*x5)*exp(3.3225*x4)*exp(7.0367*x8)*ex p(3.6237*x1)*exp(15.354*x6)*exp(7.7585*x2)；

接着，用子集模拟对上述建立好的极限状态方程求解。

计算结果如下表所示：

蒙特卡罗和子集模拟计算失效概率的对比如图4所示。通过对比可以发现，本文所提方法较传统的方法具有更高的精度和效率，大大节约了计算成本，验证了其方法的有效性。失效概率曲线图表明该可靠度计算方法与MSC方法在失效概率的计算上差别较小，其计算结果可以认为是准确的。

本文所提方法有效的解决了高维、小失效概率事件以及隐式、非线性极限状态函数的概率求解问题，可以有效求解10^-4～10^-9下的小失效概率求解问题，而且大大减少了计算成本，提高了计算效率。

Claims (1)

1.一种基于双层嵌套寻优和子集模拟的结构小失效概率计算方法，其特征在于：包括如下步骤：

步骤一：确定研究对象和其关键失效模式：

选用某型发动机低压压气机三级轮盘为关键件；共有37个榫槽，而且在榫槽处存在局部应力集中；有限元建模时轮盘材料参数设置为TC11，半空心销为3Cr13，失效模式为低循环疲劳断裂；

由式

得到压气机轮盘在0～29189.6N的脉动载荷下的疲劳强度与寿命曲线方程为

lgN_life＝lgK-mlg(S_max-703.84)

式中：N_life为疲劳载荷循环次数，σ_a为循环疲劳应力幅值，S_max为应力场强，TC11的材料参数K＝6.2784×10¹⁵，m＝4.736；

构建极限状态函数，有

式中：Z＝[K m S_max]，

为压气机轮盘性能强度，

为性能退化量，n_max为在S_max下的循环次数；现只考虑单个疲劳应力水平下的概率风险分析，即

压气机轮盘各参数详见表；

步骤二：确定影响关键件失效模式的随机变量及其分布，并建立有限元模型：

选取结构几何类随机变量、加工工艺类随机变量、载荷类随机变量和材料参数类随机变量4类随机变量；在确定好相关变量的随机统计特性后，建立关键件的有限元模型；

步骤三：选择合适的试验设计抽取样本点，并根据有限元模型求出其对应的响应：采用拉丁超立方抽样法抽取样本点，建立有限元模型，得到各样本点响应，组成样本集；

轮盘有限元模型；

采用拉丁超立方抽样法抽取样本点，建立有限元模型，得到各样本点响应，组成样本集；

步骤四：根据当前的试验设计，利用支持向量机构建初始模型：

寻找一个核函数K(s,t)使得K(x_i,x_j)＝<φ(x_i),φ(x_j)>，高维情况下多采用高斯核函数；

式中，K(s,t)是核函数，φ(x_i)和φ(x_j)是输入变量的参数；

求优化问题

的解α_i,

式中：C是惩罚因子，C＞0，核函数K(x_i·x_j)＝Φ(x_i)·Φ(x_j)，当α_i不为0时，其对应的样本点即支持向量，ε表示精度；

步骤五：针对最优化问题采用有效集进行内层优化从而构建初始模型：

设求解

其中

为核函数矩阵；

为精度控制；

为拉格朗日因子；

5.1选取初值：给定初始可行点x₀∈Rⁿ，令k＝0；

5.2解子问题：确定相应的有效集S_k＝E∪I(x_k)，求解子问题

得极小点d_k和拉格朗日乘子向量λ_k；若d_k≠0转步骤5.4；否则，转步骤5.3；

5.3检验终止准则：计算拉格朗日乘子

λ_k＝B_kg_k

其中

g_k＝Hx_k+c,

令

若(λ_k)_t≥0，则x_k是全局极小点，停算；否则，若(λ_k)_t＜0，则令S_k:＝S_k\{t}，转步骤5.2；

式中，λ_k为拉格朗日乘子向量，d_k为所得的极小点，x_k是全局极小点；

5.4确定步长α_k.令

其中

令x_k+1＝x_k+α_kd_k

式中，α_k为确定好的步长；

5.5若α_k＝1，则令S_k+1:＝S_k；否则，若α_k＜1，则令S_k+1＝S_k∪{j_k}，其中j_k满足

5.6令k＝k+1，转步骤5.1；

5.7最后输出优化后的模型；

步骤六：通过多路径粒子群寻优，对支持向量机中惩罚因子C和核函数σ两参数进行外层优化，得到优化后的支持向量机参数；

采用多路径粒子群寻优对这两参数进行寻优；

6.1.生成和初始化PSO及其各个参数，包括粒子个数，惯性权重，学习因子，最大迭代次数，并求出适应度值，确定起始个体极值p_i和全局极值p_g；

6.2.将粒子群随机分成多组，自定义每一组的c₁,c₂迭代路径；

6.3.通过多路径的循环迭代更新x、p_i和p_g：

①确定惯性权重w的取值，当w不是常数时；

②根据

更新粒子的速度

若速度中的某一维超过了V_max，则取为V_max；

③根据

更新自变量x，若x的取值超过其定义域，则在其定义域内重新初始化；

④每组粒子群根据c₁,c₂自定义的迭代路径按照不同路径更新c₁,c₂；

⑤求得适应度值，通过比较更新个体极值p_i和全局极值p_g；

6.4.判断是否满足终止条件，若不满足，则转入步骤6.3，若满足，则输出优化后的支持向量机参数，此步骤结束；

式中，

表示粒子当前的速度，

表示粒子当前位置，p_i表示个体极值，p_g表示群体极值，V_max为最大速度，c₁控制了“自我认知”部分，即粒子自身值钱的飞行经验对之后飞行方向的影响，c₂控制了“社会认知”部分，即种群中所有粒子的飞行经验对每个粒子之后飞行方向的影响；

步骤七：根据优化的惩罚因子C和核函数σ，构建最后的支持向量机回归模型，并得到最终的极限状态方程；

7.1.计算出最终的

7.2.构造非线性极限状态函数即最终的极限状态方程

式中，b为极限状态方程的参数变量，ε表示精度，K(x_i·x_j)表示核函数，当α_i不为0时，其对应的样本点即支持向量；

步骤八：利用子集模拟法对最终的极限状态方程进行失效概率求解，得到最终的失效概率；

8.1.用直接Monte Carlo模拟法产生N个服从联合概率密度函数为f_X(x)的相互独立的样本

8.2.通过功能函数g(x)得到这N个样本点对应的响应值

把N个响应值从小到大排序，记为

取第Np₀个响应值作为中间事件F₁＝{x:g(x)≤b₁}的临界值b₁，即

同时得出F₁区域失效概率的估计值P₁＝P(F₁)＝p₀；

8.3.将落在F_i-1(i＝2,3,…,m)域内的Np₀个样本作为初始马尔可夫样本点，利用MCMC模拟N个服从密度函数q(x|F_i-1)的条件样本点

8.4.通过功能函数g(x)得到这N个条件样本对应的响应值

并对响应值进行从小到大排序，记为

取第Np₀个响应值作为中间事件

F_i＝{x:g(x)≤b_i}的临界值b_i，即

同时得到F_i-1发生的条件下F_i的条件失效概率的估计值P_i＝P(F_i|F_i-1)＝p₀和F_i区域的失效概率估计值

8.5.重复步骤8.3和8.4的过程，直到某一层的功能函数值从小到大排序后的第Np₀个响应值

值小于0，则令b_m＝0，自动分层结束；统计服从密度函数q(ε|F_m-1)的条件样本点中落入失效域F中的个数N_f，则条件失效概率的估计值P_m＝P(F_m|F_m-1)＝N_f/N；

8.6.分层结束后得到失效概率估计值P_f为

式中，f_X(x)为联合概率密度函数，

为N个样本点对应的响应值，p₀为条件概率，N_f为落入失效域F中的个数，q(ε|F_m-1)为条件密度函数；

其中，在步骤6和步骤7中，设定了四组粒子群，每组粒子群c₁,c₂具有不同的迭代路径，具体路径的函数方程式如下所示；

第一组：c₁＝(-2*(t³)/M³)+2.6

c₂＝(1/M)*t+1.2

第二组：c₁＝(-2*(t³)/M³)+2.6

c₂＝(2*(t³)/M³)+0.6

第三组：c₁＝(-2.13/M)*t+2.57

c₂＝(1/M)*t+1.2

第四组：c₁＝(-2.13/M)*t+2.57

c₂＝(2*(t³)/M³)+0.6

最终结果如下所示：

pg＝2.935827145680287；1.245925317150316

Fv_last＝0.54478769385023825901960549453491

根据优化的惩罚因子C和核函数σ，构建最后的支持向量机回归模型，并得到最终的极限状态方程；接着，用子集模拟对上述建立好的极限状态方程求解。

Images