Series de Tiempo y Ciencia de Datos 2

Predicción (forecasting) de la demanda energética con machine learning

Introducción

La predicción de la demanda energética desempeña un papel fundamental en la gestión y planificación de los recursos necesarios para la generación, distribución y utilización de la energía. Predecir la demanda de energía es una tarea compleja en la que influyen factores como los patrones meteorológicos, las condiciones económicas y el comportamiento de la sociedad. Este documento muestra cómo utilizar modelos de machine learning para predecir la demanda de energía.

Series temporales y forecasting

Una serie temporal (time series) es una sucesión de datos ordenados cronológicamente, espaciados a intervalos iguales o desiguales. El proceso de forecasting consiste en predecir el valor futuro de una serie temporal, bien modelando la serie únicamente en función de su comportamiento pasado (autorregresivo) o empleando otras variables externas.

Librerías

Las librerías utilizadas en este documento son:

# Tratamiento de datos
# ==============================================================================
import numpy as np
import pandas as pd
from astral.sun import sun
from astral import LocationInfo
from skforecast.datasets import fetch_dataset
from feature_engine.datetime import DatetimeFeatures
from feature_engine.creation import CyclicalFeatures
from feature_engine.timeseries.forecasting import WindowFeatures

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
from skforecast.plot import plot_residuals
import plotly.graph_objects as go
import plotly.io as pio
import plotly.offline as poff
pio.templates.default = "seaborn"
poff.init_notebook_mode(connected=True)

{'text/html': '        <script type="text/javascript">\n        window.PlotlyConfig = {MathJaxConfig: \'local\'};\n        if (window.MathJax && window.MathJax.Hub && window.MathJax.Hub.Config) {window.MathJax.Hub.Config({SVG: {font: "STIX-Web"}});}\n        </script>\n        <script type="module">import "https://cdn.plot.ly/plotly-3.0.1.min"</script>\n        '}

plt.style.use('seaborn-v0_8-darkgrid')
plt.rcParams.update({'font.size': 8})

# Modelado y Forecasting
# ==============================================================================
import skforecast
import lightgbm
import sklearn
from lightgbm import LGBMRegressor
from sklearn.preprocessing import PolynomialFeatures
from sklearn.feature_selection import RFECV
from skforecast.recursive import ForecasterEquivalentDate, ForecasterRecursive
from skforecast.direct import ForecasterDirect
from skforecast.model_selection import TimeSeriesFold, bayesian_search_forecaster, backtesting_forecaster
from skforecast.feature_selection import select_features
from skforecast.preprocessing import RollingFeatures
from skforecast.plot import calculate_lag_autocorrelation, plot_residuals
from skforecast.metrics import calculate_coverage
import shap

# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')

Datos

Se dispone de una serie temporal de la demanda de electricidad (MW) para el estado de Victoria (Australia) desde 2011-12-31 hasta 2014-12-31. Los datos empleados en este documento se han obtenido del paquete de R tsibbledata. El set de datos contiene 5 columnas y 52608 registros completos. La información de cada columna es:

• Time: fecha y hora del registro.

• Date: fecha del registro

• Demand: demanda de electricidad (MW).

• Temperature: temperatura en Melbourne, capital de Victoria.

• Holiday: indicador si el día es festivo (vacaciones).

# Descarga de datos
# ==============================================================================
datos = fetch_dataset(name='vic_electricity', raw=True)

vic_electricity
---------------
Half-hourly electricity demand for Victoria, Australia
O'Hara-Wild M, Hyndman R, Wang E, Godahewa R (2022).tsibbledata: Diverse
Datasets for 'tsibble'. https://tsibbledata.tidyverts.org/,
https://github.com/tidyverts/tsibbledata/.
https://tsibbledata.tidyverts.org/reference/vic_elec.html
Shape of the dataset: (52608, 5)

datos.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 52608 entries, 0 to 52607
Data columns (total 5 columns):
 #   Column       Non-Null Count  Dtype  
---  ------       --------------  -----  
 0   Time         52608 non-null  object 
 1   Demand       52608 non-null  float64
 2   Temperature  52608 non-null  float64
 3   Date         52608 non-null  object 
 4   Holiday      52608 non-null  bool   
dtypes: bool(1), float64(2), object(2)
memory usage: 1.7+ MB

La columna Time está almacenada como string. Para convertirla en datetime, se emplea la función pd.to_datetime(). Una vez en formato datetime, y para hacer uso de las funcionalidades de pandas, se establece como índice. Además, dado que los datos se han registrado cada 30 minutos, se indica la frecuencia '30min'.

# Conversión del formato fecha
# ==============================================================================
datos['Time'] = pd.to_datetime(datos['Time'], format='%Y-%m-%dT%H:%M:%SZ')
datos = datos.set_index('Time')
datos = datos.asfreq('30min')
datos = datos.sort_index()
datos.head(2)

                          Demand  Temperature        Date  Holiday
Time                                                              
2011-12-31 13:00:00  4382.825174        21.40  2012-01-01     True
2011-12-31 13:30:00  4263.365526        21.05  2012-01-01     True

Uno de los primeros análisis que hay que realizar al trabajar con series temporales es verificar si la serie está completa.

