ACS6124 Laboratory I: Multisensor Systems Visakan Kadirkamanathan and Andrew Hills 2019 Aims 1. To relate theory and practice of multisensor and decision systems in application. 2. To demonstrate a comprehensive overview how different monitoring approaches are im- plemented and designed in practice. Objectives At the end, the student will be able to: 1. Appreciate the design procedure of sensor signal feature extraction. 2. Extract the informative features for health monitoring application. 3. Implement the classification algorithm based decision-making system. 4. Evaluate decision system performance for design choices. How to use MATLAB files 1. Copy all files for each lab into your working directory. 2. To get help, type help following by the name of the m-file or function. For example, if you want to know how to use a function mean, type help mean. 1 Lab A.I — Data Acquisition and Frequency Analysis Aims 1. To relate theory and practice of sensor signal analysis for health monitoring applications. 2. To demonstrate how different monitoring approaches are implemented and designed in practice. Objectives At the end, the student will be able to 1. Appreciate the design procedure of sensor signal feature extraction for health monitoring decision system. 2. Analyse sensor signal feature extraction design choices. 3. Use MATLAB commands to implement sensor signal analysis functions for health moni- toring. How to use MATLAB files Copy all files in the directory \\uosfstore.shef.ac.uk\Shared\studentacse\public\ vk\ACS6124\FeaturesLab into your working directory in U drive and add the corresponding working path to MATLAB. 2 Sensor Signal based Fault Diagnosis Features are parameters derived from the measured sensor data that robustly characterise the system it is observing, and in this laboratory, indicate the presence of faults. Objective of this lab is to obtain a feature extractor from a set of vibration sensor signals for the purposes of fault diagnosis of rotating machinery. The vibration signals of a machine normally carry the dynamic information of the machine and analysing them will help us in monitoring the condition of the machine. Here your tasks are to analyse a set of vibration sensor signals representing five different fault types and extract features from the sensor signals. The extracted features will be used in diagnosing of machine conditions using a pattern classification technique. • Test Rigs The vibration sensor data were collected from five test rigs that are used to demonstrate effects of the following built-in faults: Fault 1: Bearing Defect Fault 2: Gear mesh Fault 3: Resonance Fault 4: Imbalance Fault 5: Misalignment A brief description of each fault type is given in the Appendix. Each test rig is composed of: a rotating machinery system, a 1-axis accelerometer, a data acquisi- tion card (PCI 1711) and a PC, as shown in Figure 1. The vibration sensor data measured by the Figure 1: Machinery Vibration Monitoring System 1-axis accelerometer is captured by a PC-based data acquisition system via PCI 1711. The Real- Time Windows Target (RTWT) of MATLAB is used to facilitate the communication between PC and the A/D card. 3 Task I: Vibration Sensor Signal Analysis The purpose of this task is to perform frequency analysis and feature extraction of the vibration sensor measurement signals. • Files Used bearing.mat: contains a variable bearing which is a 50000×1 column vector of vibration measurements of the bearing-defect rig sampled over 50 sec at a frequency of 1 kHz. gearmesh.mat: contains a variable gearmesh which is a 50000 × 1 column vector of vibration measurements of the gearmesh rig sampled over 50 sec at a frequency of 1 kHz. misalignment.mat: contains a variable misalignment which is a 50000×1 column vec- tor of vibration measurements of the misalignment rig sampled over 50 sec at a frequency of 1 kHz. imbalance.mat: contains a variable imbalance which is a 50000 × 1 column vector of vibration measurements of the imbalance rig sampled over 50 sec at a frequency of 1 kHz. resonance.mat: contains a variable resonance which is a 50000 × 1 column vector of vibration measurements of the resonance rig sampled over 50 sec at a frequency of 1 kHz. To analyse the signals, follow the subsequent steps: 1. Load data files of all fault types and plot them in separate figures. Consider the vibra- tion signature in the time-domain of each fault and decide on features, which you believe distinguish different fault types. 2. Analyse signals in the frequency domain by spectral analysis. The energy distribution of the signal can be made by estimating its power spectral density (PSD). Here the Welchs method will be applied using pwelch command. Then, plot the estimated PSDs in sepa- rate figures. % Estimate the PSD of y >>[P,f]=pwelch(y,[],[],[],1000);%1000 is the sampling frequency (Hz) >>plot(f,p); In the frequency domain, decide on features, which you believe distinguish different fault types. Task II: Signal Feature Extraction This task allows you to extract appropriate features from the vibration data indicative of the fault condition. From the energy distributions, it can be seen that a large variety of features can be extracted to describe the characteristics of signals. Here we consider the energy levels which are given by the value of root mean square (RMS) of the signal. The energy levels in three frequency bands: 0–50 Hz, 50–200 Hz and 200–500 Hz will be included in the feature set. In summary, four selected features are given in Table 1. A 50000×1 sensor output vector of each fault case will be divided into (non-overlapping) 50 blocks each of length 1000. Features mentioned above are calculated based on each correspond- ing a one-second period of the measurements and providing 4 features. Features f1 to f4 of each block can be calculated by the following steps: 4 Feature Name Detail of the Feature f1 Root Mean Square (RMS) of the signal f2 RMS of low frequency: 0–50 Hz f3 RMS of middle frequency: 50–200 Hz f4 RMS of high frequency: 200–500 Hz Table 1: Features to be extracted • Normalisation: Each 1000× 1 vector is normalised for zero mean by xnormalised = x− xmean (1) where x is any 1000 × 1 data block. Let’s call the first 1000 normalised data x1, the next 1000 normalised data x2 and so on. Each xi vector is then used to extract the features f1 − f4 as follows: 1. Feature f1: calculate the PSD of xi by Welch’s method. The RMS (feature f1) is calculated by f1 = ‖xi‖2√ N (2) where N is the length of the vector xi which is 1000 in this case. ‖xi‖2 is the 2-norm or the Euclidean length of xi. The 2-norm can be obtained by using a MATLAB command norm. 2. Feature f2: signal is filtered through a low pass filter with a cut-off frequency of 50 Hz and the RMS-value is then calculated. Type of the filter here is Butterworth, which provides a maximally flat response in pass-band as well as in stop-band. Create an 11th order low pass Butterworth digital filter by butter command 1 which re- turns the filter coefficients vectors B (numerator) and A (denominator), i.e. [B,A]=butter(11,0.1); • Apply the filter to xi using filter MATLAB function: y1= filter(B,A,x1); % for x1 data vector; f2 is then obtained by calculating the RMS of the PSD of the filtered xi using (1). 3. Feature f3: Signal is filtered through a band pass filter with a cut-off frequency of 50–200 Hz and the RMS-value is then calculated. • Create a 13th order band pass Butterworth digital filer by butter command: [B,A] = butter(13,[0.1 0.4]); 1The second argument 0.1 (ωn) comes from the scaled cutoff frequency 50 Hz with 1.0 corresponding to half the sample rate, i.e. ωn = 50500 . 5 • Apply the filter to xi using filter command: y1= filter(B,A,x1); % for x1 data vector • f3 is then obtained by calculating the RMS of the PSD of the filtered xi using (1) 4. Feature f4: Signal is filtered through a high pass filter with a cut-off frequency of 200 Hz and the RMS-value is then calculated. • Create an 18th order high pass Butterworth digital filter by butter command: [B,A] = butter(18,0.
4,’high’); • Apply the filter to xi using filter command2: y1= filter(B,A,x1); % for x1 data vector Task III: Data Visualisation From Task II, you should have a set of five feature matrices each of which is composed of four features of signal energy levels in different frequency bands. As the feature dimension is four, it is impossible to visualise them and a dimension reduction is required. To this purpose, the four-dimensional feature vectors will be mapped to two dimensions by principal component analysis (PCA). First, use load to load in the feature files of fault types 1 to 5 created in Lab. Then,combine the feature matrices of fault cases 1–5 into a new matrix G: G = Fault1 Fault2 … Fault5 (3) The first 2 principal components of G are then calculated and used in the visualisation by the following MATLAB commands: c=corrcoef(G); % Calculates a correlation coefficient matrix c of G [v,d] =eig(c); % Find the eigenvectors v and the eigenvalues d of G T=[ v(:,end)’;v(:,end-1)’]; % Create the transformation matrix T from % the first two principal components z=T*G’; % Create a 2-dimensional feature vector z plot(z(1,:),z(2,:),’o’) % Scatter plot of the 2-dimensional features Is each fault type completely separate from the others? 2All filter orders are chosen by using buttord function which returns the lowest order digital Butterworth filter that loses no more than 3 dB in the passband and has at least 20 dB of attenuation in the stopband. 6 Feature Name Filter Filter Order Wn f1 Low-pass 25 Hz 7 0.05 f2 Band-pass 25 – 50 Hz 6 [0.05 0.1] f3 Band-pass 50 – 100 Hz 9 [0.1 0.2] f4 Band-pass 100 – 200 Hz 8 [0.2 0.4] f5 Band-pass 200 – 350 Hz 9 [0.4 0.7] f6 High-pass 350 Hz 16 0.7 Table 2: Features to be extracted Task IV: Optional Additional Investigation Now we consider the energy levels in six frequency bands: 0–25 Hz, 25–50 Hz, 50–100 Hz, 100– 200 Hz, 200–350 Hz and 350–500 Hz will be included in the feature set. Filter orders and Wn for each frequency band can be found in Table 2. Use these values and repeat Task II and Task III. Create and plot the 2-dimensional feature and compare with the result obtained in Task III. 7 Lab A.II — Pattern Classification Aims 1. To relate theory and practice of decision systems for health monitoring application. 2. To demonstrate how different monitoring approaches are implemented and designed in practice. Objectives At the end, the student will be able to 1. Understand the design procedure decision systems based on pattern classification ap- proaches. 2. Evaluate decision system performance for design choices on health monitoring applica- tion. How to use MATLAB files Copy all files in the directory \\uosfstore.shef.ac.uk\Shared\studentacse\public\ vk\ACS6124\ClassificationLab into your working directory in U drive and add the cor- responding working path to MATLAB. 8 Pattern Classification The task here is to classify machine conditions based on a simple classifier, namely the 1-nearest neighbour classifier. The 4-dimensional feature vectors are first transformed into 2-dimensional data to observe the data groups (this is already done in Lab A.I). The feature data are then partitioned for training and testing. The nearest neighbour classifier is implemented and its performance is estimated. Training data and test data matrices: The matrices used as the training data and the testing data will be created from the extracted features in Lab A.I. The training data are used in construction of the classifier and the test data is the data set used to evaluate the goodness or performance of the classifier. Here, the first 35 feature vectors of each fault case are used as the training data and the last 15 feature vectors are then used as the test data. Set two MATLAB variables as the training data and test data as follows: >> trainingSet = [Fault1(1:35,:); Fault2(1:35,:); Fault3(1:35,:); Fault4(1:35,:); Fault5(1:35,:)]; >> testingSet = [Fault1(36:end,:); Fault2(36:end,:);Fault3(36:end,:); Fault4(36:end,:); Fault5(36:end,:)]; By using the training data as a reference set of samples with known fault types (fault classes 1–5), an unknown sample of the test data is assigned to class of its nearest sample of the training data (see the course handout). You may also find it useful to generate training target and test target vectors corresponding to the associated fault class labels (1 for Fault 1, 2 for Fault 2 and so on) of the training data and the test data, respectively. These two target vectors can be used in the nearest neighbour classifier to assign a fault class to a sample and evaluate the classifiers performance. In this case, the training target (traintar) and the test target vectors (testtar) can be created by the following MATLAB commands >> trainingTarget=[ones(1,35) ones(1,35)*2 ones(1,35)*3 ones(1,35)*4 ones(1,35)*5]; >> testingTarget=[ones(1,15) ones(1,15)*2 ones(1,15)*3 ones(1,15)*4 ones(1,15)*5]; Write a program to classify fault conditions based on the 1-nearest neighbour classifier (1-NN) by using train as a training data set and test as a test data set and evaluate accuracy of classifi- cation (ACC) defined by %ACC = 100× Nc Na (4) where Nc is the number of all correctly classified samples and Na is the number of samples. Hint: Nc may be calculated by using find function e.g. Nc = length(find(group==testtar)) where group is a vector of the class labels (1–5) assigned to the test data , e.g. group=[1 1 2 2 2 2 5…] and testtar is the target vector of the test data. 9 Bonus Marks Extend the 1-nearest neighbours algorithm developed above to create a k-nearest neighbours solution (where k = 1, 3 and 5) and compare the classification accuracy. Template for Nearest Neighbours algorithm To assist you in this exercise, the template code is provided to create the nearest neighbours algorithm. Please note that the implementation of the algorithm is not unique, so your own solution may differ to that provided in the template code. Distance Measure The Euclidean distanc, d, between points a and b (in n-dimensional space) is defined as: d = √ (a1 − b1)2 + (a2 − b2)2 + · · ·+ (an − bn)2 Which can be simplified to just: d = √√√√ n∑ i=1 (ai − bi)2 Implemented in MATLAB as the function euc.m, this is: function distance = euc(a, b) %Euclidean Distance % Calculates the Euclidean distance between two cases which an equal % number of features. % % Author: Andrew Hills % Date: 23/02/07 if nargin ˜= 2 error(’Two input arguments required.’); return; end if ˜all(size(a) == size(b)) error(’Dimensions of inputs are not equal.’); return; end if min(size(a)) ˜= 1 error(’Input is not a vector’); return; end % Calculate the Euclidean Distance using the MATLAB’s norm function distance = norm(a – b); 10 1-Nearest Neighbours Algorithm Firstly, the variables need to be initialised: % Specify the number of training cases: numberOfTrainingCases = 35; trainingSet = [Fault1(1:numberOfTrainingCases,:); Fault2(1:numberOfTrainingCases,:); Fault3(1:numberOfTrainingCases,:); Fault4(1:numberOfTrainingCases,:); Fault5(1:numberOfTrainingCases,:)]; testingSet = [Fault1(numberOfTrainingCases+1:end,:); Fault2(numberOfTrainingCases+1:end,:); Fault3(numberOfTrainingCases+1:end,:); Fault4(numberOfTrainingCases+1:end,:); Fault5(numberOfTrainingCases+1:end,:)]; Then both the training set and testing set need to be labelled. This can be easily accom- plished by labelling Fault1 as 1, Fault2 as 2, etc. This is done for both the training and testing sets. % Note the below works because all faults are of equal lengths. numberOfTestingCases = length(Fault1) – numberOfTrainingCases; trainingTarget = [ones(1,numberOfTrainingCases), ones(1,numberOfTrainingCases)*2, ones(1,numberOfTrainingCases)*3, ones(1,numberOfTrainingCases)*4, ones(1,numberOfTrainingCases)*5]; testingTarget = [ones(1,numberOfTestingCases), ones(1,numberOfTestingCases)*2, ones(1,numberOfTestingCases)*3, ones(1,numberOfTestingCases)*4, ones(1,numberOfTestingCases)*5]; 11 Using the above initial variables, the algorithm that performs the 1-nearest neighbour search will look
as follows: % Calculate the total number of test and train classes totalNumberOfTestingCases = numberOfTestingCases * 5; totalNumberOfTrainingCases = numberOfTrainingCases * 5; % Create a vector to store assigned labels inferredLabels = zeros(1, totalNumberOfTestingCases); % This loop cycles through each unlabelled item: for unlabelledCaseIdx = 1:totalNumberOfTestingCases unlabelledCase = testingSet(unlabelledCaseIdx, :); % As any distance is shorter than infinity shortestDistance = inf; shortestDistanceLabel = 0; % Assign a temporary label % This loop cycles through each labelled item: for labelledCaseIdx = 1:totalNumberOfTrainingCases labelledCase = trainingSet(labelledCaseIdx, :); % Calculate the Euclidean distance: currentDist = euc(unlabelledCase, labelledCase); % Check the distance if currentDist < shortestDistance shortestDistance = currentDist; shortestDistanceLabel = trainingTarget(labelledCaseIdx); end end % inner loop % Assign the found label to the vector of inferred labels: inferredLabels(unlabelledCaseIdx) = shortestDistanceLabel; end % outer loop 12 Appendix — Fault Types Figure 2: Machine components and typical fault types (reproduced from [3] ) • Bearing defect: Bearing defect is a change in the condition of a bearing from the new state of installation either as a result of wear of the rollers or races- or as a result of sub surface cracking, pitting, corrosion or overheating, or cage/roller fracture of indenting due to over loading. Damage is often associated with lack of lubricating. • Gear damage: Typical problems of gear boxes are: eccentric gears, gear mesh wear, backlash, broken or chipped teeth and gear box distortion [2]. A key indicator of gear tooth wear is excitation of the Gear Natural Frequency (GNF), which is a potential vibration frequency on any machine that contains gears; equal to the number of teeth multiplied by the rotational frequency of the shaft, along with sidebands around it spaced at the running speed of the bad gear. • Resonance: Resonance occurs when a forcing frequency coincides with a system natural frequency, and can cause dramatic amplitude amplification which can result in premature or even catastrophic failure. This may be a natural frequency of the rotor but can often originate from a support frame, foundation, gearbox or drive belts. • Imbalance: Imbalance in a rotor denotes that the centre of gravity and the geometric centre of a disk are not at the same location. These two points can never be same even for a perfectly made disk, since no material is homogeneous [2]. That is why the imbalance of the rotor becomes one of the most vexing problems for rotating machine maintenance. • Misalignment: Alignment is a condition whereby machine components have the correct angular position relative to each other; either coincident, parallel, or perpendicular, according to design requirements. Poor alignment is one of the major causes of premature machine failure. 13 Figure 3: Angular and parallel misalignments (reproduced from [3]) References [1] Leeds Polytechnic, An Introduction to Plant Condition Monitoring and Its Management, 1991. [2] J.S. Rao, Vibratory Condition Monitoring of Machines, Narosa Publishing House, 2000. [3] Machine Vibration Diagnosis, http://www.vibanalysis.co.uk/. 14
