风机开裂故障预警比赛经历

背景

首先介绍一下比赛背景。
这个比赛是中电投的一个关于风机开裂故障分析的预警的比赛。

训练数据有将近5万个样本，测试数据有将近9万个样本。数据来自于SCADA采集系统。采集了10分钟之内的75个特征值得数据信息，label是一周以内风机是否会发生故障的label。

数据介绍

每个样本10分钟之内大概采集到了450条数据，一共75个特征，也就是差不多75*450个信息。最后三个特征完全没有数据，所以一开始做的时候，我们就把最后三个特征进行删除，实际上我们是对72个特征进行的数据分析。

最开始，用的是seaborn画的正常风机和不正常风机的频率分布图，比如说对于轮毂转速这个特征：

import seaborn as sns
import pandas as pd
data_file = r"D:\fan_fault\feature1.csv"
pre_process = pd.read_csv(data_file, encoding = "gbk")

pre_process = pre_process.fillna(0)
feature1_plot = pre_process["normal(0)"]

feature2_plot2 = pre_process["fault(1)"]
sns.kdeplot(feature1_plot,  shade = True)
sns.kdeplot(feature2_plot2,  shade = True)

大部分特征都是这样，没有很好地区分度。正是因为这样，也一直没有尝试出来非常好的模型。后来我们尝试用MATLAB画图，每个特征出两个图

看起来要比seaborn好一些。（后两个图和第一个不是一个特征）

我们在做数据分析这一块一个很大的问题是在于只去查了各个特征的物理含义，做了频率和频数分布图，看看是否有没有好的特征，然后就直接进入了下一步。忘了考虑是否可能会出现因为采集的问题而导致的异常值和空缺值的问题。这一点导致后面我们很多工作推到从来。

数据分析

我们从统计上来看，并没有找到很好区分度的特征，然后就考虑从物理上来找。在刘老师的建议下，我们尝试了有轮毂转速，风速为6.5m/s时，y方向振动值的特征：

依旧没有很好的区分度，对于其他风速尝试，也是如此。

我们讨论阈值、记0等方式构造新特征。就在说道记0这个新特征构造办法的时候，突然发现，大气压力这个特征，居然有0的情况。根据物理学的知识来讲，风机的大气压力是不可能为0的。然后我们才想起来，没有对数据的异常值进行处理。删除了有8万多条整行全为0的数据，导致某些文件为空，也就是这个风机没有数据信息。当然，也有某些风机是某几行为0。

除了删除空缺值，我们还对其他明显是异常的数据进行了一些数据清洗工作。因为之前我们对于数据特征数统计分析是根据未清洗的数据做的分析，所以分析的可靠性也有点问题，后面我们在做工作的时候，有了一些不必要的麻烦。

我们也做了一些相关性分析的工作，大部分特征相关性十分的高。几十个特征两两组合然后进行相关性分析，会有数千个结果，相关性分析没有办法进行下去。后来，我们就没有考虑相关性的事情。

特征工程

我们最开始尝试对前72个特征构造均值，作为基准尝试看看效果如何

import os
import pandas as pd
import numpy as np
import csv

label_file = r"C:\fan_fault\train\trainX"

train_mean = r"D:\fan_fault\train_mean_new.csv"

with open(train_mean, "a", newline = '', encoding = "utf-8") as f:
    train_mean = csv.writer(f)

    for x in range(1, 48340):
        fan_file = os.path.join(label_file, str(x) +  ".csv")
        print("程序运行进度为", x/48340)    #用该语句查看工作进度状态
        
        with open(fan_file, encoding='utf-8') as f:
            feature_read = pd.read_csv(f)
                
            #遍历打开文件的每一个特征（72），求取均值
            # a用来临时存放计算好的特征均值，外加一个label

            a = []
            for i in range(72):
                mean_num = feature_read.iloc[:, i]
                mean_num = np.array(mean_num).mean()  #生成每个特征所有数据对应的均值
                a.append(mean_num)
            
            train_mean.writerow(a)

也包括绝对值差分累计、差分均值、差分方差，用随机森林进行调参

# -*- coding: utf-8 -*-
"""
Created on Mon Aug  6 13:26:14 2018

@author: jinqiu
"""

import numpy as np
import pandas as pd
from sklearn.preprocessing import MinMaxScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
from sklearn import metrics
from sklearn.model_selection import GridSearchCV

