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

摘要

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

说明书

技术领域

本发明涉及一种基于多群体协同混沌模拟退火粒子群优化算法-支持向量机 (MCCSAPSO-SVM)的控制系统执行器故障诊断方法，属于控制系统信号处理与执行器故障诊 断技术领域。

背景技术

由于现代控制系统组成复杂，常常需要在不同环境条件下长时间、高负荷工作，这 就造成控制系统不可避免的会出现各类故障。特别是在航空航天、医疗、规模化机械生产等 领域，细微的故障有时会造成极其严重的经济损失和人员伤害的后果，因此对设备运行的 状态监测和故障诊断就成为了重要的研究课题。控制系统执行器的故障诊断可分为故障检 测和故障隔离，这两者是两种具有不同目的和方法的技术。故障检测的目的是判断系统当 前运行状态是否正常；故障隔离的目的是确定故障的类型和发生的位置。搭建一个快速、有 效、可靠的执行器故障诊断系统，其关键在于：1)准确获取能够反映控制系统当前运行状态 的特征信息；2)建立一种泛化能力强，能够处理小样本故障诊断的算法。故障特征提取的优 劣对诊断结果将起到至关重要的作用，经验模态分解(EmpiricalModeDecomposition， EMD)是一种非常适合应用于非线性、非平稳信号的特征提取方法。由于现场采集控制系统 状态信号常常存在较强的随机噪声，因此在控制系统中在应用经验模态分解前，必须先对 采样信号进行适当且必要的降噪处理。

故障诊断是一个典型的小样本问题，控制系统执行器发生故障具有突发性且难以 复现，因此可用的故障样本十分有限。基于统计学习理论的支持向量机在处理小样本、非线 性、高维问题的模式识别中具有极强的泛化能力，相较于神经网络具有所需样本少、结构更 明确、收敛速度快等优点。

在支持向量机的实际应用中，选择适合样本数据的有效核函数，能够增强支持向 量机决策函数的可解释性与鲁棒性。将多个单一核函数进行线性组合而成的组合核函数可 以弥补单一核函数对多源复杂样本数据可解释性不足的缺点。组合核函数在通过结合其构 成中各单一核函数的优势，得到性能更加优越的支持向量机的同时，仍然存在以下两个问 题：1)由于组合核函数的核参数较多，选取不当极易造成组合核函数的性能不及单一核函 数的性能；2)在已确定组合核函数所包含的单一核函数种类以及核参数的情况下，各单一 核函数所占权重的选择也会对组合核函数的性能造成影响。

通常将核函数的参数与惩罚因子等规则化参数统称为支持向量机的结构参数。目 前，对支持向量机结构参数的优化方法主要有：试凑法、交叉验证法、网格搜索法、遗传算 法、粒子群算法、人工免疫算法、蚁群算法等。其中，粒子群优化算法(ParticleSwarm Optimization，PSO)需要调整的参数较少，且具有结构简单，容易实现，收敛速度快等优势， 但是由于粒子多样性的快速消失，算法极易陷入局部极值，无法寻找出全局最优解。将混沌 思想引入PSO，不仅改善了粒子群优化算法跳出局部极值点的能力，也在一定程度上提高了 算法的收敛速度和精度，然而由于混沌PSO无法克服单独个体表现出的缺陷，从而容易造成 粒子信息交流的误判。模拟退火算法(SimulatedAnnealing，SA)利用概率突跳特性随机寻 找优化函数的全局最优解，也能够很好的避免在搜索过程中陷入局部的最优解，具有良好 的全局优化性能，但是因为初始值决定了该算法找出问题的全局最优解的性能，所以该算 法也存在一定的局限性。

发明内容

发明目的：针对上述现有技术，提出一种基于多群体协同混沌模拟退火粒子群优 化算法-支持向量机(MCCSAPSO-SVM)的控制系统执行器故障诊断方法，能够通过利用基于 多群体协同混沌模拟退火粒子群优化算法(MCCSAPSO)优化结构参数和采用组合核函数的 偏二叉树支持向量机的输出，实时有效地判断出控制系统执行器是否发生故障，并且能够 在执行器发生故障时较为准确地确定出故障类型。

技术方案：一种基于多群体协同混沌模拟退火粒子群优化算法-支持向量机 (MCCSAPSO-SVM)的控制系统执行器故障诊断方法，通过联合降噪和改进经验模态分解 (EMD)方法，对采集到的执行器输出信号进行降噪处理和特征提取；利用多群体协同混沌模 拟退火粒子群优化算法(MCCSAPSO)优化支持向量机结构参数；采用组合核函数保证了支持 向量机良好的泛化能力与学习能力；利用训练数据构建一种把一个复杂的多分类问题转化 为若干个二分类问题的偏二叉树支持向量机，降低了计算量，提高了诊断的实时性；通过偏 二叉树支持向量机的输出，能够实时有效地判断出控制系统执行器是否发生故障，并且能 够在执行器发生故障时较为准确地确定出故障类型，包括如下具体步骤：

