Percentiles of RUL at 9 cycles

During MCMC sampling, the number of samples Ns is set to 5000 (line 12). Since the sampling process takes many iterations to converge, the initial 20% of the samples are discarded in calculating the posterior distribution by command burnIn = 0.2 (line 13). In order to keep 5000 samples, Ns/(1-burnIn) = 6250 samples (line 21) are

and upper-bounds of the distribution.

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

Following is the sample output from BM.m:

(line 62–64):

end

17

for j=1:size(param,1)

b = param(1,:); s = param(2,:);

given time t (line 66).

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

is used for calculating the acceptance ratio (line 23).

time and model parameters can be an array of samples.

3.2 Bayesian parameter estimation with MCMC (lines 16–31)

In the case of the battery example, this command is equivalent to.

It is expected that the MATLAB script is saved as a file with the name of 'BM.m', which has two input arguments, para0 and weigh (line 1). The first argument, para0, is the initial sample of model parameters, and weigh is the width w of the proposal distribution in Eq. (12). The size of each array should be the same as the number of model parameters. The following is an example of calling the code in the command window of MATLAB:

#### samplResul=BM([0.011 0.02]', [0.001 0.003]');

In the above MATLAB command, para0 = [0.011 0.02]' is the initial sample of b and s, and weigh = [0.001 0.003]' is the width of proposal distribution of b and s. Since the convergence and accuracy depend on these two variables, it is suggested to try with different values. Since the users know the prior distribution of the model parameters, it is a good practice to start with the mean of the prior distribution as an initial sample. If the code ran successfully, it will return the samples of model parameters in the samplResul array and will generate two plots. The first plot is the trace of MCMC samples, and the second plot is the histogram of the RUL.

#### 3.1 Problem definition (lines 2–15, 65–67)

The problem definition means defining the degradation equation using model parameters and time/cycle. All known parameters, as well as the initial estimate of unknown parameters, parameter names, and model data, need to be defined. The problem definition consists of two parts: parameter definition and model definition. For parameter definition (lines 2–15), 'Battery' is used as a WorkName (line 3). The capacity is measured every week, so TimeUnit is 'weeks' (line 4). In line 5, time is an array of discrete times in the units of TimeUnit. Measurements and predictions will be done at these times. The relative C/1 capacity data is stored as measuData (lines 6–7), which has 10 weeks of measurement; that is, Ndata ¼ 10 (k1 in line 18). These data correspond to the first 10 times in time array. Since time starts from 0, the 10th time corresponds to 9 weeks, which is the current time; that is, tCUR ¼ 9week. The degradation will be predicted for future times that are from 10 to 50 weeks. The measurement data are generated using the MATLAB script in Section 2.1. It is noted that due to random noise, the users may experience different realizations of measurement data.

In line 8, the degradation threshold ythreshold is defined using variable 'thres' with the value of 0.7. ParamName in line 9 is the name of model parameters: model parameter 'b' and standard deviation of noise 's'. Since the parameter name will actually be used in the model, it is important to use the actual name of variables here. The number of parameters is stored in the variable 'p' in line 17. It is noted that the last parameter should be the standard deviation of noise 's'. prioDisPar stores the information of the prior distribution of model parameters, as given in

Prognostics 102: Efficient Bayesian-Based Prognostics Algorithm in MATLAB DOI: http://dx.doi.org/10.5772/intechopen.82781

Eq. (9). When a uniform prior distribution is used, each row contains the lowerand upper-bounds of the distribution.

During MCMC sampling, the number of samples Ns is set to 5000 (line 12). Since the sampling process takes many iterations to converge, the initial 20% of the samples are discarded in calculating the posterior distribution by command burnIn = 0.2 (line 13). In order to keep 5000 samples, Ns/(1-burnIn) = 6250 samples (line 21) are generated first, and then, nBurn = 1250 samples (line 30) are discarded.

Since the RUL is represented by Ns samples, the confidence interval or the prediction interval is often used to support the decision-making process. signiLevel (line 14) is the significance level of this interval in percentage. When signiLevel = 5, the code will return the lower 5 percentile, median, and the upper 5 percentile. Following is the sample output from BM.m:

> # Percentiles of RUL at 9 cycles 5prct: 18.7182, median: 20.381, 95prct: 22.1576

