Generating time series for generator ON or OFF

Question

Want to create time series for the nuclear power plant performance using the equation for Time to fail and time to repair for 8736 hours in the year so that I have the time series in hours when the generator is operating in when is not. The starting condition is that the generator is operating on the first hour.

For sure there is a more elegant solution for simulating this I'm just not able to find it. Any comment or help will be appreciated. The code is working but I think there is a better solution

TTF<-2940 # MEDIUM TIME TO FAIL(hours)
MTTR<-60 # MEDIUM TIME TO REPAIR (hours)
TTF<--MTTF*log(runif(100))# equation for Time to fail 
TTR<--MTTR*log(runif(100))# equation for Time to repair
mix<-rep(0,length(TTF)+length(TTR))
sw<-rep(0,length(TTF)+length(TTR))
for(i in 1:length(TTF)){
mix[2*i-1]<-TTF[i]
sw[2*i-1]<-1
mix[2*i]<-TTR[i]
}
cmix<-cumsum(mix)
ccmix<-cbind(cmix[1:which(cmix>8736)],sw[1:which(cmix>8736)])
ccmix[dim(ccmix)[1],1]<-8736


G1<-round(ccmix)
# transform binary values
G1[,2][G1[,2] == 1] <- 400 # 400 MW is the capacity of power plant

G1 <- cbind(G1, c(G1[1,1], diff(G1[,1])))
a1 <- rep(G1[,2], G1[,3]) ## GENERATING 8736 Values

So the desired output are 8736 values of 400 when is ON and 0 when is OFF. Finally I need to simulate 100000 years of generation.

Answer 1

I don't know of any method, but here are some small improvements:

(1) You want mix to be TTF and TTR intertwined. An easy, vectorized way to do that is:

mix <- c(rbind(TTF,TTR))

(2) You want sw to be an alternating vector of the capacity and 0. This can be vectorized as well:

cap <- 400
sw <- rep(c(cap,0), length(TTF))

(3) Next, simply find the index which is the first value above 8736 and rep that:

final <- 8736
ind <- which.max(cumsum(mix) > final)
out <- rep(sw[1:ind], round(mix)[1:ind])[1:final]

So full code would be:

MTTF <- 2940
MTTR <- 60
final <- 8736
cap <- 400
TTF <--MTTF*log(runif(100))
TTR <--MTTR*log(runif(100))

mix <- c(rbind(TTF,TTR))
sw <- rep(c(cap,0), length(TTF))
ind <- which.max(cumsum(mix) > final)
out <- rep(sw[1:ind], round(mix)[1:ind])[1:final]

Answer 2

The main problem I see is that it does not guaranty you will have data for 8736 hours. Of course one can increase the samples but that will waste a lot.

The way you describe it I think it should be written as Markov process in the form state [i] = f state[i-1]. That will get rid of all the sorting values together.

Though I am not very familiar with the probability distribution you use. To use that some complicated transformation may be required.