Archivos utilizados

• P48cierre: en formato XML, utilizaré el de cierre, es decir, el de una fecha ya pasada al día en curso. Es descargable desde la API pública de ESIOS y dispone del detalle horario de las Unidades de Programación (UP) para el día seleccionado. La definición ortodoxa es la siguiente: > Programa horario operativo cierre (P48CIERRE): Es el programa de las unidades de programación que queda establecido al finalizar el horizonte diario de programación y que contiene los programas resultantes del Programa Diario Base de Funcionamiento (PDBF), los ajustes de programa diario base de funcionamiento y de las diferentes sesiones del mercado intradiariorealizados en los mercados intradiarios, así como las modificaciones de los programas asociadas a los procesos de resoluciónsolución de restricciones técnicas y a la participación de las diferentes unidades en los servicios de regulación frecuencia-potencia del sistema y en el proceso de gestión de desvíos generación-consumo. Asimismo, incorpora los correspondientes redespachos generados como consecuencia de las declaraciones de indisponibilidad y de la comunicación de desvíos de las unidades de programación.
• Archivo que relaciona UP con sujetos de mercado: en formato csv, descargable también diariamente desde ESIOS.
• Archivo que relaciona sujeto de mercado con la descripción de la empresa o razón social: en formato xlsx, descargable desde OMIE.
• Archivo de festivos: elaborado manualmente con los calendarios oficiales.

Una vez tengo toda esta información, puedo crear una tabla base sobre la que trabajar:

library(tidyverse)
library(data.table)
library(lubridate)
library(kableExtra)
library(readxl)
library(tsibble)
library(feasts)

# Cargamos BBDD de P48
BBDD <- readRDS("/Users/pherreraariza/Documents/P48cierre/BBDD/p48cierreBBDD.rds") # %>% mutate(fecha = ymd(fecha))

# Cargamos archivo de sujetos de mercado y unidades de programación de ESIOS
sujetos <- fread("/Users/pherreraariza/Documents/P48cierre/input/export_unidades-de-programacion.csv")

# Cargamos archivo de sujetos de mercado y empresa de OMIE
empresas <- read_excel("/Users/pherreraariza/Documents/P48cierre/input/LISTA_AGENTES.xlsx")

# Cargamos festivos
festivos <- read_csv2("/Users/pherreraariza/Documents/ModeloESIOS/calendario/festivos.csv", col_types = c("cc")) %>% 
            mutate(fecha = dmy(fecha)) %>% dplyr::rename(tipo_dia = tipo)

## Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.

# Generamos tabla completa con P48, sujetos y festivos del año móvil
tabla <- BBDD %>% group_by(UP, tipo, hora, fecha, negocio) %>%
                  summarise(energia = sum(energia, na.rm = TRUE)) %>% 
                  ungroup %>%
                  arrange(fecha, hora, UP) %>%
                  left_join(sujetos, by = c("UP" = "Código de UP")) %>%
                  left_join(empresas, by = c("Sujeto del Mercado" = "CODIGO OMEL")) %>% 
                  left_join(festivos, by = "fecha") %>%
                  mutate(weekday = wday(fecha))

tabla$tipo_dia <- if_else(is.na(tabla$tipo_dia) & tabla$weekday %in% c(2,3,4,5,6), "L", "F")

El aspecto de la tabla es el siguiente:

glimpse(tabla)

