PCA故障诊断中两个关键统计变量T2T^2T2和SPESPESPE的的计算

T2T^2T2：Hotelling-T2
SPESPESPE：平方预测误差(Squared prediction error)

T2T^2T2 统计量反映了每个主成分在变化趋势和幅值上偏离模型的程度，是对模型内部化的一种度量，它可以用来对多个主元同时进行监测；SPESPESPE 统计量刻画了输入变量的测量值对主元模型的偏离程度，是对模型外部变化的一种度量。

T2T^2T2在线计算以及控制限的计算

SPESPESPE在线计算以及控制限的计算

T2T^2T2和SPESPESPE使用情况总结

T2T^2T2统计量反应的是主元空间的变化，因此不能检测到非主元变量的故障；
SPESPESPE统计量反应的是所有的变量，因此T2T^2T2统计量超限，SPESPESPE必超限（但有例外，如过程参数的变化）；
而SPESPESPE统计量衡量的是变量间相关性被改变的程度，显示异常的工况；所以当其超限时，可能是过程变量故障，也可能其它故障引起的。

根据T2T^2T2和SPESPESPE是否超限结果分析

故障使SPESPESPE和T2T^2T2统计量同时超限；
故障使SPESPESPE超限，而T2T^2T2统计量没有；
故障使T2T^2T2 统计量超限，而SPESPESPE没有；
两者都没有超限。
其中，SPESPESPE统计量对1,2,4是有效的。

两者的优缺点

T2T^2T2统计量适合来监控质量指标的变化；SPESPESPE统计量可以监控各类指标，其误报率和漏报率会少一些（针对非正态或者不平稳的过程）。

TE过程（田纳西伊斯曼过程）数据简介

具体数据见链接TE过程数据

在TE仿真平台中，可获取41个过程变量的运行数据和11个操纵变量的数据。另外，通过对进料等变量的输入控制，可以模拟得到21个故障类型，其中正常数据和异常（故障）数据的采样时间均为3min，即20次/min。
TE训练集和测试集结构：
整个TE数据集由训练集和测试集构成，TE集中的数据由22次不同的仿真运行数据构成，TE集中每个样本都有52个观测变量。d00.dat至d21.dat为训练集样本，d00_te.dat至d21_te.dat为测试集样本。前41个变量为过程变量，后11个变量为操纵变量。
d00.dat和d00_te.dat为正常工况下的样本。d00.dat训练样本是在25h运行仿真下获得的。观测数据总数为500。而d00_te.dat测试样本是在48h运行仿真下获得的，观测数据总数为960。
d01.dat至d21.dat为带有故障的训练集样本，d01_te.dat至d21_te.dat为带有故障的测试集样本。每个训练集\测试样本代表一种故障。
要值得注意的是对于带有故障的训练集样本，是在25h运行仿真下获得的。仿真开始时没有故障情况,故障是在仿真时间为1h的时候引入的。但观测数据是在引入故障后才开始采集的，即只有480个观测值。
带有故障的测试集样本是在48h运行仿真下获得的，故障在8h的时候引入，共采集960个观测值，其中前160个观测值为正常数据。

基于TE过程数据的PCA故障诊断

Python代码

read_TEdat.py：读取和处理TE过程原始.dat数据文件

import pandas as pd
import glob
import os
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来正常显示负号def read_all_data(path_test,path_train):'''读取TE过程的所有.dat数据并存人DataFrame中，输入参数为测试数据和训练数据的绝对路径'''var_name = []for i in range(1,42):var_name.append('XMEAS(' + str(i) + ')')for i in range(1,12):var_name.append('XMV(' + str(i) + ')')data_test, data_train = [], []# path_test = r'C:\Users\17253\Desktop\组内\K_shape\data\TE\test'# path_train = r'C:\Users\17253\Desktop\组内\K_shape\data\TE\train'test_join = glob.glob(os.path.join(path_test,'*.dat'))train_join = glob.glob(os.path.join(path_train,'*.dat'))for filename in test_join:data_test.append(pd.read_table(filename, sep = '\s+', header=None, engine='python', names = var_name))for filename2 in train_join:data_train.append(pd.read_table(filename2, sep = '\s+', header=None, engine='python', names = var_name))return data_test, data_train

