大作业基于PCA故障诊断.doc

大作业基于PCA故障诊断

基于主元分析（PCA）
的故障诊断

小组成员:

日期：

目录
1. 运用PCA方法的前提 2
2. PCA方法的基本理论 2
2.1 思路概述 2
2.2 基本理论 2
3. 利用PCA方法进行故障诊断的步骤 4
3.1 建立正常工况的主元模型 4
3.2 在线故障检测与诊断 4
4. PCA的局限性或优缺点 4
5. 基于TE过程的故障诊断 4
5.1 TE过程简介 4
5.2 基于PCA的故障诊断 6
5.2.1 仿真的参数设置 6
5.2.2 仿真结果 6
5.3 仿真总结 12
6. 总结 12
参考文献 13
附录 14

1. 运用PCA方法的前提
1、 样本观测相对独立
2、 潜在变量服从高斯分布
2. PCA方法的基本理论
2.1 思路概述
PCA方法是将高维过程数据投影到正交的低维子空间，并保留主要过程信息。而在几何上，把样本构成的坐标系，通过某种线性组合旋转到新的坐标空间，新的坐标轴代表了具有最大方差的方向[[] 周东华，李钢，李元．数据驱动的工业过程故障诊断技术—基于主元分析与偏最小二乘的方法[M]．北京：科学出版社，2011．
]。
2.2 基本理论
假设代表一个包含了m个传感器的测量样本，每个传感器各有n个独立采样，构造出测量数据矩阵 ，其中每一列代表一个测量变量，每一行代表一个样本。
（1） 对数据矩阵进行协方差分解，并选择主元的个数
X的协方差矩阵为 ，对其进行特征值分解，并且按照特征值的大小降序排序，如下：

其中，是一个对角阵，也是S的特征值矩阵，而且其对角线上的元素满足；是S的特征向量矩阵，维数为m x m，是的前A列，包含所有主元的信息，是余下的m-A列，包含非主元信息。
（2） 将原数据进行分解，得到主元子空间和残差子空间

对X进行特征值分解以后，X可以分解如下：

其中，，被称为主元子空间；，被称为残差子空间；，被称为得分矩阵；被称为负载矩阵，由S的前A个特征向量构成。
（3） 故障检测的两个指标或判据
下面提到的 x 是新采集的样本，维数为m x 1。
A、 SPE统计量
SPE指标衡量样本向量在残差空间的投影的变化

其中，表示置信度为的控制限。
常用的计算公式如下：

其中，，为X的协方差矩阵的特征值，为标准正态分布在置信度为下的阈值。
B、 统计量
Hotelling's 统计量衡量样本向量在主元空间的变化

其中，，为置信度为的控制限。
控制限的常用计算方法如下：

其中，是带有A和n-A个自由度、置信度为的F分布值。
（4） 计算贡献率
下面提到的 x 是新采集的样本，维数为m x 1。 P是负载矩阵。
最常用的贡献图的统计量是SPE和。
基于SPE的贡献图定义如下：

其中，表示每个变量对SPE统计量的贡献值，，表示单位矩阵I的第i列。
基于的贡献图的定义如下：

