4.3 Results

f ykj<sup>θ</sup> � � <sup>¼</sup> <sup>1</sup>

Fault Detection, Diagnosis and Prognosis

where ζ<sup>k</sup> ¼

m and C.

assumption as

ing code:

code:

22

4.2 Modifying the code

ykζ<sup>k</sup>

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ln 1 þ s=~ykÞ

r h � i

ffiffiffiffiffi <sup>2</sup><sup>π</sup> <sup>p</sup> exp � <sup>1</sup>

Eq. (19) with the likelihood function in Eq. (18).

2

2

ln yk � η<sup>k</sup> ζk � �<sup>2</sup> " #

and <sup>η</sup><sup>k</sup> <sup>¼</sup> ln <sup>~</sup>ykÞ � <sup>0</sup>:<sup>5</sup> <sup>ζ</sup><sup>k</sup> ð Þ<sup>2</sup> �

Also, the prior distribution of each parameter is assumed as a normal distribution as f mð Þ� <sup>N</sup> <sup>4</sup>; <sup>0</sup>:2<sup>2</sup> � �, <sup>f</sup>ð Þ� ln<sup>C</sup> <sup>N</sup> �22:33; <sup>0</sup>:5<sup>2</sup> � �, and f sðÞ� <sup>N</sup> <sup>5</sup> � <sup>10</sup>�4; <sup>1</sup> � <sup>10</sup>�<sup>4</sup> � �<sup>2</sup> � �

f <sup>Θ</sup>ð Þ¼ θ f mð Þ� fð Þ� lnC f sð Þ (19)

Therefore, the joint prior distribution can be obtained from the independence

The posterior distribution can be obtained by multiplying the prior distribution in

For the crack growth example, the code in Appendix needs to be changed as follows. First, the problem definition part in lines 2–15 is replaced with the follow-

%===== PROBLEM DEFINITION 1 (Required Parameters)================== WorkName='Crack'; %work results are saved by WorkName TimeUnit='cycles'; % time unit name time=(0:50:3600)'; % time including both at measurement and at prediction

1.14 1.20 1.22 1.37 1.21 1.25 1.25 1.36 1.30 1.32 1.48 ...

thres=0.043; % threshold - critical value ParamName=['m'; 'C'; 's']; % model parameters' name to be estimated prioDisPar=[4 0.2; -22.33 .5; 5E-4 1E-4]; % parameter prior distributions Ns=10000; % num. of samples for MCMC simulation burnIn=0.2; % burn-in fraction signiLevel=2.5; % significance level for C.I. and P.I. %===============================================================

A total of Ndata ¼ 25 measurement data are provided up to tCUR ¼ 1, 200 cycles, and the degradation is predicted until tend ¼ 3, 600 cycles. Since the prior distributions are assumed a normal distribution, the first column of prioDisPar is the mean,

Next, the model definition part in lines 65–67 is replaced with the following

degraModel=(t.\*exp(C).\*coef.\*(dsig\*sqrt(pi)).^m + a0.^coef).^(1./coef); loca=imag(degraModel)�=0; degraModel(loca)=real(degraModel(loca)); %===============================================================

%===== PROBLEM DEFINITION 2 (model equation)=======================

measuData=[1.03 1.00 0.96 1.14 1.13 1.10 1.15 1.15 1.19 1.19 1.14 ...

1.52 1.47 1.59]'\*0.01;

and the second column is the standard deviation.

a0=0.01; dsig=75; coef=1-m/2;

deviation and mean of a lognormal distribution, respectively. In the above equation, ~ykð Þθ is the model prediction from Eq. (17) at time tk with given model parameters

, k ¼ 1, …, Ndata (18)

are the standard

.

Similar to the battery example, the MATLAB code can be used to plot the trace of MCMC sampling and the histogram of model parameters and RUL. Using a similar code provided in Section 3.4, the degradation trend could also be plotted. For example, Figure 5 shows the predicted degradation with the true degradation. Even if the median (1553 cycles) has a relatively large error with the true RUL (1709 cycles), the 95% confidence interval covers the true RUL.

An important information that was not discussed before is the correlation between model parameters. It is well known that the Paris-Erdogan model parameters, m and C, are strongly correlated [20]. Therefore, it would be beneficial to plot the MCMC samples in the parameter space. The following MATLAB script plots the MCMC samples of the parameters.

