Context

In order to illustrate the package SMM we introduce a practical example taken from Barbu and Limnios (2008), where a semi-Markov modeling is appropriate. Let ℱ be a textile factory. To avoid river pollution, the factory waste is treated in the unit 𝒰 before being thrown in the river ℛ. In order for the factory not to be stopped if a failure occurs in the treatment unit 𝒰, the waste is stocked in a tank 𝒯. If 𝒰 is repaired before 𝒯 is full, then the factory starts to work again normally and we suppose that the tank is emptied instantaneously (before another failure of 𝒰). Otherwise, the whole factory has to stop and a certain time is necessary to restart it. The figure below describes the waste treatment system of the textile factory.

Waste treatment for a textile factory

As the output of the factory is variable, the filling time of the tank is a random variable. Moreover, if we suppose that the filling time is deterministic, we can still make the assumption that we have a random time that has a Dirac distribution.

We propose a semi-Markov modelization for the evolution of the described system. Let E = {"1", "2", "3"} be the set of possible states of the system, where we set:

• State 1 - everything is working, that is, ℱ and 𝒰 are in good state and 𝒯 is empty;
• State 2 - 𝒰 has failed, but 𝒯 is not yet full, so the factory is still working;
• State 3 - 𝒰 has failed and 𝒯 is full, so the factory is not working.

The possible transitions between states are given in the following figure:

Three-state discrete-time semi-Markov system modelization

The system is defined by:

• The initial distribution α = (α₁, α₂, α₃);

• The transition matrix p of the embedded Markov chain (J_n)_{n ∈ ℕ} $\boldsymbol{p} = \begin{pmatrix} 0 & 1 & 0 \\ a & 0 & b \\ 1 & 0 & 0 \end{pmatrix}$,

• The conditional sojourn time distributions $\boldsymbol{f(k)} = \begin{pmatrix} 0 & f_{12}(k) & 0 \\ f_{21}(k) & 0 & f_{23}(k) \\ f_{31}(k) & 0 & 0 \end{pmatrix}, \ k \in \mathbb{N}$,

  where,

  • f₁₂ is the distribution of the failure time of the treatment unit 𝒰;

  • f₂₁ is the distribution of the repairing time of 𝒰;

  • f₂₃ is the distribution of the filling time of the tank 𝒯;

  • f₃₁ is the distribution of time needed to restart the entire factory ℱ, after it has been shut down.

This is an irreducible discrete-time semi-Markov system, or, equivalently, a Markov renewal system. We chose the following sojourn time distributions:

• f₁₂ is the geometric distribution on ℕ^* of parameter p, 0 < p < 1, i.e., f₁₂(k) = p(1 − p)^k − 1, k ∈ ℕ^*.

• f₂₁ = W_q1, b1, f₂₃ = W_q2, b2, f₃₁ = W_q3, b3 are discrete-time first-type Weibull distributions (Nakagawa and Osaki (1975)), defined by: W_q, b(0) = 0, W_q, b(k) = q^{(k − 1)b} − q^kb, k ∈ ℕ^*

The choice of discrete-time Weibull distribution is motivated by its use in reliability theory, due to its flexibility in modeling failure rates. In fact, our main purpose in this document is to compute and estimate the main reliability indicators of this system.

Numerical application

Let us continue the example of the textile factory presented above. We take the initial distribution α = (1, 0, 0), the transition matrix of the embedded Markov chain (J_n)_{n ∈ ℕ}, $$ \boldsymbol{p} = \begin{pmatrix} 0 & 1 & 0 \\ 0.95 & 0 & 0.05 \\ 1 & 0 & 0 \end{pmatrix} $$ and the following conditional sojourn time distributions:

• f₁₂ is the geometric distribution on ℕ^* with parameter p = 0.8;

• f₂₁ = W_q1, b1, f₂₃ = W_q2, b2, f₃₁ = W_q3, b3 are discrete-time Weibull distributions with parameters: q₁ = 0.3, b₁ = 0.5, q₂ = 0.5, b₂ = 0.7, q₃ = 0.6, b₃ = 0.9.

states <- c("1", "2", "3") # State space

alpha <- c(1, 0, 0) # Initial distribution

p <- matrix(data = c(0, 1, 0, 
                     0.95, 0, 0.05, 
                     1, 0, 0), nrow = 3, byrow = TRUE) # Transition matrix

distr <- matrix(c(NA, "geom", NA, 
                  "dweibull", NA, "dweibull", 
                  "dweibull", NA, NA), 
                nrow = 3, ncol = 3, byrow = TRUE) # Distribution matrix

param1 <- matrix(c(NA, 0.8, NA, 
                   0.3, NA, 0.5,
                   0.6, NA, NA), 
                 nrow = 3, ncol = 3, byrow = TRUE)

