...
 
Commits (6)
# Created by https://www.gitignore.io/api/emacs,matlab,jupyternotebook
### Emacs ###
# -*- mode: gitignore; -*-
*~
\#*\#
/.emacs.desktop
/.emacs.desktop.lock
*.elc
auto-save-list
tramp
.\#*
# Org-mode
.org-id-locations
*_archive
# flymake-mode
*_flymake.*
# eshell files
/eshell/history
/eshell/lastdir
# elpa packages
/elpa/
# reftex files
*.rel
# AUCTeX auto folder
/auto/
# cask packages
.cask/
dist/
# Flycheck
flycheck_*.el
# server auth directory
/server/
# projectiles files
.projectile
projectile-bookmarks.eld
# directory configuration
.dir-locals.el
# saveplace
places
# url cache
url/cache/
# cedet
ede-projects.el
# smex
smex-items
# company-statistics
company-statistics-cache.el
# anaconda-mode
anaconda-mode/
### JupyterNotebook ###
.ipynb_checkpoints
*/.ipynb_checkpoints/*
__pycache__
# Remove previous ipynb_checkpoints
# git rm -r .ipynb_checkpoints/
#
### Matlab ###
##---------------------------------------------------
## Remove autosaves generated by the Matlab editor
## We have git for backups!
##---------------------------------------------------
# Windows default autosave extension
*.asv
# OSX / *nix default autosave extension
*.m~
# Compiled MEX binaries (all platforms)
*.mex*
# Simulink Code Generation
slprj/
# Session info
octave-workspace
# Simulink autosave extension
*.autosave
# End of https://www.gitignore.io/api/emacs,matlab,jupyternotebook
\ No newline at end of file
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import scipy.io as sio
def IsolabilityMatrix(fsm):
"""Compute isolability matrix based on a fault signature matrix"""
nf = fsm.shape[1]
im = np.ones((nf, nf), dtype=np.int)
for ri in fsm:
im[np.ix_(ri>0, ri==0)] = 0
return im
def DiagnosesAndConfusionMatrix(data, residx=None):
"""Compute consistency based diagnoses and corresponding confusion matrix based on a dataset."""
if isinstance(residx, type(None)):
dx = SingleFaultIsolability(data['res'], data['fsm'])
else:
dx = SingleFaultIsolability(data['res'][:, residx], data['fsm'][residx, :])
nf = data['fsm'].shape[1]
C = np.zeros((nf, nf))
for fi in range(nf):
for fj in range(nf):
fjIdx = data['mode']==fj
C[fj, fi] = np.sum(dx[fjIdx, fi])/np.sum(fjIdx)
return dx, C
def SingleFaultIsolability(res, fsm):
"""Compute sinlg fault consistency based diagnoses."""
N = res.shape[0]
nf = fsm.shape[1]
dx = np.zeros((N, nf))
for k, rk in enumerate(res):
alarm = (rk>0.5)
if np.any(alarm):
dx[k] = np.all(fsm[alarm], axis=0)
else:
dx[k] = np.hstack((True, np.zeros(nf-1)))
return dx
def make_colormap(seq):
"""Return a LinearSegmentedColormap
seq: a sequence of floats and RGB-tuples. The floats should be increasing
and in the interval (0,1).
"""
seq = [(None,) * 3, 0.0] + list(seq) + [1.0, (None,) * 3]
cdict = {'red': [], 'green': [], 'blue': []}
for i, item in enumerate(seq):
if isinstance(item, float):
r1, g1, b1 = seq[i - 1]
r2, g2, b2 = seq[i + 1]
cdict['red'].append([item, r1, r2])
cdict['green'].append([item, g1, g2])
cdict['blue'].append([item, b1, b2])
return mcolors.LinearSegmentedColormap('CustomMap', cdict)
summer_cmap = make_colormap([
(0.000, 0.500, 0.400),
(0.016, 0.508, 0.400),
(0.032, 0.516, 0.400),
(0.048, 0.524, 0.400),
(0.063, 0.532, 0.400),
(0.079, 0.540, 0.400),
(0.095, 0.548, 0.400),
(0.111, 0.556, 0.400),
(0.127, 0.563, 0.400),
(0.143, 0.571, 0.400),
(0.159, 0.579, 0.400),
(0.175, 0.587, 0.400),
(0.190, 0.595, 0.400),
(0.206, 0.603, 0.400),
(0.222, 0.611, 0.400),
(0.238, 0.619, 0.400),
(0.254, 0.627, 0.400),
(0.270, 0.635, 0.400),
(0.286, 0.643, 0.400),
(0.302, 0.651, 0.400),
(0.317, 0.659, 0.400),
(0.333, 0.667, 0.400),
(0.349, 0.675, 0.400),
(0.365, 0.683, 0.400),
(0.381, 0.690, 0.400),
(0.397, 0.698, 0.400),
(0.413, 0.706, 0.400),
(0.429, 0.714, 0.400),
(0.444, 0.722, 0.400),
(0.460, 0.730, 0.400),
(0.476, 0.738, 0.400),
(0.492, 0.746, 0.400),
(0.508, 0.754, 0.400),
(0.524, 0.762, 0.400),
(0.540, 0.770, 0.400),
(0.556, 0.778, 0.400),
(0.571, 0.786, 0.400),
(0.587, 0.794, 0.400),
(0.603, 0.802, 0.400),
(0.619, 0.810, 0.400),
(0.635, 0.817, 0.400),
(0.651, 0.825, 0.400),
(0.667, 0.833, 0.400),
(0.683, 0.841, 0.400),
(0.698, 0.849, 0.400),
(0.714, 0.857, 0.400),
(0.730, 0.865, 0.400),
(0.746, 0.873, 0.400),
(0.762, 0.881, 0.400),
(0.778, 0.889, 0.400),
(0.794, 0.897, 0.400),
(0.810, 0.905, 0.400),
(0.825, 0.913, 0.400),
(0.841, 0.921, 0.400),
(0.857, 0.929, 0.400),
(0.873, 0.937, 0.400),
(0.889, 0.944, 0.400),
(0.905, 0.952, 0.400),
(0.921, 0.960, 0.400),
(0.937, 0.968, 0.400),
(0.952, 0.976, 0.400),
(0.968, 0.984, 0.400),
(0.984, 0.992, 0.400),
(1.000, 1.000, 0.400)])
def PlotConfusionMatrix(C):
"""Plot a confusion matrix in a suitable colormap."""
nf = C.shape[0]
plt.imshow(C, cmap=summer_cmap)
for fi in range(nf):
for fj in range(nf):
plt.text(fi, fj, '%.1f' % (C[fj, fi]*100), ha='center', va='center', color='k')
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_xticks(np.arange(C.shape[1]+1)-.5, minor=True)
ax.set_yticks(np.arange(C.shape[0]+1)-.5, minor=True)
ax.tick_params(which="minor", bottom=False, left=False)
ax.grid(which="minor", color="w", linestyle='-', linewidth=3)
def BoxOff(*argin):
"""Remove (ugly) box around plots."""
if len(argin)>0:
ax=argin[0]
else:
ax=plt.gca();
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
def loadmat(filename):
'''
this function should be called instead of direct spio.loadmat
as it cures the problem of not properly recovering python dictionaries
from mat files. It calls the function check keys to cure all entries
which are still mat-objects
'''
data = sio.loadmat(filename, struct_as_record=False, squeeze_me=True)
return _check_keys(data)
def _check_keys(dict):
'''
checks if entries in dictionary are mat-objects. If yes
todict is called to change them to nested dictionaries
'''
for key in dict:
if isinstance(dict[key], sio.matlab.mio5_params.mat_struct):
dict[key] = _todict(dict[key])
return dict
def _todict(matobj):
'''
A recursive function which constructs from matobjects nested dictionaries
'''
dict = {}
for strg in matobj._fieldnames:
elem = matobj.__dict__[strg]
if isinstance(elem, sio.matlab.mio5_params.mat_struct):
dict[strg] = _todict(elem)
else:
dict[strg] = elem
return dict
\ No newline at end of file
This source diff could not be displayed because it is too large. You can view the blob instead.
% 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
% 4. Licenstext?
% 5. Uppdatera readme
% 6. Kommentarer om version av matlab och figurer. Klura ut vilken version
% som använts i artikeln.
......@@ -15,7 +12,6 @@
% Consistency Based Diagnosis Using Machine Learning Models. In Safeprocess
% 2018, Warsaw, Polen.
%
% Some parts of the code requires a path to the Fault Diagnosis Toolbox.
clc
clear
......@@ -30,14 +26,15 @@ addpath src
% pretreatment of the data has been performed:
%
% 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.
% 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
% 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
% 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
% modes will be greater than one.
load data.mat
......@@ -66,16 +63,14 @@ clear absdata
no_residuals = size(data.res, 2);
all_residuals = 1:no_residuals; % select all residuals
plotDecisionMatrix(all_residuals, data);
clear no_residuals ans
clear no_residuals
%% Plot Fig 6: Structural isolability of all available residuals
load model.mat;
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');
ylabel('Injected fault', 'FontWeight', 'bold');
xlabel('Diagnsoed fault', 'FontWeight', 'bold');
clear ans
%% Create thresholded residuals
thdata = data;
......@@ -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
figure(17);
CBDPlots(thdata, model, all_residuals);
clear all_residuals ans
CBDPlots(thdata, all_residuals);
%% Train a random forest of 100 classification trees using the all data.
trees = 100;
num_trees = 100;
rng(1); % For reproducibility
t = templateTree('PredictorSelection', 'interaction-curvature');
fprintf('Building random forest classifier ... ');
tic;
tic;
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);
clear t trees
clear t num_trees
%% Plot Fig 8: Random forest classification performance using all tests on out-of-bag data.
Y = oobPredict(Mdl);
......@@ -128,7 +122,7 @@ for k = 1:length(b)
end
colormap('Summer')
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.
figure(19);
......@@ -140,7 +134,7 @@ clear oobErr
%% Estimate residual importance by permuting out-of-bag observations.
fprintf('Compute variable importance for the random forest (this may take a while) ... ');
tic;
tic;
importance = oobPermutedPredictorImportance(Mdl);
fprintf(' Finished in %.2f sec.\n', toc);
......@@ -165,7 +159,7 @@ 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);
IsolabilityAnalysisFSM(data.fsm(:, 2:8), data.modes(2:8), 'permute', 0);
stim2 = sum(((2.^(0:6)).*structuralIsolationMatrix), 2);
probabilityMaximumIsolation = zeros(42, no_modes-1);
......@@ -185,7 +179,7 @@ for i = 1:length(residualRanking)
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_residuals);
C{i} = CBDPlots(data, 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));
......@@ -198,7 +192,7 @@ clear structuralIsolationMatrix i j
% and aggregated isolation error as a function of selected tests.
figure(21);
plot([mean_md, fa, isolationerror], 'LineWidth', 2)
hold
hold on
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')
......@@ -208,7 +202,8 @@ legend('Missed detection probability', 'False alarm probability', ...
xlabel('Number of selected residuals','FontWeight', 'bold')
ylabel('Probability', 'FontWeight', 'bold')
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.
figure(22);
......@@ -218,7 +213,6 @@ set(l, 'Interpreter', 'none')
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:
......@@ -228,7 +222,7 @@ 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);
C{i} = CBDPlots(data, selected_residuals);
title(['No of tests: ' num2str(selected_no_of_residuals(i))]);
end
clear C i
......@@ -238,9 +232,8 @@ 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)), ...
IsolabilityAnalysisFSM(data.fsm(selected_residuals, (2:8)), data.modes(2:8), ...
'permute', 0); % full str isolability
title(['No of tests: ' num2str(selected_no_of_residuals(i))])
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
[dx] = EvaluateResiduals(data, FSM42, resIdx);
%% Evaluate selected set of residuals
[dx] = EvaluateResiduals(data, FSM42, resIdx);
%% Confusion matrix computation
% C(i,j) = P(fj diagnos|fi injected fault)
nf = numel(data.modes);
C = zeros(nf,nf);
%% Confusion matrix computation
% C(i,j) = P(fj diagnos|fi injected fault)
nf = numel(data.modes);
C = zeros(nf,nf);
for fi = 1:nf
for fj = 1:nf
fjIdx = (data.mode==fj-1);
C(fj,fi) = sum(dx(fjIdx,fi))/sum(fjIdx);
for fi = 1:nf
for fj = 1:nf
fjIdx = (data.mode==fj-1);
C(fj,fi) = sum(dx(fjIdx,fi))/sum(fjIdx);
end
end
end
%% Plot confusion matrix in a 3D-bar plot
if nargin==4
figure(i); clf;
end
%% Plot confusion matrix in a 3D-bar plot
b=bar3(C);
set(gca, 'XTickLabel', data.modes, 'TickLabelInterpreter', 'none')
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);
set(gca, 'XTickLabel', data.modes, 'TickLabelInterpreter', 'none')
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')
for finjIdx=1:nf
for fdiagIdx=1:nf
text(fdiagIdx,finjIdx,C(finjIdx,fdiagIdx)+0.05,...
sprintf('%.1f',C(finjIdx,fdiagIdx)*100),'HorizontalAlignment', 'center')
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
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)
......@@ -52,8 +60,8 @@ function [dx] = EvaluateResiduals(data, FSM, resIdx)
% Run basic consistency based fault isolation strategy (assumes threshold
% = 1)
fprintf('Running basic consistency based fault isolation ... ');
tic;
% fprintf('Running basic consistency based fault isolation ... ');
%tic;
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