步骤1)将离散含噪信号f(k)(信号长度为N)先进行中值滤波处理得到，滤除 可能的脉冲噪声；

步骤2)对中值滤波处理后的进行改进小波阀值降噪处理：

步骤2.1)将进行第j层(从j＝1开始)小波分解，其中j的确定方式为：对小波 分解得到的高频细节系数d_j(k)进行式(1)自相关系数λ计算，若满足自由度l的x²分布，则对 d_j继续进行j+1层分解，直到不满足x²分布为止；式中，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·{\nu }^{-{\left(\mu \right|W|-\delta )}^2},& \left|W\right|\le \delta \end{array}\right)---\left(2\right)$

步骤2.3)对处理后的各小波系数进行小波逆变换，得到降噪后的真实信号的估计

步骤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,...,\eta ),$则提取出的特征向量为$[E_1/\underset{i=1}{\overset{\eta }{\Sigma }}{E}_i,E_2/\underset{i=1}{\overset{\eta }{\Sigma }}{E}_i,...,E_{\eta }/\underset{i=1}{\overset{\eta }{\Sigma }}{E}_i];$

步骤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的最值，通过调节引入 的收敛因子δ，可以加快粒子群的收敛速度；

步骤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错误分类 数；

Ψ＝100r0.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)/(E_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 \left({\mathrm{wv}}_{Mi}^t+c_1{r}_1\left(p_{Mibest}-x_{Mi}^t\right)+{\mathrm{\phi c}}_2{r}_2\left(g_{Mbest}-x_{Mi}^t\right)+\left(1-\phi \right)c_3{r}_3\left(g_{Sbest}-x_{Si}^t\right)\right)\\ 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及适应值变化量ΔΨ_i＝ Ψ′_i-Ψ_i；若ΔΨ_i＜0，或ΔΨ_i＞0时exp(-ΔΨ/T)在区间[0，1]上，则进行降温操作T←K_TT， 否则温度不变；

步骤5.7)终止条件：当满足温度降至T₀，或达到最大迭代次数t_max时停止迭代，否 则返回步骤5.3；

步骤6)构造组合核函数支持向量机：

步骤6.1)构造式(12)形式的组合核函数，式中，d为多项式核函数的阶数，σ为高斯 核函数的核半径，ρ为权重；

K＝ρK_poly+(1-ρ)K_RBF＝ρ(xx_i+1)^d+(1-ρ)exp(-‖x-x_i‖²/σ²)(12)

步骤6.2)构造组合核函数支持向量机的决策函数可表示为式(13)，式中，{x_i，y_i} 为训练样本，l是训练样本个数，K(x_i，x_j)是给定的核函数，α_i和b为根据训练样本得出的参 数；

$f\left(x\right)=\mathrm{sgn}\left(\underset{i=1}{\overset{l}{\Sigma }}{\alpha }_i{y}_i\left({\mathrm{\rho K}}_{poly}\left(x_i,x_j\right)+\left(1-\rho \right)K_{RBF}\left(x_i,x_j\right)\right)+b\right)---\left(13\right)$

步骤6.3)将组合核函数支持向量机的优化模型(14)转化为式(15)，式中，C为惩罚 因子，ξ_i为松弛变量；

$\begin{array}{cc}\underset{\rho ,\left(1-\rho \right),f_{poly},f_{RBF},b,{\xi }_i}{\min }& \left(\frac{1}{\rho }\left|\right|f_{poly}|{|}^2\right)+\left(\frac{1}{\left(1-\rho \right)}\left|\right|f_{RBF}|{|}^2\right)+C\underset{i}{\overset{l}{\Sigma }}{\xi }_i\\ s.t.& y_i\left(f_{poly}+f_{RBF}+b\right)\ge 1-{\xi }_i\\ & {\xi }_i\ge 0,0\le \rho \le 1\end{array}---\left(14\right)$

$\begin{array}{cc}\underset{\rho ,\left(1-\rho \right),f_{poly},f_{RBF},b,{\xi }_i}{\min }& \frac{1}{2}\left(\frac{1}{\rho }\left|\right|f_{poly}|{|}^2\right)+\left(\frac{1}{\left(1-\rho \right)}\left|\right|f_{RBF}|{|}^2\right)+C\underset{i}{\overset{l}{\Sigma }}{\xi }_i\\ s.t.& y_i\left(f_{poly}+f_{RBF}+b\right)\ge 1-{\xi }_i\\ & {\xi }_i\ge 0,0\le \rho \le 1\end{array}---\left(15\right)$

