Visualización de bilaterales físicos en el sistema eléctrico

Descarga del archivo I90

Este archivo tipo excel contiene multitud de información sobre energía transaccionada, tanto a nivel de ofertas en el mercado mayorista, como de ejecución de contratos bilaterales, los cuales son trazables a través del código de contrato que identifica a cada uno. Un vistazo general al índice de este archivo muestra el tipo de información que contiene:

Concretamente, la hoja 27 contiene la información que busco:

library(RCurl)
library(plyr)
library(tidyverse)
library(lubridate)
library(httr)
library(reshape2)
library(data.table)
library(readxl)
library(viridis)
library(kableExtra)

Fijo las fechas inicio y fin para descargar mediante un loop cada I90 del mes de enero de 2020, concretamente la hoja 27, y exportarlo como csv para analizarlo mejor (la descarga no es especialmente rápida, pero en 1 minuto lo tendremos):

# Fechas de publicación del I90 diario
fecha_final <-  ymd("20200101")
fecha_inicial <- ymd("20200131")
dias <- seq(fecha_final, fecha_inicial, 1)

# Loop para descargar los I90 diarios
lapply(dias, function(x) {
  
  # URL
  uridownlaod <- paste0("https://api.esios.ree.es/archives/34/download?start_date=", x, "T00:00:00+00:00&end_date=", x, "T23:59:59+00:00")
  
  # Tratamos respuesta y contenido de XML a dataframe
  message(paste0("Descargando I90DIA del día ", x))
  #response <- GET(uridownlaod, httpheader = httpheader)
  temp <- tempfile()
  download.file(uridownlaod ,temp)  
  unzip(temp, paste0("I90DIA_", year(x), format.Date(x, "%m"), format.Date(x, "%d"), ".xls"), exdir = "temp")
  I90dia27 <- read_xls(paste0("temp/I90DIA_", year(x), format.Date(x, "%m"), format.Date(x, "%d"), ".xls"), sheet = "I90DIA27", skip = 3)

  # Exportamos
  write_csv2(I90dia27, path = paste0("/Users/pherreraariza/Documents/energychisquared/post_inputs/10_post/i90/I90DIA_", x, ".csv"))
  message(paste0("Generado I90 del día ", x))
  unlink(temp)
  
})

Se leen los ficheros y se genera una única tabla:

# Lectura de ficheros
lista_I90_bilaterales <- list.files("/Users/pherreraariza/Documents/energychisquared/post_inputs/10_post/i90", pattern = ".csv", full.names = TRUE)

BBDD_I90 <- rbindlist(use.names = TRUE, fill = TRUE, lapply(lista_I90_bilaterales, function(x) {
                              I90_data <- read_csv2(x, locale = locale(decimal_mark = "."))
                              I90_data$año <- substr(x, 80, 83)
                              I90_data$mes <- substr(x, 85, 86)
                              I90_data$dia <- substr(x, 88, 89)
                              I90_data$fecha <- ymd(paste0(I90_data$año, I90_data$mes, I90_data$dia))
                              I90_data <- I90_data %>% select(-año, -mes, -dia)
                            }))

Los datos en bruto tienen este aspecto:

glimpse(BBDD_I90)