For the model definition, the degradation model ~y tð Þ ; b in Eq. (1) is defined in lines 65–67. In this equation, t is the time of measurement, and b is the model parameter as defined in line 9. The model equation needs to be defined in such a way that component-by-component operations are possible. This is because the time and model parameters can be an array of samples.

#### 3.2 Bayesian parameter estimation with MCMC (lines 16–31)

In the Bayesian parameter estimation process, the posterior distribution is expressed in terms of the product of the prior distribution and the likelihoods of all measured data. In the MATLAB code BM.m, the degradation model and the posterior distribution are calculated in the function BMappl (lines 60–76). First, the parameter samples in param are assigned to the variables using the eval function (line 62–64):

for j=1:size(param,1) eval([ParamName(j,:) '=param(j,:);']); end

In the case of the battery example, this command is equivalent to.

b = param(1,:);

3. MATLAB implementation

Fault Detection, Diagnosis and Prognosis

command window of MATLAB:

3.1 Problem definition (lines 2–15, 65–67)

realizations of measurement data.

16

quent sections with an example of battery degradation.

In this section, MATLAB implementation of prognostics using the Bayesian method is discussed. In the following explanation, 'line' or 'lines' in a parenthesis indicated the number of the line of the code in Appendix. The code is divided into three parts: (1) Problem Definition (lines 2–15, 65–67) (2) Bayesian parameter estimation and MCMC sampling (lines 16–29, 60–76) (3) Post-processing for displaying results (lines 40–57). Only the Problem Definition part needs to be changed for different applications. Detailed explanations are given in the subse-

It is expected that the MATLAB script is saved as a file with the name of 'BM.m',

which has two input arguments, para0 and weigh (line 1). The first argument, para0, is the initial sample of model parameters, and weigh is the width w of the proposal distribution in Eq. (12). The size of each array should be the same as the number of model parameters. The following is an example of calling the code in the

samplResul=BM([0.011 0.02]', [0.001 0.003]');

trace of MCMC samples, and the second plot is the histogram of the RUL.

In the above MATLAB command, para0 = [0.011 0.02]' is the initial sample of b and s, and weigh = [0.001 0.003]' is the width of proposal distribution of b and s. Since the convergence and accuracy depend on these two variables, it is suggested to try with different values. Since the users know the prior distribution of the model parameters, it is a good practice to start with the mean of the prior distribution as an initial sample. If the code ran successfully, it will return the samples of model parameters in the samplResul array and will generate two plots. The first plot is the

The problem definition means defining the degradation equation using model parameters and time/cycle. All known parameters, as well as the initial estimate of unknown parameters, parameter names, and model data, need to be defined. The problem definition consists of two parts: parameter definition and model definition. For parameter definition (lines 2–15), 'Battery' is used as a WorkName (line 3). The capacity is measured every week, so TimeUnit is 'weeks' (line 4). In line 5, time is an array of discrete times in the units of TimeUnit. Measurements and predictions will be done at these times. The relative C/1 capacity data is stored as measuData (lines 6–7), which has 10 weeks of measurement; that is, Ndata ¼ 10 (k1 in line 18). These data correspond to the first 10 times in time array. Since time starts from 0, the 10th time corresponds to 9 weeks, which is the current time; that is, tCUR ¼ 9week. The degradation will be predicted for future times that are from 10 to 50 weeks. The measurement data are generated using the MATLAB script in Section 2.1. It is noted that due to random noise, the users may experience different

In line 8, the degradation threshold ythreshold is defined using variable 'thres' with

the value of 0.7. ParamName in line 9 is the name of model parameters: model parameter 'b' and standard deviation of noise 's'. Since the parameter name will actually be used in the model, it is important to use the actual name of variables here. The number of parameters is stored in the variable 'p' in line 17. It is noted that the last parameter should be the standard deviation of noise 's'. prioDisPar stores the information of the prior distribution of model parameters, as given in

```
s = param(2,:);
```
This is why the ParamName in line 9 must have the same name with the actual variable. Then, these parameters are used to calculate the degradation model with given time t (line 66).

