%% Anomaly detection in an Industrial Machine using 3-axis vibration measurements
%
% This example shows how to detect anomalies in vibration data using different techniques such as machine learning and deep learning. We work with vibration data collected from an indutrial machine and start by extracting features from the raw measurements using the Diagnostic Feature Designer App. The selected features are used to train each model under consideration. Each trained model is then used to identify whether the machine is operating under normal or anomalous condition.
% Dataset
% The dataset contains 3-axis vibration measurements from an industrial machine under both anomalous condition with a faulty component and normal condition with component replaced. Data for each axis is stored in a separate column and each file contains 7000 measurements. Load the data from the MathWorks supportfiles site (this is a large dataset, ~70 MB) and plot a sample of the normal and anomalous dataset.

load("MachineData.mat")
trainData

% To better understand the data, visualize it in normal and anomalous conditions. Pick the fourth member of the ensemble for the plot and you can see that the data for the two conditions look noticeably different.

ensMember = 4;
plotVibrationData(trainData, ensMember)

% Extracting features with Diagnostic Feature Designer App
% Using raw data for training machine learning models is not very efficient. The Diagnostic Feature Designer App lets you interactively explore and preprocess your data, extract time- and frequency-domain features, then rank the features to determine which will be most effective. You can then export a function to extract the selected features from your dataset programmatically. Open the DFD app:

diagnosticFeatureDesigner

% Press the New Session button, select the trainData as the source and then set label as a condition variable.
%
% You can use the Diagnostic Feature Designer to iterate on the features and rank them. A histogram view for all generated features is created to visualize the distributions separated by labels for various features extracted from Ch1. Note that the histograms shown below are from a much larger dataset, so the separation is easier to visualize. The dataset used so far is a subset of the original data, so your results might be different.
%
%
%
% Use the top four ranked features for each channel.
% Ch1 : Crest Factor, Kurtosis, RMS, Std
% Ch2 : Mean, RMS, Skewness, Std
% Ch3 : Crest Factor, SINAD, SNR, THD
% Export a function to generate the features from the Diagnostic Feature designer app and save it with the name generateFeatures. This function extracts the top 4 relevant features from each channel in the entire dataset from the commandline.

trainFeatures = generateFeatures(trainData);

%% Train Models for Anomaly Detection
% The dataset we looked at so far is only a small subset of a much larger dataset to illustrate the process of feature extraction and selection. Training your algorithm on all available data yields the best performance. To this end, load the 12 features that were previously extracted from the larger dataset of 17,642 signals.

load("FeatureEntire.mat")
head(featureAll)

%
% Use cvpartition to partition data into a training set and an independent test set.

rng(0) % set for reproducibility
idx = cvpartition(featureAll.label, 'holdout', 0.1);
feature_train = featureAll(idx.training, :);
feature_val = featureAll(idx.test, :);

% trueAnomaliesTest = renamecats(feature_val.label,["Before","After"],["Anomaly","Healthy"]);
trueAnomaliesTest = feature_val.label;
tmp_Train = feature_train(feature_train.label=='Normal', :);

%
%% Detect Anomalies with One-class SVM
% Support Vector Machines are powerful classifiers, and their variant, which trains on only one-class, models the "normal" data. This model works well for identifying abnormalities as being "far" from the normal data. Select the training data for normal condition and train the SVM model using fitcsvm function.

mdlSVM = fitcsvm(tmp_Train, 'label', 'Standardize', true, 'OutlierFraction', 0);

% Validate the trained SVM model by using test data which contains both normal and anomalous data.

tmp_all = feature_val(:, 2:end);
[~,scoreSVM] = predict(mdlSVM,tmp_all);
isanomalySVM = scoreSVM<0;
predSVM = categorical(isanomalySVM, [1, 0], ["Anomaly", "Normal"]);

confusionchart(trueAnomaliesTest, predSVM, Title="Anomaly Detection with One-class SVM", Normalization="row-normalized");

% From the confusion matrix, it is seen that the one-class SVM does pretty well. Only 0.3% of anomalous samples are misclassified as normal and about 0.9% of normal data is misclassified as anomalous.
%
%% Detect Anomalies with Isolation Forest
% The decision trees of an isolation forest isolate each observation in a leaf. How many decisions a sample passes through to get to its leaf is a measure of how complicated it was to isolate it from the others. The average depth of trees for a specific sample is used as their anomaly score, and returned by the iforest. Once again, train the isolation forest model on normal data only.