## Rows: 6,021
## Columns: 30
## $ `Unidad de Programación` <chr> "AGRIC01", "ALBCC01", "ALBCC01", "ALBCC…
## $ `Tipo Oferta`            <dbl> 8, 8, 8, 8, 1, 1, 1, 1, 8, 1, 1, 1, 8, …
## $ `Nº contrato`            <dbl> 1984, 2018, 2019, 2025, 1929, 1994, 192…
## $ Hora                     <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Total                    <dbl> -541.0, 3.8, 7.4, 2.4, 12624.0, 7776.0,…
## $ `00-01`                  <dbl> -18.0, 0.1, 0.3, 0.1, 526.0, 324.0, 525…
## $ `01-02`                  <dbl> -17.0, 0.1, 0.3, 0.1, 526.0, 324.0, 525…
## $ `02-03`                  <dbl> -16.0, 0.1, 0.3, 0.1, 526.0, 324.0, 525…
## $ `03-04`                  <dbl> -16.0, 0.1, 0.3, 0.1, 526.0, 324.0, 525…
## $ `04-05`                  <dbl> -16.0, 0.1, 0.2, 0.1, 526.0, 324.0, 525…
## $ `05-06`                  <dbl> -17.0, 0.1, 0.2, 0.1, 526.0, 324.0, 525…
## $ `06-07`                  <dbl> -22.0, 0.2, 0.3, 0.1, 526.0, 324.0, 525…
## $ `07-08`                  <dbl> -25.0, 0.2, 0.3, 0.1, 526.0, 324.0, 525…
## $ `08-09`                  <dbl> -27.0, 0.2, 0.3, 0.1, 526.0, 324.0, 525…
## $ `09-10`                  <dbl> -27.0, 0.2, 0.3, 0.1, 526.0, 324.0, 525…
## $ `10-11`                  <dbl> -27.0, 0.2, 0.3, 0.1, 526.0, 324.0, 525…
## $ `11-12`                  <dbl> -27.0, 0.2, 0.3, 0.1, 526.0, 324.0, 525…
## $ `12-13`                  <dbl> -27.0, 0.2, 0.3, 0.1, 526.0, 324.0, 525…
## $ `13-14`                  <dbl> -25.0, 0.2, 0.3, 0.1, 526.0, 324.0, 525…
## $ `14-15`                  <dbl> -23.0, 0.2, 0.3, 0.1, 526.0, 324.0, 525…
## $ `15-16`                  <dbl> -22.0, 0.1, 0.3, 0.1, 526.0, 324.0, 525…
## $ `16-17`                  <dbl> -25.0, 0.1, 0.3, 0.1, 526.0, 324.0, 525…
## $ `17-18`                  <dbl> -26.0, 0.1, 0.3, 0.1, 526.0, 324.0, 525…
## $ `18-19`                  <dbl> -27.0, 0.1, 0.3, 0.1, 526.0, 324.0, 525…
## $ `19-20`                  <dbl> -26.0, 0.2, 0.4, 0.1, 526.0, 324.0, 525…
## $ `20-21`                  <dbl> -24.0, 0.2, 0.4, 0.1, 526.0, 324.0, 525…
## $ `21-22`                  <dbl> -22.0, 0.2, 0.4, 0.1, 526.0, 324.0, 525…
## $ `22-23`                  <dbl> -20.0, 0.2, 0.4, 0.1, 526.0, 324.0, 525…
## $ `23-24`                  <dbl> -19.0, 0.2, 0.3, 0.1, 526.0, 324.0, 525…
## $ fecha                    <date> 2020-01-01, 2020-01-01, 2020-01-01, 20…

Data wrangling

Una vez se dispone de los datos, se deben separar los comprados y los vendedores, para luego casar los contratos por su ID y agruparlos por mes para su posterior visualización:

## Tratamos compardores y vendedores del bilateral
compradores <- BBDD_I90 %>% filter(Total < 0) %>% 
                            select(`Unidad de Programación`, `Tipo Oferta`, fecha, `Nº contrato`) %>%
                            rename(UP_compradora = `Unidad de Programación`,
                                   Tipo_oferta_compradora = `Tipo Oferta`)

vendedores <- BBDD_I90 %>% filter(Total > 0) %>% 
                           select(`Unidad de Programación`, `Tipo Oferta`, Total, fecha, `Nº contrato`) %>%
                           rename(UP_vendedora = `Unidad de Programación`,
                                  Tipo_oferta_vendedora = `Tipo Oferta`)

## Casamos contratos
bilaterales <- compradores %>% left_join(vendedores, by = c("fecha", "Nº contrato"))

bilaterales_mensual <- bilaterales %>% mutate(año = year(fecha),
                                              mes = month(fecha)) %>% 
                                       group_by(UP_compradora, Tipo_oferta_compradora, UP_vendedora, Tipo_oferta_vendedora, año, mes) %>%
                                       summarise(energia = sum(Total, na.rm = T))

glimpse(bilaterales_mensual)

## Rows: 117
## Columns: 7
## Groups: UP_compradora, Tipo_oferta_compradora, UP_vendedora, Tipo_oferta_vendedora, año [117]
## $ UP_compradora          <chr> "AGRIC01", "AQUIC01", "ATLSC01", "CATRC01…
## $ Tipo_oferta_compradora <dbl> 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,…
## $ UP_vendedora           <chr> "NEXUC01", "EFGC01", "NEXUC01", "EALBC01"…
## $ Tipo_oferta_vendedora  <dbl> 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,…
## $ año                    <int> 2020, 2020, 2020, 2020, 2020, 2020, 2020,…
## $ mes                    <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ energia                <dbl> 15522.3, 950.6, 29749.1, 171.5, 1736.3, 1…

Visualización de bilaterales

Aquí viene el verdadero reto de implementar una visualización cómoda en R. Necesito un visualización que me pueda aportar información tanto del flujo de energía mensual transaccionada entre unidades de programación, como de su volumen. Una buena opción, flexible (y algo compleja al principio de utilizar) es a través de diagramas de cuerdas. Para cada ejemplo, aislaré la unidad a analizar tanto de las UP vendedoras como de las compradoras, para ver el flujo de energía de dicha unidad.

# Input de datos para el diagrama de cuerdas de ENDESA
adjacency_list_bilaterales_ENDESA <- bilaterales_mensual %>% 
                                              ungroup() %>% 
                                              filter(energia > 1000) %>%
                                              filter(UP_vendedora == "ENDE01" | UP_compradora == "ENDE01") %>%
                                                        select(UP_vendedora, UP_compradora, energia)

library(chorddiag)

# Creación de la matriz de datos
m1 <- acast(adjacency_list_bilaterales_ENDESA, UP_vendedora ~ UP_compradora, value.var = "energia", sum)

