model-lab/PCA_Test_offline.py

# -*- coding: utf-8 -*-
"""
Created on Sun Feb 28 10:04:26 2016
PCA source code————最新更新————————————————————————————
@author: liudiwei
"""


import numpy as np
import json
import pymssql
import requests
import jenkspy
import gc
import pyodbc
from recon import Lars, recon_fault_diagnosis_r, recon_fault_diagnosis_r_l, recon_fault_diagnosis_r_c
import config


"""
参数：
	- XMat：传入的是一个numpy的矩阵格式，行表示样本数，列表示特征    
	- k：表示取前k个特征值对应的特征向量
返回值：
	- finalData：参数一指的是返回的低维矩阵，对应于输入参数二
	- reconData：参数二对应的是移动坐标轴后的矩阵
"""


def min_pos(X):
    X[X <= 0] = np.max(X)
    m = np.min(X)
    re = np.where(X == m)
    min_i = re[0]
    min_j = re[1]
    if m < 0:
        m = 0
    return m, min_i, min_j

class MSSQL:
    def __init__(self,host,user,pwd,database):
        self.host = host
        self.user = user
        self.pwd = pwd
        self.db = database

    def __GetConnect(self):
        """
        得到连接信息
        返回: conn.cursor()
        """
        if not self.db:
            raise(NameError,"没有设置数据库信息")
        self.conn = pymssql.connect(host=self.host,user=self.user,password=self.pwd,database=self.db,charset="utf8")
        cur = self.conn.cursor()
        if not cur:
            raise(NameError,"连接数据库失败")
        else:
            return cur

    def ExecQuery(self,sql):
        """
        执行查询语句
        返回的是一个包含tuple的list，list的元素是记录行，tuple的元素是每行记录的字段
        """
        cur = self.__GetConnect()
        cur.execute(sql)
        resList = cur.fetchall()

        #查询完毕后必须关闭连接
        self.conn.close()
        return resList

    def ExecNonQuery(self,sql):
        """
        执行非查询语句

        调用示例：
            cur = self.__GetConnect()
            cur.execute(sql)
            self.conn.commit()
            self.conn.close()
        """
        cur = self.__GetConnect()
        cur.execute(sql)
        self.conn.commit()
        self.conn.close()

def get_model_by_ID(model_id):
    ms = MSSQL(host="120.26.116.243", user="root", pwd="powerSIS#123", database="alert")
    resList = ms.ExecQuery("SELECT Model_info FROM model_cfg where \"model_id\"="+str(model_id))
    #return json.loads(resList[0][0])["para"]
    return json.loads(resList[0][0])


def get_model_by_id(model_id):
    try:
        conn = pyodbc.connect(
            r"DRIVER={ODBC Driver 17 for SQL Server};SERVER=%s;DATABASE=alert;UID=sa;PWD=powerSIS#123" % config._SQL_IP)  # 连接数据库
    except pyodbc.Error:
        conn = pyodbc.connect(
            r"DRIVER={SQL SERVER NATIVE CLIENT 10.0};SERVER=%s;DATABASE=alert;UID=sa;PWD=powerSIS#123" % config._SQL_IP)  # 连接数据库
    cursor = conn.cursor()  # 获得操作的游标
    cursor.execute(f"SELECT Model_info FROM model_cfg where model_id={model_id}")
    res_list = cursor.fetchall()  # 获取查询的结果
    conn.commit()  # 提交执行
    cursor.close()  # 关闭游标
    conn.close()  # 关闭数据库连接
    return json.loads(res_list[0][0])


def get_model_by_id_and_version(model_id, version):
    try:
        conn = pyodbc.connect(
            r"DRIVER={ODBC Driver 17 for SQL Server};SERVER=%s;DATABASE=alert;UID=sa;PWD=powerSIS#123" % config._SQL_IP)  # 连接数据库
    except pyodbc.Error:
        conn = pyodbc.connect(
            r"DRIVER={SQL SERVER NATIVE CLIENT 10.0};SERVER=%s;DATABASE=alert;UID=sa;PWD=powerSIS#123" % config._SQL_IP)  # 连接数据库
    cursor = conn.cursor()  # 获得操作的游标
    cursor.execute(f"SELECT Model_info FROM [alert].[dbo].[Model_Version] where model_id={model_id} and version='{version}'")
    res_list = cursor.fetchall()  # 获取查询的结果
    conn.commit()  # 提交执行
    cursor.close()  # 关闭游标
    conn.close()  # 关闭数据库连接
    return json.loads(res_list[0][0])