#数据导入、检查空缺值
data = pd.read_csv(r'D:\next\8_19\train_data.csv',encoding = "gbk")
label = pd.read_csv(r"D:\next\8_19\train_label.csv")
data.info()
data.notnull().sum(axis=0)/data.shape[0]
train = data.iloc[:,:-1]
label = label.iloc[:,-1]

#数据标准化
scaler = MinMaxScaler()
train = scaler.fit(train).transform(train)

#单个分类器
clf = RandomForestClassifier(random_state=14)
f1 = cross_val_score(clf, train, label, scoring='f1')
print("f1:{0:.1f}%".format(np.mean(f1)*100))

#调参
parameter_space = {
        'n_estimators':range(10,200,10),
        'max_depth':range(1,10),
        'min_samples_split':range(2,10),
        }
clf = RandomForestClassifier(random_state=14)
grid = GridSearchCV(clf,parameter_space,scoring='f1', n_jobs = 6)
grid.fit(train,label)
print("f1:(0:.1f)%".format(grid.best_score_*100))
print(grid.best_estimator_)

#调参后的分类器
new_clf = RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',
            max_depth=7, max_features='auto', max_leaf_nodes=None,
            min_impurity_decrease=0.0, min_impurity_split=None,
            min_samples_leaf=1, min_samples_split=7,
            min_weight_fraction_leaf=0.0, n_estimators=10, n_jobs=1,
            oob_score=False, random_state=14, verbose=0, warm_start=False)
print("f1:{0:.1f}%".format(np.mean(f1)*100))

测试集输出预测结果

#数据标准化
scaler = MinMaxScaler()
train = scaler.fit(train).transform(train)
test = scaler.fit(test).transform(test)

#训练分类器
clf = RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',
            max_depth=8, max_features='auto', max_leaf_nodes=None,
            min_impurity_decrease=0.0, min_impurity_split=None,
            min_samples_leaf=1, min_samples_split=5,
            min_weight_fraction_leaf=0.0, n_estimators=20, n_jobs=5,
            oob_score=False, random_state=14, verbose=0, warm_start=False)
clf = clf.fit(train, label)
#预测结果
pre = clf.predict(test)

#测试结果文件写入
import csv

label = r"D:/fan_fault/label.csv"

with open(label, "w", newline = '', encoding = "utf-8") as f:
    label = csv.writer(f)

    for x in range(len(pre)):
        label.writerow(pre[x:x+1])

测试效果来看，并没有取得十分理想的效果。

之后想起来没有考虑数据清洗，之前的工作全部推倒从来。我们在此期间，也调查整理关于这72个特征和风机叶片开裂故障的原因，发现从文献上没有找到和叶片开裂故障有很好相关性的特征，我们对于特征进行一些改造，也没有很好的区分效果。后期我们咨询过金风科技的工程师和阜特科技的工程师，他们说这些特征中也没有十分强相关性的特征。我们就考虑到特征之间两两交叉相乘看看效果。

# 交叉特征有2844列，分别是自身的平方和相互交叉，最后求均值方差，最后三个特征不单独再生成交叉特征

import os
import pandas as pd
import numpy as np
import csv
from sklearn.preprocessing import PolynomialFeatures

label_file = r"F:\User\Xinyuan Huang\train_labels.csv"
fan_folder = r"F:\User\Xinyuan Huang"
read_label = pd.read_csv(label_file)

cross_var = r"F:\User\Xinyuan Huang\CaiJi\Feature_crosses\cross_var.csv"

with open(cross_var, "a", newline = '', encoding = "utf-8") as f:
    cross_var = csv.writer(f)
    
    # 该for循环用于定位要打开的文件
    for x in range(len(read_label)-1):
        column1 = str(read_label["f_id"][x:x+1])      #遍历DataFrame第一列的f_id标签下面的每一个数
        column2 = str(read_label["file_name"][x:x+1]) #遍历DataFrame第二列的file_name标签下面的每一个数
        column3 = str(read_label["ret"][x:x+1])       #遍历DataFrame第三列的ret标签下面的每一个数

        f_id = column1.split()[1]  #第一行的文件所对应的f_id进行切片操作，获取对应的数字
        # 对f_id进行补0操作
        f_id = f_id.zfill(3) # 比如2补成002，所以这里写3
        file_name = column2.split()[1] #第一行的文件所对应的file_name
        label = column3.split()[1]     #第一行文件所对应的ret

        fan_file = os.path.join(fan_folder, "train", f_id,  file_name)
        print("程序运行进度为", x/(len(read_label)-1))    #用该语句查看工作进度状态
        


        # 打开相应的fan_file文件进行读取操作
        with open(fan_file, encoding='utf-8') as f:
            dataset = pd.read_csv(f)
            
            #数据集名称为dataset
            poly = PolynomialFeatures(degree=2, include_bias=False, interaction_only=False)
            X_ploly = poly.fit_transform(dataset)


            data_ploly = pd.DataFrame(X_ploly, columns=poly.get_feature_names())

            new_data = data_ploly.ix[:,75:-6]


            #ploly_mean,ploly_var为交叉特征均值方差
            ploly_mean = np.mean(new_data)
            ploly_var = np.var(ploly_mean)

            ploly_var = list(ploly_var)
            ploly_var.append(label)

          
            cross_var.writerow(ploly_var)

