Lieral a
library(markovchain)
## Package: markovchain
## Version: 0.8.6
## Date: 2021-05-17
## BugReport: https://github.com/spedygiorgio/markovchain/issues
# X = # rollos en la prensa en t
# Y = # rollos en recubrimiento en t
# Z = # rollos en satinado en t
#Capacidades
capPre<-3
capRec<-2
capSat<-2
#Tasas de servicio
tasaPre<-0.104
tasaRec<-1/8
tasaSat<-1/5
#Tasa de llegada
tasaLlegada<-3*1/60
#Estados de conjunta
estados<-c()
for (i in 0:capPre) {
for (j in 0:capRec) {
for (k in 0:capSat) {
estados<-c(estados,paste(i,j,k,sep = ","))
}
}
}
estados
## [1] "0,0,0" "0,0,1" "0,0,2" "0,1,0" "0,1,1" "0,1,2" "0,2,0" "0,2,1" "0,2,2"
## [10] "1,0,0" "1,0,1" "1,0,2" "1,1,0" "1,1,1" "1,1,2" "1,2,0" "1,2,1" "1,2,2"
## [19] "2,0,0" "2,0,1" "2,0,2" "2,1,0" "2,1,1" "2,1,2" "2,2,0" "2,2,1" "2,2,2"
## [28] "3,0,0" "3,0,1" "3,0,2" "3,1,0" "3,1,1" "3,1,2" "3,2,0" "3,2,1" "3,2,2"
# Matriz de probabilidad --------------------------------------------------
matQ<-matrix(0,nrow = length(estados), ncol=length(estados),dimnames = list(estados,estados))
for (fila in estados) {
for (columna in estados) {
i<-as.numeric(strsplit(fila,",")[[1]][1])
j<-as.numeric(strsplit(fila,",")[[1]][2])
k<-as.numeric(strsplit(fila,",")[[1]][3])
iprima<-as.numeric(strsplit(columna,",")[[1]][1])
jprima<-as.numeric(strsplit(columna,",")[[1]][2])
kprima<-as.numeric(strsplit(columna,",")[[1]][3])
#Rollo entra al sistema
if (iprima==i+1&jprima==j & kprima==k & i<capPre) {
matQ[fila,columna]<-tasaLlegada
}
#Rollo sale de prensa a recubrimiento
if (iprima==i-1&jprima==j+1 & kprima==k & i>0 & j<capRec) {
matQ[fila,columna]<-tasaPre
}
#Rollo sale de recubrimiento a alisado
if (iprima==i &jprima==j-1 & kprima==k+1 & j>0 & k<capSat) {
matQ[fila,columna]<-tasaRec*min(j,2)
}
#Rollo sale de alisado
if (iprima==i &jprima==j & kprima==k-1 & k>0 ) {
matQ[fila,columna]<-tasaSat*11/12
}
}
}
diag(matQ)<-(-rowSums(matQ))
rowSums(matQ)
## 0,0,0 0,0,1 0,0,2 0,1,0 0,1,1
## 0.000000e+00 1.387779e-17 1.387779e-17 1.387779e-17 1.387779e-17
## 0,1,2 0,2,0 0,2,1 0,2,2 1,0,0
## 1.387779e-17 1.387779e-17 1.387779e-17 1.387779e-17 0.000000e+00
## 1,0,1 1,0,2 1,1,0 1,1,1 1,1,2
## -2.775558e-17 -2.775558e-17 -2.775558e-17 -2.775558e-17 -2.775558e-17
## 1,2,0 1,2,1 1,2,2 2,0,0 2,0,1
## 1.387779e-17 1.387779e-17 1.387779e-17 0.000000e+00 -2.775558e-17
## 2,0,2 2,1,0 2,1,1 2,1,2 2,2,0
## -2.775558e-17 -2.775558e-17 -2.775558e-17 -2.775558e-17 1.387779e-17
## 2,2,1 2,2,2 3,0,0 3,0,1 3,0,2
## 1.387779e-17 1.387779e-17 0.000000e+00 1.387779e-17 1.387779e-17
## 3,1,0 3,1,1 3,1,2 3,2,0 3,2,1
## 1.387779e-17 1.387779e-17 1.387779e-17 0.000000e+00 0.000000e+00
## 3,2,2
## 0.000000e+00
#CRear cadena de markov
cadena<-new(Class="ctmc",states=estados,generator=matQ)
Lieral B
# Literal B ---------------------------------------------------------------
# Tiempo de ocupacion -----------------------------------------------------
#1 Construir EMC
EMC<-matrix(0,nrow = length(estados),ncol = length(estados),dimnames = list(estados,estados))
for(i in estados){
for(j in estados){
if(i!=j){
EMC[i,j]<-matQ[i,j]/(-matQ[i,i])
}
}
}
rowSums(EMC)
## 0,0,0 0,0,1 0,0,2 0,1,0 0,1,1 0,1,2 0,2,0 0,2,1 0,2,2 1,0,0 1,0,1 1,0,2 1,1,0
## 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1,1,1 1,1,2 1,2,0 1,2,1 1,2,2 2,0,0 2,0,1 2,0,2 2,1,0 2,1,1 2,1,2 2,2,0 2,2,1
## 1 1 1 1 1 1 1 1 1 1 1 1 1
## 2,2,2 3,0,0 3,0,1 3,0,2 3,1,0 3,1,1 3,1,2 3,2,0 3,2,1 3,2,2
## 1 1 1 1 1 1 1 1 1 1
EMC2<-generatorToTransitionMatrix(cadena@generator)
#2. Hallar n equivalente
#Elemnetos importantes
t<-7*60*60
alpha<-rep(0,length(estados))
names(alpha)<-estados
alpha["0,0,0"]<-1
vectorR<-(-diag(matQ))
vectorL<-(1/vectorR)
#Inicializar
E_t<-alpha%*%vectorL
aux<-EMC
k<-1
#Hallar n equivalente
while(E_t<t){
k=k+1
E_t<-E_t+alpha%*%aux%*%vectorL
aux<-aux%*%EMC
}
n<-k-1
n
## [1] 4692
#3.Hallar matriz de ocupacion M de EMC en n
M<-matrix(0,ncol=length(estados),nrow=length(estados),dimnames=list(estados,estados))
diag(M)<-1
aux<-EMC
for (i in 1:n) {
M<-M+aux
aux<-aux%*%EMC
}
#Mumero de transiciones/veces en el estado
M["0,0,0","3,2,2"]
## [1] 0.6023113
#4. Tiempo en segundos que dura en 3,2,2
M["0,0,0","3,2,2"]*vectorL["3,2,2"]
## 3,2,2
## 3.285334
(alpha%*%M*vectorL)["3,2,2"]
## [1] NA
rowSums(alpha%*%M*vectorL)
## [1] 25201.85
Literal C y D
# Literal C numero de rollos en estado estable ----------------------------
probaLP<-steadyStates(cadena)
NumeroRollos<-0
contador<-0
for (e in estados) {
contador<-contador+1
i<-as.numeric(strsplit(e,split = ",")[[1]][1])
j<-as.numeric(strsplit(e,split = ",")[[1]][2])
k<-as.numeric(strsplit(e,split = ",")[[1]][3])
NumeroRollos<-NumeroRollos+(i+j+k)*probaLP[contador]
}
NumeroRollos
## [1] 1.455402+0i
# Literal D: Tiempo primera pasada ----------------------------------------
ExpectedTime(cadena,which(estados=="1,1,1"),which(estados=="3,2,2"))
## [1] 44224.58