If measurement data (measuData) is empty, then BMappl only calculates the degradation model with given parameter samples at a given time t. This corresponds to propagating the degradation model using parameter samples to the future time for prognostics. If measurement data are provided (lines 71–74), then BMappl calculates the values of the posterior distribution at the parameter samples. In this case, the time should be given as an array of Ndata components from the start time to the current time. The posterior distribution in Eq. (11) (line 74) is the multiplication of the prior distribution (line 71) in Eq. (9) with the likelihood of all measured data (lines 72–73) in Eq. (10). The calculated posterior distribution from BMappl is used for calculating the acceptance ratio (line 23).

MCMC sampling using the MH algorithm starts with the initial sample that is provided by the users (line 19) and the value of the posterior PDF (line 20). In the loop of MCMC sampling (lines 21–29), a candidate sample is randomly generated from a uniform distribution, centered at the current sample and the width of �w (line 22). The value of posterior PDF for the candidate sample is also calculated (line 23). If the acceptance ratio in Eq. (13) is greater than a randomly generated number u, then the candidate sample is accepted as a new sample (lines 25–26). Otherwise, the current sample and its posterior PDF are kept as a new sample and PDF. Once the MCMC sampling loop is over, the first 20% of the samples are discarded as a burn-in process (line 31). At the end of Bayesian parameter estimation, samplResul array contains Ns samples of model parameters.

### 3.3 Remaining useful life prediction (lines 32–39)

Once the samples of model parameters are obtained based on the posterior distribution, they can be used to find the RUL, which is the time when the degradation prediction reaches the threshold. First, for all samples in samplResul, the degradation is predicted in the future times between tCUR ¼ k1 and tend.

for k=1:length(time(k1:end)) [degrPreCon(k,:),�]=BMappl(samplResul,ParamName,time(k1-1+k),[],[]); end

Once the degradations in the future times are calculated, MATLAB function interp1 is used to find the time when ~y tð Þ¼ EOL; b ythreshold. Once tEOL is found, the RUL can be calculated using Eq. (3).

```
for i=1:Ns
 RUL(i)=interp1(degrPreCon(:,i),time(k1:end),thres,'pchip') - time(k1);
end
```
The option 'pchip' in interp1 uses a shape-preserving piecewise cubic interpolation, which preserves C1 -continuity.

### 3.4 Postprocessing (lines 40–58)

In the postprocessing stage, the results given in samples are interpreted in terms of statistical quantities or in the form of graphs. First, in the RUL array, those components that have an infinite life should be removed (line 41). Then, the confidence intervals of [5%, median, 95%] are calculated from the RUL array and stored in rulPerce (line 42–43).

The MATLAB code plots two figures. The first figure plots the trace of MCMC samples (lines 45–50) as shown in Figure 6(a). This trace shows the quality of MCMC samples. If samples are distributed randomly and symmetrically around the mean with a constant bound, then it means that the samples are stabilized and well represent the posterior distribution. If the trace of samples shows an irregular behavior as shown in Figure 6(b), the samples may not represent the posterior PDF properly. This can happen when the initial sample was far away from the mean and when the width of the proposal distribution is too narrow or too wide.

The second plot is the histogram of RUL (lines 51–53) as shown in Figure 7. At the end of the code, the confidence intervals of [5%, median, 95%] are printed on the command window (lines 54–56). All variables are saved in the computer file so that they can be loaded to the memory for further analysis (line 57). The name of the saved database is "WorkName at tCUR.mat". For example, in the battery case, the saved database name is "Battery at 9.mat".

Although the MATLAB code plots two figures, it is possible that the users can plot different figures using the saved database. After calling the BM.m function, the saved database has to be loaded to the memory using the following commands:

Trace of samples from MCMC sampling (a) proper samples and (b) improper samples.

Prognostics 102: Efficient Bayesian-Based Prognostics Algorithm in MATLAB

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

clear; clc; load('Battery at 9.mat')

Histogram of the remaining useful life.

Figure 6.

Figure 7.

19

Prognostics 102: Efficient Bayesian-Based Prognostics Algorithm in MATLAB DOI: http://dx.doi.org/10.5772/intechopen.82781

Figure 6.