交叉相乘之后的文件有数千个特征，生成的文件有将近2G大小。考虑到服务器性能不高，计算旷日持久。不对特征进行筛选，直接进行交叉之后跑算法这条路被我们放弃了。

后来阜特科技的杨工帮我们筛选了一些比较重要的特征，我们在此基础之上进行了一些特征交叉和重要性排序的操作，特征缩小到了几百个（包含交叉、均值、方差等，经过重要性排序），然后用它来跑的模型。

特征里面有一些特征是离散特征，对于这些特征我们进行单独处理，进行离散化。比如说偏航要求值总共有3个值分别是1,2,3。我们对其进行离散化处理，一个特征就变成了三个特征。每个样本统计出现这三个特征的频率。

import os
import pandas as pd
import numpy as np
import csv

label_file = r"E:\8_19\testX_csv"

train_mean = r"E:\8_19\disperse\discrete56.csv"

with open(train_mean, "a", newline = '', encoding = "utf-8") as f:
    train_mean = csv.writer(f)

    for x in range(1, 451):
        fan_file = os.path.join(label_file, str(x) +  ".csv")
#        print("程序运行进度为", x/451)    #用该语句查看工作进度状态
        
        with open(fan_file, encoding='utf-8') as f:
            feature_read = pd.read_csv(f, header = None)
            
            num1 = 0
            num2 = 0
            num3 = 0
            
            a = []
            
            for x in range(len(feature_read)):
                if feature_read[55][x] == 0:
                    num1 = num1+1
                if feature_read[55][x] == 1:
                    num2 = num2+1
                if feature_read[55][x] == 2:
                    num3 = num3+1

            num1 = num1/len(feature_read)
            num2 = num2/len(feature_read)
            num3 = num3/len(feature_read)

            a.append(num1)
            a.append(num2)
            a.append(num3)

            train_mean.writerow(a)

算法

我们最后主要用的算法是Xgboost，期间也尝试过LightGBM，因为算力不够的原因，没有办法尝试一些算法（包括杨工说的SVM以及深度学习的想法），最后主要用Xgboost进行调参，直接一起调参的话算力不够，我们是单个调参以及两两调参组合的形式进行参数调整。

from xgboost import XGBClassifier
import xgboost as xgb

import pandas as pd 
import numpy as np

from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import StratifiedKFold

from sklearn.metrics import log_loss
from sklearn.preprocessing import MinMaxScaler


#数据导入、检查空缺值
data = pd.read_csv(r'D:\next\8_19\train_data.csv',encoding = "gbk")
label = pd.read_csv(r"D:\next\8_19\train_label.csv")
test = pd.read_csv(r"D:\next\8_19\test_data.csv", encoding = "gbk")
train = data.iloc[:,:-1]
label = label.iloc[:,-1]

X_train = train
y_train = label

#数据标准化
scaler = MinMaxScaler()
train = scaler.fit(train).transform(train)
test = scaler.fit(test).transform(test)

#交叉验证
kfold = StratifiedKFold(n_splits=5, shuffle=True, random_state=3)

param_test1 = {
 'max_depth':list(range(3,10,1)),
 'min_child_weight':list(range(1,6,1))
}
gsearch1 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=400, max_depth=5,
 min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,
 objective= 'binary:logistic', nthread=4, scale_pos_weight=1, seed=27), 
 param_grid = param_test1, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch1.fit(X_train,y_train)
gsearch1.grid_scores_, gsearch1.best_params_, gsearch1.best_score_

调完这个参数之后，把最好的输出结果拿出来放在下一个参数调优里进行调整：

aram_test1 = {
 'learning_rate':[i/100.0 for i in range(6,14,2)]
}
gsearch1 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=400, max_depth=6,
 min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,
 objective= 'binary:logistic', nthread=6, scale_pos_weight=1, seed=27), 
 param_grid = param_test1, scoring='roc_auc',n_jobs=-1,iid=False, cv=5)
