核主元分析（Kernel principal component analysis ，KPCA）在降维、特征提取以及...

Colbie

; O% e4 v; r' z7 |+ N
核主元分析（Kernel principal component analysis ，KPCA）在降维、特征提取以及故障检测中的应用。主要功能有：（1）训练数据和测试数据的非线性主元提取（降维、特征提取）
5 p+ O7 x' X8 X  Q" d7 H( m8 ~) v（2）SPE和T2统计量及其控制限的计算
) X$ L  X7 g$ k: C（3）故障检测
' M' a5 k0 p! S: W/ C1 d! s: ]
参考文献：
2 \0 q; W* |, D9 h4 MLee J M, Yoo C K, Choi S W, et al. Nonlinear process monitoring using kernel principal component analysis[J]. Chemical engineering science, 2004, 59(1) : 223-234.
9 z' y5 V" P4 \  H7 w
1. KPCA的建模过程（故障检测）：
2 g! ~4 h/ G2 T: r4 }1 C% R3 [（1）获取训练数据（工业过程数据需要进行标准化处理）
（2）计算核矩阵
' g  R. t' h4 J7 k, p6 Z4 ?/ c（3）核矩阵中心化
（4）特征值分解
4 s4 m6 P2 g! ]9 Z) k0 z（5）特征向量的标准化处理
% Q) g9 D8 |2 K* r) E" o（6）主元个数的选取
6 X2 H; \+ o1 K  t# ^  d8 m+ F  |（7）计算非线性主成分（即降维结果或者特征提取结果）
1 l2 d4 Y& o5 D（8）SPE和T2统计量的控制限计算

• function model = kpca_train(X,options)
• % DESCRIPTION
• % Kernel principal component analysis (KPCA)
• %
• %       mappedX = kpca_train(X,options)
• %
• % INPUT
• %   X            Training samples (N*d)
• %                N: number of samples
• %                d: number of features
• %   options      Parameters setting
• %
• % OUTPUT
• %   model        KPCA model
• %
• %
• % Created on 9th November, 2018, by Kepeng Qiu.
• % number of training samples
• L = size(X,1);
• % Compute the kernel matrix
• K = computeKM(X,X,options.sigma);
• % Centralize the kernel matrix
• unit = ones(L,L)/L;
• K_c = K-unit*K-K*unit+unit*K*unit;
• % Solve the eigenvalue problem
• [V,D] = eigs(K_c/L);
• lambda = diag(D);
• % Normalize the eigenvalue
• V_s = V ./ sqrt(L*lambda)';
• % Compute the numbers of principal component
• % Extract the nonlinear component
• if options.type == 1 % fault detection
•     dims = find(cumsum(lambda/sum(lambda)) >= 0.85,1, 'first');
• else
•     dims = options.dims;
• end
• mappedX  = K_c* V_s(:,1:dims) ;
• % Store the results
• model.mappedX =  mappedX ;
• model.V_s = V_s;
• model.lambda = lambda;
• model.K_c = K_c;
• model.L = L;
• model.dims = dims;
• model.X = X;
• model.K = K;
• model.unit = unit;
• model.sigma = options.sigma;
• % Compute the threshold
• model.beta = options.beta;% corresponding probabilities
• [SPE_limit,T2_limit] = comtupeLimit(model);
• model.SPE_limit = SPE_limit;
• model.T2_limit = T2_limit;
• end
  

