## ---------------------------------------------------------------------------------
#  Demo del paquete "queueing"
#  Ejemplos tomados de "Teoría de Colas" (Miguel Miranda)
#  Daniel Rigou (2019)
## ---------------------------------------------------------------------------------
# Modificar el directorio de trabajo empleado en la sentencia siguiente
setwd("~/datos/R/R_Archivos de Trabajo/Queueing")

## Carga el paquete "queueing". Requiere instalación previa.
library(queueing)

# La documentación del paquete "queueing" (Canadilla-Montoya) puede encontrarse en:
# https://www.rdocumentation.org/packages/queueing/versions/0.2.11

## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014") 


## Ejemplo 2.1  M/M/1

lambda<-5  # (clientes/hora)
mu<-6      # (clientes/hora)


Ejemplo_2_1 <- NewInput.MM1(lambda, mu, n=10)
CheckInput(Ejemplo_2_1)

model <- QueueingModel(Ejemplo_2_1)
print(summary(model), digits=4)
Report(model)

# Grafico las probabilidades de estado
s<-seq(from=0, to=10)
plot(s, model$Pn,type="p", main="Probabilidad de Estado", xlab="Cantidad de Clientes en el Sistema", 
     ylab="P(n)")
grid(col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)

## a) probabilidad de no esperar p(0) (notar que el primer elemento del vector Pn lleva
#     el subíndice 1

model$Pn[1]

## b) porcentaje de utilización del canal (RO)

model$RO

## c) probabilidad de que un cliente que llega al sistema encuentre que hay
## más de tres personas esperando en cola. Es igual a la probabilidad de que haya 
## cuatro o más en el sistema

(model$RO)^4
1-sum(model$Pn[1:4])

## d) Probabilidad de que encuentre que hay menos de dos personas en el sistema

sum(model$Pn[1:2])

## e) Longitud promedio del sistema

model$L

## f) Longitud promedio de la cola

model$Lq

## g) Ingreso esperado ($/hora)

lambda * 50

## h) Tiempo promedio de espera en cola (hora)

model$Wq

## i) Tiempo promedio de permanencia dentro del sistema (hora)

model$W

## j) Probabilidad de que un cliente permanezca en el sistema más de 5 minutos (0.83333 h)

1-model$FW(0.083333)

## k) Probabilidad de que un cliente espere en cola más de 3 minutos (0.05 h)

1-model$FWq(0.05)

## l) Velocidad óptima de atención (clientes/hora)

lambda+sqrt(lambda*20/15)


# Gráfico de la probabilidad de los tiempos de espera en el sistema y en la cola
# Función de distribución de W y Wq
# w: probabilidad de demorar menos de t en el sistema
# wq:probabilidad de demorar menos de t en la cola

gTitle <- "Función de distribución de W y Wq"
fw     <- model$FW
fwq    <- model$FWq
nn     <- 6
ty     <- "l"
ylab   <- "FW(t), FWq(t)"
xlab   <- "t"
cols   <- c("black", "red")
leg    <- c("FW(t)", "FWq(t)")