步骤6.4)利用Largrange乘数法和KKT条件，计算其对偶问题(16)；

$\begin{array}{cc}\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}\left(x_i,x_j\right)+\left(1-\rho \right)K_{RBF}\left(x_i,x_j\right)\right)\\ s.t.& \underset{i=1}{\overset{l}{\Sigma }}{\alpha }_i{y}_i=0,0\le {\alpha }_i\le C,i=1,\mathrm{...},l\\ & 0\le \rho \le 1\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\left({\mathrm{\rho K}}_{poly}\left(x_i,x_j\right)+\left(1-\rho \right)K_{RBF}\left(x_i,x_j\right)\right)+b^{*}\right)---\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的输 出值实时有效地判断出控制系统执行器是否发生故障，以及准确判定执行器发生的具体故 障类型。

有益效果：一种基于多群体协同混沌模拟退火粒子群优化算法-支持向量机 (MCCSAPSO-SVM)的控制系统执行器故障诊断方法，通过联合降噪和改进经验模态分解 (EMD)方法，对采集到的执行器输出信号进行降噪处理和特征提取；利用多群体协同混沌模 拟退火粒子群优化算法(MCCSAPSO)优化支持向量机结构参数；采用组合核函数保证了支持 向量机良好的泛化能力与学习能力；利用训练数据构建一种把一个复杂的多分类问题转化 为若干个二分类问题的偏二叉树支持向量机，降低了计算量，提高了诊断的实时性；通过偏 二叉树支持向量机的输出，能够实时有效地判断出控制系统执行器是否发生故障，并且能 够在执行器发生故障时较为准确地确定出故障类型。具有如下具体优点：

①该支持向量机采用高斯核函数与多项式核函数的线性组合，既具有全局核函数 良好的泛化能力，也具有局部核函数良好的学习能力；

②支持向量机的结构参数优化采用多群体协同方法，有效克服了单一粒子信息交 流可能造成的误判问题；

③在每个粒子的速度和位置更新过程中加入模拟退火思想，粒子群迭代后的适应 值既接受优化解，也一定概率的接受恶化解，这样可以使PSO从局部极值区域跳出，并调整 模拟退火算法的温度；

④通过偏二叉树支持向量机的输出，不仅能够实时有效地判断出控制系统执行器 是否发生故障，而且能够在执行器发生故障时较为准确地确定出故障类型；

⑤本发明提高了故障诊断的针对性和准确性，实现了控制系统执行器故障的自动 化方式诊断，同时对噪声干扰具有较强的鲁棒性。

本发明所提方法作为一种控制系统执行器故障诊断的改进方法，具有一定的应用 意义，易于实现，实时性好，准确性高，能够有效提高控制系统安全性且可操作性强，节省时 间，效率更高，可广泛应用于各种控制系统的执行器故障诊断中。

附图说明

图1是本发明方法的流程图；

图2是控制系统执行器故障诊断框图；

图3是特征提取与故障诊断模块图；

图4是偏二叉树支持向量机故障诊断工作原理图；

图5是Quanser公司研制的用以研究四旋翼直升机控制的实验装置Qball-X4四旋 翼直升机；

图6是电机正常及在各种故障状态下的样本数据图；

图7是卡死故障状态下的电机转速图；

图8是SVM故障诊断结果图；

图9是恒偏差故障状态下的电机转速图；

图10是SVM故障诊断结果图。

具体实施方式

下面结合附图对本发明做更进一步的解释。

如图1所示，一种基于多群体协同混沌模拟退火粒子群优化算法-支持向量机 (MCCSAPSO-SVM)的控制系统执行器故障诊断方法，通过联合降噪和改进经验模态分解 (EMD)方法，对采集到的执行器输出信号进行降噪处理和特征提取；利用多群体协同混沌模 拟退火粒子群优化算法(MCCSAPSO)优化支持向量机结构参数；采用组合核函数保证了支持 向量机良好的泛化能力与学习能力；利用训练数据构建一种把一个复杂的多分类问题转化 为若干个二分类问题的偏二叉树支持向量机，降低了计算量，提高了诊断的实时性；通过偏 二叉树支持向量机的输出，能够实时有效地判断出控制系统执行器是否发生故障，并且能 够在执行器发生故障时较为准确地确定出故障类型，包括如下具体步骤：

