# @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_Temp(): matrixD = None matrixL = None healthyResidual = None normalDataBallTree = None def __init__(self): self.model = None def _get_by_id(self, table_name, ids): """Get data from MySQL database by IDs""" df_res = [] engine = create_engine('mysql+pymysql://root:admin123456@106.120.102.238:10336/energy_data_prod') for id in ids: table_name=windcode+'_minute' lastday_df_sql = f"SELECT * FROM {table_name} where id = {id}" df = pd.read_sql(lastday_df_sql, engine) df_res.append(df) return df_ress 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 flag.append(si) return flag if __name__=='__main__': start_time = time.time() myMSET = MSET_Temp() # 1、计算各子系统的健康度(子系统包括：发电机组、传动系统(直驱机组无齿轮箱、无数据)、机舱系统、变流器系统、电网环境、辅助系统（无数据)) # 1.1、发电机组健康度评分 # Title = pd.read_excel(r'/Users/xmia/Documents/code/Temp_Diag/34_QITAIHE.xlsx', header=None, nrows=1, usecols=[12, 13, 14]) # df_D = pd.read_excel(r'/Users/xmia/Documents/code/Temp_Diag/34_QITAIHE.xlsx', # usecols=[0, 12, 13, 14], parse_dates=True) # 读取温度指标：齿轮箱油温12、驱动侧发电机轴承温度13、非驱动侧发电机轴承温度14 df = pd.read_csv('/Users/xmia/Desktop/ZN/华电中光/471_QTH1125/2024/#34.csv', usecols=[ 'wind_turbine_number', 'time_stamp', 'main_bearing_temperature', 'gearbox_oil_temperature', 'generatordrive_end_bearing_temperature', 'generatornon_drive_end_bearing_temperature' ], parse_dates=['time_stamp']) df = df[df['wind_turbine_number'] == 'WOG01312'] df = df.drop(columns=['wind_turbine_number']) df_D = df[df['time_stamp'] > '2024-11-15'] cols = df_D.columns df_D[cols] = df_D[cols].apply(pd.to_numeric, errors='coerce') df_D = df_D.dropna() # 删除空行和非数字项 x_date = df_D.iloc[len(df_D) // 2 + 1:, 0] #获取时间序列 x_date = pd.to_datetime(x_date) df_Temp = df_D.iloc[:, 1:] #获取除时间列以外的数据 df_Temp_values = df_Temp.values df_Temp_values = np.array(df_Temp_values) [m, n] = np.shape(df_Temp_values) # Residual = [] flag_Spart_data = [] for i in range(0, n): df_Temp_i = df_Temp_values[:, i] trainDataSet_data = df_Temp_i[0:len(df_Temp_i) // 2] testDataSet_data = df_Temp_i[len(df_Temp_i) // 2 + 1:] trainDataSet_data = np.reshape(trainDataSet_data, (len(trainDataSet_data), 1)) testDataSet_data = np.reshape(testDataSet_data, (len(testDataSet_data), 1)) myMSET.genDLMatrix(trainDataSet_data, 60, 5) # stateResidual = self.calcStateResidual(testDataSet_data) flag_data = myMSET.calcSPRT(testDataSet_data, np.array(1), decisionGroup=1) # Residual.append(stateResidual) flag_Spart_data.append(flag_data) flag_Spart_data = np.array(flag_Spart_data) flag_Spart_data = flag_Spart_data.T Temp1 = flag_Spart_data[:,0] Temp2 = flag_Spart_data[:,1] Temp3 = flag_Spart_data[:,2] Temp1_lable = "gearbox_oil_temperature" Temp2_lable = "generatordrive_end_bearing_temperature" Temp3_lable = "generatornon_drive_end_bearing_temperature" print(x_date) # alarmedFlag = np.array([[i, Temp1[i]] for i, x in enumerate(Temp1) if x > 0.99]) # Flag中选出大于0.99的点 plt.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文标签 plt.close('all') plt.subplot(311) plt.plot(x_date, Temp1, 'b-o', label=Temp1_lable) plt.ylabel(Temp1_lable) plt.xlabel('时间') # plt.scatter(alarmedFlag1[:, 0], alarmedFlag1[:, 2], marker='x', c='r') plt.subplot(312) plt.plot(x_date, Temp2) plt.ylabel(Temp2_lable) plt.xlabel('时间') #plt.scatter(alarmedFlag[:, 0], alarmedFlag[:, 1], marker='x', c='r') plt.subplot(313) plt.plot(x_date, Temp3) plt.ylabel(Temp3_lable) plt.xlabel('时间') #plt.scatter(alarmedFlag[:, 0], alarmedFlag[:, 1], marker='x', c='r') plt.show() print(flag_Spart_data) print(np.shape(flag_Spart_data)) 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() '''