# Verificar que un índice temporal está completo
# ==============================================================================
fecha_inicio = datos.index.min()
fecha_fin = datos.index.max()
date_range_completo = pd.date_range(start=fecha_inicio, end=fecha_fin, freq=datos.index.freq)
print(f"Índice completo: {(datos.index == date_range_completo).all()}")

Índice completo: True

print(f"Filas con valores ausentes: {datos.isnull().any(axis=1).mean()}")

Filas con valores ausentes: 0.0

# Completar huecos en un índice temporal
# ==============================================================================
# datos.asfreq(freq='30min', fill_value=np.nan)

Aunque los datos se encuentran en intervalos de 30 minutos, el objetivo es crear un modelo capaz de predecir la demanda eléctrica a nivel horario, por lo que se tienen que agregar los datos. Este tipo de transformación es muy sencillas si se combina el índice DatetimeIndex de pandas y su método resample().

Es muy importante utilizar correctamente los argumentos closed=‘left’ y label=‘right’ para no introducir en el entrenamiento información a futuro (leakage)). Supóngase que se dispone de valores para las 10:10, 10:30, 10:45, 11:00, 11:12 y 11:30. Si se quiere obtener el promedio horario, el valor asignado a las 11:00 debe calcularse utilizando los valores de las 10:10, 10:30 y 10:45; y el de las 12:00, con el valor de las 11:00, 11:12 y 11:30.

Para el valor promedio de las 11:00 no se incluye el valor puntual de las 11:00 por que, en la realidad, en ese momento exacto no se dispone todavía del valor.

# Agregado en intervalos de 1H
# ==============================================================================
# Se elimina la columna Date para que no genere error al agregar.
datos = datos.drop(columns='Date')
datos = (
    datos
    .resample(rule="h", closed="left", label="right")
    .agg({
        "Demand": "mean",
        "Temperature": "mean",
        "Holiday": "mean",
    })
)
datos

                          Demand  Temperature  Holiday
Time                                                  
2011-12-31 14:00:00  4323.095350       21.225      1.0
2011-12-31 15:00:00  3963.264688       20.625      1.0
2011-12-31 16:00:00  3950.913495       20.325      1.0
2011-12-31 17:00:00  3627.860675       19.850      1.0
2011-12-31 18:00:00  3396.251676       19.025      1.0
...                          ...          ...      ...
2014-12-31 09:00:00  4069.625550       21.600      0.0
2014-12-31 10:00:00  3909.230704       20.300      0.0
2014-12-31 11:00:00  3900.600901       19.650      0.0
2014-12-31 12:00:00  3758.236494       18.100      0.0
2014-12-31 13:00:00  3785.650720       17.200      0.0

[26304 rows x 3 columns]

El set de datos empieza el 2011-12-31 14:00:00 y termina el 2014-12-31 13:00:00. Se descartan los primeros 10 y los últimos 13 registros para que empiece el 2012-01-01 00:00:00 y termine el 2014-12-30 23:00:00. Además, para poder optimizar los hiperparámetros del modelo y evaluar su capacidad predictiva, se dividen los datos en 3 conjuntos, uno de entrenamiento, uno de validación y otro de test.

# Separación datos train-val-test
# ==============================================================================
datos = datos.loc['2012-01-01 00:00:00':'2014-12-30 23:00:00', :].copy()
fin_train = '2013-12-31 23:59:00'
fin_validacion = '2014-11-30 23:59:00'
datos_train = datos.loc[: fin_train, :].copy()
datos_val   = datos.loc[fin_train:fin_validacion, :].copy()
datos_test  = datos.loc[fin_validacion:, :].copy()

print(f"Fechas train      : {datos_train.index.min()} --- {datos_train.index.max()}  (n={len(datos_train)})")

Fechas train      : 2012-01-01 00:00:00 --- 2013-12-31 23:00:00  (n=17544)

print(f"Fechas validacion : {datos_val.index.min()} --- {datos_val.index.max()}  (n={len(datos_val)})")

Fechas validacion : 2014-01-01 00:00:00 --- 2014-11-30 23:00:00  (n=8016)

print(f"Fechas test       : {datos_test.index.min()} --- {datos_test.index.max()}  (n={len(datos_test)})")

Fechas test       : 2014-12-01 00:00:00 --- 2014-12-30 23:00:00  (n=720)

Exploración gráfica

La exploración gráfica de series temporales es una forma eficaz de identificar tendencias, patrones y estacionalidad. Esto, a su vez, ayuda a orientar la selección del modelo de forecasting más adecuado.

Gráfico de la serie temporal

Serie temporal completa

# Gráfico interactivo de la serie temporal
# ==============================================================================
fig = go.Figure()
fig.add_trace(go.Scatter(x=datos_train.index, y=datos_train['Demand'], mode='lines', name='Train'))

fig.add_trace(go.Scatter(x=datos_val.index, y=datos_val['Demand'], mode='lines', name='Validation'))

fig.add_trace(go.Scatter(x=datos_test.index, y=datos_test['Demand'], mode='lines', name='Test'))

fig.update_layout(
    title  = 'Demanda eléctrica horaria',
    xaxis_title="Fecha",
    yaxis_title="Demanda (MWh)",
    legend_title="Partición:",
    width=800,
    height=400,
    margin=dict(l=20, r=20, t=35, b=20),
    legend=dict(orientation="h", yanchor="top", y=1, xanchor="left", x=0.001)
)

#fig.update_xaxes(rangeslider_visible=True)
fig.show()