基于Bagging集成策略和多元状态估计的风电机组齿轮箱状态监测

2020-08-03赵劲松王梓齐刘长良

科学技术与工程 2020年20期

赵劲松，王梓齐，刘长良,2*

(1.华北电力大学控制与计算机工程学院,保定 071003；2.新能源电力系统国家重点实验室(华北电力大学),北京 102206)

中国风电机组的主要机电设备中,齿轮箱的故障率位列第四且维修费用和时间也相对较高[1-2]。因此,有必要实时监测风电机组齿轮箱的运行状态,及早发现潜在故障并提前安排维护,从而降低运维成本。

监控与数据采集(supervisory control and data acquisition, SCADA)系统被大量应用于风电场中,主要功能是储存历史运行数据和在线采集实时运行数据。相比其他常见的油液、振动分析等,基于SCADA运行数据的状态监测方法无须新加装金属磨粒、加速度传感器等,在经济性上的优势相对更大，近年来受到了较为广泛的关注[3]。

多元状态估计(multivariate state estimate technique, MSET)[4]是根据正常运行数据构成记忆矩阵,进而对设备当前运行数据进行估计的一种基于SCADA数据的状态监测方法,有着无须训练参数、精度高、泛化能力强等优势。近年来,MSET在风电机组齿轮箱[5-6]、火电机组磨煤机[7]和引风机[8]等设备中得到了一定应用。MSET对记忆矩阵的依赖大,要想进一步提高精度,往往需要对记忆矩阵进行扩展。但是,由于每个计算周期均需要执行时间复杂度较高的矩阵乘法运算,扩展记忆矩阵会对MSET的实时性造成较大影响[9]。

集成学习是一种通用的机器学习模型结构,即用多个差异化的个体学习器对同一问题进行训练,再将各学习器的结果综合后输出。对于相同的机器学习算法,相比建立单一的全局学习器,集成学习模型通常有着更强的泛化能力。近年来,集成学习已在风电功率预测[10-11]、齿轮箱故障诊断[12]等相关领域得到了一定应用。

为解决记忆矩阵较大时,常规MSET实时计算速度慢的问题,提出了一种基于Bagging集成策略和MSET的状态监测方法。首先基于Bagging对数量较多的训练样本进行随机抽样,构造多个记忆矩阵规模较小的MSET个体学习器,最后将各MSET的结果平均后输出。随着个体学习器数量的增多,集成MSET对训练样本的利用率增大、估计精度提高且计算时间仅为线性增长,能够在对精度影响不大的前提下,有效改善算法的实时性。

以某2 MW风电机组SCADA系统的实际运行数据为算例,基于正常运行数据将常规与集成MSET的性能进行了对比；进一步基于故障数据,并结合统计过程控制技术设计预警阈值,使用集成MSET实现了对齿轮箱实际故障的预警。

1 算法基础及分析

1.1 多元状态估计(MSET)

MSET的核心是基于正常运行数据构造的记忆矩阵,通过记忆矩阵对设备当前的运行状态进行估计:若设备处于正常状态,则对当前运行数据的估计精度高；若设备可能发生潜在故障,则估计精度相对较低。MSET的原理和实现过程如下。

某设备共有n个参数参与状态监测,将某一时刻这些参数的测量值记为观测向量X=[x1,x2, …,xn]T。在设备正常运行时,采集的m个观测向量X构造记忆矩阵D如下(要求m个观测向量尽可能覆盖设备所有的运行工况)：

(1)

MSET的输入为当前时刻的观测向量Xo,输出为观测向量Xo的估计值Xe。对于每个输入Xo,均需要计算一个对应的权值向量W=[w1,w2, …,wm]T,使得记忆矩阵D和权值向量W的乘积等于输出Xe。具体公式如下:

Xe=DW=w1X1+w2X2+…+wmXm

(2)

权值向量W的计算公式如下(推导过程见文献[13]):

W=(DT⊗D)-1·(DT⊗Xo)

(3)

