首先介紹一下比賽背景。這個比賽是中電投的一個關于風機開裂故障分析的預警的比賽。
02 資料介紹
每個樣本10分鐘之内大概采集到了450條資料，一共75個特征，也就是差不多75*450個資訊。最後三個特征完全沒有資料，是以一開始做的時候，我們就把最後三個特征進行删除，實際上我們是對72個特征進行的資料分析。最開始，用的是seaborn畫的正常風機和不正常風機的頻率分布圖，比如說對于輪毂轉速這個特征：
1import seaborn as snsimport pandas as pd
2data_file = r"D:\fan_fault\feature1.csv"
3pre_process = pd.read_csv(data_file, encoding = "gbk")
4
5pre_process = pre_process.fillna(0)
6feature1_plot = pre_process["normal(0)"]
7
8feature2_plot2 = pre_process["fault(1)"]
9sns.kdeplot(feature1_plot, shade = True)
10sns.kdeplot(feature2_plot2, shade = True)
大部分特征都是這樣，沒有很好地區分度。正是因為這樣，也一直沒有嘗試出來非常好的模型。後來我們嘗試用MATLAB畫圖，每個特征出兩個圖
看起來要比seaborn好一些（後兩個圖和第一個不是一個特征）。我們在做資料分析這一塊一個很大的問題是在于隻去查了各個特征的實體含義，做了頻率和頻數分布圖，看看是否有沒有好的特征，然後就直接進入了下一步。忘了考慮是否可能會出現因為采集的問題而導緻的異常值和空缺值的問題。這一點導緻後面我們很多工作推到從來。
03 資料分析
我們從統計上來看，并沒有找到很好區分度的特征，然後就考慮從實體上來找。在老師的建議下，我們嘗試了有輪毂轉速，風速為6.5m/s時，y方向振動值的特征：
依舊沒有很好的區分度，對于其他風速嘗試，也是如此。
我們讨論門檻值、記0等方式構造新特征。就在說道記0這個新特征構造辦法的時候，突然發現，大氣壓力這個特征，居然有0的情況。根據實體學的知識來講，風機的大氣壓力是不可能為0的。然後我們才想起來，沒有對資料的異常值進行處理。删除了有8萬多條整行全為0的資料，導緻某些檔案為空，也就是這個風機沒有資料資訊。當然，也有某些風機是某幾行為0。除了删除空缺值，我們還對其他明顯是異常的資料進行了一些資料清洗工作。因為之前我們對于資料特征數統計分析是根據未清洗的資料做的分析，是以分析的可靠性也有點問題，後面我們在做工作的時候，有了一些不必要的麻煩。我們也做了一些相關性分析的工作，大部分特征相關性十分的高。幾十個特征兩兩組合然後進行相關性分析，會有數千個結果，相關性分析沒有辦法進行下去。後來，我們就沒有考慮相關性的事情。
04 特征工程
我們最開始嘗試對前72個特征構造均值，作為基準嘗試看看效果如何。
1import os
2import pandas as pd
3import numpy as np
4import csv
5
6label_file = r"C:\fan_fault\train\trainX"
7train_mean = r"D:\fan_fault\train_mean_new.csv"
8
9with open(train_mean, "a", newline = '', encoding = "utf-8") as f:
10 train_mean = csv.writer(f)
11
12 for x in range(1, 48340):
13 fan_file = os.path.join(label_file, str(x) + ".csv")
14 print("程式運作進度為", x/48340) #用該語句檢視工作進度狀态
15
16 with open(fan_file, encoding='utf-8') as f:
17 feature_read = pd.read_csv(f)
18 #周遊打開檔案的每一個特征（72），求取均值
19 # a用來臨時存放計算好的特征均值，外加一個label
20
21 a = []
22 for i in range(72):
23 mean_num = feature_read.iloc[:, i]
24 mean_num = np.array(mean_num).mean()
25 #生成每個特征所有資料對應的均值
26 a.append(mean_num)
27
28 train_mean.writerow(a)
也包括絕對值差分累計、差分均值、差分方差，用随機森林進行調參。
1# -*- coding: utf-8 -*-"""
2
4import pandas as pd
5from sklearn.preprocessing import MinMaxScaler
6from sklearn.ensemble import RandomForestClassifier
7from sklearn.model_selection import cross_val_scorefrom sklearn
8import metrics
9from sklearn.model_selection import GridSearchCV
10
11#資料導入、檢查空缺值
12data = pd.read_csv(r'D:\next\8_19\train_data.csv',encoding = "gbk")
13label = pd.read_csv(r"D:\next\8_19\train_label.csv")
14data.info()
15data.notnull().sum(axis=0)/data.shape[0]
16train = data.iloc[:,:-1]
17label = label.iloc[:,-1]
18
19#資料标準化
20scaler = MinMaxScaler()
21train = scaler.fit(train).transform(train)
22
23#單個分類器
24clf = RandomForestClassifier(random_state=14)
25f1 = cross_val_score(clf, train, label, scoring='f1')
26print("f1:{0:.1f}%".format(np.mean(f1)*100))
28#調參
29parameter_space = {
30 'n_estimators':range(10,200,10),
31 'max_depth':range(1,10),
32 'min_samples_split':range(2,10),
33 }
34clf = RandomForestClassifier(random_state=14)
35grid = GridSearchCV(clf,parameter_space,scoring='f1', n_jobs = 6)
36grid.fit(train,label)
37print("f1:(0:.1f)%".format(grid.best_score_*100))
38print(grid.best_estimator_)
39
40#調參後的分類器
41new_clf = RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',
42 max_depth=7, max_features='auto', max_leaf_nodes=None,
43 min_impurity_decrease=0.0, min_impurity_split=None,
44 min_samples_leaf=1, min_samples_split=7,
45 min_weight_fraction_leaf=0.0, n_estimators=10, n_jobs=1,
46 oob_score=False, random_state=14, verbose=0,warm_start=False)
47print("f1:{0:.1f}%".format(np.mean(f1)*100))
測試集輸出預測結果：
1#資料标準化
2scaler = MinMaxScaler()
3train = scaler.fit(train).transform(train)
4test = scaler.fit(test).transform(test)
6#訓練分類器
7clf = RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',
8 max_depth=8, max_features='auto', max_leaf_nodes=None,
9 min_impurity_decrease=0.0, min_impurity_split=None,
10 min_samples_leaf=1, min_samples_split=5,
11 min_weight_fraction_leaf=0.0, n_estimators=20, n_jobs=5,
12 oob_score=False, random_state=14, verbose=0, warm_start=False)
13clf = clf.fit(train, label)
14#預測結果
15pre = clf.predict(test)
16
17#測試結果檔案寫入
18import csv
19
20label = r"D:/fan_fault/label.csv"
21
22with open(label, "w", newline = '', encoding = "utf-8") as f:
23 label = csv.writer(f)
24
25 for x in range(len(pre)):
26 label.writerow(pre[x:x+1])
測試效果來看，并沒有取得十分理想的效果。
之後想起來沒有考慮資料清洗，之前的工作全部推倒從來。我們在此期間，也調查整理關于這72個特征和風機葉片開裂故障的原因，發現從文獻上沒有找到和葉片開裂故障有很好相關性的特征，我們對于特征進行一些改造，也沒有很好的區分效果。後期我們咨詢過金風科技的工程師和阜特科技的工程師，他們說這些特征中也沒有十分強相關性的特征。我們就考慮到特征之間兩兩交叉相乘看看效果。
1# 交叉特征有2844列，分别是自身的平方和互相交叉，最後求均值方差，最後三個特征不單獨再生成交叉特征
3import os
5import numpy as np
6import csv
7from sklearn.preprocessing import PolynomialFeatures
9label_file = r"F:\User\Xinyuan Huang\train_labels.csv"
10fan_folder = r"F:\User\Xinyuan Huang"
11read_label = pd.read_csv(label_file)
12
13cross_var = r"F:\User\Xinyuan Huang\CaiJi\Feature_crosses\cross_var.csv"
14
15with open(cross_var, "a", newline = '', encoding = "utf-8") as f:
16 cross_var = csv.writer(f)
17
18 # 該for循環用于定位要打開的檔案
19 for x in range(len(read_label)-1):
20 column1 = str(read_label["f_id"][x:x+1])
21 #周遊DataFrame第一列的f_id标簽下面的每一個數
22 column2 = str(read_label["file_name"][x:x+1])
23 #周遊DataFrame第二列的file_name标簽下面的每一個數
24 column3 = str(read_label["ret"][x:x+1])
25 #周遊DataFrame第三列的ret标簽下面的每一個數
26
27 f_id = column1.split()[1]
28 #第一行的檔案所對應的f_id進行切片操作，擷取對應的數字
29 # 對f_id進行補0操作
30 f_id = f_id.zfill(3)
31 # 比如2補成002，是以這裡寫3
32 file_name = column2.split()[1]
33 #第一行的檔案所對應的file_name
34 label = column3.split()[1]
35 #第一行檔案所對應的ret
36
37 fan_file = os.path.join(fan_folder, "train", f_id, file_name)
38 print("程式運作進度為", x/(len(read_label)-1))
39 #用該語句檢視工作進度狀态
40
41 # 打開相應的fan_file檔案進行讀取操作
42 with open(fan_file, encoding='utf-8') as f:
43 dataset = pd.read_csv(f)
44 #資料集名稱為dataset
45 poly = PolynomialFeatures(degree=2, include_bias=False,interaction_only=False)
46 X_ploly = poly.fit_transform(dataset)
47 data_ploly = pd.DataFrame(X_ploly, columns=poly.get_feature_names())
48
49 new_data = data_ploly.ix[:,75:-6]
50
51 #ploly_mean,ploly_var為交叉特征均值方差
52 ploly_mean = np.mean(new_data)
53 ploly_var = np.var(ploly_mean)
54
55 ploly_var = list(ploly_var)
56 ploly_var.append(label)
57
58 cross_var.writerow(ploly_var)
交叉相乘之後的檔案有數千個特征，生成的檔案有将近2G大小。考慮到伺服器性能不高，計算曠日持久。不對特征進行篩選，直接進行交叉之後跑算法這條路被我們放棄了。
後來阜特科技的楊工幫我們篩選了一些比較重要的特征，我們在此基礎之上進行了一些特征交叉和重要性排序的操作，特征縮小到了幾百個（包含交叉、均值、方差等，經過重要性排序），然後用它來跑的模型。
特征裡面有一些特征是離散特征，對于這些特征我們進行單獨處理，進行離散化。比如說偏航要求值總共有3個值分别是1,2,3。我們對其進行離散化處理，一個特征就變成了三個特征。每個樣本統計出現這三個特征的頻率。
6label_file = r"E:\8_19\testX_csv"
8train_mean = r"E:\8_19\disperse\discrete56.csv"
9
10with open(train_mean, "a", newline = '', encoding = "utf-8") as f:
11 train_mean = csv.writer(f)
13 for x in range(1, 451):
14 fan_file = os.path.join(label_file, str(x) + ".csv")
15# print("程式運作進度為", x/451) #用該語句檢視工作進度狀态
17 with open(fan_file, encoding='utf-8') as f:
18 feature_read = pd.read_csv(f, header = None)
20 num1 = 0
21 num2 = 0
22 num3 = 0
23
24 a = []
25
26 for x in range(len(feature_read)):
27 if feature_read[55][x] == 0:
28 num1 = num1+1
29 if feature_read[55][x] == 1:
30 num2 = num2+1
31 if feature_read[55][x] == 2:
32 num3 = num3+1
33
34 num1 = num1/len(feature_read)
35 num2 = num2/len(feature_read)
36 num3 = num3/len(feature_read)
37
38 a.append(num1)
39 a.append(num2)
40 a.append(num3)
41
42 train_mean.writerow(a)
05 算法
我們最後主要用的算法是Xgboost，期間也嘗試過LightGBM，因為算力不夠的原因，沒有辦法嘗試一些算法（包括楊工說的SVM以及深度學習的想法），最後主要用Xgboost進行調參，直接一起調參的話算力不夠，我們是單個調參以及兩兩調參組合的形式進行參數調整。
1from xgboost import XGBClassifier
2import xgboost as xgb
3
4import pandas as pd
6
7from sklearn.model_selection import GridSearchCV
8from sklearn.model_selection import StratifiedKFold
10from sklearn.metrics import log_loss
11from sklearn.preprocessing import MinMaxScaler
13
14#資料導入、檢查空缺值
15data = pd.read_csv(r'D:\next\8_19\train_data.csv',encoding = "gbk")
16label = pd.read_csv(r"D:\next\8_19\train_label.csv")
17test = pd.read_csv(r"D:\next\8_19\test_data.csv", encoding = "gbk")
18train = data.iloc[:,:-1]
19label = label.iloc[:,-1]
21X_train = train
22y_train = label
24#資料标準化
25scaler = MinMaxScaler()
26train = scaler.fit(train).transform(train)
27test = scaler.fit(test).transform(test)
28
29#交叉驗證
30kfold = StratifiedKFold(n_splits=5, shuffle=True, random_state=3)
31
32param_test1 = {
33 'max_depth':list(range(3,10,1)),
34 'min_child_weight':list(range(1,6,1))
35}
36gsearch1 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=400, max_depth=5,
37 min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,
38 objective= 'binary:logistic', nthread=4, scale_pos_weight=1, seed=27),
39 param_grid = param_test1, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
40gsearch1.fit(X_train,y_train)
41gsearch1.grid_scores_, gsearch1.best_params_, gsearch1.best_score_
調完這個參數之後，把最好的輸出結果拿出來放在下一個參數調優裡進行調整：
1aram_test1 = {
2 'learning_rate':[i/100.0 for i in range(6,14,2)]
3}
4gsearch1 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=400, max_depth=6,
5 min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,
6 objective= 'binary:logistic', nthread=6, scale_pos_weight=1, seed=27),
7 param_grid = param_test1, scoring='roc_auc',n_jobs=-1,iid=False, cv=5)
8gsearch1.fit(X_train,y_train)
9gsearch1.grid_scores_, gsearch1.best_params_, gsearch1.best_score_
11param_test1 = {
12 'subsample':[0.8, 0.9]
13}
14gsearch1 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=310, max_depth=6,
15 min_child_weight=1, gamma=0, subsample=0.9, colsample_bytree=0.8,
16 objective= 'binary:logistic', nthread=6, scale_pos_weight=1, seed=27),
17 param_grid = param_test1, scoring='roc_auc',n_jobs=-1,iid=False, cv=5)
18gsearch1.fit(X_train,y_train)
19gsearch1.grid_scores_, gsearch1.best_params_, gsearch1.best_score_
還可以調整好幾個參數，最後結果輸出：
1import xgboost as xgb
2dtrain=xgb.DMatrix(X_train ,label=y_train)
3dtest=xgb.DMatrix(test)
5params={
6 'objective': 'binary:logistic',
7 'max_depth':6,
8 'subsample':0.8,
9 'colsample_bytree':0.8,
10 'min_child_weight':1,
11 'seed':27,
12 'nthread':6,
13 'learning_rate':0.1,
14 'n_estimators':292,
15 'gamma':0,
16 'scale_pos_weight':1}
18watchlist = [(dtrain,'train')]
20bst=xgb.train(params,dtrain,num_boost_round=100,evals=watchlist)
22ypred=bst.predict(dtest)
24import csv
26test_label = r"D:\next\8_20\test_label_new.csv"
27with open(test_label, "a", newline = '', encoding = "utf-8") as f:
28 test_label = csv.writer(f)
29
30 for x in range(len(ypred)):
31 a = []
32 if ypred[x] < 0.5:
33 a.append(0)
34 test_label.writerow(a)
35 else:
36 a.append(1)
37 test_label.writerow(a)
即使是單個調參和兩兩調參，對于我們而言，計算速度還是太慢，我們為此也嘗試了Hyperopt方法。通俗的說，我們用的是擲骰子方法，也就是，在一個劃定參數區域内随機的擲骰子，哪個參數被擲成幾，我們就用這個數來訓練模型。最後傳回一個擲的最好的參數的結果。這個方法有很大的局限性，一是結果的随機性，二是很容易局部收斂。但是，如果用來粗糙的驗證一個特征構造的好壞，也不失為一個不錯的方法。
1# -*- coding: utf-8 -*-
2"""
3Created on Fri May 18 14:09:06 2018
6"""
8import numpy as np
9import pandas as pd
10from sklearn.preprocessing import MinMaxScaler
11import xgboost as xgb
12from random import shuffle
13from xgboost.sklearn import XGBClassifier
14from sklearn.cross_validation import cross_val_score
15import pickle
16import time
17from hyperopt import fmin, tpe, hp,space_eval,rand,Trials,partial,STATUS_OK
18import random
20data = pd.read_csv(r'D:\next\select_data\new_feature.csv', encoding = "gbk").values
21label = pd.read_csv(r'D:\next\select_data\new_label.csv').values
22labels = label.reshape((1,-1))
23label = labels.tolist()[0]
25minmaxscaler = MinMaxScaler()
26attrs = minmaxscaler.fit_transform(data)
28index = range(0,len(label))
29random.shuffle(label)
30trainIndex = index[:int(len(label)*0.7)]
31print (len(trainIndex))
32testIndex = index[int(len(label)*0.7):]
33print (len(testIndex))
34attr_train = attrs[trainIndex,:]
35print (attr_train.shape)
36attr_test = attrs[testIndex,:]
37print (attr_test.shape)
38label_train = labels[:,trainIndex].tolist()[0]
39print (len(label_train))
40label_test = labels[:,testIndex].tolist()[0]
41print (len(label_test))
42print (np.mat(label_train).reshape((-1,1)).shape)
43
44
45def GBM(argsDict):
46 max_depth = argsDict["max_depth"] + 5
47# n_estimators = argsDict['n_estimators'] * 5 + 50
48 n_estimators = 627
49 learning_rate = argsDict["learning_rate"] * 0.02 + 0.05
50 subsample = argsDict["subsample"] * 0.1 + 0.7
51 min_child_weight = argsDict["min_child_weight"]+1
52
53 print ("max_depth:" + str(max_depth))
54 print ("n_estimator:" + str(n_estimators))
55 print ("learning_rate:" + str(learning_rate))
56 print ("subsample:" + str(subsample))
57 print ("min_child_weight:" + str(min_child_weight))
58
59 global attr_train,label_train
60
61 gbm = xgb.XGBClassifier(nthread=6, #程序數
62 max_depth=max_depth, #最大深度
63 n_estimators=n_estimators, #樹的數量
64 learning_rate=learning_rate, #學習率
65 subsample=subsample, #采樣數
66 min_child_weight=min_child_weight, #孩子數
67
68 max_delta_step = 50, #50步不降則停止
69 objective="binary:logistic")
70
71 metric = cross_val_score(gbm,attr_train,label_train,cv=3, scoring="f1", n_jobs = -1).mean()
72 print (metric)
73 return -metric
74
75space = {"max_depth":hp.randint("max_depth",15),
76 "n_estimators":hp.quniform("n_estimators",100,1000,1), #[0,1,2,3,4,5] -> [50,]
77 #"learning_rate":hp.quniform("learning_rate",0.01,0.2,0.01), #[0,1,2,3,4,5] -> 0.05,0.06
78 #"subsample":hp.quniform("subsample",0.5,1,0.1),#[0,1,2,3] -> [0.7,0.8,0.9,1.0]
79 #"min_child_weight":hp.quniform("min_child_weight",1,6,1), #
80
81 #"max_depth":hp.randint("max_depth",15),
82 # "n_estimators":hp.randint("n_estimators",10), #[0,1,2,3,4,5] -> [50,]
83 "learning_rate":hp.randint("learning_rate",6), #[0,1,2,3,4,5] -> 0.05,0.06
84 "subsample":hp.randint("subsample",3),#[0,1,2,3] -> [0.7,0.8,0.9,1.0]
85 "min_child_weight":hp.randint("min_child_weight",2)
86
87 }
88algo = partial(tpe.suggest,n_startup_jobs=1)
89best = fmin(GBM,space,algo=algo,max_evals=50) #max_evals表示想要訓練的最大模型數量，越大越容易找到最優解
90
91print (best)
92print (GBM(best))
06 最終結果
我們首先把資料進行分類處理。對于那些空缺值的資料，我們直接給label為1（表示異常），沒有道理，對于空缺值的處理隻能摸獎。在分析訓練樣本的分布的時候，我們還發現有一些門檻值的特征，就是那些特征大于某些值或者小于某些值之後故障風機要明顯比正常風機多出很多，這樣，我們可以用門檻值判斷直接給label，剩下不能判斷的，我們再用算法進行判斷。然後最後時間比較緊，門檻值部分的沒有做，除去了空缺值之後，其他的全部用算法進行判斷。
楊工告訴我們，應該做分析的時候分析“輪毂轉速”大于3的資料，因為風機工作才可以檢測出來異常，如果風機不工作，是很難判斷的。但是因為時間比較緊，對于訓練集我們就沒有進行這樣的操作，于是測試集也沒有進行這樣的劃分。全部都一起塞進算法裡了。
總結一下。我們對于特征進行交叉和重要性排序，綜合考慮楊工說的重要特征和算法回報的重要特征排序。最後生成一百多個特征的特征檔案用來訓練。（訓練樣本是經過資料清洗過後的樣本）
測試集分為兩部分，一部分是空缺值，直接标1，另一部分放算法裡出結果。
07 總結
首先最大的一個坑是開始沒有做資料清洗的工作，後來發現了之後從新來了一遍。再後來是楊工和我們說應該分析工作風機，拿工作風機來進行訓練。如果這樣做的話又要推到從來，當時時間已經十分緊張了，心有餘而力不足。對于比賽或者說資料分析工作來說，資料的了解是第一位的。否則很可能會做不少無用功。有的時候，受限于專業背景，我們很難充分的了解資料和業務場景，這時候應該向專業人士進行請教，把這些工作都做好之後再進行資料分析要好很多。
其次，提高自己代碼和算法的能力。既要懂算法有能撸出一手好代碼，這樣才能提高效率。我寫代碼寫的太慢，十分制約我的想法實作速度。算法不太懂，不能很好地參與到算法讨論環節。
另外，版本控制十分重要。我們每天的在實作新的想法，檔案很多，也很亂。經常出現剛發一個檔案，過一段時間就找不到那個檔案或者忘了有沒有發那個檔案。
原文釋出時間為：2018-09-15
本文作者：Einstellung
本文來自雲栖社群合作夥伴“
Python愛好者社群”，了解相關資訊可以關注“
”。