Assuming the time between demand periods is exponentially distributed, the probability that machine $n$ is in a demand period at some time $t$ in the future is
$$p_n=\frac{m_n}{60*h_n}$$
The probability that the three machines are in a demand state of, e.g., {1, 1, 0} (where 0 = no demand, and 1 = demand) is
$$p_1p_2(1-p_3)$$
Calculating the full set of probabilities in R:
Brute-force with expand.grid
library(data.table)
library(Rfast) # for the colprods function
f <- c(20, 25, 40)
h <- c(1, 2, 1)
m <- c(10, 20, 5)
p <- m/(60*h)
states <- t(expand.grid(setNames(rep(list(0:1), length(f)), paste0("n=", seq_along(f)))))
setorder(
cbind(
t(states),
data.frame(
total_demand = colSums(states*f),
probability = colprods(states*p + (1 - states)*(1 - p))
)
), total_demand
)[]
#> n=1 n=2 n=3 total_demand probability
#> 1 0 0 0 0 0.636574074
#> 2 1 0 0 20 0.127314815
#> 3 0 1 0 25 0.127314815
#> 5 0 0 1 40 0.057870370
#> 4 1 1 0 45 0.025462963
#> 6 1 0 1 60 0.011574074
#> 7 0 1 1 65 0.011574074
#> 8 1 1 1 85 0.002314815
Iterative approach
If there are many machines whose flow rates share a common factor (e.g., the flow rates are all integer values), a more efficient algorithm would be to compute the total flow rate by iteratively adding machines.
library(data.table)
f <- c(20, 25, 40)
h <- c(1, 2, 1)
m <- c(10, 20, 5)
p <- cbind(m/60/h, 1 - m/60/h)
dt <- data.table(total_demand = c(0, f[1]), probability = p[1,])
for (i in 2:length(f)) {
dt <- dt[
, .(
probability = c(outer(probability, p[1,])),
total_demand = c(total_demand, total_demand + f[i])
)
][, .(probability = sum(probability)), total_demand]
}
dt
#> total_demand probability
#> 1: 0 0.636574074
#> 2: 20 0.127314815
#> 3: 25 0.127314815
#> 4: 45 0.025462963
#> 5: 40 0.057870370
#> 6: 60 0.011574074
#> 7: 65 0.011574074
#> 8: 85 0.002314815
Using the FFT
For the problem in the question, all flow rates are integer multiples of 5, so the total flow rate can take values only from $\{0,5,10,...,(\sum{f_n})=85\}$, so only 18 bins are needed. If there are a large number of machines whose flow rates do not share a common factor, the number of bins should be selected based on the desired precision, since the output will be discretized into bins ranging from 0 to the maximum possible flow rate.
library(data.table)
f <- c(20, 25, 40)
h <- c(1, 2, 1)
m <- c(10, 20, 5)
n.bins <- 18L
x.hat <- rep(1, n.bins)
p <- m/60/h
idx <- 1 + round((n.bins - 1)/(sumf <- sum(f))*f)
for (i in seq_along(f)) {
y <- numeric(n.bins)
y[idx[i]] <- p[i]
y[1] <- 1 - p[i]
x.hat <- fft(y)*x.hat
}
data.table(
total_demand = seq(0, sumf, length.out = n.bins),
probability = zapsmall(Re(fft(x.hat, TRUE)))/n.bins
)[probability > 0]
#> total_demand probability
#> 1: 0 0.636574056
#> 2: 20 0.127314833
#> 3: 25 0.127314833
#> 4: 40 0.057870389
#> 5: 45 0.025462944
#> 6: 60 0.011574056
#> 7: 65 0.011574056
#> 8: 85 0.002314833