gsearch1.fit(X_train,y_train)
gsearch1.grid_scores_, gsearch1.best_params_, gsearch1.best_score_

param_test1 = {
 'subsample':[0.8, 0.9]
}
gsearch1 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=310, max_depth=6,
 min_child_weight=1, gamma=0, subsample=0.9, colsample_bytree=0.8,
 objective= 'binary:logistic', nthread=6, scale_pos_weight=1, seed=27), 
 param_grid = param_test1, scoring='roc_auc',n_jobs=-1,iid=False, cv=5)
gsearch1.fit(X_train,y_train)
gsearch1.grid_scores_, gsearch1.best_params_, gsearch1.best_score_

还可以调整好几个参数，最后结果输出：

import xgboost as xgb
dtrain=xgb.DMatrix(X_train ,label=y_train)
dtest=xgb.DMatrix(test)

params={
    'objective': 'binary:logistic',
    'max_depth':6,
    'subsample':0.8,
    'colsample_bytree':0.8,
    'min_child_weight':1,
    'seed':27,
    'nthread':6,
    'learning_rate':0.1,
    'n_estimators':292,
    'gamma':0,
    'scale_pos_weight':1}

watchlist = [(dtrain,'train')]

bst=xgb.train(params,dtrain,num_boost_round=100,evals=watchlist)

ypred=bst.predict(dtest)

import csv

test_label = r"D:\next\8_20\test_label_new.csv"
with open(test_label, "a", newline = '', encoding = "utf-8") as f:
    test_label = csv.writer(f)
    
    for x in range(len(ypred)):
        a = []
        if ypred[x] < 0.5:
            a.append(0)
            test_label.writerow(a)
        else:
            a.append(1)
            test_label.writerow(a)

即使是单个调参和两两调参，对于我们而言，计算速度还是太慢，我们为此也尝试了Hyperopt方法。通俗的说，我们用的是掷骰子方法，也就是，在一个划定参数区域内随机的掷骰子，哪个参数被掷成几，我们就用这个数来训练模型。最后返回一个掷的最好的参数的结果。这个方法有很大的局限性，一是结果的随机性，二是很容易局部收敛。但是，如果用来粗糙的验证一个特征构造的好坏，也不失为一个不错的方法。

# -*- coding: utf-8 -*-
"""
Created on Fri May 18 14:09:06 2018

@author: jinqiu
"""

import numpy as np
import pandas as pd
from sklearn.preprocessing import MinMaxScaler
import xgboost as xgb
from random import shuffle
from xgboost.sklearn import XGBClassifier
from sklearn.cross_validation import cross_val_score
import pickle
import time
from hyperopt import fmin, tpe, hp,space_eval,rand,Trials,partial,STATUS_OK
import random

data = pd.read_csv(r'D:\next\select_data\new_feature.csv', encoding = "gbk").values
label = pd.read_csv(r'D:\next\select_data\new_label.csv').values
labels = label.reshape((1,-1))
label = labels.tolist()[0]

minmaxscaler = MinMaxScaler()
attrs = minmaxscaler.fit_transform(data)

index = range(0,len(label))
random.shuffle(label)
trainIndex = index[:int(len(label)*0.7)]
print (len(trainIndex))
testIndex = index[int(len(label)*0.7):]
print (len(testIndex))
attr_train = attrs[trainIndex,:]
print (attr_train.shape)
attr_test = attrs[testIndex,:]
print (attr_test.shape)
label_train = labels[:,trainIndex].tolist()[0]
print (len(label_train))
label_test = labels[:,testIndex].tolist()[0]
print (len(label_test))
print (np.mat(label_train).reshape((-1,1)).shape)


def GBM(argsDict):
    max_depth = argsDict["max_depth"] + 5
#    n_estimators = argsDict['n_estimators'] * 5 + 50
    n_estimators = 627
    learning_rate = argsDict["learning_rate"] * 0.02 + 0.05
    subsample = argsDict["subsample"] * 0.1 + 0.7
    min_child_weight = argsDict["min_child_weight"]+1

    print ("max_depth:" + str(max_depth))
    print ("n_estimator:" + str(n_estimators))
    print ("learning_rate:" + str(learning_rate))
    print ("subsample:" + str(subsample))
    print ("min_child_weight:" + str(min_child_weight))

    global attr_train,label_train

    gbm = xgb.XGBClassifier(nthread=6,    #进程数
                            max_depth=max_depth,  #最大深度
                            n_estimators=n_estimators,   #树的数量
                            learning_rate=learning_rate, #学习率
                            subsample=subsample,      #采样数
                            min_child_weight=min_child_weight,   #孩子数

                            max_delta_step = 50,  #50步不降则停止
                            objective="binary:logistic")

    metric = cross_val_score(gbm,attr_train,label_train,cv=3, scoring="f1", n_jobs = -1).mean()
    print (metric)
    return -metric

