|
|
@@ -0,0 +1,418 @@
|
|
|
+# @Time : 2019-6-19
|
|
|
+# @Author : wilbur.cheng@SKF
|
|
|
+# @FileName: baseMSET.py
|
|
|
+# @Software: a package in skf anomaly detection library to realize multivariate state estimation technique
|
|
|
+# this package has included all the functions used in MSET, similarity calculation, state memory matrix D,
|
|
|
+# normal state matrix L, residual calculation, and the basic SPRT (Sequential Probability Ratio Test)
|
|
|
+# function for binary hypothesis testing.
|
|
|
+#
|
|
|
+
|
|
|
+import numpy as np
|
|
|
+import pandas as pd
|
|
|
+import matplotlib.pyplot as plt
|
|
|
+from sklearn.neighbors import BallTree
|
|
|
+import openpyxl
|
|
|
+import time
|
|
|
+#from scikit-learn.neighbors import BallTree
|
|
|
+
|
|
|
+class MSET():
|
|
|
+
|
|
|
+ matrixD = None
|
|
|
+ matrixL = None
|
|
|
+ healthyResidual = None
|
|
|
+
|
|
|
+ normalDataBallTree = None
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+ def __init__(self):
|
|
|
+ self.model = None
|
|
|
+
|
|
|
+
|
|
|
+ def calcSimilarity(self, x, y, m = 'euc'):
|
|
|
+ """
|
|
|
+ Calcualte the similartity of two feature vector
|
|
|
+ :param x: one of input feature list
|
|
|
+ :param y: one of input feature list
|
|
|
+ :param m: method of the similarity calculation method, default is the Eucilidean distance named 'euc';
|
|
|
+ a city block distance function is used when m set to "cbd'.
|
|
|
+ :return: the two feature similarity, float, range (0,1)
|
|
|
+ """
|
|
|
+
|
|
|
+ if (len(x) != len(y)):
|
|
|
+ return 0
|
|
|
+
|
|
|
+ if (m == 'cbd'):
|
|
|
+ dSimilarity = [1/(1+np.abs(p-q)) for p,q in zip(x,y)]
|
|
|
+ dSimilarity = np.sum(dSimilarity)/len(x)
|
|
|
+ else:
|
|
|
+ dSimilarity = [np.power(p-q,2) for p, q in zip(x, y)]
|
|
|
+ dSimilarity = 1/(1+np.sqrt(np.sum(dSimilarity)))
|
|
|
+
|
|
|
+ return dSimilarity
|
|
|
+
|
|
|
+
|
|
|
+ def genDLMatrix(self, trainDataset, dataSize4D=100, dataSize4L=50):
|
|
|
+ """
|
|
|
+ automatically generate the D and L matrix from training data set, assuming the training data set is all normal
|
|
|
+ state data.
|
|
|
+ :param trainDataset: 2D array, [count of data， length of feature]
|
|
|
+ :param dataSize4D: minimized count of state for matrix D
|
|
|
+ :param dataSize4L: minimized count of state for matrix L
|
|
|
+ :return: 0 if successful, -1 if fail
|
|
|
+ """
|
|
|
+ [m,n] = np.shape(trainDataset)
|
|
|
+
|
|
|
+ if m < dataSize4D + dataSize4L:
|
|
|
+ print('training dataset is too small to generate the matrix D and L')
|
|
|
+ return -1
|
|
|
+
|
|
|
+ self.matrixD = []
|
|
|
+ selectIndex4D = []
|
|
|
+ # Step 1: add all the state with minimum or maximum value in each feature into Matrix D
|
|
|
+
|
|
|
+ for i in range(0, n):
|
|
|
+ feature_i = trainDataset[:,i]
|
|
|
+ minIndex = np.argmin(feature_i)
|
|
|
+ maxIndex = np.argmax(feature_i)
|
|
|
+
|
|
|
+ self.matrixD.append(trainDataset[minIndex, :].tolist())
|
|
|
+ selectIndex4D.append(minIndex)
|
|
|
+ self.matrixD.append(trainDataset[maxIndex, :].tolist())
|
|
|
+ selectIndex4D.append(maxIndex)
|
|
|
+
|
|
|
+ # Step 2: iteratively add the state with the maximum average distance to the states in selected matrixD
|
|
|
+ while(1):
|
|
|
+
|
|
|
+ if (len(selectIndex4D) >= dataSize4D):
|
|
|
+ break
|
|
|
+
|
|
|
+ # Get the free state list
|
|
|
+ freeStateList = list(set(np.arange(0, len(trainDataset))) - set(selectIndex4D))
|
|
|
+
|
|
|
+ # Calculate the average dist of each state in free to selected state in matrix D
|
|
|
+ distList = []
|
|
|
+ for i in freeStateList:
|
|
|
+ tmpState = trainDataset[i]
|
|
|
+ tmpDist = [1-self.calcSimilarity(x, tmpState) for x in self.matrixD]
|
|
|
+ distList.append(np.mean(tmpDist))
|
|
|
+
|
|
|
+ # select the free state with largest average distance to states in matrixD, and update matrixD
|
|
|
+ selectId = freeStateList[distList.index(np.max(distList))]
|
|
|
+ self.matrixD.append(trainDataset[selectId, :].tolist())
|
|
|
+ selectIndex4D.append(selectId)
|
|
|
+ #format matrixD from list to array
|
|
|
+ self.matrixD = np.array(self.matrixD)
|
|
|
+ self.normalDataBallTree = BallTree(self.matrixD, leaf_size=4, metric = lambda i,j: 1-self.calcSimilarity(i,j))
|
|
|
+
|
|
|
+
|
|
|
+ # Step 3. select remaining state for matrix L
|
|
|
+ #index4L = list(set(np.arange(0, len(trainDataset))) - set(selectIndex4D))
|
|
|
+ #self.matrixL = trainDataset[index4L, :]
|
|
|
+ # consider the limited dataset, using all the train data to matrix L
|
|
|
+ self.matrixL = trainDataset
|
|
|
+
|
|
|
+ # Calculate the healthy Residual from matrix L
|
|
|
+ lamdaRatio = 1e-3
|
|
|
+ [m, n] = np.shape(self.matrixD)
|
|
|
+ self.DDSimilarity = np.array([[1-self.calcSimilarity(x,y) for x in self.matrixD] for y in self.matrixD] + lamdaRatio*np.eye(n))
|
|
|
+ self.DDSimilarity = 1/self.DDSimilarity
|
|
|
+
|
|
|
+ #self.healthyResidual = self.calcResidual(self.matrixL)
|
|
|
+ self.healthyResidual = self.calcResidualByLocallyWeightedLR(self.matrixL)
|
|
|
+
|
|
|
+
|
|
|
+ return 0
|
|
|
+
|
|
|
+ def calcResidualByLocallyWeightedLR(self, newStates):
|
|
|
+ """
|
|
|
+ find the K-nearest neighbors for each input state, then calculate the estimated state by locally weighted average.
|
|
|
+ :param newStates: input states list
|
|
|
+ :return: residual R_x
|
|
|
+ """
|
|
|
+ [m,n] = np.shape(newStates)
|
|
|
+ est_X = []
|
|
|
+ # print(newStates)
|
|
|
+ for x in newStates:
|
|
|
+ (dist, iList) = self.normalDataBallTree.query([x], 20, return_distance=True)
|
|
|
+ weight = 1/(dist[0]+1e-1)
|
|
|
+ weight = weight/sum(weight)
|
|
|
+ eState = np.sum([w*self.matrixD[i] for w,i in zip(weight, iList[0])])
|
|
|
+ est_X.append(eState)
|
|
|
+
|
|
|
+ est_X = np.reshape(est_X, (len(est_X),1))
|
|
|
+ # print(est_X)
|
|
|
+ # print(newStates)
|
|
|
+ return est_X - newStates
|
|
|
+
|
|
|
+ def calcStateResidual(self, newsStates):
|
|
|
+ stateResidual = self.calcResidualByLocallyWeightedLR(newsStates)
|
|
|
+ return stateResidual
|
|
|
+
|
|
|
+ def calcSPRT(self, newsStates, feature_weight, alpha=0.1, beta=0.1, decisionGroup=5):
|
|
|
+ """
|
|
|
+ anomaly detection based Wald's SPRT algorithm, refer to A.Wald, Sequential Analysis,Wiley, New York, NY, USA, 1947
|
|
|
+ :param newsStates: input states list
|
|
|
+ :param feature_weight: the important weight for each feature, Normalized and Nonnegative
|
|
|
+ :param alpha: prescribed false alarm rate， 0 < alpha < 1
|
|
|
+ :param beta: prescribed miss alarm rate, 0 < beta < 1
|
|
|
+ :param decisionGroup: length of the test sample when the decision is done
|
|
|
+
|
|
|
+ :return: anomaly flag for each group of states, 1:anomaly, -1:normal, (-1:1): unable to decision
|
|
|
+ """
|
|
|
+
|
|
|
+ # Step 1. transfer the raw residual vector to dimension reduced residual using feature weight
|
|
|
+ #stateResidual = self.calcResidual(newsStates)
|
|
|
+ stateResidual = self.calcResidualByLocallyWeightedLR(newsStates)
|
|
|
+ # print(stateResidual)
|
|
|
+ weightedStateResidual = [np.dot(x, feature_weight) for x in stateResidual]
|
|
|
+ # print(weightedStateResidual)
|
|
|
+ weightedHealthyResidual = [np.dot(x, feature_weight) for x in self.healthyResidual]
|
|
|
+
|
|
|
+ '''
|
|
|
+ plt.subplot(211)
|
|
|
+ plt.plot(weightedHealthyResidual)
|
|
|
+ plt.xlabel('Sample index')
|
|
|
+ plt.ylabel('Healthy State Residual')
|
|
|
+ plt.subplot(212)
|
|
|
+ plt.plot(weightedStateResidual)
|
|
|
+ plt.xlabel('Sample index')
|
|
|
+ plt.ylabel('All State Residual')
|
|
|
+ plt.show()
|
|
|
+ '''
|
|
|
+ # Calculate the distribution of health state residual
|
|
|
+ mu0 = np.mean(weightedHealthyResidual)
|
|
|
+ sigma0 = np.std(weightedHealthyResidual)
|
|
|
+
|
|
|
+
|
|
|
+ #sigma1 = np.std(weightedStateResidual)
|
|
|
+
|
|
|
+ lowThres = np.log(beta/(1-alpha))
|
|
|
+ highThres = np.log((1-beta)/alpha)
|
|
|
+
|
|
|
+ flag = []
|
|
|
+ for i in range(0, len(newsStates)-decisionGroup+1):
|
|
|
+ # For each group to calculate the new state residual distribution
|
|
|
+ # Then check the hypothesis testing results
|
|
|
+ mu1 = np.mean(weightedStateResidual[i:i+decisionGroup])
|
|
|
+ si = np.sum(weightedStateResidual[i:i+decisionGroup])*(mu1-mu0)/sigma0**2 - decisionGroup*(mu1**2-mu0**2)/(2*sigma0**2)
|
|
|
+
|
|
|
+ if (si > highThres):
|
|
|
+ si = highThres
|
|
|
+ if (si < lowThres):
|
|
|
+ si = lowThres
|
|
|
+
|
|
|
+ if (si > 0):
|
|
|
+ si = si/highThres
|
|
|
+ if (si < 0):
|
|
|
+ si = si/lowThres
|
|
|
+ si = 100-si*100
|
|
|
+ flag.append(si)
|
|
|
+
|
|
|
+ return flag
|
|
|
+
|
|
|
+ def CRITIC_prepare(self, data, flag=1): # 计算权重前的数据预处理
|
|
|
+ '''
|
|
|
+
|
|
|
+ :param data: 输入数据，类型是DataFrame
|
|
|
+ :param flag: flag=0特征正向归一化，flag=1特征负向归一化
|
|
|
+ :return:返回 DataFrame数据
|
|
|
+ '''
|
|
|
+ data_columns = data.columns.values
|
|
|
+ maxnum = np.max(data, axis=0)
|
|
|
+ minnum = np.min(data, axis=0)
|
|
|
+ if flag == 0: # 正向指标归一化计算
|
|
|
+ Y = (data - minnum) * 1.0 / (maxnum - minnum)
|
|
|
+ if flag == 1: # 负向指标归一化计算
|
|
|
+ Y = (maxnum - data) / (maxnum - minnum)
|
|
|
+
|
|
|
+ # 对ln0处理
|
|
|
+ Y0 = np.array(Y * 1.0)
|
|
|
+ Y0[np.where(Y0 == 0)] = 0.00001
|
|
|
+ Y0 = pd.DataFrame(Y0, columns=data_columns)
|
|
|
+ return Y0
|
|
|
+
|
|
|
+ def CRITIC(self, data): # CRITIC客观法计算子系统中各测点的权重
|
|
|
+ '''
|
|
|
+
|
|
|
+ :param data: 归一化预处理之后的DataFrame数据
|
|
|
+ :return: 返回权重Series以及按指标排序后的得分项
|
|
|
+ '''
|
|
|
+ n, m = data.shape
|
|
|
+ s = np.std(data, axis=0)
|
|
|
+ r = np.corrcoef(data, rowvar=False)
|
|
|
+ a = np.sum(1 - r, axis=1)
|
|
|
+ c = s * a
|
|
|
+ w = c / np.sum(c)
|
|
|
+ # score = np.round(np.sum(data * w, axis=1), 6)
|
|
|
+ # data['score'] = score
|
|
|
+ # data.sort_values(by=['score'], ascending=False, inplace=True)
|
|
|
+ return w
|
|
|
+
|
|
|
+ def calcHealthy(self, df_data): # 计算子系统的健康度
|
|
|
+ cols = df_data.columns
|
|
|
+ df_data[cols] = df_data[cols].apply(pd.to_numeric, errors='coerce')
|
|
|
+ df_data = df_data.dropna() # 删除空行和非数字项
|
|
|
+
|
|
|
+ W_prepare_data = self.CRITIC_prepare(df_data)
|
|
|
+ Weights_data = self.CRITIC(W_prepare_data)
|
|
|
+
|
|
|
+
|
|
|
+ df_data_values = df_data.values
|
|
|
+ df_data_values = np.array(df_data_values)
|
|
|
+ [m, n] = np.shape(df_data_values)
|
|
|
+ # Residual = []
|
|
|
+ flag_Spart_data = []
|
|
|
+ for i in range(0, n):
|
|
|
+ df_data_i = df_data_values[:, i]
|
|
|
+ trainDataSet_data = df_data_i[0:len(df_data_i) // 2]
|
|
|
+ testDataSet_data = df_data_i[len(df_data_i) // 2 + 1:]
|
|
|
+ trainDataSet_data = np.reshape(trainDataSet_data, (len(trainDataSet_data), 1))
|
|
|
+ testDataSet_data = np.reshape(testDataSet_data, (len(testDataSet_data), 1))
|
|
|
+
|
|
|
+ self.genDLMatrix(trainDataSet_data, 60, 5)
|
|
|
+ # stateResidual = self.calcStateResidual(testDataSet_data)
|
|
|
+ flag_data = self.calcSPRT(testDataSet_data, np.array(1), decisionGroup=1)
|
|
|
+ # Residual.append(stateResidual)
|
|
|
+ flag_Spart_data.append(flag_data)
|
|
|
+ # weights = np.array([1/3, 1/3, 1/3])
|
|
|
+ W_data = np.array(Weights_data)
|
|
|
+ W_data = W_data.reshape(-1, 1)
|
|
|
+ flag_Spart_data = np.array(flag_Spart_data)
|
|
|
+ flag_Spart_data = flag_Spart_data.T
|
|
|
+ Score_data = np.dot(flag_Spart_data, W_data)
|
|
|
+ Health_data = np.mean(Score_data)
|
|
|
+
|
|
|
+ return Weights_data, Health_data
|
|
|
+
|
|
|
+ def ahp(self, matrix): #层次分析法ahp(主观)计算风机各子系统的权重
|
|
|
+ # 计算特征值和特征向量
|
|
|
+ eigenvalue, eigenvector = np.linalg.eig(matrix)
|
|
|
+ # 提取最大特征值对应的特征向量，并归一化
|
|
|
+ max_eigenvalue_index = np.argmax(eigenvalue)
|
|
|
+ max_eigenvalue = eigenvalue[max_eigenvalue_index]
|
|
|
+ weight = eigenvector[:, max_eigenvalue_index]
|
|
|
+ weight = weight / np.sum(weight)
|
|
|
+ weight = weight.real
|
|
|
+ return weight, max_eigenvalue
|
|
|
+
|
|
|
+ def consistency_check(self, matrix, weight): #层次分析法ahp的一致性检验
|
|
|
+ n = len(matrix)
|
|
|
+ # 计算一致性指标 CI
|
|
|
+ lambda_max = np.sum(np.dot(matrix, weight) / weight) / n
|
|
|
+ CI = (lambda_max - n) / (n - 1)
|
|
|
+
|
|
|
+ # 随机一致性指标 RI 的值，这里我们假设矩阵大小不超过10
|
|
|
+ RI_list = [0, 0, 0.58, 0.90, 1.12, 1.24, 1.32, 1.41, 1.45, 1.49, 1.52, 1.54, 1.56, 1.58, 1.59, 1.61, 1.63, 1.64,
|
|
|
+ 1.65]
|
|
|
+ RI = RI_list[n - 1]
|
|
|
+
|
|
|
+ # 计算一致性比例 CR
|
|
|
+ CR = CI / RI
|
|
|
+ return CI, CR
|
|
|
+
|
|
|
+if __name__=='__main__':
|
|
|
+
|
|
|
+ start_time = time.time()
|
|
|
+ myMSET = MSET()
|
|
|
+ # 1、计算各子系统的健康度(子系统包括：发电机组、传动系统(直驱机组无齿轮箱、无数据)、机舱系统、变流器系统、电网环境、辅助系统（无数据))
|
|
|
+ # 1.1、发电机组健康度评分
|
|
|
+ df_Gen = pd.read_excel(r'F26_SEP_ZHANGYAOXIAN.xlsx',
|
|
|
+ usecols=[10, 11, 12, 17]) # 读取发电机组指标：发电机U相温度10、发电机V相温度11、发电机W相温度12、发电机轴承温度17
|
|
|
+ W_Gen,Health_Gen = myMSET.calcHealthy(df_Gen)
|
|
|
+ print(W_Gen)
|
|
|
+ print('发电机组健康度:', Health_Gen)
|
|
|
+
|
|
|
+ # 1.2、机舱系统健康度评分
|
|
|
+ df_Nac = pd.read_excel(r'F26_SEP_ZHANGYAOXIAN.xlsx',
|
|
|
+ usecols=[23,24,41,42]) # 读取机舱系统指标：塔筒前后振动23、塔筒左右振动24、机舱位置41、机舱温度42
|
|
|
+ W_Nac, Health_Nac = myMSET.calcHealthy(df_Nac)
|
|
|
+ print(W_Nac)
|
|
|
+ print('机舱系统健康度:', Health_Nac)
|
|
|
+
|
|
|
+ # 1.3、变流器系统健康度评分
|
|
|
+ df_Con = pd.read_excel(r'F26_SEP_ZHANGYAOXIAN.xlsx',
|
|
|
+ usecols=[18,19]) # 读取变流器系统指标：变流器有功功率19、变流器冷却介质温度18
|
|
|
+ W_Con, Health_Con = myMSET.calcHealthy(df_Con)
|
|
|
+ print(W_Con)
|
|
|
+ print('变流器系统健康度:', Health_Con)
|
|
|
+
|
|
|
+ # 1.4、电网环境健康度评分
|
|
|
+ df_Grid = pd.read_excel(r'F26_SEP_ZHANGYAOXIAN.xlsx',
|
|
|
+ usecols=[32,38,64,65,66]) # 读取电网环境指标：有功功率38、无功功率32、电网A相电流64、电网B相电流65、电网C相电流66
|
|
|
+ W_Grid, Health_Grid = myMSET.calcHealthy(df_Grid)
|
|
|
+ print(W_Grid)
|
|
|
+ print('电网环境健康度:', Health_Grid)
|
|
|
+
|
|
|
+ # 2、计算各子系统的权重
|
|
|
+ # 输入判断矩阵（发电机组、机舱系统、变流器系统、电网环境）
|
|
|
+ matrix_subsys = np.array([
|
|
|
+ [1, 2, 3, 4],
|
|
|
+ [1/2, 1, 2, 3],
|
|
|
+ [1/3, 1/2, 1, 2],
|
|
|
+ [1/4, 1/3, 1/2, 1],
|
|
|
+ ])
|
|
|
+ # 计算判断矩阵的权重和最大特征值
|
|
|
+ weight_subsys, max_eigenvalue_subsys = myMSET.ahp(matrix_subsys)
|
|
|
+ print(weight_subsys)
|
|
|
+
|
|
|
+ # 检查一致性
|
|
|
+ CI1, CR1 = myMSET.consistency_check(matrix_subsys, weight_subsys)
|
|
|
+
|
|
|
+ # 3、计算整机的健康度
|
|
|
+ Score_subsys = np.array([Health_Gen, Health_Nac, Health_Con, Health_Grid])
|
|
|
+ weight_subsys = weight_subsys.reshape(-1, 1)
|
|
|
+ Score_WT = np.dot(Score_subsys, weight_subsys)
|
|
|
+ print('整机健康度:', Score_WT)
|
|
|
+
|
|
|
+ end_time = time.time()
|
|
|
+ execution_time = end_time - start_time
|
|
|
+ print(f"Execution time: {execution_time} seconds")
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+#
|
|
|
+
|
|
|
+ '''
|
|
|
+ f = open("speed_vib.txt", "r")
|
|
|
+ data1 = f.read()
|
|
|
+ f.close()
|
|
|
+ data1 = data1.split('\n')
|
|
|
+ rpm = [(row.split('\t')[0]).strip() for row in data1]
|
|
|
+ vib = [(row.split('\t')[-1]).strip() for row in data1]
|
|
|
+
|
|
|
+ # print(rpm)
|
|
|
+ rpm = np.array(rpm).astype(np.float64)
|
|
|
+ vib = np.array(vib).astype(np.float64)
|
|
|
+
|
|
|
+ #vib = [(x-np.mean(vib))/np.std(vib) for x in vib]
|
|
|
+ #print(vib)
|
|
|
+ trainDataSet = [vib[i] for i in range(0,100) if vib[i] < 5]
|
|
|
+ trainDataSet = np.reshape(trainDataSet,(len(trainDataSet),1))
|
|
|
+ testDataSet = np.reshape(vib, (len(vib),1))
|
|
|
+ '''
|
|
|
+
|
|
|
+
|
|
|
+ # Title = pd.read_csv(r'F1710001001.csv', header=None, nrows=1, usecols=[36,37,38], encoding='gbk')
|
|
|
+
|
|
|
+ '''
|
|
|
+ alarmedFlag = np.array([[i, flag[i]] for i, x in enumerate(flag) if x > 0.99]) # Flag中选出大于0.99的点
|
|
|
+ plt.close('all')
|
|
|
+ plt.subplot(211)
|
|
|
+ plt.plot(testDataSet)
|
|
|
+ plt.ylabel('Vibration')
|
|
|
+ plt.xlabel('Sample index')
|
|
|
+ plt.subplot(212)
|
|
|
+ plt.plot(flag)
|
|
|
+ plt.ylabel('SPRT results')
|
|
|
+ plt.xlabel('Sample index')
|
|
|
+ plt.scatter(alarmedFlag[:,0], alarmedFlag[:,1], marker='x',c='r')
|
|
|
+
|
|
|
+ plt.show()
|
|
|
+ '''
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+
|
lcb
1 maand geleden