## Observations: 43,125,312
## Variables: 20
## $ UP                   <chr> "AAYUC01", "ABA1", "ABA2", "ABO1", "ABO2", …
## $ tipo                 <chr> "C", "V", "V", "V", "V", "C", "C", "V", "C"…
## $ hora                 <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ fecha                <date> 2016-01-01, 2016-01-01, 2016-01-01, 2016-0…
## $ negocio              <chr> "mercado", "mercado", "mercado", "mercado",…
## $ energia              <dbl> 0.2, 0.1, 0.1, 245.0, 424.0, 0.0, 0.0, 168.…
## $ `Descripción corta`  <chr> NA, "ABA1", "ABA2", "ABOÑO 1", "ABOÑO 2", "…
## $ `Descripción larga`  <chr> NA, "ABARAN", "ABARAN 2", "C.T. ABOÑO 1", "…
## $ `Potencia máxima MW` <chr> NA, "0,3", "0,4", "341,7", "561,8", "0,1", …
## $ `Código EIC`         <chr> NA, "18WABA1-123456-K", "18WABA2-123456-A",…
## $ `Tipo de producción` <chr> NA, "Hidráulica UGH", "Hidráulica UGH", "Hu…
## $ Negocio              <chr> NA, "Venta", "Venta", "Venta", "Venta", "Co…
## $ `Zona de Regulación` <chr> NA, "SZR", "SZR", "HC", "HC", "SZR", "SZR",…
## $ `Sujeto del Mercado` <chr> NA, "BARER", "BARER", "HC G", "HC G", "ABOU…
## $ `Tipo de UP`         <chr> NA, "Generación", "Generación", "Generación…
## $ `CODIGO MINISTERIO`  <chr> NA, "BARER", "BARER", "HC G", "HC G", "ABOU…
## $ DESCRIPCIÓN          <chr> NA, "BARBO ENERGÍAS RENOVABLES S.A.", "BARB…
## $ `TIPO AGENTE`        <chr> NA, "GENERACIÓN", "GENERACIÓN", "GENERACIÓN…
## $ tipo_dia             <chr> "F", "F", "F", "F", "F", "F", "F", "F", "F"…
## $ weekday              <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6…

La variable negocio corresponde con el nodo del XML del P48 llamado TipoNegocio, que puede tomar dos valores: Z21 y Z24. El primero corresponde a energía a mercado y el segundo a bilateral. Sin embargo, el último valor Z24 data del 12 de junio de 2018.

Volúmenes de energía programados: a mercado y bilateral

Una parte de la energía programada no es comprada o vendida en el mercado OMIE, sino transaccionada entre dos unidades que han de declarar el volumen intercambiado comercialmente ante REE a un precio confidencial que únicamente las dos entidades conocen, según contrato establecido.

Por tanto, no puedo saber con el P48 entre qué UP se intercambia la energía, pero sí qué cuota representa y agruparla por tecnologías:

# Genero tabla sumarizada de energía por tipo de negocio
bilaterales <- tabla %>% group_by(negocio, tipo, año = year(fecha)) %>% summarise(energia = sum(energia, na.rm = TRUE))

# Grafico
bilaterales %>% ggplot(aes(x = año, y = energia, fill = tipo)) + geom_bar(stat="identity", position=position_dodge()) + facet_grid(. ~ negocio) + scale_fill_brewer(palette="Paired") + theme_minimal()

# Tabla
dcast(bilaterales, negocio + tipo ~ año, sum) %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))

Existe una discontinuidad de los datos a partir de 2018. Tras consultarlo con el servidio digame de REE, parece que debido a la entrada en funcionamiento del Mercado Intradiario Continuo (MIC), el valor que idenficaba los bilaterales se dejo de publicar, pasando a todo a asignarse a mercado (cuando esto es falso, ya que según el I90 que publica OMIE, obviamente sigue programándose energía por bilateales).

Ránking de energía gestionada durante el año móvil

Un primer punto de partida sería analizar las comercializadoras no incumbentes y graficar un ranking de la energía programada. Dicha energía no tiene por qué ser la medida real de sus puntos de suministro, ya que la reducción de la penalización por desvíos entre programa y medida puede venir por estimar correctamente el sentido de los mismos, siendo el operador del sistema el encargado de comprar o vender la energía deficitaria o excedentaria.

La agrupación más sencilla consiste en hacerla por sujeto de mercado, aunque dará lugar a una clasificación con ciertas particularidades, al no contemplar la energía total de un grupo o holding empresarial (caso de Energya VM, por citar algún ejemplo):