% [mdlIF,~,scoreTrainIF] = iforest(trainHealthy,'ContaminationFraction',0.05);
[mdlIF,~,scoreTrainIF] = iforest(tmp_Train(:,2:13),'ContaminationFraction',0.01);

% Validate the trained isolation forest model by using test data which contains both normal and anomalous data. Visualize the performance of this model by using a Confusion Chart.

[isanomalyIF,scoreTestIF] = isanomaly(mdlIF,tmp_all);
predIF = categorical(isanomalyIF, [1, 0], ["Anomaly", "Normal"]);

confusionchart(trueAnomaliesTest,predIF,Title="Anomaly Detection with Isolation Forest",Normalization="row-normalized");

% On this data, the isolation forest doesn't do as well as the one-class SVM, but it's not your typical anomaly detection with so many "anomalies" in this almost balanced data set. The reason for this is that the training data contained normal data only while the test data contains about 30% anomalous data. This indicates that the isolation forest model is better suited when the distribution of the anomalies and normal data is closer.
%
%% Detect anomalies with LSTM Autoencoder Network
% Autoencoders are a type of neural network that learn a compressed representation of unlabeled data. LSTM autoencoders are a variant of this network, which can learn a compressed representation of sequence data. This LSTM autoencoder will be trained with normal data only and the trained network will be used to identify when the sequence does not look like normal data is encountered.

X = extractLabeledData(feature_train, ...
   "Label", "Normal");

% Then, construct the biLSTM autoencoder network and set the training options.

featureDimension = 1;

% Define biLSTM network layers
layers = [ sequenceInputLayer(featureDimension, 'Name', 'in')
   bilstmLayer(16, 'Name', 'bilstm1')
   reluLayer('Name', 'relu1')
   bilstmLayer(32, 'Name', 'bilstm2')
   reluLayer('Name', 'relu2')
   bilstmLayer(16, 'Name', 'bilstm3')
   reluLayer('Name', 'relu3')
   fullyConnectedLayer(featureDimension, 'Name', 'fc')
   regressionLayer('Name', 'out') ];

% Set Training Options
options = trainingOptions('adam', ...
   'Plots', 'training-progress', ...
   'MiniBatchSize', 500,...
   'MaxEpochs',50);

% Now we're ready to train the model. .

net = trainNetwork(X, X, layers, options);

%
%% Visualize Model Behavior and Error on Validation Data
% Extract and visualize a sample each from Anomalous and Normal condition. Note that the plots below are comparing the error in the values of each of the 12 features (indicated on the x-axis). In this sample, features 10, 11 and 12 do not reconstruct well for the anomalous input and thus have a high error value. This is used to identify an anomaly.

testNormal = {feature_val(1200, 2:end).Variables};
testAnomaly = {feature_val(200, 2:end).Variables};

% Predict decoded signal for both
decodedNormal = predict(net,testNormal);
decodedAnomaly = predict(net,testAnomaly);

% Visualize
visualizeModelBehavior(testNormal, testAnomaly, decodedNormal, decodedAnomaly)

% Extract features for normal and anomalous data. We can see in the figure that the reconstruction error for the anomalous data is significantly higher than the normal data. This makes sense, since the autoencoder was trained on the normal data, so it will better reconstruct similar signals.

% Extract data Before maintenance
XTestAnomaly = extractLabeledData(feature_val, ...
   "Label", "Anomaly");

% Predict output before maintenance and calculate error
yHatAnomaly = predict(net, XTestAnomaly);
errorAnomaly = calculateError(XTestAnomaly, yHatAnomaly);

% Extract data after maintenance
XTestNormal = extractLabeledData(feature_val, ...
   "Label", "Normal");

% Predict output after maintenance and calculate error
yHatNormal = predict(net, XTestNormal);
errorNormal = calculateError(XTestNormal, yHatNormal);

visualizeError(errorAnomaly, errorNormal);

%
%% Identify Anomalies
% Calculate error on full validation data.

XTestAll = extractLabeledData(feature_val, ...
   "Label", "All");

yHatAll = predict(net, XTestAll);
errorAll = calculateError(XTestAll, yHatAll);

% Define an anomaly as a point which has reconstruction error 0.5 times larger than the mean across all observations. This threshold was determined through experimentation and can be changed as required.

