introducción:
Estamos modelando una estación de reproceso. Las piezas llegan con un periodo fijo de 20 minutos. Cada pieza puede presentar hasta 3 defectos. El número de defectos por pieza sigue una distribución Binomial con media 2.4. Cada defecto requiere un tiempo de reparación aleatorio, modelado por una Exponencial con tasa 0.2. El objetivo es estimar el instante en que se termina de procesar la pieza número 200, simulando 10 réplicas independientes.
“Por cada réplica generamos 200 piezas. Para cada pieza primero muestreamos el número de defectos K ~ Binomial(3, 0.8). Si K>0, generamos K tiempos exponenciales (independientes) y sumamos para obtener el tiempo de servicio de la pieza. Las piezas llegan en tiempos 0, 20, 40, … y el inicio de servicio para la pieza i es el máximo entre su tiempo de llegada y el fin de servicio de la pieza i-1. Guardamos el instante de finalización de la pieza 200. Repetimos 10 veces y analizamos la distribución de esos tiempos.”
# Parámetros
n_piezas <- 200
n_replicas <- 10
n_binom <- 3
p_binom <- 2.4 / n_binom # 0.8
rate_defecto <- 0.2 # tasa exponencial (por minuto)
intervalo_llegada <- 20 # minutos entre piezas
# Simulación
sim_results <- tibble(replica = integer(), tiempo_total = double())
for(rep in 1:n_replicas){
# número de defectos por pieza
defects <- rbinom(n_piezas, size = n_binom, prob = p_binom)
# tiempos de servicio por pieza (suma de exponenciales si hay defectos)
service_times <- numeric(n_piezas)
for(i in 1:n_piezas){
k <- defects[i]
if(k > 0) service_times[i] <- sum(rexp(k, rate = rate_defecto)) else service_times[i] <- 0
}
# tiempos de llegada determinísticos
arrival_times <- (0:(n_piezas-1)) * intervalo_llegada
start_service <- numeric(n_piezas)
end_service <- numeric(n_piezas)
for(i in 1:n_piezas){
if(i == 1) start_service[i] <- arrival_times[i] else start_service[i] <- max(arrival_times[i], end_service[i-1])
end_service[i] <- start_service[i] + service_times[i]
}
tiempo_total <- end_service[n_piezas]
sim_results <- sim_results %>% add_row(replica = rep, tiempo_total = tiempo_total)
}
sim_results
## # A tibble: 10 × 2
## replica tiempo_total
## <int> <dbl>
## 1 1 3985.
## 2 2 4000.
## 3 3 3989.
## 4 4 3980.
## 5 5 4001.
## 6 6 3996.
## 7 7 4017.
## 8 8 3995.
## 9 9 3988.
## 10 10 3985.
# Estadísticos
summary_tab <- sim_results %>% summarise(media = mean(tiempo_total), sd = sd(tiempo_total))
summary_tab
## # A tibble: 1 × 2
## media sd
## <dbl> <dbl>
## 1 3994. 10.8
Interpretación:
“Tras ejecutar la simulación con 10 réplicas, la media del tiempo total para procesar las 200 piezas fue 3993.84 minutos y la desviación estándar fue 10.78 minutos. Esto equivale a aproximadamente 66.56 horas en promedio. Estas cifras implican que, bajo los supuestos hechos (Binomial(3,0.8) y tiempos exponenciales de reparación), el sistema completa 200 piezas en alrededor de 66.56 horas. La pequeña desviación estándar indica que el tiempo total es relativamente consistente entre réplicas.”
library(ggplot2)
p1 <- ggplot(sim_results, aes(x = tiempo_total)) +
geom_histogram(bins = 10, fill = "steelblue", alpha = 0.7) +
labs(title = "Distribución del tiempo total (200 piezas, 10 réplicas)",
x = "Tiempo total (minutos)", y = "Frecuencia") +
geom_vline(aes(xintercept = mean(tiempo_total)), color = "red", linetype = "dashed") +
annotate("text", x = mean(sim_results$tiempo_total), y = Inf, label = paste0("Media = ", round(mean(sim_results$tiempo_total),2)),
vjust = 2, hjust = 1.1, color = "red", size = 3.5)
p1
Puntos clave:
introducción:
En este punto simulamos la operación de un camión que se carga,
viaja para descargar y regresa. Los parámetros son: carga 30 ± 10 min
(Normal), descarga 20 ± 5 min (Normal) y traslado exponencial con media
40 min (rate = 1/40). Solo hay un muelle de carga (un camión a la vez).
Primero simulamos 10 horas con 5 réplicas. Luego calculamos un intervalo
de confianza para el número de viajes en 24 horas usando 2000
réplicas.
set.seed(4040)
horizon_10h <- 10 * 60 # minutos
nrep_10h <- 5
rnorm_trunc <- function(n, mean, sd){
x <- rnorm(n, mean = mean, sd = sd)
x[x < 1] <- 1
x
}
simulate_once_10h <- function(horizon){
t <- 0
viajes <- 0
while(TRUE){
# carga
t <- t + rnorm_trunc(1, 30, 10)
if(t > horizon) break
# traslado ida
t <- t + rexp(1, rate = 1/40)
if(t > horizon) break
# descarga
t <- t + rnorm_trunc(1, 20, 5)
if(t > horizon) break
viajes <- viajes + 1
# regreso
t <- t + rexp(1, rate = 1/40)
if(t > horizon) break
}
viajes
}
viajes_10h <- sapply(1:nrep_10h, function(i) simulate_once_10h(horizon_10h))
viajes_10h
## [1] 3 2 5 3 5
interpretación:
“En la simulación de 10 horas con 5 réplicas obtuvimos los siguientes números de viajes por réplica: 3, 2, 5, 3, 5. Esto ilustra la variabilidad en un periodo acotado; la media es 3.6 viajes en 10 horas.”
qplot(x = factor(1:length(viajes_10h)), y = viajes_10h, geom = "col") +
labs(title = "Viajes completados en 10 horas (5 réplicas)", x = "Réplica", y = "Viajes")
set.seed(2023)
horizon_24h <- 24 * 60
nrep_ci <- 2000
simulate_day <- function(horizon){
t <- 0
viajes <- 0
while(TRUE){
t <- t + rnorm_trunc(1, 30, 10)
if(t > horizon) break
t <- t + rexp(1, rate = 1/40)
if(t > horizon) break
t <- t + rnorm_trunc(1, 20, 5)
if(t > horizon) break
viajes <- viajes + 1
t <- t + rexp(1, rate = 1/40)
if(t > horizon) break
}
viajes
}
# ejecutar replicaciones
viajes_24h <- replicate(nrep_ci, simulate_day(horizon_24h))
# estadísticas
mean_24h <- mean(viajes_24h)
sd_24h <- sd(viajes_24h)
ci_emp <- quantile(viajes_24h, probs = c(0.025, 0.975))
list(mean = mean_24h, sd = sd_24h, ci_95_emp = ci_emp)
## $mean
## [1] 11.014
##
## $sd
## [1] 1.448742
##
## $ci_95_emp
## 2.5% 97.5%
## 8 14
interpretación:
“La simulación sugiere que un camión realiza en promedio 11.01
viajes por día, con desviación estándar 1.45. El intervalo empírico del
95% (percentiles 2.5% y 97.5%) es
[r ci_emp[1], r ci_emp[2]] viajes por día. Notamos que,
debido a la variabilidad de cargas, traslados y descargas, existe
incertidumbre y la probabilidad de no alcanzar exactamente 10 viajes en
un día no es nula.”
hist(viajes_24h, breaks = 20, main = "Histograma: viajes en 24 horas (simulación)",
xlab = "Viajes por camión en 24h")
abline(v = mean_24h, col = "red", lwd = 2)
text(mean_24h + 1, max(hist(viajes_24h, plot=FALSE)$counts)*0.9, paste0("Media=", round(mean_24h,2)))
plot(density(viajes_24h), main = "Densidad: viajes por día", xlab = "Viajes en 24h")
abline(v = mean_24h, col = "red")
“Si la empresa requiere al menos 10 entregas por día, debemos mirar la media y el riesgo: mi simulación indica una media de 11.01 viajes y un intervalo empírico de 8 a 14. Para garantizar con alta probabilidad ≥10 entregas, recomiendo:”
Efecto de reducir el tiempo medio de carga.
# Este bloque está desactivado por defecto (eval=FALSE) para no alargar el KNIT.
# Si desea correrlo en vivo, cambie eval=TRUE.
library(purrr)
horizon_day <- 24*60
means <- c(30, 24, 20)
res <- tibble(mean_load = means) %>%
mutate(mean_viajes = map_dbl(mean_load, function(mu){
set.seed(100 + round(mu))
nrep <- 300
mean(replicate(nrep, {
t <- 0; viajes <- 0
while(TRUE){
t <- t + max(rnorm(1, mu, 10), 1)
if(t > horizon_day) break
t <- t + rexp(1, rate = 1/40)
if(t > horizon_day) break
t <- t + max(rnorm(1, 20, 5), 1)
if(t > horizon_day) break
viajes <- viajes + 1
t <- t + rexp(1, rate = 1/40)
if(t > horizon_day) break
}
viajes
}))
}))
res