param2 <- matrix(c(NA, NA, NA, 
                   0.5, NA, 0.7,
                   0.9, NA, NA), 
                 nrow = 3, ncol = 3, byrow = TRUE)

parameters <- array(c(param1, param2), c(3, 3, 2))

factory <- smmparametric(states = states, init = alpha, ptrans = p, 
                         type.sojourn = "fij", distr = distr, param = parameters)

M <- 10000
seq <- simulate(object = factory, nsim = M)

estimate <- fitsmm(sequences = seq, states = states, type.sojourn = "fij")
## Some transitions from state i to state j are not observed: (i = "3" to j = "2"), (i = "1" to j = "3")
## The probabilities of the initial state(s) "2", "3" are 0.

print(x = estimate$ptrans, digits = 2)
##      1 2     3
## 1 0.00 1 0.000
## 2 0.95 0 0.048
## 3 1.00 0 0.000

plot(x = estimate, i = "2", j = "3", type = "l", col = "blue")

lines(x = 1:estimate$kmax, y = ddweibull(x = 1:estimate$kmax, q = 0.5, beta = 0.7), 
       col = "red", pch = "x")

legend(x = "topright", 
       legend = c("True value", "Estimate"), 
       col = c("red", "blue"), lty = c(1, 1))

k <- 300
upstates <- c("1", "2") # Working states of the semi-Markov system

trueReliab <- reliability(x = factory, k = k, upstates = upstates)
estReliab <- reliability(x = estimate, k = k, upstates = upstates)

plot(x = 0:k, y = trueReliab[, 1], type = "l", cex = 2.5, ylim = c(0, 1), 
     col = "red", main = "Reliability", xlab = "k", ylab = "R(k)")

lines(x = estReliab[, 1], col = "blue")
lines(x = estReliab[, 3], lty = 4, col = "blue")
lines(x = estReliab[, 4], lty = 4, col = "blue")
legend(x = "topright", 
       legend = c("True value", "Estimated value", "95% confidence interval"), 
       col = c("red", "blue", "blue"), lty = c(1, 1, 4))

trueAvail <- availability(x = factory, k = k, upstates = upstates)
estAvail <- availability(x = estimate, k = k, upstates = upstates)

plot(x = 0:k, y = trueAvail[, 1], type = "l", cex = 2.5, ylim = c(0.95, 1), 
     col = "red", main = "Availability", xlab = "k", ylab = "A(k)")

lines(x = estAvail[, 1], col = "blue")
lines(x = estAvail[, 3], lty = 4, col = "blue")
lines(x = estAvail[, 4], lty = 4, col = "blue")
legend(x = "topright", 
       legend = c("True value", "Estimated value", "95% confidence interval"), 
       col = c("red", "blue", "blue"), lty = c(1, 1, 4))

trueBMP <- failureRate(x = factory, k = k, upstates = upstates)
estBMP <- failureRate(x = estimate, k = k, upstates = upstates)

plot(x = 0:k, y =  trueBMP[, 1], type = "l", cex = 2.5, ylim = c(0, 0.025), 
     col = "red", main = "BMP-failure rate", xlab = "k", ylab = bquote(lambda(k)))

lines(x = estBMP[, 1], col = "blue")
lines(x = estBMP[, 3], lty = 4, col = "blue")
lines(x = estBMP[, 4], lty = 4, col = "blue")
legend(x = "topright", 
       legend = c("True value", "Estimated value", "95% confidence interval"), 
       col = c("red", "blue", "blue"), lty = c(1, 1, 4))

trueRG <- failureRate(x = factory, k = k, upstates = upstates, failure.rate = "RG")
estRG <- failureRate(x = estimate, k = k, upstates = upstates, failure.rate = "RG")

plot(x = 0:k, y =  trueRG[, 1], type = "l", cex = 2.5, ylim = c(0, 0.03), 
     col = "red", main = "RG-failure rate", xlab = "k", ylab = "r(k)")

lines(x = estRG[, 1], col = "blue")
lines(x = estRG[, 3], lty = 4, col = "blue")
lines(x = estRG[, 4], lty = 4, col = "blue")
legend(x = "topright", 
       legend = c("True value", "Estimated value", "95% confidence interval"), 
       col = c("red", "blue", "blue"), lty = c(1, 1, 4))

trueMTTF <- mttf(x = factory, upstates = upstates)
estMTTF <- mttf(x = estimate, upstates = upstates)

print(trueMTTF)
##       mttf   sigma2
## 1 67.19099 307414.6
## 2 65.94099 307372.0
print(estMTTF)
##       mttf   sigma2      lwr    upper
## 1 69.06976 330291.4 57.80790 80.33161
## 2 67.83568 330252.4 56.57449 79.09687

trueMTTR <- mttr(x = factory, upstates = upstates)
estMTTR <- mttr(x = estimate, upstates = upstates)