Reliability Growth: Method

Reliability Growth: Duane Plot, Crow-AMSAA Model, Goodness of Fit, Bounds and Projection

What reliability growth tracking does

A development test of a new system runs under a test-analyze-and-fix (TAAF) program: failures are logged with the cumulative test time at which they occur, their causes are found, design changes are made, and testing continues on the changed configuration. If the fixes work, the time between failures lengthens as the test proceeds. Reliability growth tracking fits a model to that lengthening to answer three management questions (MIL-HDBK-189C section 6.2): what is the reliability of the configuration now on test (the demonstrated MTBF), is it growing at the planned rate, and when will it reach the requirement if the program continues as it is. The data are the N successive cumulative test times T₁ < T₂ < … < T_N at which failures occurred, summed over every system on test, and the total accumulated time T. Nothing is purged: early failures that were fixed are as much a part of the growth record as late ones, and leaving them out destroys the trend the model measures.

Duane's postulate

Duane (1964) examined the development records of several General Electric systems and found that the cumulative failure rate N(T)/T, plotted against the cumulative test time T on log-log paper, lay on a straight line with a negative slope. Equivalently the cumulative MTBF m_c = T/N(T) rises along a straight line of slope α, the growth rate: ln m_c = ln b + α ln T, where b is the cumulative MTBF at T = 1. The tool fits that line by least squares to the points (T_i, T_i/i), which is the ReliaSoft RGA procedure, and reports α, b and the coefficient of determination. The expected number of failures is then N(T) = T^1−α/b, and its derivative, the instantaneous failure intensity (1 − α)/(b T^α), shows that the current MTBF is the cumulative value divided by (1 − α): the instantaneous MTBF line lies parallel to the cumulative line, shifted up by the factor 1/(1 − α). The regression supplies standard errors: σ² = SSE/(n − 2), SE(α) = σ/√S_xx and SE(ln b) = σ √(Σ(ln T_i)²/(n S_xx)), with t-based bounds on α and b and normal bounds on the fitted line. Duane's slopes ranged from 0.24 to 0.65; MIL-HDBK-189C Table II gives a historical mean of 0.34 (range 0.23 to 0.53) for time- or distance-based systems. The Duane method is deterministic: it makes no allowance for the randomness of the failure times, so it cannot say whether an observed slope differs from zero by more than chance, and its lines pass through the origin on linear paper (zero MTBF at the start of test). Those are the gaps the AMSAA model fills.

The Crow-AMSAA model

Crow (1974) observed that Duane's power law is the mean value function of a non-homogeneous Poisson process (NHPP) with the Weibull-form intensity ρ(t) = λ β t^β−1, so that the expected number of failures by time t is E[N(t)] = λ t^β and the number of failures in any interval is Poisson with mean λ(t₂^β − t₁^β). The shape β governs the trend: β < 1 is reliability growth (falling intensity), β = 1 a homogeneous Poisson process (constant rate, the exponential case, no growth) and β > 1 deterioration; the growth rate is α = 1 − β. The scale λ depends on the unit of time. The model is a Weibull process, not a Weibull distribution: it describes the arrival of failures of one repairable system being improved, and the statistics of the Weibull life distribution do not apply to it. Within one test phase, where fixes are incorporated as the test proceeds, the model describes the growth; across phases with delayed fixes the planning models of the handbook are the right tool, and a jump in the intensity between phases shows as a kink on the log-log plot and a failed goodness-of-fit test.

Maximum likelihood for individual failure times

For N failures at cumulative times T₁ … T_N in a test ending at T* the log-likelihood is N ln λ + N ln β − λ T*^β + (β − 1) Σ ln T_i, and its maximum (MIL-HDBK-189C 6.2.2.8.1) is β̂ = N/(N ln T* − Σ ln T_i) and λ̂ = N/T*^β̂, which equates the expected failures at T* to the observed N. For a time-terminated test T* is the planned end time T (the test may also simply be in progress with data to date); for a failure-terminated test T* is the last failure time T_N, and the same formula applies with the last term cancelling. The intensity at the end of test is ρ(T) = λ β T^β−1 = β (N/T): the model takes the average rate N/T and reduces it by the factor β, as if the fraction 1 − β of the failures had been designed out, which is the sense in which the estimator is "self-purging" although every failure is used. Its reciprocal m(T) = 1/ρ(T) is the demonstrated MTBF, the MTBF the configuration on test would show if no further changes were made; the cumulative MTBF T/N is β times it. The estimate β̂ is biased high in small samples; the unbiased estimate is β̄ = (N − 1)β̂/N for a time-terminated test and (N − 2)β̂/(N − 1) for a failure-terminated test. The handbook uses β̄ in the goodness-of-fit statistic and β̂ elsewhere, and so does the tool (ReliaSoft reports β̂ by default).