```
clear; clc; load('Crack at 1200.mat');
[m, C]=meshgrid((3.4:0.009:4.3)',(-24:0.025:-21.5));
for i=1:101; for j=1:101
  para=[m(i,j); C(i,j); 0.0005];
  [,post(i,j)]=BMappl(para,ParamName,time(1:k1),measuData,prioDisPar);
end; end
plot(samplResul(1,:),samplResul(2,:),'.','color',[0.5 0.5 0.5]); hold on
contour(m,C,real(post));
plot(3.8,log(1.5E-10),'pk','markersize',14);
```
Figure 9 shows the MCMC samples in the parameter space along with the exact value of the parameters (star marker). The figure also shows the contour of the joint posterior PDF based on the grid method. It can be observed that the MCMC samples represent the joint posterior PDF well, and the joint PDF covers the true parameter values. However, due to a strong correlation between the two Paris-Erdogan model parameters, the joint PDF shows a narrow but long tail. Any combination of model parameters along the correlation line may yield a similar damage growth trend.

22 para1(:,1)=para0+weigh.\*(2\*rand(p,1)-1); % sample from proposal dist. 23 [,jPdf1]=BMappl(para1,ParamName,time(1:k1),measuData,prioDisPar); 24 if rand<(jPdf1/jPdf0) && jPdf1>0 % acceptance criterion

Prognostics 102: Efficient Bayesian-Based Prognostics Algorithm in MATLAB

28 sampl(:,i)=para0; % new MCMC sample

30 nBurn=Ns/(1-burnIn)-Ns; % No. of effective MCMC samples 31 samplResul=sampl(:,nBurn+1:end); % Final Sampling results

33 for k=1:length(time(k1:end)) % degradation prediction 34 [degrPreCon(k,:),]=BMappl(samplResul,ParamName,time(k1-1+k),[],[]);

37 for i=1:Ns % RUL prediction 38 RUL(i)=interp1(degrPreCon(:,i),time(k1:end),thres,'pchip') - time(k1);

42 perceValue=[50 signiLevel 100-signiLevel]; % median & confi-

43 rulPerce=prctile(RUL,perceValue); %percentiles of RUL

51 figure(2); set(gca,'fontsize',14); hist(RUL,30); % RUL histogram

53 titleName=['at ' num2str(time(k1)) ' ' TimeUnit]; title(titleName) 54 fprintf( '\n # Percentiles of RUL at %g cycles \n', time(k1))

55 fprintf('\n %gprct: %g, median: %g, %gprct: %g \n' , perceValue(2), ... 56 rulPerce(2), rulPerce(1), perceValue(3), rulPerce(3))

57 Name=[WorkName ' at ' num2str(time(k1)) '.mat']; save(Name); % save work

60 function [degraModel, poste]=BMappl(param,ParamName,t,measuData,

65 %===== PROBLEM DEFINITION 2 (model equation)====================

67 %===========================================================

71 prior=prod(unifpdf(param,prioDisPar(:,1),prioDisPar(:,2))); % prior 72 likel=(1./s).^length(measuData) ... % likelihood

52 xlim([min(RUL) max(RUL)]); xlabel(['RUL' ' (' TimeUnit ')']);

61 % Evaluate the degradation model or posterior PDF

35 degraPredi(k,:)=degrPreCon(k,:)+normrnd(0,samplResul(end,:));

25 para0=para1; 26 jPdf0=jPdf1;

32 %%% RUL prediction

DOI: http://dx.doi.org/10.5772/intechopen.82781

40 %%% POST-PROCESSING

47 plot(samplResul(j,:)); 48 ylabel(ParamName(j,:)); 49 title('MCMC sample trace');

41 Index=isnan(RUL); RUL(Index)=[];

45 for j=1:p % plotting MCMC sample trace 46 subplot(p,1,j); % for all model parameters

27 end

29 end

36 end

39 end

dence intervals

44 figure(1);

50 end

58 end 59 %

64 end

prioDisPar)

69 poste=0; 70 else

25

62 for j=1:size(param,1)

66 degraModel=exp(-b.\*t);

68 if isempty(measuData)

63 eval([ParamName(j,:) '=param(j,:);']);

Figure 9. Correlation between the two Paris-Erdogan model parameters.