def pca(model,LockVariable, Data_origin):
    Data = (Data_origin - model["Train_X_mean"]) / model["Train_X_std"]
    featValue = np.array(model["featValue"])  # 训练数据的特征值
    featVec = np.array(model["featVec"])  # 训练数据的特征向量
    k = (model["K"])  # 主元个数
    # selectVec = np.array(model["selectVec"])
    selectVec = featVec[:, 0:k]####自己选择的,取前k个特征向量
    # index = np.argsort(-np.array(featValue))  # 按照featValue进行从大到小排序
    featValue_sort = featValue#[index]  # 排序后的特征值
    '''
    featValue, featVec = np.linalg.eig(model["COV"])  # 求解协方差矩阵的特征值和特征向量

    index = np.argsort(-featValue)  # 按照featValue进行从大到小排序
    featValue_sort =featValue[index]#排序后的特征值
    selectVec = np.matrix(featVec[:,index[:int(model["K"])]])  # 所以这里需要进行转置P
    '''
    ############----------*********-SPE-**************----------########################
    numbel_variable = featValue.shape[0]#取特征值的个数
    # LockVariable="3,5"

    C_ = np.eye(numbel_variable) - np.dot(selectVec, selectVec.T)#生成numbel_variable阶单位矩阵
    X_SPE = C_.T
    D_SPE = C_
    DIAG_SPE = np.eye(numbel_variable)#生成numbel_variable阶单位矩阵

    '''
    # ************************调用LARS*******************************
    t = 50000
    lamuta = 1
    limit_line = model["QCUL_99"]
    beta_path=[]
    for i in range(Data.shape[0]):
        Y=Data[i,:]
        beta, mse=Lars(X_SPE, Y, D_SPE, DIAG_SPE, t, limit_line, lamuta)
        beta_end=abs(beta[-1,:])
        jenk=jenkspy.jenks_breaks(beta_end,5)
        limt=(jenk[1]+jenk[2])/2
        index=np.where(beta_end>limt)[0]        
        beta_path.append(beta[-1,:])
    '''
    ############----------*********-T2-**************----------########################
    DIAG_T2 = np.linalg.pinv(np.diag(featValue_sort[:int(model["K"])]))
    D_T2 = selectVec.copy()
    X_T2 = np.dot(D_T2, np.linalg.cholesky(DIAG_T2)).T

    ############----------*********-综合指标-**************----------########################
    II = featValue_sort.copy()
    II[:int(model["K"])] = II[:int(model["K"])] * model["T2CUL_99"]
    II[int(model["K"]):] = model["QCUL_99"]
    DIAG_Fai = np.linalg.pinv(np.diag(II))
    D_Fai = featVec.copy()
    X_Fai = np.dot(D_Fai, np.linalg.cholesky(DIAG_Fai)).T
    # ************************调用LARS*******************************
    t = 50000
    lamuta = 1
    #limit_line = model["Kesi_99"]/np.sqrt(numbel_variable)#修改
    limit_line = model["Kesi_99"]
    beta_path = []
    SPE_list = []
    FAI_list=[]
    paraState = np.zeros([np.array(Data_origin).shape[0], np.array(Data_origin).shape[1]])
    if Data.shape[1] >= 12:
        para_length = 3
    elif 12 > Data.shape[1] >= 7:
        para_length = 2
    else:
        para_length = 1
    Y = None
    plots_matrix = []  # 贡献图法的矩阵
    plots_index = []  # 贡献图法的index
    for i in range(Data.shape[0]):
        Y = Data[i, :]  # 测试数据的每一行
        #########*********************计算SPE******************************
        SPE_line = np.dot(Y, C_).dot(Y.T)
        SPE_list.append(SPE_line)
        #########################计算综合指标##########
        FAI_list.append(np.dot(Y.T, D_Fai).dot(DIAG_Fai).dot(D_Fai.T).dot(Y))
        # **************计算LARS***************
        beta, mse = Lars(X_Fai, Y, D_Fai, DIAG_Fai, t, limit_line,LockVariable)
        beta_end = abs(beta[-1, :])
        pi=len(beta_end)
        if pi>7:
            jenk = jenkspy.jenks_breaks(beta_end, 5)
        else:
            jenk = jenkspy.jenks_breaks(beta_end, 2)
        limt = (jenk[1] + jenk[2]) / 2
        index = np.where(beta_end > 0)[0]
        if len(index) > para_length:
            # res = recon_fault_diagnosis_r_c(Y, D_Fai @ DIAG_Fai @ D_Fai.T, limit_line, list(zip(index, beta_end[index])), model,
            #                               True,  X_SPE @ X_SPE.T, rbc=None)
            res = recon_fault_diagnosis_r_c(Y, D_Fai @ DIAG_Fai @ D_Fai.T, limit_line,
                                            list(zip(index, beta_end[index])), model,
                                            True, X_SPE @ X_SPE.T, LockVariable, selectVec, rbc=None)
            if not isinstance(res[0], list):
                if res[1] == "plot":
                    # beta[-1, :] = res[0]
                    plots_matrix.append(res[0])
                    plots_index.append(i)
                else:
                    beta[-1, :], index = res[0].T, res[1]
            # beta[-1, :], index = res[0].T, res[1]
        elif len(index) <= para_length and len(index) != 0:
            res = recon_fault_diagnosis_r_l(Y, D_Fai @ DIAG_Fai @ D_Fai.T, index)
            beta[-1, :], index = res[0].T, res[1]
        paraState[i, index] = 1
        beta_new=beta[-1, :]*paraState[i,:]
        beta_path.append(beta_new)
    del Y
    gc.collect()
    beta_path = np.array(beta_path)
    ################-------------------------------------------------------------###############
    #finalData = np.dot(Data - beta_path, selectVec).dot(selectVec.T)
    finalData=Data - beta_path
    reconData = np.add(np.multiply(finalData, model["Train_X_std"]), model["Train_X_mean"])  # 重构值
    if len(plots_matrix) != 0:
        reconData[plots_index] = plots_matrix
    errorData = Data_origin - reconData  # 偏差值
    # cos检验值
    R = 0
    for index in range(0, reconData.shape[1]):
        vector1 = Data_origin[:, index]
        vector2 = np.array(reconData)[:, index]
        R += np.dot(vector1, vector2.T) / (np.sqrt(np.sum(vector1 ** 2)) * np.sqrt(np.sum(vector2 ** 2)))
    R /= reconData.shape[1]
    items = [('reconData', reconData.tolist())
        , ('errorData', errorData.tolist()), ('R', R.tolist()), ('SPE', SPE_list),('FAI', FAI_list),
             ('paraState', paraState.tolist())]
    result = json.dumps(dict(items))  # json.dumps(result)
    return result


