Provision spares with confidence

[shg4421]

Last year I asked a question in

https://www.eng-tips.com/viewthread.cfm?qid=437516

about a function to give the number of spares required for a lifetime purchase given only the number of operating hours and failures to date, the operating hours remaining, and the confidence level. -- something I thought would be pretty common, but apparently not. In the end, I wrote a VBA function for Excel to do just that via numerical integration:

= NumSpares(HrsSoFar, nFail, HrsToGo, ConfMin, [Stats])

E.g.,

=NumSpares(100000, 0, 150000, 95%) returns 1

=NumSpares(100000, 0, 150000, 95%, 1) returns

1 98.4% 27,109 3,949,789[/tt]

... which means 1 spare gives 98.4% confidence and shows the MTBF integration limits.

=NumSpares(100000, 0, 150000, 95%, 2) returns

0 88.9% 27,109 3,949,789
1 98.4% 27,109 3,949,789

... which means 0 spares gives 88.9% confidence and the second line is as described above. The output lists all results from 0 spares to the number that achieves the required confidence.

I'm happy to post the code if anyone want to give it a test drive and compare the results to any alternative method.

 

[IRstuff]

Sure, I'd be interested, but I'm bothered by the seeming lack of an input for the MTBF itself, although you apparently do something with MTBF.

TTFN (ta ta for now)
I can do absolutely anything. I'm an expert!

https://www.youtube.com/watch?v=BKorP55Aqvg

faq731-376 forum1529 Entire Forum list

http://www.eng-tips.com/forumlist.cfm

 

[shg4421]

The true MTBF is unknowable; that was the point. The code integrates over the 95% confidence interval (the constant [tt]MTBFConfInt[/tt] in the first function).

Here's the code:

Function NumSpares(HrsSoFar As Double, nFail As Long, _
                   HrsToGo As Double, ConfMin As Double, _
                   Optional iStats As Long) As Variant
  ' shg 2018, 2019-0706

  ' Calculates the number of spares required such that an LRU that has had
  ' nFail failures in HrsSoFar will not exceed NumSpares failures in HrsToGo
  ' (beyond the failures that have already occurred) with confidence ConfMin.

  ' HrsSoFar and nFail are assumed as time-truncated (vs failure-truncated) values

  ' iStats  Return
  ' ------  ---------------------------------------------------------------
  '   0     number of spares required
  '   1     numbers of spares, resulting confidence, and integration limits
  '   2     spares {0,1,2,...}, resulting confidence, and integration limits

  Const m           As Long = 100       ' number of MTBF steps
  Dim r             As Double           ' MTBF geometric ratio

  Dim i             As Long             ' index to all arrays
  Dim conf          As Double           ' confidence of having nSpare or fewer failures in HrsToGo

  Const MTBFConfInt As Double = 0.95    ' MTBF Confidence Interval
  Dim MTBFMin       As Double           ' low limit of integration
  Dim MTBFMax       As Double           ' high limit of integration
  Dim adMTBF(0 To m) As Double          ' geometric progression of MTBFs spanning MTBFMin to MTBFMax

  Dim adPrb1(0 To m) As Double          ' prob of LRU having exactly nFail failures in HrsSoFar at this MTBF
  Dim dPrb1Int      As Double           ' numerical (trapezoidal) integral of adPrb1 over adMTBF

  Dim adPrb2(0 To m) As Double          ' prob of LRU having nSpare failures in HrsToGo at this MTBF

  Dim col           As Collection       ' collection of results

  If HrsSoFar <= 0# Then
    NumSpares = "HrsSoFar > 0!"

  ElseIf HrsToGo <= 0# Then
    NumSpares = "HrsToGo > 0!"

  ElseIf nFail < 0& Then
    NumSpares = "nFail >= 0!"

  ElseIf ConfMin <= 0# Or ConfMin >= 1# Then
    NumSpares = "0 < ConfMin < 1!"

  Else
    ' build an array of MTBFs and the probabilities of having
    ' nFail failures in HrsSoFar at those MTBFs

    Set col = New Collection

    With WorksheetFunction
      If nFail = 0 Then
        MTBFMin = ConfidenceLimit(HrsSoFar, nFail, MTBFConfInt, 4)
        ' there is no uper limit, so pretend there was one failure
        MTBFMax = ConfidenceLimit(HrsSoFar, 1, MTBFConfInt, 5)

      Else
        MTBFMin = ConfidenceLimit(HrsSoFar, nFail, MTBFConfInt, 4)
        MTBFMax = ConfidenceLimit(HrsSoFar, nFail, MTBFConfInt, 5)
      End If

      r = (MTBFMax / MTBFMin) ^ (1# / m)

      For i = 0 To m
        adMTBF(i) = MTBFMin * r ^ i
        adPrb1(i) = .Poisson_Dist(nFail, HrsSoFar / adMTBF(i), False)
      Next i

      ' calculate the integral
      dPrb1Int = dDefInt(adMTBF, adPrb1)

      NumSpares = -1

      Do While conf < ConfMin
        NumSpares = NumSpares + 1

        ' calculate the product and integrate
        For i = 0 To m
          adPrb2(i) = adPrb1(i) * .Poisson_Dist(NumSpares, HrsToGo / adMTBF(i), False)
        Next i
        conf = conf + dDefInt(adMTBF, adPrb2) / dPrb1Int

        col.Add Item:=VBA.Array(NumSpares, conf, MTBFMin, MTBFMax)
      Loop
    End With

    Select Case iStats
      Case 0  ' just the required spares
        NumSpares = col.Item(col.Count)(0)
      Case 1  ' required spares, confidence, and mtbf limits
        NumSpares = col.Item(col.Count)
      Case 2  ' 0 to required spares and associated confidence
        Dim avOut   As Variant
        ReDim avOut(1 To col.Count)
        For i = 1 To col.Count
          avOut(i) = col.Item(i)
        Next i
        NumSpares = avOut
      Case Else
        NumSpares = Array("iStats {0|1|2}!", "", "", "")
    End Select
  End If
End Function

Function DefIntegral(vY As Variant, vX As Variant) As Variant
  ' shg 2018

  ' UDF wrapper for dDefInt

  Dim adX()         As Double
  Dim adY()         As Double

  If Not bMake1D(vX, adX, 1) Then
    DefIntegral = "X!"

  ElseIf Not bMake1D(vY, adY, 1) Then
    DefIntegral = "Y!"

  ElseIf UBound(adX) <> UBound(adY) Then
    DefIntegral = "XY!"

  Else
    DefIntegral = dDefInt(adX, adY)

  End If
End Function

Function dDefInt(adX() As Double, adY() As Double) As Double
  Dim i             As Long

  ' shg 2018
  ' VBA only

  ' returns the trapezoidal integral of adY over adX

  For i = LBound(adX) + 1 To UBound(adX)
    dDefInt = dDefInt + (adX(i) - adX(i - 1)) * (adY(i) + adY(i - 1)) / 2
  Next i
End Function

Function bMake1D(av As Variant, ad() As Double, Optional iBase As Long = 0) As Boolean
  ' shg 2014
  ' VBA only

  ' Returns a 1D iBase-based array of the values in av, which can be a
  ' column or row vector, a 1D array, or a scalar

  Dim rArea         As Range
  Dim cell          As Range
  Dim nOut          As Long
  Dim i             As Long
  Dim iLB1          As Long
  Dim iUB1          As Long
  Dim iLB2          As Long
  Dim iUB2          As Long

  If TypeOf av Is Range Then av = av.Value2

  Select Case NumDim(av)
    Case 0
      Select Case VarType(av)
        Case vbDouble, vbLong, vbInteger, vbSingle, vbByte
          ReDim ad(iBase To iBase)
          ad(iBase) = av
          bMake1D = True
      End Select

    Case 1
      iLB1 = LBound(av)
      iUB1 = UBound(av)

      ReDim ad(iBase To iBase + iUB1 - iLB1)
      For i = iLB1 To iUB1
        If VarType(av(i)) = vbDouble Then
          ad(i - iLB1 + iBase) = av(i)
        Else
          GoTo Oops
        End If
      Next i
      bMake1D = True

    Case 2
      iLB1 = LBound(av, 1)
      iUB1 = UBound(av, 1)
      iLB2 = LBound(av, 2)
      iUB2 = UBound(av, 2)

      If iUB1 <> iLB1 And iUB2 <> iLB2 Then
        ' one dimension must be a  single element
        GoTo Oops

      ElseIf iLB2 = iUB2 Then
        ' column vector
        ReDim ad(iBase To iBase + iUB1 - iLB1)

        For i = iLB1 To iUB1
          Select Case VarType(av(i, iLB2))
            Case vbDouble, vbLong, vbInteger, vbSingle, vbByte
              ad(i - iLB1 + iBase) = av(i, iLB2)
            Case Else
              GoTo Oops
          End Select
        Next i
        bMake1D = True

      Else    ' iLB1 = iUB1
        ' row vector
        ReDim ad(iBase To iBase + iUB2 - iLB2)

        For i = iLB2 To iUB2
          Select Case VarType(av(iLB1, i))
            Case vbDouble, vbLong, vbInteger, vbSingle, vbByte
              ad(i - iLB2 + iBase) = av(iLB1, i)
            Case Else
              GoTo Oops
          End Select
        Next i
        bMake1D = True

      End If
  End Select

Oops:
End Function

Private Function NumDim(av As Variant) As Long
  Dim i             As Long

  If TypeOf av Is Range Then
    If av.Count = 1 Then NumDim = 1 Else NumDim = 2

  ElseIf IsArray(av) Then
    On Error GoTo Done

    For NumDim = 0 To 6000
      i = LBound(av, NumDim + 1)
    Next NumDim
Done:
    Err.Clear
  End If
End Function

Function ConfidenceLimit(HrsSoFar As Double, nFail As Long, _
                         conf As Double, iType As Long) As Variant
  ' shg 2019

  ' Returns the specified MTBF confidence limit
  ' [URL unfurl="true"]https://reliabilityanalyticstoolkit.appspot.com/confidence_limits_exponential_distribution[/URL]
  ' [URL unfurl="true"]https://www.itl.nist.gov/div898/handbook/apr/section4/apr451.htm[/URL]

  ' iType       limit               data
  ' -----   ---------------   -----------------
  '   0     lower one-sided   failure-truncated
  '   1     lower two-sided   failure truncated
  '   2     upper two-sided   failure truncated
  '   3     lower one-sided   time-truncated
  '   4     lower two-sided   time-truncated
  '   5     upper two-sided   time-truncated

  Dim alpha         As Double

  If nFail < 0& Or conf <= 0# Or conf >= 1# Then
    ConfidenceLimit = CVErr(xlErrNum)

  Else
    alpha = 1# - conf

    With WorksheetFunction
      Select Case iType Mod 6
        Case 0  ' lower, one-sided, failure-truncated
          If nFail = 0 Then
            ConfidenceLimit = CVErr(xlErrNum)
          Else
            ConfidenceLimit = 2# * HrsSoFar / .ChiSq_Inv_RT(alpha, 2# * nFail)
          End If

        Case 1  ' lower, two-sided, failure truncated
          If nFail = 0 Then
            ConfidenceLimit = CVErr(xlErrNum)
          Else
            ConfidenceLimit = 2# * HrsSoFar / .ChiSq_Inv_RT(alpha / 2#, 2# * nFail)
          End If

        Case 2  ' upper, two-sided, failure truncated
          If nFail = 0 Then
            ConfidenceLimit = CVErr(xlErrNum)
          Else
            ConfidenceLimit = 2# * HrsSoFar / .ChiSq_Inv_RT(1# - alpha / 2#, 2# * nFail)
          End If

        Case 3  ' lower, one-sided, time-truncated
          ConfidenceLimit = 2# * HrsSoFar / .ChiSq_Inv_RT(alpha, 2# * nFail + 2#)

        Case 4  ' lower, two-sided, time-truncated
          ConfidenceLimit = 2# * HrsSoFar / .ChiSq_Inv_RT(alpha / 2#, 2# * nFail + 2#)

        Case 5  ' upper, two-sided, time-truncated
          If nFail = 0 Then
            ConfidenceLimit = "(none)"
          Else
            ConfidenceLimit = 2# * HrsSoFar / .ChiSq_Inv_RT(1 - alpha / 2#, 2# * nFail)
          End If
      End Select
    End With
  End If
End Function

 

[IRstuff]

The calculated upper bound seems a bit high to me. Shouldn't 95% reliability at 100,000 hrs result in an MTBF of 1.95 million hr, just from -100000/ln(0.95)?

3949789 hr MTBF results in exp(-100000/3949789) = 97.5%, which looks like that's because of the two-side choice

TTFN (ta ta for now)
I can do absolutely anything. I'm an expert!

https://www.youtube.com/watch?v=BKorP55Aqvg

faq731-376 forum1529 Entire Forum list

http://www.eng-tips.com/forumlist.cfm

 

[shg4421]

3949789 hr MTBF results in exp(-100000/3949789) = 97.5%, which looks like that's because of the two-side choice

Right; it's a 95% two-sided confidence interval -- 2.5% above and below. The function [tt]ConfidenceLimit[/tt] computes one- or two-sided limits for time- or failure-truncated data.

Thanks for looking.

 

[IRstuff]

I guess the question is whether that's appropriate; seems to me that reliability is always one-sided in a practical sense, since we never really care about the lower side, as that's often buried by the infant mortality.

TTFN (ta ta for now)
I can do absolutely anything. I'm an expert!

https://www.youtube.com/watch?v=BKorP55Aqvg

faq731-376 forum1529 Entire Forum list

http://www.eng-tips.com/forumlist.cfm

 

[shg4421]

The function integrations over a ratiometric scale from the low limit to the high limit. I reckon it could integrate from a smaller value, but not zero. I don't think it would change the answers out to several decimal places, but it's a valid point, thanks.

I'd be very interested if you compared and posted results with any other source you have.