! ]/ H: P; R, m% `6 o! B6 h. i" z
2. KPCA的测试过程：
1 Q- b8 ~0 \6 B, {7 _9 }（1）获取测试数据（工业过程数据需要利用训练数据的均值和标准差进行标准化处理）
$ i* ?3 k" C# M2 A! ?# u（2）计算核矩阵
0 n6 H$ g# u0 Q; G) ~（3）核矩阵中心化
（4）计算非线性主成分（即降维结果或者特征提取结果）
（5）SPE和T2统计量的计算
1 U% o. T1 o4 b! E1 l$ }8 z2 q. n: D

• function [SPE,T2,mappedY] = kpca_test(model,Y)
• % DESCRIPTION
• % Compute the T2 statistic, SPE statistic,and the nonlinear component of Y
• %
• %       [SPE,T2,mappedY] = kpca_test(model,Y)
• %
• % INPUT
• %   model       KPCA model
• %   Y           test data
• %
• % OUTPUT
• %   SPE         the SPE statistic
• %   T2          the T2 statistic
• %   mappedY     the nonlinear component of Y
• %
• % Created on 9th November, 2018, by Kepeng Qiu.
• % Compute Hotelling's T2 statistic
• % T2 = diag(model.mappedX/diag(model.lambda(1:model.dims))*model.mappedX');
• % the number of test samples
• L = size(Y,1);
• % Compute the kernel matrix
• Kt = computeKM(Y,model.X,model.sigma );
• % Centralize the kernel matrix
• unit = ones(L,model.L)/model.L;
• Kt_c = Kt-unit*model.K-Kt*model.unit+unit*model.K*model.unit;
• % Extract the nonlinear component
• mappedY = Kt_c*model.V_s(:,1:model.dims);
• % Compute Hotelling's T2 statistic
• T2 = diag(mappedY/diag(model.lambda(1:model.dims))*mappedY');
• % Compute the squared prediction error (SPE)
• SPE = sum((Kt_c*model.V_s).^2,2)-sum(mappedY.^2 ,2);
• end
  6 _' P9 ~! {; x) K

. N8 Q9 B& z* l

3. demo1: 降维、特征提取
(1) 源代码
2 C7 T8 ~; v: T

• % Demo1: dimensionality reduction or feature extraction
• % ---------------------------------------------------------------------%
• clc
• clear all
• close all
• addpath(genpath(pwd))
• % 4 circles
• load circledata
• %
• X = circledata;
• for i = 1:4
•     scatter(X(1+250*(i-1):250*i,1),X(1+250*(i-1):250*i,2))
•     hold on
• end
• % Parameters setting
• options.sigma = 5;   % kernel width
• options.dims  = 2;   % output dimension
• options.type  = 0;   % 0:dimensionality reduction or feature extraction
•                      % 1:fault detection
• options.beta  = 0.9; % corresponding probabilities (for ault detection)
• options.cpc  = 0.85; % Principal contribution rate (for ault detection)
• % Train KPCA model
• model = kpca_train(X,options);
• figure
• for i = 1:4
•     scatter(model.mappedX(1+250*(i-1):250*i,1), ...
•         model.mappedX(1+250*(i-1):250*i,2))
•     hold on
• end
  

/ F1 h" e3 j" x5 A- [; @0 C2 j（2）结果 （分别为原图和特征提取后的图）
. h0 L/ q0 s5 f0 O) c

4. demo2: 故障检测（需要调节核宽度、主元贡献率和置信度等参数来提高故障检测效果）
（1）源代码
- |7 B2 H( k# \5 a2 h" g

• % Demo2: Fault detection
• % X: training samples
• % Y: test samples
• % Improve the peRFormance of fault detection by adjusting parameters
• % 1. options.sigma = 16;   % kernel width
• % 2. options.beta          % corresponding probabilities
• % 3. options.cpc  ;        % principal contribution rate
• % ---------------------------------------------------------------------%
• clc
• clear all
• close all
• addpath(genpath(pwd))
• %
• X = rand(200,10);
• Y = rand(100,10);
• Y(20:40,: ) = rand(21,10)+3;
• Y(60:80,: ) = rand(21,10)*3;
• % Normalization (if necessary)
• % mu = mean(X);
• % st = std(X);
• % X = zscore(X);
• % Y = bsxfun(@rdivide,bsxfun(@minus,Y,mu),st);
• % Parameters setting
• options.sigma = 16;   % kernel width
• options.dims  = 2;   % output dimension
• options.type  = 1;   % 0:dimensionality reduction or feature extraction
•                      % 1:fault detection
• options.beta  = 0.9; % corresponding probabilities (for ault detection)
• options.cpc  = 0.85; % principal contribution rate (for ault detection)
• % Train KPCA model
• model = kpca_train(X,options);
• % Test a new sample Y (vector of matrix)
• [SPE,T2,mappedY] = kpca_test(model,Y);
• % Plot the result
• plotResult(model.SPE_limit,SPE);
• plotResult(model.T2_limit,T2);
  

（2）结果（分别是SPE统计量和T2统计量的结果图）
( [5 h+ w7 P2 l' S" c$ x

" q1 @9 g, J1 P# S) C附件是基于KPCA的降维、特征提取和故障检测程序源代码。如有错误的地方请指出，谢谢。

游客，如果您要查看本帖隐藏内容请回复

Demyar

看看楼主的代码。

szm2000

fanichicl

基于KPCA的降维、特征提取和故障检测程序源代码。