loop of MCMC sampling (lines 21–29), a candidate sample is randomly generated from a uniform distribution, centered at the current sample and the width of �w (line 22). The value of posterior PDF for the candidate sample is also calculated (line 23). If the acceptance ratio in Eq. (13) is greater than a randomly generated number u, then the candidate sample is accepted as a new sample (lines 25–26). Otherwise, the current sample and its posterior PDF are kept as a new sample and PDF. Once the MCMC sampling loop is over, the first 20% of the samples are discarded as a burn-in process (line 31). At the end of Bayesian parameter estimation, samplResul

Once the samples of model parameters are obtained based on the posterior distribution, they can be used to find the RUL, which is the time when the degradation prediction reaches the threshold. First, for all samples in samplResul, the

[degrPreCon(k,:),�]=BMappl(samplResul,ParamName,time(k1-1+k),[],[]);

Once the degradations in the future times are calculated, MATLAB function interp1 is used to find the time when ~y tð Þ¼ EOL; b ythreshold. Once tEOL is found, the

RUL(i)=interp1(degrPreCon(:,i),time(k1:end),thres,'pchip') - time(k1);


The option 'pchip' in interp1 uses a shape-preserving piecewise cubic interpola-

In the postprocessing stage, the results given in samples are interpreted in terms

of statistical quantities or in the form of graphs. First, in the RUL array, those components that have an infinite life should be removed (line 41). Then, the confidence intervals of [5%, median, 95%] are calculated from the RUL array and

The MATLAB code plots two figures. The first figure plots the trace of MCMC samples (lines 45–50) as shown in Figure 6(a). This trace shows the quality of MCMC samples. If samples are distributed randomly and symmetrically around the mean with a constant bound, then it means that the samples are stabilized and well represent the posterior distribution. If the trace of samples shows an irregular behavior as shown in Figure 6(b), the samples may not

represent the posterior PDF properly. This can happen when the initial sample was far away from the mean and when the width of the proposal distribution is too

The second plot is the histogram of RUL (lines 51–53) as shown in Figure 7. At the end of the code, the confidence intervals of [5%, median, 95%] are printed on the command window (lines 54–56). All variables are saved in the computer file so that they can be loaded to the memory for further analysis (line 57). The name of the saved database is "WorkName at tCUR.mat". For example, in the battery case,

degradation is predicted in the future times between tCUR ¼ k1 and tend.

array contains Ns samples of model parameters.

for k=1:length(time(k1:end))

Fault Detection, Diagnosis and Prognosis

RUL can be calculated using Eq. (3).

3.4 Postprocessing (lines 40–58)

stored in rulPerce (line 42–43).

the saved database name is "Battery at 9.mat".

narrow or too wide.

18

end

end

for i=1:Ns

tion, which preserves C1

3.3 Remaining useful life prediction (lines 32–39)

Trace of samples from MCMC sampling (a) proper samples and (b) improper samples.

Figure 7. Histogram of the remaining useful life.

Although the MATLAB code plots two figures, it is possible that the users can plot different figures using the saved database. After calling the BM.m function, the saved database has to be loaded to the memory using the following commands:

clear; clc; load('Battery at 9.mat')

In addition to the trace of samples shown in Figure 6, it is possible to plot the histogram of model parameters using the same samples. The following commands plot the histograms of all model parameters.

20.38 weeks, which is close to the true RUL. In addition, the 90% confidence interval is about 3.4 weeks; that is, the uncertainty in the prediction is about 17%.

Prognostics 102: Efficient Bayesian-Based Prognostics Algorithm in MATLAB

The code can be modified easily by the users for various applications. In this section, an example of crack growth is used to explain how the MATLAB code BM.m

In fatigue crack growth, the failure criterion is given in terms of the crack size. Therefore, it would be appropriate to use the size of crack as a degradation feature. In this case, the degradation feature monotonically increases, while it was monotonically decreased for the battery example. Assuming that a through-the-thickness center crack exists in an infinite plate under mode-I loading condition, the rate of

<sup>d</sup><sup>N</sup> <sup>¼</sup> <sup>C</sup>ð Þ <sup>Δ</sup><sup>K</sup> <sup>m</sup> (16)

<sup>1</sup>�<sup>m</sup> 2 0

<sup>2</sup> ð Þ, s <sup>¼</sup> <sup>0</sup>:0005m to the true crack sizes. The

2�m