步骤1)将离散含噪信号f(k)(信号长度为N)先进行中值滤波处理得到滤除 可能的脉冲噪声；

步骤2)对中值滤波处理后的进行改进小波阀值降噪处理：

步骤2.1)将进行第j层(从j＝1开始)小波分解，其中j的确定方式为：对小波 分解得到的高频细节系数d_j(k)进行式(1)自相关系数λ计算，若满足自由度l的x²分布，则对 d_j继续进行j+1层分解，直到不满足x²分布为止；式中，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·{\nu }^{-{\left(\mu \right|W|-\delta )}^2},& \left|W\right|\le \delta \end{array}\right)---\left(2\right)$

步骤2.3)对处理后的各小波系数进行小波逆变换，得到降噪后的真实信号的估计

步骤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,...,\eta ),$则提取出的特征向量为$[E_1/\underset{i=1}{\overset{\eta }{\Sigma }}{E}_i,E_2/\underset{i=1}{\overset{\eta }{\Sigma }}{E}_i,...,E_{\eta }/\underset{i=1}{\overset{\eta }{\Sigma }}{E}_i];$

步骤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的最值，通过调节引入 的收敛因子δ，可以加快粒子群的收敛速度；

步骤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错误分类 数；

Ψ＝100r0.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 \left({\mathrm{wv}}_{Mi}^t+c_1{r}_1\left(p_{Mibest}-x_{Mi}^t\right)+{\mathrm{\phi c}}_2{r}_2\left(g_{Mbest}-x_{Mi}^t\right)+\left(1-\phi \right)c_3{r}_3\left(g_{Sbest}-x_{Si}^t\right)\right)\\ 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及适应值变化量ΔΨ_i＝ Ψ′_i-Ψ_i；若ΔΨ_i＜0，或ΔΨ_i＞0时exp(-ΔΨ/T)在区间[0，1]上，则进行降温操作T←K_TT， 否则温度不变；

步骤5.7)终止条件：当满足温度降至T₀，或达到最大迭代次数t_max时停止迭代，否 则返回步骤5.3；

步骤6)构造组合核函数支持向量机：

步骤6.1)构造式(12)形式的组合核函数，式中，d为多项式核函数的阶数，σ为高斯 核函数的核半径，ρ为权重；

K＝ρK_poly+(1-ρ)K_RBF＝ρ(xx_i+1)^d+(1-ρ)exp(-‖x-x_i‖²/σ²)(12)

步骤6.2)构造组合核函数支持向量机的决策函数可表示为式(13)，式中，{x_i，y_i} 为训练样本，l是训练样本个数，K(x_i，x_j)是给定的核函数，α_i和b为根据训练样本得出的参 数；

$f\left(x\right)=\mathrm{sgn}\left(\underset{i=1}{\overset{l}{\Sigma }}{\alpha }_i{y}_i\left({\mathrm{\rho K}}_{poly}\left(x_i,x_j\right)+\left(1-\rho \right)K_{RBF}\left(x_i,x_j\right)\right)+b\right)---\left(13\right)$

步骤6.3)将组合核函数支持向量机的优化模型(14)转化为式(15)，式中，C为惩罚 因子，ξ_i为松弛变量；

$\begin{array}{cc}\underset{\rho ,\left(1-\rho \right),f_{poly},f_{RBF},b,{\xi }_i}{\min }& \left(\frac{1}{\rho }\left|\right|f_{poly}|{|}^2\right)+\left(\frac{1}{\left(1-\rho \right)}\left|\right|f_{RBF}|{|}^2\right)+C\underset{i}{\overset{l}{\Sigma }}{\xi }_i\\ s.t.& y_i\left(f_{poly}+f_{RBF}+b\right)\ge 1-{\xi }_i\\ & {\xi }_i\ge 0,0\le \rho \le 1\end{array}---\left(14\right)$

$\begin{array}{cc}\underset{\rho ,\left(1-\rho \right),f_{poly},f_{RBF},b,{\xi }_i}{\min }& \frac{1}{2}\left(\frac{1}{\rho }\left|\right|f_{poly}|{|}^2\right)+\left(\frac{1}{\left(1-\rho \right)}\left|\right|f_{RBF}|{|}^2\right)+C\underset{i}{\overset{l}{\Sigma }}{\xi }_i\\ s.t.& y_i\left(f_{poly}+f_{RBF}+b\right)\ge 1-{\xi }_i\\ & {\xi }_i\ge 0,0\le \rho \le 1\end{array}---\left(15\right)$

步骤6.4)利用Largrange乘数法和KKT条件，计算其对偶问题(16)；