# Ranking de comercializadoras no incumbentes por energía total acumulada del año móvil
ranking_comerc_libres_no_incb <- tabla %>% mutate(energia = if_else(tipo == "C", energia, energia * (-1))) %>%
                                   filter(`Tipo de producción` == "Comercializadores mercado libre") %>% 
                                   filter(!`Sujeto del Mercado` %in% c("ENDCO", "IBCLI", "GNSE", "GNCO", "NATGA", "VIEE", "HC G", "HCENE")) %>% 
                                   filter(fecha >= last(fecha) - years(1)) %>%
                                   group_by(`Sujeto del Mercado`, DESCRIPCIÓN) %>%
                                   summarise(energia = sum(energia, na.rm = TRUE)/1000) %>%
                                   arrange(desc(energia)) %>% 
                                   ungroup() %>%
                                   mutate(sujeto = paste0(`Sujeto del Mercado`, " - ", DESCRIPCIÓN)) %>%
                                   select(sujeto, energia)

# Dotplot
ranking_comerc_libres_no_incb %>% top_n(50) %>%
                    arrange(energia) %>%
                    mutate(sujeto = factor(sujeto, sujeto)) %>%
                    ggplot(aes(x = sujeto, y = energia)) + 
                                geom_point(col="orange", size = 1) +   # Draw points
                                geom_segment(aes(x = sujeto, 
                                                 xend = sujeto, 
                                                 y = min(energia), 
                                                 yend = max(energia)), 
                                             linetype="dashed", 
                                             size=0.1) +   # Draw dashed lines
                                labs(title="Ranking comercializadoras no incumbentes", 
                                     subtitle="Energía programada", 
                                     caption="source: ESIOS",
                                     y = "GWh",
                                     x = "Sujeto de mercado") +  
                                coord_flip() +   
                                theme_light() +
                                theme(panel.grid.major.y = element_blank(),
                                    panel.border = element_blank(),
                                    axis.ticks.y = element_blank(),
                                    axis.text=element_text(size=5)) +
                                scale_y_continuous(labels = scales::comma)

ranking_comerc_libres_no_incb %>% top_n(50) %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))

Durante el año móvil (tomado desde el 2 de agosto de 2019) se diferencian varios grupos destacables:

• Fortia, con unca cartera de cerca de 10 TWh/año.
• Acciona y Fenie: con alrededor de 5 TWh/año.
• Energya VM, Axpo, Nexus, Aldro y Engie: alrededor de 3 TWh/año.
• Finalmente una gran cola con una distribución más suave que engloba al resto.

Evolución anual del top 10

Teniendo datos desde 2016, se puede ver la evolución mensual hasta julio de 2019 (incluido) para ver tendencias del top 10 de comercializadoras de mercado libre:

# Selecciono el top 10 del ranking anterior
top10 <- ranking_comerc_libres_no_incb %>% top_n(10) %>% select(sujeto) %>% pull()

# Creo la tabla con el año y mes para cada sujeto
evolucion <- tabla %>% mutate(energia = if_else(tipo == "C", energia, energia * (-1))) %>%
                       filter(`Tipo de producción` == "Comercializadores mercado libre" & fecha < ymd("2019-08-01")) %>% 
                       group_by(`Sujeto del Mercado`, DESCRIPCIÓN, año = year(fecha), mes = month(fecha)) %>%
                       summarise(energia = sum(energia, na.rm = TRUE)/1000) %>%
                       arrange(desc(energia)) %>% 
                       ungroup() %>%
                       mutate(sujeto = paste0(`Sujeto del Mercado`, " - ", DESCRIPCIÓN),
                              añomes = ymd(paste(año, mes, "01"))) %>%
                       filter(sujeto %in% top10) %>% 
                       select(añomes, sujeto, energia)

