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_m^l &= \sum_{m=0}^{l-1}{\frac{\Gamma(cm+1)\Gamma(cl-cm+1)}{\Gamma(m+1)\Gamma(l-m+1)}} \\ \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.
This can all be implemented in Excel (use =EXP(GAMMALN(X))
for $\Gamma(X)$). I would suggest creating $J \times N$ grids (where $J$ and $N$ are your choice for how high to sum over $j$ and $n$) for $\Pr(N(T_m)=n)$ and calculation of $E(N(T_m))$, drawing from vectors created for $\alpha_m^n$ and $\alpha_j^n$.
I'm afraid all the legwork is left for you to do, but hopefully this is a sufficiently comprehensive answer to get you on your way!
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:
- Define base case $A(0, j) = \frac{\Gamma(cj+1)}{\Gamma(j+1)}$
- Define recursive case $A(n+1, j) = \sum_{m=n}^{j-1}{A(n, m)\frac{\Gamma(cj-cm+1)}{\Gamma(j-m+1)}}$
- 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