clear
clc
% close all;
%% PARAMETERS
% variable selection
% s=[1,2,3,4,5,6,9,10,11,13,14,16,18,19,21,22];
s=[1:22,42:52];
% s=[1:52];
% time-lag
q = 2;
% number of dynamic factors
a = 6;
% number of AR model
p = 4;
% confidence limit
alpha = 0.99;
%% OFFLINE MODELING
load d00_te.mat
X0 = d00_te(:, s);
% load N00.mat;
% X0 = N00;
% augmented matrix
[n, m] = size(X0);
% standardization
[Xref, Xmean, Xstd] = zscore(X0);
% structured dynamic pca
[R, P, W, T] = sdpca(Xref, q, a);
n = size(T, 1);
for i = 1:n-p
tL = [];
for j = 1:p
tL = [tL, T(i+j-1, :)];
end
Ti(i, :) = tL;
end
To = T(p+1:end, :);
C = inv(Ti'*Ti)*Ti'*To;
Vref = To - Ti*C;
invSv = inv(Vref'*Vref/(n-p-1));
% pca on residual
Eref = Xref - T*P';
[resvec, ressco, resvals] = pca(Eref);
numPC = cpv(resvals, 0.85);
numPC = 16;
Pref = resvec(:, 1:numPC);
invS = inv(diag(resvals(1:numPC)));
% control limits
T2c_d = a*(n-1)/(n-a) * finv(alpha, a, n-a);
T2c_s = numPC*(n-1)/(n-numPC) * finv(alpha, numPC, n-numPC);
Fref = Eref - ressco(:, 1:numPC)*Pref';
Qref = diag(Fref * Fref');
Qmean = mean(Qref); Qvar = var(Qref);
g = Qvar/Qmean/2; h = 2*Qmean^2/Qvar;
Qc = g*chi2inv(alpha, h);
%% ONLINE MONITORING
load d10_te.mat
X1 = d10_te(:, s);
% load N04.mat;
% X1 = N04;
Xcrt = autoscale(X1, Xmean, Xstd);
Td = Xcrt*R;
n = size(Td, 1);
for i = 1:n-p
tL = [];
for j = 1:p
tL = [tL, Td(i+j-1, :)];
end
Ti_crt(i, :) = tL;
end
To_crt = Td(p+1:end, :);
Vcrt = To_crt - Ti_crt*C;
T2d = diag(Vcrt* invSv * Vcrt');
Ecrt = Xcrt - Td*P';
Tcrt = Ecrt*Pref;
T2s = diag(Tcrt*invS*Tcrt');
Fcrt = Ecrt - Tcrt*Pref';
Q = diag(Fcrt * Fcrt');
%% PLOTS
figure
N = size(T2d, 1);
plot(1:N, T2d, 'k-');
hold on
plot(1:N, repmat(T2c_d, N, 1), 'k-');
ylabel('T^2');
hold off
axis tight
figure
N = size(T2s, 1);
subplot(2,1,1)
plot(1:N, T2s, 'k-');
hold on
plot(1:N, repmat(T2c_s, N, 1), 'k-');
ylabel('T^2');
hold off
axis tight
subplot(2,1,2)
plot(1:N, Q, 'k-');
hold on
plot(1:N, repmat(Qc, N, 1), 'k-');
ylabel('Q');
xlabel('Samples')
hold off
axis tight
%% TYPE I/II ERROR
nf = 160;
% % Normal Data Testing
% k = 0;
% for i = 1:N-p
% if T2d(i) > T2c_d
% k = k+1;
% end
% end
% T2d_TypeI = 100*k/(N-p)
%
% k = 0;
% for i = 1:N
% if T2s(i) > T2c_s
% k = k+1;
% end
% end
% T2s_TypeI = 100*k/N
% k = 0;
% for i = 1:N
% if Q(i) > Qc
% k = k+1;
% end
% end
% Q_TypeI = 100*k/N
% Type II error (missing alarm rates)
k = 0;
for i = nf+1:N-p
if T2d(i) <= T2c_d
k = k+1;
end
end
T2d_TypeII = 100*k/(N-p-nf)
k = 0;
for i = nf+1:N
if T2s(i) <= T2c_s
k = k+1;
end
end
T2s_TypeII = 100*k/(N-nf)
k = 0;
for i = nf+1:N
if Q(i) <= Qc
k = k+1;
end
end
Q_TypeII = 100*k/(N-nf)
clc
% close all;
%% PARAMETERS
% variable selection
% s=[1,2,3,4,5,6,9,10,11,13,14,16,18,19,21,22];
s=[1:22,42:52];
% s=[1:52];
% time-lag
q = 2;
% number of dynamic factors
a = 6;
% number of AR model
p = 4;
% confidence limit
alpha = 0.99;
%% OFFLINE MODELING
load d00_te.mat
X0 = d00_te(:, s);
% load N00.mat;
% X0 = N00;
% augmented matrix
[n, m] = size(X0);
% standardization
[Xref, Xmean, Xstd] = zscore(X0);
% structured dynamic pca
[R, P, W, T] = sdpca(Xref, q, a);
n = size(T, 1);
for i = 1:n-p
tL = [];
for j = 1:p
tL = [tL, T(i+j-1, :)];
end
Ti(i, :) = tL;
end
To = T(p+1:end, :);
C = inv(Ti'*Ti)*Ti'*To;
Vref = To - Ti*C;
invSv = inv(Vref'*Vref/(n-p-1));
% pca on residual
Eref = Xref - T*P';
[resvec, ressco, resvals] = pca(Eref);
numPC = cpv(resvals, 0.85);
numPC = 16;
Pref = resvec(:, 1:numPC);
invS = inv(diag(resvals(1:numPC)));
% control limits
T2c_d = a*(n-1)/(n-a) * finv(alpha, a, n-a);
T2c_s = numPC*(n-1)/(n-numPC) * finv(alpha, numPC, n-numPC);
Fref = Eref - ressco(:, 1:numPC)*Pref';
Qref = diag(Fref * Fref');
Qmean = mean(Qref); Qvar = var(Qref);
g = Qvar/Qmean/2; h = 2*Qmean^2/Qvar;
Qc = g*chi2inv(alpha, h);
%% ONLINE MONITORING
load d10_te.mat
X1 = d10_te(:, s);
% load N04.mat;
% X1 = N04;
Xcrt = autoscale(X1, Xmean, Xstd);
Td = Xcrt*R;
n = size(Td, 1);
for i = 1:n-p
tL = [];
for j = 1:p
tL = [tL, Td(i+j-1, :)];
end
Ti_crt(i, :) = tL;
end
To_crt = Td(p+1:end, :);
Vcrt = To_crt - Ti_crt*C;
T2d = diag(Vcrt* invSv * Vcrt');
Ecrt = Xcrt - Td*P';
Tcrt = Ecrt*Pref;
T2s = diag(Tcrt*invS*Tcrt');
Fcrt = Ecrt - Tcrt*Pref';
Q = diag(Fcrt * Fcrt');
%% PLOTS
figure
N = size(T2d, 1);
plot(1:N, T2d, 'k-');
hold on
plot(1:N, repmat(T2c_d, N, 1), 'k-');
ylabel('T^2');
hold off
axis tight
figure
N = size(T2s, 1);
subplot(2,1,1)
plot(1:N, T2s, 'k-');
hold on
plot(1:N, repmat(T2c_s, N, 1), 'k-');
ylabel('T^2');
hold off
axis tight
subplot(2,1,2)
plot(1:N, Q, 'k-');
hold on
plot(1:N, repmat(Qc, N, 1), 'k-');
ylabel('Q');
xlabel('Samples')
hold off
axis tight
%% TYPE I/II ERROR
nf = 160;
% % Normal Data Testing
% k = 0;
% for i = 1:N-p
% if T2d(i) > T2c_d
% k = k+1;
% end
% end
% T2d_TypeI = 100*k/(N-p)
%
% k = 0;
% for i = 1:N
% if T2s(i) > T2c_s
% k = k+1;
% end
% end
% T2s_TypeI = 100*k/N
% k = 0;
% for i = 1:N
% if Q(i) > Qc
% k = k+1;
% end
% end
% Q_TypeI = 100*k/N
% Type II error (missing alarm rates)
k = 0;
for i = nf+1:N-p
if T2d(i) <= T2c_d
k = k+1;
end
end
T2d_TypeII = 100*k/(N-p-nf)
k = 0;
for i = nf+1:N
if T2s(i) <= T2c_s
k = k+1;
end
end
T2s_TypeII = 100*k/(N-nf)
k = 0;
for i = nf+1:N
if Q(i) <= Qc
k = k+1;
end
end
Q_TypeII = 100*k/(N-nf)