基于粒子群和支持向量机的控制系统执行器故障诊断方法

摘要

本发明公开了一种基于多群体协同混沌模拟退火粒子群优化算法-支持向量机(Multi-swarm Cooperative Chaos Simulated Annealing Particle Swarm Optimization-Support Vector Machine，MCCSAPSO-SVM)的控制系统执行器故障诊断方法。通过联合降噪和改进经验模态分解(Empirical Mode Decomposition，EMD)方法，对采集到的执行器输出信号进行特征提取；利用MCCSAPSO优化支持向量机结构参数；采用组合核函数保证了SVM良好的泛化能力与学能力；利用训练数据构建偏二叉树SVM，把一个复杂的多分类问题转化为若干个二分类问题。本发明方法对控制系统易获得的常见状态信号进行处理，通过偏二叉树SVM的输出，能够实时有效地判断出控制系统执行器是否发生故障，并且能够在执行器发生故障时较为准确地确定出故障类型。本发明用于高精度控制系统的实时故障诊断。

主权项

一种基于多群体协同混沌模拟退火粒子群优化算法一支持向量机(MCCSAPSO‑SVM)的控制系统执行器故障诊断方法，其特征在于：通过联合降噪和改进经验模态分解(EMD)方法，对采集到的执行器输出信号进行降噪处理和特征提取；利用多群体协同混沌模拟退火粒子群优化算法(MCCSAPSO)优化支持向量机结构参数；采用组合核函数保证了支持向量机良好的泛化能力与学习能力；利用训练数据构建一种把一个复杂的多分类问题转化为若干个二分类问题的偏二叉树支持向量机，降低了计算量，提高了诊断的实时性；通过偏二叉树支持向量机的输出，能够实时有效地判断出控制系统执行器是否发生故障，并且能够在执行器发生故障时较为准确地确定出故障类型，包括如下具体步骤：步骤1)将离散含噪信号f(k)(信号长度为N)先进行中值滤波处理得到滤除可能的脉冲噪声；步骤2)对中值滤波处理后的进行改进小波阀值降噪处理：步骤2.1)将进行第j层(从j＝1开始)小波分解，其中j的确定方式为：对小波分解得到的高频细节系数d_j(k)进行式(1)自相关系数λ计算，若满足自由度l的χ²分布，则对d_j继续进行j+1层分解，直到不满足χ²分布为止；式中，d_j(k)高频细节系数，为d_j(k)的平均值；${\lambda }_i=\frac{\underset{k}{\overset{N-i}{\Sigma }}[d_j\left(k\right)-{\bar{d}}_j][d_j(k+i)-{\bar{d}}_j]}{\underset{k}{\overset{m-i}{\Sigma }}{[d_j\left(k\right)-{\bar{d}}_j]}^2},i=1,2,...,l---\left(1\right)$步骤2.2)将分解得到的小波系数W_j(k)带入式(2)加以处理；式中，W为信号小波变换后的小波系数W_j(k)，μ(μ＞0)，v(v＞1)，p(p∈[0，1])，q(q≥0)均为可调参数，W_δ为经过阀值降噪后的小波系数，为阀值，σ为噪声标准差；$W_{\delta }=\left\{\begin{array}{cc}\mathrm{sgn}\left(W\right)\left(\right|W|-\frac{2p}{1+e^{q\left(\right|W|-\delta )}}\delta ),& \left|W\right|\ge \delta \\ W·v^{-{\left(\mu \right|W|-\delta )}^2},& \left|W\right|\le \delta \end{array}\right.---\left(2\right)$步骤23)对处理后的各小波系数进行小波逆变换，得到降噪后的真实信号的估计步骤3)对联合降噪后的信号进行EMD分解，通过比较各本征模态函数(IMF)与信号的相关系数ρ与预设阀值ζ的大小，剔除max(ρ)后的所有小于阀值的IMF即“伪分量”；其中，ρ通过式(3)计算，σ[·]为标准差，ζ＝0.1max(ρ)；$\rho =\frac{\underset{k=1}{\overset{N}{\Sigma }}\hat{s}\left(k\right)·{\mathrm{imf}}_i\left(k\right)}{\sigma \left[\hat{s}\left(k\right)\right]·\sigma \left[{\mathrm{imf}}_i\left(k\right)\right]}---\left(3\right)$步骤4)对剔除“伪分量”后的η个IMF分量，计算各IMF分量c_i(k)的能量$E_i=\underset{k=1}{\overset{N}{\Sigma }}{c_i}^2\left(k\right)$(i＝1，2，...，η)，则提取出的特征向量为步骤5)利用多群体协同混沌模拟退火粒子群优化算法(MCCSAPSO)优化组合核函数支持向量机结构参数：步骤5.1)构建粒子群优化算法迭代公式(4)，该迭代公式能够使得粒子群在迭代初期具备较强的搜索能力，能搜索到较大的解空间，并具有不断搜索新区域的能力，而在迭代后期，算法快速收敛到全局最优区域；式中，c₁、c₂为学习因子，r₁、r₂为介于[0，1]之间的随机数，p_ibest＝(p_i1，p_i2，...，p_iD)为粒子历史最优适应值，g_best＝(g₁，g₂，...，g_D)为全局最优适应值，w为惯性权重，t_max为最大迭代次数，w_min，w_max是惯性权重w的最值，通过调节引入的收敛因子δ，可以加快粒子群的收敛速度；$\left\{\begin{array}{c}{v}_i^{t+1}=\delta ({\mathrm{wv}}_i^t+c_1{r}_1(p_{ibest}-x_i^t)+c_2{r}_2(g_{best}-x_i^t))\\ x_i^{t+1}=x_i^t+v_i^{t+1}\\ w=w_{\min }+\left[(t_{\max }-t)(w_{\max }-w_{\min })\right]/t_{max}\\ \delta =2/|2-(c_1+c_2)-\sqrt{{(c_1+c_2)}^2-4(c_1+c_2)}\end{array}\right.---\left(4\right)$步骤5.2)初始化：分别对主粒子群和从粒子群的粒子群规模L、维度D、粒子i的位置x_i＝(x_i1，x_i2，...，x_iD)、速度v_i＝(v_i1，v_i2，...，v_iD)、惯性权重w、学习因子c₁，c₂、当前温度T、结束温度T₀、退火速度K_T、最大迭代次数t_max初始化；步骤5.3)从群更新：计算从粒子群每个粒子的当前适应值Ψ_i，并和自身最优适应值p_ibest比较，更新自身最优适应值p_ibest和全局最优适应值g_best；根据式(4)分别各自更新粒子的速度和位置；其中，选取适应度函数Ψ为式(5)，即适应度函数是以最大化形式进行目标优化的分类器正确率的倍数，如式(6)；式中，为SVM_j正确分类数，为SVM_j错误分类数；Ψ＝100r 0.1≤r≤1      (5)$max\Psi (C,\sigma ,d,\rho )=100r=100\times \underset{j=1}{\overset{Q}{\Sigma }}\frac{N_j^c}{N_j^c+N_j^w}---\left(6\right)$步骤5.4)混沌优化：先将g_best映射到Logistic方程[0，1]定义域上，如式(7)：R＝(g_best‑R_min)/(R_max‑R_min)      (7)式中，R_max、R_min为g_best的上下限。接着对当前序列中的最优粒子位置和速度进行P次迭代后产生混沌序列为R＝(R₁，R₂，…，R_P)，通过式(8)可以得到如式(9)所示的最优粒子序列；g_best，p＝R_min+R_n(R_max‑R_min)，n＝1，2，…，P   (8)$g_{best}^{*}=(g_{best,1}^{*},g_{best,2}^{*},...,g_{best,p}^{*})---\left(9\right)$步骤5.5)主群更新：在每一代主群选取从群体中最好的粒子，并根据从群的经验进行状态更新，其速度和位置更新方程为式(10)：$\left\{\begin{array}{c}{v}_{Mi}^{t+1}=\delta ({\mathrm{wv}}_{Mi}^t+c_1{r}_1(p_{Mibest}-x_{Mi}^t)+{\mathrm{\phi c}}_2{r}_2(g_{Mbest}-x_{Mi}^t)+(1-\phi )c_3{r}_3(g_{Sbest}-x_{Si}^t))\\ x_{Mi}^{t+1}=x_{Mi}^t+v_{Mi}^{t+1}\end{array}\right.---\left(10\right)$式中，M为主群，S为从群，c₁为学习因子，r₃为介于[0，1]之间的随机数，φ为满足式(11)的迁移因子，g_Mbest和g_Sbest分别为主群和从群中的全局最优适应值，G_Mbest和G_Sbest分别为由g_Mbest和g_Sbest确定的适应值；$\phi =\left\{\begin{array}{cc}0,& G_{Sbest}<G_{Mbest}\\ 0.5,& G_{Sbest}=G_{Mbest}\\ 1,& G_{Sbest}>G_{Mbest}\end{array}\right.---\left(11\right)$步骤5.6)退火优化：计算每个粒子更新后的适应值及适应值变化量若ΔΨ_i＜0，或ΔΨ_i＞0时exp(‑ΔΨ/T)在区间[0，1]上，则进行降温操作T←K_TT，否则温度不变；步骤5.7)终止条件：当满足温度降至T₀，或达到最大迭代次数t_max时停止迭代，否则返回步骤5.3；步骤6)构造组合核函数支持向量机：步骤6.1)构造式(12)形式的组合核函数，式中，d为多项式核函数的阶数，σ为高斯核函数的核半径，ρ为权重；$K={\mathrm{\rho K}}_{poly}+(1-\rho )K_{RBF}=\rho {({\mathrm{xx}}_i+1)}^d+(1-\rho )\exp (-||x-x_i|{|}^2/{\sigma }^2)---\left(12\right)$步骤6.2)构造组合核函数支持向量机的决策函数可表示为式(13)，式中，{x_i，y_i}为训练样本，l是训练样本个数，K(x_i，x_j)是给定的核函数，α_i和b为根据训练样本得出的参数；$f\left(x\right)=\mathrm{sgn}\left(\underset{j=1}{\overset{l}{\Sigma }}{\alpha }_i{y}_i\right({\mathrm{\rho K}}_{poly}\left(x_i,x_j\right)+\left(1-\rho \right)K_{RBF}\left(x_j,x_j\right))+b)---\left(13\right)$步骤6.3)将组合核函数支持向量机的优化模型(14)转化为式(15)，式中，C为惩罚因子，ξ_i为松弛变量；$\begin{array}{c}\underset{\rho ,(1-\rho ),f_{poly},f_{RBF},b,{\xi }_i}{\min }\left(\frac{1}{\rho }\right|\left|f_{poly}\right|{|}^2)+\left(\frac{1}{(1-\rho )}\right|\left|f_{RBF}\right|{|}^2)+C\underset{i}{\overset{l}{\Sigma }}{\xi }_i\\ \begin{array}{cc}s.t.& y_i(f_{poly}+f_{RBF}+b)\ge 1-{\xi }_i\\ & {\xi }_i\ge 0,0\le \rho \le 1\end{array}\end{array}---\left(14\right)$$\begin{array}{c}\underset{\rho ,(1-\rho ),f_{poly},f_{RBF},b,{\xi }_i}{\min }\frac{1}{2}[\left(\frac{1}{\rho }\right|\left|f_{poly}\right|{|}^2)+\left(\frac{1}{(1-\rho )}\right|\left|f_{RBF}\right|{|}^2)]+C\underset{i}{\overset{l}{\Sigma }}{\xi }_i\\ \begin{array}{cc}s.t.& y_i(f_{poly}+f_{RBF}+b)\ge 1-{\xi }_i\\ & {\xi }_i\ge 0,0\le \rho \le 1\end{array}\end{array}---\left(15\right)$步骤6.4)利用Largrange乘数法和KKT条件，计算其对偶问题(16)；$\begin{array}{c}\underset{\alpha }{\max }\underset{i=1}{\overset{l}{\Sigma }}{\alpha }_i-\frac{1}{2}\underset{i=1}{\overset{l}{\Sigma }}\underset{j=1}{\overset{l}{\Sigma }}{\alpha }_i{\alpha }_j{y}_i{y}_j\left({\mathrm{\rho K}}_{poly}\right(x_i,x_j)+(1-\rho \left)K_{RBF}\right(x_i,x_j\left)\right)\\ \begin{array}{cc}s.t.& \underset{i=1}{\overset{l}{\Sigma }}{\alpha }_i{y}_i=0,0\le {\alpha }_i\le C,i=1,...,l\\ & 0\le \rho \le 1\end{array}\end{array}---\left(16\right)$步骤6.5)令解得最优解为则最优分类函数为式(17)；$f\left(x\right)=\mathrm{sgn}\left(\underset{i=1}{\overset{l}{\Sigma }}{\alpha }_i^{*}{y}_i\right({\mathrm{\rho K}}_{poly}\left(x_i,x_j\right)+\left(1-\rho \right)K_{RBF}\left(x_i,x_j\right))+b^{*})---\left(17\right)$步骤7)构建偏二叉树支持向量机分类器，针对含有正常状态和(Q‑1)种故障状态的诊断问题，需训练(Q‑1)个支持向量机分类器：步骤7.1)构建第一个分类器SVM₁是对正常状态样本与其余所有(Q‑1)种故障样本进行分类，将正常状态样本作为一类，分类器输出为+1，其余所有(Q‑1)种故障样本的分类器输出为‑1；步骤7.2)构建第二个分类器SVM₂是对故障1样本与剩余的其他(Q‑2)种故障样本进行分类，将故障1样本作为一类，分类器输出为+1，剩余的其他(Q‑2)种故障样本的分类器输出为‑1；步骤7.3)以此类推，直到第(Q‑1)个分类器SVM_(Q‑1)是对故障(Q‑2)样本与故障(Q‑1)样本进行分类，故障(Q‑2)样本分类器输出为+1，故障(Q‑1)样本分类器输出为‑1；步骤8)将测试数据应用于训练好的(Q‑1)个支持向量机分类器，根据各个SVM的输出值实时有效地判断出控制系统执行器是否发生故障，以及准确判定执行器发生的具体故障类型。