关键词
叶根螺栓 时间序列 LSTM
项目背景
以下直接引用了某些论文的内容，帮助大家了解背景。
螺栓连接是风力发电机组装配中的重要装配方式，几乎涉及到风力发电机组的所有部件。因此，螺栓的选用和强度校核是风力发电机组可靠性的重要保证。随着我国风电事业的跨越式发展，伴随着风力发电成本不断下降，风电机组的价格也越来越低，各大风电设备总装企业的价格战已经进行到了白热化阶段。如何在降低成本的情况下，保证风电机组的质量，成为各大风电企业面临的重要问题。
--
现阶段，我国风电机组的螺栓失效问题已经在连接塔筒法兰的高强度螺栓上有所体现。主要失效形式为：安装麦抢带发生滑丝、扭断、屈服、甚至拉断等现象;设备运行过程中发生螺栓断裂，威胁机组运行，严重者甚至造成风力发电机组倒塌。
--
螺栓断裂与以下几种因素有关：
- 螺栓的质量
- 螺栓的预紧力矩
- 螺栓的强度
- 螺栓的疲劳强度
数据及问题描述
该项目是想基于风机scada系统数据通过机器学习的方式对叶根螺栓的断裂进行提前预测，从而避免重大事故和经济损失。
拿到的样本数据中包括33组风机3月、4月、5月、6月共四个月的数据，其中部分数据有少量缺失。故障风机为17#、12#、6#，其他均为正常风机。根据实际数据情况与建模需求选取全部故障风机数据、随机选取适量的正常数据进行建模。具体数据使用情况如下（参与分析风机均已去除停机数据）：
数据使用情况
解决方案
- 经过对数据的深入考察发现，所给数据中，与变桨系统相关的变量及部分温度特征与叶根螺栓故障有明显的相关性，如平均桨叶角度、最高Topbox温度、最高控制柜温度等。理论上基于足够的数据可以训练识别故障的机器学习模型。
- 考虑到叶根螺栓断裂在时间上应该存在某种应力累积的过程，达到一定程度时则会发生断裂。因此尝试基于时间序列的LSTM深度学习模型。为了预测可能发生叶根螺栓断裂的风机大概的断裂时间，考虑搭建两个模型，第一个模型识别是否为可能发生故障的风机并标识蕴含故障模式的数据与正常数据，第二个模型对故障风机预测发生断裂时间。具体思路如下：
- 首先去除停机部分数据，仅考察风机运行数据；
- 然后针对第一阶段的模型，对故障风机停机时间点之前的数据逐条贴标签1，取时间窗口为10天，对故障点t0后的正常运行数据和随机选取的正常风机数据贴标签0，分别表示叶根螺栓故障发生和未发生。
- 对上述数据进行滑窗处理，为使得参与训练数据中正负样本均衡，令正常数据的滑窗步长大于故障数据的滑窗步长。基于上述数据训练一个LSTM分类模型。
- 针对第二阶段的故障模型，对故障风机停机时间点之前的数据逐条贴标签，以递增的整数表示距离该时间点的由近及远，取时间窗口为10天，为提高模型预测效果，相邻数据之间的标签为间隔为10的整数。例如对于10分钟级的数据，某条数据的标签为60则表示距离故障发生时间t0为1小时左右。
- 同样进行滑窗处理，一个滑窗内有144条数据。训练LSTM回归模型。
实施过程
特征处理
考察样本数据，剔除无关变量，选取与叶根螺栓断裂故障相关的特征。经过分析发现如下特征与叶根螺栓断裂故障有一定相关性：平均桨叶角度(deg)、最小桨叶角度(deg)、最大桨叶角度(deg)、最大机舱加速度(g)、最高电机绕组温度(℃)、最高Topbox温度(℃)、最高机舱温度(℃)、最高控制柜温度(℃)、最高变桨电机1温度(℃)、最高变桨电机2温度(℃)、最高变桨电机3温度(℃)、最高变桨柜1的柜体温度(℃)、最高变桨柜2的柜体温度(℃)、最高变桨柜3的柜体温度(℃)、平均变流器功率(kW)、最高变桨柜1备电柜温度(℃)、最高变桨柜2备电柜温度(℃)、最高变桨柜3备电柜温度(℃)。
为便于观察，以下绘制了部分特征在一定时间范围内的变化趋势对比：
topbox温度变化
平均桨叶角度
最高机舱温度
变桨柜1柜体温度
经过对所有特征的遍历分析，选取上述特征或其组合参与模型训练，对最高变桨电机1温度(℃)、最高变桨电机2温度(℃)、最高变桨电机3温度(℃)，最高变桨柜1的柜体温度(℃)、最高变桨柜2的柜体温度(℃)、最高变桨柜3的柜体温度(℃)，最高变桨柜1备电柜温度(℃)、最高变桨柜2备电柜温度(℃)、最高变桨柜3备电柜温度(℃)等特征进行合并（分别取其均值作为新变量代替三个近似特征）。
基于LSTM的故障分类模型
- 模型的结构上，LSTM层设置cell_size 为13，其中有dropout机制，会随机忘记前几层的Cell中神经元，该层设置共13个神经元作为LSTM的输出。DeepLearnning全连接层共设置3层，其中第一层有18个神经元，激活函数为ReLU；第二层有15个神经元，激活函数为ReLU；第三层有12个神经元，激活函数也为ReLU。最终隐层的输出进入输出层，输出层共有1个神经元，输出0/1数值标签。
-
训练方面，采取128Batch Size大小的数据分批投入。由于已在数据滑窗时均衡了正负样本比例，这里直接按4:1的比例随机划分训练集和测试集。
训练过程中的loss下降趋势如下图所示：
损失函数
经过测试模型最终的综合准确率为96%，故障召回率95%。如下图：
report
基于LSTM的故障时间预测模型
- 该模型为回归模型，模型的结构上，LSTM层设置cell_size 为20，num_size为3，DeepLearnning全连接层,同样设置3层，其中第一层有20个神经元，激活函数为ReLU；第二层有15个神经元，激活函数为ReLU；第三层有11个神经元，激活函数也为ReLU。学习率为0.01。
- 最终隐层的输出进入输出层，输出层共有1个神经元，输出距离故障停机时间点的整数值标签，激活函数为Linear线性激活函数，构造出基于LSTM神经网络结构的回归模型。
-
训练方面，采取64Batch Size大小的数据分批投入。Loss的计算采用均方误差。其模型架构如下图所示：
模型结构
Loss值变化如下，模型loss最终收敛于232.86，对于一个Batch Size（10080条数据）来说已经很低，训练效果比较理想。
loss
将实际数据与预测数据对比如下图，对于叶根螺栓故障来说该预测误差是可以接受的。
result
核心代码
这里贴出部分代码(关键代码)
- lstm网络实现（基于tensorflow1.2.0）
# -*- coding: utf-8 -*-
"""
Created on Wed Sep 13 13:20:45 2018
@author: xuanlei
"""
import pandas as pd
import tensorflow as tf
import numpy as np
#==============================================================================
# Batch Normalization
#==============================================================================
def batch_norm_layer(x, train_phase, scope_bn):
with tf.variable_scope(scope_bn):
beta = tf.Variable(tf.constant(0.0, shape=[x.shape[-1]]), name='beta', trainable=True)
gamma = tf.Variable(tf.constant(1.0, shape=[x.shape[-1]]), name='gamma', trainable=True)
axises = np.arange(len(x.shape) - 1)
batch_mean, batch_var = tf.nn.moments(x, axises, name='moments')
ema = tf.train.ExponentialMovingAverage(decay=0.5)
def mean_var_with_update():
ema_apply_op = ema.apply([batch_mean, batch_var])
with tf.control_dependencies([ema_apply_op]):
return tf.identity(batch_mean), tf.identity(batch_var)
mean, var = tf.cond(train_phase, mean_var_with_update,
lambda: (ema.average(batch_mean), ema.average(batch_var)))
normed = tf.nn.batch_normalization(x, mean, var, beta, gamma, 1e-3)
return normed
#==============================================================================
# RNN Structure
#==============================================================================
class LSTMRNN():
#initial setting
def __init__(self, n_steps, input_size, output_size, cell_size, h1_size, h2_size, h3_size, LR,num_size, batch_size):
self.n_steps = n_steps
self.input_size = input_size
self.output_size = output_size
self.cell_size = cell_size
self.num_size = num_size
self.h1_size = h1_size
self.h2_size = h2_size
self.h3_size = h3_size
self.batch_size = batch_size
self.LR = LR
self.num_size = num_size
with tf.name_scope('inputs'):
self.xs = tf.placeholder(tf.float32, [None, n_steps, input_size], name='xs')
self.ys = tf.placeholder(tf.float32, [None, output_size], name='ys')
self.keep_prob = tf.placeholder(tf.float32, name='keep_prob')
self.train_phase = tf.placeholder(tf.bool, name='train_phase')
with tf.name_scope('in_hidden'):
self.add_input_layer()
with tf.name_scope('LSTM_Cell'):
self.add_cell_layer()
with tf.name_scope('hidden_1'):
self.add_h1_layer()
with tf.name_scope('hidden_2'):
self.add_h2_layer()
with tf.name_scope('hidden_3'):
self.add_h3_layer()
with tf.name_scope('out_hidden'):
self.add_output_layer()
with tf.name_scope('cost'):
self.compute_cost()
with tf.name_scope('train'):
self.train_op = tf.train.AdamOptimizer(learning_rate=self.LR).minimize(self.cost)
def add_input_layer(self):
with tf.name_scope('input_layer'):
l_in_x = tf.reshape(self.xs,[-1,self.input_size], name='x_input')
#Ws_in = tf.Variable(tf.truncated_normal([self.input_size, self.cell_size], mean=1, stddev=0.5))
Ws_in = tf.get_variable("W", shape=[self.input_size, self.cell_size],initializer=tf.contrib.layers.xavier_initializer())
bs_in = tf.Variable(tf.zeros([self.cell_size,])+0.01)
l_in_y = tf.matmul(l_in_x,Ws_in)+bs_in
self.l_in_y = tf.reshape(l_in_y,[-1,self.n_steps,self.cell_size],name='cell_input')
def add_cell_layer(self):
with tf.name_scope('LSTM_layer'):
lstm_cell = tf.contrib.rnn.DropoutWrapper(tf.contrib.rnn.BasicLSTMCell(self.cell_size, forget_bias=1.0, state_is_tuple=True), output_keep_prob=self.keep_prob)
lstm_cells = tf.contrib.rnn.MultiRNNCell([lstm_cell]*self.num_size, state_is_tuple=True)
self.cells_init_state = lstm_cells.zero_state(self.batch_size,dtype=tf.float32)
self.cells_outputs, self.cells_final_state = tf.nn.dynamic_rnn(lstm_cells, self.l_in_y, initial_state=self.cells_init_state, time_major=False)
tf.summary.histogram('cells_outputs',self.cells_outputs)
tf.summary.histogram('cells_final_state',self.cells_final_state)
def add_h1_layer(self):
with tf.name_scope('h1_layer'):
h1_x = tf.reshape(self.cells_outputs, [-1,self.cell_size])
#Ws_h1 = tf.Variable(tf.truncated_normal([self.cell_size, self.h1_size], mean=3, stddev=1))
Ws_h1 = tf.get_variable("W1", shape=[self.cell_size, self.h1_size],initializer=tf.contrib.layers.xavier_initializer())
bs_h1 = tf.Variable(tf.zeros([self.h1_size,])+0.01)
non_bn_h1 = tf.nn.relu(tf.matmul(h1_x,Ws_h1)+bs_h1)
non_bn_h1 = tf.matmul(h1_x,Ws_h1)+bs_h1
self.h1_y = batch_norm_layer(non_bn_h1, train_phase=self.train_phase, scope_bn='bn_h1')
def add_h2_layer(self):
with tf.name_scope('h2_layer'):
h2_x = tf.reshape(self.h1_y, [-1,self.h1_size])
#Ws_h2 = tf.Variable(tf.truncated_normal([self.h1_size, self.h2_size], mean=3, stddev=2))
Ws_h2 = tf.get_variable("W2", shape=[self.h1_size, self.h2_size],initializer=tf.contrib.layers.xavier_initializer())
bs_h2 = tf.Variable(tf.zeros([self.h2_size,])+0.01)
non_bn_h2 = tf.nn.relu(tf.matmul(h2_x,Ws_h2)+bs_h2)
non_bn_h2 = tf.matmul(h2_x,Ws_h2)+bs_h2
self.h2_y = batch_norm_layer(non_bn_h2, train_phase=self.train_phase, scope_bn='bn_h2')
def add_h3_layer(self):
with tf.name_scope('h3_layer'):
h3_x = tf.reshape(self.h2_y, [-1,self.h2_size])
#Ws_h3 = tf.Variable(tf.truncated_normal([self.h2_size, self.h3_size], mean=3, stddev=2))
Ws_h3 = tf.get_variable("W3", shape=[self.h2_size, self.h3_size],initializer=tf.contrib.layers.xavier_initializer())
bs_h3 = tf.Variable(tf.zeros([self.h3_size,])+0.01)
non_bn_h3 = tf.nn.relu(tf.matmul(h3_x,Ws_h3)+bs_h3)
non_bn_h3 = tf.matmul(h3_x,Ws_h3)+bs_h3
self.h3_y = batch_norm_layer(non_bn_h3, train_phase=self.train_phase, scope_bn='bn_h3')
def add_output_layer(self):
layer_name='output_layer'
with tf.name_scope('output_layer'):
l_out_x = tf.reshape(self.h3_y,[-1,self.h3_size],name = 'y_input')
#Ws_out = tf.Variable(tf.truncated_normal([self.h3_size, self.output_size], mean=3, stddev=1))
Ws_out = tf.get_variable("W4", shape=[self.h3_size, self.output_size],initializer=tf.contrib.layers.xavier_initializer())
bs_out = tf.Variable(tf.zeros([self.output_size,]))
#self.pred = tf.matmul(l_out_x,Ws_out)+bs_out
self.pred = tf.nn.sigmoid(tf.matmul(l_out_x,Ws_out)+bs_out)
#self.pred = tf.matmul(l_out_x,Ws_out)+bs_out
#self.pred = tf.nn.softmax(tf.matmul(l_out_x,Ws_out)+bs_out)
tf.summary.histogram('w', Ws_out)
tf.summary.histogram('b', bs_out)
tf.summary.histogram('out', self.pred)
#===============================================================================
# 交叉熵
#===============================================================================
def compute_cost(self):
with tf.name_scope('loss'):
self.cost = -tf.reduce_sum(self.ys*tf.log(self.pred+0.001)+(1-self.ys)*tf.log(1-self.pred+0.001))#可以调节权重控制样本不平衡的数据训练
tf.summary.scalar('result_cost', self.cost)
小结
- 值得指出的是故障记录中的故障发生时间并非实际的叶根螺栓断裂时间而是发现时间，尽管真正的断裂时间对模型更有意义但实际中难以捕捉，因此若进一步提升模型的实际效果，可将贴标签的基准时间往前顺延。
- 尽管样本数据中存在与叶根螺栓断裂故障相关的特征，模型效果较为理想，但并不能排除其他故障或风机自身特性的影响。要训练一个稳健的、高鲁棒性的模型还需要更多的故障数据。