space = {"max_depth":hp.randint("max_depth",15),
         "n_estimators":hp.quniform("n_estimators",100,1000,1),  #[0,1,2,3,4,5] -> [50,]
        #"learning_rate":hp.quniform("learning_rate",0.01,0.2,0.01),  #[0,1,2,3,4,5] -> 0.05,0.06
         #"subsample":hp.quniform("subsample",0.5,1,0.1),#[0,1,2,3] -> [0.7,0.8,0.9,1.0]
         #"min_child_weight":hp.quniform("min_child_weight",1,6,1), #
         
         #"max_depth":hp.randint("max_depth",15),
        # "n_estimators":hp.randint("n_estimators",10),  #[0,1,2,3,4,5] -> [50,]
         "learning_rate":hp.randint("learning_rate",6),  #[0,1,2,3,4,5] -> 0.05,0.06
         "subsample":hp.randint("subsample",3),#[0,1,2,3] -> [0.7,0.8,0.9,1.0]
         "min_child_weight":hp.randint("min_child_weight",2)

        }
algo = partial(tpe.suggest,n_startup_jobs=1)
best = fmin(GBM,space,algo=algo,max_evals=50) #max_evals表示想要训练的最大模型数量，越大越容易找到最优解

print (best)
print (GBM(best))

最终结果

我们首先把数据进行分类处理。对于那些空缺值的数据，我们直接给label为1（表示异常），没有道理，对于空缺值的处理只能摸奖。在分析训练样本的分布的时候，我们还发现有一些阈值的特征，就是那些特征大于某些值或者小于某些值之后故障风机要明显比正常风机多出很多，这样，我们可以用阈值判断直接给label，剩下不能判断的，我们再用算法进行判断。然后最后时间比较紧，阈值部分的没有做，除去了空缺值之后，其他的全部用算法进行判断。

杨工告诉我们，应该做分析的时候分析“轮毂转速”大于3的数据，因为风机工作才可以检测出来异常，如果风机不工作，是很难判断的。但是因为时间比较紧，对于训练集我们就没有进行这样的操作，于是测试集也没有进行这样的划分。全部都一起塞进算法里了。

总结一下。我们对于特征进行交叉和重要性排序，综合考虑杨工说的重要特征和算法反馈的重要特征排序。最后生成一百多个特征的特征文件用来训练。（训练样本是经过数据清洗过后的样本）
测试集分为两部分，一部分是空缺值，直接标1，另一部分放算法里出结果。

总结

首先最大的一个坑是开始没有做数据清洗的工作，后来发现了之后从新来了一遍。再后来是杨工和我们说应该分析工作风机，拿工作风机来进行训练。如果这样做的话又要推到从来，当时时间已经十分紧张了，心有余而力不足。对于比赛或者说数据分析工作来说，数据的理解是第一位的。否则很可能会做不少无用功。有的时候，受限于专业背景，我们很难充分的理解数据和业务场景，这时候应该向专业人士进行请教，把这些工作都做好之后再进行数据分析要好很多。

其次，提高自己代码和算法的能力。既要懂算法有能撸出一手好代码，这样才能提高效率。我写代码写的太慢，十分制约我的想法实现速度。算法不太懂，不能很好地参与到算法讨论环节。

另外，版本控制十分重要。我们每天的在实现新的想法，文件很多，也很乱。经常出现刚发一个文件，过一段时间就找不到那个文件或者忘了有没有发那个文件。

最后，想说声感谢。在这个比赛过程中得到了刘老师、杨工以及金风科技的工程师们的支持和帮助。特别感谢我的队友们，很多工作都是在他们的帮助之下完成的。leader金指明方向并提供算法，小妹处理数据，牛哥从家里回来的第二天就投入到工作中。金芳虽然去了哈工大，但她也十分关心我们比赛的进程，并常常体意见。我在整个比赛过程中只做了一点点微不足道的贡献，谢谢诸位大佬们。相信我们经过这次比赛以及总结和分享之后将来可以做的更好。