curve(fw, from=0, to=nn, type=ty, ylab=ylab, xlab=xlab, col=cols[1], main=gTitle, las=1)
curve(fwq, from=0, to=nn, type=ty, col=cols[2], add=T)
legend("bottomright", leg, lty=c(1, 1), col=cols)
grid(col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")


## Ejemplo 2.2  M/M/1

lambda=2.5            # (clientes/hora)
n<-c(1:8)             # n: cantidad de operarios asignados (parámetro)
Ts<-0.5*exp(-0.1*n)   # Tiempo de servicio en función de la cantidad de operarios (horas)
mu<-1/Ts              # (clientes/hora)
Cnop<-30              # Costo ómnibus no operativo ($/hora)
Sala<-15              # Salario de cada operario ($/hora)

L<-rep(NA, 8)         # inicializo vector L

for(i in 1:8){
  if(mu[i]>lambda){
          Ejemplo_2_2 <- NewInput.MM1(lambda, mu=mu[i])
          model <- QueueingModel(Ejemplo_2_2)
          L[i]<-model$L
  }
}


# Función objetivo
Z<-Cnop*L+Sala*n

tabla<-data.frame(n, Ts, mu, L, Z)
tabla

# Número óptimo de operarios
which.min(Z)

# Costo óptimo
min(Z, na.rm=TRUE)

# Gráfico del Costo Total como función de la cantidad de Operarios
plot(n, Z, type="b", xlab="Cantidad de Operarios", ylab="Costo Total ($/hora)", main="Costo Total en
     función de la Cantidad de Operarios", col="blue", las=1)
grid(col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")


## Ejemplo 2.3  M/M/1

lambda<-12      # (PC/mes)
mu1<-15         # (PC/mes)
mu2<-20         # (PC/mes)
Sala1<-4000     # Costo operario ($/mes)
Sala2<-6000     # Costo operario + Costo ayudante ($/mes)
Cnop<-1000      # Costo por computadora no operativa ($/mes)

Ejemplo_2_3_1 <- NewInput.MM1(lambda, mu=mu1, n=10)
CheckInput(Ejemplo_2_3_1)
model1 <- QueueingModel(Ejemplo_2_3_1)

Ejemplo_2_3_2 <- NewInput.MM1(lambda, mu=mu2, n=10)
CheckInput(Ejemplo_2_3_2)
model2 <- QueueingModel(Ejemplo_2_3_2)

print(summary(model1), digits=4)
print(summary(model2), digits=4)

Z1<-Cnop*model1$L+Sala1
Z1

Z2<-Cnop*model2$L+Sala2
Z2


# Grafico del Costo Total como función del Lucro Cesante
LC<-seq(from=500, to=1500, by=100)
CT1<-LC*model1$L+4000
CT2<-LC*model2$L+6000

cols   <- c("blue", "red")
yscale<-c(5500,10500)

plot(LC, CT1, type="b", ylim = yscale, ylab="Costo Total", xlab="Lucro Cesante", col=cols[1],
     main="Costo Total en función del Lucro Cesante", las=1)
par(new=TRUE)
plot(LC, CT2, type="b", ylim = yscale, pch=2, yaxt="n", xaxt="n", ylab=" ", xlab=" ", col=cols[2])

leg<-c("Actual", "Alternativa")
pchs<-c(1,2)
legend("bottomright", leg, lty=c(1, 1), col=cols, pch=pchs)
abline(v=1000, col="grey")
grid(col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")


# Ejemplo 2.4   M/M/1/10

lambda<-12   # (clientes/hora)
mu<-10       # (clientes/hora)
k<-10        # Capacidad del sistema (cola + canal de atención)   
Iserv<-30    # Ingreso por cada servicio ($/servicio)

Ejemplo_2_4 <- NewInput.MM1K(lambda, mu, k)
CheckInput(Ejemplo_2_4)

model <- QueueingModel(Ejemplo_2_4)
print(summary(model), digits=4)


# a) Probabilidad de no esperar
model$P[1]

# b) Porcentaje de ocupación del canal
1-model$P[1]

# c) Probabilidad de que un cliente que llega al sistema encuentre que hay tres o más
# personas
1-model$P[1]-model$P[2]-model$P[3]

# d) Probabilidad de encontrar menos de dos personas en el sistema
model$P[1]+model$P[2]

# e) Longitud promedio del sistema
model$L

# f) Longitud promedio de la cola
model$Lq

# g) Ingreso monetario esperado

#lambda_ingreso<-lambda * (1-model$Pn[11])
#lambda_ingreso

Ingreso<-Iserv*lambda * (1-model$Pn[11])
Ingreso

# h) Tiempo promedio de espera en cola
model$Wq

# i) Tiempo promedio de permanencia dentro del sistema
model$W


# Gráfico de la tasa de ingreso efectivo al sistema en función de la tasa de arribos lambda

lambda_aux<-seq(from=1, to=20, by=1)
lambda_ingreso<-rep(0, 20)

for(i in 1:20){
    Ejemplo_2_4 <- NewInput.MM1K(lambda_aux[i], mu, k)
    model_aux <- QueueingModel(Ejemplo_2_4)
    lambda_ingreso[i]<-lambda_aux[i]*(1-model_aux$Pn[11])
  
}

plot(lambda_aux, lambda_ingreso, type="l", xlab="Tasa de arribos (clientes/hora)", 
     ylab="Tasa de ingreso efectivo (clientes/hora)",
     main="Tasa de Ingreso Efectivo al Sistema en función de la Tasa de Arribos", 
     col="blue", las=1)

grid(col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)

## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014") 


# Ejemplo 2.5   M/M/1/4

lambda<-6     # (clientes/hora)
mu<-8         # (clientes/hora)
k<-4          # Capacidad del sistema (cola + canal de atención)   
Iserv<-100    # Ingreso por cada servicio ($/servicio)

Ejemplo_2_5 <- NewInput.MM1K(lambda, mu, k)
CheckInput(Ejemplo_2_5)

model <- QueueingModel(Ejemplo_2_5)
print(summary(model), digits=4)
Report(model)

n<-seq(from=0, to=4)
tabla<-data.frame(n, model$Pn)
tabla

# Verificamos que la tasa de ingreso sea igual a la de salida
Qing<-lambda*(1-model$Pn[5])
Qing

Qsal<-mu*(1-model$P[1])
Qsal

# Tasa de rechazos
Qrec<-lambda*model$Pn[5]
Qrec
# debe ser igual a:
lambda-Qsal

# Lucro cesante
LC<-Iserv*Qrec
LC

# Longitud del sistema
model$L

# Longitud de la cola
model$Lq

# Tiempo en cola
model$Wq

# Tiempo de permanencia en el sistema
model$W
# debe ser igual al tiempo en cola más el tiempo de servicio
model$Wq+1/mu



## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014") 


# Ejemplo 2.6   M/M/1/k

lambda<-2           # (clientes/día)
mu<-2               # (clientes/día)
CostoCola<-10       # Costo de cada posición de cola por día
IServ<-100          # Ingreso por cada servicio ($/servicio)
  
Qrec<-rep(0,10)     # Inicializo vector para clientes rechazados
Z<-rep(0,10)        # Iniciliazo vector para funcional

for(k in 1:10){
  Ejemplo_2_6 <- NewInput.MM1K(lambda, mu, k)
  model_aux <- QueueingModel(Ejemplo_2_6)
  Qrec[k]=lambda*model_aux$Pn[k]                  # Cantidad de clientes rechazados
  Z[k]=CostoCola*(k-1)+IServ*Qrec[k]              # Funcional: Costo de posiciones en cola +
                                                  # costo de oportunidad de ingresos perdidos
}

# Tamaño óptimo del sistema
which.min(Z)

k<-seq(from=1,to=10)
tabla<-data.frame(k,Z)
tabla                                             # Tamaño del Sistema y Funcional

plot(k, Z, type="b", col="blue", xlab="Tamaño del Sistema", ylab="Costo Total ($/día)",
     main="Costo de Oportunidad de los Clientes Rechazados + Costo de lugares en cola", las=1)
grid(col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")


# Ejemplo 3.1   M/M/4/

lambda<-12          # (clientes/hora)
mu<-6               # (clientes/hora)
c<-4                # Cantidad de canales de servicio
IServ<-50           # Ingreso por cada servicio ($/servicio)

Ejemplo_3_1 <- NewInput.MMC(lambda, mu, c, n=10)
CheckInput.i_MMC(Ejemplo_3_1)
model <- QueueingModel(Ejemplo_3_1)

print(summary(model), digits=4)
Report(model)

# a) Cantidad mínima de canales
floor(lambda/mu)+1

# b) Probabilidad de no esperar. P(n<=3)
sum(model$Pn[1:4])

# c) Porcentaje de utilización del canal
model$RO          # en este paquete, RO=lambda/(mu.c) 

# d) Probabilidad de que haya menos de dos personas en el sistema: P(0)+P(1)
model$Pn[1]+model$Pn[2]

# e) Número promedio de clientes esperando ser atendidos
model$Lq

# f) Número promedio de clientes en el sistema
model$L

# g) Ingreso de caja esperado ($/hora)
IServ * Throughput(model)

# h) Tiempo promedio de permanencia de un cliente dentro del sistema (horas)
model$W


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")


# Ejemplo 3.2   M/M/C/

lambda<-5           # (camiones/hora)
mu<-2               # (camiones/hora)
Cnop<-50            # Costo por camión no operativo ($/hora)
Ccan<-30            # Costo operativo por canal ($/hora)

# Rango de alternativas para c (cantidad de canales)
rango<-4
c<-seq(from=floor(lambda/mu)+1, to=floor(lambda/mu)+rango)
Z<-rep(0,rango)     # Iniciliazo vector para funcional

# Resolución del modelo para diferentes valores de c (cantidad de canales)
for(i in 1:rango){
    Ejemplo_3_2 <- NewInput.MMC(lambda, mu, c[i])
    model_aux <- QueueingModel(Ejemplo_3_2)
    Z[i]<- Cnop * model_aux$Lq + Ccan * c[i]
}

# Número óptimo de canales de atención
copt<-c[which.min(Z)]
copt

# Costo total óptimo
Zopt<-Z[which.min(Z)]
Zopt

# Calculo el modelo para la cantidad óptima de canales
Ejemplo_3_2 <- NewInput.MMC(lambda, mu, c=copt, n=10)
model <- QueueingModel(Ejemplo_3_2)
print(summary(model), digits=4)
Report(model)

# Gráfico del Costo Total en función de la cantidad de canales de atención
plot(c, Z, type="b", col="blue", xlab="Cantidad de canales de atención", ylab="Costo Total ($/hora)",
     main="Costo Total en función de la cantidad de canales de atención", 
     xaxp=c(c[1],c[rango],c[rango]-c[1]), las=1)
grid(col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")


# Ejemplo 3.3   M/M/3/ versus 3 sistemas M/M/1

lambda1<-10         # (clientes/hora)
lambda2<-30         # (clientes/hora)
mu<-12              # (clientes/hora)


# 3 Sistemas M/M/1 en paralelo
Ejemplo_3_3 <- NewInput.MM1(lambda1, mu)
model_1 <- QueueingModel(Ejemplo_3_3)
print(summary(model_1), digits=4)

# Tiempo en cola
Wq(model_1)                          # es equivalente a: model_1$Wq
# Tiempo en el sistema
W(model_1)


# Sistema M/M/3
Ejemplo_3_3 <- NewInput.MMC(lambda2, mu, c=3)
model_2 <- QueueingModel(Ejemplo_3_3)
print(summary(model_2), digits=4)

# Tiempo en cola
Wq(model_2)
# Tiempo en el sistema
W(model_2)


# ........................
# Análisis complementario: M/M/1 con tasa de servicio tres veces más rápida

lambda3<-lambda2
mu3<-3*mu

Ejemplo_3_3 <- NewInput.MM1(lambda3, mu3)
model_3 <- QueueingModel(Ejemplo_3_3)
print(summary(model_3), digits=4)

# Tiempo en cola
Wq(model_3)
# Tiempo en el sistema
W(model_3)


# ........................
# Resumen
W<-c(W(model_1),W(model_2),W(model_3))
Wq<-c(Wq(model_1),Wq(model_2),Wq(model_3))
tabla<-data.frame(W,Wq)
row.names(tabla)<-c("3xMM1", "MM3", "MM1 3xmu")
colnames(tabla)<-c("Tiempo en el Sistema", "Tiempo en Cola")

tabla

# Plot de las tres alternativas
plot(tabla, ylim=c(0.1,0.5), xlim=c(0.1,0.6), main="Comparación del los tiempos en cola y 
     en el sistema para los tres modelos", las=1)
text(tabla, labels = row.names(tabla), pos=3)
grid(col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")


# Ejemplo 3.4   M/M/2/4/

lambda<-12         # (clientes/hora)
mu<-6              # (clientes/hora)
c<-2               # cantidad de canales de atención
k<-4               # capacidad del sistema

# Resolución del modelo
Ejemplo_3_4 <- NewInput.MMCK(lambda, mu, c, k)
model <- QueueingModel(Ejemplo_3_4)
CheckInput.i_MMCK(Ejemplo_3_4)
print(summary(model), digits=4)
Report(model)

# Cantidad de clientes rechazados, e ingreso de clientes en función de la velocidad de atención
Qrec<-rep(0,15)
Qing<-rep(0,15)
for(i in 1:15){
  Ejemplo_3_4 <- NewInput.MMCK(lambda, i, c, k)
  model_aux <- QueueingModel(Ejemplo_3_4)
  Qrec[i]<-lambda-Throughput(model_aux)
  Qing[i]<-Throughput(model_aux)
}

# Gráfico de la cantidad de clientes rechazados en función de la velocidad de atención
plot(Qrec, type="b", col="blue", xlab="Velocidad de atención (clientes/hora)", 
     ylab="Clientes rechazados (clientes/hora)",
     main="Clientes rechazados en función de la velocidad de atención", las=1)
grid(col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)
abline(v=mu, col="grey")

# Gráfico del ingreso de clientes en función de la velocidad de atención
plot(Qing, type="b", col="blue", xlab="Velocidad de atención (clientes/hora)", 
     ylab="Ingreso de clientes (clientes/hora)",
     main="Ingreso de clientes en función de la velocidad de atención", las=1)
grid(col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)
abline(v=mu, col="grey")


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")


# Ejemplo 4.1   M/M/1/(5) (Población finita)

lambda<-2          # (clientes/hora)
mu<-3              # (clientes/hora)
k<-5               # capacidad del sistema


# Resolución del modelo
Ejemplo_4_1 <- NewInput.MM1KK(lambda, mu, k)
model <- QueueingModel(Ejemplo_4_1)
CheckInput.i_MM1KK(Ejemplo_4_1)
print(summary(model), digits=4)
Report(model)


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")

# Ejemplo 4.3   M/M/3/(10) (Población finita)

lambda<-0.25       # (pedidos/hora)
mu<-2              # (servicios/hora)
c<-3               # cantidad de canales de atención
k<-10              # capacidad del sistema
m<-10              # tamaño de la población

# Resolución del modelo
Ejemplo_4_3 <- NewInput.MMCKM(lambda, mu, c, k, m, method=0)
model <- QueueingModel(Ejemplo_4_3)
CheckInput.i_MMCKM(Ejemplo_4_3)
print(summary(model), digits=4)
Report(model)

D<-sum(model$Pn[4:10])
D


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")

# Ejemplo 4.4   M/M/2/(5) (Población finita)

lambda<-0.5        # (pedidos/hora)
mu<-1              # (servicios/hora)
c<-2               # cantidad de canales de atención
k<-5               # capacidad del sistema
m<-5               # tamaño de la población

# Resolución del modelo
Ejemplo_4_4 <- NewInput.MMCKM(lambda, mu, c, k, m)
model <- QueueingModel(Ejemplo_4_4)
CheckInput.i_MMCKM(Ejemplo_4_4)
print(summary(model), digits=4)
Report(model)


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")

# Ejemplo 5.4  Red abierta (Open Jackson Network)

# Matriz de transiciones entre nodos
m<-c(0,0.3,0.3,0.4,0,0,0,0.5,0,0,0,1,0,0,0,0)
prob<-matrix(data=m, nrow=4, ncol=4, byrow=TRUE)
prob

# Definición de las características de cada nodo
# Nota 1: los valores de lambda son los arribos desde el exterior de la red al nodo
# si no hay arribos desde el exterior de la red, entonces lambda=0 para ese nodo
n1 <- NewInput.MM1(lambda=10, mu=15)
n2 <- NewInput.MM1(lambda=15, mu=20)
n3 <- NewInput.MM1(lambda=0, mu=4)
n4 <- NewInput.MM1(lambda=0, mu=20)

# Resolución del modelo
ojn1 <- NewInput.OJN(prob, n1, n2, n3, n4)
model <- QueueingModel(ojn1)

print(summary(model), digits=4)
Report(model)


# Nota 2: los tiempos en cada nodo que brinda el reporte son para el total de
# los clientes (los que pasan por el nodo y los que no pasan).
# La suma de estos tiempos es igual al tiempo promedio en el sistema: W=Wk[1]+Wk[2]+Wk[3]+Wk[4]
# Los tiempos en el nodo para los clientes que pasan por el nodo pueden 
# calcularse por Little

WN<-model$Lk/model$Throughputk
WN


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")

# Ejemplo 5.7  Red abierta con Reciclado

# Matriz de transiciones entre nodos
m<-c(0,0.2,0,0.7,0,0.3,0,1,0)
prob<-matrix(data=m, nrow=3, ncol=3, byrow=TRUE)
prob

# Definición de las características de cada nodo
# Nota 1: los valores de lambda son los arribos desde el exterior de la red al nodo
# si no hay arribos desde el exterior de la red, entonces lambda=0 para ese nodo
n1 <- NewInput.MM1(lambda=10, mu=14)
n2 <- NewInput.MM1(lambda=0, mu=4)
n3 <- NewInput.MM1(lambda=0, mu=2)

# Resolución del modelo
ojn1 <- NewInput.OJN(prob, n1, n2, n3)
model <- QueueingModel(ojn1)

print(summary(model), digits=4)
Report(model)

WN<-model$Lk/model$Throughputk
WN

# Probabilidad de que el sistema esté vacio
Pvacio<-model$Pn[1]*model$Pn[2]*model$Pn[3]
Pvacio

## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")

# Ejemplo 5.8  Red cerrada

m<-c(0,1,1,0)
prob<-matrix(data=m, nrow=2, ncol=2, byrow=TRUE)
prob

n<-9                                   # número de clientes en el sistema
z<-0                                   # think time of the client ?
operational<-FALSE                     # los valores de la matriz prob son probabilidades

n1 <- NewInput.MM1(lambda=0, mu=10)    # subsistema 1
n2 <- NewInput.MM1(lambda=0, mu=9.6)   # subsistema 2

# Resolución del modelo
cjn1 <- NewInput.CJN(prob, n, z, operational, 0, 0.001, n1, n2)
model<-QueueingModel(cjn1)

print(summary(model), digits=4)
Report(model)


## ---------------------------------------------------------------------------------
rm(list=ls())
dev.off()
cat("\014")

# Ejemplo 5.9  Red cerrada (Comparar los resultados con el método aproximado del libro)

m<-c(0,0.3,0.7,1,0,0,1,0,0)
prob<-matrix(data=m, nrow=3, ncol=3, byrow=TRUE)
prob

n<-20                                  # número de clientes en el sistema
z<-0                                   # think time of the client ?
operational<-FALSE                     # los valores de la matriz prob son probabilidades

n1 <- NewInput.MM1(lambda=0, mu=10)    # subsistema 1
n2 <- NewInput.MM1(lambda=0, mu=4)     # subsistema 2
n3 <- NewInput.MM1(lambda=0, mu=7)     # subsistema 3

# Resolución del modelo
cjn1 <- NewInput.CJN(prob, n, z, operational, 0, 0.001, n1, n2, n3)
model<-QueueingModel(cjn1)

print(summary(model), digits=4)
Report(model)

# ...............
# Comparación de los resultados con el modelo aproximado para distintos
# valores del número de clientes en el sistema

# Análisis exacto con el modelo
cjn1_20<- NewInput.CJN(prob, n=20, z, operational, 0, 0.001, n1, n2, n3)
cjn1_100<- NewInput.CJN(prob, n=100, z, operational, 0, 0.001, n1, n2, n3)
cjn1_200<- NewInput.CJN(prob, n=200, z, operational, 0, 0.001, n1, n2, n3)
cjn1_400<- NewInput.CJN(prob, n=400, z, operational, 0, 0.001, n1, n2, n3)
model_20<-QueueingModel(cjn1_20)
model_100<-QueueingModel(cjn1_100)
model_200<-QueueingModel(cjn1_200)
model_400<-QueueingModel(cjn1_400)

lambda1<-c(model_20$Throughputk[1],model_100$Throughputk[1],
           model_200$Throughputk[1],model_400$Throughputk[1])

# ..............
# Análisis aproximado despejando lambda1 de L1+L2+L3
# n=20
func <- function(x) {
  x/(10-x) + 0.3*x/(4-0.3*x) + 0.7*x/(7-0.7*x) -20
}
lambda1_20<-uniroot(func, c(0,10))

#n=100
func <- function(x) {
  x/(10-x) + 0.3*x/(4-0.3*x) + 0.7*x/(7-0.7*x) -100
}
lambda1_100<-uniroot(func, c(0,10))

# n=200
func <- function(x) {
  x/(10-x) + 0.3*x/(4-0.3*x) + 0.7*x/(7-0.7*x) -200
}
lambda1_200<-uniroot(func, c(0,10))

# n=400
func <- function(x) {
  x/(10-x) + 0.3*x/(4-0.3*x) + 0.7*x/(7-0.7*x) -400
}
lambda1_400<-uniroot(func, c(0,10))

lambda1_aprox<-c(lambda1_20$root, lambda1_100$root, lambda1_200$root,lambda1_400$root )

error<-(lambda1-lambda1_aprox)/lambda1*100
error

# Gráfico del error porcentual para lambda1
x<-c(20,100,200,400)
plot(x,error, xlab="cantidad de clientes en el sistema",
     ylab="error porcentual", main="Error porcentual para lambda 1 con el modelo aproximado",
     type="b", col="blue", las=1)
grid(col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)


## ---------------------------------------------------------------------------------