Degradation Model in MATLAB

Question

I am trying (using MATLAB) to generate the following image (from the section 4.1 Numerical Examples) of Wu Tian Chen research article 'Condition-based Maintenance Optimization Using Neural Network-based Health Condition Prediction':

The model is derived from the equation: $$ L(t)=θ′+β′t+ε(t) $$ where $θ′$ is a normal random variable with mean $5$ and variance $1{^2}$, and $β′$ also be a normal random variable with mean $5$ and variance $1.5{^2}$. $ε(t)=σW(t)$ is a centered Brownian motion such that the mean of $ε(t)$ is zero and the variance of $ε(t)$ is $σ{^2}t$.The parameters are set according to the article.

The MATLAB code I have written is the following:

% Initialization:
Ts          = 0;
Te          = 150;
Tn          = 101;
mu0         = 5;
sigma0      = 1;
mu1         = 5;
sigma1      = 1.5;
sigma       = .5;
D           = 500; % failure threshold
paths = 50; % Number of paths in accordance with the article and graph

figure, hold on, box on, grid off

t = linspace(Ts, Te, Tn)';

for i = 1:paths 
   % Generate θ':
    teta1 = randn() * sigma0 + mu0;
   % Generate β':
    beta1 = (randn() * sigma1 + mu1) .* t;

   % Generate Brownian path
    dW = sqrt(Te / Tn) * randn(Tn, 1);
    W = cumsum(dW, 1); % cumulative sum
    e_t = sigma * W;

   L = teta1 + beta1 + e_t;   
   plot(t, L);
   xlim([Ts, Te])
   ylim([0, 600])

end

% Draws threshold
plot([Ts, Te], [D, D], 'k', 'LineWidth', 2.5)

title('Degradation Signal', 'FontWeight', 'bold', 'FontSize', 14);
xlabel('Time (day)', 'FontSize', 12);
ylabel('Amplitude', 'FontSize', 12);

The output of this code is different from the original one, but I don't understand why. In have already applied the changes from A. Donda, but the plot is still different, due to the abscence of fluctuations.

Thank you

Answer 1

First, your code for the linear part is wrong. As far as I understand, $\theta'$ and $\beta'$ vary randomly over instantiations, but are constant over time. To implement this, you have to write

% Generate θ':
theta1 = randn() * sigma0 + mu0;
% Generate β':
beta1 = (randn() * sigma1 + mu1) .* t;

Secondly, there are several strange things in your code to generate a Brownian path: There's an $\exp()$ without motivation, you use the constant 0.5 instead of the defined parameter sigma, and you have an additional factor 0.5 in the code for dW. IMHO, that section should look like this:

dW = sqrt(Te / Tn) * randn(Tn, 1);
W = cumsum(dW);
e_t = sigma * W;

There'd be also a problem here if Ts is not zero, so I'd recommend to just remove that parameter and replace it by 0 everywhere.

With this modified code, the result looks like this:

which is still not what you aim at. However, I think the parameters you use are not consistent with your plot anyway. For one, if the variance of $\beta'$ is supposed to be 1.5, then you have to write

sigma1 = sqrt(1.5);

since $\sigma_1$ is the standard deviation. Moreover, if $\sigma = 0.5$ for the Brownian motion, then the variance of $\epsilon(150)$ is 37.5, and the standard deviation 6.12. The fluctuations in your plot however are quite a bit larger.

Hope this helps.