Improving the results from SVM classification of faults (2025)

mdl = 'Boost_converter';

open_system(mdl)

% Set fault parameters

%**************************************************************************

% Default parameters for simulation

% Bondwire RC 1e-4

% ESR is 1e-6

% Thermal Resistance 5th element in Foster Model [9.5e-3 0.125039

% 0.132857 0.256654 "0.021551" 2.1e-3] - 0.021551 relates to solder between

% layers

%**************************************************************************

%Bondwire RC failire at 3e-3

Bondwire_fault_set=linspace(1e-4,3e-3,10);

% Electrolytic Capacitor 220uF failure at ESR of 0.5 Ohms

ESR_fault_set=linspace(1e-6,0.5,10);

% Foster Thermal Resistance paramters, failure at 0.2 increase of 5th

% element in Foster model

Thermal_Resistance_set=linspace(0.021551,0.2,10);

IGBT_RC=[mdl '/IGBT Full1']; 'RC';

%**************************************************************************

% set number of simulation runs here default is 20, this gives a data set

% of 80.

nPerGroup = 20; % 125 Number of elements in each fault group

%**************************************************************************

% No fault simulations

Bondwire = repmat(Bondwire_fault_set(1),nPerGroup,1);

ESR = repmat(ESR_fault_set(1),nPerGroup,1);

Thermal_Resistance = repmat(Thermal_Resistance_set(1),nPerGroup,1);

% Single fault simulations

idx = ceil(10*rand(nPerGroup,1));

