Commit 87134f51 authored by Erik Frisk's avatar Erik Frisk

Removed toolbox dependency and revised code slightly

parent 5619803b
% TODO % TODO
% %
% 1. Ta bort beroendet av toolbox
% 2. Gör så att alla celler kan köras (om föregående celler har körts)
% 3. Flytta viss kod till förklarande funktioner % 3. Flytta viss kod till förklarande funktioner
% 4. Licenstext?
% 5. Uppdatera readme % 5. Uppdatera readme
% 6. Kommentarer om version av matlab och figurer. Klura ut vilken version % 6. Kommentarer om version av matlab och figurer. Klura ut vilken version
% som använts i artikeln. % som använts i artikeln.
...@@ -15,7 +12,6 @@ ...@@ -15,7 +12,6 @@
% Consistency Based Diagnosis Using Machine Learning Models. In Safeprocess % Consistency Based Diagnosis Using Machine Learning Models. In Safeprocess
% 2018, Warsaw, Polen. % 2018, Warsaw, Polen.
% %
% Some parts of the code requires a path to the Fault Diagnosis Toolbox.
clc clc
clear clear
...@@ -30,14 +26,15 @@ addpath src ...@@ -30,14 +26,15 @@ addpath src
% pretreatment of the data has been performed: % pretreatment of the data has been performed:
% %
% 1. 100 samples has been removed after a mode switch to reduce the % 1. 100 samples has been removed after a mode switch to reduce the
% influence on mode transitions. % influence from mode transitions.
% 2. Data has been sorted according to mode. % 2. Data has been sorted according to mode.
% 3. Equal number of samples from each mode has been selected. % 3. Data has been balanced, i.e., equal number of samples from each mode
% has been selected.
% 4. A mean value of 1000 samples has been computed in order to reduce % 4. A mean value of 1000 samples has been computed in order to reduce
% noise an reduce the size of the data set. % noise and reduce the size of the data set.
% 5. Each residual has been normalized based on fault free data and data % 5. Each residual has been normalized based on fault free data and data
% from decoupled faults. Normalization has been performed by taking the % from decoupled faults. Normalization has been performed by taking the
% residual absolut value and then divide it by its 99.5 %-percentil. In % residual absolut value and then divide it by its 99.5 %-percentile. In
% this way 0.5% of the residual samples in the fault free and decoupled % this way 0.5% of the residual samples in the fault free and decoupled
% modes will be greater than one. % modes will be greater than one.
load data.mat load data.mat
...@@ -66,16 +63,14 @@ clear absdata ...@@ -66,16 +63,14 @@ clear absdata
no_residuals = size(data.res, 2); no_residuals = size(data.res, 2);
all_residuals = 1:no_residuals; % select all residuals all_residuals = 1:no_residuals; % select all residuals
plotDecisionMatrix(all_residuals, data); plotDecisionMatrix(all_residuals, data);
clear no_residuals ans clear no_residuals
%% Plot Fig 6: Structural isolability of all available residuals %% Plot Fig 6: Structural isolability of all available residuals
load model.mat;
figure(16); figure(16);
model.IsolabilityAnalysisFSM(data.fsm(:, 2:8), 'permute', 0); IsolabilityAnalysisFSM(data.fsm(:, 2:8), data.modes(2:8), 'permute', 0);
title('Structural isolability matrix using all available tests'); title('Structural isolability matrix using all available tests');
ylabel('Injected fault', 'FontWeight', 'bold'); ylabel('Injected fault', 'FontWeight', 'bold');
xlabel('Diagnsoed fault', 'FontWeight', 'bold'); xlabel('Diagnsoed fault', 'FontWeight', 'bold');
clear ans
%% Create thresholded residuals %% Create thresholded residuals
thdata = data; thdata = data;
...@@ -83,20 +78,19 @@ thdata.res = abs(data.res) >= 1; ...@@ -83,20 +78,19 @@ thdata.res = abs(data.res) >= 1;
%% Plot Fig 7: Confusion matrix obtained by consistency-based diagnosis (CBD) using all thresholded residuals %% Plot Fig 7: Confusion matrix obtained by consistency-based diagnosis (CBD) using all thresholded residuals
figure(17); figure(17);
CBDPlots(thdata, model, all_residuals); CBDPlots(thdata, all_residuals);
clear all_residuals ans
%% Train a random forest of 100 classification trees using the all data. %% Train a random forest of 100 classification trees using the all data.
trees = 100; num_trees = 100;
rng(1); % For reproducibility rng(1); % For reproducibility
t = templateTree('PredictorSelection', 'interaction-curvature'); t = templateTree('PredictorSelection', 'interaction-curvature');
fprintf('Building random forest classifier ... '); fprintf('Building random forest classifier ... ');
tic; tic;
Mdl = fitcensemble(single(thdata.res), thdata.mode, 'Method', 'bag',... Mdl = fitcensemble(single(thdata.res), thdata.mode, 'Method', 'bag',...
'NumLearningCycles', trees, 'Learners', t); 'NumLearningCycles', num_trees, 'Learners', t);
fprintf(' Finished in %.2f sec.\n', toc); fprintf(' Finished in %.2f sec.\n', toc);
clear t trees clear t num_trees
%% Plot Fig 8: Random forest classification performance using all tests on out-of-bag data. %% Plot Fig 8: Random forest classification performance using all tests on out-of-bag data.
Y = oobPredict(Mdl); Y = oobPredict(Mdl);
...@@ -128,7 +122,7 @@ for k = 1:length(b) ...@@ -128,7 +122,7 @@ for k = 1:length(b)
end end
colormap('Summer') colormap('Summer')
view(0, 90) view(0, 90)
clear Y Ival nf C b fdiagIdx finjIdx k zdata thdata clear Y Ival nf C b fdiagIdx finjIdx k zdata
%% Plot Fig 9: Random forest out-of-bag classification error as a function of the number of trees. %% Plot Fig 9: Random forest out-of-bag classification error as a function of the number of trees.
figure(19); figure(19);
...@@ -140,7 +134,7 @@ clear oobErr ...@@ -140,7 +134,7 @@ clear oobErr
%% Estimate residual importance by permuting out-of-bag observations. %% Estimate residual importance by permuting out-of-bag observations.
fprintf('Compute variable importance for the random forest (this may take a while) ... '); fprintf('Compute variable importance for the random forest (this may take a while) ... ');
tic; tic;
importance = oobPermutedPredictorImportance(Mdl); importance = oobPermutedPredictorImportance(Mdl);
fprintf(' Finished in %.2f sec.\n', toc); fprintf(' Finished in %.2f sec.\n', toc);
...@@ -165,7 +159,7 @@ no_modes = numel(data.modes); % m ...@@ -165,7 +159,7 @@ no_modes = numel(data.modes); % m
no_of_samples_per_mode = no_samples/no_modes; % k no_of_samples_per_mode = no_samples/no_modes; % k
structuralIsolationMatrix = ... structuralIsolationMatrix = ...
model.IsolabilityAnalysisFSM(data.fsm(:, 2:8),'permute', 0); IsolabilityAnalysisFSM(data.fsm(:, 2:8), data.modes(2:8), 'permute', 0);
stim2 = sum(((2.^(0:6)).*structuralIsolationMatrix), 2); stim2 = sum(((2.^(0:6)).*structuralIsolationMatrix), 2);
probabilityMaximumIsolation = zeros(42, no_modes-1); probabilityMaximumIsolation = zeros(42, no_modes-1);
...@@ -185,7 +179,7 @@ for i = 1:length(residualRanking) ...@@ -185,7 +179,7 @@ for i = 1:length(residualRanking)
j*no_of_samples_per_mode)==stim2(j))/no_of_samples_per_mode; j*no_of_samples_per_mode)==stim2(j))/no_of_samples_per_mode;
end end
% Compute confusion matrix of CBD using the i best residuals % Compute confusion matrix of CBD using the i best residuals
C{i} = CBDPlots(data, model, selected_residuals); C{i} = CBDPlots(data, selected_residuals);
% Compute performance indicators % 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; 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)); mean_md(i) = mean(C{i}(2:end, 1));
...@@ -198,7 +192,7 @@ clear structuralIsolationMatrix i j ...@@ -198,7 +192,7 @@ clear structuralIsolationMatrix i j
% and aggregated isolation error as a function of selected tests. % and aggregated isolation error as a function of selected tests.
figure(21); figure(21);
plot([mean_md, fa, isolationerror], 'LineWidth', 2) plot([mean_md, fa, isolationerror], 'LineWidth', 2)
hold hold on
selected_no_of_residuals = [12 14 26 27]; selected_no_of_residuals = [12 14 26 27];
plot(selected_no_of_residuals,isolationerror(selected_no_of_residuals), ... plot(selected_no_of_residuals,isolationerror(selected_no_of_residuals), ...
'x', 'LineWidth', 2, 'MarkerSize', 8, 'Color', 'k') 'x', 'LineWidth', 2, 'MarkerSize', 8, 'Color', 'k')
...@@ -208,7 +202,8 @@ legend('Missed detection probability', 'False alarm probability', ... ...@@ -208,7 +202,8 @@ legend('Missed detection probability', 'False alarm probability', ...
xlabel('Number of selected residuals','FontWeight', 'bold') xlabel('Number of selected residuals','FontWeight', 'bold')
ylabel('Probability', 'FontWeight', 'bold') ylabel('Probability', 'FontWeight', 'bold')
grid on grid on
clear mean_md isolationerror selected_no_of_residuals hold off
clear selected_no_of_residuals
%% Plot Fig 12: Consistency based isolation performance per mode as a function of number of selected residuals. %% Plot Fig 12: Consistency based isolation performance per mode as a function of number of selected residuals.
figure(22); figure(22);
...@@ -218,7 +213,6 @@ set(l, 'Interpreter', 'none') ...@@ -218,7 +213,6 @@ set(l, 'Interpreter', 'none')
xlabel('Number of selected residuals', 'FontWeight', 'bold') xlabel('Number of selected residuals', 'FontWeight', 'bold')
ylabel('Probability of maximum isolation performance', 'FontWeight', 'bold') ylabel('Probability of maximum isolation performance', 'FontWeight', 'bold')
grid on grid on
clear probabilityMaximumIsolation fa l
%% Plot Fig 13: Consistency-based diagnosis results obtained when running different numbers k of most important residuals. %% 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: % Plots the cases when using 10, 12, 26, and 27 residuals:
...@@ -228,7 +222,7 @@ C = cell(length(selected_no_of_residuals), 1); ...@@ -228,7 +222,7 @@ C = cell(length(selected_no_of_residuals), 1);
for i = 1:length(selected_no_of_residuals) for i = 1:length(selected_no_of_residuals)
subplot(2, 2, i); subplot(2, 2, i);
selected_residuals = residualRanking(1:selected_no_of_residuals(i)); selected_residuals = residualRanking(1:selected_no_of_residuals(i));
C{i} = CBDPlots(data, model, selected_residuals); C{i} = CBDPlots(data, selected_residuals);
title(['No of tests: ' num2str(selected_no_of_residuals(i))]); title(['No of tests: ' num2str(selected_no_of_residuals(i))]);
end end
clear C i clear C i
...@@ -238,9 +232,8 @@ figure(24); ...@@ -238,9 +232,8 @@ figure(24);
for i = 1:length(selected_no_of_residuals) for i = 1:length(selected_no_of_residuals)
subplot(2, 2, i) subplot(2, 2, i)
selected_residuals = residualRanking(1:selected_no_of_residuals(i)); selected_residuals = residualRanking(1:selected_no_of_residuals(i));
model.IsolabilityAnalysisFSM(data.fsm(selected_residuals, (2:8)), ... IsolabilityAnalysisFSM(data.fsm(selected_residuals, (2:8)), data.modes(2:8), ...
'permute', 0); % full str isolability 'permute', 0); % full str isolability
title(['No of tests: ' num2str(selected_no_of_residuals(i))]) title(['No of tests: ' num2str(selected_no_of_residuals(i))])
end end
%clear selected_no_of_residuals selected_residuals residualRanking i ans
function [C] = CBDPlots(data,model,resIdx,i) function C = CBDPlots(data, resIdx)
% CBDPlots Plot result from consistency based diagnosis based on residuals
%
% C = CBDPlots(data, resIdx)
%
% Inputs:
% data - residual data structure
% resIdx - Index to residuals to evaluate
%
% Outputs:
% C - Confusion matrix
%
FSM42 = data.fsm(:,[2:end]); FSM42 = data.fsm(:, 2:end);
%% Evaluate selected set of residuals %% Evaluate selected set of residuals
[dx] = EvaluateResiduals(data, FSM42, resIdx); [dx] = EvaluateResiduals(data, FSM42, resIdx);
%% Confusion matrix computation %% Confusion matrix computation
% C(i,j) = P(fj diagnos|fi injected fault) % C(i,j) = P(fj diagnos|fi injected fault)
nf = numel(data.modes); nf = numel(data.modes);
C = zeros(nf,nf); C = zeros(nf,nf);
for fi = 1:nf for fi = 1:nf
for fj = 1:nf for fj = 1:nf
fjIdx = (data.mode==fj-1); fjIdx = (data.mode==fj-1);
C(fj,fi) = sum(dx(fjIdx,fi))/sum(fjIdx); C(fj,fi) = sum(dx(fjIdx,fi))/sum(fjIdx);
end
end end
end
%% Plot confusion matrix in a 3D-bar plot %% Plot confusion matrix in a 3D-bar plot
if nargin==4 b=bar3(C);
figure(i); clf; set(gca, 'XTickLabel', data.modes, 'TickLabelInterpreter', 'none')
end set(gca, 'YTickLabel', data.modes, 'TickLabelInterpreter', 'none')
xlabel('Diagnosed fault', 'FontWeight', 'bold');
ylabel('Injected fault', 'FontWeight', 'bold');
zlabel('P(fi diag|fj)');
title('Fault Isolation Performance Matrix')
b=bar3(C); for finjIdx=1:nf
set(gca, 'XTickLabel', data.modes, 'TickLabelInterpreter', 'none') for fdiagIdx=1:nf
set(gca, 'YTickLabel', data.modes, 'TickLabelInterpreter', 'none') text(fdiagIdx,finjIdx,C(finjIdx,fdiagIdx)+0.05,...
xlabel('Diagnosed fault', 'FontWeight', 'bold'); sprintf('%.1f',C(finjIdx,fdiagIdx)*100),'HorizontalAlignment', 'center')
ylabel('Injected fault', 'FontWeight', 'bold'); end
zlabel('P(fi diag|fj)');
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')
end end
for k = 1:length(b)
zdata = b(k).ZData;
b(k).CData = zdata;
b(k).FaceColor = 'interp';
end
colormap('Summer')
view(0,90)
end end
for k = 1:length(b)
zdata = b(k).ZData;
b(k).CData = zdata;
b(k).FaceColor = 'interp';
end
colormap('Summer')
view(0,90)
function [dx] = EvaluateResiduals(data, FSM, resIdx) function [dx] = EvaluateResiduals(data, FSM, resIdx)
...@@ -52,8 +60,8 @@ function [dx] = EvaluateResiduals(data, FSM, resIdx) ...@@ -52,8 +60,8 @@ function [dx] = EvaluateResiduals(data, FSM, resIdx)
% Run basic consistency based fault isolation strategy (assumes threshold % Run basic consistency based fault isolation strategy (assumes threshold
% = 1) % = 1)
fprintf('Running basic consistency based fault isolation ... '); % fprintf('Running basic consistency based fault isolation ... ');
tic; %tic;
dx = SingleFaultIsolation(selectedRes, FSM(resIdx,:)); dx = SingleFaultIsolation(selectedRes, FSM(resIdx,:));
fprintf(' Finished in %.2f sec.\n', toc); %fprintf(' Finished in %.2f sec.\n', toc);
end
function [im,df,ndf] = IsolabilityAnalysisFSM(fsm, f, varargin)
% IsolabilityAnalysisFSM Perform structural single fault isolability analysis of a Fault Signature Matrix (FSM)
%
% [im,df,ndf] = IsolabilityAnalysisFSM(fsm, options )
%
% With no output arguments, then the command plots the isolability
% analysis results.
%
% Options are key/value pairs
%
% Inputs:
% fsm - Fault signature matriux
%
% Key Value
% permute If true, permute the fault variables such that the
% isolability matrix gets a block structure for easier
% interpretation when plotted. Does not affect the output
% argument im, only the plot (default true)
%
%
% Outputs:
% im - Isolability matrix, im(i,j)=1 if fault i can be isolated
% from fault j, 0 otherwise
% df - Detectable faults
% ndf - Non-detectable faults
%
p = inputParser;
p.addOptional('permute',true);
p.parse( varargin{:} );
opts = p.Results;
% Compute isolability matrix
nf = size(fsm,2);
nr = size(fsm,1);
im = ones( nf, nf );
for k=1:nr
im(fsm(k,:)>0,fsm(k,:)==0)=0;
end
% Compute detectable and non-detectable faults
ndf = f(any(fsm,1)==0);
df = f(any(fsm,1)>0);
% Plot
if nargout==0
if opts.permute
[p,q] = dmperm(im);
else
p = 1:numel(f);
q = p;
end
spy(im(p,q), 40)
set(gca,'XTick', 1:nf);
set(gca,'YTick', 1:nf);
if verLessThan('matlab', '8.4')
set(gca,'XTickLabel', f(p));
set(gca,'YTickLabel', f(p));
else
set(gca,'XTickLabel', f(p), 'TickLabelInterpreter','none');
set(gca,'YTickLabel', f(p), 'TickLabelInterpreter','none');
end
xlabel('')
end
end
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment