CN112861289B - 基于IMM-KEnPF的风机桨距系统故障诊断方法 - Google Patents
基于IMM-KEnPF的风机桨距系统故障诊断方法 Download PDFInfo
- Publication number
- CN112861289B CN112861289B CN202110254185.9A CN202110254185A CN112861289B CN 112861289 B CN112861289 B CN 112861289B CN 202110254185 A CN202110254185 A CN 202110254185A CN 112861289 B CN112861289 B CN 112861289B
- Authority
- CN
- China
- Prior art keywords
- model
- particle
- fault
- probability
- particles
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 69
- 238000003745 diagnosis Methods 0.000 title claims abstract description 60
- 239000002245 particle Substances 0.000 claims abstract description 150
- 238000009826 distribution Methods 0.000 claims abstract description 41
- 238000001914 filtration Methods 0.000 claims abstract description 37
- 238000005070 sampling Methods 0.000 claims abstract description 37
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 23
- 230000002452 interceptive effect Effects 0.000 claims abstract description 22
- 238000004364 calculation method Methods 0.000 claims abstract description 12
- 230000006870 function Effects 0.000 claims description 28
- 239000011159 matrix material Substances 0.000 claims description 26
- 230000003993 interaction Effects 0.000 claims description 24
- 238000012952 Resampling Methods 0.000 claims description 14
- 230000008569 process Effects 0.000 claims description 13
- 238000012546 transfer Methods 0.000 claims description 10
- 230000004927 fusion Effects 0.000 claims description 9
- 230000007704 transition Effects 0.000 claims description 9
- 238000005457 optimization Methods 0.000 claims description 7
- 238000009825 accumulation Methods 0.000 claims description 4
- 238000010606 normalization Methods 0.000 claims description 4
- 230000003044 adaptive effect Effects 0.000 claims description 3
- 230000009191 jumping Effects 0.000 claims description 3
- 238000005315 distribution function Methods 0.000 claims description 2
- 101150049349 setA gene Proteins 0.000 claims description 2
- 239000000126 substance Substances 0.000 claims description 2
- 230000009466 transformation Effects 0.000 claims description 2
- 230000002776 aggregation Effects 0.000 claims 2
- 238000004220 aggregation Methods 0.000 claims 2
- 230000015556 catabolic process Effects 0.000 abstract description 6
- 238000006731 degradation reaction Methods 0.000 abstract description 6
- 238000004088 simulation Methods 0.000 description 6
- 238000012986 modification Methods 0.000 description 5
- 230000004048 modification Effects 0.000 description 5
- 238000013461 design Methods 0.000 description 4
- 230000006872 improvement Effects 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 238000013016 damping Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000013178 mathematical model Methods 0.000 description 2
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000002405 diagnostic procedure Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 235000003642 hunger Nutrition 0.000 description 1
- 230000000116 mitigating effect Effects 0.000 description 1
- 230000006855 networking Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000037351 starvation Effects 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
- 238000013024 troubleshooting Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/17—Mechanical parametric or variational design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/25—Design optimisation, verification or simulation using particle-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/06—Wind turbines or wind farms
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- General Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Computer Hardware Design (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Testing And Monitoring For Control Systems (AREA)
Abstract
本发明公开了基于IMM‑KEnPF的风机桨距系统故障诊断方法，涉及故障诊断技术领域。针对变桨距系统的非线性非高斯的特点引入粒子滤波算法，从准确性和实时性两个方面对其进行了改进。首先引入集合卡尔曼滤波(EnKF)优化PF的建议分布，缓解PF的粒子退化问题，提高PF的估计精度；其次引入KLD采样，根据后验估计分布与真实分布的KLD来自适应的调整所需要采样的粒子数目，提高PF的计算效率和实时性，最后将交互式多模型方法和粒子滤波算法结合，提出了基于IMM‑KEnPF的风机桨距系统故障诊断方法。
Description
技术领域
本发明涉及故障诊断技术领域，特别是涉及基于IMM-KEnPF的风机桨距系统故障诊断方法。
背景技术
风电机组中的变桨距系统是故障常发部件之一，故障会影响风机功率输出的稳定性，并会导致转子、转轴和轴承上的负载不平衡，给风机的安全运行带来隐患。因此，对风机变桨距系统的故障诊断方法进行研究，对于保障风电机组的安全稳定运行具有重要作用。
相对于风电机组的其他子系统，变桨距系统故障诊断方法的研究成果较少。ZhangYoumin最先提出了基于IMM的故障诊断的基本框架，该方法对于系统参数或者结构会发生变化的系统也可以较好地完成自适应估计。胡艳艳等提出了网络化异步IMM融合滤波算法，将未知的故障幅值作为增广的系统状态变量，根据模型概率检测和定位故障，同时得到故障幅值和系统状态的联合估计。王剑等将无迹卡尔曼滤波与交互式多模型相结合，提出了IMM-UKF故障诊断算法。GADSDEN S A等使用时变平滑边界层(VBL)来改进平滑变结构滤波器(SVSF)，并与交互式多模型结合设计出基于IMM-SVVF-VBL的故障诊断算法。但是UKF和SVSF处理非线性的能力有限，诊断性能有待进一步提升。Zhengjiang提出将PF与交互式多模型方法结合，设计出一种适用于强非线性系统的故障诊断方法。但PF中存在的粒子退化和样本匮乏导致的估计精度下降以及计算耗时问题。同时，交互式多模型需要根据所建立的模型集，集成多个滤波器，模型集中的子模型越多所需要的滤波器就越多，计算负荷会急剧增加。
针对以上问题，目前的研究中主要使用群智能算法对粒子后验分布进行优化，但是群智能优化容易陷入局部最优，且寻优迭代过程会耗费许多计算时间，而使用次优滤波器直接优化建议分布改善重要性采样过程则更直观和高效。其中常见的次优滤波器是无迹卡尔曼滤波(UKF)。Merwe等将无迹卡尔曼滤波(Unscented Kalman Filter,UKF)引入粒子滤波，根据最新的观测值通过无迹卡尔曼滤波来预测状态值的变化，以此来优化粒子滤波的建议分布，缓解粒子退化现象。但是UKF中的Sigma点的计算时间会随着维度和粒子增加而急剧增加。Delft G V等设计出一种选取随机点，用来预测状态后验分布作为粒子滤波的建议分布，称作集合粒子滤波(Ensemble Particle Filter,EnPF)，有着更好的处理系统非线性的能力和计算效率。当粒子数量增加时，使用次优滤波器优化粒子滤波的方法同时也增加了算法的计算复杂度，会严重影响算法的实时性，如何在提高粒子滤波算法估计精度的同时从算法层面提高算法的实时性仍是一个研究难题。
发明内容
本发明提供的基于IMM-KEnPF的风机桨距系统故障诊断方法，可以解决现有技术中存在的问题。
本发明提供了基于IMM-KEnPF的风机桨距系统故障诊断方法，包括以下步骤：
步骤1，输入交互；风机变桨距系统的多种故障模型构成故障模型集，故障模型集中的故障模型之间随机跳变切换，前一时刻的故障模型到下一时刻的故障模型之间的切换服从故障模型转移概率；基于交互式多模型IMM，利用上一时刻的各故障模型概率和故障模型转移概率计算故障模型之间的交互概率，再使用交互概率融合各粒子滤波器的上一时刻的滤波信息，得到各故障模型在上一时刻初始的估计交互值和对应的误差协方差矩阵；
步骤2，并行滤波；利用步骤1得到的各故障模型在上一时刻初始的估计交互值和对应的误差协方差矩阵，为各个粒子滤波器生成粒子集，按照集合卡尔曼滤波EnKF优化的建议分布以及变桨距系统的观测值，对粒子集中的粒子进行重要性采样，得到粒子的预测集合；在重要性采样过程中，利用库尔贝克-莱布勒散度KLD采样算法自适应调整采样粒子的数量；
得到预测集合后，根据变桨距系统的观测值计算粒子权值，对粒子残值利用残差重采样方法进行重采样，得到新的粒子集合，将新的粒子集合中粒子的粒子均值、各个粒子滤波器状态估计的误差协方差矩阵、估计残差和残差协方差矩阵作为输出；
步骤3，模型概率更新；根据步骤2输出的估计残差和残差协方差矩阵建立各故障模型滤波输出的似然函数，使用该似然函数计算各个故障模型的模型概率；
步骤4，融合估计；得到各个故障模型的模型概率后，将其与步骤2输出的粒子均值和各个粒子滤波器状态估计的误差协方差矩阵进行融合，得到估计结果；
步骤5，故障诊断；基于步骤3得到的各个故障模型的模型概率，利用决策函数实现变桨距系统的故障诊断。
本发明中的基于IMM-KEnPF的风机桨距系统故障诊断方法，针对变桨距系统的非线性非高斯的特点引入粒子滤波算法，然后从准确性和实时性两个方面对其进行了改进，最后将交互式多模型方法和粒子滤波算法结合。首先引入集合卡尔曼滤波(EnKF)优化PF的建议分布，缓解PF的粒子退化问题，提高PF的估计精度；其次引入KLD采样，根据后验估计分布与真实分布的KLD来自适应的调整所需要采样的粒子数目，提高PF的计算效率和实时性。最后将改进的PF与IMM相结合，在风电机组仿真平台上实现了变桨距系统多传感器故障的诊断。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明中诊断方法的流程图；
图2为仿真实验框架；
图3为模型概率的变化曲线；
图4为故障诊断结果示意图；
图5为状态估计结果示意图；
图6为本发明方法和对比方法的计算时间对比。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
在介绍本发明的方法前，首先对风机变桨距系统的故障进行介绍。
将风机变桨距系统考虑成一个活塞伺服系统，进而将其建模成一个二阶传递函数模型：
式中，ζ为变桨距系统的阻尼系数，ωn为变桨距系统的自然频率。阻尼系数和自然频率的变化可以反应变桨距系统的动态特性，通过实验辨识参数得到变桨距系统正常工作时ζ＝0.6，ωn＝11.11rad/s。
根据上述的基准模型对不同的故障建立相应的数学模型。将变桨距系统的数学模型从传递函数的形式转化成状态空间的形式：
式中，f(·)为变桨距系统的状态转移函数，h(·)为观测函数，wk-1为k-1时刻的过程噪声，vk为k时刻的观测噪声，fs,k为k时刻出现的传感器故障。为了模拟传感器可能会出现的多种故障，并能在特定的时刻注入故障，本发明为fs,k做了如表1所示的定义。
表1传感器故障模型
参照图1，本发明提供了基于IMM-KEnPF的风机桨距系统故障诊断方法，该方法包括以下五个步骤：
步骤1，输入交互。风机变桨距系统的多种故障模型构成故障模型集，故障模型集中的故障模型之间是随机跳变切换的，前一时刻的故障模型到下一时刻的故障模型之间的切换服从故障模型转移概率。基于交互式多模型IMM，利用上一时刻的各故障模型概率和故障模型转移概率计算故障模型之间的交互概率，再使用交互概率融合各粒子滤波器的上一时刻的滤波信息，得到各故障模型在上一时刻初始的估计交互值和对应的误差协方差矩阵，实现粒子滤波器的初始化。
步骤2，并行滤波。利用步骤1得到的各故障模型在上一时刻初始的估计交互值和对应的误差协方差矩阵，为各个粒子滤波器生成粒子集，按照集合卡尔曼滤波EnKF优化的建议分布以及变桨距系统的观测值，对粒子集中的粒子进行重要性采样，得到粒子的预测集合。在重要性采样过程中，利用库尔贝克-莱布勒散度KLD采样算法自适应调整采样粒子的数量。
得到预测集合后，根据变桨距系统的观测值计算粒子权值，对粒子残值利用残差重采样方法进行重采样，得到新的粒子集合，将新的粒子集合中粒子的均值、各个粒子滤波器状态估计的误差协方差矩阵、估计残差和残差协方差矩阵作为输出。
步骤3，模型概率更新。根据步骤2输出的估计残差和残差协方差矩阵建立各故障模型滤波输出的似然函数，使用该似然函数计算各个故障模型的模型概率。
步骤4，融合估计。得到各个故障模型的模型概率后，将其与步骤2输出的粒子均值和各个粒子滤波器状态估计的误差协方差矩阵进行融合，得到融合估计结果。
步骤5，故障诊断。基于步骤3得到的各个故障模型的模型概率，利用决策函数实现变桨距系统的工作模式诊断。
下面对每个步骤进行详细说明。
步骤1，由多个故障模型建立的故障模型集为M＝{M1,M2,…,Mr}，并用如式(3)所示的混杂系统状态空间模型表示：
式中：Mk-1和Mk分别是与变桨距系统工作模式相匹配的前一时刻和当前时刻的故障模型，交互式多模型IMM中假设混杂系统的离散模式之间的切换是随机跳变的，为一阶马尔可夫过程，前一时刻的故障模型i到下一时刻故障模型j之间的切换服从：
P{Mk＝Mj|Mk-1＝Mi}＝πij,i,j＝1,2,…r (4)
式中，πij为模型转移概率，且
是决定每一时刻粒子滤波器之间进行状态估计交互的参数，一般根据经验和先验知识预设。
将交互式多模型IMM与粒子滤波KEnPF进行结合时，交互环节的实施有模型层面的、粒子层面的以及前两者相结合的三种方式。现有文献对以上三种交互的实施方式做了对比实验，结果表明第一种方式的目标跟踪的性能最好，因此，本发明在故障诊断设计上也采取该方式。
在步骤1迭代周期开始时，利用上一时刻的各故障模型概率和故障模型转移概率计算故障模型之间的交互概率，再使用交互概率融合各粒子滤波器的上一时刻的滤波信息进行融合，实现滤波器的初始化，交互的过程如式(5)-(9)所示：
式中：i,j∈M，
为k-1时刻第i个模型粒子滤波器的估计值，
为k-1时刻第i个模型粒子滤波器的估计误差协方差矩阵，
为k-1时刻第j个模型在迭代周期初始时的估计交互值，
为与其对应的误差协方差矩阵，
为k-1时刻模型i的模型概率，
为k-1时刻模型i与模型j的交互概率，πij,k-1为k-1时刻模型i与模型j之间的模型转移概率矩阵。
步骤2，首先，利用前一步的交互结果
和
为各个粒子滤波器生成粒子集
N为粒子数。
其中，rand表示随机函数。
随后，进入改进粒子滤波KEnPF的流程，详细流程如下：
1、初始化。在k＝0时刻，从先验概率分布中随机采样N个权值大小相同的粒子
粒子权值
权值累加和初值S＝0。
2、定义背景集合。为第i个粒子随机采样n个集合点，构成背景集合
3、更新背景集合。
随机采样n个初始样本集合
作为上一时刻的背景集合
在预测阶段，利用状态方程进行背景集合的更新得到新一时刻的背景集合
如下式所示：
根据式(12)计算背景集合的均值，并用式(13)和式(14)代替方差计算。
通过下式计算卡尔曼增益：
式中，Rk为k时刻系统的观测误差协方差矩阵。然后，用式(16)根据观测信息zk和计算所得的卡尔曼增益更新背景集合得到新的集合
称作分析集合：
最后根据下式计算分析集合的均值和方差：
使用式(11)-(16)更新背景集合得到分析集合
并用式(17)-(18)计算分析集合的均值
和方差
4、重要性采样和权值更新。从EnKF优化的建议分布，即上述得到的各个粒子滤波器生成粒子集中采集粒子，得到粒子的预测集合：
将标准粒子滤波器中的状态转移函数作为建议分布替换，使得每个粒子
的重要性采样过程，都以EnKF更新得到的分析集合
作为建议分布：
由于充分利用最新时刻的观测信息，获得的建议分布将更合理，可改善粒子滤波的粒子退化现象和样本匮乏，更易于应对系统状态的突变。
对预测集合中的第i个粒子进行重要性采样
然后更新粒子的权值并累加：
S＝S+ωi,k (22)
5、更新自适应粒子数目。判断xi,k是否落入空子集Xk，若是则k＝k+1，当i＞nmin，nmin表示最小采样数，则用式(27)更新nkld，当i＜nmin，i＝i+1。
6、当i＜nkld，则返回2，否则跳出2到5的循环。
7、权值归一化和重采样。使用权值累加和进行权值归一化：
得到粒子集
对其进行重采样并将权值赋为1/N，得到粒子集
8、加权输出估计结果：
KLD采样的方法为：
任意两个分布p与q之间的KLD由式(25)定义：
利用KLD调整粒子采样规模，就是要使粒子的采样数目满足真实后验分布与估计后验分布之间的KLD小于ε的概率为1-δ。假设有n个粒子从k个离散的空间中随机地取出，用X＝(X1,…,Xk)代表从每个空间中取出粒子的数目，X服从多项式分布，也就是X～M(n,p)，其中p＝(p1,…,pk)是从对应空间中取出粒子的概率。使用采样的n个粒子能得到概率p的最大似然估计是
可以看出所采样的粒子数越大，其越接近于真实的后验分布，因此可以基于此进行推得所需要的最优粒子数：
其中，
表示自由度为k-1的卡方分布，1-δ为卡方分布的分位数。为了根据式(26)确定nkld，根据Wilson-Hilferty变换可以获得卡方分布分位点的估计，进而最终可以得到：
式中，z1-δ是标准正态分布的上1-δ分位值，其中的参数设定：δ＝0.01，ε＝0.15。
根据观测值计算残差，用残差计算粒子权值，并归一化。
式中：zk为k时刻的观测值，利用计算得到的粒子权值根据残差重采样方法进行重采样，得到新的粒子集合
将粒子的均值作为输出：
最后，分别计算各个粒子滤波器状态估计的误差协方差矩阵
估计残差
和残差协方差矩阵
步骤3，使用各故障模型滤波输出的似然函数来计算故障模型概率，如式(35)所示。
式中，似然函数根据各滤波器估计结果的残差
和残差方差
进行计算，由于观测噪声通常服从高斯分布，故似然函数选为标准正态分布函数：
故障模型概率在交互式多模型故障诊断算法中起到关键性的作用，除了利用其对滤波器的滤波结果进行融合输出，也需要被用来进行故障模型的匹配进行诊断。在系统处于工作正常的情况下，与正常模式对应故障模型的故障模型概率应该接近于1，当工作模式发生变化，或发生故障时，对应的匹配故障模型的故障模型概率也就进阶地上升并接近于1，而其他故障模型的故障模型概率下降到接近于0。
步骤4，在获得故障模型概率之后，使用故障模型概率对各滤波器的滤波结果进行加权融合，得到IMM-KEnPF的估计结果，如(37)-(39)所示。
步骤5，对于故障模型概率所蕴含的模式信息，可以利用式(40)-(41)的决策函数实现工作模式的识别和诊断。当某一个故障模型与当前系统匹配，而其他故障模型都失配时，说明系统正处于对应的故障模式当中。
式中，μ′T为预先设定的诊断阈值，阈值小则诊断延迟小，但可能会因为模型竞争导致误诊，为提高诊断率而将阈值增大则会增加诊断的时间延迟。
仿真实验
为了运用所提出的基于IMM-KEnPF方法解决变桨距系统的故障诊断问题并验证所提出改进的效果，本发明建立了如图2所示的实验框架，先通过Matlab/Simulation搭建的风机仿真平台采集变桨距系统含有故障的运行数据，再将数据传给交互式多模型进行故障诊断和状态估计。
首先，根据建立的故障模型集，对应设计4个KEnPF滤波器，通过各滤波的滤波结果进行故障模型概率的计算，从而进行变桨系统的故障诊断和状态估计，最后将此方法与基于MMAE-PF以及IMM-PF的故障诊断方法进行了对比。实验参数设置如下：仿真时间为T＝6s，采样间隔ΔT＝0.01，粒子滤波器粒子数N＝50，过程噪声和观测噪声分别为w～Γ(0.1,0.1)，v～N(0,0.001)，KEnPF中最小采样粒子数nmin＝20，EnKF集合样本数n＝10。桨距角估计值初始值x0＝10，初始协方差矩阵方差P0＝diag[10 10]，诊断阈值μ'T设置为9.9。初始故障模型概率为μ0＝(0.7,0.1,0.1,0.1)T，初始故障模型转移概率设置为：
实验硬件为Intel(R)Corei5-8250U CPU@1.6GHz，内存8G，操作系统为Windows10，编程软件为Matlab R2017b。
方法改进前后的故障模型概率变化对比如图3所示。可以看到在200ΔT处故障模型切换后的故障模型概率发生了反转，故障模型0的故障模型概率上升到接近1，其他故障模型的概率下降到接近0。而改进后的故障诊断算法相较标准IMM-PF算法，由于KEnPF改进的估计性能，故障模型概率也更接近于1。由于诊断和融合估计都要使用故障模型概率，这样的效果有利于提高算法诊断的准确度和状态估计的精度。从图4改进后的诊断结果可以看出，正确识别故障的数量明显增加，故障模型失配得到了减少；从图5改进后的状态估计结果可以看出桨距角也能更好地逼近真实值。
图4中，从上到下分别是故障的预设模式序列、基于MMAE-PF的故障诊断结果、基于标准的IMM-PF故障诊断结果、基于本发明的IMM-KEnPF的故障诊断结果。可以看出基于MMAE-PF方法的故障诊断结果相较于其他两种方法有较多的误报和漏报，基于IMM-PF的故障诊断结果中，除了在故障模式发生切换时存在一定的时间延迟，整个过程有较少的漏报，诊断结果整体上与预设的模式序列相吻合。
图5是MMAE-PF、IMM-PF和IMM-KEnPF三种算法的自适应估计结果，从中可以看出MMAE-PF的估计偏离实际桨距角的程度较大，而IMM-PF大部分区间内的桨距角估计较为准确，误差较小，只是在500～600ΔT发生传感器卡死时，桨距角估计都发生了振荡。这是由于输入值在变化，而滤波器无法获得相应的输出值的更新。图6显示了IMM-PF与IMM-KEnPF的整个故障诊断方法计算的时间对比，可以看出相较于固定粒子数的IMM-PF方法，IMM-KEnPF方法由于具有采用KLD采样对粒子数进行削减，而在计算时间上有大幅的节省，并且随着粒子数的提升，效果越明显。
从图3～6的实验结果可以看出，本发明提出的基于IMM-KEnPF的变桨距传感器的多故障诊断方法总体上是有效的，而在此基础上，对交互式多模型和滤波器的改进进一步提升了故障诊断和状态估计性能。为综合评价这两种算法的性能，选取了以下三个指标：
(1)正确诊断(CDID)：诊断结果与实际模式匹配的个数；
(2)均方根误差(RMSE)：如式(42)所示的交互式多模型的融合估计的均方根误差；
式中：xl,n为第l次迭代第n个采样点的实际值，
为第l次迭代第n个采样点的IMM估计值，由于卡死阶段无法有效估计状态，本发明选取前500个采样点计算均方根误差。
(3)计算时间(runtime):每个故障诊断算法运行周期的计算时间。本发明选取粒子数为50的情况下统计计算时间。
将改进前和改进后的故障诊断方法各独立运行50次蒙特卡洛实验，计算各次的平均值，得到的各性能指标统计结果如表2所示。
表2 IMM-KEnPF与IMM-PF性能对比
从表中的数据可以看出，改进后的故障诊断方法提高了18.22个正确诊断数，提升了近3％，降低了56.17％的估计误差，并且计算时间减少了72.2％。此统计结果表明本发明所提出的IMM-KEnPF的故障诊断方法从诊断准确性、状态估计精度以及实时性都得到了较好的提升。
尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。
Claims (5)
1.基于IMM-KEnPF的风机桨距系统故障诊断方法，其特征在于，包括以下步骤：
步骤1，输入交互；风机变桨距系统的多种故障模型构成故障模型集，故障模型集中的故障模型之间随机跳变切换，前一时刻的故障模型到下一时刻的故障模型之间的切换服从故障模型转移概率；基于交互式多模型IMM，利用上一时刻的各故障模型概率和故障模型转移概率计算故障模型之间的交互概率，再使用交互概率融合各粒子滤波器的上一时刻的滤波信息，得到各故障模型在上一时刻初始的估计交互值和对应的误差协方差矩阵；
步骤2，并行滤波；利用步骤1得到的各故障模型在上一时刻初始的估计交互值和对应的误差协方差矩阵，为各个粒子滤波器生成粒子集，按照集合卡尔曼滤波EnKF优化的建议分布以及变桨距系统的观测值，对粒子集中的粒子进行重要性采样，得到粒子的预测集合；在重要性采样过程中，利用库尔贝克-莱布勒散度KLD采样算法自适应调整采样粒子的数量；
得到预测集合后，根据变桨距系统的观测值计算粒子权值，对粒子残值利用残差重采样方法进行重采样，得到新的粒子集合，将新的粒子集合中粒子的粒子均值、各个粒子滤波器状态估计的误差协方差矩阵、估计残差和残差协方差矩阵作为输出；
步骤3，模型概率更新；根据步骤2输出的估计残差和残差协方差矩阵建立各故障模型滤波输出的似然函数，使用该似然函数计算各个故障模型的模型概率；
步骤4，融合估计；得到各个故障模型的模型概率后，将其与步骤2输出的粒子均值和各个粒子滤波器状态估计的误差协方差矩阵进行融合，得到融合估计结果；
步骤5，故障诊断；基于步骤3得到的各个故障模型的模型概率，利用决策函数实现变桨距系统的故障诊断；
步骤1具体包括：
由多个故障模型建立的故障模型集为M＝{M1，M2，…，Mr}，用如式(1)所示的混杂系统状态空间模型表示：
式中：f(·)为变桨距系统的状态转移函数，h(·)为观测函数，wk-1为k-1时刻的过程噪声，vk为k时刻的观测噪声，Mk-1和Mk分别是与变桨距系统工作模式相匹配的前一时刻和当前时刻的故障模型，交互式多模型IMM中假设混杂系统的离散模式之间的切换随机跳变，为一阶马尔可夫过程，前一时刻的故障模型i到下一时刻故障模型j之间的切换服从：
P{Mk＝M1|Mk-1＝M1|＝πij，i，j＝1，2，…r (2)
式中，πij为模型转移概率；
在迭代周期开始时，利用上一时刻的各故障模型概率和故障模型转移概率计算故障模型之间的交互概率，再使用交互概率融合各粒子滤波器的上一时刻的滤波信息进行融合，实现滤波器的初始化，交互的过程如式(3)-(7)所示：
式中：i，j∈M，
为k-1时刻第i个模型粒子滤波器的估计值，
为k-1时刻第i个模型粒子滤波器的估计误差协方差矩阵，
为k-1时刻第j个模型在迭代周期初始时的估计交互值，
为与其对应的误差协方差矩阵，
为k-1时刻模型i的模型概率，
为k-1时刻模型i与模型j的交互概率，πij，k-1为k-1时刻模型i与模型j之间的模型转移概率矩阵；
步骤2具体包括：
子步骤1，利用步骤1的交互结果
和
为各个粒子滤波器生成粒子集
N为粒子数；
其中，rand表示随机函数；
子步骤2，进入改进粒子滤波KEnPF的流程：
子步骤21，初始化；在k＝0时刻，从先验概率分布中随机采样N个权值大小相同的粒子
粒子权值
权值累加和初值S＝0；
子步骤22，定义背景集合；为第i个粒子随机采样n个集合点，构成背景集合
子步骤23，更新背景集合；随机采样n个初始样本集合
作为上一时刻的背景集合
在预测阶段，利用状态方程进行背景集合的更新得到新一时刻的背景集合
如下式所示：
根据式(10)计算背景集合的均值，并用式(11)和式(12)代替方差计算：
通过下式计算卡尔曼增益：
Kk＝Pk 1HT(HPk 1HT+Rk)-1 (13)
式中，Rk为k时刻系统的观测误差协方差矩阵；然后，用式(14)根据观测信息zk和计算所得的卡尔曼增益更新背景集合得到新的集合
称作分析集合：
最后，根据下式计算分析集合的均值和方差：
子步骤24，重要性采样和权值更新；从EnKF优化的建议分布，即上述得到的各个粒子滤波器生成粒子集中采集粒子，得到粒子的预测集合：
将标准粒子滤波器中的状态转移函数作为建议分布替换，使得每个粒子
的重要性采样过程，都以EnKF更新得到的分析集合
作为建议分布：
对预测集合中的第i个粒子进行重要性采样
然后更新粒子的权值并累加：
S＝S+ωi，k (20)
子步骤25，更新自适应粒子数目；判断xi，k是否落入空子集Xk，若是则k＝k+1，当i＞nmin，nmin表示最小采样数，则更新KLD采样粒子数nk/d，当i＜nmin，i＝i+1；
子步骤26，当i＜nk/d，则返回子步骤2，否则跳出子步骤2到5的循环；
子步骤27，权值归一化和重采样；使用权值累加和进行权值归一化：
得到粒子集
对其进行重采样并将权值赋为1/N，得到粒子集
子步骤28，加权输出估计结果：
子步骤3，根据观测值计算残差，用残差计算粒子权值，并归一化；
式中：zk为k时刻的观测值，利用计算得到的粒子权值根据残差重采样方法进行重采样，得到新的粒子集合
将粒子的均值作为输出：
子步骤4，分别计算各个粒子滤波器状态估计的误差协方差矩阵
估计残差
和残差协方差矩阵
2.如权利要求1所述的基于IMM-KEnPF的风机桨距系统故障诊断方法，其特征在于，KLD采样的方法为：
任意两个分布p与q之间的KLD由式(30)定义：
利用KLD调整粒子采样规模，就是要使粒子的采样数目满足真实后验分布与估计后验分布之间的KLD小于ε的概率为1-δ；假设有n个粒子从k个离散的空间中随机地取出，用X＝(X1，...，Xk)代表从每个空间中取出粒子的数目，X服从多项式分布，也就是X～M(n，p)，其中p＝(p1，...，pk)是从对应空间中取出粒子的概率；使用采样的m个粒子能得到概率p的最大似然估计是
基于此进行推得所需要的KLD采样粒子数：
其中，
表示自由度为k-1的卡方分布，1-δ为卡方分布的分位数，为了根据式(31)确定nk/d，根据Wilson-Hilferty变换获得卡方分布分位点的估计，进而最终得到：
式中，z1-δ是标准正态分布的上1-δ分位值。
3.如权利要求1所述的基于IMM-KEnPF的风机桨距系统故障诊断方法，其特征在于，步骤3具体包括：
使用各故障模型滤波输出的似然函数来计算故障模型概率，如式(33)所示：
式中，似然函数根据各滤波器估计结果的残差
和残差方差
进行计算，由于观测噪声通常服从高斯分布，故似然函数选为标准正态分布函数：
4.如权利要求3所述的基于IMM-KEnPF的风机桨距系统故障诊断方法，其特征在于，步骤4具体包括：
在获得故障模型概率之后，使用故障模型概率对各滤波器的滤波结果进行加权融合，得到IMM-KEnPF的估计结果，如(37)-(39)所示：
5.如权利要求3所述的基于IMM-KEnPF的风机桨距系统故障诊断方法，其特征在于，步骤5具体包括：
对于故障模型概率所蕴含的模式信息，利用式(40)-(41)的决策函数实现工作模式的识别和诊断；当某一个故障模型与当前系统匹配，而其他故障模型都失配时，说明系统正处于对应的故障模式当中：
式中，μ′T为预先设定的诊断阈值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110254185.9A CN112861289B (zh) | 2021-03-09 | 2021-03-09 | 基于IMM-KEnPF的风机桨距系统故障诊断方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110254185.9A CN112861289B (zh) | 2021-03-09 | 2021-03-09 | 基于IMM-KEnPF的风机桨距系统故障诊断方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112861289A CN112861289A (zh) | 2021-05-28 |
| CN112861289B true CN112861289B (zh) | 2021-10-29 |
Family
ID=75994947
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110254185.9A Active CN112861289B (zh) | 2021-03-09 | 2021-03-09 | 基于IMM-KEnPF的风机桨距系统故障诊断方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112861289B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113887055A (zh) * | 2021-10-11 | 2022-01-04 | 西安因联信息科技有限公司 | 一种基于运行参数的离心风机性能退化评估方法及系统 |
| CN117406589A (zh) * | 2023-10-08 | 2024-01-16 | 哈尔滨工业大学 | 一种针对模型未知的机动目标交互式平滑变结构滤波方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102725520A (zh) * | 2009-12-21 | 2012-10-10 | 维斯塔斯风力系统集团公司 | 具有用于执行对风力涡轮发电机的预测控制的控制方法和控制器的风力涡轮机 |
| CN106599368A (zh) * | 2016-11-14 | 2017-04-26 | 浙江大学 | 基于改进粒子提议分布和自适应粒子重采样的FastSLAM方法 |
| CN108681614A (zh) * | 2018-03-07 | 2018-10-19 | 南京航空航天大学 | 基于改进高斯粒子滤波的涡扇发动机突变故障诊断方法 |
| CN110597203A (zh) * | 2019-09-09 | 2019-12-20 | 兰州理工大学 | 一种基于多gpu并行crpf的故障诊断方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11486925B2 (en) * | 2020-05-09 | 2022-11-01 | Hefei University Of Technology | Method for diagnosing analog circuit fault based on vector-valued regularized kernel function approximation |
-
2021
- 2021-03-09 CN CN202110254185.9A patent/CN112861289B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102725520A (zh) * | 2009-12-21 | 2012-10-10 | 维斯塔斯风力系统集团公司 | 具有用于执行对风力涡轮发电机的预测控制的控制方法和控制器的风力涡轮机 |
| CN106599368A (zh) * | 2016-11-14 | 2017-04-26 | 浙江大学 | 基于改进粒子提议分布和自适应粒子重采样的FastSLAM方法 |
| CN108681614A (zh) * | 2018-03-07 | 2018-10-19 | 南京航空航天大学 | 基于改进高斯粒子滤波的涡扇发动机突变故障诊断方法 |
| CN110597203A (zh) * | 2019-09-09 | 2019-12-20 | 兰州理工大学 | 一种基于多gpu并行crpf的故障诊断方法 |
Non-Patent Citations (1)
| Title |
|---|
| 基于修正IMM 的风机变桨系统故障诊断方法;王进花 等;《北京航空航天大学学报》;20200831;第46卷(第8期);第1460-1466页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112861289A (zh) | 2021-05-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112861289B (zh) | 基于IMM-KEnPF的风机桨距系统故障诊断方法 | |
| CN112686464A (zh) | 短期风电功率预测方法及装置 | |
| CN112990556A (zh) | 一种基于Prophet-LSTM模型的用户用电能耗预测方法 | |
| CN105046045B (zh) | 一种基于贝叶斯组合的Web服务QoS预测方法 | |
| CN114399032B (zh) | 一种电能表计量误差预测方法及系统 | |
| Lever et al. | Modelling policies in mdps in reproducing kernel hilbert space | |
| CN111709454B (zh) | 一种基于最优copula模型的多风电场出力聚类评估方法 | |
| CN111237988A (zh) | 地铁车载空调机组控制方法及系统 | |
| CN104036328A (zh) | 自适应风电功率预测系统及预测方法 | |
| CN112819238A (zh) | 基于混沌鸡群优化算法的短期风电功率预测方法 | |
| Chowdhury et al. | Adaptive regulatory genes cardinality for reconstructing genetic networks | |
| CN115659665A (zh) | 一种基于野马优化算法的电压暂降源信号降噪与识别的方法 | |
| CN113449919A (zh) | 一种基于特征和趋势感知的用电量预测方法及系统 | |
| CN114166509A (zh) | 一种电机轴承故障预测方法 | |
| CN116303786B (zh) | 一种基于多维数据融合算法的区块链金融大数据管理系统 | |
| CN116722653A (zh) | 一种电力系统动态状态检测方法及系统 | |
| CN112149896A (zh) | 一种基于注意力机制的机械设备多工况故障预测方法 | |
| CN115718880A (zh) | 一种复杂装备退化阶段的预测方法 | |
| CN115616333A (zh) | 一种配电网线损预测方法及系统 | |
| CN115660038A (zh) | 基于误差因素和改进moea/d-sas的多阶段集成短期负荷预测 | |
| CN115459982A (zh) | 一种电力网络虚假数据注入攻击检测方法 | |
| CN114021465A (zh) | 基于深度学习的电力系统鲁棒状态估计方法及系统 | |
| Husaini et al. | The effect of network parameters on pi-sigma neural network for temperature forecasting | |
| CN116595883B (zh) | 数值反应堆实时在线系统状态修正方法 | |
| CN117910550B (zh) | 一种基于深度学习的无油超高速离心压缩机自动优化系统 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20240527 Address after: 731100 No. 18 Xinzhuang Lane, Suonan Town, Dongxiang Autonomous County, Linxia Hui Autonomous Prefecture, Gansu Province Patentee after: Shengshi Hengrui Construction Group Co.,Ltd. Country or region after: China Address before: 730050, No. 287 Lan Ping Road, Qilihe District, Gansu, Lanzhou Patentee before: LANZHOU University OF TECHNOLOGY Country or region before: China |