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.

• flodel: I feel you are going to introduce a bias by assuming that the generator is up on January 1st of every year. So instead of running this same code 100000 times, I'd suggest you run it once for 8736 * 100000 hours, then split that long vector by year (or put that vector into a matrix with 8736 rows.) This way you are only requiring that the generator be up on January 1st of the 1st year, a much weaker assumption.

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]

• kelamahim: For sure more elegant solution
• kelamahim: and much faster
• kelamahim: Will simulate to compare the two approaches if the ON and OFF times are the same
• slamballais: I tried it just now. They give the same results, except for some rounding error. You used cumsum(mix) and then you rounded that, whereas I simply took mix and rounded that. The result is that sometimes your fails/repairs are 1 hour earlier/later. Other than that, the results are identical, and the sum of capacity is still the same (try all.equal(sum(a1),sum(out)).
• kelamahim: Exactly! Just tried out and the results are basically the same. Thank you for the simplification.

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.