un1 <- unique(sort(c(colnames(m1), rownames(m1))))
m2 <- matrix(0, length(un1), length(un1), dimnames = list(un1, un1))
m2[row.names(m1), colnames(m1)] <- m1


# Construcción del diagrama de cuerdas
n <- dim(m2)[1]

groupColors <- viridis_pal()(n)           # para los grupos de colores
groupColors <- substr(groupColors, 0, 7)  # eliminamos las transparencias ('FF')

plot_endesa <- chorddiag(m2, type = "directional", groupColors = groupColors)
plot_endesa

library(circlize)

## ========================================
## circlize version 0.4.8
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: http://jokergoo.github.io/circlize_book/book/
## 
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization 
##   in R. Bioinformatics 2014.
## ========================================

# Tabla de colores
col <-  rand_color(nrow(adjacency_list_bilaterales_ENDESA))

# Creación del diagrama de cuerdas
chordDiagram(adjacency_list_bilaterales_ENDESA, annotationTrack = c("grid", "axis"), preAllocateTracks = 1, col = col,
             directional = 1, direction.type = c("diffHeight", "arrows"),
    link.arr.type = "big.arrow", diffHeight = -0.01)
  circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
  xlim = get.cell.meta.data("xlim")
  xplot = get.cell.meta.data("xplot")
  ylim = get.cell.meta.data("ylim")
  sector.name = get.cell.meta.data("sector.index")
    circos.text(mean(xlim), ylim[1], sector.name, facing = "clockwise",
                niceFacing = TRUE, adj = c(-0.5, 0.5), cex = 0.5)
}, bg.border = NA)

# Limpieza de las parametrizaciones
circos.clear()

# Tabla
adjacency_list_bilaterales_ENDESA %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))

# Input de datos para el diagrama de cuerdas de STATKRAFT
adjacency_list_bilaterales_STATKRAFT <- bilaterales_mensual %>% 
                                              ungroup() %>% 
                                              filter(energia > 1000) %>%
                                              filter(UP_vendedora == "GSMG1" | UP_compradora == "GSMG1") %>%
                                                        select(UP_vendedora, UP_compradora, energia)

# Tabla de colores
col <-  rand_color(nrow(adjacency_list_bilaterales_STATKRAFT))

# Creación del diagrama de cuerdas
chordDiagram(adjacency_list_bilaterales_STATKRAFT, annotationTrack = c("grid", "axis"), preAllocateTracks = 1, col = col,
             directional = 1, direction.type = c("diffHeight", "arrows"),
    link.arr.type = "big.arrow", diffHeight = -0.01)
  circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
  xlim = get.cell.meta.data("xlim")
  xplot = get.cell.meta.data("xplot")
  ylim = get.cell.meta.data("ylim")
  sector.name = get.cell.meta.data("sector.index")
    circos.text(mean(xlim), ylim[1], sector.name, facing = "clockwise",
                niceFacing = TRUE, adj = c(-0.5, 0.5), cex = 0.5)
}, bg.border = NA)

# Limpieza de las parametrizaciones
circos.clear()

# Tabla
adjacency_list_bilaterales_STATKRAFT %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))

# Input de datos para el diagrama de cuerdas de STATKRAFT
adjacency_list_bilaterales_NEXUS <- bilaterales_mensual %>% 
                                              ungroup() %>% 
                                              filter(energia > 1000) %>%
                                              filter(UP_vendedora == "NEXUC01" | UP_compradora == "NEXUC01") %>%
                                                        select(UP_vendedora, UP_compradora, energia)

# Tabla de colores
col <-  rand_color(nrow(adjacency_list_bilaterales_NEXUS))

# Creación del diagrama de cuerdas
chordDiagram(adjacency_list_bilaterales_NEXUS, annotationTrack = c("grid", "axis"), preAllocateTracks = 1, col = col,
             directional = 1, direction.type = c("diffHeight", "arrows"),
    link.arr.type = "big.arrow", diffHeight = -0.01)
  circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
  xlim = get.cell.meta.data("xlim")
  xplot = get.cell.meta.data("xplot")
  ylim = get.cell.meta.data("ylim")
  sector.name = get.cell.meta.data("sector.index")
    circos.text(mean(xlim), ylim[1], sector.name, facing = "clockwise",
                niceFacing = TRUE, adj = c(-0.5, 0.5), cex = 0.5)
}, bg.border = NA)

# Limpieza de las parametrizaciones
circos.clear()

# Tabla
adjacency_list_bilaterales_NEXUS %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))

Y para qué tanto lío…

Para disponer de información de la competencia, por supuesto. Cierto es que este tipo de análisis de seguimiento siempre ocurren a posteriori (90 días ni más ni menos), pero permite ver por qué un sujeto liquida un volumen de energía mayor al programa final de su unidad de compra (NEXUS), cómo se abastecen las incumbentes (ENDESA) o si los PPA’s firmados están en ejecución o no (STATKRAFT).