式(3)中：⊗为非线性运算符,常选取为两向量间的欧几里得距离。

将式(3)代入式(2)中,输出Xe的计算公式如下:

Xe=D·(DT⊗D)-1·(DT⊗Xo)

(4)

MSET不需要实现基于训练样本进行参数训练,但在计算每一个输出Xe时,都需要重新基于记忆矩阵D计算权值向量W。因此,记忆矩阵D对MSET的估计精度和计算速度有着决定性的影响,具体分析如下。

(1)记忆矩阵D的规模较小时,MSET的计算速度较快,但对设备正常运行状态空间的覆盖有限,估计精度相对较低，状态监测的准确性较差。

(2)记忆矩阵D的规模较大时,MSET的估计精度相对较高，准确性较好。但在计算权重W和输出Xe时,需要分别进行矩阵乘法运算[矩阵(DT⊗D)-1和(DT⊗Xo)相乘、矩阵D和W相乘],而矩阵乘法的时间复杂度一般为o(n3)左右[14]。因此,扩大记忆矩阵D会使计算速度大幅度增加,同样也会占用更多的存储空间。

根据上述分析,从原理上MSET的估计精度和实时性之间存在较大的矛盾,有必要针对该问题进行深入研究,在保证一定估计精度的同时尽可能提高MSET的实时计算速度。

1.2 集成MSET和Bagging集成策略

易证,计算1次n阶矩阵乘法所花费的时间,要大于计算k(k≥2)次n/k阶矩阵乘法所花费的时间。将上述结论推广至MSET并结合集成学习思想,可由多个记忆矩阵规模较小的MSET模型作为个体学习器,得到在计算速度上优于常规方法、估计精度与常规方法相当的集成MSET模型,从而改善算法的实时性。

集成学习的常见构造策略有Bagging、Boosting和Stacking等[15]。其中,Boosting的个体学习器是串行关系,即根据前一个学习器的训练结果调整权重后训练下一个学习器；Stacking要求个体学习器是异质的,即不同学习器所使用的算法不同,均很难与MSET进行结合。

Bagging集成策略的主要思想是对原始训练样本进行多次随机抽样,并行生成多组训练样本以训练多个个体学习器[16],能够与MSET进行方便地结合。Bagging集成策略的具体结构和实现流程如图1所示。

图1 Bagging集成学习策略Fig.1 The Bagging ensemble learning strategy

由图1可知,Bagging集成策略主要分为训练和测试两个阶段。

(1)训练阶段:首先需要确定个体学习器的个数n和对应的训练样本容量(一般小于等于原始训练样本个数)；之后,对原始训练样本进行n次随机抽样,得到n组训练样本并训练得到n个基于相同算法的个体学习器。

(2)测试阶段:基于测试样本计算n个个体学习器的输出,对于分类问题采取投票方式，对于回归问题采取算数平均方法,将n个输出综合后作为集成学习模型整体输出。

基于Bagging集成策略建立的集成MSET模型,其估计精度和计算速度受每个MSET个体学习器的记忆矩阵规模和个体学习器的个数决定:记忆矩阵规模越大、个体学习器的个数越多,集成MSET模型的估计精度越高、计算速度越慢；在记忆矩阵规模确定的情况下,集成MSET的计算时间与个体学习器的个数呈线性关系。

2 基于Bagging和MSET的状态监测

2.1 算法流程

针对风电机组齿轮箱这一研究对象,基于Bagging和MSET的状态监测方法的具体实现流程如图2所示。

由图2可知,算法主要分为3部分：数据预处理、集成模型建立和设计预警阈值。

图2 状态监测算法的实现流程Fig.2 Implementation process of the condition monitoring method

(1)数据预处理的第一步是对SCADA历史数据中属性缺失、参数异常的数据进行删除,以获取代表机组正常运行状态的数据。第二步是确定构成MSET观测向量X的运行参数,首先确定能够直接反映齿轮箱运行状态的状态变量,通常只选取1个,记作y；其次确定与状态变量密切相关的其他运行参数即辅助变量,通常选取多个,记作xi。第三步是将数据分为训练集和验证集,训练集用来建立集成MSET模型、样本数量相对较多,验证集用来设计预警阈值、样本数量相对较少。

