基于PLS的TE过程故障检测

Alep__

已于 2024-06-25 10:00:36 修改

阅读量585

点赞数 8

文章标签：人工智能 python 机器学习

于 2024-05-20 09:48:13 首次发布

本文链接：https://blog.csdn.net/Alep__/article/details/139054942

版权

PLS常用来进行关键性能指标的预测，并且可以通过计算SPE和T2统计量来检测系统故障。PLS常见的算法有很多种，这里介绍在故障检测中最常用的非线性偏最小二乘法。

数据选取

这里选用的是TE数据，选取过程变量XMEAS（1-22），操纵变量XMV（1-11）作为预测变量X，选取过程变量XMEAS（35-36）作为响应变量Y。在进行PLS之前，要进行标准化操作。X为960*33维矩阵，960为样本数，33为变量数。Y为960*2维矩阵，960为样本数，2为变量数。训练集和测试集维数相同。

链接：https://pan.baidu.com/s/1yHgZRqy6yyk9wsm6igqUSw
提取码：odsz

MATLAB

clear all;
close all;
clc;
%置信度
alpha = 0.05;

%数据导入
X_Y_train = matfile('X_Y_train.mat');
X_train = X_Y_train.X_train;
Y_train = X_Y_train.Y_train;
X_Y_test = matfile('X_Y_test_IDV01.mat');
X_test = X_Y_test.X_test;
Y_test = X_Y_test.Y_test;

%数据预处理
[X_train_z,X_train_mean,X_train_std] = zscore(X_train);
[Y_train_z,Y_train_mean,Y_train_std] = zscore(Y_train);
X_test  = (X_test-X_train_mean)./(X_train_std);

%样本数与变量数
[n,m] = size(X_train);

u = Y_train_z(:,1);

num_c = 10;%主元选取
W_T=[];Q_T=[];T=[];P_T=[];
for k = 1:num_c
    u_old = randn(1,m);
    for i = 1:100
        w_T = (u' * X_train_z)/(u' *u);
        w_T = w_T/norm(w_T);
        t = (X_train_z * w_T')/(w_T * w_T');
        q_T = (t' * Y_train_z) /(t' * t);
        u = (Y_train_z * q_T') / (q_T * q_T' );
        if norm(u-u_old)<1e-6
            break;
        end
        u_old = u;
    end
    p_T = (t' * X_train_z)/(t'*t);
    T=[T,t];
    Q_T = [Q_T;q_T];
    W_T = [W_T;w_T];
    P_T = [P_T;p_T];
    X_train_z = X_train_z -t*p_T;
    Y_train_z = Y_train_z - t*q_T;
end

R = [];  % 初始化 R 矩阵
R_c = eye(m);  % 初始化 R_cc 为 m 维单位矩阵

for i = 1:num_c
    if i > 1
        R_c = R_c * (eye(m) - W_T(i-1,:)' * P_T(i-1,:));
    end
    R = [R, R_c * W_T(i,:)'];  % 直接计算并追加到 R
end

%预测Y
Beta = W_T'*inv(P_T*W_T') * Q_T;
Y_pred = (X_test * Beta) .* Y_train_std +Y_train_mean;

%T2控制限
T2_limit = calculateT2Limit(n, num_c, alpha);
%T2统计量
T2 = X_test * R * inv((T' * T)/(n-1)) * R' * X_test';
T2 = diag(T2);

%SPE控制限
SPE_limit_s = ((eye(m) - P_T' * R')*X_train_z')' * ((eye(m) - P_T' * R')*X_train_z');
SPE_limit_s = diag(SPE_limit_s);

mean_SPE_limit = mean(SPE_limit_s);%均值
var_SPE_limit = var(SPE_limit_s);%方差

g =var_SPE_limit/(2*mean_SPE_limit);
h = (2 * mean_SPE_limit * mean_SPE_limit)/var_SPE_limit;

SPE_limit = g * chi2inv(1-alpha,h);

%SPE统计量
SPE = ((eye(m) - P_T' * R')*X_test')' * ((eye(m) - P_T' * R')*X_test');
SPE = diag(SPE);

%绘图真实值/预测值
figure;
plot(Y_test(:,1), 'bo-');
hold on;
plot(Y_pred(:,1), 'r*-');
hold on;
plot(Y_test(:,2), 'ko-');
hold on;
plot(Y_pred(:,2), 'g*-');
legend('X35','X35预测值','X36','X36预测值');
xlabel('样本编号');
ylabel('响应值');
title('PLS回归预测');

%绘制T2/SPE
figure;
subplot(2,1,1);
plot(T2);
hold on;
line([1 n], [T2_limit T2_limit], 'Color', 'red');
title('测试集 Hotelling''s T² 控制图');
xlabel('样本编号');
ylabel('T² 值');

subplot(2,1,2);
plot(SPE);
hold on;
line([1 n], [SPE_limit SPE_limit], 'Color', 'red');
title('测试集 SPE 控制图');
xlabel('样本编号');
ylabel('SPE 值');