TA的每日心情 | 开心
2019-11-20 15:05 |
|---|
签到天数: 2 天 [LV.1]初来乍到
|
EDA365欢迎您登录！
您需要 登录 才可以下载或查看，没有帐号？注册
x
本帖最后由 Colbie 于 2020-3-18 09:58 编辑 9 ~: ^; t; \, G6 X! C1 r x
: i& n7 j/ _1 s9 A ]3 x, d; q
核主元分析（Kernel principal component analysis ，KPCA）在降维、特征提取以及故障检测中的应用。主要功能有：（1）训练数据和测试数据的非线性主元提取（降维、特征提取）
6 l& _, t& z- A9 ~（2）SPE和T2统计量及其控制限的计算
9 f4 I0 x/ t; O0 ^: G( j（3）故障检测& |3 q: n7 B) H$ E: e0 m+ m
1 F1 g+ t, T" A5 }! Y8 D* e参考文献：) E( Z8 g+ q K4 X' _% ^$ u
Lee 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.
6 B* p! I6 I- |1 n2 z# s: q: A5 q* _; h& J! { T5 `, U0 \0 e
1. KPCA的建模过程（故障检测）：* G& W5 y/ y: |' a P G
（1）获取训练数据（工业过程数据需要进行标准化处理）
" }9 q' ~: z" b+ O. H) n（2）计算核矩阵
' u& k# ~; @' u/ Q, K6 ~! G（3）核矩阵中心化
" B# h. W8 G: h: Z3 }# x（4）特征值分解
0 {* v" r/ c. O（5）特征向量的标准化处理1 l% ?0 i$ M8 m7 q; z
（6）主元个数的选取+ ^7 M9 F' q5 p% H/ l4 M8 P+ i
（7）计算非线性主成分（即降维结果或者特征提取结果）3 A3 o& S" v V/ ]: l' F" @
（8）SPE和T2统计量的控制限计算
0 n1 y. _3 }* R" X- 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
" Z* d) o6 z9 f" S2 ^) j0 ~% ^: @
5 n6 s0 k8 a {1 I1 l
0 A9 @% C; q# e' v
2. KPCA的测试过程：
: H4 ~) g; E! |2 q$ P- [（1）获取测试数据（工业过程数据需要利用训练数据的均值和标准差进行标准化处理）
- C1 S# [# g4 z! I0 Q（2）计算核矩阵+ ^2 _; n( V3 U& w" j9 j9 A# @
（3）核矩阵中心化& m7 |0 t9 l+ Y/ W" k
（4）计算非线性主成分（即降维结果或者特征提取结果）
0 T) [: m/ w5 ^" X（5）SPE和T2统计量的计算
) q) Q$ w+ }& @$ Q2 A9 W- 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
3 j7 V' R0 e! J
3 l6 S. g( z- i3 X8 @3 s" m! Q+ F$ O4 O' @
6 \! D! C. ^2 V* {* q; P" L5 N3 n
3. demo1: 降维、特征提取
) x7 P; D& @! X) I( B(1) 源代码7 d4 {9 m/ y' a; w% K+ {
- % 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
* n; f& ?' x4 f0 r
( ^8 |/ x# r% s& K! ~$ {
, x$ ^. {0 Q2 b0 J（2）结果 （分别为原图和特征提取后的图）) O4 y, F% x/ d: Z) {! _+ B
2 P D' R7 m* p" [8 Y
8 H% r& P. _' t/ p2 W0 q8 W+ T4. demo2: 故障检测（需要调节核宽度、主元贡献率和置信度等参数来提高故障检测效果）
$ H) t2 _! K$ p2 ]* j# S9 \（1）源代码9 j) Q1 _. m( {; k& T: c
- % 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);
$ Z9 T, i/ l9 w- I4 G, V+ J% T
! c/ r+ [$ }! O; f$ k- Y5 \
. T3 O6 S% F v0 c（2）结果（分别是SPE统计量和T2统计量的结果图）
) i% O0 t9 `$ k/ V& M6 H
4 R$ O" I0 [. [! R2 j- y4 b" n B; Q
v" Y* d" Z' y7 [附件是基于KPCA的降维、特征提取和故障检测程序源代码。如有错误的地方请指出，谢谢。# _% k! ?2 z- y3 A" f
0 \6 b0 {4 s) k- d |
|
/1 
发表于 2020-3-18 09:56
收藏
淘帖
支持！
反对！