This is an attempt to model the breakdown events of a system as a progressive series of failure events. This model provides an opportunity to create a simulated failure log. This in turn can be used test a hypothesis that suggests that the order of failures within a log could be used to determine the edges of a network of failures. If proven valid, service logs could be used to provide insights that would improve the service and maintenance of systems.

In this experimental model, a virtual motor is running and is subject to wear and tear that lead to motor failure. Minor problems such as loss of lubricant and accumulation of dust give way to more serious conditions as loss of bearings and overheating which in the end result in total motor failure. In the majority of cases, symptoms are treated and the motor is restored to a more normal running condition. The transition between states can be seen in the following graphic.

## Warning: package 'DiagrammeR' was built under R version 3.5.3

Simplified network of events leading to breakdown

For the sake of this study the description was removed to produce the following rendering of this problem network.

transitions = rbind(
#    dA , dB,  dC,  dD,  dE,  dF,  dG,  dH,  dI,  dJ,  dK,  dL,  dM
   c(0.5, 0.3, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), # sA
   c(0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), # sB
   c(0.5, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0), # sC
   c(0.5, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.0), # sD
   c(0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), # sE
   c(1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), # sF
   c(0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0), # sG
   c(0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.0, 0.1, 0.0, 0.0), # sH
   c(0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0), # sI
   c(0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), # sJ
   c(0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5), # sK
   c(0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), # sL
   c(0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) # sM
)

The time interval between failures is determined by Poisson distribution and the state transition is selected as per a uniform distribution. These relationships were used to develop a simulated dataset of the failure event time log.

machines = c(10000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
history = rbind(machines, machines %*% transitions)
for (i in 2:24) {
  history = rbind(history, history[i,] %*% transitions)
}
rownames(history) = c(1:25)
colnames(history) = letters[1:13]

barplot(t(history), main="Machine states",
  xlab="Number of Epocs", col=rainbow(15),
  legend = colnames(history))

If we represent the number of machines in each state as the vector \(s\) and the transition matrix as \(T\), steady state of this network was defined mathematically as \[s \times T = s\] For the values used in this experiment, steady state achieved after 25 epochs.

State Description Subsequent Step Steady state counts
A Start A, B, G 506.009
B Squeak A, C 151.803
C Rubbing A, D, H 75.901
D Lost Bearing A, E, L 7.590
E Frozen Axel A, F 0.759
F Motor Stop A 12.018
G Dust Build Up A, H 101.202
H Overheating A, 80.961
I Melted Wires A, I, K 32.385
J Electrical Short A, F 16.193
K Worn Belt A, F 8.096
L Broken Axel A, F 3.0360539
M Broken Belt A, F 4.0480720 
hist(rpois(10000,12), xlab="Interval between events", col="blue", breaks=30,
     main="Distribution of Interval")

Establishing timing sequence of events

The reduced version of this network shown below was used for the initial series of experiments.

t = rbind(
#    dA , dB,  dC,  dD,  dG
   c(0.5, 0.3, 0.0, 0.0, 0.2), # sA
   c(0.5, 0.0, 0.5, 0.0, 0.0), # sB
   c(0.5, 0.0, 0.0, 0.5, 0.0), # sC
   c(1.0, 0.0, 0.0, 0.0, 0.0), # sD
   c(0.5, 0.0, 0.0, 0.5, 0.0)) # sG
s1 = c(10000,0,0,0,0)
for (i in 1:30)
{
  s1 = s1 %*% t
  cat(i,s1,"\n")
}
## 1 5000 3000 0 0 2000 
## 2 5000 1500 1500 1000 1000 
## 3 5500 1500 750 1250 1000 
## 4 5625 1650 750 875 1100 
## 5 5437.5 1687.5 825 925 1125 
## 6 5462.5 1631.25 843.75 975 1087.5 
## 7 5487.5 1638.75 815.625 965.625 1092.5 
## 8 5482.812 1646.25 819.375 954.0625 1097.5 
## 9 5477.031 1644.844 823.125 958.4375 1096.562 
## 10 5479.219 1643.109 822.4219 959.8438 1095.406 
## 11 5479.922 1643.766 821.5547 958.9141 1095.844 
## 12 5479.457 1643.977 821.8828 958.6992 1095.984 
## 13 5479.35 1643.837 821.9883 958.9336 1095.891 
## 14 5479.467 1643.805 821.9186 958.9398 1095.87 
## 15 5479.47 1643.84 821.9024 958.8942 1095.893 
## 16 5479.447 1643.841 821.92 958.8979 1095.894 
## 17 5479.449 1643.834 821.9205 958.907 1095.889 
## 18 5479.454 1643.835 821.9171 958.905 1095.89 
## 19 5479.452 1643.836 821.9173 958.9034 1095.891 
## 20 5479.452 1643.836 821.918 958.904 1095.89 
## 21 5479.452 1643.836 821.9179 958.9043 1095.89 
## 22 5479.452 1643.836 821.9178 958.9041 1095.89 
## 23 5479.452 1643.836 821.9178 958.9041 1095.89 
## 24 5479.452 1643.836 821.9178 958.9041 1095.89 
## 25 5479.452 1643.836 821.9178 958.9041 1095.89 
## 26 5479.452 1643.836 821.9178 958.9041 1095.89 
## 27 5479.452 1643.836 821.9178 958.9041 1095.89 
## 28 5479.452 1643.836 821.9178 958.9041 1095.89 
## 29 5479.452 1643.836 821.9178 958.9041 1095.89 
## 30 5479.452 1643.836 821.9178 958.9041 1095.89
colnames(s1) = c("รค","b","c","d","g")
bp = barplot(s1,col="blue", main="Steady state of events")
text(bp,s1 + 200, format(round(s1)),xpd=TRUE)

The steady state is achieved in 22 epochs.

intervals = rpois(10000,12)
l = rep('a',10000)
times = data.frame(intervals,l)
times[1,1] = 0
for (i in 2:10000) {
  times[i,1] = times[i-1,1] + times[i,1]
}

selector = runif(10000,0,10)
b = times[selector < 3,]
b[,1] = b[,1] + rpois(length(b[,1]),5)
b[,2] = 'b'

g = times[selector > 8,]
g[,1] = g[,1] + rpois(length(g[,1]),5)
g[,2] = 'g'

selector = runif(length(g),0,10)
d1 = g[selector < 5,]
d1[,1] = d1[,1] + rpois(length(d1[,1]),5)
d1[,2] = 'd'

selector = runif(length(b[,1]),0,10)
c = b[selector < 5,]
c[,1] = c[,1] + rpois(length(c[,1]),5)
c[,2] = 'c'

selector = runif(length(c[,1]),0,10)
d = c[selector < 5,]
d[,1] = d[,1] + rpois(length(d[,1]),5)
d[,2] = 'd'

log = rbind(times,b,c,d,d1,g)
log = log[order(log[,1]),]
write(log[,2],file = "failurelog.txt")

entries = log[5001:17000,]
tab = table(entries[,2])
bp = barplot(tab,col="blue", main="Number of events in the simulated log")
text(bp,tab + 400, format(round(tab)),xpd=TRUE)

chisq.test(rbind(round(tab),s1))
## 
##  Pearson's Chi-squared test
## 
## data:  rbind(round(tab), s1)
## X-squared = 0.8609, df = 4, p-value = 0.9301