Health_Evaluation/balltreeMSET.PY

# @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()
    '''