%% Matlab Routine to Perform a Step-Stress Testing Example
% Coded by Reuel Smith 2015-2017
% v. MATLAB R2015b through 2017a
% ========================================================================
% Example Problem 4.19
% Consider a step-stress test of cable insulation. The test was run to
% estimate life at a design stress of 400 volts/mm. The data in Table 4.24
% illustrate the data analysis methods. Each specimen was first held for
% 10 minutes each at 5 kV, 10 kV, 15 kV, and 20 kV before it went into step
% 1 at 26 kV. In steps 1 through 10, a specimen has the same hold time at
% each voltage (15 minutes, 1 hour, 4 hours, or 16 hours). Table 4.24
% shows the step number and the total time on test when a specimen broke
% down, and its insulation thickness (used to calculate its stress as the
% voltage divided by the thickness). Assume that this tested cable 1 is
% comparable to another cable 2 tested before, so we can use the MLE result
% of the test of cable 2 as the prior information. Assume that cable 2 was
% tested and the data shown in Table 4.25 were obtained. Update the models
% parameters using the Bayesian estimation.
%
% Solution:
% Follow the methodology outlined in Section 4.16.3 with the second
% likelihood function where all the previous steps are converted to the
% final step n.
% ========================================================================
clear
clc
% ========================================================================
% Initialize data and constants
% ========================================================================
% Definitions of different variables used in this code
% x(1): A
% x(2): p
% x(3): beta
% Voltage pattern with step number on the left column and voltage in
% kilovolts on the right
voltpattern = [1 26;2 28.5;3 31;4 33.4;5 36;6 38.5;7 41;8 43.5;9 46;10 48.5];
V = 1000.*voltpattern(:,2); % voltage in volts
% Total number of steps
STEPS = 10;
% Cable data:
% First Column - The hold time for each step. If the final step is 5, then
% the total time elapsed will be Hold Time x 5
% Second Column - Final Step which is the step that testing ends either
% because of specimen failure or for censoring
% Third Column - Total Time-to-Failure. Times-to-failure that are given as
% negative are treated as censored data or a runtime without failure.
% Fourth Column - Thickness of the cable insulation (mm).
% Cable data for first one (will be for updating of priors)
cable1 = [15 5 102 27;...
15 5 113 27;...
15 5 113 27;...
60 6 -370 29.5;...
60 6 345 29.5;...
60 6 -345 28;...
240 6 1249 29;...
240 6 1333 29;...
240 6 -1333 29;...
240 5 1096.9 29;...
240 6 1250.8 30;...
240 5 1097.9 29;...
960 3 2460.9 30;...
960 3 -2460.9 30;...
960 3 2700.4 30;...
960 4 2923.9 30;...
960 2 1160 30;...
960 3 1962.9 30;...
960 1 -363.9 30;...
960 1 -898.4 30;...
960 5 4142.1 30];
% Cable data for second one (will be used to obtain the priors
cable2 = [15 5 114 27;...
15 6 119 27;...
15 6 116 27;...
15 6 117 27;...
60 7 423 29.5;...
60 7 -423 28;...
240 6 1401.5 29;...
240 6 1261.7 29;...
240 6 1419.4 29;...
240 6 1364.2 29;...
960 2 1398.1 30;...
960 3 2691.7 30;...
960 3 2000 30;...
960 3 -2000 30;...
960 3 2149.6 30;...
960 2 1429.7 30;...
960 2 1502.7 30;...
960 2 1440.9 30;...
960 4 3349.1 30;...
960 4 -3349.1 30];
% The initial time in minutes. Before step 1, each specimen was held 10
% minutes each at 5, 10, 15, and 20 kV, so 40 minutes passed before the
% start of the step-stress testing.
t0 = 40;
% ========================================================================
% Separation of Failure and Survivor Data for Cable 1
% ========================================================================
% This loop takes the existing cable data and computes the stress levels of
% each data (in V/mm) and the true start time based on the initial testing
% time of 40 minutes. The true start time is computed from based on the
% final step n t_(n-1)
for i=1:length(cable1(:,1))
% Stress levels of Cable 1 in V/mm (i)
Si_1(i,1)=(V(cable1(i,2)))/cable1(i,4);
% Final Step Stress levels of Cable 1 in V/mm
SF_1(i,1)=V(STEPS)/cable1(i,4);
% Previous time on Cable 1 in minutes (i-1)
if cable1(i,2)>1
tim1_1(i,1)=t0+(cable1(i,2)-1)*(cable1(i,1));
else
tim1_1(i,1)=t0;
end
timF_1(i,1)=t0+(STEPS-1)*(cable1(i,1));
end
cable1 = [cable1 tim1_1 timF_1 Si_1 SF_1];
% This loop separates the failure data from the survivor data for
% preparation of the Bayesian analysis step
i_F=1;i_S=1;
for i=1:length(cable1(:,1))
if cable1(i,3)>0
cable1F(i_F,:) = cable1(i,:);
i_F = i_F + 1;
else
cable1S(i_S,:) = abs(cable1(i,:));
i_S = i_S + 1;
end
end
% ========================================================================
% Separation of Failure and Survivor Data for Cable 2
% ========================================================================
% This loop takes the existing cable data and computes the stress levels of
% each data (in V/mm) and the true start time based on the initial testing
% time of 40 minutes.
for i=1:length(cable2(:,1))
% Stress levels of Cable 2 in V/mm
Si_2(i,1)=(V(cable2(i,2)))/cable2(i,4);
% Final Step Stress levels of Cable 2 in V/mm
SF_2(i,1)=V(STEPS)/cable2(i,4);
% Previous time on Cable 2 in minutes (i-1)
if cable2(i,2)>1
tim1_2(i,1)=t0+(cable2(i,2)-1)*(cable2(i,1));
else
tim1_2(i,1)=t0;
end
timF_2(i,1)=t0+(STEPS-1)*(cable2(i,1));
end
cable2 = [cable2 tim1_2 timF_2 Si_2 SF_2];
% This loop separates the failure data from the survivor data for
% preparation of the Bayesian analysis step
i_F=1;i_S=1;
for i=1:length(cable2(:,1))
if cable2(i,3)>0
cable2F(i_F,:) = cable2(i,:);
i_F = i_F + 1;
else
cable2S(i_S,:) = abs(cable2(i,:));
i_S = i_S + 1;
end
end
% ========================================================================
% Initialize Formulae for evaluation
% ========================================================================
% SETUP FOR TAU (tau1 through tau9 for stress levels 2 through 10. The
% tau values calculated follow the equation:
%
% tau_(i-1)=(t_(i-1)-t_(i-2)+tau_(i-2) ) (V_(i-1)/V_i )^p
%
% The (V(i-1)/V(i))^p term is defined as the inverse of the acceleration
% factor (AF) between stress V(i-1) and the next stress V(i) under the
% Inverse Power Relation (IPR) for life (or 1/AF = 1/[(V(i)/V(i-1))^p]).
% This value will have to be changed in the tau statements when and if
% assuming a different life model. Parameters being sought must all be
% given an x(*) value:
% Example: IPR - x(1)=A, x(2)=p, x(3)=beta
% ========================================================================
% Acceleration factors: Change (V(i-1)/V(i))^p statement in tau for
% different models:
%
% Inverse Power Law: (V(i-1)/V(i))^p
% Eqn(4.69), V-stress, p-parameter
%
% Arrhenius: exp(a*(1/T(i) - 1/T(i-1)))
% Eqn(4.40), T-temperature, a-slope
%
% Eyring: (V(i-1)/V(i))*exp(a*(1/T(i) - 1/T(i-1)))
% Eqn(4.55), V-stress, T-temperature, a-slope
%
% Dual Stress: exp(a*(1/T(i) - 1/T(i-1)) + b*(1/V(i) - 1/V(i-1)))
% Eqn(4.85), T-temperature, V-non-thermal stress, a-slope 1, b-slope 2
%
% Power Exponential: ((V(i-1)/V(i))^n)*exp(a*(1/T(i) - 1/T(i-1)))
% Eqn(4.97), V-stress, T-temperature, a-slope 1, n-slope 2
% ========================================================================
% Adjusted time for failures
adjtimeC1 = @(x) (cable1F(:,3) - cable1F(:,5)).*((cable1F(:,8)./cable1F(:,7)).^x(2));
adjtimeC2 = @(x) (cable2F(:,3) - cable2F(:,5)).*((cable2F(:,8)./cable2F(:,7)).^x(2));
% Adjusted time for survivals
adjtimeC1S = @(x) (cable1S(:,3) - cable1S(:,5)).*((cable1S(:,8)./cable1S(:,7)).^x(2));
adjtimeC2S = @(x) (cable2S(:,3) - cable2S(:,5)).*((cable2S(:,8)./cable2S(:,7)).^x(2));
% ====================================================================
% SETUP FOR MH SAMPLING
% ====================================================================
% number of randomly generated samples
nsamples=20000;
%this value can remain as it is even if model parameters or model change
K = 1000;
% This parameter controls the size of the new markov chain which omitts M-1
% out of M values. This will curb the effect of autocorrelation
M = 200;
% number of parameters in likelihood equation: x(1): A, p, and beta
n=3;
% Initialize parameter guess
init_param=[0.000000001 1 1];
% sigma for the proposed PDF, an nxn identity matrix
prop_sig=eye(n);
% define proposal distribution (proppdf) and random number generator
% (proprnd). The standard deviation for each parameter in the proprnd is
% to be chosen such that the updated guess covers a reasonable range of
% estimates. These should be small if the initial estimate is a close one.
proppdf = @(x,y) mvnpdf(x,y,prop_sig);
proprnd = @(y) [normrnd(y(1),.000001),normrnd(y(2),.1),normrnd(y(3),.1)];
% Definition of the Likelihood
% In this example, the log-likelihood of Eqn(4.145) is used for the
% Bayesian updating. Often, the product of an initial likelihood estimate
% will be zero, which MHSAMPLE cannot process. In those instances it is
% recommended to process the log-likelihoodxpriors. This is made up of the
% log of Eqns(4.142) which is the Weibull PDF and reliability function of
% an IPR in place of the Weibull alpha parameter. For using other
% Life-stress models and/or distributions, adjust the equation manually and
% enter them into the next three lines of code.
pdf1C1 = @(x,Si) log(x(1).*x(3))+x(2).*Si+(x(3)-1).*log(adjtimeC1(x).*x(1).*Si.^x(3))-(adjtimeC1(x).*x(1).*Si.^x(2));
pdf1C2 = @(x,Si) log(x(1).*x(3))+x(2).*Si+(x(3)-1).*log(adjtimeC2(x).*x(1).*Si.^x(3))-(adjtimeC2(x).*x(1).*Si.^x(2));
R1 = @(x,Si,tr) -((tr.*x(1).*Si.^x(2)).^x(3));
% define the Log-likelihoodxPriors
% The prior distributions for A, p, and beta are assumed to be
% non-informative, therefore they are treated as uniform distributions.
% Parameters A and beta are known to be positive so the uniform
% distribution for them are between 0 and 100. Parameter p is variable
% between -100 and 100 since it can be either positive or negative. These
% priors multiplied by the sum of the log-likelihood to form the
% Log-likelihoodxPriors term that undergoes Bayesian updating.
pdfC2 = @(x) unifpdf(x(1),0,1000).*unifpdf(x(2),-1000,1000).*unifpdf(x(3),0,1000).*(sum(pdf1C2(x,cable2F(:,8)))+sum(R1(x,cable2S(:,8),adjtimeC2S(x))));
% MH-sampling routine
[result,accept] = mhsample(init_param,nsamples,'pdf',pdfC2,'proppdf',proppdf,'proprnd',proprnd,'burnin',K,'thin',M);
% Alternatively the BURN-IN, and THIN value may be left out for a faster,
% but less precise MH-sampling procedure
% [result,accept] = mhsample(init_param,nsamples,'pdf',pdfC2,'proppdf',proppdf,'proprnd',proprnd);
% Fit the posterior parameters A, p, and beta to a uniform parameter. This
% can be fit to other distributions based on the histogram, but for this
% example, it is assumed uniform.
[A_min,A_max] = unifit(result(:,1));
[p_min,p_max] = unifit(result(:,2));
[beta_min,beta_max] = unifit(result(:,3));
prior1 = lognfit(result(:,1));
[prior2m,prior2s] = normfit(result(:,2));
prior3 = lognfit(result(:,3));
% Documents the upper and lower bounds of the three parameters
priortable = [A_min,A_max;p_min,p_max;beta_min,beta_max];
% ====================================================================
% POSTERIOR PARAMETER UPDATE
% ====================================================================
SD1 = std(result(:,1));
SD2 = std(result(:,2));
SD3 = std(result(:,3));
% redefine the Log-likelihoodxPriors
% This time, the prior distributions for A, p, and beta are based on the
% posterior parameter estimates from the Bayesian estimation of Cable 2.
% These priors multiplied by the sum of the log-likelihood to form the
% Log-likelihoodxPriors term that undergoes Bayesian updating. The new
% random number generator (proprnd2) is also based on the new priors.
% pdfC1 = @(x) unifpdf(x(1),A_min,A_max).*unifpdf(x(2),p_min,p_max).*unifpdf(x(3),beta_min,beta_max).*(sum(pdf1C1(x,cable1F(:,8)))+sum(R1(x,cable1S(:,8),adjtimeC1S(x))));
pdfC1 = @(x) lognpdf(x(1),prior1(1),prior1(2)).*normpdf(x(2),prior2m,prior2s).*lognpdf(x(3),prior3(1),prior3(2)).*(sum(pdf1C1(x,cable1F(:,8)))+sum(R1(x,cable1S(:,8),adjtimeC1S(x))));
proprnd2 = @(y) [normrnd(y(1),SD1),normrnd(y(2),SD2),normrnd(y(3),SD3)];
init_param2 = result(1,:);
% MH-sampling routine
[result2,accept2] = mhsample(init_param2,nsamples,'pdf',pdfC1,'proppdf',proppdf,'proprnd',proprnd,'burnin',K,'thin',M);
% Alternatively the BURN-IN, and THIN value may be left out for a faster,
% but less precise MH-sampling procedure
% [result2,accept2] = mhsample(init_param2,nsamples,'pdf',pdfC1,'proppdf',proppdf,'proprnd',proprnd2);
% Fit the new posterior parameters A, p, and beta to a uniform parameter.
% This can be fit to other distributions based on the histogram, but for
% this example, it is assumed uniform.
[A_min2,A_max2] = unifit(result2(:,1));
[p_min2,p_max2] = unifit(result2(:,2));
[beta_min2,beta_max2] = unifit(result2(:,3));
posteriortable = [A_min2,A_max2;p_min2,p_max2;beta_min2,beta_max2];
figure(1)
subplot(1,3,1)
hist([result(:,1) result2(:,1)],100)
xlabel('Parameter A')
subplot(1,3,2)
hist([result(:,2) result2(:,2)],100)
xlabel('Parameter p')
subplot(1,3,3)
hist([result(:,3) result2(:,3)],100)
xlabel('Parameter \beta')
legend('Posterior - Cable 2','Posterior - Cable 1')
set(gcf,'color','w');
Top