diagnosis_pca.py：用PCA算法进行故障诊断

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from scipy import stats
from scipy.stats import norm, chi2
from sklearn.preprocessing import StandardScaler
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来正常显示负号def t2_online(x, p, v):'''p：特征向量组成的降维矩阵，负载矩阵x：在线样本，shape为m*1v：特征值由大至小构成的对角矩阵'''T_2 = np.dot(np.dot((np.dot((np.dot(x.T, p)), np.linalg.inv(v))), p.T), x)return T_2def spe_online(x, p):'''p：特征向量组成的降维矩阵，负载矩阵x：在线样本，shape为m*1'''I = np.eye(len(x))spe = np.dot(np.dot(x.T, I - np.dot(p, p.T)), x)# Q_count = np.linalg.norm(np.dot((I - np.dot(p_k, p_k.T)), test_data_nor), ord=None, axis=None, keepdims=False)  #二范数计算方法return spedef pca_control_limit(Xtrain, ratio = 0.95, confidence = 0.99):'''计算出T2和SPE统计量'''pca = PCA(n_components = ratio)pca.fit(Xtrain)evr = pca.explained_variance_ratio_ev = pca.explained_variance_ # 方差，相当于X的协方差的最大的前几个特征值n_com = pca.n_componentsp = (pca.components_).T # 负载矩阵v = np.diag(ev) # 特征值组成的对角矩阵v_all = PCA(n_components = Xtrain.shape[1]).fit(Xtrain).explained_variance_p_all = (PCA(n_components = Xtrain.shape[1]).fit(Xtrain).components_).Tk = len(evr)n_sample = pca.n_samples_##T统计量阈值计算coe = k* (n_sample - 1) * (n_sample + 1) / ((n_sample - k) * n_sample)t_limit = coe * stats.f.ppf(confidence, k, (n_sample - k))##SPE统计量阈值计算theta1 = np.sum((v_all[k:]) ** 1)theta2 = np.sum((v_all[k:]) ** 2)theta3 = np.sum((v_all[k:]) ** 3)h0 = 1 - (2 * theta1 * theta3) / (3 * (theta2 ** 2))c_alpha = norm.ppf(confidence)spe_limit = theta1 * ((h0 * c_alpha * ((2 * theta2) ** 0.5)/ theta1 + 1 + theta2 * h0 * (h0 - 1) / (theta1 ** 2)) ** (1 / h0))return t_limit, spe_limit, p, v, v_all, k, p_alldef pca_model_online(X, p, v):t_total = []q_total = []for x in range(np.shape(X)[0]):data_in = X[x]t = t2_online(data_in, p, v)q = spe_online(data_in, p)t_total.append(t)q_total.append(q)return t_total, q_totaldef figure_control_limit(X, t_limit, spe_limit, t_total, q_total):## 画控制限的图plt.figure(2, figsize=(12, 7))ax1 = plt.subplot(2, 1, 1)plt.plot(t_total)plt.plot(np.ones((len(X))) * t_limit, 'r', label='95% $T^2$ control limit')# ax1.set_ylim(0,100)# plt.xlim(0,100)ax1.set_xlabel(u'Samples')ax1.set_ylabel(u'Hotelling $T^2$ statistic')plt.legend()ax2 = plt.subplot(2, 1, 2)plt.plot(q_total)plt.plot(np.ones((len(X))) * spe_limit, 'r', label='95% spe control limit')# ax1.set_ylim(0,30)# plt.xlim(0,100)ax2.set_xlabel(u'Samples')ax2.set_ylabel(u'SPE statistic')plt.legend()plt.show()def Contribution_graph(test_data, trian_data, index, p, p_all, v_all, k, t_limit):# 贡献图index = 160#1.确定造成失控状态的得分test_data = fault02_testdata_mean = data_mean = np.mean(Xtrain_nor, 0)data_std = np.std(Xtrain_nor, 0)test_data_submean = np.array(test_data - data_mean)test_data_norm = np.array((test_data - data_mean) / data_std)t = test_data_norm[index,:].reshape(1,test_data.shape[1])S = np.dot(t,p[:,:])r = []for i in range(k):if S[0,i]**2/v_all[i] > t_limit/k:r.append(i)print(r)#2.计算每个变量相对于上述失控得分的贡献cont = np.zeros([len(r),test_data.shape[1]])for i in range(len(r)):for j in range(test_data.shape[1]):cont[i,j] = S[0,i]/v_all[r[i]]*p_all[r[i],j]*test_data_submean[index,j]if cont[i,j] < 0:cont[i,j] = 0#3.计算每个变量对T的总贡献a = cont.sum(axis = 0)#4.计算每个变量对Q的贡献I = np.eye(test_data.shape[1])e = (np.dot(test_data_norm[index,:],(I - np.dot(p, p.T))))**2##画图plt.rcParams['font.sans-serif']=['SimHei']plt.rcParams['axes.unicode_minus'] = Falsefont1 = {'family' : 'SimHei','weight' : 'normal','size'   : 23,}plt.figure(2,figsize=(16,9))ax1=plt.subplot(2,1,1)plt.bar(range(test_data.shape[1]),a)plt.xlabel(u'变量号',font1)plt.ylabel(u'T^2贡献率 %',font1)plt.legend()plt.showax1=plt.subplot(2,1,2)plt.bar(range(test_data.shape[1]),e)plt.xlabel(u'变量号',font1)plt.ylabel(u'Q贡献率 %',font1)plt.legend()plt.show()