def get_history_value(points,time,interval):
    #url="http://192.168.1.201:8080/openPlant/getMultiplePointHistorys"
    url=f"http://{config._EXA_IP}:9000/exawebapi/exatime/GetSamplingValueArrayFloat"
    headers = {"Content-Type": "application/json;charset=utf-8"}#,"token":get_token()
    point_array = points.split(",")
    time_span = time.split(";")
    value_array = []
    for item in point_array:
        for time_piece in time_span:
            st = time_piece.split(",")[0]
            et = time_piece.split(",")[1]
            para = {"ItemName": item, "StartingTime": st, "TerminalTime": et, "SamplingPeriod": interval}
            response = requests.get(url, headers=headers, params=para)
            value = eval(str(response.text).replace("\"","").replace("null","0"))
            value_group = []
            for row in value:
                value_group.append(row[1])
        value_array.append(value_group)
    return np.transpose(np.array(value_array))

# 根据数据集data.txt


if __name__ == "__main__":
    #lifan 调试模型计算引擎
    jsonstr = '{"Model_id":764,"version":"v-2021-04-02 09:43:50","Test_Data":[[129.3936,0.8944824,152.4081,119.2822,0.4589844]],"Target_Data":[[]]}'
    jsonstr = json.loads(jsonstr)
    model_id = jsonstr["Model_id"]
    version = jsonstr["version"]
    res = get_model_by_id_and_version(model_id, version)
    filename = res["algorithm"]
    Data = jsonstr["Test_Data"]
    if filename == "PCA":
        model = res["para"]["Model_info"]
        lock = []
        point_info = res["pointInfo"]
        for i in range(len(point_info)):
            try:
                if point_info[i]["lock"]:
                    lock.append(i)
            except:
                continue
        result = pca(model, lock, np.array(Data))
    print('aaa')