How to write this cost equation?

Question

I have a mean time to failure variable that is Weibull-distributed with shape parameter less than $1$. Whenever a failure occurs, corrective action of cost $c_c$ is performed. From time to time, preventive maintenance of cost $c_m$ is performed, which resets the time in the Weibull distribution.

My question is, how do I write an annual cost equation using this?

This is what I have so far. My expected number of failures in a year is $\frac{365}{\mu}$, where $\mu$ is the mean of the Weibull. Therefore, the cost for corrective action is:

$$c_c * \frac{365}{\mu}$$

When I add in the preventive maintenance, it becomes something like:

$$c_c * \frac{365}{\mu} + c_m * f_m$$

Where $f_m$ is the frequency of doing preventive maintenance in a year.

But the problem is the frequency of preventive maintenance affects the number of times corrective action is done (in this case it increases it). So a more correct formula would be:

$$c_c * fcn(f_m) + c_m * f_m$$

Where $fcn(f_m) = \frac{365}{\mu} + something$, so for every increase in $f_m$, the function value increases. Obviously, $something$ should have $f_m$ in there somewhere, but I don't know where. The idea is that every increase in frequency of preventive maintenance increases the expected number of times corrective maintenance has to be done, since the function is being reset to time zero and there is a large infant mortality in this distribution.

I tried thinking about it in terms of just probability, but it's difficult. Say preventive maintenance is done some time $t_p$ after corrective maintenance. Then the probability it will have broken down in that time is $Pr(t<T=t_p)$ and the probability for the Weibull time to failure resets to time zero.

But then I'm stuck, because how can I model the corrective actions in probability terms? Plus, I don't know how to convert that probability into frequency. I dont think it's as simple as $f = \frac{365}{t_p}$ because preventive maintenance happens on fixed intervals with respect to the previous preventive maintenance, not with respect to the previous corrective action, so $t_p$ always varies. And since the occurrence of one preventive maintenance resets the time to failure function, then all corrective actions are shifted.

Any help or advice would be appreciated!

Answer 1

Here is a suggested solution. Note that if failure times followed an exponential distribution then we would have a Poisson count process for the number of failures between each preventative maintenance action. Since failure times do not follow an exponential distribution but a Weibull distribution we need to look elsewhere.

A counting process with Weibull distribution for time between events will not be possible to place in closed form because the convolution of the PDF and CDF of the Weibull distribution has no closed form.

I found this paper:

Bradlow, Eric and Fader, Peter and Adrian, Moshe and McShane, Blakeley B., Count Models Based on Weibull Interarrival Times (January 2006). Available at SSRN: http://ssrn.com/abstract=729886 or http://dx.doi.org/10.2139/ssrn.729886

Which suggests that the Taylor series expansion may be used and the authors identify recursion formulae which are useful.

Assume the Weibull distribution is parameterised by $\lambda$ and $c$ (as in the paper) such that the CDF is $(1-\exp(-\lambda t^c))$ (if you have another parameterisation you can convert as necessary).

The authors calculate first Taylor expansions for the CDF and PDF:

$$ \begin{align} F(t) &= \sum_{j=1}^\infty{\frac{(-1)^{j+1}(\lambda t^c)^j}{\Gamma(j+1)}} \\ f(t) &= \sum_{j=1}^\infty{\frac{(-1)^{j+1}cj\lambda^j t^{cj-1}}{\Gamma(j+1)}} \end{align} $$

They then define the following:

$$ \begin{align} \alpha_j^0 &= \frac{\Gamma(cj+1)}{\Gamma(j+1)}, \quad j=0,1,2,\ldots \\ \alpha_j^{n+1} &= \sum_{m=n}^{j-1}{\alpha_m^n \frac{\Gamma(cj-cm+1)}{\Gamma(j-m+1)}} \end{align} $$

And use these for the following results:

$$ \begin{align} \Pr(N(t)=n) &= \sum_{j=n}^\infty{\frac{(-1)^{j+n}(\lambda t^c)^j \alpha_j^n}{\Gamma(cj+1)}} \\ E(N(t)) &= \sum_{n=1}^\infty{\sum_{j=n}^\infty{\frac{n(-1)^{j+n}(\lambda t^c)^j \alpha_j^n}{\Gamma(cj+1)}}} \end{align} $$

SO! How does this help you?

Well, you need to calculate $E(N(T_m))$ where $T_m$ is the interval between preventative maintenance to give you the expected number of corrective maintenance actions for each preventative maintenance action. You then scale up by how many preventative maintenance actions there are in a year and apply the relevant costs to corrective and preventative maintenance.

You can sensibly sum $n$ only up to whatever is an extremely unlikely number of failures within time $T_m$. Choosing increasing values of $j$ will result in progressively more accurate estimation, so this should be summed up as appropriate for your desired accuracy.

EDIT

As stated in a comment below, the notation employed by the authors is a little confusing, here is a clarification with some of my own notation:

1. Define base case $A(0, j) = \frac{\Gamma(cj+1)}{\Gamma(j+1)}$
2. Define recursive case $A(n+1, j) = \sum_{m=n}^{j-1}{A(n, m)\frac{\Gamma(cj-cm+1)}{\Gamma(j-m+1)}}$
3. Obtain results:
   • $C_n(t) = \Pr(N(t)=n) = \sum_{j=n}^\infty{\frac{(-1)^{j+n}(\lambda t^c)^j A(n, j)}{\Gamma(cj+1)}}$
   • $E[N(t)] = \sum_{n=0}^\infty{n C_n(t)}$

I have not yet coded this as an Excel spreadsheet (though it is doubtless possible), instead here is Python v3 code:

from math import exp
from scipy.special import gamma, gammaln

rate = 1.5
shape = 0.8

Acache = dict()

def A(n, j):
    if not (n, j) in Acache:
        if n == 0:
            Acache[(n, j)] = exp(gammaln(shape*j+1) - gammaln(j+1))
        else:
            Acache[(n, j)] = sum( A(n-1, m) * exp(gammaln(shape*(j-m)+1) - gammaln(j-m+1)) for m in range(n-1, j) )
    return Acache[(n, j)]

def Cn(t, n, J=20):
    return sum( ( (-1)**(j+n) * (rate * t**shape)**j * A(n, j) )/gamma(shape*j+1) for j in range(n, n+J) )

def E(t, J=20, N=20):
    return sum( (n * Cn(t, n, J)) for n in range(0, N) )

C = list(Cn(1, n) for n in range(0, 20))

print('n', 'Cn(1)', sep='\t')
for n in range(0, 20):
    print(n, C[n], sep='\t')

print(r'\sum_{n=0}^{19}{C_n(1)} =', sum(C))

print('E[N(1)] =', E(1))

Which gives the following output:

n   Cn(1)
0   0.223130160148
1   0.285693698536
2   0.232154412273
3   0.142130392306
4   0.0707688290882
5   0.0299340192311
6   0.0110641801411
7   0.00364513700374
8   0.00108625766843
9   0.000296143808589
10  7.45318952329e-05
11  1.74439294931e-05
12  3.81998309182e-06
13  7.86741813662e-07
14  1.53063575444e-07
15  2.82382329525e-08
16  4.95652671252e-09
17  8.30172205633e-10
18  1.33028122034e-10
19  2.0441696595e-11
\sum_{n=0}^{19}{C_n(1)} = 0.999999999996
E[N(1)] = 1.71339144029