TE_mian：主函数

import read_TEdat
import diagnosis_pca as pca
from sklearn.preprocessing import StandardScaler
path_test = r'C:\Users\17253\Desktop\组内\K_shape\data\TE\test'
path_train = r'C:\Users\17253\Desktop\组内\K_shape\data\TE\train'
data_test, data_train = read_TEdat.read_all_data(path_test, path_train)
fault02_train, nor_train = data_train[2], data_train[0]
fault02_test, nor_test = data_test[2], data_test[0]#数据标准化
scaler = StandardScaler().fit(nor_train)
Xtrain_nor = scaler.transform(nor_train)
Xtest_nor = scaler.transform(nor_test)
Xtrain_fault = scaler.transform(fault02_train)
Xtest_fault = scaler.transform(fault02_test)
# PCA
t_limit, spe_limit, p, v, v_all, k, p_all = pca.pca_control_limit(Xtrain_nor)# 在线监测
t2, spe = pca.pca_model_online(Xtest_fault, p, v)pca.figure_control_limit(Xtest_fault, t_limit, spe_limit, t2, spe)
pca.Contribution_graph(fault02_test, nor_train, 600, p, p_all, v_all, k, t_limit)
print(i)

结果分析

本次实验选择了故障2的数据作为测试数据d02_te.dat，训练数据使用训练集中的正常数据d00.dat。故障2的具体故障表现为组分B含量发生变化，A/C进料流量比不变（流4），类型为阶跃型。

T2T^2T2和SPESPESPE统计量监测

最后测试集的T2T^2T2和SPESPESPE统计量如图：

可以看出，在160个样本左右诊断出故障，与实际情况相符（具体的故障数据说明参考TE过程数据）

故障贡献图

一般变量的贡献越大，说明受其影响越大。
相关代码如下(diagnosis.Contribution_graph)：

def Contribution_graph(test_data, trian_data, index, p, p_all, v_all, k, t_limit):# 贡献图index = 160#1.确定造成失控状态的得分test_data = fault02_testdata_mean = data_mean = np.mean(Xtrain_nor, 0)data_std = np.std(Xtrain_nor, 0)test_data_submean = np.array(test_data - data_mean)test_data_norm = np.array((test_data - data_mean) / data_std)t = test_data_norm[index,:].reshape(1,test_data.shape[1])S = np.dot(t,p[:,:])r = []for i in range(k):if S[0,i]**2/v_all[i] > t_limit/k:r.append(i)print(r)#2.计算每个变量相对于上述失控得分的贡献cont = np.zeros([len(r),test_data.shape[1]])for i in range(len(r)):for j in range(test_data.shape[1]):cont[i,j] = S[0,i]/v_all[r[i]]*p_all[r[i],j]*test_data_submean[index,j]if cont[i,j] < 0:cont[i,j] = 0#3.计算每个变量对T的总贡献a = cont.sum(axis = 0)#4.计算每个变量对Q的贡献I = np.eye(test_data.shape[1])e = (np.dot(test_data_norm[index,:],(I - np.dot(p, p.T))))**2##画图plt.rcParams['font.sans-serif']=['SimHei']plt.rcParams['axes.unicode_minus'] = Falsefont1 = {'family' : 'SimHei','weight' : 'normal','size'   : 23,}plt.figure(2,figsize=(16,9))ax1=plt.subplot(2,1,1)plt.bar(range(test_data.shape[1]),a)plt.xlabel(u'变量号',font1)plt.ylabel(u'T^2贡献率 %',font1)plt.legend()plt.showax1=plt.subplot(2,1,2)plt.bar(range(test_data.shape[1]),e)plt.xlabel(u'变量号',font1)plt.ylabel(u'Q贡献率 %',font1)plt.legend()plt.show()