# Grafico
p <- evolucion %>% ggplot(aes(x = añomes, y = energia, fill = sujeto)) + 
                    geom_area(stat = "identity") + facet_wrap(sujeto ~ ., scales="free_y") + 
                    theme_minimal() + 
                    theme(legend.position="none") #, strip.text.y = element_blank())

p

Existen varias tipologías de tendencias a la vista de los datos en esta selección. Existen comercializadoras, como Fortia, Cide, con una clara estacionalidad de tipo mensual, donde la tendencia es plana. Otras, por otra parte, poseen una estacionalidad muy diluida pero con tendencias muy acentuadas, como Aldro y Fenie. Pero se observa en general que cada una tiene sus propios drivers y su propia evolución.

Estacionalidad

A continuación, para el análisis de las estacionalidad se puede graficar la energía programada para cada año. Al no estar normalizada o escalada, también tiene el efecto de crecimiento neto para cada compañía.

# Convierto tabla 'evolucion' en un time-series para utilizar la función de graficación de la estacionalidad
estacionalidad <- evolucion %>% mutate(indice = yearmonth(añomes)) %>% 
                                as_tsibble(index = indice, key = sujeto) %>% 
                                select(-añomes)

# Grafico
estacionalidad %>% gg_season(energia, labels = "both") +
                            ylab("GWh") +
                            ggtitle("Estacionalidad por sujeto")

Existe una estacionalidad variable, desde curvas semiplanas con tendencia lineal positiva, hasta curvas en forma de W, con picos en enero, julio/agosto y diciembre, propias de carteras domésticas.

Curvas de programa operativo por mes y año

Finalmente, se puede visualizar la curva horaria, distinguiendo entre días laborales y festivos + fines de semana. De esta manera se llega a la máxima granularidad a la hora de extraer posibles conclusiones.

# Verifico que hay una única UP para cada entrada de sujeto de mercado de tipo comercializador
top10_UP <- tabla %>% select(`Sujeto del Mercado`, UP, `Tipo de producción`, DESCRIPCIÓN) %>%
                       filter(`Tipo de producción` == "Comercializadores mercado libre") %>% 
                       mutate(sujeto = paste0(`Sujeto del Mercado`, " - ", DESCRIPCIÓN)) %>%
                       filter(sujeto %in% top10) %>%
                       select(sujeto, UP) %>%
                       distinct()

top10_UP %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))

# Tabla de curvas horarias de programa para las UP del top 10
curvas <-  tabla %>% mutate(energia = if_else(tipo == "C", energia, energia * (-1))) %>%
                     filter(UP %in% top10_UP$UP & fecha < ymd("2019-08-01")) %>%
                     group_by(fecha, UP, hora, tipo_dia) %>%
                     summarise(energia = sum(energia, na.rm = TRUE))

ALDRO

El primer impacto visual viene dado por el fuerte crecimiento en energía (cartera) que presenta desde el Q4 de 2018, con una curva horaria típica de las tarifas 3.0A, con un pico a mediodía y forma decreciente sin pico nocturno. Además, excepto por los meses de octubre y noviembre de 2018, las curvas son muy estables.

curvas %>% filter(UP == "ALDRC01") %>% ggplot(aes(x=hora, y=energia, group = fecha, color = tipo_dia)) + geom_line(alpha=0.2) + facet_grid(year(fecha) ~ month(fecha)) + ggtitle("ALDRO")

AUDAX

Se observa un perfil de curva típicamente de empresa, bastante estable en general, exceptuando los meses de agosto y diciembre. Al igual que Aldro, hay una clara diferenciación entre festivos y laborables.

curvas %>% filter(UP == "AUDAC01") %>% ggplot(aes(x=hora, y=energia, group = fecha, color = tipo_dia)) + geom_line(alpha=0.2) + facet_grid(year(fecha) ~ month(fecha)) + ggtitle("AUDAX")

ENERGYA VM