(17)

<sup>π</sup><sup>a</sup> <sup>p</sup> is the range of stress intensity factor, and <sup>Δ</sup><sup>σ</sup> is the stress

<sup>π</sup> � � <sup>p</sup> <sup>m</sup> <sup>þ</sup> <sup>a</sup>

h i <sup>2</sup>

fatigue crack growth can be expressed using the Paris-Erdogan model as

da

where a is the half crack size, N is the number of cycles, m and C are model

range. It is assumed that time is the number of fatigue loading cycles. For the consistent notation, the crack size and cycles are replaced with ~y ¼ a and t ¼ N in the following explanation. Since the degradation model requires the crack size as a function of time and model parameters, Eq. (16) can be integrated to obtain the

> 2 � � Δσ ffiffiffi

<sup>m</sup>true <sup>¼</sup> <sup>3</sup>:8 and <sup>C</sup>true <sup>¼</sup> <sup>1</sup>:<sup>5</sup> � <sup>10</sup>�10. The measured crack size data are then generated

θ ¼ f g m; ln ð Þ C ; s . In the Paris-Erdogan model, the y-intercept C is very small but changes its magnitude by several orders. Therefore, it would be better to identify logarithm of C. For RUL calculation, the critical crack size is determined as 0.043 m. For the Bayesian method, it is necessary to define the prior distribution and the

likelihood function. In the battery example, it was assumed that noise in data follows a normal distribution. However, when the distribution type of measurement noise is unknown, it is possible that the likelihood function might be different from the true noise distribution. The same is true for the prior/initial distribution. Therefore, it would be a good exercise to study the effect of different distribution types by changing the MATLAB codes. In this example, the lognormal distribution

the true crack size data are generated at every 50 cycles using Eq. (17) with

measured crack size data are used to identify three model parameters,

The system is under fatigue loading with the range of stress being Δσ ¼ 75MPa at each cycle. It is assumed that the health monitoring is performed every 50 cycles to measure the crack size yk until the current time tCUR ¼ 1, 200 cycles, and the initial size of the crack is a<sup>0</sup> ¼ 0:01m. Similar to the battery degradation example, the measurement data are simulated by adding random noise to the true crack size. First,

<sup>~</sup>y tð Þ¼ ; <sup>m</sup>;<sup>C</sup> tC <sup>1</sup> � <sup>m</sup>

4. Application to crack growth prognostics

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

4.1 Model definition: crack growth

parameters, <sup>Δ</sup><sup>K</sup> <sup>¼</sup> <sup>Δ</sup><sup>σ</sup> ffiffiffiffiffi

following degradation model:

by adding Gaussian noise ε � N 0; s

is employed for the likelihood function as:

21

can be modified.

```
figure; % histogram of parameters
for j=1:p;
 subplot(1,p,j);
 hist(samplResul(j,:),30);
end
```
When the true degradation model is known, it is possible to compare the predicted degradation with the true one. The following MATLAB script plots the median and confidence intervals of the predicted degradation along with the true degradation and the threshold.

```
figure;
degraTrue=exp(-0.012.*time);
degraPI=prctile(degraPredi',perceValue)';
plot(time,degraTrue,'k'); hold on;
plot(time(1:k1),measuData,'*b');
plot(time(k1:end),degraPI(:,1),'–r');
plot(time(k1:end),degraPI(:,2:3),':r');
plot([0 time(end)],[thres thres],'g');
xlabel('weeks'); ylabel('Relative C/1 capacity');
```
The degradation curves up to 50 weeks are shown in Figure 8 for the battery example. The true degradation with b ¼ 0:012 is shown with the black curve. The red curves show 5, 50 (median) and 95 percentiles of the predicted degradation, which are caused by signiLevel = 5 (line 14). The plot also shows the threshold (green line) and measurement data (blue asterisk marks). Based on the true model, the end of life of the battery is tEOL ¼ 29:72 weeks, and thus, the true RUL should be tRUL ¼ 20:72 weeks. The prediction shows that the median of RUL is

Figure 8. Comparison of the predicted degradation with the true degradation.

20.38 weeks, which is close to the true RUL. In addition, the 90% confidence interval is about 3.4 weeks; that is, the uncertainty in the prediction is about 17%.