Bondwire = [Bondwire; Bondwire_fault_set(idx)'];

ESR = [ESR;repmat(ESR_fault_set(1),nPerGroup,1)];

Thermal_Resistance= [Thermal_Resistance;repmat(Thermal_Resistance_set(1),nPerGroup,1)];

idx = ceil(10*rand(nPerGroup,1));

Bondwire = [Bondwire; repmat(Bondwire_fault_set(1),nPerGroup,1)];

ESR = [ESR;ESR_fault_set(idx)'];

Thermal_Resistance = [Thermal_Resistance;repmat(Thermal_Resistance_set(1),nPerGroup,1)];

idx = ceil(10*rand(nPerGroup,1));

Bondwire = [Bondwire; repmat(Bondwire_fault_set(1),nPerGroup,1)];

ESR = [ESR;repmat(ESR_fault_set(1),nPerGroup,1)];

Thermal_Resistance = [Thermal_Resistance;Thermal_Resistance_set(idx)'];

%

% %

% % %**************************************************************************

% % % Double fault simulations

% % idxA = ceil(10*rand(nPerGroup,1));

% % idxB = ceil(10*rand(nPerGroup,1));

% % Bondwire = [Bondwire; Bondwire_fault_set(idxA)'];

% % ESR = [ESR;ESR_fault_set(idxB)'];

% % Thermal_Resistance = [Thermal_Resistance;repmat(Thermal_Resistance_set(1),nPerGroup,1)];

% % idxA = ceil(10*rand(nPerGroup,1));

% % idxB = ceil(10*rand(nPerGroup,1));

% % Bondwire = [Bondwire; Bondwire_fault_set(idxA)'];

% % ESR = [ESR;repmat(ESR_fault_set(1),nPerGroup,1)];

% % Thermal_Resistance = [Thermal_Resistance;Thermal_Resistance_set(idxB)'];

% % idxA = ceil(10*rand(nPerGroup,1));

% % idxB = ceil(10*rand(nPerGroup,1));

% % Bondwire = [Bondwire; repmat(Bondwire_fault_set(1),nPerGroup,1)];

% % ESR = [ESR;ESR_fault_set(idxA)'];

% % Thermal_Resistance = [Thermal_Resistance;Thermal_Resistance_set(idxB)'];

% % %**************************************************************************

% % % Triple fault simulations

% % idxA = ceil(10*rand(nPerGroup,1));

% % idxB = ceil(10*rand(nPerGroup,1));

% % idxC = ceil(10*rand(nPerGroup,1));

% % Bondwire = [Bondwire; Bondwire_fault_set(idxA)'];

% % ESR = [ESR; ESR_fault_set(idxB)'];

% % Thermal_Resistance = [Thermal_Resistance; Thermal_Resistance_set(idxC)'];

% %**************************************************************************

%Use the fault parameter combinations to create Simulink.SimulationInput objects.

%For each simulation input ensure that the random seed is set differently to generate different results.

IGBT=[mdl '/IGBT Full1'];'RC';

CAP=[mdl '/Capacitor'];'r';

IGBT2=[mdl '/IGBT Full1'];'Rth_Foster(5)';

Noise_Source=[mdl '/Voltage Source'];'user_seed';

for ct = numel(Bondwire):-1:1

simInput(ct) = Simulink.SimulationInput(mdl);

simInput(ct) = setVariable(simInput(ct),'IGBT',Bondwire(ct));

simInput(ct) = setVariable(simInput(ct),'CAP',ESR(ct));

simInput(ct) = setVariable(simInput(ct),'IGBT2',Thermal_Resistance(ct));

simInput(ct) = setVariable(simInput(ct),'Noise_Source',ct-1);

end

%Use the generateSimulationEnsemble function to run the simulations defined by the Simulink.

%SimulationInput objects defined above and store the results in a local sub-folder.

%Then create a simulationEnsembleDatastore from the stored results.

% Run the simulation 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

% [ok,e] = generateSimulationEnsemble(simInput(1:10),fullfile(pwd,'Data'),'UseParallel',true);

% [ok,e] = generateSimulationEnsemble(simInput(1:80),fullfile('.','Data'),'UseParallel',true);

% [ok,e] = generateSimulationEnsemble(simInput,fullfile('.','Data'),'UseParallel',true);

% Use above for full data set

runAll = true;

if runAll

[ok,e] = generateSimulationEnsemble(simInput,fullfile('.','Data'),'UseParallel',true);

else

[ok,e] = generateSimulationEnsemble(simInput(1:10),fullfile('.','Data')); %#ok<UNRCH>

end

%**************************************************************************

% Important - Run from here once data has been generated!!!

% *************************************************************************

ens = simulationEnsembleDatastore(fullfile(pwd,'Data'));

% %

ens.DataVariables

% Configure the ensemble so that the read only returns the variables of interest

% Call the preprocess function

ens.SelectedVariables = ["Vce", "SimulationInput"];

reset(ens);

data = read(ens);

[Vce,VceSpectrum,VceFrequencies,faultValues] = preprocess2(data);

%

% %Plot the voltage and voltage spectrum. The plotted data is for a fault free condition.

% % Figure with nominal

subplot(211);

plot(Vce.Time,Vce.Data);

subplot(212);

semilogx(VceFrequencies,pow2db(VceSpectrum));

xlabel('Hz')

%Configure the ensemble with data variables for the condition indicators and condition variables for fault variable values.

%Then call the extractCI function to compute the features, and use the writeToLastMemberRead command to add the feature and fault variable values to the ensemble.

%The extractCI function is separate.

%

ens.DataVariables = [ens.DataVariables; ...

"fPeak"; "pLow"; "pMid"; "pHigh"; "pKurtosis"; ...

"qMean"; "qVar"; "qSkewness"; "qKurtosis"; ...

"qPeak2Peak"; "qCrest"; "qRMS"; "qMAD"; "qCSRange"];

ens.ConditionVariables = ["Bondwire_Fault","ESR_Fault","Thermal_Fault"];

%

feat = extractCI2(Vce,VceSpectrum,VceFrequencies);

dataToWrite = [faultValues, feat];

writeToLastMemberRead(ens,dataToWrite{:})

%

%The above code preprocesses and computes the condition indicators for the first member of the simulation ensemble.

%Repeat this for all the members in the ensemble using the ensemble hasdata command.

%To get an idea of the simulation results under different fault conditions plot every hundredth element of the ensemble.

%

% %Figure with nominal and faults

figure,

subplot(211);

lVce = plot(Vce.Time,Vce.Data,'LineWidth',2);

subplot(212);

lVceSpectrum = semilogx(VceFrequencies,pow2db(VceSpectrum),'LineWidth',2);

xlabel('Hz')

ct = 1;

lColors = get(lVce.Parent,'ColorOrder');

lIdx = 2;

%

% % Loop over all members in the ensemble, preprocess

% % and compute the features for each member

while hasdata(ens)

%

% % Read member data

data = read(ens);

%

% % Preprocess and extract features from the member data

[Vce,VceSpectrum,VceFrequencies,faultValues] = preprocess2(data);

feat = extractCI2(Vce,VceSpectrum,VceFrequencies);

% % Add the extracted feature values to the member data

dataToWrite = [faultValues, feat];

writeToLastMemberRead(ens,dataToWrite{:})

%

% % Plot member signal and spectrum for every 10th member

if mod(ct,10) == 0

line('Parent',lVce.Parent,'XData',Vce.Time,'YData',Vce.Data, ...

'Color', lColors(lIdx,:));

line('Parent',lVceSpectrum.Parent,'XData',VceFrequencies,'YData',pow2db(VceSpectrum), ...

'Color', lColors(lIdx,:));

if lIdx == size(lColors,1)

lIdx = 1;

else

lIdx = lIdx+1;

end

end

ct = ct + 1;

end

% %Configure the simulation ensemble to read the condition indicators, and

% %use the tall and gather commands to load all the condition indicators and fault variable values into memory

%

% % Get data to design a classifier.

reset(ens)

ens.SelectedVariables = ["fPeak","pLow","pMid","pHigh","pKurtosis",...

"qMean","qVar","qSkewness","qKurtosis",...

"qPeak2Peak","qCrest","qRMS","qMAD","qCSRange",...

"Bondwire_Fault","ESR_Fault","Thermal_Fault"];

idxLastFeature = 14;

%

% % Load the condition indicator data into memory

data = gather(tall(ens));

data(1:10,:)

%data(:,:)

%

% % Convert the fault variable values to flags

data.BondwireFlag = data.Bondwire_Fault > 1e-4;

data.ESRFlag = data.ESR_Fault > 1e-6;

data.ThermalFlag = data.Thermal_Fault > 0.021551;

data.CombinedFlag =data.BondwireFlag+2*data.ESRFlag+4*data.ThermalFlag;

% %Create a classifier that takes as input the condition indicators and returns the combined fault flag.

% %Train a support vector machine that uses a 2nd order polynomial kernel.

% %Use the cvpartition command to partition the ensemble members into a set for training and a set for validation.

%

rng('default') % for reproducibility

predictors = data(:,1:idxLastFeature);

response = data.CombinedFlag;

cvp = cvpartition(size(predictors,1),'KFold',5);

%

% % Create and train the classifier

template = templateSVM(...

'KernelFunction', 'polynomial', ...

'PolynomialOrder', 2, ...

'KernelScale', 'auto', ...

'BoxConstraint', 1, ...

'Standardize', true);

combinedClassifier = fitcecoc(...

predictors(cvp.training(1),:), ...

response(cvp.training(1),:), ...

'Learners', template, ...

'Coding', 'onevsone', ...

'ClassNames', [0; 1; 2; 3; 4; 5; 6; 7]);

% %Check the performance of the trained classifier using the validation data and plot the results on a confusion plot.

%

% % Check performance by computing and plotting the confusion matrix

actualValue = response(cvp.test(1),:);

predictedValue = predict(combinedClassifier, predictors(cvp.test(1),:));

confdata = confusionmat(actualValue,predictedValue);

figure,

labels = {'None','Capacitor','Thermal','Bondwire','Bondwire & Capacitor', ...

'Thermal & Bondwire','Thermal & Capacitor','All'};

h = heatmap(confdata)

% ,...

% 'YLabel', 'Actual IGBT fault', ...

% 'YDisplayLabels', labels, ...

% 'XLabel', 'Predicted IGBT fault', ...

% 'XDisplayLabels', labels, ...

% 'ColorbarVisible','off');

%

% % Compute the overall accuracy of the classifier

sum(diag(confdata))/sum(confdata(:))

%

% % Compute the accuracy of the classifier at predicting

% % that there is a fault

1-sum(confdata(2:end,1))/sum(confdata(:))

function [Vce,VceSpectrum,VceFrequencies,faultValues] = preprocess2(data)

% Remove the 1st 0.1 seconds of the Voltage signal to eliminate transient

tMin = seconds(0.1);

Vce = data.Vce{1};

Vce = Vce(Vce.Time >= tMin,:);

Vce.Time = Vce.Time - Vce.Time(1);

% Ensure the Voltage is sampled at a uniform sample rate

Vce = retime(Vce,'regular','linear','TimeStep',seconds(5e-5)); %5e-5 is 20KHz sample rate

% Remove the mean from the Vce and compute the Vce spectrum

fA = Vce;

fA.Data = fA.Data-mean(fA.Data);

%

[VceSpectrum,VceFrequencies] = pspectrum(fA);

% 'FrequencyLimits',[2 2500]

% % Find the values of the fault variables from the SimulationInput

simin = data.SimulationInput{1};

vars = {simin.Variables.Name};

idx = strcmp(vars,'IGBT');

Bondwire_Fault = simin.Variables(idx).Value;

idx = strcmp(vars,'CAP');

ESR_Fault = simin.Variables(idx).Value;

idx = strcmp(vars,'IGBT2');

Thermal_Fault = simin.Variables(idx).Value;

%

% % Collect the fault values in a cell array

faultValues = {...

'Bondwire_Fault', Bondwire_Fault, ...

'ESR_Fault', ESR_Fault, ...

'Thermal_Fault', Thermal_Fault};

end

function ci = extractCI2(Vce,VceSpectrum,VceFrequencies)

% % Find the frequency of the peak magnitude in the power spectrum.

pMax = max(VceSpectrum);

fPeak = VceFrequencies(VceSpectrum==pMax);

% % Compute the power in the low frequency range Hz.

fRange = VceFrequencies >=10 & VceFrequencies <= 100;

pLow = sum(VceSpectrum(fRange));

% % Compute the power in the mid frequency range Hz.

fRange = VceFrequencies >= 100 & VceFrequencies <= 600;

pMid = sum(VceSpectrum(fRange));

%

% % Compute the power in the high frequency range Hz.

fRange = VceFrequencies >= 600;

pHigh = sum(VceSpectrum(fRange));

%

% % Find the frequency of the spectral kurtosis peak

[pKur,fKur] = pkurtosis(VceSpectrum);

pKur = fKur(pKur == max(pKur));

%

% % Compute the Vce cumulative sum range.

csVce = cumsum(Vce.Data);

csVceRange = max(csVce)-min(csVce);

%

% % Collect the feature and feature values in a cell array.

% % Add flow statistic (mean, variance, etc.) and common signal

% % characteristics (rms, peak2peak, etc.) to the cell array.

ci = {...

'qMean', mean(Vce.Data), ...

'qVar', var(Vce.Data), ...

'qSkewness', skewness(Vce.Data), ...

'qKurtosis', kurtosis(Vce.Data), ...

'qPeak2Peak', peak2peak(Vce.Data), ...

'qCrest', peak2rms(Vce.Data), ...

'qRMS', rms(Vce.Data), ...

'qMAD', mad(Vce.Data), ...

'qCSRange',csVceRange,...

'fPeak', fPeak, ...

'pLow', pLow, ...

'pMid', pMid, ...

'pHigh', pHigh, ...

'pKurtosis', pKur(1)};

end