$\begin{array}{cc}\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}\left(x_i,x_j\right)+\left(1-\rho \right)K_{RBF}\left(x_i,x_j\right)\right)\\ s.t.& \underset{i=1}{\overset{l}{\Sigma }}{\alpha }_i{y}_i=0,0\le {\alpha }_i\le C,i=1,\mathrm{...},l\\ & 0\le \rho \le 1\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\left({\mathrm{\rho K}}_{poly}\left(x_i,x_j\right)+\left(1-\rho \right)K_{RBF}\left(x_i,x_j\right)\right)+b^{*}\right)---\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的输 出值实时有效地判断出控制系统执行器是否发生故障，以及准确判定执行器发生的具体故 障类型。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人 员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应 视为本发明的保护范围。

如图2、3所示，一种基于多群体协同混沌模拟退火粒子群优化算法-支持向量机 (MCCSAPSO-SVM)的控制系统执行器故障诊断方法的具体实现方式：就特征提取与故障诊断 模块而言，包括作为上位机的计算机1，计算机1中装有多通道数据采集卡2、降噪模块3、特 征量提取模块4和支持向量机故障诊断模块5。由于需要对控制系统的执行器部分进行故障 诊断，因此就将特征提取与支持向量机故障诊断装置并入控制系统的执行器部分旁。其中， 偏二叉树支持向量机故障诊断装置的工作原理如图4。

下面以实际案例仿真说明实施方案的有效性。

采用由加拿大Quanser公司研制的Qball-X4四旋翼直升机飞行控制系统执行器作 为应用研究对象。如图5所示，上位机中包含故障诊断模块，与四旋翼直升机通过WiFi进行 通信，四旋翼中包括嵌入式计算机和数据采集卡。该飞行器存在六维度变量即(X，Y，Z，ψ，θ， φ)，其中X，Y，Z为位置变量，ψ为偏航角，θ为俯仰角，φ为滚转角。本文以飞行器向X轴方向 平飞时的电机转速作为测试对象，分别在电机正常(normalstate)、卡死(lockinplace， LIP)、饱和(hardoverfault，HOF)、恒偏差(constantdeviation)、恒增益(constant gain)故障等状态下，采集电机转速信号进行故障诊断。采样频率为1kHz，采样点为1000。图 6为电机正常及在各种故障状态下的样本数据。

首先，采用联合降噪和改进EMD方法，小波函数选择db5小波，分解层数确定为5层， 计算相关系数ρ如表1并与预设阀值ζ＝0.1max(ρ)比较，剔除伪分量并提取出正常及各故障 状态下样本的特征向量。接着，通过4个SVM分类器进行具体故障诊断。第一个分类器SVM₁是 对正常样本与其余所有故障样本进行分类，将正常样本作为一类，分类器输出为+1，其余所 有故障样本的分类器输出为-1；第二个分类器SVM₂是对卡死故障样本与剩余的其他故障样 本进行分类，将卡死故障样本作为一类，分类器输出为+1，剩余的其他故障样本的分类器输 出为-1；第3个分类器SVM₃是对饱和故障样本与剩余的其他两种故障样本进行分类，饱和故 障样本的分类器输出为+1，剩余的其他两种故障样本的分类器输出为-1；第4个分类器SVM₄是对恒增益故障样本与恒偏差故障样本进行分类，恒偏差故障样本的分类器输出为+1，恒 增益故障样本的分类器输出为-1。

其次，通过MCCSAPSO算法优化SVM结构参数。其中，粒子群主群数为1，从群数为3， 粒子群规模为L＝30，最大迭代次数t_max＝500，学习因子c₁＝2.05，c₂＝2.05，c₂＝0.8，权重 系数w_min＝0.2，w_max＝0.9，初始温度T＝100，结束温度T₀＝0.01，退火速度K_T＝0.9。

最后，利用训练好的故障诊断支持向量机，对故障进行辨识。由于篇幅所限，仅以 卡死和恒偏差故障辨识为例。

卡死故障发生在第550个采样点并一直持续到采样结束，诊断结果如图7、8所示。 从图8中可以看出，当故障数据通过SVM₁-SVM₄时，SVM₂能够立即给出准确的诊断结果，且其 余SVM不存在误判的情况。

恒偏差故障发生在第550个采样点，诊断结果如图9、10所示。从图10中可以看出， 当发生此故障时，SVM₄能够迅速作出诊断同时持续10个采样点，综合4个SVM得出与实际情 况完全一致的诊断结果。

从以上的故障诊断仿真结果可以看出，支持向量机分类方法能够对传感器故障进 行在线诊断，并具有能够识别故障类型和快速、准确性较高的特点。

表1相关系数ρ