Perfil muy parecido al de Audax. La única diferencia aparente radica en la evolución de la altura de la curva.

curvas %>% filter(UP == "CENTC01") %>% ggplot(aes(x=hora, y=energia, group = fecha, color = tipo_dia)) + geom_line(alpha=0.2) + facet_grid(year(fecha) ~ month(fecha)) + ggtitle("ENERGYA VM")

CIDE

Este perfil tiene poca diferenciación entre festivo y laborable, y presenta un doble pico de consumo a lo largo del día, lo que indica que se trata eminentemente de una cartera de tipo doméstico.

curvas %>% filter(UP == "CHCC01") %>% ggplot(aes(x=hora, y=energia, group = fecha, color = tipo_dia)) + geom_line(alpha=0.2) + facet_grid(year(fecha) ~ month(fecha)) + ggtitle("CIDE")

ENGIE

Este perfil es casi plano, con diferenciación entre festivo y laborable acentuado en enero y diciembre durante 2016 y que se va dispersando en los siguientes años. Se trata principalmente de una cartera industrial, con procesos estables, continuos que generan este tipo de curvas.

curvas %>% filter(UP == "EBLCO01") %>% ggplot(aes(x=hora, y=energia, group = fecha, color = tipo_dia)) + geom_line(alpha=0.2) + facet_grid(year(fecha) ~ month(fecha)) + ggtitle("ENGIE")

ACCIONA

Esta unidad es de las que más dispersión presenta, pese a tratarse de una cartera eminentemente industrial como la anterior. La curva deja de ser plana para pasar a incorporar una zona valle con un menor consumo de energía.

curvas %>% filter(UP == "EGEDC01") %>% ggplot(aes(x=hora, y=energia, group = fecha, color = tipo_dia)) + geom_line(alpha=0.2) + facet_grid(year(fecha) ~ month(fecha)) + ggtitle("ACCIONA")

AXPO

Perfil similar a Engie, estable excepto los meses de agosto, donde paran los procesos productivos.

curvas %>% filter(UP == "EGLEC01") %>% ggplot(aes(x=hora, y=energia, group = fecha, color = tipo_dia)) + geom_line(alpha=0.2) + facet_grid(year(fecha) ~ month(fecha)) + ggtitle("AXPO")

FENIE

Fenie presenta uno de los perfiles con más casuísticas: tiene un crecimiento mensual muy acentuado entre 2017 y 2018, para después volver a programar menos energía durante el primer semestre de 2019. Además, la forma de la curva es un mix entre las de doméstico y empresa, donde posiblemente los industriales simplemente eleven o disminuyan la altura de la misma. Llama la atención la dispersión atípica de los meses de enero y febrero, que podría deberse a finalización de contratos de industriales de consumo plano.

curvas %>% filter(UP == "FENC01") %>% ggplot(aes(x=hora, y=energia, group = fecha, color = tipo_dia)) + geom_line(alpha=0.2) + facet_grid(year(fecha) ~ month(fecha)) + ggtitle("FENIE")

FORTIA

Pefil muy marcado con elevado consumo durante las primeras horas de la noche, para bajar casi un 50% después. Llama la atención que durante los festivos se programe más energía que durante los días laborables, a diferencia de lo que se ha visto en el resto de casos.

curvas %>% filter(UP == "FORTC01") %>% ggplot(aes(x=hora, y=energia, group = fecha, color = tipo_dia)) + geom_line(alpha=0.2) + facet_grid(year(fecha) ~ month(fecha)) + ggtitle("FORTIA")

NEXUS

Nexus presenta un perfil de empresa acentuado, muy estable y con poca dispersión, y con clara diferenciación entre festivo y laborable.

curvas %>% filter(UP == "NEXUC01") %>% ggplot(aes(x=hora, y=energia, group = fecha, color = tipo_dia)) + geom_line(alpha=0.2) + facet_grid(year(fecha) ~ month(fecha)) + ggtitle("NEXUS")