Home > Codi > VRE_fault_rec.m

VRE_fault_rec

PURPOSE ^

-------------------------------------------------------------------------

SYNOPSIS ^

function [f, xiMV, Xrec] = VRE_fault_rec(X, PCA_Struct)

DESCRIPTION ^

 -------------------------------------------------------------------------
 |            [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.                            |
 -------------------------------------------------------------------------

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Wed 12-Sep-2012 13:03:54 by m2html © 2005