From f0c9988d2ced1c6e7495a11292117d7d0411516b Mon Sep 17 00:00:00 2001 From: Erik Frisk Date: Tue, 19 Jun 2018 23:24:46 +0200 Subject: [PATCH] Slight cleanup of code and added SingleFaultIsolation.m to src directory --- README.md | 2 +- code/main.m | 154 +++++++++++++++++--------------- code/src/SingleFaultIsolation.m | 14 +++ 3 files changed, 97 insertions(+), 73 deletions(-) create mode 100644 code/src/SingleFaultIsolation.m diff --git a/README.md b/README.md index e0625ce..d6fae27 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ This repository contains data and code to generate results from the paper _Residual Selection for Consistency Based Diagnosis Using Machine Learning Models_ -by Erik Frisk and Mattias Krysander, Department of Electrical Engineering, Linköping Univeristy, Sweden \ No newline at end of file +by Erik Frisk and Mattias Krysander, Department of Electrical Engineering, Linköping University, Sweden diff --git a/code/main.m b/code/main.m index d86569f..844a08c 100644 --- a/code/main.m +++ b/code/main.m @@ -1,6 +1,6 @@ % This script produces the results shown in: % -% Erik Frisk and Mattias Krysander (2018). Residaul Selection for +% Erik Frisk and Mattias Krysander (2018). Residual Selection for % Consistency Based Diagnosis Using Machine Learning Models. In Safeprocess % 2018, Warsaw, Polen. % @@ -52,47 +52,48 @@ plotResiduals(absdata) clear absdata %% Plot Fig. 5: Probability of residual alarm given a given mode -no_residuals = size(data.res,2); -all_residuals =[1:no_residuals]; % select all residuals -plotDecisionMatrix(all_residuals,data); +no_residuals = size(data.res, 2); +all_residuals = 1:no_residuals; % select all residuals +plotDecisionMatrix(all_residuals, data); clear no_residuals ans %% Plot Fig 6: Structural isolability of all available residuals load model.mat; -figure; -model.IsolabilityAnalysisFSM(data.fsm(:,[2:8]),'permute',0); +figure(16); +model.IsolabilityAnalysisFSM(data.fsm(:, 2:8), 'permute', 0); title('Structural isolability matrix using all available tests'); -ylabel('Injected fault','FontWeight', 'bold'); -xlabel('Diagnsoed fault','FontWeight', 'bold'); +ylabel('Injected fault', 'FontWeight', 'bold'); +xlabel('Diagnsoed fault', 'FontWeight', 'bold'); clear ans + %% Create thresholded residuals thdata = data; -thdata.res = abs(data.res)>=1; +thdata.res = abs(data.res) >= 1; %% Plot Fig 7: Confusion matrix obtained by consistency-based diagnosis (CBD) using all thresholded residuals -figure; -CBDPlots(thdata,model,all_residuals); +figure(17); +CBDPlots(thdata, model, all_residuals); clear all_residuals ans %% Train a random forest of 100 classification trees using the all data. trees = 100; rng(1); % For reproducibility -t = templateTree('PredictorSelection','interaction-curvature'); +t = templateTree('PredictorSelection', 'interaction-curvature'); fprintf('Building random forest classifier ... '); tic; -Mdl = fitcensemble(single(thdata.res),thdata.mode,'Method','bag',... - 'NumLearningCycles',trees,'Learners',t); +Mdl = fitcensemble(single(thdata.res), thdata.mode, 'Method', 'bag',... + 'NumLearningCycles', trees, 'Learners', t); fprintf(' Finished in %.2f sec.\n', toc); clear t trees %% Plot Fig 8: Random forest classification performance using all tests on out-of-bag data. Y = oobPredict(Mdl); -Ival = confusionmat(categorical(thdata.mode),categorical(Y)); +Ival = confusionmat(categorical(thdata.mode), categorical(Y)); nf = 8; -C = Ival./(sum(Ival')'*ones(1,nf)); +C = Ival./(sum(Ival, 2)*ones(1, nf)); -figure; +figure(18); b = bar3(C); title('Confusion matrix - Validation data set') set(gca, 'XTickLabel', data.modes, 'TickLabelInterpreter', 'none') @@ -104,8 +105,9 @@ title('Fault Isolation Performance Matrix') for finjIdx=1:nf for fdiagIdx=1:nf - text(fdiagIdx,finjIdx,C(finjIdx,fdiagIdx)+0.05,... - sprintf('%.1f',C(finjIdx,fdiagIdx)*100),'HorizontalAlignment', 'center') + text(fdiagIdx, finjIdx, C(finjIdx, fdiagIdx) + 0.05,... + sprintf('%.1f', C(finjIdx,fdiagIdx)*100), ... + 'HorizontalAlignment', 'center') end end for k = 1:length(b) @@ -114,15 +116,15 @@ for k = 1:length(b) b(k).FaceColor = 'interp'; end colormap('Summer') -view(0,90) +view(0, 90) clear Y Ival nf C b fdiagIdx finjIdx k zdata thdata %% Plot Fig 9: Random forest out-of-bag classification error as a function of the number of trees. -figure; -oobErr = oobLoss(Mdl,'mode','cumulative'); -plot(oobErr,'LineWidth',2) -xlabel( 'Number of grown trees' , 'FontWeight', 'bold'); -ylabel( 'Out-of-bag classification error' , 'FontWeight', 'bold'); +figure(19); +oobErr = oobLoss(Mdl, 'mode', 'cumulative'); +plot(oobErr, 'LineWidth', 2) +xlabel('Number of grown trees', 'FontWeight', 'bold'); +ylabel('Out-of-bag classification error', 'FontWeight', 'bold'); clear oobErr %% Estimate residual importance by permuting out-of-bag observations. @@ -132,13 +134,14 @@ importance = oobPermutedPredictorImportance(Mdl); fprintf(' Finished in %.2f sec.\n', toc); %% Plot Fig 10: Residual importance -[imp,residualRanking] = sort(importance,'descend'); +[imp,residualRanking] = sort(importance, 'descend'); -figure; -plot(imp,'LineWidth',2) -xticks([1:length(residualRanking)]) +figure(20); +plot(imp,'LineWidth', 2) +xticks(1:length(residualRanking)) axis([0 43 -0.4 1.6]) -set(gca, 'XTickLabel', residualRanking, 'TickLabelInterpreter', 'none', 'FontSize',8) +set(gca, 'XTickLabel', residualRanking, 'TickLabelInterpreter', 'none', ... + 'FontSize', 8) title('Out-of-Bag Permuted Predictor Importance Estimates'); ylabel('Estimates','FontWeight', 'bold'); xlabel('Predictors','FontWeight', 'bold'); @@ -146,80 +149,87 @@ clear imp %% Performance evaluation of selected tests using consistency based diagnosis. -no_samples = size(data.res,1); % n +no_samples = size(data.res, 1); % n no_modes = numel(data.modes); % m no_of_samples_per_mode = no_samples/no_modes; % k -structuralIsolationMatrix = model.IsolabilityAnalysisFSM(data.fsm(:,[2:8]),'permute',0); +structuralIsolationMatrix = ... + model.IsolabilityAnalysisFSM(data.fsm(:, 2:8),'permute', 0); + +stim2 = sum(((2.^(0:6)).*structuralIsolationMatrix), 2); +probabilityMaximumIsolation = zeros(42, no_modes-1); -stim2 = sum(((2.^[0:6]).*structuralIsolationMatrix)'); -probabilityMaximumIsolation = zeros(42,no_modes-1); -C ={}; +C = cell(length(residualRanking), 1); +isolationerror = zeros(length(residualRanking), 1); +mean_md = zeros(length(residualRanking), 1); +fa = zeros(length(residualRanking), 1); for i = 1:length(residualRanking) - selected_residauls = residualRanking(1:i); % select the i best residuals - dx = SingleFaultIsolation(data.res(:,selected_residauls),... - data.fsm(selected_residauls,[2:8])); % apply CBD. + selected_residuals = residualRanking(1:i); % select the i best residuals + dx = SingleFaultIsolation(data.res(:, selected_residuals),... + data.fsm(selected_residuals, (2:8))); % apply CBD. % Compute for each mode the probability of maximum isolation - dx2 = sum(((2.^[0:6]).*dx(:,[2:end]))'); + dx2 = sum(((2.^(0:6)).*dx(:, (2:end))), 2); for j = 1:no_modes-1 % not NF - probabilityMaximumIsolation(i,j) = sum(dx2([1:no_of_samples_per_mode]+... + probabilityMaximumIsolation(i,j) = sum(dx2((1:no_of_samples_per_mode)+... j*no_of_samples_per_mode)==stim2(j))/no_of_samples_per_mode; end % Compute confusion matrix of CBD using the i best residuals - C{i} = CBDPlots(data,model,selected_residauls); + C{i} = CBDPlots(data, model, selected_residuals); % Compute performance indicators isolationerror(i) = sum(sum(abs(structuralIsolationMatrix-C{i}(2:8,2:8)./((1-C{i}(2:end,1))*ones(1,7)))))/7^2; - mean_md(i) = mean(C{i}(2:end,1)); - fa(i) = 1-C{i}(1,1); + mean_md(i) = mean(C{i}(2:end, 1)); + fa(i) = 1-C{i}(1, 1); end -clear stim2 dx dx2 no_samples no_modes no_of_samples_per_mode selected_residauls C +clear stim2 dx dx2 no_samples no_modes no_of_samples_per_mode selected_residuals C clear structuralIsolationMatrix i j %% Plot Fig 11: False alarm probability, mean missed detection probability, % and aggregated isolation error as a function of selected tests. -figure; -plot([mean_md;fa;isolationerror]','LineWidth',2) +figure(21); +plot([mean_md, fa, isolationerror], 'LineWidth', 2) hold -selected_no_of_residauls = [12 14 26 27]; -plot(selected_no_of_residauls,isolationerror(selected_no_of_residauls),'x',... - 'LineWidth',2,'MarkerSize',8,'Color','k') +selected_no_of_residuals = [12 14 26 27]; +plot(selected_no_of_residuals,isolationerror(selected_no_of_residuals), ... + 'x', 'LineWidth', 2, 'MarkerSize', 8, 'Color', 'k') -legend('Missed detection probability','False alarm probability','Fault isolation errors') -xlabel('Number of selected residauls','FontWeight', 'bold') -ylabel('Probability','FontWeight', 'bold') +legend('Missed detection probability', 'False alarm probability', ... + 'Fault isolation errors') +xlabel('Number of selected residuals','FontWeight', 'bold') +ylabel('Probability', 'FontWeight', 'bold') grid on -clear mean_md isolationerror selected_no_of_residauls +clear mean_md isolationerror selected_no_of_residuals %% Plot Fig 12: Consistency based isolation performance per mode as a function of number of selected residuals. -figure; -plot([1-fa; probabilityMaximumIsolation']','LineWidth',2) +figure(22); +plot([1-fa, probabilityMaximumIsolation], 'LineWidth', 2) l = legend(data.modes{1:end}); set(l, 'Interpreter', 'none') -xlabel('Number of selected residuals','FontWeight', 'bold') -ylabel('Probability of maximum isolation performance','FontWeight', 'bold') +xlabel('Number of selected residuals', 'FontWeight', 'bold') +ylabel('Probability of maximum isolation performance', 'FontWeight', 'bold') grid on clear probabilityMaximumIsolation fa l %% Plot Fig 13: Consistency-based diagnosis results obtained when running different numbers k of most important residuals. % Plots the cases when using 10, 12, 26, and 27 residuals: -selected_no_of_residauls = [10 12 26 27]; -figure; -C ={}; -for i = 1:length(selected_no_of_residauls) - subplot(2,2,i); - selected_residauls = residualRanking(1:selected_no_of_residauls(i)); - C{i} = CBDPlots(data,model,selected_residauls); - title(['No of tests: ' num2str(selected_no_of_residauls(i))]); +selected_no_of_residuals = [10 12 26 27]; +figure(23); +C = cell(length(selected_no_of_residuals), 1); +for i = 1:length(selected_no_of_residuals) + subplot(2, 2, i); + selected_residuals = residualRanking(1:selected_no_of_residuals(i)); + C{i} = CBDPlots(data, model, selected_residuals); + title(['No of tests: ' num2str(selected_no_of_residuals(i))]); end clear C i %% Plot structural isolability using the k best residuals. Compare with previous plot. -figure; -for i = 1:length(selected_no_of_residauls) - subplot(2,2,i) - selected_residauls = residualRanking(1:selected_no_of_residauls(i)); - model.IsolabilityAnalysisFSM(data.fsm(selected_residauls,[2:8]),'permute',0); % full str isolability - title(['No of tests: ' num2str(selected_no_of_residauls(i))]) +figure(24); +for i = 1:length(selected_no_of_residuals) + subplot(2, 2, i) + selected_residuals = residualRanking(1:selected_no_of_residuals(i)); + model.IsolabilityAnalysisFSM(data.fsm(selected_residuals, (2:8)), ... + 'permute', 0); % full str isolability + title(['No of tests: ' num2str(selected_no_of_residuals(i))]) end -%clear selected_no_of_residauls selected_residauls residualRanking i ans +%clear selected_no_of_residuals selected_residuals residualRanking i ans diff --git a/code/src/SingleFaultIsolation.m b/code/src/SingleFaultIsolation.m new file mode 100644 index 0000000..7673f5a --- /dev/null +++ b/code/src/SingleFaultIsolation.m @@ -0,0 +1,14 @@ +function isol=SingleFaultIsolation( res, FSM ) + N = size(res,1); + nf = size(FSM,2); + + isol = zeros(N,nf + 1); + for k=1:N + alarms = abs(res(k,:))>=1; + if any(alarms) + isol(k,:) = [false all(FSM(alarms,:),1)]; + else + isol(k,:) = [true boolean(zeros(1,nf))]; + end + end +end \ No newline at end of file -- 2.22.0