贡献图

PCA针对TE过程（田纳西伊斯曼过程）数据的故障监测处理相关推荐

Tennessee Eastman（TE）田纳西-伊斯曼仿真平台应用试验与分析
上次说到典型相关分析,在故障诊断邻域用及,下面来实操. 这里是几个数据集下载及相关文档的Link: MATLAB代码:http://depts.washington.edu/control/LARRY ...
弗朗西斯科·罗曼·阿拉尔孔·苏亚雷斯!!伊斯科!!
目录弗朗西斯科·罗曼·阿拉尔孔·苏亚雷斯弗朗西斯科·罗曼·阿拉尔孔·苏亚雷斯弗朗西斯科·罗曼·阿拉尔孔·苏亚雷斯(Francisco Román Alarcón Suarez),昵称伊斯科(Is ...
荷兰帝斯曼收购初创科技公司尚善无氟光伏背板被指有利回收再利用
不久前,在荷兰拥有110多年悠久历史的帝斯曼集团,宣布全资收购了我国一科技型初创企业--苏州尚善新材料科技有限公司. 据记者了解,促使帝斯曼下定决心完成这一交易的原因之一,正是尚善所具备的APE高性能 ...
NBA名人堂之-朱利叶斯·欧文|埃尔文·海耶斯|多尔夫·谢伊斯|鲍勃·佩蒂特|戴夫·冰
朱利叶斯·欧文英文名:Julius Erving 出生:1950年2月22日身高:2.01米 NBA生涯:1976-1983年司职:前锋主要荣誉: ①1983年获NBA总冠军: ②5次入选NB ...
什么叫小米粒吸尘机器人_【诉说心声】德国克林斯曼全自动吸小米粒小石头怎么样?老铁诉说使用全过程 | 智能扫地机器人评测...
[诉说心声]德国克林斯曼全自动吸小米粒小石头怎么样?老铁诉说使用全过程操作简单千帕吸力电池耐用上个月在市场上面了解了下怎么样,当时超出了我的预算,没有入手,这个月活动,果断入手了,对于很多人最关 ...
pca针对初学者_针对初学者和专家的12酷数据科学项目创意
pca针对初学者 The domain of Data Science brings with itself a variety of scientific tools, processes, alg ...
集训感言：来自祖斯曼
今天已经是5月27号了, 已经陆陆续续写了27天的题,尽管困难的题目基本都是看着英雄哥的代码理解来解决的,但是收获还是很大的,对于很多细节上的处理,以及我不曾学习过的语法,怎么说呢,加入星球,每天训练 ...
训练局限玻尔斯曼机(Training Restricted Boltzmann Machines)
本栏目(Training Restricted Boltzmann Machines)包括内容的介绍.局限玻尔斯曼机与对比差异(CD)的回顾.使用CD时怎样纠正统计.迷你块大小(mini-ba ...
佳明、杜斯曼集团等巨头遭勒索攻击；Emotet利用新技术实施攻击
聚焦源代码安全,网罗国内外最新资讯! 编译:奇安信代码卫士团队勒索攻击层出不穷. 两大跨国巨头遭勒索攻击最近,知名 GPS 导航设备及运动穿戴设备制造商佳明 (Garmin) 官方证实称遭到勒索 ...

PCA针对TE过程（田纳西伊斯曼过程）数据的故障监测处理