Partiendo de los datos de movimientos de BiciMAD realizamos dos modelos, uno para predicir la demanda por estación y hora del día y otro para predecir la demanda por estación y día en los próximos meses.
En esta ocasión, no hemos considerado los datos de marzo a junio de 2020 ya que la crisis del Coronavirus ha cambiado mucho el comportamiento normal de los usuarios de BiciMAD. En la presentación de la serie hemos visto como durante la segunda quincena de marzo y hasta finales de abril, no ha habido uso de BiciMAD. En mayo y junio todavía existían bastantes restricciones horarias en cuanto a la movilidad de la población. Todas estas circunstancias nos han llevado a prescindir de ese periodo temporal.
# Librerías utilizadas durante todo el proyecto
library(ndjson)
library(stringr)
library(chron)
library(dplyr)
library(rjson)
library(geosphere)
library(caret)
library(ranger)
library(data.table)
library(ggplot2)
library(lubridate)
library(Hmisc)
library(utils)# Directorio de trabajo
setwd("C:/Users/rodsp/Desktop/DataScience/DATASETS/TFM/MOVIMIENTOS TOTAL")
# Leer datos
movimientos <- as.data.frame(fread("movimientos2019-jun2020.csv"))
summary(movimientos)## _id.$oid user_day_code user_type ageRange
## Length:5062186 Length:5062186 Min. :1.000 Min. :1.000
## Class :character Class :character 1st Qu.:1.000 1st Qu.:4.000
## Mode :character Mode :character Median :1.000 Median :4.000
## Mean :1.148 Mean :4.239
## 3rd Qu.:1.000 3rd Qu.:5.000
## Max. :3.000 Max. :6.000
## idplug_base idunplug_base idplug_station idunplug_station
## Min. : 1.00 Min. : 1.00 Min. : 1.00 Min. : 1.0
## 1st Qu.: 5.00 1st Qu.: 5.00 1st Qu.: 51.00 1st Qu.: 51.0
## Median :12.00 Median :12.00 Median : 96.00 Median : 96.0
## Mean :12.13 Mean :12.11 Mean : 99.84 Mean : 99.9
## 3rd Qu.:19.00 3rd Qu.:19.00 3rd Qu.:149.00 3rd Qu.:149.0
## Max. :30.00 Max. :30.00 Max. :219.00 Max. :219.0
## unplug_hourTime travel_time zip_code fe_year
## Length:5062186 Min. : 0 Length:5062186 Min. :2019
## Class :character 1st Qu.: 462 Class :character 1st Qu.:2019
## Mode :character Median : 700 Mode :character Median :2019
## Mean : 1491 Mean :2019
## 3rd Qu.: 1071 3rd Qu.:2019
## Max. :19272889 Max. :2020
## fe_month fe_day fe_wday fe_hour
## Min. : 1.00 Min. : 1.00 Min. :1.000 Min. : 0.00
## 1st Qu.: 3.00 1st Qu.: 8.00 1st Qu.:3.000 1st Qu.: 9.00
## Median : 5.00 Median :14.00 Median :4.000 Median :14.00
## Mean : 5.83 Mean :15.08 Mean :4.066 Mean :13.67
## 3rd Qu.: 9.00 3rd Qu.:23.00 3rd Qu.:6.000 3rd Qu.:18.00
## Max. :12.00 Max. :31.00 Max. :7.000 Max. :23.00
## fe_season
## Length:5062186
## Class :character
## Mode :character
##
##
##
# Nos quedamos solo con los datos hasta febrero 2020
# Conviertimos a fecha_hms
movimientos$unplug_hourTime <- ymd_hms(movimientos$unplug_hourTime)
movimientos$fecha <- as.Date(movimientos$unplug_hourTime, format="%d/%m/%Y")
movimientos <- subset(movimientos, movimientos$fecha < '2020-03-01')
# Eliminamos la variable unplug_hourTime ya que da la misma info que fecha
movimientos$unplug_hourTime <- NULLA continuación, exploramos y limpiamos las variables de este dataset creado con los datos de movimientos de BiciMAD desde enero de 2019 a febrero de 2020.
Las variables de id no son relevantes ya que son códigos generamos por BiciMAD para cada viaje, pero no nos proporciona información de los usuarios ni de las estaciones, por lo que las eliminamos.
Para nuestra predicción hemos optado por eliminar los usuarios tipo 3 ya que no son viajes propiamente dichos sino movimientos que realizan los empleados de BiciMAD para el mantimiento de las bicis.
Estas variables relativas al usuario tampoco las consideramos relevantes para el modelo, pues lo que estamos prediciendo son número de viajes por estación a cada hora del día sin saber previamente la información de los usuarios para cada uno de esos viajes.
Recategorizamos esta variable como factor y cambiamos por su frecuencia.
Consideramos que los viajes de menos de 60 sg son viajes erroneos, por ejemplo, porque la bici no funciona correctamente y se vuelve a dejar en la misma estación y los viajes superiores a 7200 sg, por ejemplo, porque no se ha contabilizado bien el tiempo (suele pasar y tenemos experiencia en ello). Se crea un nuevo dataframe llamado movimientos_sint3_time.
Por lo tanto, ahora tenemos viajes con los siguientes valores:
Después de filtrar la duración de viajes con la que trabajaremos, eliminamos la variable travel time de nuestro data frame, pues es información de usuarios.
# Eliminación de viajes cortos y largos
movimientos_sint3_time <- subset(movimientos_sint3, movimientos_sint3$travel_time > 60)
movimientos_sint3_time <- subset(movimientos_sint3_time, movimientos_sint3_time$travel_time < 7200)
# Visualización
summary(movimientos_sint3_time$travel_time)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 61.0 454.0 667.0 841.5 971.0 7199.0
## [1] "user_type" "idplug_base" "idplug_station" "idunplug_station"
## [5] "travel_time" "fe_year" "fe_month" "fe_day"
## [9] "fe_wday" "fe_hour" "fecha" "fe2_season"
El servicio BiciMAD funciona en el centro de la ciudad de Madrid. Dado que es un servicio al aire libre, consideramos relevante estudiar si la metereologia del día influye en su uso (parece que sí según los datos de meses y estaciones vistos en el EDA).
Descargamos estos datos de la Agencia Estatal de Metereología (Aemet) que también cuenta con una base de datos pública a través de su API. Hemos descargado los datos de clima de los meses que nos interesan para la estación de Retiro ya que es la más centrica y acorde al area de influencia de BiciMAD.
Para la lectura de datos utilizamos la función read.csv2 ya que los decimales los configura con una coma en lugar de punto, con esta función asume que esas comas son puntos y trata la variable como numérica (sino lo hace como character).
# Directorio de trabajo
setwd("C:/Users/rodsp/Desktop/DataScience/DATASETS/TFM/MOVIMIENTOS TOTAL")
# Lectura de datos
clima <- read.csv2("aemet_ene2019-feb2020.csv")
head(clima)## 'data.frame': 425 obs. of 19 variables:
## $ fecha : chr "01/01/2019" "02/01/2019" "03/01/2019" "04/01/2019" ...
## $ indicativo : int 3195 3195 3195 3195 3195 3195 3195 3195 3195 3195 ...
## $ nombre : chr "MADRID, RETIRO" "MADRID, RETIRO" "MADRID, RETIRO" "MADRID, RETIRO" ...
## $ provincia : chr "MADRID" "MADRID" "MADRID" "MADRID" ...
## $ altitud : int 667 667 667 667 667 667 667 667 667 667 ...
## $ tmed : num 6.6 6 7.3 5.9 5.4 6.6 7.1 7.1 7.6 4.7 ...
## $ prec : chr "0" "0" "0" "0" ...
## $ tmin : num 1.4 0.4 2.8 0.4 -0.1 -0.5 0.5 1.2 2 0.6 ...
## $ horatmin : chr "7:30" "6:30" "7:00" "8:05" ...
## $ tmax : num 11.7 11.7 11.8 11.4 10.8 13.8 13.7 13 13.2 8.8 ...
## $ horatmax : chr "15:10" "15:15" "15:20" "15:10" ...
## $ dir : int 6 6 5 6 9 2 8 34 3 3 ...
## $ velmedia : num 0.3 0.6 0.6 0.3 0.6 1.1 0.6 0.8 1.4 2.8 ...
## $ racha : num 5.3 7.8 7.8 3.6 3.1 3.3 5 3.6 9.7 10 ...
## $ horaracha : chr "0:00" "23:20" "0:10" "5:40" ...
## $ presMax : num 954 952 951 953 955 ...
## $ horaPresMax: chr "9" "0" "10" "Varias" ...
## $ presMin : num 952 948 949 950 952 ...
## $ horaPresMin: chr "5" "16" "5" "5" ...
Primero, quitamos las variables no relevantes para el ejercicio: el indicativo, el nombre de la estación, la provincia y la altitud ya que son siempre iguales para la estación Retiro.
Dirección y racha del viento tienen valores nulos que imputamos por aleatorio.
Para poder unir esta dataframe al anterior de movimiento lo hacemos con la variable fecha, para ello, debe tener el mismo formato en ambos dataframes.
# Quitamos ciudad, provincia etc y las de hora
clima_final <- clima[ , -c(2:5, 9, 11, 15, 16, 17, 19)]
# Imputamos dir y racha
clima_final$dir <- impute(clima_final$dir, "random")
clima_final$racha <- impute(clima_final$racha, "random")
# Cambiamos a fecha la fecha del clima
clima_final$fecha <- as.Date(clima_final$fecha, format="%d/%m/%Y")Con las precipitaciones creamos una nueva variable binaria: - 1: llovió - 0: no llovió
La dirección del viento la modificamos a intervalos abiertos por la derecha según Rosa del viento (según este método se podrían hacer más o menos intervalos, hemos cogido el intermedio): - Norte: 0-45º - Noreste: 45-90º - Este: 90-135º - Sureste: 135-180º - Sur: 180-225º - Suroeste: 225-270º - Oeste: 270-335º - Noroeste: 335-360º En nuestro caso solo hay valores hasta el 99 por lo que solo habrá viento Norte, Noreste y Este que hemos llamado respectivamente 1, 2 y 3
La velocidad del viento la modificamos en tramos (está medido en m/s): - Flojo: 0 - 5.4 que serán los valores 1 - Fresco: 5.5 - 17.1 que serán los valores 2 - Temporal > 17.2 que serán los valores 3 (aunque no habrá ninguno)
Estas ideas sobre la dirección y velocidad del viento han sido extraidas de esta página: https://es.windfinder.com/wind/windspeed.htm
clima_final$fe_prec_bin <- factor(replace(clima_final$prec, which(clima_final$prec != 0), 1))
clima_final$fe_dir <- cut(clima_final$dir, breaks = c(-1, 45, 90, 135), labels = c("1", "2", "3")) # no ponemos más intervalos porque hemos visto en el summary que el máximo es 99
clima_final$fe_velmedia <- cut(clima_final$velmedia, breaks = c(-1, 5.4, 17.1, 50), labels = c("1", "2", "3"))Estas variables las unimos al dataframe anterior de movimientos a través de la variable fecha, creando un nuevo dataframe que llamaremos mov_con_clima.
#Unión de ambos df
mov_con_clima <-merge(movimientos_sint3_time, clima_final, by = "fecha")
head(mov_con_clima)Cargamos los datos de las estaciones que también tratamos anteriormente, en este caso, utilizaremos las variables de localización de la estación, su id para poder unirlo al otro dataframe y el número de bases con las que cuenta.
Como estamos prediciendo el número de viajes desde la estación origen, consideramos estos valores para ella (idunplug_station).
Este conjunto de datos nos permite añadir más información a nuestro data frame, especialmente localizar en qué parte de la ciudad se situa cada estación.
# Directorio de trabajo
setwd("C:/Users/rodsp/Desktop/DataScience/DATASETS/TFM/MOVIMIENTOS TOTAL")
# Leer datos
estaciones <- as.data.frame(fread("estaciones.csv"))
summary(estaciones)## activate name reservations_count light
## Min. :1 Length:214 Min. :0 Min. :0.000
## 1st Qu.:1 Class :character 1st Qu.:0 1st Qu.:2.000
## Median :1 Mode :character Median :0 Median :2.000
## Mean :1 Mean :0 Mean :1.925
## 3rd Qu.:1 3rd Qu.:0 3rd Qu.:2.000
## Max. :1 Max. :0 Max. :3.000
## total_bases free_bases number longitude
## Min. :12.0 Min. : 0.00 Length:214 Min. :-3.725
## 1st Qu.:24.0 1st Qu.:10.00 Class :character 1st Qu.:-3.704
## Median :24.0 Median :12.00 Mode :character Median :-3.695
## Mean :23.9 Mean :11.63 Mean :-3.694
## 3rd Qu.:24.0 3rd Qu.:13.00 3rd Qu.:-3.683
## Max. :30.0 Max. :24.00 Max. :-3.656
## no_available address latitude dock_bikes
## Min. :0.000000 Length:214 Min. :40.39 Min. : 0.00
## 1st Qu.:0.000000 Class :character 1st Qu.:40.41 1st Qu.:10.00
## Median :0.000000 Mode :character Median :40.42 Median :12.00
## Mean :0.004673 Mean :40.42 Mean :11.37
## 3rd Qu.:0.000000 3rd Qu.:40.43 3rd Qu.:13.00
## Max. :1.000000 Max. :40.47 Max. :19.00
## id
## Min. : 1.00
## 1st Qu.: 57.25
## Median :112.50
## Mean :111.60
## 3rd Qu.:165.75
## Max. :219.00
# Nos quedamos con las variables que nos interesan y creamos el df
estaciones_origen <- estaciones[ , c(5,8,11,13)]
# Cambiamos los nombres
names(estaciones_origen)[4] <- "idunplug_station"
names(estaciones_origen)[1] <- "total_bases_origen"
names(estaciones_origen)[2] <- "longitude_origen"
names(estaciones_origen)[3] <- "latitude_origen"Habitualmente encontramos las coordenadas geográficas representadas por Latitud y Longitud separadas por una coma, por ejemplo: 40.4168961º, -3.7024255º. Creamos una nueva variable que así lo represente. Después, las cambiamos por su frecuencia para que entren correctamente en el modelo.
# Creamos coordenadas
estaciones_origen$coordenadas_origen <- paste(estaciones_origen$latitude_origen, "º", ", ", estaciones_origen$longitude_origen, "º", sep = "")
# Cambiamos por frecuencia. Antes nos aseguramos que es data.table
estaciones_origen <- as.data.table(estaciones_origen)
estaciones_origen[ , fe_coordenadas_origen := .N , by = .(coordenadas_origen)]
# Eliminamos original
estaciones_origen$coordenadas_origen <- NULLTambién creamos una variable que situe cada estación según lo lejos que esté del centro. Tomamos como referencia el Kilometro cero de la Puerta del Sol y calculamos la distancia en km que hay desde ese punto a cada estación.
estaciones_origen <- estaciones_origen%>% rowwise() %>%
mutate(dist_centro_origen = distHaversine(c(latitude_origen, longitude_origen), c(40.41664,-3.7059988))/1000)
head(estaciones_origen)Y unimos al dataframe con el que veníamos trabajando haciendo un merge tomando como referencia la estación de desanclaje.
Nuestro objetivo es predecir el número de viajes que se realizarán cada hora del día por estación. Para ello, creamos nuestra variable objetivo o target con el número de viajes que se realizan en cada estación en cada hora del día y le unimos el df que hemos trabajado.
# Contamos los viajes por estación por día y hora
target <- mov_con_clima_estaciones%>%count(idunplug_station, fecha, fe_hour, sort=TRUE)
# Nos quedamos con las filas únicas.
mov_con_clima_estaciones<-distinct(mov_con_clima_estaciones)
# Y unimos el df original
mov_con_clima_estaciones_target <-merge(target, mov_con_clima_estaciones, by = c("idunplug_station", "fecha", "fe_hour"))
# Cambiamos nombre a la target
names(mov_con_clima_estaciones_target)[4] <- "target"
# Visualización df final con la target
head(mov_con_clima_estaciones_target)## [1] 0
## 'data.frame': 1203308 obs. of 25 variables:
## $ idunplug_station : int 1 1 1 1 1 1 1 1 1 1 ...
## $ fecha : Date, format: "2019-01-01" "2019-01-01" ...
## $ fe_hour : int 12 15 16 17 18 19 20 22 23 9 ...
## $ target : int 1 1 2 2 9 1 1 1 2 1 ...
## $ fe_year : int 2019 2019 2019 2019 2019 2019 2019 2019 2019 2019 ...
## $ fe_month : int 1 1 1 1 1 1 1 1 1 1 ...
## $ fe_day : int 1 1 1 1 1 1 1 1 1 1 ...
## $ fe_wday : int 3 3 3 3 3 3 3 3 3 3 ...
## $ fe2_season : int 1582753 1582753 1582753 1582753 1582753 1582753 1582753 1582753 1582753 1582753 ...
## $ tmed : num 6.6 6.6 6.6 6.6 6.6 6.6 6.6 6.6 6.6 6.6 ...
## $ prec : chr "0" "0" "0" "0" ...
## $ tmin : num 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 ...
## $ tmax : num 11.7 11.7 11.7 11.7 11.7 11.7 11.7 11.7 11.7 11.7 ...
## $ dir : int 6 6 6 6 6 6 6 6 6 6 ...
## $ velmedia : num 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 ...
## $ racha : num 5.3 5.3 5.3 5.3 5.3 5.3 5.3 5.3 5.3 5.3 ...
## $ presMin : num 952 952 952 952 952 ...
## $ fe_prec_bin : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ fe_dir : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ fe_velmedia : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ total_bases_origen : int 30 30 30 30 30 30 30 30 30 30 ...
## $ longitude_origen : num -3.7 -3.7 -3.7 -3.7 -3.7 ...
## $ latitude_origen : num 40.4 40.4 40.4 40.4 40.4 ...
## $ fe_coordenadas_origen: int 1 1 1 1 1 1 1 1 1 1 ...
## $ dist_centro_origen : num 0.473 0.473 0.473 0.473 0.473 ...
A continuación hacemos el mismo proceso, pero para crear una nueva variable target que usaremos en el segundo modelo. En este caso el número de viajes por estación y día.
target2<-mov_con_clima_estaciones_target%>%group_by(idunplug_station,fecha)%>%summarise(n=sum(target)) # Agrupamos por fecha.
mov_con_clima_estaciones2<-mov_con_clima_estaciones[,-7] # Quitamos la variable "fe_hour"
mov_con_clima_estaciones2<-distinct(mov_con_clima_estaciones2) # Filtramos solo por las filas únicas.
mov_dia <-merge(target2, mov_con_clima_estaciones2, by = c("idunplug_station", "fecha")) # Juntamos la variable target con el data frame.
names(mov_dia)[3]<-"target"Ahora que tenemos nuestro data frame pulido y con nuevas variables que consideramos importantes, procedemos a hacer un modelo predictivo del número de viajes según la estación, el día y la hora.
Hacemos una partición de 70% train y 30% test.
Vamos a hacer un modelo de random forest, utilizando la implementación de la librería ranger.
modelo_ranger1 <- ranger(
formula = target~.,
data = data_train,
num.trees= 100,
importance = 'impurity',
verbose = TRUE)## Growing trees.. Progress: 22%. Estimated remaining time: 1 minute, 49 seconds.
## Growing trees.. Progress: 47%. Estimated remaining time: 1 minute, 9 seconds.
## Growing trees.. Progress: 71%. Estimated remaining time: 38 seconds.
## Growing trees.. Progress: 96%. Estimated remaining time: 5 seconds.
## Ranger result
##
## Call:
## ranger(formula = target ~ ., data = data_train, num.trees = 100, importance = "impurity", verbose = TRUE)
##
## Type: Regression
## Number of trees: 100
## Sample size: 842317
## Number of independent variables: 24
## Mtry: 4
## Target node size: 5
## Variable importance mode: impurity
## Splitrule: variance
## OOB prediction error (MSE): 4.921925
## R squared (OOB): 0.4466778
Este modelo ranger de tipo regresión, con 100 árboles, nos da un R squared (OOB) de 0.44. Pero este resultado es sobre el conjunto train. Debemos probar el modelo con el conjunto test, para comprobar la precisión del modelo.
Analizando la importancia que le da el modelo a las variables, destacan la hora y la estación de origen. Las variables de geolocalización también tienen bastante relevancia, así como las climáticas.
rang_ranger1 <- as.data.frame(importance(modelo_ranger1))
rang_ranger1$vars <- rownames(rang_ranger1)
names(rang_ranger1)[1]<- c('Importance')
imp_var<-rang_ranger1 %>%
arrange(desc(Importance))
imp_var_top<-head(imp_var,15)
ggplot(imp_var_top, aes(x=reorder(vars,-Importance), y=Importance))+
geom_bar(stat="identity",color="darkgreen", fill="lightgreen")+
theme(axis.text.x=element_text(size=12, angle=45, color="blue",face="bold.italic"))+
labs(x = "Variables",y = "Importancia")Predicción sobre test:
# Coeficiente de determinación de la predicción TEST
R2 <- 1 - (sum((data_test$target-pred_mod1$predictions)^2)/sum((data_test$target-mean(data_test$target))^2))
# Coeficientee de determinación ajustado.
n<-length(pred_mod1$predictions)
k<-pred_mod1$num.independent.variables
R2A<- 1-((n-1)/(n-k-1))*(1-R2)
R2## [1] 0.4514743
## [1] 0.4514378
El R2 ajustado en test es 0.45. Un resultado muy parecido al R squared OOB del modelo. No tenemos un R2 muy alto, pero sí un modelo que nos complementará la información que obtengamos en el siguiente.
Este modelo nos permitirá adelantarnos a la demanda de las estaciones de Bicimad para cada hora y día del año. Sabiendo cuales van a ser los picos de demanda durante el día, se optimizarían las estaciones ya que los usuarios nunca encontrarían una estación vacía.
En este caso, la partición la hacemos sobre nuestro data frame de número de viajes por día.
set.seed(1234)
trainIndex2 <- createDataPartition(mov_dia$target, p=0.7, list=FALSE)
data_train2 <- mov_dia[trainIndex2,]
data_test2 <- mov_dia[-trainIndex2,]Hacemos un modelo random forest (ranger), pero en este caso con 200 árboles. Tenemos menos datos por lo que el esfuerzo computacional es menor.
modelo_ranger2 <- ranger(
formula = target~.,
data = data_train2,
num.trees= 200,
importance = 'impurity',
verbose = TRUE)
modelo_ranger2## Ranger result
##
## Call:
## ranger(formula = target ~ ., data = data_train2, num.trees = 200, importance = "impurity", verbose = TRUE)
##
## Type: Regression
## Number of trees: 200
## Sample size: 52842
## Number of independent variables: 23
## Mtry: 4
## Target node size: 5
## Variable importance mode: impurity
## Splitrule: variance
## OOB prediction error (MSE): 240.8185
## R squared (OOB): 0.8151926
Predicción sobre test:
Vemos que la importancia de las variables no es muy diferente al modelo anterior excepto por la mayor relevancia de la longitud y la distancia de las estaciones al centro. Tampoco está la hora de desanclaje porque no entraba en el modelo (predicción por día).
imp2<-as.data.frame(modelo_ranger2$variable.importance)
imp2<-cbind(rownames(imp2), imp2$`modelo_ranger2$variable.importance`)
imp2<-as.data.frame(imp2)
imp2$V2<-as.numeric(imp2$V2)
imp2<-imp2[order(-imp2$V2),]
colnames(imp2)[1:2]<-c("Variable", "Importancia")
imp2# Coeficiente de determinación de la predicción TEST
R2 <- 1 - (sum((data_test2$target-pred_mod2$predictions)^2)/sum((data_test2$target-mean(data_test2$target))^2))
# Coeficientee de determinación ajustado.
n<-length(pred_mod2$predictions)
k<-pred_mod2$num.independent.variables
R2A<- 1-((n-1)/(n-k-1))*(1-R2)
R2## [1] 0.8150124
## [1] 0.8148243
El coeficiente de determinación del segundo modelo es muy similar al R squared OOB, 0,81. Tenemos una buena medida de ajuste.
Con este modelo no podemos hacer una predicción tan detallada de la demanda como hicimos en el primero (demanda por hora), en cambio su R2 es mucho más alto. Este modelo permitiría a las empresas reubicar las bicicletas para que haya siempre suficiente stock en las estaciones, y así satisfacer la demanda del día.