%% Weibull Plotting (Single) Example
% Coded by Reuel Smith 2015-2017
% v. MATLAB R2015b through 2017a
% ========================================================================
% Example Problem 5.1
% Five specimens of a new corrosion-resistant material are tested for 240
% hours in a highly corrosive environment. The density of the material is
% 7.6 g/cm3, and the exposed surface area of each specimen is 4.3 cm2. At
% the end of the test period, the measured weight losses in mg were 11.1,
% 10.4, 12.1, 11.4, and 9.8 (assume no measurement error). If a
% degradation of 1 mm or more results in a structural failure, predict the
% failure times for the five specimens. Based on these failure times,
% determine the probability distribution that best represents the life of
% the material.
%
% Solution:
% TTF = 70659, 75415, 64820, 68800, and 80033
% The proposed probability solutions for this problem are Weibull and
% Lognormal. This routine will process the Weibull solution for the
% parameters alpha and beta, and the Weibull fit probability plot
% ========================================================================
clc
clear
% ========================================================================
% PART 1: Initialize data and constants
% ========================================================================
% sorting the times to failure in ascending order (CHANGE NUMBERS FOR
% DIFERENT DATA)
TTF = sort([70659 75415 64820 68800 80033])';
% set the length of the TTF vector or the number of failures
n = length(TTF);
% 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)
i = 1:n; % uncensored data points from 1 to n
F = 100.*((i'-0.375)./(n+0.25));
% Weibull scale conversion of the Failure Percent data
FB = log(log(1./(1-F./100)));
% reliability for probability plot
R = (100-F)./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 = polyfit(log(TTF),log(log(1./R)),1);
beta = wblparm(1);alpha = exp(-wblparm(2)/beta);
disp(['Weibull parameters: alpha = ',num2str(alpha),' hours and beta = ',num2str(beta)])
% ========================================================================
% 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 = wblparm(2);
ttfc = exp([(log(log(1./(1-0.001)))-intercept)/beta (log(log(1./(1-0.999)))-intercept)/beta]);
% Computes the upper and lower bound for the TTF axis in terms of log-time
signs1 = [floor(log10(ttfc(1))):1:ceil(log10(ttfc(2)))];
figure(1)
semilogx(TTF,FB,'b*',ttfc,fcB,'b-');
xlabel('Time (hours)')
ylabel('Percent Failure')
grid on
legend('Data','Weibull best fit line','Location','Northwest')
set(gca,'XTickLabel',{num2str(10.^signs1')})
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'})
xlim([10^signs1(1) 10^signs1(end)])
ylim([log(log(1./(1-0.001))) log(log(1./(1-0.999)))])
set(gcf,'color','w');
Top