Grouped data

When failures are found at inspections and only the number in each interval is known (at least three intervals, failures in at least two), the likelihood is the product of the Poisson probabilities of the interval counts, and β̂ is the root of Σ F_i [(t_i^β ln t_i − t_i−1^β ln t_i−1)/(t_i^β − t_i−1^β) − ln t_K] = 0 with t₀ ln t₀ = 0 (MIL-HDBK-189C 6.2.3.1), solved by bisection; λ̂ = N/t_K^β̂ as before. The intensity and MTBF functions follow, but because growth under grouped data happens configuration by configuration the handbook also reports the MTBF of the last group, M_K = (t_K − t_K−1)/E_K, the last interval's length over its expected failures. The fit is tested by chi-square: the expected failures E_i = λ(t_i^β − t_i−1^β) are pooled over adjacent intervals until each is at least five, and Σ(O − E)²/E is compared with the chi-square distribution on K_R − 2 degrees of freedom.

Several systems

When several prototypes are tested concurrently the model needs the cumulative test time on all systems at each failure. If the systems run on their own clocks, each starting at zero, the tool builds the equivalent single system (ESS, ReliaSoft RGA): a failure of one system at its clock time t is placed at Σ_j min(t, T_j), the operating time of every system at that instant (a system that has already stopped contributes its end time), and the ESS ends at the sum of the end times. The handbook's two-prototype example is already in that form: its "cumulative hours" column is the sum of the two prototypes' clocks at each failure. If the systems' operating times at each failure are recorded directly, add them and paste the sums as cumulative times.

Goodness of fit: the Cramér-von Mises statistic

MIL-HDBK-189 Appendix C (20.2.1.3) tests the hypothesis that the NHPP with the power-law intensity describes the data with a Cramér-von Mises statistic computed from the unbiased shape: C_M² = 1/(12M) + Σ_i=1..M [(T_i/T)^β̄ − (2i − 1)/(2M)]², with M = N and T the end time for a time-terminated test, M = N − 1 and T = T_N for a failure-terminated test (the last failure fixes the end and is left out). Under the model the transformed times (T_i/T)^β are the order statistics of a uniform sample, so the statistic measures their departure from the uniform plotting positions. The model is rejected when the statistic exceeds the critical value of Table C-II for the chosen significance level (0.20, 0.15, 0.10, 0.05 or 0.01), indexed by M; the tool interpolates between the tabulated rows as the handbook's examples do and uses the M = 100 row beyond 100. A larger significance level is the stricter test. The handbook names two usual causes of rejection: several failures logged at the same cumulative time (group the data instead) and a discontinuity in the intensity (a batch of design changes, a change of test conditions), which shows as a kink on the cumulative-failure plot; the remedy is to fit the phase after the jump on its own, with the times measured from the jump.

Is there growth at all?

Before the shape is interpreted, the handbook's trend test (MIL-HDBK-189 5.3.5) asks whether the intensity changed at all: under a constant rate the statistic 2 Σ ln(T*/T_i) = 2N/β̂ is chi-square with 2N degrees of freedom (2(N − 1) for a failure-terminated test). Values in the lower tail mean the later intervals between failures are longer than chance allows, growth; values in the upper tail mean deterioration. The tool prints the one-sided probability and the verdict at the 10 % level. The same question is asked of the average failure rates in equal intervals of the test (the counts in K intervals against the constant rate, chi-square with K − 1 degrees of freedom), and those average rates are drawn in Figure 3 under the fitted intensity, the handbook's recommended picture: the staircase should fall and the curve should run through it.

Confidence bounds