(2)集成MSET模型的建立过程与Bagging集成策略的训练过程基本一致。由于MSET不需要训练参数,通过随机抽样确定训练样本后,构造n个记忆矩阵D即完成了集成MSET模型的建立。

2.2 基于统计过程控制技术的预警阈值设计

统计过程控制(statistical process control, SPC)技术是一种基于统计理论的过程监测方法,近年来在风电机组状态监测的预警阈值设计中得到了一定应用[17-18],具体的阈值设计方法如下。

设随机变量X～N(μ,σ2),则根据正态分布相关理论,随机变量X落在(μ-2σ,μ+2σ]和(μ-3σ,μ+3σ]区间内的概率如下:

P(μ-2σ

(5)

P(μ-3σ

(6)

若X的取值长期超出这两个区间,可以认为过程受到了异常因素的影响,导致随机变量的X概率分布出现了显著变化。

(7)

式(7)中：ei为验证集上状态变量的预测残差；m为验证集的样本数。

(8)

预警阈值构造完成后,结合状态变量的预测残差,即可初步实现风电机组齿轮箱的状态监测、故障预警及可视化。

3 算例及分析

3.1 数据来源及预处理

数据来源为福建省某陆上风电场的一台双馈式机组。机组的切入风速为4 m/s,切出风速为25 m/s,额定功率为2 MW,SCADA系统的采样周期为10 min。该机组于2016年7月13日10:20发生齿轮箱故障导致非计划停运,故障发生前SCADA系统未提前发出报警信号。

从历史运行数据库中,导出2月13日10:20至齿轮箱故障前的运行数据。剔除这些数据中运行参数缺失、有功功率小于等于零、风速小于4 m/s或大于25 m/s的数据点,并进一步基于拉依达准则去除运行参数异常的数据点,最终用余下的10 000组数据参与接下来的实验。

从运行参数中,选择风速、发电机转速、环境温度、有功功率和齿轮箱轴承温度共5个参数构成观测向量X,选择齿轮箱轴承温度作为齿轮箱状态监测的状态变量y。参与实验的10 000组数据中,上述运行参数的变化范围如表1所示。

表1 运行参数的变化范围Table 1 Range of the operation parameters

为了去除量纲不同对结果造成的影响,将运行参数进行min-max归一化:

(9)

式(9)中：xij为参数j的第i个数据点；maxxj和minxj分别为参数j的最大值和最小值；x′ij为归一化后的数据。

min-max归一化后,以1～5000号样本作为训练集,用于建立常规和集成MSET模型；5001～6000号作为验证集,用来对比验证常规和集成MSET的性能并计算预警阈值；5001～10000号作为测试集(包含验证集),用来进行齿轮箱故障的状态监测实验。

3.2 状态监测方法的性能对比实验

分别以训练集中时间标签靠后的1 000、2 000～5 000组样本,构造常规MSET的记忆矩阵并基于验证集进行测试。使用均方根误差RMSE作为误差函数以评价估计精度,计算公式为

(10)

程序的运行环境为MATLAB R2018a,CPU型号为Intel i5-7500,内存为8 G。常规MSET在验证集上的测试结果如图3所示。图3中,计算时间的含义为计算1 000组验证集所花费的时间。

图3 常规MSET方法的测试结果Fig.3 Test results of conventional MSET method

由图3可知,从估计精度的角度,随着记忆矩阵规模的扩大,常规MSET的误差有一定下降,但矩阵规模大于2 000后降幅较小,说明扩展记忆矩阵能够改善MSET的估计精度,但矩阵达到一定规模后,继续扩展对精度的改善较为有限；从计算时间的角度,随着记忆矩阵规模的扩大,常规MSET的计算时间大幅度增长且为非线性增长,说明较大规模的记忆矩阵会严重影响MSET的实时性。

下面进行常规和集成MSET的对比实验。确定MSET个体学习器的记忆矩阵规模为1 000,分别建立由1～5个个体学习器构成的集成MSET模型。由于Bagging集成策略中训练样本的选取具有随机性,故在验证集上进行50次重复实验。50次实验所得到的RMSE的箱形图如图4所示。图4中,横轴记忆矩阵规模的含义为MSET个体学习器的记忆矩阵规模之和。

图4 集成MSET方法50次重复实验的测试结果Fig.4 Test results of 50 repetitive experiments of ensemble MSET method

由图4可知,随着个体学习器个数的增多,集成MSET模型的误差逐渐减小且分布范围逐渐收窄,说明集成MSET能够保证一定的估计精度,且个体学习器的数目越多,估计精度相对更高，在重复测试中的表现更稳定。尽管集成MSET可能需要进行多次重复测试,才能得到估计精度相对较高的模型,但该过程属于离线训练阶段,不会对在线计算的实时性产生影响。

不同记忆矩阵规模下,常规MSET和集成MSET的估计精度和计算时间如表2所示。表2中,t为计算1 000组验证集所花费的时间,集成MSET的t为50次重复实验的平均计算时间；Q3、Q2、Q1分别为集成MSET的RMSE在50次重复实验中的上四分位数、中位数和下四分位数。

表2 集成MSET和常规MSET的对比Table 2 Comparison of ensemble MSET and conventional MSET

由表2可知,随着个体学习器个数的增多,集成MSET的计算时间呈线性增长；当个体学习器为3

个时,RMSE的中位数Q2基本与矩阵规模为2 000的常规MSET相当,但计算时间仅为常规方法的60%,说明在估计精度接近的前提下,集成MSET的计算速度相比常规方法有了较大改善。

3.3 齿轮箱实际故障实验

使用5 000组测试集数据,验证集成MSET的故障预警能力。选择集成MSET的个体学习器个数为3,总的记忆矩阵规模为3 000,模型在验证集上的RMSE为0.003 6。

图5 集成MSET方法在测试集上的残差Fig.5 Residual of ensemble MSET method on test set

由图5可知,第2500点前,残差总体的幅值相对较小、超过预警阈值的数据相对较少,可以认为齿轮箱的运行状态相对正常；第2500点后,残差的波动程度相对剧烈、超过预警阈值的数据明显增多,说明齿轮箱的运行状态发生了较大程度的改变。

通过图5能够大致判断齿轮箱的运行情况,但由于预测残差具有较强的随机性,提供给运行人员的信息较为杂乱、直观性差。为改善以上问题,结合滑动窗口法和预警阈值,设计了一种名为滑动窗口异常率的状态监测指标,具体如下。

对于预警阈值T1和T2,分别定义一级异常率R1和二级异常率R2:

(11)

式(11)中：M为滑动窗口的长度；N1和N2为当前滑动窗口中,残差超过预警阈值T1和T2的数据个数。

取滑动窗口的长度M=1 000,计算测试集上的滑动窗口异常率如图6所示。

图6 测试集上的滑动窗口异常率Fig.6 Sliding window exception rate on test set

由图6可知,第1 500个滑动窗口前,异常率相对较小,R1小于6.5%、R2小于2.5%；第1 500～2 500个滑动窗口,异常率不断上升直至达到峰值,R1和R2高达26.0%和15.8%,说明齿轮箱可能发生了早期故障,导致运行状态发生了显著变化；第2 500个滑动窗口后,异常率的幅值和波动情况均有一定程度的回落,但仍显著高于正常运行状态,R1维持在12.5%以上、R2维持在5.2%以上。根据异常率的变化曲线,可在第2 500个滑动窗口即故障发生前10 d左右,判定齿轮箱发生了较严重的早期故障。

相比图5的预测残差序列,图6所示的滑动窗口异常率曲线平滑度较好，能够提供清晰、直观的状态监测信息,更有利于运行人员的分析和决策。