其中，（因为这里是参考书籍，所以与实际程序并不相符），其余同SPE贡献图的定义。
当检测到故障后，贡献图中较大的变量被认为是可能造成故障的变量。但需要具有过程背景
知识的人员确定最终的故障原因。
3. 利用PCA方法进行故障诊断的步骤
3.1 建立正常工况的主元模型
Step 1 将正常样本数据进行标准化，变换为均值为0，方差为1的标准数据集
Step 2 对Step 1中的标准数据集，建立PCA主元模型，提取主元
Step 3 计算Step 1中的标准数据集的PCA模型的统计量及相应的控制限
3.2 在线故障检测与诊断
Step 1 在线采集数据，从采样中获得新的数据x，并进行标准化
Step 2 对标准化后的数据，计算统计量和SPE统计量，监控其数据是否超过正常状态的控制限。若没有超限，重复Step 1，否则进入Step 3。
Step 3 计算每个过程变量对统计量和SPE统计量的贡献率，贡献率最大的变量就是可能引起故障的变量。
4. PCA的局限性或优缺点
PCA方法理论基础简单，适合处理多变量统计问题，而且完全基于系统传感器的测量数据，工程应用容易实现。
但是其未考虑变量间的相关性，可能会丢失一些信息；另外，对于非正态分布的数据，非线性的数据，不能很好的进行故障诊断；而且其诊断出来的故障物理意义不明确，难于解释。
5. 基于TE过程的故障诊断
5.1 TE过程简介
TE过程模型，是根据实际化工过程的建立的模型，其已经被广泛的用作进行控制与监控研究的基准过程（Benchmark Process），图5-1展示了这个过程的流程图。该过程共包含5个主要单元（反应器、冷凝器、压缩机、分离器、汽提塔）。TE过程是通过4种反应物产生2中产物，同时还有一种惰性产物和一种副产物，共8种成分。由于专利的原因，他们分别被命名为A,B,C,D,E,F,G,H，具体的工艺运行过程可详见蒋浩天等的文章或者书籍。
Downs和Vogel[[] J.J.DOWNS，E.F.VOGEL．A plant-wide industrial process control problem[J]．Computers&Chemical Engineering，1993,17(3)：245-255．
]在论文中提到有6中操作模式，本文只考虑基本工况模式。由于TE过程是开环不稳定的，所以许多学者对其控制方法进行了许多研究，如McAvoy[[] T.J.MCAVOY ，N.YE．base control for the Tennessee Eastman problem[J]．Computers&Chemical Engineering，1998,18(5)：383-413．
]提出的基本控制（Base Control）,Ricker N L[[] N.L.RICKER．Decentralized control of the Tennessee Eastman Challenge Process[J]．Journal of Process Control，1996,6(4)：205-221．
]提出的离散控制策略以及Lyman.P.R提出的厂级控制策略。本文采用离散控制策略，包含19个PI控制器和11个设定点。该过程共有53个变量，其中41个是过程变量，12个是操纵变量，详见参考文献，故障类型见表5-1。采样时间为3分钟，这充分考虑了该过程闭环时间常数为2小时的因素。整体过程并未采用先进的质量控制策略，因此可以认为这是一个稳定的系统。

图5-1

表5-1 TE故障汇总
变量符号
过程变量
类型
IDV(0)
正常操作
无
IDV(1)
A/C加料比率，B成分不变
阶跃
IDV(2)
B成分，A/C进料比不变
阶跃
IDV(3)
D的进口温度
阶跃
IDV(4)
反应器冷却水的入口温度
阶跃
IDV(5)
冷凝器冷却水的入口温度
阶跃
IDV(6)
A进料损失
阶跃
IDV(7)
C存在压力损失-可用性降低
阶跃
IDV(8)
A、B、C进料成分
随机变量
IDV(9)
D的进料温度
随机变量
IDV(10)
C的进料温度
随机变量
IDV(11)
反应器冷却水的入口温度
随机变量
IDV(12)
冷凝器冷却水的入口温度
随机变量
IDV(13)
反应动态
随机变量
IDV(14)
反应器冷却水阀门
粘住
IDV(15)
冷凝器冷却水阀门
粘住
IDV(16)
未知
未知
IDV(17)
未知
未知
IDV(18)
未知
未知
IDV(19)
未知
未知
IDV(20)
未知
未知

5.2 基于PCA的故障诊断
5.2.1 仿真的参数设置
本文利用的数据是从TE过程产生的数据，采样时间为3min，所以共1000个采样点可用，本文利用前900个。数据集都经过平滑滤波和归一化。
过程变量即被监控的变量为变量1~变量16。
故障类为故障1、故障4、故障17，分别是阶跃型，瞬间突变型和未知类型。这样的选择有利于测试该方法的有效性。
对于控制限，统计量和SPE统计量的置信度均为99%。
5.2.2 仿真结果
1）建立PCA模型
Step 1 利用累计方差贡献率的方法选取主元
利用正常数据集的前800个采样点进行训练，得到的结果如下所示：

图 5-2

由图5-2可知，在累计贡献率为85%时，主元个数为10；所以本文的累计贡献率选为85%。
箱线图如下

图5-3
Step 2 T^2统计量和SPE统计量监控
结果如下：

图5-4

2）在线监控
故障1的诊断结果：

图5-5

图5-6

图5-7
故障4的诊断结果：

图5-8

图5-9

图5-10

故障17的诊断结果：

图5-11

图5-12

图5-13