Two families are reported. The Fisher-matrix bounds (ReliaSoft RGA Appendix C, after Nelson) invert the observed information, the negative matrix of second derivatives of the log-likelihood at the estimate, to get Var(λ̂), Var(β̂) and their covariance; each quantity of interest, being positive, is treated as log-normal, so the bounds are the estimate times exp(∓z √Var/estimate) with z the standard normal quantile at (1 + C)/2, and the variances of the cumulative and instantaneous intensity and MTBF and of N(t) follow by the delta method from their partial derivatives in β and λ. These are large-sample bounds, optimistic with few failures. Crow's bounds are exact for the model: on β, χ²_2N/(2(N − 1)) times β̂ for a time-terminated test and N χ²_2(N−1)/(2(N − 1)(N − 2)) times β̂ for a failure-terminated one; on λ and on the cumulative intensity, the chi-square interval of a Poisson count, χ²_{α/2, 2N}/2T below and χ²_{1−α/2, 2N+2}/2T (2N for failure-terminated) above; and on the instantaneous MTBF the coefficients L and U of MIL-HDBK-189 Tables C-I and C-III, L m̂ ≤ m ≤ U m̂. Crow (1977) derived those tables from the distribution of the estimate: for a time-terminated test the lower coefficient Π = 4N²/x² solves H(x|N) = (1 − γ)/2 and the upper solves H(x|N − 1) = (1 + γ)/2, where H(x|k) is the ratio of the first k terms of the series of the modified Bessel function I₁(x) to the whole function; for a failure-terminated test the coefficient c solves G(n²/c | n) = (1 ∓ γ)/2 with G(μ|n) = ∫ e^−x xⁿ⁻²/(n − 2)! P[Poisson(μ/x) ≤ n − 1] dx. The tool evaluates both equations numerically for the entered confidence rather than reading the tables, and reproduces their entries to the printed digits (checks in the test script); above N = 400 it uses the handbook's large-sample forms (1 ± z/√(2N))⁻². Because N is discrete the handbook intervals are slightly conservative.

Projection

With the fit in hand the growth curve can be extended: the instantaneous MTBF m(t) = t^1−β/(λβ) reaches a target m* at t = (λ β m*)^1/(1−β), and the expected cumulative failures at a future time are λ t^β, so the additional failures beyond those seen are λ t^β − N. The Duane line gives the same figures from b and α, t = [(1 − α) m*/b]^1/α. These assume that the factors driving the growth, the test environment, the rate of finding and fixing causes, the effectiveness of the fixes, stay what they were; they are an extrapolation of the current phase, not a forecast of the next one, and the handbook warns against extending a curve far beyond the data. A requirement that needs several times the accumulated test time calls for a new planning curve rather than a longer extrapolation. Projecting the reliability that a set of delayed fixes will deliver is a different problem, treated by the Crow extended projection model (MIL-HDBK-189C 6.3), which is not implemented here.

Reading the results

Figure 1 is the cumulative-failure plot on log-log axes: the data should scatter about a straight line, and a kink is a change of growth rate. Figure 2 is the Duane plot proper, the cumulative MTBF points with the cumulative and instantaneous lines of both models; the Fisher bounds on the Crow-AMSAA instantaneous MTBF widen to the left because the model has a singularity at t = 0. Figure 3 overlays the fitted intensity on the average failure rates, the handbook's first check. The bounds table gives the demonstrated MTBF with its handbook interval, the number that goes into the test report; the Fisher and Crow intervals should be of similar width when N is in the tens, and the Crow interval is the one to quote with few failures. The worked example on the page is the two-prototype system of MIL-HDBK-189 Appendix C (also Table VI of MIL-HDBK-189C): 27 failures in 300 hours, β̂ = 0.716, intensity 0.325 t^−0.284, demonstrated MTBF 15.5 hours with the 90 % interval 9.9 to 26.1 hours, β̄ = 0.690, Cramér-von Mises 0.091 against the critical value 0.218 at the 5 % level. (The text of MIL-HDBK-189C 6.2.2.9 quotes F = 45, β̂ = 0.826 and 8.07 hours for that table; those figures belong to a different data set than the 27 times it prints, which give the 1981 values.)

Limits

• One test phase with continuous usage (time, miles, cycles); the discrete (one-shot) tracking model, the planning models and the Crow extended projection model of MIL-HDBK-189C are not implemented.
• Failure times must be distinct in principle; repeated times are accepted but noted, and grouped data are the handbook's answer when they are many.
• Gap analysis (an interval of unreliable reporting) and a change of slope are not fitted automatically: fit the phase after the jump separately, with its times measured from the jump.
• Several systems are combined by the equivalent single system, which assumes they share the configuration and the fixes; the common-β test across systems is not offered.
• The Fisher bounds are large-sample; the Crow coefficients are exact for N ≤ 400 and the handbook's asymptotic form beyond.