library(markovchain)
## Package: markovchain
## Version: 0.8.6
## Date: 2021-05-17
## BugReport: https://github.com/spedygiorgio/markovchain/issues
library(expm)
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.1.3
##
## Attaching package: 'expm'
## The following object is masked from 'package:Matrix':
##
## expm
# Variables ---------------------------------------------------------------
#X(t)= el n?mero de unidades en la zona de prensas de presi?n en el tiempo t
#Y(t)= el n?mero de unidades en la zona de recubrimiento en el tiempo t
#Z(t)= el n?mero de unidades en la zona de alisado en el tiempo t
#Espacio de estados
estados1<-c(0:3)
estados2<-c(0:2)
estados3<-c(0:2)
#Crear espacio de estados conjunto
estados<-c()
for (i in estados1) {
for (j in estados2) {
for (k in estados3) {
estados<-c(estados,paste(i,j,k,sep = ","))
}
}
}
estados2 <- apply(expand.grid(estados1,estados2,estados3),1,paste,collapse=",")
# Matriz de transiciones --------------------------------------------------
mat<-matrix(0,nrow = length(estados),ncol = length(estados),dimnames = list(estados,estados))
for (fila in estados) {
for (columna in estados) {
#Se realiza la division de los componentes, seran el origen
#i=as.numeric(strsplit(fila,",")[[1]][1]) #Opcion 2
i=as.numeric(unlist(strsplit(fila,","))[1])
j=as.numeric(unlist(strsplit(fila,","))[2])
k=as.numeric(unlist(strsplit(fila,","))[3])
#Se realiza la division de los componentes, seran el destino
ip=as.numeric(unlist(strsplit(columna,","))[1])
jp=as.numeric(unlist(strsplit(columna,","))[2])
kp=as.numeric(unlist(strsplit(columna,","))[3])
if (i<3 & ip==i+1 & jp==j & kp==k) { #Llegada a la etapa 3
mat[fila,columna]=0.05
}
if (i>0 & j<2 & ip==i-1 & jp==j+1 & kp==k) { #Sale de prensa a recubrimiento
mat[fila,columna]=0.104
}
if (j>0 & k<2 & ip==i & jp==j-1 & kp==k+1) { #Sale de recubrimiento a alisado
mat[fila,columna]=0.125*j
}
if (k>0 & ip==i & jp==j & kp==k-1) { #Sale de alisado
mat[fila,columna]=(11/12)*0.2
}
}
}
diag(mat)<--rowSums(mat)
rowSums(mat)
## 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
CMTC<-new(Class = "ctmc",states=estados,generator=mat)
# C Numero de rollos en estado estable --------------------------------------
is.CTMCirreducible(CMTC)
## [1] TRUE
probaLP<-as.numeric(steadyStates(CMTC))
names(probaLP)<-estados
Rollos<-0
for (a in estados) {
#Se realiza la division de los componentes
i=as.numeric(unlist(strsplit(a,","))[1])
j=as.numeric(unlist(strsplit(a,","))[2])
k=as.numeric(unlist(strsplit(a,","))[3])
Rollos<-Rollos+(i+j+k)*probaLP[a]
}
Rollos<-as.numeric(Rollos)
Rollos
## [1] 1.455402
# Tiempos de primera pasada CMTC ------------------------------------------
#Se hallar? este tiempo volviendo el estado a evaluar final (3,2,2) absorbente
alpha<-c(rep(0,length(estados)))
names(alpha)<-estados
alpha["1,1,1"]<-1
matAbsorbente<-mat
matAbsorbente['3,2,2',]<-0
estadosTransitoriosU<-estados[-which(estados=='3,2,2')]
matU<-matAbsorbente[-which(estados=='3,2,2'),-which(estados=='3,2,2')]
dimnames(matU)<-list(estadosTransitoriosU,estadosTransitoriosU)
#Calcular tiempo antes de a absorbcion como -u^-1
matTiemposAntesAbsorcion<- -solve(matU)
rowSums(matTiemposAntesAbsorcion)["1,1,1"] #Tiempo de primera pasada en segundos
## 1,1,1
## 44224.58
#Opcion 2
ExpectedTime(CMTC,which(estados=="1,1,1"),which(estados=='3,2,2'))
## [1] 44224.58