%% Routine for Multi-Stress Lognormal Plotting Example
% Coded by Reuel Smith 2015-2017
% v. MATLAB R2015b through 2017a
% ========================================================================
% Example Figure 4.7
% Consider the set of data provided in Table 4.1 that represents the ALT
% data of a component in which temperature is the stress-causing agent of
% failure. Further assume that it has already been shown that the scatter
% in times to failure data listed in Table 4.1 can be best modeled by the
% lognormal distribution, and that the Arrhenius model (exponential) is the
% most appropriate life-stress relationship for this particular ALT.
%
% Temperature # of Tested Units Recorded Failure Times (hrs)
% =========== ================= ======================================
% 150 ?C 3 2350 2560 2980
% 200 ?C 9 220 250 330 370 380 460 460 510 610
%
% Solution:
% This routine will generate the Lognormal probability plots for the two
% stress levels 150 ?C (423 ?K) and 200 ?C (472 ?K)
% ========================================================================
clc
clear
% ========================================================================
% PART 1: Initialize data and constants
% ========================================================================
% sorting the times to failure in ascending order
TTF1 = [2350;2560;2980];
TTF2 = [220;250;330;370;380;460;460;510;610];
% set the length of the TTF vector or the number of failures
n1 = length(TTF1);
n2 = length(TTF2);
% CDF or unreliability for probability plot at Stress Levels 1 and 2
i1 = 1:n1; % uncensored data points from 1 to n1
i2 = 1:n2; % uncensored data points from 1 to n2
% CDF or Failure Percent for probability plot under the Kimball probability
% plotting position: (i - 0.375)/(n + 0.25). Other plotting positions may
% be assumed:
% Midpoint - (i - 0.5)/n
% Mean - i/(n + 1)
% Median - (i - 0.3)/(n + 0.4)
F1 = 100.*((i1'-0.375)/(n1+0.25));
F2 = 100.*((i2'-0.375)/(n2+0.25));
% Lognormal scale conversion of the Failure Percent data
F1B = icdf('Normal',F1./100,0,1);
F2B = icdf('Normal',F2./100,0,1);
% reliability for probability plot
R1 = (100-F1)./100;
R2 = (100-F2)./100;
% Scale conversion to Lognormal plot scale
% Upper and lower bounds of the Percent Failure axis in percent
fc = [.1 99.9];
% Lognormal scale conversion of the upper and lower bounds of the Percent
% Failure axis
fcB = icdf('Normal',fc./100,0,1);
% 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];
% Lognormal scale conversion of Percent Failure ticks cited between the
% lower and upper bounds
Pticks = icdf('Normal',Pticks1./100,0,1);
% ========================================================================
% PART 2: Solve for Lognormal Parameters log mean and log sigma
% ========================================================================
% Use a polynomial fit (polyfit) to draw the best fit line (fc) through the
% known times to failure TTF from F(t) = 1% to F(t) = 99%
% polynomial fit for LN(TTF) and F
pb_1 = polyfit(log(TTF1),F1,1);
pb_2 = polyfit(log(TTF2),F2,1);
% the log mean of the data
meant_1 = (50-pb_1(2))/pb_1(1);
meant_2 = (50-pb_2(2))/pb_2(1);
% the 84% value for LN(t)
t84_1 = (84-pb_1(2))/pb_1(1);
t84_2 = (84-pb_2(2))/pb_2(1);
% the log-sigma
sigmat_1 = t84_1-meant_1;
sigmat_2 = t84_2-meant_2;
disp(['Lognormal parameters at 150 C (423 K): log mean = ',num2str(meant_1),' and log sigma = ',num2str(sigmat_1)])
disp(['Lognormal parameters at 200 C (473 K): log mean = ',num2str(meant_2),' and log sigma = ',num2str(sigmat_2)])
% ========================================================================
% PART 3: Lognormal probability plot
% ========================================================================
% Plot points at 0.1% and 99.9%
% These have a different polynomial pair than that which was used to find
% the mean and SD
pb_1B = polyfit([meant_1 meant_1+sigmat_1],[icdf('Normal',0.5,0,1) icdf('Normal',0.84,0,1)],1);
pb_2B = polyfit([meant_2 meant_2+sigmat_2],[icdf('Normal',0.5,0,1) icdf('Normal',0.84,0,1)],1);
ttfc_1 = exp([(icdf('Normal',0.001,0,1)-pb_1B(2))/pb_1B(1) (icdf('Normal',0.999,0,1)-pb_1B(2))/pb_1B(1)]);
ttfc_2 = exp([(icdf('Normal',0.001,0,1)-pb_2B(2))/pb_2B(1) (icdf('Normal',0.999,0,1)-pb_2B(2))/pb_2B(1)]);
figure(1)
semilogx(TTF1,F1B,'b*',ttfc_1,fcB,'b-',TTF2,F2B,'rs',ttfc_2,fcB,'r-');
hold on
points1 = semilogx(TTF1,F1B,'b*',TTF2,F2B,'rs');
hold off
xlabel('Time (hours)')
ylabel('Percent Failure')
grid on
set(gca,'YTick',Pticks,'YTickLabel',{'0.10';'0.20';' ';' ';'0.50';' ';' ';' ';' ';'1.00';'2.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',{'100';'1000';'10000'}) % CHANGE SCALE MANUALLY OR COMMENT OUT
legend(points1,'Stress Level 1 - 423 deg K','Stress Level 2 - 473 deg K','Location','Northwest')
ylim([icdf('Normal',0.001,0,1) icdf('Normal',0.999,0,1)])
set(gcf,'color','w');
Top