%% Matlab Routine to Perform MLE on a Step-Stress Testing Example
% Coded by Reuel Smith 2015-2017
% v. MATLAB R2015b through 2017a
% ========================================================================
% Example Problem 4.17
% For the data given in Table 4.8,
% ========================================================================
% Step Stress Test Time Test Result (time of failure occurrence)
% (volts) (hrs)
% ==== ======= ========= ==========================================
% 1 2 250 0 fails in 11
% 2 3 100 3 fails in 11 (30; 60; 80)
% 3 4 20 3 fails in 8 (2; 10; 16)
% 4 5 10 3 fails in 5 (1; 4; 8)
% 5 6 10 2 fails in 2 (1; 5)
% 6 7 10 0 fails in 0
% ========================================================================
% stresses are
%
% S_1=2, 0<t<=250
% S_2=3, 250<t<=350
% S_3=4, 350<t<=370
% S_4=5, 370<t<=380
% S_5=6, 380<t<=390
% S_6=7, 390<t<=Inf
%
% Using the log-likelihood function in Equation (4.139) estimate the
% parameters beta, K, and n.
% ========================================================================
% Solution:
% The derivatives of the log-likelihood L, dL/dbeta=0, dL/dK=0, and dL/dn=0
% are found by an online derivative calculator and the values are input
% into the MATLAB routine. The 0 values of all three parameter derivatives
% (or approximation to 0) are found % through MATLAB's FSOLVE command, that
% solves for the zero value of a set of functions in a vector. Online
% derivative calculators include:
% https://www.derivative-calculator.net/
% http://www.wolframalpha.com/calculators/derivative-calculator/
% https://www.symbolab.com/solver/first-derivative-calculator
% ========================================================================
clc
clear
% ========================================================================
% Initialize data and constants
% ========================================================================
% Definitions of different variables used in this code
% x(1): beta parameter
% x(2): K parameter
% x(3): n parameter
% Stress Levels 1 through 6 in Volts
V = [2 3 4 5 6 7];
% The end time of each of the Stress Levels 1 through 5 (in hours)
tend = [250 350 370 380 390];
% Times to failure
% These times have already been defined with respect to Step 1. That is
% Step 2 failure time 30 hrs is defined as 280 hrs with respect to the 250
% hour test time of Step 1.
TTF = [280 310 330 352 360 366 371 374 378 381 385]';
% ========================================================================
% The next step is to define the the failure and survival times as well as
% failure and survival stress levels for each of the six step stress
% levels. Note that all instances where the failure time vector (TTF) and
% survival time vector (TTS) is left as 0, are not added to the operation
% because the log-likelihood only counts the times failed and succeded
% through all the steps.
% ========================================================================
% STEP 1: 2V; 250hr test time; 0 fails in 11
% Failure Times in hours (0)
TTF1 = 0;
% Survival times in hours (11)
TTS1 = tend(1).*ones(11,1);
% Failure Stress (V)
VF1 = V(1);
% Survival Stress (V)
VS1 = V(1).*ones(11,1);
% STEP 2: 3V; 100hr test time; 3 fails in 11 (30; 60; 80)
% Failure Times in hours (3)
TTF2 = TTF(1:3);
% Survival times in hours (8)
TTS2 = tend(2).*ones(8,1);
% Failure Stress (V)
VF2 = V(2).*ones(3,1);
% Survival Stress (V)
VS2 = V(2).*ones(8,1);
% STEP 3: 4V; 20hr test time; 3 fails in 8 (2; 10; 16)
% Failure Times in hours (3)
TTF3 = TTF(4:6);
% Survival times in hours (5)
TTS3 = tend(3).*ones(5,1);
% Failure Stress (V)
VF3 = V(3).*ones(3,1);
% Survival Stress (V)
VS3 = V(3).*ones(5,1);
% STEP 4: 5V; 10hr test time; 3 fails in 5 (1; 4; 8)
% Failure Times in hours (3)
TTF4 = TTF(7:9);
% Survival times in hours (2)
TTS4 = tend(4).*ones(2,1);
% Failure Stress (V)
VF4 = V(4).*ones(3,1);
% Survival Stress (V)
VS4 = V(4).*ones(2,1);
% STEP 5: 6V; 10hr test time; 2 fails in 2 (1; 5)
% Failure Times in hours (2)
TTF5 = TTF(10:end);
% Survival times in hours (8)
TTS5 = 0;
% Failure Stress (V)
VF5 = V(5).*ones(2,1);
% Survival Stress (V)
VS5 = V(5);
% STEP 6: 7V; 10hr test time; 0 fails in 0
% Failure Times in hours (0)
TTF6 = 0;
% Survival times in hours (0)
TTS6 = 0;
% Failure Stress (V)
VF6 = V(6);
% Survival Stress (V)
VS6 = V(6);
% ========================================================================
% SETUP FOR TAU (tau1 through tau5 for stress levels 2 through 6
% This establishes the equivalent time (tau_(i-1)) of step i-1 if the item
% is operated at step stress level i. Tau_1 follows Equation 4.135 in the
% textbook while all other tau values follow Equation 4.136.
%
% 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): beta parameter, x(2): K parameter, and x(3): n parameter
% ========================================================================
% 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
% ========================================================================
% tau for step 2
tau1 = @(x) (tend(1)).*(V(1)./V(2)).^x(3);
% tau for step 3
tau2 = @(x) (tend(2)-tend(1)+tau1(x)).*(V(2)/V(3)).^x(3);
% tau for step 4
tau3 = @(x) (tend(3)-tend(2)+tau2(x)).*(V(3)/V(4)).^x(3);
% tau for step 5
tau4 = @(x) (tend(4)-tend(3)+tau3(x)).*(V(4)/V(5)).^x(3);
% tau for step 6
tau5 = @(x) (tend(5)-tend(4)+tau4(x)).*(V(5)/V(6)).^x(3);
% ========================================================================
% SETUP FOR ADJUSTED TIME VECTORS
% The procedure for setting up the adjusted failure time vector ti, and the
% adjusted survival time vector tj follows the expression given in Equation
% 4.134. The vector tvect is the t_(i-1) value, the end time of step i-1,
% and tau is the tau_(i-1) value. Each failure time and survival time is
% assigned a t_(i-1) and tau_(i-1) according to the step.
tvecti = [tend(1).*ones(3,1);tend(2).*ones(3,1);tend(3).*ones(3,1);tend(4).*ones(2,1)];
tvectj = [zeros(11,1);tend(1).*ones(8,1);tend(2).*ones(5,1);tend(3).*ones(2,1)];
taui = @(x) [tau1(x);tau1(x);tau1(x);tau2(x);tau2(x);tau2(x);tau3(x);tau3(x);tau3(x);tau4(x);tau4(x)];
tauj = @(x) [zeros(11,1);tau1(x);tau1(x);tau1(x);tau1(x);tau1(x);tau1(x);tau1(x);tau1(x);tau2(x);tau2(x);tau2(x);tau2(x);tau2(x);tau3(x);tau3(x)];
ti = @(x) [TTF2;TTF3;TTF4;TTF5] - tvecti + taui(x);
tj = @(x) [TTS1;TTS2;TTS3;TTS4] - tvectj + tauj(x);
% The failure and survival stress levels are assigned a vector of the same
% sizes as the adjusted time vectors for the MLE operation.
Vi = [VF2;VF3;VF4;VF5];
Vj = [VS1;VS2;VS3;VS4];
% ========================================================================
% SETUP FOR DERIVATIVES
% The derivatives defined here were computed in an online derivative
% calculator. When converting your own equations, take care to group the
% failure and survival stresses and times together for summation of each.
% Derivative dL/dbeta
deriv1 = @(x) -((sum((x(2).*Vj.^x(3).*tj(x)).^x(1).*log(x(2).*Vj.^x(3).*tj(x)))+sum((x(2).*Vi.^x(3).*ti(x)).^x(1).*log(x(2).*Vi.^x(3).*ti(x))-log(x(2).*Vi.^x(3).*ti(x)))).*x(1)-1)./x(1);
% Derivative dL/dK
deriv2 = @(x) -(x(1).*(sum((Vj.^x(3).*tj(x).*x(2)).^x(1))+sum((Vi.^x(3).*ti(x).*x(2)).^x(1))-1))./x(2);
% Derivative dL/dn
deriv3 = @(x) -x(1).*(sum(log(Vj).*(x(2).*Vj.^x(3).*tj(x)).^x(1))+sum(log(Vi).*(x(2).*Vi.^x(3).*ti(x)).^x(1)-log(Vi)));
% Derivative vector [dL/dbeta;dL/dK;dL/dn]
F1 = @(x) [deriv1(x);deriv2(x);deriv3(x)];
% Initial estimate for parameters
% Take care to find an initial estimate that coincides with the assigned
% parameter bounds. For example, beta can only be positive so select a
% positive initial guess. Note however that there may be multiple guesses
% that produce a solution that may or may not be close to zero. A
% reasonable first guess may be acquired by using the Plotting Step-Stress
% method defined in section 4.16.1. Randomly estimating x0 will almost
% always result in a poor estimate.
x0 = [5 1e-6 4];
% x0 = [5 1 4];
% PARAMHAT solve
% This is the initial calculation for the MLE parameter estimates using
% FSOLVE on the derivative vector. Sometimes the warning message "fsolve
% stopped because it exceeded the function evaluation limit,
% options.MaxFunEvals = 300 (the default value)." will appear because
% there may not be an exact solution where the derivative vector is
% completely zero. In those cases, run the FSOLVE routine as a loop until
% the parameter estimate PARAMHAT converges to a solution.
% initialize PARAMHAT estimates
paramhatv = x0;
for i = 1:100
if i == 1
paramhat = fsolve(F1,x0);
else
paramhat = fsolve(F1,paramhat);
end
paramhatv(i+1,:) = paramhat;
end
% Display the derivative vector based on the PARAMHAT estimate
F1(paramhat)
disp(['MLE Estimates: beta = ',num2str(paramhat(1)),', K = ',num2str(paramhat(2)), ', and n = ', num2str(paramhat(3))])
% ========================================================================
% Plot the iterations of the MLE parameter estimates to determine
% convergance
figure(1)
subplot(3,1,1)
plot(paramhatv(:,1),'b+-')
xlabel('iteration')
ylabel('\beta parameter')
subplot(3,1,2)
plot(paramhatv(:,2),'rs-')
xlabel('iteration')
ylabel('K parameter')
subplot(3,1,3)
plot(paramhatv(:,3),'g^-')
xlabel('iteration')
ylabel('n parameter')
set(gcf,'color','w');
Top