Problema a resover
El problema que se plantea es optimizar la integración de baterías que funcionen como respaldo en el sistema eléctrico mexicano. Suponiendo que estas baterías se cargan de manera directa desde la red y entregan energía cuando el sistema enfrenta situaciones críticas, es decir, cuando la generación disponible es baja, asegurando así el cumplimiento de la demanda.
Un punto clave a destacar en el mercado eléctrico es que la demanda siempre debe igualar a la oferta. Cualquier desviación de esta igualdad podría provocar el colapso del sistema. Para evitar esto, existe un operador, el CENACE (Centro Nacional de Control de Energía), encargado de asegurar que esta condición se cumpla en todo momento.
Con la creciente incorporación de energías renovables en los sistemas energéticos a nivel global, la incorporación de reservas se ha vuelto cada vez más necesaria. Sin embargo, en México, la regulación para su implementación aún está en revisión. En este proyecto, se realiza un análisis utilizando datos públicos disponibles sobre precios y condiciones del sistema, con el objetivo de evaluar la viabilidad de integrar estas baterías en el sistema eléctrico mexicano.
Supuestos
Para nuestro anláisis supondremos que tenemos uan batería que carga en dos horas completamente y descarga completamente en dos horas (Batería de 0.5C de carga y 0.5C de descarga). Para fines prácticos supondremos que no hay una tasa de eficiencia ni un mínimo de carga. La batería que entrega 1 MWh ya que MWh es la unidad de medida a la que se determina el precio.
Un MW no es significativo en la demanda del sistema eléctrico.
Vemamos la siguiente gráfica
print(descriptiva_region)## # A tibble: 9 × 6
## region demanda_mean demanda_sd demanda_median demanda_max demanda_min
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BCA 1944. 576. 1749. 3585. 936.
## 2 BCS 382. 109. 352. 669. 206.
## 3 CEN 6577. 917. 6693. 8468. 111.
## 4 NES 7230. 1572. 7121. 10978. 164.
## 5 NOR 3228. 1106. 2907. 5763. 65.5
## 6 NTE 3613. 779. 3473. 5440. 127.
## 7 OCC 8996. 1226. 9111. 12134. 234.
## 8 ORI 6905. 920. 6883. 9362. 218.
## 9 PEN 1918. 397. 1928. 2948. 43.9Lo primero que debemos saber sobre la determinación de precios de electricidad es que se varían de forma horaria y de acuerdo a la ubicación.
En la siguiente gráfica se muestran los precios promedio por región,
según el mapa previamente presentado. Los datos utilizados provienen de
la base de datos publicada por CENACE. En general, podemos observar que
los precios tienden a aumentar durante la noche y suelen ser más bajos
por la mañana, aunque este comportamiento puede variar dependiendo de la
región.
`
Para el primero definimos el Profit Score como:
\[ \left( P_{\text{max}_1} + P_{\text{max}_2} \right) - \left( P_{\text{min}_1} + P_{\text{min}_2} \right) \]
Dado que nuestro objetivo es maximizar la ganancia que puede obtener la batería en un día, asumimos que la batería se carga durante las dos horas con el precio más bajo (al cual se le cobra la electricidad) y se descarga durante las dos horas con el precio más alto (al cual se le paga por la electricidad que entrega).
En la primera base de datos que analizamos, contamos con información sobre los precios de electricidad para el año 2024. A continuación, realizaremos un pequeño análisis descriptivo de la información.
print(descriptiva_region)## # A tibble: 9 × 6
## region ps_mean ps_sd ps_median ps_max ps_min
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 BCA 590. 1701. 258. 17657. 9.13
## 2 BCS 3418. 2656. 2668. 18263. 496.
## 3 CEN 3669. 4441. 2025. 21707. 90.2
## 4 NES 3307. 4297. 1655. 21603. 70.7
## 5 NOR 3209. 4167. 1389. 21493. 123.
## 6 NTE 3603. 4516. 1873. 23399. 114.
## 7 OCC 3599. 4426. 1905. 22099. 109.
## 8 ORI 4068. 4706. 2569. 23383. 108.
## 9 PEN 6210. 5689. 4642. 30836. 116.Utilizamos una analis ANOVA para determinar que region en media es significativamente mayor que las demás porque esta sería la región que nos estaría dado el mayor beneficio ya que la batería se puede llevar un diferencial mayor.
summary(comparaciones_tukey)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = ps ~ region, family = poisson, data = ANOVA_DATA)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## BCS - BCA == 0 2.457319 0.003417 719.048 <1e-06 ***
## CEN - BCA == 0 2.195245 0.003468 632.932 <1e-06 ***
## NES - BCA == 0 2.056354 0.003494 588.476 <1e-06 ***
## NOR - BCA == 0 2.072488 0.003489 594.036 <1e-06 ***
## NTE - BCA == 0 2.204240 0.003465 636.193 <1e-06 ***
## OCC - BCA == 0 2.167926 0.003473 624.183 <1e-06 ***
## ORI - BCA == 0 2.323966 0.003447 674.169 <1e-06 ***
## PEN - BCA == 0 2.912624 0.003375 862.950 <1e-06 ***
## CEN - BCS == 0 -0.262074 0.001442 -181.748 <1e-06 ***
## NES - BCS == 0 -0.400965 0.001503 -266.702 <1e-06 ***
## NOR - BCS == 0 -0.384831 0.001490 -258.193 <1e-06 ***
## NTE - BCS == 0 -0.253079 0.001433 -176.586 <1e-06 ***
## OCC - BCS == 0 -0.289393 0.001454 -199.089 <1e-06 ***
## ORI - BCS == 0 -0.133353 0.001390 -95.927 <1e-06 ***
## PEN - BCS == 0 0.455305 0.001201 379.216 <1e-06 ***
## NES - CEN == 0 -0.138891 0.001616 -85.957 <1e-06 ***
## NOR - CEN == 0 -0.122757 0.001604 -76.542 <1e-06 ***
## NTE - CEN == 0 0.008995 0.001551 5.801 <1e-06 ***
## OCC - CEN == 0 -0.027319 0.001570 -17.406 <1e-06 ***
## ORI - CEN == 0 0.128721 0.001511 85.189 <1e-06 ***
## PEN - CEN == 0 0.717379 0.001339 535.872 <1e-06 ***
## NOR - NES == 0 0.016135 0.001659 9.724 <1e-06 ***
## NTE - NES == 0 0.147886 0.001608 91.970 <1e-06 ***
## OCC - NES == 0 0.111572 0.001626 68.609 <1e-06 ***
## ORI - NES == 0 0.267612 0.001570 170.480 <1e-06 ***
## PEN - NES == 0 0.856270 0.001405 609.577 <1e-06 ***
## NTE - NOR == 0 0.131752 0.001596 82.557 <1e-06 ***
## OCC - NOR == 0 0.095438 0.001614 59.122 <1e-06 ***
## ORI - NOR == 0 0.251477 0.001557 161.476 <1e-06 ***
## PEN - NOR == 0 0.840135 0.001391 604.051 <1e-06 ***
## OCC - NTE == 0 -0.036314 0.001561 -23.256 <1e-06 ***
## ORI - NTE == 0 0.119726 0.001503 79.678 <1e-06 ***
## PEN - NTE == 0 0.708384 0.001329 532.922 <1e-06 ***
## ORI - OCC == 0 0.156040 0.001522 102.517 <1e-06 ***
## PEN - OCC == 0 0.744698 0.001351 551.129 <1e-06 ***
## PEN - ORI == 0 0.588658 0.001283 458.908 <1e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)Haciendo los ajustes llegamos a que todas las medias a nivel región son diferentes entre si a excepción de la Occidental, Norte y Central. Por lo menos la región Peninsular es significativamente más grande que las demás.
Como se determina el precio en mercado
Las plantas de generación se ordenan de acuerdo con su costo variable (principalmente costos de combustibles), de menor a mayor. Se da preferencia a las plantas con costos más bajos, como las de energía renovable (solar, eólica) y las hidroeléctricas.
El precio que se paga a los generadores es el que corresponde al punto de intersección entre la demanda y la oferta, conocido como el Precio Marginal Local.
- Principio de Minimización de Costos:
El objetivo es cubrir la demanda eléctrica al menor costo posible mientras se garantizan las condiciones de confiabilidad y estabilidad del sistema.
- Precio Marginal Local (PML):
El costo del despacho económico se refleja en el PML, que varía según la ubicación y las restricciones del sistema. Este precio incluye:
- Costo de energía.
- Congestión en la red.
- Pérdidas técnicas.
- Restricciones Operativas:
Aunque el despacho sigue el orden de mérito, se consideran factores como:
- Límites de capacidad de transmisión.
- Restricciones técnicas de las plantas.
- Necesidades de estabilidad del sistema.
Entonces que puede afectar el precio: cambios en la demanda, cambios en generación disponible, salidas de generación por mantenimeinto y paros forzados.
Adicionalmente, existen dos tipos de despacho: el MDA (Mercado Diario de Energía), donde los precios se determinan a través de pronósticos, y el MTR (Mercado de Tiempo Real), donde los precios se fijan con base en datos reales. En el MDA, los precios se publican ex post, mientras que en el MTR, los precios se determinan de forma anticipada para el día siguiente.
La normativa actual establece que las baterías se despachan en tiempo real, lo cual representa una ventaja competitiva, ya que las condiciones del día de adelanto son conocidas con antelación y se pueden usar para pronosticar el precio en el MTR.
Colinealidad
## [1] 0.1584971## [1] "No hay evidencia fuerte de colinealidad."## [1] 2.4643434 1.4939231 0.7860672 0.6713825 0.3535499 0.2307340## [1] "No hay colinealidad según los valores propios."Regresión
receta<- recipe(pml_mtr~ pml_mda+region+hora+forzadas+pronostico_mwh, data = datos_sin_outliers) %>%
step_dummy(region) %>%
step_normalize(all_numeric_predictors()) %>%
prep()
receta_juice<-juice(receta)
# Ajustamos modelo
modelo1 <- linear_reg() %>%
set_mode("regression") %>%
set_engine("lm")
lm_fit1 <- train(pml_mtr ~ ., data = receta_juice, method = "lm")##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6447.1 -264.5 -90.1 105.7 8729.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 996.896 2.254 442.219 < 2e-16 ***
## pml_mda 511.361 2.967 172.339 < 2e-16 ***
## hora 84.356 2.398 35.184 < 2e-16 ***
## forzadas 24.976 6.817 3.664 0.000249 ***
## pronostico_mwh 301.272 7.711 39.069 < 2e-16 ***
## region_BCS 572.280 3.878 147.589 < 2e-16 ***
## region_CEN -54.422 5.310 -10.248 < 2e-16 ***
## region_NES -99.769 8.188 -12.185 < 2e-16 ***
## region_NOR 10.865 3.254 3.339 0.000843 ***
## region_NTE 15.790 3.323 4.752 2.02e-06 ***
## region_OCC -130.635 6.628 -19.708 < 2e-16 ***
## region_ORI -58.388 5.575 -10.473 < 2e-16 ***
## region_PEN 179.290 3.133 57.219 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 607.5 on 72611 degrees of freedom
## Multiple R-squared: 0.6817, Adjusted R-squared: 0.6817
## F-statistic: 1.296e+04 on 12 and 72611 DF, p-value: < 2.2e-16Regresión Lógistica
Aplicamos dos regresiones logísticas para determinar, entre las horas del día, cuáles tienen mayor probabilidad de ser las de precios más altos y cuáles las de precios más bajos.
Primero, asignamos una etiqueta a los datos históricos: asignamos 1 si la hora correspondió a una de las dos horas de precio más alto en tiempo real y 0 si no fue así. De manera análoga, realizamos el mismo proceso para identificar las horas con precios más bajos.
data_completa<-data
lista_fecha<-as.list(unique(data_completa$fecha))
lista_zona<-as.list(unique(data_completa$region))
data_variable<-NULL
data_logic<-NULL
for(i in lista_fecha){
for(j in lista_zona){
data_variable=data_completa%>%filter(fecha==i, region==j)
data_variable=data_variable%>% arrange(hora)%>%
mutate(indicador_hora_max=ifelse(pml_mtr ==max_pnd_1|pml_mtr ==max_pnd_2,1,0),
indicador_hora_min=ifelse(pml_mtr ==min_pnd_1|pml_mtr==min_pnd_2,1,0))
data_logic<-rbind(data_logic,data_variable)
}
# print(i)
}
data_logic<-data_logic%>%mutate(indicador_hora_max=as.factor(indicador_hora_max),
indicador_hora_min=as.factor(indicador_hora_min),
region=as.factor(region))Se necesita una sub-base de entrenamiento y una para las pruebas 80% entrenamiento y 20%.
# set.seed(1234)
# data_split <- initial_split(data_logic, prop = 0.8)
# data_train <- training(data_split)
# data_test <- testing(data_split)
data_train<-data_logic
data_test<-data_logicGeneramos el modelo de regresión logística para el máximo:
receta_logic<- recipe(indicador_hora_max~ pml_mda+fecha+hora+region+pronostico_mwh+forzadas, data = data_train) %>%
step_rm(fecha)%>%
step_dummy(all_nominal_predictors()) %>%
# step_normalize(all_numeric_predictors())%>%
prep()
receta_logic_juice<-juice(receta_logic)
logistic_model <- logistic_reg() %>%
set_engine("glm")
logistic_fit1 <- parsnip::fit(logistic_model, indicador_hora_max ~ ., receta_logic_juice)
predicciones <- predict(logistic_fit1,receta_logic_juice,type="prob")
head(predicciones)## # A tibble: 6 × 2
## .pred_0 .pred_1
## <dbl> <dbl>
## 1 0.992 0.00772
## 2 0.991 0.00888
## 3 0.990 0.0103
## 4 0.989 0.0113
## 5 0.987 0.0132
## 6 0.985 0.0154logistic_p_test <- predict(logistic_fit1, receta_logic_juice,type="prob") %>%
bind_cols(receta_logic_juice) Algoritmo de asignación
Dado que una regresión logística asigna probabilidades dentro de un rango, es posible que un día tenga más de dos horas asignadas o ninguna. Por lo tanto, nos centramos únicamente en aquellas horas del día que tienen la mayor probabilidad de ser las de precios más altos..
# lista_regiones<-list("region_PEN","region_ORI")
lista_dias<-as.list(unique(data_test$fecha))
data_variable<-NULL
data_predict<-NULL
for(i in lista_dias){
data_variable=data_test%>%filter(fecha==i)
# receta aplicada al subconjunto
receta_logic<- recipe(indicador_hora_max~ pml_mda+fecha+hora+region+pronostico_mwh+forzadas, data = data_variable)%>%
step_rm(fecha)%>%
step_dummy(all_nominal_predictors()) %>%
# step_normalize(all_numeric_predictors())%>%
prep()
receta_logic_juice<-juice(receta_logic)
logistic_p_test <- predict(logistic_fit1, receta_logic_juice,type="prob") %>%
bind_cols(receta_logic_juice) %>%mutate(fecha=i, indicador_hora_max=as.numeric(indicador_hora_max))
data_predict<-rbind(data_predict,logistic_p_test)
# print(i)
}
data_predict<-data_predict%>%gather(key="region",value="indicador",8:15)%>%filter(indicador==1)%>%
mutate(region=substr(region,8,nchar(region)))%>%dplyr::select(-indicador)
# Calculamos la proba
data_predict<-data_predict%>%group_by(fecha,region) %>%
mutate(max_pred_1 = max(.pred_1))%>%
ungroup()
segundo_valor_max <- data_predict %>% group_by(fecha,region) %>%
arrange(desc(.pred_1)) %>%
slice(2)%>%dplyr::select(fecha,region,.pred_1)%>%
rename(max_pred_2=.pred_1)
data_predict<-merge(data_predict,segundo_valor_max,by=c("fecha","region"))lista_fecha<-as.list(unique(data_predict$fecha))
lista_zona<-as.list(unique(data_predict$region))
data_variable<-NULL
data_predict_max<-NULL
for(i in lista_fecha){
for(j in lista_zona){
data_variable=data_predict%>%filter(fecha==i, region==j)
data_variable=data_variable%>% arrange(hora)%>%
mutate(indicador_hora_max_predict=ifelse(.pred_1 ==max_pred_1|.pred_1 ==max_pred_2,1,0))
data_predict_max<-rbind(data_predict_max,data_variable)
}
# print(i)
}Lo mismo pero horas mínimas
Generamos una regresión logística para el precio mínimo
receta_logic<- recipe(indicador_hora_min~ pml_mda+fecha+hora+region+pronostico_mwh+forzadas, data = data_train) %>%
step_rm(fecha)%>%
step_dummy(all_nominal_predictors()) %>%
# step_normalize(all_numeric_predictors())%>%
prep()
receta_logic_juice<-juice(receta_logic)
logistic_model <- logistic_reg() %>%
set_engine("glm")
logistic_fit1 <- parsnip::fit(logistic_model, indicador_hora_min ~ ., receta_logic_juice)
predicciones <- predict(logistic_fit1,receta_logic_juice,type="prob")
head(predicciones)## # A tibble: 6 × 2
## .pred_0 .pred_1
## <dbl> <dbl>
## 1 0.926 0.0743
## 2 0.917 0.0828
## 3 0.920 0.0804
## 4 0.852 0.148
## 5 0.861 0.139
## 6 0.868 0.132logistic_p_test <- predict(logistic_fit1, receta_logic_juice,type="prob") %>%
bind_cols(receta_logic_juice) # lista_regiones<-list("region_PEN","region_ORI")
lista_dias<-as.list(unique(data_test$fecha))
data_variable<-NULL
data_predict<-NULL
for(i in lista_dias){
data_variable=data_test%>%filter(fecha==i)
# receta aplicada al subconjunto
receta_logic<- recipe(indicador_hora_min~ pml_mda+fecha+hora+region+pronostico_mwh+forzadas, data = data_variable)%>%
step_rm(fecha)%>%
step_dummy(all_nominal_predictors()) %>%
# step_normalize(all_numeric_predictors())%>%
prep()
receta_logic_juice<-juice(receta_logic)
logistic_p_test <- predict(logistic_fit1, receta_logic_juice,type="prob") %>%
bind_cols(receta_logic_juice) %>%mutate(fecha=i, indicador_hora_min=as.numeric(indicador_hora_min))
data_predict<-rbind(data_predict,logistic_p_test)
# print(i)
}
data_predict<-data_predict%>%gather(key="region",value="indicador",8:15)%>%filter(indicador==1)%>%
mutate(region=substr(region,8,nchar(region)))%>%dplyr::select(-indicador)
# Calculamos la proba
data_predict<-data_predict%>%group_by(fecha,region) %>%
mutate(min_pred_1 = max(.pred_1))%>%
ungroup()
segundo_valor_min<- data_predict %>% group_by(fecha,region) %>%
arrange(desc(.pred_1)) %>%
slice(2)%>%dplyr::select(fecha,region,.pred_1)%>%
rename(min_pred_2=.pred_1)
data_predict<-merge(data_predict,segundo_valor_min,by=c("fecha","region"))
data_predict<-data_predict%>%dplyr::select(fecha,hora,region,min_pred_1,min_pred_2,indicador_hora_min,.pred_1,.pred_0)lista_fecha<-as.list(unique(data_predict$fecha))
lista_zona<-as.list(unique(data_predict$region))
data_variable<-NULL
data_predict_min<-NULL
for(i in lista_fecha){
for(j in lista_zona){
data_variable=data_predict%>%filter(fecha==i, region==j)
data_variable=data_variable%>% arrange(hora)%>%
mutate(indicador_hora_min_predict=ifelse(.pred_1 ==min_pred_1|.pred_1 ==min_pred_2,1,0))
data_predict_min<-rbind(data_predict_min,data_variable)
}
# print(i)
}Ganancia
data_predict<-merge(data_predict_max,data_predict_min,by=c("fecha","hora","region"))
data_pml_mtr<-data%>%dplyr::select(fecha,hora,region,pml_mtr)
data_predict<-merge(data_pml_mtr,data_predict,by=c("fecha","hora","region"))
ingreso<-data_predict%>%filter(indicador_hora_max_predict==1)%>%
group_by(region,fecha)%>%
summarise(ingreso=sum(pml_mtr))
costo<-data_predict%>%filter(indicador_hora_min_predict==1)%>%
group_by(region,fecha)%>%
summarise(egreso=sum(pml_mtr))
resultado<-merge(ingreso,costo,by=c("fecha","region"))%>%mutate(utilidad=ingreso-egreso)%>%filter(utilidad>0)%>%group_by(region)%>%summarise(utilidad=sum(utilidad))
resultado_esperado<-data%>%group_by(region,fecha)%>%filter(region!="BCA")%>%summarise(utilidad=mean(ps))%>%group_by(region)%>%
summarise(utilidad=sum(utilidad))El resultado que ubieramos obtenido si le atinamos al total de los valores mayores
print(resultado_esperado)## # A tibble: 8 × 2
## region utilidad
## <chr> <dbl>
## 1 BCS 1247405.
## 2 CEN 1335616.
## 3 NES 1203639.
## 4 NOR 1167910.
## 5 NTE 1311365.
## 6 OCC 1310157.
## 7 ORI 1480779.
## 8 PEN 2260365.Contra el resultado que ubieramos obtenido con nuestro modelo
print(resultado)## # A tibble: 8 × 2
## region utilidad
## <chr> <dbl>
## 1 BCS 413422.
## 2 CEN 467504.
## 3 NES 412488.
## 4 NOR 385929.
## 5 NTE 456784.
## 6 OCC 446216.
## 7 ORI 516687.
## 8 PEN 1050689.Conslusión y comentarios
Además del diferencial de precios, la batería recibe un pago por un concepto conocido como Potencia, que es un ingreso anual otorgado por estar disponible durante las horas críticas, las cuales coinciden con las horas de precio máximo. Este pago varía anualmente, pero generalmente ronda los 4,600,000 MXN por MW de capacidad instalada.
En nuestro modelo, asumimos un 50% de éxito y contamos con una capacidad de 1 MW, lo que nos permite estimar un ingreso anual de 2,300,000 MXN por este concepto. La vida útil de la batería es de aproximadamente 20 años, y su costo de inversión se encuentra entre 12,000,000 y 20,000,000 MXN, por lo que consideramos un costo promedio de inversión de 16,000,000 MXN. Esto nos brinda una Tasa Interna de Retorno (TIR) aproximada del 20%, con un periodo de recuperación de la inversión de alrededor de 5 años.
Cabe destacar que los precios en tiempo real presentan dificultades para su estimación. En el modelo presentado, únicamente se utilizó información pública. Sin embargo, con la adición de datos provenientes de proveedores externos, áreas seguras del CENACE, y la aplicación de algoritmos más sofisticados, se podría alcanzar una rentabilidad máxima de hasta 43%.
El esquema de despacho de baterías utilizado en nuestra modelación aún no se encuentra operativo en México. Actualmente, las disposiciones para su implementación están en análisis; sin embargo, el esquema propuesto presenta un despacho similar al utilizado en nuestro algoritmo y se estima que comenzará a operar en febrero de 2025.
Ingresar rápidamente a este mercado representa una ventaja competitiva significativa, ya que permitiría aprovechar el diferencial de precios antes de que las baterías tengan un impacto considerable en la estimación de precios.