5.3 仿真总结
由图5-5到图5-13可以看出，PCA方法可以有效的检测出各种类型的故障，并且根据贡献图可以初步判断引发故障的变量。
当然，截止目前已有许多改进的PCA方法，如MPCA，MSPCA，基于权重的PCA，KPCA，DPCA等等，总之，该技术方兴未艾，正处于蓬勃发展之中。
6. 总结
I、统计量使用情况总结
首先，统计量反应的是主元空间的变化，因此不能检测到非主元变量的故障；SPE统计量反应的是所有的变量，因此统计量超限，SPE必超限（但有例外，如过程参数的变化）；但是SPE统计量本身主要代表噪声，所以当其超限时，可能是非主元变量故障，也可能噪声引起的。
II、使用以上两种统计量会出现以下几种情况：
（1） 故障使SPE和统计量同时超限；
（2） 故障使SPE超限，而统计量没有；
（3） 故障使统计量超限，而SPE没有；
（4） 两者都没有超限。
其中，SPE统计量对（1）、（2）、（4）是有效的。
III、两者的优缺点
统计量适合来监控质量指标的变化；SPE统计量对应的误报率和漏报率会少一些（针对非正态的或者不平稳的过程）。
IV、引入箱线图的意义
由于箱线图可以反映历史数据的统计分布情况，可以用来判断数据的异常值；本文引入箱线图是为了辅助贡献图的方法，进行故障分离或者定位。

参考文献

18 / 19
附录
MATLAB源程序
% 作者，北京化工大学，控制科学与工程，张静静
% 参考文献：周东华的书籍《数据驱动的工业过程故障诊断
% 技术--基于主元分析与偏最小二乘的方法》.
% 一些变量的说明：
% x 是 n*m 的矩阵,其中 n 是一个测量变量的不同时刻的采样值,m 是有m个传感器；即行是一个样本
% x,test_x 是归一化后的值
% N 是测试数据集的行数

% P 是负载矩阵,大小为m*kk，kk为主元的个数
% T_test 是T^2统计量式子左边的计算结果，SPE是 e'*e 的结果
% T_alpha 是T^2统计量的控制限，SPE_alpha 是SPE的控制限
% contri4 是每个变量对SPE的贡献值, contri6 是 contri4 贡献率转换成百分比的形式
% t2 是每个变量对T^2统计量的贡献值, t4 是 t2 贡献率转换成百分比的形式
% monitortime表示在该时刻绘制的贡献图，因为贡献图都是某一时刻的
%%
clear all
clc
%% 载入数据
%*************构造训练样本，测试样本**************
simout2=xlsread('output_normal',2);
dataset=simout2(1:5:4500,1:16); %正常数据集
x=dataset(1:800,:);
% test_x=x;
% N=size(x,1);
% load output_fault4;
load output_fault17;
dataset2=simout;
% dataset2=xlsread('output_fault1',1);
test_x=dataset2(1:5:4500,1:16);
N=size(test_x,1);
monitortime=200;
%% 调用自编的PCA函数
[P,T_test,SPE,T_alpha,SPE_alpha,contri6,t4,x,test_x,Latent]=PCA(x,test_x,N,monitortime);
%% 绘制必要的曲线图形
%********************************************
figure(1);
subplot(211),
x_row=1:N;
T_alpha=T_alpha*(ones(size(x_row)));
plot(x_row,T_alpha,'r-','LineWidth',2.5),hold on
plot(T_test,'b-'); grid on; xlabel('样本点');ylabel('T^2统计量');
subplot(212),
x_row=1:N;
SPE_alpha=SPE_alpha*(ones(size(x_row)));
plot(x_row,SPE_alpha,'r-','LineWidth',2.5),hold on
plot(SPE,'b-'); grid on; xlabel('样本点');ylabel('SPE统计量');
% 误报率
count=0;
for i=1:N
if T_test(i)>T_alpha(i)
count=count+1;
end
end
false_alarm=count/N
count2=0;
for i=1:N
if SPE(i)>SPE_alpha(i)
count2=count2+1;
end
end
false_alarm2=count2/N
figure(2);
subplot(121),bar(t4),title('T^2统计量贡献图');
subplot(122),bar(contri6),title('SPE统计量贡献图');
figure(3)
boxplot(test_x)

