------------------------------------------------------------------------- | [f, xiMV, Xrec] = VRE_fault_rec(X, PCA_Struct) | ------------------------------------------------------------------------- ----------------- | FUNCTIONALITY | ------------------------------------------------------------------------- | This function aim is to first identify the fault direction affecting | | each faulty observation. Then, the method computes the fault | | magnitude in every variable presenting the abnormal behaviour, and | | finally, reconstructing each faulty observation. To do so, the method | | employs the reconstruction-based contributions (Alcala and Qin, 2009) | | to identify the variables with an abnormal behaviour. Then, the | | fault magnitude and fault reconstruction are computed using the | | formulae in (Dunia and Qin, 1996). | ------------------------------------------------------------------------- -------------------- | INPUT PARAMETERS | ------------------------------------------------------------------------- | X: Array or matrix of double values whose observations | | (rows) have to be reconstructed (if their square | | prediction error surpasses the statistical limit of the | | PCA model). | | PCA_Struct: Structure that contains the information required to | | project new observations to a given principal component | | analysis (PCA) model. This structure must have, at least: | | -> mx: array containing the mean for each variable | | (column) to subtract to X. | | -> std: array containing the standard deviation of each | | variables (column) in X. This value is used to rescalate | | (divide) each variable. | | -> P: matrix containing the weights (loadings) to project | | each observation (row) in X to the PCA model. | | -> SPE_lim: statistical limit for the SPE. This value is | | used to determine if an observation (row) in X is faulty | | (greater than SPE_lim) or not (equal or lower). | ------------------------------------------------------------------------- --------------------- | OUTPUT PARAMETERS | ------------------------------------------------------------------------- | f: Array containing the fault magnitude for each faulty | | observation in X. A non-faulty observation has f set to 0, | | while any value different value indicates that the | | observation is faulty. | | xiMV: Array or matrix (if X was a matrix) used to reconstruct each | | observation (row) in X. Each row in xiMV is the fault | | direction that affects the corresponding observation in X. A | | non-faulty observation has all elements set to 0, while a | | faulty observation has the corresponding elements affected by | | the fault set to a value different to 0. This value | | represents the unitary vector representing the projection of | | the fault for the variable (column) in question. | | Xrec: Reconstructed values for the faulty observations that | | surpassed the SPE limit. Those observations that were not | | considered as faulty have the same values, while the faulty | | ones present only variations on the variables affected by the | | the corresponding fault direction. | -------------------------------------------------------------------------
0001 % ------------------------------------------------------------------------- 0002 % | [f, xiMV, Xrec] = VRE_fault_rec(X, PCA_Struct) | 0003 % ------------------------------------------------------------------------- 0004 % ----------------- 0005 % | FUNCTIONALITY | 0006 % ------------------------------------------------------------------------- 0007 % | This function aim is to first identify the fault direction affecting | 0008 % | each faulty observation. Then, the method computes the fault | 0009 % | magnitude in every variable presenting the abnormal behaviour, and | 0010 % | finally, reconstructing each faulty observation. To do so, the method | 0011 % | employs the reconstruction-based contributions (Alcala and Qin, 2009) | 0012 % | to identify the variables with an abnormal behaviour. Then, the | 0013 % | fault magnitude and fault reconstruction are computed using the | 0014 % | formulae in (Dunia and Qin, 1996). | 0015 % ------------------------------------------------------------------------- 0016 % -------------------- 0017 % | INPUT PARAMETERS | 0018 % ------------------------------------------------------------------------- 0019 % | X: Array or matrix of double values whose observations | 0020 % | (rows) have to be reconstructed (if their square | 0021 % | prediction error surpasses the statistical limit of the | 0022 % | PCA model). | 0023 % | PCA_Struct: Structure that contains the information required to | 0024 % | project new observations to a given principal component | 0025 % | analysis (PCA) model. This structure must have, at least: | 0026 % | -> mx: array containing the mean for each variable | 0027 % | (column) to subtract to X. | 0028 % | -> std: array containing the standard deviation of each | 0029 % | variables (column) in X. This value is used to rescalate | 0030 % | (divide) each variable. | 0031 % | -> P: matrix containing the weights (loadings) to project | 0032 % | each observation (row) in X to the PCA model. | 0033 % | -> SPE_lim: statistical limit for the SPE. This value is | 0034 % | used to determine if an observation (row) in X is faulty | 0035 % | (greater than SPE_lim) or not (equal or lower). | 0036 % ------------------------------------------------------------------------- 0037 % --------------------- 0038 % | OUTPUT PARAMETERS | 0039 % ------------------------------------------------------------------------- 0040 % | f: Array containing the fault magnitude for each faulty | 0041 % | observation in X. A non-faulty observation has f set to 0, | 0042 % | while any value different value indicates that the | 0043 % | observation is faulty. | 0044 % | xiMV: Array or matrix (if X was a matrix) used to reconstruct each | 0045 % | observation (row) in X. Each row in xiMV is the fault | 0046 % | direction that affects the corresponding observation in X. A | 0047 % | non-faulty observation has all elements set to 0, while a | 0048 % | faulty observation has the corresponding elements affected by | 0049 % | the fault set to a value different to 0. This value | 0050 % | represents the unitary vector representing the projection of | 0051 % | the fault for the variable (column) in question. | 0052 % | Xrec: Reconstructed values for the faulty observations that | 0053 % | surpassed the SPE limit. Those observations that were not | 0054 % | considered as faulty have the same values, while the faulty | 0055 % | ones present only variations on the variables affected by the | 0056 % | the corresponding fault direction. | 0057 % ------------------------------------------------------------------------- 0058 function [f, xiMV, Xrec] = VRE_fault_rec(X, PCA_Struct) 0059 0060 %% PARAMETER CHECK 0061 0062 % First of all, we check if the user passed all required arguments to 0063 % proceed with the function. 0064 if nargin ~= 2 0065 % - We inform the user that we need 2 input parameters. 0066 error('VRE_fault_rec:notEnoughInputs', ... 0067 'ERROR: The method requires two input parameters!'); 0068 end 0069 0070 % Next, we check if the first input (X) is a numeric array or matrix 0071 if ~isnumeric(X) 0072 % - We inform the user that X must be a numeric array or matrix. 0073 error('VRE_fault_rec:wrongType', ... 0074 'ERROR: X must be a numeric array or matrix!'); 0075 end 0076 0077 % Finally, we check if the second input parameter is a structure... 0078 if ~isstruct(PCA_Struct) 0079 % - We inform the user that PCA_Struct must be a Matlab structure 0080 % object. 0081 error('VRE_fault_rec:wrongType', ... 0082 'ERROR: PCA_Struct must be a structure!'); 0083 else 0084 % - If PCA_Struct is a structure, it has to have somespecific fields: 0085 % -- Mean (mx): 0086 if ~isfield(PCA_Struct, 'mx') 0087 error('VRE_fault_rec:missingField', ... 0088 'ERROR: PCA_Struct must have the field mx!'); 0089 end 0090 % -- Standard deviation (std): 0091 if ~isfield(PCA_Struct, 'std') 0092 error('VRE_fault_rec:missingField', ... 0093 'ERROR: PCA_Struct must have the field std!'); 0094 end 0095 % -- Loading matrix (P): 0096 if ~isfield(PCA_Struct, 'P') 0097 error('VRE_fault_rec:missingField', ... 0098 'ERROR: PCA_Struct must have the field P!'); 0099 end 0100 % -- Square prediction limit (SPE_lim): 0101 if ~isfield(PCA_Struct, 'SPE_lim') 0102 error('VRE_fault_rec:missingField', ... 0103 'ERROR: PCA_Struct must have the field SPE_lim!'); 0104 end 0105 end 0106 0107 %% INITIALISATIONS. 0108 0109 % The first step will consist on auto scaling the given data matrix (X) 0110 % using the mean and standard deviations obtained to build the PCA model. 0111 X_as = standardise_to(X, PCA_Struct.mx, PCA_Struct.std); 0112 0113 % Next, we initialise all output parameters: 0114 % - Fault magnitude (f): 0115 f = zeros(size(X, 1), 1); 0116 % - Fault direction to use when reconstructing each faulty observation 0117 % (xiMV): 0118 xiMV = zeros(size(X)); 0119 % - Values of the reconstructed faulty observations (Xrec): 0120 Xrec = X; 0121 0122 % With the standardised data, we project it into the projection space of 0123 % the PCA model (T): 0124 T = X_as * PCA_Struct.P; 0125 0126 % We compute the projection error (E): 0127 E = X_as - T * PCA_Struct.P'; 0128 0129 % And the values of the squared prediction error (SPE): 0130 SPE = sum(E .^ 2, 2); 0131 0132 % Next, we compare the SPE values with the statistical limit for this 0133 % index in order to identify the faulty observations. 0134 SPE_fd = SPE > PCA_Struct.SPE_lim; 0135 0136 % Next, we create the fault direction matrix. In this case, we are only 0137 % interested in the reconstruction of individual sensor faults, therefore, 0138 % the resulting fault direction matrix is the identity matrix. 0139 Xi = eye(size(X, 2)); 0140 0141 % Next, we project the fault direction matrix into the residual subspace. 0142 Xi_rs = (eye(size(X, 2)) - PCA_Struct.P * PCA_Struct.P') * Xi; 0143 0144 % The last step before computing the fault magnitude of each faulty 0145 % observation is to compute the norm of each fault direction in the 0146 % residual subpsace. 0147 % - To do so, we initialise the vector to store the norm 0148 % of each direction: 0149 normXi_rs = zeros(1, size(Xi_rs, 2)); 0150 % - And for each fault direction... 0151 for i = 1:length(normXi_rs) 0152 % -- We compute the norm of the i-th fault direction. 0153 normXi_rs(i) = norm(Xi_rs(:, i)); 0154 end 0155 0156 %% FAULT MAGNITUDE AND INIDIVIDUAL VARIABLE FAULT IDENTIFICATION. 0157 0158 % Using all the information we have from the previous steps (fault 0159 % directions projected into the residual subspace, residual matrix of the 0160 % observations and norm of each fault direction) we compute the fault 0161 % magnitude for each faulty observation. 0162 % - First, we initialise the matrix that will contain the values of the 0163 % fault projected in each possible fault direction. 0164 fMat = zeros(size(X)); 0165 % - And then, we compute the fault magnitude for each faulty observation. 0166 fMat(SPE_fd, :) = (Xi_rs' * X_as(SPE_fd, :)')' ./ ... 0167 repmat(normXi_rs .^ 2, sum(SPE_fd), 1); 0168 0169 % We also compute the reconstruction-based contributions in order to grant 0170 % that at least individual variable faults can be correctly diagnosed. 0171 RBC_Struct = compute_RBC(X_as, PCA_Struct); 0172 0173 % Next, the first fault direction candidate is the one that presents a 0174 % higher RBC to the SPE index. 0175 [~, SPE_RBC_sorted] = sort(RBC_Struct.SPE.RBC, 2, 'descend'); 0176 0177 % We reconstruct the fault directions using the fault direction associated 0178 % to the maximum RBC value for each faulty observation. 0179 % - First, we create the variable to store the different fault directions 0180 % to used when reconstructing... 0181 Xi_rec1 = zeros(size(X)); 0182 % - And now, we indicate the variables to reconstruct (only one at a time) 0183 % for each faulty observation. 0184 Xi_rec1(sub2ind(size(Xi_rec1), (1:size(X, 1))', SPE_RBC_sorted(:, 1))) = 1; 0185 % - Next, we indicate the fault magnitude of the fault projected in the 0186 % first candidate direction. 0187 f_rec1 = fMat(sub2ind(size(fMat), (1:size(X, 1))', SPE_RBC_sorted(:, 1))); 0188 % - And with these two values, we proceed to the reconstruction of the 0189 % faulty observations (considering only individual variable faults). 0190 X_rec_as1 = X_as - (Xi_rec1 .* repmat(f_rec1, 1, size(X, 2))); 0191 0192 % Next, we project the reconstructed faulty observations to the calibration 0193 % model. 0194 T_rec1 = X_rec_as1 * PCA_Struct.P; 0195 0196 % We compute the residual matrix for the reconstructed observations. 0197 E_rec1 = X_rec_as1 - T_rec1 * PCA_Struct.P'; 0198 0199 % And we compute the SPE for the reconstructed observations. 0200 SPE_rec1 = sum(E_rec1 .^ 2, 2); 0201 0202 % Next, we compare the new values of the SPE with the SPE statistical 0203 % control limit in order to find those observations whose misbehaviour has 0204 % been corrected. 0205 SPE_rec1_fd = SPE_rec1 > PCA_Struct.SPE_lim; 0206 0207 % Those observations whose SPE value after being reconstructed are below 0208 % the statistical control limit have been correctly reconstructed. 0209 % Therefore, we can copy these values to the corresponding variables. 0210 % - First, we locate the positions to overwrite as those that were 0211 % initially over the SPE limit and once reconstructed are below the SPE 0212 % limit. 0213 overwrite = xor(SPE_fd, SPE_rec1_fd); 0214 % - Now, we desnormalise the values of the reconstructed faulty 0215 % observations. 0216 Xrec(overwrite, :) = X_rec_as1(overwrite, :) .* ... 0217 repmat(PCA_Struct.std, sum(overwrite), 1) + ... 0218 repmat(PCA_Struct.mx, sum(overwrite), 1); 0219 0220 % We also update the values of the fault magnitude to reconstruct those 0221 % positions: 0222 f(overwrite) = f_rec1(overwrite); 0223 0224 % And the fault direction to reconstruct those observations. 0225 xiMV(overwrite, :) = Xi_rec1(overwrite, :); 0226 0227 %% MULTIVARIABLE FAULT TREATEMENT. 0228 0229 % In this part of the code we will deal with multivariable faults ocurring 0230 % on the studied process. In the previous section we reconstructed the 0231 % individual variable faults, so in this section we will focus on those 0232 % observations for which the SPE value surpasses the SPE limit even after 0233 % reconstructing one of the variables. 0234 0235 % The first step is to know how faults propagate from one variable to the 0236 % rest in the residual space (since we are detecting faults using the SPE). 0237 % Initially, we compute this value considering a unit fault, and then we 0238 % will compute the expected value considering the fault magnitude estimated 0239 % in the initial fault direction reconstruction. 0240 % - The unit RBC propagation of a fault depends solely on the relation 0241 % between the residual covariance between each pair of variables and its 0242 % own residual variance. So first, we compute the residual covariance 0243 % matrix. 0244 Crs = eye(size(X, 2)) - PCA_Struct.P * PCA_Struct.P'; 0245 % - Now, we can compute the RBC propagation of a unit fault magnitude for 0246 % each pair of variables. 0247 cProp = (Crs .^ 2) ./ repmat(diag(Crs)', size(X, 2), 1); 0248 % - Next, we compute how we had expected the contributions of the faulty 0249 % observation should an individual sensor fault had happened. 0250 expectedProp = cProp(SPE_RBC_sorted(:, 1), :) .* ... 0251 (repmat(fMat(sub2ind(size(fMat), (1:size(fMat, 1))', ... 0252 SPE_RBC_sorted(:, 1))), 1, size(cProp, 2)) .^ 2); 0253 % - And now, we compute the difference between the expected value and the 0254 % real value we've found. 0255 difProp = RBC_Struct.SPE.RBC - expectedProp; 0256 % - Finally, we sort the differences from higher (the next candidate) to 0257 % lower (the one we've reconstructed). 0258 [~, mvCandidatesSorted] = sort(difProp, 2, 'descend'); 0259 0260 % Now comes the iterative part to identify the fault affecting each 0261 % observation. 0262 % - First, we initialise the variable that will control the loop. 0263 loop = sum(SPE_rec1_fd) > 0 ; 0264 % - Next, we initialise the variable that indicates how many variables have 0265 % to be used to reconstruct each faulty observation... 0266 up2 = 1; 0267 % - We also create a variable to store the previous values of the SPE in 0268 % order to determine which observations can be reconstructed in the current 0269 % iteration. 0270 SPE_ant_fd = SPE_rec1_fd; 0271 % - We also store the values of the previously reconstructed values (in the 0272 % auto scaled space). 0273 X_rec_as_ant = X_rec_as1; 0274 % - While we have to loop... 0275 while loop 0276 % -- We initialise the fault direction matrix to reconstruct the faulty 0277 % observations. 0278 Xi_rec_i = zeros(size(X)); 0279 0280 % -- Next, we form the list of the candidates to use when 0281 % reconstructing each faulty observation. This uses the first candidate 0282 % and the set of variables with the highest difference between the 0283 % expected value and the real RBC. 0284 % --- To do so, first, we copy the fault magnitude values used for the 0285 % first reconstruction. 0286 Xi_rec_i(sub2ind(size(Xi_rec_i), find(SPE_ant_fd), ... 0287 SPE_RBC_sorted(SPE_ant_fd, 1))) = ... 0288 fMat(sub2ind(size(Xi_rec_i), find(SPE_ant_fd), ... 0289 SPE_RBC_sorted(SPE_ant_fd, 1))); 0290 % --- Next, we find the positions of the faulty observations. 0291 fPos_i = find(SPE_ant_fd); 0292 % --- Now, for each faulty observation... 0293 for i = 1:length(fPos_i) 0294 % ---- We copy all the values of the fault direction candidates for 0295 % the i-th faulty observation. 0296 Xi_rec_i(fPos_i(i), mvCandidatesSorted(fPos_i(i), 1:up2)) = ... 0297 fMat(fPos_i(i), mvCandidatesSorted(fPos_i(i), 1:up2)); 0298 end 0299 % % --- Next, we add the rest of reconstruction candidates, whose 0300 % % difference between the expected value and the current RBC values. 0301 % Xi_rec_i(sub2ind(size(Xi_rec_i), find(SPE_ant_fd), ... 0302 % mvCandidatesSorted(SPE_ant_fd, 1:up2))) = ... 0303 % fMat(sub2ind(size(Xi_rec_i), find(SPE_ant_fd), ... 0304 % mvCandidatesSorted(SPE_ant_fd, 1:up2))); 0305 0306 % -- We compute the norm of each row, which is equivalent to computing 0307 % the global fault magnitude affecting the observations. 0308 % --- First of all, we initialise the variable to ones, since we will 0309 % be dividing all values by this value. 0310 f_i = ones(size(X, 1), 1); 0311 % --- Next, we compute the norm of all faulty observations in the 0312 % previous iteration. 0313 f_i(SPE_ant_fd) = sqrt(sum(Xi_rec_i(SPE_ant_fd, :) .^ 2, 2)); 0314 0315 % -- Finally, we divide each observation's fault direction by the fault 0316 % magnitude affecting the whole process. 0317 Xi_rec_i = Xi_rec_i ./ repmat(f_i, 1, size(Xi_rec_i, 2)); 0318 0319 % -- Next, we proceed to reconstruct the faulty observations using both 0320 % the multivariable fault directions and the global fault magnitude. 0321 % --- First, we reconstructed values in the previous observation... 0322 X_rec_as_i = X_rec_as_ant; 0323 % --- Now, we reconstruct the multivariate faulty observations. 0324 X_rec_as_i(SPE_ant_fd, :) = X_rec_as_i(SPE_ant_fd, :) - ... 0325 (Xi_rec_i(SPE_ant_fd, :) .* repmat(f_i(SPE_ant_fd), 1, ... 0326 size(X, 2))); 0327 % --- We project the reconstructed observations into the PCA model. 0328 T_i = X_rec_as_i * PCA_Struct.P; 0329 % --- We also compute the residual matrix of the observations. 0330 E_i = X_rec_as_i - T_i * PCA_Struct.P'; 0331 % --- And the SPE values. 0332 SPE_i = sum(E_i .^ 2, 2); 0333 % --- Which observations are still faulty? 0334 SPE_i_fd = SPE_i > PCA_Struct.SPE_lim; 0335 0336 % -- Next, we look for which observations have been correctly 0337 % reconstructed in this iteration. 0338 overwrite_i = xor(SPE_ant_fd, SPE_i_fd); 0339 0340 % -- And we update the values of both the fault magnitude and fault 0341 % direction for the correctly reconstructed observations in this loop. 0342 % --- Fault magnitude. 0343 f(overwrite_i) = f_i(overwrite_i); 0344 % --- Fault direction to use when reconstructing. 0345 xiMV(overwrite_i, :) = Xi_rec_i(overwrite_i, :); 0346 % --- Reconstructed observations denormalised. 0347 Xrec(overwrite_i, :) = X_rec_as_i(overwrite_i, :) .* ... 0348 repmat(PCA_Struct.std, sum(overwrite_i), 1) + ... 0349 repmat(PCA_Struct.mx, sum(overwrite_i), 1); 0350 0351 % -- Finally, we check if we have to loop (there are still some 0352 % observations that have to be reconstructed) or not (all faulty 0353 % observations were correctly reconstructed), or the observations 0354 % cannot be reconstructed (all variables used to reconstruct). 0355 loop = sum(SPE_i_fd) > 0 && up2 < (size(xiMV, 2) - 1); 0356 0357 % -- If we have to loop another time... 0358 if loop 0359 % --- Then, we have to update some variables for the next 0360 % iteration. 0361 % ---- Position to use for the next reconstruction step. 0362 up2 = up2 + 1; 0363 % ---- Faulty observations found on the previous iteration after 0364 % reconstruction. 0365 SPE_ant_fd = SPE_i_fd; 0366 % ---- Previously reconstructed values (on the auto scaled space). 0367 X_rec_as_ant = X_rec_as_i; 0368 end 0369 end