thresh = 0.5;
anomalies = errorAll > thresh*mean(errorAll);

visualizeAnomalies(anomalies, errorAll, feature_val);

%
%% Summary
% Anomaly detection was performed on timeseries signals from a three-axis accelerometer on an industrial machine. The top four features were selected and extracted for data in each channel. These features were used to detect anomalies using one-class SVM model, isolationForest model and a LSTM autoencoder network. For each case, the trained model's performance is evaluated using a confusion matrix.
%
%% Helper functions
%
% Calculate Error

function E = calculateError(X, Y)

E = zeros(length(X),1);
for i = 1:length(X)
   E(i,:) = sqrt(sum((Y{i} - X{i}).^2));
end

end

% Visualize Error

function visualizeError(errorBefore, errorAfter)

figure("Color", "W")
tiledlayout("flow")

nexttile
plot(1:length(errorBefore), errorBefore, 'LineWidth',1.5), grid on
title(["Before Maintenance", ...
   sprintf("Mean Error: %.2f\n", mean(errorBefore))])
xlabel("Observations")
ylabel("Reconstruction Error")
ylim([0 15])

nexttile
plot(1:length(errorAfter), errorAfter, 'LineWidth',1.5), grid on,
title(["After Maintenance", ...
   sprintf("Mean Error: %.2f\n", mean(errorAfter))])
xlabel("Observations")
ylabel("Reconstruction Error")
ylim([0 15])

end

% Visualize Anomalies

function visualizeAnomalies(anomalies, errorAll, featureTest)

anomalyIdx = find(anomalies);
anomalyErr = errorAll(anomalies);

predAE = categorical(anomalies, [1, 0], ["Anomaly", "Normal"]);
% trueAE = renamecats(featureTest.label,["Before","After"],["Anomaly","Healthy"]);
trueAE = featureTest.label;
acc = numel(find(trueAE == predAE))/numel(predAE)*100;
figure;
t = tiledlayout("flow");
title(t, "Test Accuracy: " + round(mean(acc),2) + "%");
nexttile
hold on
plot(errorAll)
plot(anomalyIdx, anomalyErr, 'x')
hold off
ylabel("Reconstruction Error")
xlabel("Observation")
legend("Error", "Candidate Anomaly")

nexttile
confusionchart(trueAE,predAE)

end

% Visualize Model Behavior

function visualizeModelBehavior(normalData, abnormalData, decodedNorm, decodedAbNorm)
%visualizeModelBehavior Visualize model behavior on sample validation data

arguments
   normalData(1,1) cell
   abnormalData(1,1) cell
   decodedNorm(1,1) cell
   decodedAbNorm(1,1) cell
end

figure("Color", "W")
tiledlayout("flow")

nexttile()
hold on
colororder('default')
yyaxis left
plot(normalData{:})
plot(decodedNorm{:},":","LineWidth",1.5)
hold off
title("Normal Input")
grid on
ylabel("Feature Value")
yyaxis right
stem(abs(normalData{:} - decodedNorm{:}))
ylim([0 2])
ylabel("Error")
legend(["Input", "Decoded","Error"],"Location","southwest")

nexttile()
hold on
yyaxis left
plot(abnormalData{:})
plot(decodedAbNorm{:},":","LineWidth",1.5)
hold off
title("Abnormal Input")
grid on
ylabel("Feature Value")
yyaxis right
stem(abs(abnormalData{:} - decodedAbNorm{:}))
ylim([0 2])
ylabel("Error")
legend(["Input", "Decoded","Error"],"Location","southwest")

end

function X = extractLabeledData(featureTable, options)
%EXTRACTLABELEDDATA Extract data from normal or abnormal operating conditions
% and re-format to support input to autoencoder network

arguments
   featureTable table
   options.Label(1,1) string {mustBeMember(options.Label,["Normal", "Anomaly", "All"])} = "All"
end

% Select data with label After
if options.Label == "All"
   Xtemp = featureTable(:, 2:end).Variables;
else
   tF = featureTable.label == options.Label;
   Xtemp = featureTable(tF, 2:end).Variables;
end

% Arrange data into cells
X = cell(length(Xtemp),1);
for i = 1:length(Xtemp)
   X{i,:} = Xtemp(i,:);
end

end

% Copyright 2021 The MathWorks, Inc.