function [P,T_test,SPE,T_alpha,SPE_alpha,contri6,t4,x,test_x,Latent]=PCA(x,test_x,N,monitortime)
%%
%**************对数据进行标准化*************
% 对训练数据进行按列归一化
row=size(x,1);
column=size(x,2);
mean_x=mean(x); % 每一列的均值
var_x=var(x); % 每一列的方差
for i=1:row
for j=1:column
if var_x(j)~=0 % 这里栽过跟头
x(i,j)=(x(i,j)-mean_x(j))/sqrt(var_x(j));
end
end
end
% 对测试数据进行归一化
row1=size(test_x,1);
column1=size(test_x,2);
mean_test_x=mean(test_x); % 每一列的均值
var_test_x=var(test_x); % 每一列的方差
for i=1:row1
for j=1:column1
if var_test_x(j)~=0
test_x(i,j)=(test_x(i,j)-mean_test_x(j))/sqrt(var_test_x(j));
end
end
end
%%
filter=3;
for i=1:column
for j=2:row-1
x(j,i)=sum(x(j-1:j+1,i))/filter;
end
end
for i=1:column1
for j=2:row1-1
test_x(j,i)=sum(test_x(j-1:j+1,i))/filter;
end
end
%%
%**********对标准化后的数据进行协方差分解***********
S=x'*x/(row-1);
[e,d,M]=svd(S);
% [e,d]=eig(S); % e是特征向量，d是特征值
% d=diag(d);
% d=d(end:-1:1); % 保证特征值从大到小排序
% e=fliplr(e);
%%
%*****************选择主元**************************
% 这里利用累计方差贡献率的方法
Latent=diag(d); % 将d从矩阵形式变换为列向量的形式
SUM=0;
for i=1:length(Latent)
SUM=SUM+Latent(i);
percent=SUM/sum(Latent);
if(percent>=0.85)
kk=i;
break;
end
end
disp('主元个数为');disp(kk);
disp('***********************************************');
%*********************构造 P 矩阵*********************
P=e(:,1:kk); % P 的维数是m*kk

%%
%****************计算T^2统计量，并确定控制限*****************
% 确定控制限
alpha1=0.99; % 置信度
n=row;
F=finv(alpha1,kk,(n-kk));
T_alpha=kk*(n-1)*(n+1)*F/(n*(n-kk)); % F 指 F(kk,n-kk),alpha1的F分布
% 计算T^2统计量
D=d(1:kk,1:kk);
T_test=zeros(N,1);
for i=1:N
T_test(i)=test_x(i,:)*P*inv(D)*P'*test_x(i,:)';
end
disp('**********************************************');
%****************计算SPE统计量，并确定控制限*****************
% 确定控制限
m=column; % 样本点的个数
theta1=sum(Latent(kk+1:m));
theta2=sum(Latent(kk+1:m).^2);
theta3=sum(Latent(kk+1:m).^3);
h0=1-2*theta1*theta3/(3*theta1^2);

alpha2=0.99; % 置信度
C=norminv(alpha2,0,1);
SPE_alpha=theta1*(((C*sqrt(2*theta2*h0^2))/theta1)+1+(theta2*h0*(h0-1))/theta1^2)^(1/h0);

% 计算SPE统计量
SPE=zeros(N,1);
for i=1:N
SPE(i)=(test_x(i,:)'-P*P'*test_x(i,:)')'*(test_x(i,:)'-P*P'*test_x(i,:)');
% SPE(i)=(test_x(i,:)-test_x(i,:)*P*P')*(test_x(i,:)-test_x(i,:)*P*P')';
end
%% 利用贡献图确定故障的样本点的位置
%****************************************
% 下面基于T^2的贡献图法
t2=zeros(1,m); kexi=eye(m);
t1=P*D*P';
for i=1:m
t2(i)=test_x(monitortime,:)*t1*kexi(:,i)*kexi(:,i)'*test_x(monitortime,:)';
end
t3=sum(t2);
t4=t2/t3; % 求百分比
% 下面利用SPE的贡献图法
contri1=eye(m);
contri2=contri1-P*P';
contri3=zeros(1,m);
for i=1:m
contri3(i)=contri1(:,i)'*contri2*test_x(200,:)';
% 这一句中测试集的行一般大于列的个数，否则出错
end
contri4=contri3.^2;
contri5=sum(contri4);
contri6=contri4/contri5; % 得到contri4已经可以了，这一步是将贡献率转换成百分比的形式