using simulink to generate fault data -凯发k8网页登录
this example shows how to use a simulink® model to generate fault and healthy data. the fault and healthy data is used to develop a condition monitoring algorithm. the example uses a transmission system and models a gear tooth fault, a sensor drift fault and shaft wear fault.
transmission system model
the transmission casing model uses simscape™ driveline™ blocks to model a simple transmission system. the transmission system consists of a torque drive, drive shaft, clutch, and high and low gears connected to an output shaft.
mdl = 'pdmtransmissioncasing';
open_system(mdl)
the transmission system includes a vibration sensor that is monitoring casing vibrations. the casing model translates the shaft angular displacement to a linear displacement on the casing. the casing is modelled as a mass spring damper system and the vibration (casing acceleration) is measured from the casing.
open_system([mdl '/casing'])
fault modelling
the transmission system includes fault models for vibration sensor drift, gear tooth fault, and shaft wear. the sensor drift fault is easily modeled by introducing an offset in the sensor model. the offset is controlled by the model variable sdrift, note that sdrift=0 implies no sensor fault.
open_system([mdl '/vibration sensor with drift'])
the shaft wear fault is modeled by a variant subsystem. in this case the subsystem variants change the shaft damping but the variant subsystem could be used to completely change the shaft model implementation. the selected variant is controlled by the model variable shaftwear, note that shaftwear=0 implies no shaft fault.
open_system([mdl '/shaft'])
open_system([mdl,'/shaft/output shaft'])
the gear tooth fault is modelled by injecting a disturbance torque at a fixed position in the rotation of the drive shaft. the shaft position is measured in radians and when the shaft position is within a small window around 0 a disturbance force is applied to the shaft. the magnitude of the disturbance is controlled by the model variable toothfaultgain, note that toothfaultgain=0 implies no gear tooth fault.
simulating fault and healthy data
the transmission model is configured with variables that control the presence and severity of the three different fault types, sensor drift, gear tooth, and shaft wear. by varying the model variables, sdrift, toothfaultgain, and shaftwear, you can create vibration data for the different fault types. use an array of simulink.simulationinput objects to define a number of different simulation scenarios. for example, choose an array of values for each of the model variables and then use the ndgrid function to create a simulink.simulationinput for each combination of model variable values.
toothfaultarray = -2:2/10:0; % tooth fault gain values sensordriftarray = -1:0.5:1; % sensor drift offset values shaftweararray = [0 -1]; % variants available for drive shaft conditions % create an n-dimensional array with combinations of all values [toothfaultvalues,sensordriftvalues,shaftwearvalues] = ... ndgrid(toothfaultarray,sensordriftarray,shaftweararray); for ct = numel(toothfaultvalues):-1:1 % create a simulink.simulationinput for each combination of values siminput = simulink.simulationinput(mdl); % modify model parameters siminput = setvariable(siminput,'toothfaultgain',toothfaultvalues(ct)); siminput = setvariable(siminput,'sdrift',sensordriftvalues(ct)); siminput = setvariable(siminput,'shaftwear',shaftwearvalues(ct)); % collect the simulation input in an array gridsimulationinput(ct) = siminput; end
similarly create combinations of random values for each model variable. make sure to include the 0 value so that there are combinations where only subsets of the three fault types are represented.
rng('default'); % reset random seed for reproducibility toothfaultarray = [0 -rand(1,6)]; % tooth fault gain values sensordriftarray = [0 randn(1,6)/8]; % sensor drift offset values shaftweararray = [0 -1]; % variants available for drive shaft conditions %create an n-dimensional array with combinations of all values [toothfaultvalues,sensordriftvalues,shaftwearvalues] = ... ndgrid(toothfaultarray,sensordriftarray,shaftweararray); for ct=numel(toothfaultvalues):-1:1 % create a simulink.simulationinput for each combination of values siminput = simulink.simulationinput(mdl); % modify model parameters siminput = setvariable(siminput,'toothfaultgain',toothfaultvalues(ct)); siminput = setvariable(siminput,'sdrift',sensordriftvalues(ct)); siminput = setvariable(siminput,'shaftwear',shaftwearvalues(ct)); % collect the simulation input in an array randomsimulationinput(ct) = siminput; end
with the array of simulink.simulationinput objects defined use the generatesimulationensemble function to run the simulations. the generatesimulationensemble function configures the model to save logged data to file, use the timetable format for signal logging and store the simulink.simulationinput objects in the saved file. the generatesimulationensemble function returns a status flag indicating whether the simulation completed successfully.
the code above created 110 simulation inputs from the gridded variable values and 98 simulation inputs from the random variable values giving 208 total simulations. running these 208 simulations in parallel can take as much as two hours or more on a standard desktop and generates around 10gb of data. an option to only run the first 10 simulations is provided for convenience.
% run the simulations and create an ensemble to manage the simulation results if ~exist(fullfile(pwd,'data'),'dir') mkdir(fullfile(pwd,'data')) % create directory to store results end runall = true; if runall [ok,e] = generatesimulationensemble([gridsimulationinput, randomsimulationinput], ... fullfile(pwd,'data'),'useparallel', true); else [ok,e] = generatesimulationensemble(gridsimulationinput(1:10), fullfile(pwd,'data')); %#ok<*unrch> end
[21-jul-2022 15:36:27] checking for availability of parallel pool... starting parallel pool (parpool) using the 'local' profile ... connected to the parallel pool (number of workers: 6). [21-jul-2022 15:37:34] starting simulink on parallel workers... [21-jul-2022 15:38:00] configuring simulation cache folder on parallel workers... [21-jul-2022 15:38:01] transferring base workspace variables used in the model to parallel workers... [21-jul-2022 15:38:01] loading model on parallel workers... [21-jul-2022 15:38:20] running simulations... [21-jul-2022 15:38:55] completed 1 of 208 simulation runs [21-jul-2022 15:38:56] completed 2 of 208 simulation runs [21-jul-2022 15:38:57] completed 3 of 208 simulation runs [21-jul-2022 15:38:58] completed 4 of 208 simulation runs [21-jul-2022 15:38:59] completed 5 of 208 simulation runs [21-jul-2022 15:39:00] completed 6 of 208 simulation runs [21-jul-2022 15:39:09] completed 7 of 208 simulation runs [21-jul-2022 15:39:09] completed 8 of 208 simulation runs [21-jul-2022 15:39:10] completed 9 of 208 simulation runs [21-jul-2022 15:39:11] completed 10 of 208 simulation runs [21-jul-2022 15:39:12] completed 11 of 208 simulation runs [21-jul-2022 15:39:13] completed 12 of 208 simulation runs [21-jul-2022 15:39:21] completed 13 of 208 simulation runs [21-jul-2022 15:39:22] completed 14 of 208 simulation runs [21-jul-2022 15:39:23] completed 15 of 208 simulation runs [21-jul-2022 15:39:23] completed 16 of 208 simulation runs [21-jul-2022 15:39:24] completed 17 of 208 simulation runs [21-jul-2022 15:39:25] completed 18 of 208 simulation runs [21-jul-2022 15:39:34] completed 19 of 208 simulation runs [21-jul-2022 15:39:35] completed 20 of 208 simulation runs [21-jul-2022 15:39:35] completed 21 of 208 simulation runs [21-jul-2022 15:39:36] completed 22 of 208 simulation runs [21-jul-2022 15:39:37] completed 23 of 208 simulation runs [21-jul-2022 15:39:38] completed 24 of 208 simulation runs [21-jul-2022 15:39:47] completed 25 of 208 simulation runs [21-jul-2022 15:39:47] completed 26 of 208 simulation runs [21-jul-2022 15:39:48] completed 27 of 208 simulation runs [21-jul-2022 15:39:49] completed 28 of 208 simulation runs [21-jul-2022 15:39:50] completed 29 of 208 simulation runs [21-jul-2022 15:39:51] completed 30 of 208 simulation runs [21-jul-2022 15:40:00] completed 31 of 208 simulation runs [21-jul-2022 15:40:00] completed 32 of 208 simulation runs [21-jul-2022 15:40:01] completed 33 of 208 simulation runs [21-jul-2022 15:40:02] completed 34 of 208 simulation runs [21-jul-2022 15:40:03] completed 35 of 208 simulation runs [21-jul-2022 15:40:05] completed 36 of 208 simulation runs [21-jul-2022 15:40:13] completed 37 of 208 simulation runs [21-jul-2022 15:40:14] completed 38 of 208 simulation runs [21-jul-2022 15:40:15] completed 39 of 208 simulation runs [21-jul-2022 15:40:16] completed 40 of 208 simulation runs [21-jul-2022 15:40:17] completed 41 of 208 simulation runs [21-jul-2022 15:40:18] completed 42 of 208 simulation runs [21-jul-2022 15:40:26] completed 43 of 208 simulation runs [21-jul-2022 15:40:27] completed 44 of 208 simulation runs [21-jul-2022 15:40:27] completed 45 of 208 simulation runs [21-jul-2022 15:40:28] completed 46 of 208 simulation runs [21-jul-2022 15:40:29] completed 47 of 208 simulation runs [21-jul-2022 15:40:30] completed 48 of 208 simulation runs [21-jul-2022 15:40:38] completed 49 of 208 simulation runs [21-jul-2022 15:40:39] completed 50 of 208 simulation runs [21-jul-2022 15:40:40] completed 51 of 208 simulation runs [21-jul-2022 15:40:41] completed 52 of 208 simulation runs [21-jul-2022 15:40:42] completed 53 of 208 simulation runs [21-jul-2022 15:40:43] completed 54 of 208 simulation runs [21-jul-2022 15:40:52] completed 55 of 208 simulation runs [21-jul-2022 15:40:55] completed 56 of 208 simulation runs [21-jul-2022 15:40:56] completed 57 of 208 simulation runs [21-jul-2022 15:40:57] completed 58 of 208 simulation runs [21-jul-2022 15:40:58] completed 59 of 208 simulation runs [21-jul-2022 15:40:59] completed 60 of 208 simulation runs [21-jul-2022 15:41:09] completed 61 of 208 simulation runs [21-jul-2022 15:41:10] completed 62 of 208 simulation runs [21-jul-2022 15:41:11] completed 63 of 208 simulation runs [21-jul-2022 15:41:12] completed 64 of 208 simulation runs [21-jul-2022 15:41:13] completed 65 of 208 simulation runs [21-jul-2022 15:41:14] completed 66 of 208 simulation runs [21-jul-2022 15:41:20] completed 67 of 208 simulation runs [21-jul-2022 15:41:21] completed 68 of 208 simulation runs [21-jul-2022 15:41:22] completed 69 of 208 simulation runs [21-jul-2022 15:41:22] completed 70 of 208 simulation runs [21-jul-2022 15:41:23] completed 71 of 208 simulation runs [21-jul-2022 15:41:24] completed 72 of 208 simulation runs [21-jul-2022 15:41:33] completed 73 of 208 simulation runs [21-jul-2022 15:41:33] completed 74 of 208 simulation runs [21-jul-2022 15:41:34] completed 75 of 208 simulation runs [21-jul-2022 15:41:35] completed 76 of 208 simulation runs [21-jul-2022 15:41:36] completed 77 of 208 simulation runs [21-jul-2022 15:41:37] completed 78 of 208 simulation runs [21-jul-2022 15:41:45] completed 79 of 208 simulation runs [21-jul-2022 15:41:46] completed 80 of 208 simulation runs [21-jul-2022 15:41:47] completed 81 of 208 simulation runs [21-jul-2022 15:41:48] completed 82 of 208 simulation runs [21-jul-2022 15:41:49] completed 83 of 208 simulation runs [21-jul-2022 15:41:49] completed 84 of 208 simulation runs [21-jul-2022 15:41:57] completed 85 of 208 simulation runs [21-jul-2022 15:41:58] completed 86 of 208 simulation runs [21-jul-2022 15:41:59] completed 87 of 208 simulation runs [21-jul-2022 15:42:00] completed 88 of 208 simulation runs [21-jul-2022 15:42:01] completed 89 of 208 simulation runs [21-jul-2022 15:42:02] completed 90 of 208 simulation runs [21-jul-2022 15:42:10] completed 91 of 208 simulation runs [21-jul-2022 15:42:11] completed 92 of 208 simulation runs [21-jul-2022 15:42:12] completed 93 of 208 simulation runs [21-jul-2022 15:42:13] completed 94 of 208 simulation runs [21-jul-2022 15:42:14] completed 95 of 208 simulation runs [21-jul-2022 15:42:15] completed 96 of 208 simulation runs [21-jul-2022 15:42:23] completed 97 of 208 simulation runs [21-jul-2022 15:42:23] completed 98 of 208 simulation runs [21-jul-2022 15:42:24] completed 99 of 208 simulation runs [21-jul-2022 15:42:25] completed 100 of 208 simulation runs [21-jul-2022 15:42:26] completed 101 of 208 simulation runs [21-jul-2022 15:42:27] completed 102 of 208 simulation runs [21-jul-2022 15:42:35] completed 103 of 208 simulation runs [21-jul-2022 15:42:36] completed 104 of 208 simulation runs [21-jul-2022 15:42:37] completed 105 of 208 simulation runs [21-jul-2022 15:42:38] completed 106 of 208 simulation runs [21-jul-2022 15:42:38] completed 107 of 208 simulation runs [21-jul-2022 15:42:39] completed 108 of 208 simulation runs [21-jul-2022 15:42:47] completed 109 of 208 simulation runs [21-jul-2022 15:42:48] completed 110 of 208 simulation runs [21-jul-2022 15:42:50] completed 111 of 208 simulation runs [21-jul-2022 15:42:51] completed 112 of 208 simulation runs [21-jul-2022 15:42:52] completed 113 of 208 simulation runs [21-jul-2022 15:42:53] completed 114 of 208 simulation runs [21-jul-2022 15:43:01] completed 115 of 208 simulation runs [21-jul-2022 15:43:02] completed 116 of 208 simulation runs [21-jul-2022 15:43:03] completed 117 of 208 simulation runs [21-jul-2022 15:43:04] completed 118 of 208 simulation runs [21-jul-2022 15:43:05] completed 119 of 208 simulation runs [21-jul-2022 15:43:06] completed 120 of 208 simulation runs [21-jul-2022 15:43:13] completed 121 of 208 simulation runs [21-jul-2022 15:43:14] completed 122 of 208 simulation runs [21-jul-2022 15:43:14] completed 123 of 208 simulation runs [21-jul-2022 15:43:15] completed 124 of 208 simulation runs [21-jul-2022 15:43:15] completed 125 of 208 simulation runs [21-jul-2022 15:43:16] completed 126 of 208 simulation runs [21-jul-2022 15:43:28] completed 127 of 208 simulation runs [21-jul-2022 15:43:29] completed 128 of 208 simulation runs [21-jul-2022 15:43:30] completed 129 of 208 simulation runs [21-jul-2022 15:43:31] completed 130 of 208 simulation runs [21-jul-2022 15:43:32] completed 131 of 208 simulation runs [21-jul-2022 15:43:33] completed 132 of 208 simulation runs [21-jul-2022 15:43:38] completed 133 of 208 simulation runs [21-jul-2022 15:43:40] completed 134 of 208 simulation runs [21-jul-2022 15:43:41] completed 135 of 208 simulation runs [21-jul-2022 15:43:42] completed 136 of 208 simulation runs [21-jul-2022 15:43:43] completed 137 of 208 simulation runs [21-jul-2022 15:43:44] completed 138 of 208 simulation runs [21-jul-2022 15:43:51] completed 139 of 208 simulation runs [21-jul-2022 15:43:52] completed 140 of 208 simulation runs [21-jul-2022 15:43:53] completed 141 of 208 simulation runs [21-jul-2022 15:43:54] completed 142 of 208 simulation runs [21-jul-2022 15:43:55] completed 143 of 208 simulation runs [21-jul-2022 15:43:56] completed 144 of 208 simulation runs [21-jul-2022 15:44:03] completed 145 of 208 simulation runs [21-jul-2022 15:44:05] completed 146 of 208 simulation runs [21-jul-2022 15:44:06] completed 147 of 208 simulation runs [21-jul-2022 15:44:07] completed 148 of 208 simulation runs [21-jul-2022 15:44:08] completed 149 of 208 simulation runs [21-jul-2022 15:44:09] completed 150 of 208 simulation runs [21-jul-2022 15:44:16] completed 151 of 208 simulation runs [21-jul-2022 15:44:17] completed 152 of 208 simulation runs [21-jul-2022 15:44:18] completed 153 of 208 simulation runs [21-jul-2022 15:44:18] completed 154 of 208 simulation runs [21-jul-2022 15:44:19] completed 155 of 208 simulation runs [21-jul-2022 15:44:20] completed 156 of 208 simulation runs [21-jul-2022 15:44:28] completed 157 of 208 simulation runs [21-jul-2022 15:44:29] completed 158 of 208 simulation runs [21-jul-2022 15:44:29] completed 159 of 208 simulation runs [21-jul-2022 15:44:31] completed 160 of 208 simulation runs [21-jul-2022 15:44:32] completed 161 of 208 simulation runs [21-jul-2022 15:44:33] completed 162 of 208 simulation runs [21-jul-2022 15:44:42] completed 163 of 208 simulation runs [21-jul-2022 15:44:42] completed 164 of 208 simulation runs [21-jul-2022 15:44:43] completed 165 of 208 simulation runs [21-jul-2022 15:44:44] completed 166 of 208 simulation runs [21-jul-2022 15:44:45] completed 167 of 208 simulation runs [21-jul-2022 15:44:46] completed 168 of 208 simulation runs [21-jul-2022 15:44:54] completed 169 of 208 simulation runs [21-jul-2022 15:44:55] completed 170 of 208 simulation runs [21-jul-2022 15:44:56] completed 171 of 208 simulation runs [21-jul-2022 15:44:57] completed 172 of 208 simulation runs [21-jul-2022 15:44:57] completed 173 of 208 simulation runs [21-jul-2022 15:44:58] completed 174 of 208 simulation runs [21-jul-2022 15:45:07] completed 175 of 208 simulation runs [21-jul-2022 15:45:08] completed 176 of 208 simulation runs [21-jul-2022 15:45:08] completed 177 of 208 simulation runs [21-jul-2022 15:45:09] completed 178 of 208 simulation runs [21-jul-2022 15:45:10] completed 179 of 208 simulation runs [21-jul-2022 15:45:11] completed 180 of 208 simulation runs [21-jul-2022 15:45:20] completed 181 of 208 simulation runs [21-jul-2022 15:45:21] completed 182 of 208 simulation runs [21-jul-2022 15:45:22] completed 183 of 208 simulation runs [21-jul-2022 15:45:23] completed 184 of 208 simulation runs [21-jul-2022 15:45:24] completed 185 of 208 simulation runs [21-jul-2022 15:45:25] completed 186 of 208 simulation runs [21-jul-2022 15:45:32] completed 187 of 208 simulation runs [21-jul-2022 15:45:33] completed 188 of 208 simulation runs [21-jul-2022 15:45:34] completed 189 of 208 simulation runs [21-jul-2022 15:45:35] completed 190 of 208 simulation runs [21-jul-2022 15:45:36] completed 191 of 208 simulation runs [21-jul-2022 15:45:37] completed 192 of 208 simulation runs [21-jul-2022 15:45:45] completed 193 of 208 simulation runs [21-jul-2022 15:45:45] completed 194 of 208 simulation runs [21-jul-2022 15:45:46] completed 195 of 208 simulation runs [21-jul-2022 15:45:47] completed 196 of 208 simulation runs [21-jul-2022 15:45:48] completed 197 of 208 simulation runs [21-jul-2022 15:45:49] completed 198 of 208 simulation runs [21-jul-2022 15:45:57] completed 199 of 208 simulation runs [21-jul-2022 15:45:58] completed 200 of 208 simulation runs [21-jul-2022 15:45:59] completed 201 of 208 simulation runs [21-jul-2022 15:46:00] completed 202 of 208 simulation runs [21-jul-2022 15:46:01] completed 203 of 208 simulation runs [21-jul-2022 15:46:02] completed 204 of 208 simulation runs [21-jul-2022 15:46:09] completed 205 of 208 simulation runs [21-jul-2022 15:46:10] completed 206 of 208 simulation runs [21-jul-2022 15:46:11] completed 207 of 208 simulation runs [21-jul-2022 15:46:12] completed 208 of 208 simulation runs [21-jul-2022 15:46:12] cleaning up parallel workers...
the generatesimulationensemble ran and logged the simulation results. create a simulation ensemble to process and analyze the simulation results using the simulationensembledatastore command.
ens = simulationensembledatastore(fullfile(pwd,'data'));processing the simulation results
the simulationensembledatastore command created an ensemble object that points to the simulation results. use the ensemble object to prepare and analyze the data in each member of the ensemble. the ensemble object lists the data variables in the ensemble and by default all the variables are selected for reading.
ens
ens =
simulationensembledatastore with properties:
datavariables: [6×1 string]
independentvariables: [0×0 string]
conditionvariables: [0×0 string]
selectedvariables: [6×1 string]
readsize: 1
nummembers: 208
lastmemberread: [0×0 string]
files: [208×1 string]
ens.selectedvariables
ans = 6×1 string
"simulationinput"
"simulationmetadata"
"tacho"
"vibration"
"xfinal"
"xout"
for analysis only read the vibration and tacho signals and the simulink.simulationinput. the simulink.simulationinput has the model variable values used for simulation and is used to create fault labels for the ensemble members. use the ensemble read command to get the ensemble member data.
ens.selectedvariables = ["vibration" "tacho" "simulationinput"]; data = read(ens)
data=1×3 table
vibration tacho simulationinput
___________________ ___________________ ______________________________
{40272×1 timetable} {40272×1 timetable} {1×1 simulink.simulationinput}
extract the vibration signal from the returned data and plot it.
vibration = data.vibration{1};
plot(vibration.time,vibration.data)
title('vibration')
ylabel('acceleration')
the first 10 seconds of the simulation contains data where the transmission system is starting up; for analysis discard this data.
idx = vibration.time >= seconds(10); vibration = vibration(idx,:); vibration.time = vibration.time - vibration.time(1);
the tacho signal contains pulses for the rotations of the drive and load shafts but analysis, and specifically time synchronous averaging, requires the times of shaft rotations. the following code discards the first 10 seconds of the tacho data and finds the shaft rotation times in tachopulses.
tacho = data.tacho{1};
idx = tacho.time >= seconds(10);
tacho = tacho(idx,:);
plot(tacho.time,tacho.data)
title('tacho pulses')
legend('drive shaft','load shaft') % load shaft rotates more slowly than drive shaft
idx = diff(tacho.data(:,2)) > 0.5; tachopulses = tacho.time(find(idx) 1)-tacho.time(1)
tachopulses = 8×1 duration
2.8543 sec
6.6508 sec
10.447 sec
14.244 sec
18.04 sec
21.837 sec
25.634 sec
29.43 sec
the simulink.simulationinput.variables property contains the values of the fault parameters used for the simulation, these values allow us to create fault labels for each ensemble member.
vars = data.simulationinput{1}.variables;
idx = strcmp({vars.name},'sdrift');
if any(idx)
sf = abs(vars(idx).value) > 0.01; % small drift values are not faults
else
sf = false;
end
idx = strcmp({vars.name},'shaftwear');
if any(idx)
sv = vars(idx).value < 0;
else
sv = false;
end
if any(idx)
idx = strcmp({vars.name},'toothfaultgain');
st = abs(vars(idx).value) < 0.1; % small tooth fault values are not faults
else
st = false
end
faultcode = sf 2*sv 4*st; % a fault code to capture different fault conditionsthe processed vibration and tacho signals and the fault labels are added to the ensemble to be used later.
sdata = table({vibration},{tachopulses},sf,sv,st,faultcode, ...
'variablenames',{'vibration','tachopulses','sensordrift','shaftwear','toothfault','faultcode'}) sdata=1×6 table
vibration tachopulses sensordrift shaftwear toothfault faultcode
___________________ ______________ ___________ _________ __________ _________
{30106×1 timetable} {8×1 duration} true false false 1
ens.datavariables = [ens.datavariables; "tachopulses"];the ensemble conditionvariables property can be used to identify the variables in the ensemble that contain condition or fault label data. set the property to contain the newly created fault labels.
ens.conditionvariables = ["sensordrift","shaftwear","toothfault","faultcode"];
the code above was used to process one member of the ensemble. to process all the ensemble members the code above was converted to the function preparedata and using the ensemble hasdata command a loop is used to apply preparedata to all the ensemble members. the ensemble members can be processed in parallel by partitioning the ensemble and processing the ensemble partitions in parallel.
reset(ens) runlocal = false; if runlocal % process each member in the ensemble while hasdata(ens) data = read(ens); adddata = preparedata(data); writetolastmemberread(ens,adddata) end else % split the ensemble into partitions and process each partition in parallel n = numpartitions(ens,gcp); parfor ct = 1:n subens = partition(ens,n,ct); while hasdata(subens) data = read(subens); adddata = preparedata(data); writetolastmemberread(subens,adddata) end end end
plot the vibration signal of every 10th member of the ensemble using the hasdata and read commands to extract the vibration signal.
reset(ens) ens.selectedvariables = "vibration"; figure, ct = 1; while hasdata(ens) data = read(ens); if mod(ct,10) == 0 vibration = data.vibration{1}; plot(vibration.time,vibration.data) hold on end ct = ct 1; end hold off title('vibration signals') ylabel('acceleration')
analyzing the simulation data
now that the data has been cleaned and pre-processed the data is ready for extracting features to determine the features to use to classify the different faults types. first configure the ensemble so that it only returns the processed data.
ens.selectedvariables = ["vibration","tachopulses"];
for each member in the ensemble compute a number of time and spectrum based features. these include signal statistics such as signal mean, variance, peak to peak, non-linear signal characteristics such as approximate entropy and lyapunov exponent, and spectral features such as the peak frequency of the time synchronous average of the vibration signal, and the power of the time synchronous average envelope signal. the analyzedata function contains a full list of the extracted features. by way of example consider computing the spectrum of the time synchronous averaged vibration signal.
reset(ens) data = read(ens)
data=1×2 table
vibration tachopulses
___________________ ______________
{30106×1 timetable} {8×1 duration}
vibration = data.vibration{1};
% interpolate the vibration signal onto periodic time base suitable for fft analysis
np = 2^floor(log(height(vibration))/log(2));
dt = vibration.time(end)/(np-1);
tv = 0:dt:vibration.time(end);
y = retime(vibration,tv,'linear');
% compute the time synchronous average of the vibration signal
tp = seconds(data.tachopulses{1});
vibrationtsa = tsa(y,tp);
figure
plot(vibrationtsa.time,vibrationtsa.tsa)
title('vibration time synchronous average')
ylabel('acceleration')
% compute the spectrum of the time synchronous average
np = numel(vibrationtsa);
f = fft(vibrationtsa.tsa.*hamming(np))/np;
frtsa = f(1:floor(np/2) 1); % tsa frequency response
wtsa = (0:np/2)/np*(2*pi/seconds(dt)); % tsa spectrum frequencies
mtsa = abs(frtsa); % tsa spectrum magnitudes
figure
semilogx(wtsa,20*log10(mtsa))
title('vibration spectrum')
xlabel('rad/s')the frequency that corresponds to the peak magnitude could turn out to be a feature that is useful for classifying the different fault types. the code below computes the features mentioned above for all the ensemble members (running this analysis can take up to an hour on a standard desktop. optional code to run the analysis in parallel using the ensemble partition command is provided.) the names of the features are added to the ensemble data variables property before calling writetolastmemberread to add the computed features to each ensemble member.
reset(ens) ens.datavariables = [ens.datavariables; ... "sigmean";"sigmedian";"sigrms";"sigvar";"sigpeak";"sigpeak2peak";"sigskewness"; ... "sigkurtosis";"sigcrestfactor";"sigmad";"sigrangecumsum";"sigcorrdimension";"sigapproxentropy"; ... "siglyapexponent";"peakfreq";"highfreqpower";"envpower";"peakspeckurtosis"]; if runlocal while hasdata(ens) data = read(ens); adddata = analyzedata(data); writetolastmemberread(ens,adddata); end else % split the ensemble into partitions and analyze each partition in parallel n = numpartitions(ens,gcp); parfor ct = 1:n subens = partition(ens,n,ct); while hasdata(subens) data = read(subens); adddata = analyzedata(data); writetolastmemberread(subens,adddata) end end end
selecting features for fault classification
the features computed above are used to build a classifier to classify the different fault conditions. first configure the ensemble to read only the derived features and the fault labels.
featurevariables = analyzedata('getfeaturenames');
ens.datavariables = [ens.datavariables; featurevariables];
ens.selectedvariables = [featurevariables; ens.conditionvariables];
reset(ens)gather the feature data for all the ensemble members into one table.
featuredata = gather(tall(ens))
consider the sensor drift fault. use the fscnca command with all the features computed above as predictors and the sensor drift fault label (a true false value) as the response. the fscnca command returns weights for each feature and features with higher weights have higher importance in predicting the response. for the sensor drift fault the weights indicate that two features are significant predictors (the signal cumulative sum range and the peak frequency of the spectral kurtosis) and the remaining features have little impact.
idxresponse = strcmp(featuredata.properties.variablenames,'sensordrift'); idxlastfeature = find(idxresponse)-1; % index of last feature to use as a potential predictor featureanalysis = fscnca(featuredata{:,1:idxlastfeature},featuredata.sensordrift); featureanalysis.featureweights idxselectedfeature = featureanalysis.featureweights > 0.1; classifysd = [featuredata(:,idxselectedfeature), featuredata(:,idxresponse)]
a grouped histogram of the cumulative sum range gives us insight into why this feature is a significant predictor for the sensor drift fault.
figure histogram(classifysd.sigrangecumsum(classifysd.sensordrift),'binwidth',5e3) xlabel('signal cumulative sum range') ylabel('count') hold on histogram(classifysd.sigrangecumsum(~classifysd.sensordrift),'binwidth',5e3) hold off legend('sensor drift fault','no sensor drift fault')
the histogram plot shows that the signal cumulative sum range is a good featured for detecting sensor drift faults though an additional feature is probably needed as there may be false positives when the signal cumulative sum range is below 0.5 if just the signal cumulative sum range is used to classify sensor drift.
consider the shaft wear fault. in this case the fscnca function indicates that there are 3 features that are significant predictors for the fault (the signal lyapunov exponent, peak frequency, and the peak frequency in the spectral kurtosis), choose these to classify the shaft wear fault.
idxresponse = strcmp(featuredata.properties.variablenames,'shaftwear');
featureanalysis = fscnca(featuredata{:,1:idxlastfeature},featuredata.shaftwear);
featureanalysis.featureweights
idxselectedfeature = featureanalysis.featureweights > 0.1;
classifysw = [featuredata(:,idxselectedfeature), featuredata(:,idxresponse)]the grouped histogram for the signal lyapunov exponent shows why that feature alone is not a good predictor.
figure histogram(classifysw.siglyapexponent(classifysw.shaftwear)) xlabel('signal lyapunov exponent') ylabel('count') hold on histogram(classifysw.siglyapexponent(~classifysw.shaftwear)) hold off legend('shaft wear fault','no shaft wear fault')
the shaft wear feature selection indicates multiple features are needed to classify the shaft wear fault, the grouped histogram confirms this as the most significant feature (in this case the lyapunov exponent) has a similar distribution for both faulty and fault free scenarios indicating that more features are needed to correctly classify this fault.
finally consider the tooth fault, the fscnca function indicates that there are 3 features primary that are significant predictors for the fault (the signal cumulative sum range, the signal lyapunov exponent and the peak frequency in the spectral kurtosis). choosing those 3 features to classify the tooth fault results in a classifier that has poor performance. instead, use the 6 most important features.
idxresponse = strcmp(featuredata.properties.variablenames,'toothfault'); featureanalysis = fscnca(featuredata{:,1:idxlastfeature},featuredata.toothfault); [~,idxselectedfeature] = sort(featureanalysis.featureweights); classifytf = [featuredata(:,idxselectedfeature(end-5:end)), featuredata(:,idxresponse)] figure histogram(classifytf.sigrangecumsum(classifytf.toothfault)) xlabel('signal cumulative sum range') ylabel('count') hold on histogram(classifytf.sigrangecumsum(~classifytf.toothfault)) hold off legend('gear tooth fault','no gear tooth fault')
using the above results a polynomial svm to classify gear tooth faults. split the feature table into members that are used for training and members for testing and validation. use the training members to create a svm classifier using the fitcsvm command.
rng('default') % for reproducibility cvp = cvpartition(size(classifytf,1),'kfold',5); % randomly partition the data for training and validation classifiertf = fitcsvm(... classifytf(cvp.training(1),1:end-1), ... classifytf(cvp.training(1),end), ... 'kernelfunction','polynomial', ... 'polynomialorder',2, ... 'kernelscale','auto', ... 'boxconstraint',1, ... 'standardize',true, ... 'classnames',[false; true]);
use the classifier to classify the test points using the predict command and check the performance of the predictions using a confusion matrix.
% use the classifier on the test validation data to evaluate performance actualvalue = classifytf{cvp.test(1),end}; predictedvalue = predict(classifiertf, classifytf(cvp.test(1),1:end-1)); % check performance by computing and plotting the confusion matrix confdata = confusionmat(actualvalue,predictedvalue); h = heatmap(confdata, ... 'ylabel', 'actual gear tooth fault', ... 'ydisplaylabels', {'false','true'}, ... 'xlabel', 'predicted gear tooth fault', ... 'xdisplaylabels', {'false','true'}, ... 'colorbarvisible','off');
the confusion matrix indicates that the classifier correctly classifies all non-fault conditions but incorrectly classifies one expected fault condition as not being a fault. increasing the number of features used in the classifier can help improve the performance further.
summary
this example walked through the workflow for generating fault data from simulink, using a simulation ensemble to clean up the simulation data and extract features. the extracted features were then used to build classifiers for the different fault types.
see also
related topics
select a web site
choose a web site to get translated content where available and see local events and offers. based on your location, we recommend that you select: .
you can also select a web site from the following list:
americas
- (español)
- (english)
- (english)
europe
- (english)
- (english)
- (deutsch)
- (español)
- (english)
- (français)
- (english)
- (italiano)
- (english)
- (english)
- (english)
- (deutsch)
- (english)
- (english)
- switzerland
- (english)
asia pacific
- (english)
- (english)
- (english)
- 中国
- (日本語)
- (한국어)