%% Routine for Multi-Stress Weibull Plotting Example
% Coded by Reuel Smith 2015-2017
% v. MATLAB R2015b through 2017a
% ========================================================================
% Example Problem 4.7
% Consider an accelerated test where temperature is the accelerating
% variable and 20 units were tested to failure assuming complete failures
% (no censoring). Eight units were tested at 406 ?K, and six units each at
% 436 ?K and 466 ?K, with times to failure tabulated below. Assuming an
% Arrhenius-Weibull life-stress relationship, find the parameters of the
% model, using both the plotting method and MLE method and compare the
% results. Find the expected life at the use level temperature of 353 ?K
% and compare the results.
%
% Temperature Time to Failure (hrs)
% =========== =====================================
% 406 K 248 456 528 731 813 965 972 1528
% 436 K 164 176 289 319 386 459
% 466 K 92 105 155 184 219 235
%
% Solution:
% This routine will only cover the Weibull multi-plot at the three
% different stress levels.
% ========================================================================
clc
clear
% ========================================================================
% PART 1: Initialize data and constants
% ========================================================================
% sorting the times to failure in ascending order (CHANGE NUMBERS AND
% AMOUNT FOR DIFFERENT DATA)
TTF1 = [248 456 528 731 813 965 972 1528]';
TTF2 = [164 176 289 319 386 459]';
TTF3 = [92 105 155 184 219 235]';
% set the length of the TTF vector or the number of failures
n1 = length(TTF1);
n2 = length(TTF2);
n3 = length(TTF3);
% CDF or unreliability for probability plot at Stress Levels 1, 2, and 3
i1 = 1:n1; % uncensored data points from 1 to n1
i2 = 1:n2; % uncensored data points from 1 to n2
i3 = 1:n3; % uncensored data points from 1 to n3
% Median plotting position (i - 0.3)/(n + 0.4). Other plotting positions
% may be assumed:
% Kimball - (i - 0.375)/(n + 0.25)
% Midpoint - (i - 0.5)/n
% Mean - i/(n + 1)
F1 = 100.*((i1'-0.3)./(n1+0.4));
F2 = 100.*((i2'-0.3)./(n2+0.4));
F3 = 100.*((i3'-0.3)./(n3+0.4));
% Weibull scale conversion of the Failure Percent data
F1B = log(log(1./(1-F1./100)));
F2B = log(log(1./(1-F2./100)));
F3B = log(log(1./(1-F3./100)));
% reliability for probability plot
R1 = (100-F1)./100;
R2 = (100-F2)./100;
R3 = (100-F3)./100;
% Scale conversion to Weibull plot scale
% Upper and lower bounds of the Percent Failure axis in percent
fc = [.1 99.9];
% Weibull scale conversion of the upper and lower bounds of the Percent
% Failure axis
fcB = log(log(1./(1-fc./100)));
% Percent Failure ticks cited between the lower and upper bounds
Pticks1 = [0.1:0.1:1 2:1:10 20:10:90 95 99 99.9];
% Weibull scale conversion of Percent Failure ticks cited between the lower
% and upper bounds
Pticks = log(log(1./(1-Pticks1./100)));
% ========================================================================
% PART 2: Solve for Weibull Parameters alpha and beta
% ========================================================================
% The Weibull parameters alpha and beta may be computed by assuming a
% relationship for the reliability function R(t) = exp[-(t/alpha)^beta] or
% 1/R(t) = exp[(t/alpha)^beta]. The relation translates to:
%
% ln[ln(1/R(t))] = beta*ln(t) - beta*ln(alpha)
%
% The POLYFIT function is used to fit this linear relationship and extract
% the Weibull parameters alpha and beta. See Reliability Engineering and
% Risk Analysis: A Practical Guide Section 3.3.2.2. for method.
wblparm_1 = polyfit(log(TTF1),log(log(1./R1)),1);
wblparm_2 = polyfit(log(TTF2),log(log(1./R2)),1);
wblparm_3 = polyfit(log(TTF3),log(log(1./R3)),1);
beta_1 = wblparm_1(1);alpha_1 = exp(-wblparm_1(2)/beta_1);
beta_2 = wblparm_2(1);alpha_2 = exp(-wblparm_2(2)/beta_2);
beta_3 = wblparm_3(1);alpha_3 = exp(-wblparm_3(2)/beta_3);
% The average beta value
betamean = mean([beta_1 beta_2 beta_3]);
disp(['Weibull parameters at 406 K: alpha = ',num2str(alpha_1),' hours and beta = ',num2str(beta_1)])
disp(['Weibull parameters at 436 K: alpha = ',num2str(alpha_2),' hours and beta = ',num2str(beta_2)])
disp(['Weibull parameters at 466 K: alpha = ',num2str(alpha_3),' hours and beta = ',num2str(beta_3)])
disp(['Average beta = ',num2str(betamean)])
% ========================================================================
% PART 3: Weibull probability plot
% ========================================================================
% Plot points at 0.1% and 99.9%
% These have a different polynomial pair than that which was used to find
% the alpha and beta
% ========================================================================
intercept_1 = wblparm_1(2);
intercept_2 = wblparm_2(2);
intercept_3 = wblparm_3(2);
ttfc_1 = exp([(log(log(1./(1-0.001)))-intercept_1)/beta_1 (log(log(1./(1-0.999)))-intercept_1)/beta_1]);
ttfc_2 = exp([(log(log(1./(1-0.001)))-intercept_2)/beta_2 (log(log(1./(1-0.999)))-intercept_2)/beta_2]);
ttfc_3 = exp([(log(log(1./(1-0.001)))-intercept_3)/beta_3 (log(log(1./(1-0.999)))-intercept_3)/beta_3]);
figure(1)
semilogx(TTF1,F1B,'b*',TTF2,F2B,'rs',TTF3,F3B,'gp',ttfc_1,fcB,'b-',ttfc_2,fcB,'r-',ttfc_3,fcB,'g-');
hold on
points1 = semilogx(TTF1,F1B,'b*',TTF2,F2B,'rs',TTF3,F3B,'gp');
hold off
xlabel('Time (hours)')
ylabel('Percent Failure')
grid on
set(gca,'YTick',Pticks,'YTickLabel',{'0.10';'0.20';'0.30';' ';'0.50';' ';' ';' ';' ';'1.00';'2.00';'3.00';' ';'5.00';' ';' ';' ';' ';'10.00';'20.00';'30.00';'40.00';'50.00';'60.00';'70.00';'80.00';'90.00';'95.00';'99.00';'99.90'})
set(gca,'XTickLabel',{'10';'100';'1000';'10000'}) % CHANGE SCALE MANUALLY OR COMMENT OUT
legend(points1,'Stress Level 1 - 406 deg K','Stress Level 2 - 436 deg K','Stress Level 3 - 466 deg K','Location','Northwest')
ylim([log(log(1./(1-0.001))) log(log(1./(1-0.999)))])
Top