rm(list=ls())

options(digits=4)

QueuingFullCapacity <- function(lamda, mu, maxNumOfCustomer) {
  rho <- lamda / mu
  
  sumTemp <- 0
  for (i in 0:3) {
    sumTemp <- sumTemp + (rho ^ i) * (factorial(maxNumOfCustomer) / factorial(maxNumOfCustomer - i))
  }

  P0 <- 1/sumTemp
  Lq <- maxNumOfCustomer - ((lamda + mu) * (1 - P0)) / lamda
  L <- Lq + (1 - P0)
  
  avgLamda <- (maxNumOfCustomer - L) * lamda
  
  W <- L / avgLamda
  Wq <- Lq / avgLamda
  
  ans <- c(rho, P0, L, W, Lq, Wq)
  ans
}

queuing <- function(maxNumOfCustomer, numOfCustomer, rho, P0) {
  Pn <- P0 * (rho ^ numOfCustomer) * factorial(maxNumOfCustomer) / factorial(maxNumOfCustomer - numOfCustomer)
  Pn
}

avgNumOfCustomer <- 1/6
numOfRequestCanHold <- 1/2

ans <- QueuingFullCapacity(lamda=avgNumOfCustomer, mu=numOfRequestCanHold, maxNumOfCustomer=3)
P0 <- ans[2]

print(sprintf("維修員閒置機率=%.4f", P0))
print(sprintf("服務系統內平均客戶數=%.4f", ans[3]))
print(sprintf("服務系統內平均處理時間=%.4f", ans[4]))
print(sprintf("等候線內平均客戶數=%.4f", ans[5]))
print(sprintf("等候線內平均處理時間=%.4f", ans[6]))

for (i in 1:3) {
  Pn <- queuing(maxNumOfCustomer=3, numOfCustomer=i, rho=ans[1], P0=P0)  
  print(sprintf("設備 %d 部維修機率=%.4f", i, Pn))
}

Pn <- queuing(maxNumOfCustomer=3, numOfCustomer=1, rho=ans[1], P0=P0)
print(sprintf("設備等候維修機率=%.4f", 1 - P0 - Pn))

[1] 維修員閒置機率=0.3462
[2] 服務系統內平均客戶數=1.0385
[3] 服務系統內平均處理時間=3.1765
[4] 等候線內平均客戶數=0.3846
[5] 等候線內平均處理時間=1.1765
[6] 設備 1 部維修機率=0.3462
[7] 設備 2 部維修機率=0.2308
[8] 設備 3 部維修機率=0.0769
[9] 設備等候維修機率=0.3077