Si te gusta Skforecast , ayúdanos dándonos una estrella en GitHub! ⭐️
Predicción (forecasting) de la demanda eléctrica con Python
Más sobre ciencia de datos: cienciadedatos.net
- Forecasting series temporales con Python y Scikit Learn
- Predicción (forecasting) de visitas a página web con machine learning
- Forecasting series temporales con gradient boosting: Skforecast, XGBoost, LightGBM y CatBoost
- Predicción del precio de Bitcoin con Python
- Intervalos de predicción en modelos de forecasting
- Multi-series forecasting
- Predicción de demanda intermitente
Introducción¶
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.
Cuando se trabaja con series temporales, raramente se quiere predecir solo el siguiente elemento de la serie ($t_{+1}$), sino todo un intervalo futuro ($t_{+1}$), ..., ($t_{+n}$)) o un punto alejado en el tiempo ($t_{+n}$). Existen varias estrategias que permiten generar este tipo de predicciones múltiples, la librería skforecast recoge las siguientes para series temporales univariantes:
Forecasting multi-step recursivo: dado que, para predecir el momento $t_{n}$ se necesita el valor de $t_{n-1}$, y $t_{n-1}$ se desconoce, se aplica un proceso recursivo en el que, cada nueva predicción, se basa en la anterior. Este proceso se conoce como predicción recursiva o predicción multipaso recursiva y puede generarse fácilmente con las clases
ForecasterAutoregyForecasterAutoregCustom.Forecasting multi-step directo: este método consiste en entrenar un modelo diferente para cada valor futuro (step) del horizonte de predicción. Por ejemplo, para predecir los 5 valores siguientes de una serie temporal, se entrenan 5 modelos diferentes, uno para cada step. De este modo, las predicciones son independientes entre sí. Todo este proceso se automatiza en la clase
ForecasterAutoregDirect.
En este documento se muestra un ejemplo de cómo utilizar métodos de forecasting para predecir la demanda eléctrica a nivel horario. En concreto, se hace uso de skforecast, una librería que contiene las clases y funciones necesarias para adaptar cualquier modelo de regresión de scikit-learn a problemas de forecasting.
Para casos de uso más detallados visitar skforecast-examples
Caso de uso¶
Se dispone de una serie temporal con la demanda eléctrica (MW) del estado de Victoria (Australia) desde 2011-12-31 al 2014-12-31. Se pretende generar un modelo de forecasting capaz de predecir la demanda energética del día siguiente a nivel horario.
# Tratamiento de datos
# ==============================================================================
import numpy as np
import pandas as pd
# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
plt.style.use('seaborn-v0_8-darkgrid')
# Modelado y Forecasting
# ==============================================================================
from sklearn.linear_model import Ridge
from lightgbm import LGBMRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_absolute_error
from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.ForecasterAutoregDirect import ForecasterAutoregDirect
from skforecast.model_selection import grid_search_forecaster
from skforecast.model_selection import backtesting_forecaster
from sklearn.preprocessing import OneHotEncoder
from sklearn.compose import make_column_transformer
from sklearn.compose import make_column_selector
# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')
Datos¶
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
# ==============================================================================
url = ('https://raw.githubusercontent.com/JoaquinAmatRodrigo/skforecast/master/' +
'data/vic_elec.csv')
datos = pd.read_csv(url, sep=',')
datos.info()
La columna Time se ha almacenado 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)
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
# ==============================================================================
(datos.index == pd.date_range(start=datos.index.min(),
end=datos.index.max(),
freq=datos.index.freq)).all()
print(f"Filas con missing values: {datos.isnull().any(axis=1).mean()}")
# 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 de tipo temporal 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. La columna Holiday
# no genera error ya que es booleana y se trata como 0-1.
datos = datos.drop(columns='Date')
datos = datos.resample(rule='H', closed='left', label ='right').mean()
datos
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)})")
print(f"Fechas validacion : {datos_val.index.min()} --- {datos_val.index.max()} (n={len(datos_val)})")
print(f"Fechas test : {datos_test.index.min()} --- {datos_test.index.max()} (n={len(datos_test)})")
Exploración gráfica¶
Cuando se quiere generar un modelo de forecasting, es importante representar los valores de la serie temporal. Esto permite identificar patrones tales como tendencias y estacionalidad.
Serie temporal completa
# Gráfico serie temporal
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 3.5))
datos_train.Demand.plot(ax=ax, label='entrenamiento', linewidth=1)
datos_val.Demand.plot(ax=ax, label='validación', linewidth=1)
datos_test.Demand.plot(ax=ax, label='test', linewidth=1)
ax.set_title('Demanda eléctrica')
ax.legend();
El gráfico anterior muestra que la demanda eléctrica tiene estacionalidad anual. Se observa un incremento centrado en el mes de Julio y picos de demanda muy acentuados entre enero y marzo.
Sección de la serie temporal
Debido a la varianza de la serie temporal, no es posible apreciar con un solo gráfico el posible patrón intradiario.
# Gráfico serie temporal con zoom
# ==============================================================================
zoom = ('2013-05-01 14:00:00','2013-06-01 14:00:00')
fig = plt.figure(figsize=(8, 4))
grid = plt.GridSpec(nrows=8, ncols=1, hspace=0.6, wspace=0)
main_ax = fig.add_subplot(grid[1:3, :])
zoom_ax = fig.add_subplot(grid[5:, :])
datos.Demand.plot(ax=main_ax, c='black', alpha=0.5, linewidth=0.5)
min_y = min(datos.Demand)
max_y = max(datos.Demand)
main_ax.fill_between(zoom, min_y, max_y, facecolor='blue', alpha=0.5, zorder=0)
main_ax.set_xlabel('')
datos.loc[zoom[0]: zoom[1]].Demand.plot(ax=zoom_ax, color='blue', linewidth=1)
main_ax.set_title(f'Demanda electricidad: {datos.index.min()}, {datos.index.max()}', fontsize=10)
zoom_ax.set_title(f'Demanda electricidad: {zoom}', fontsize=10)
zoom_ax.set_xlabel('')
plt.subplots_adjust(hspace=1)
Al aplicar zoom sobre la serie temporal, se hace patente una clara estacionalidad semanal, con consumos más elevados durante la semana laboral (lunes a viernes) y menor en los fines de semana. Se observa también que existe una clara correlación entre el consumo de un día con el del día anterior.
Estacionalidad anual, semanal y diaria
# Gráfico boxplot para estacionalidad anual
# ==============================================================================
fig, ax = plt.subplots(figsize=(7, 2.5))
datos['mes'] = datos.index.month
datos.boxplot(column='Demand', by='mes', ax=ax,)
datos.groupby('mes')['Demand'].median().plot(style='o-', linewidth=0.8, ax=ax)
ax.set_ylabel('Demand')
ax.set_title('Distribución demanda por mes')
fig.suptitle('');
Se observa que hay una estacionalidad anual, con valores de demanda (mediana) superiores en los meses de Junio, Julio y Agosto, y con elevados picos de demanda en los meses de Noviembre, Diciembre, Enero, Febrero y Marzo.
# Gráfico boxplot para estacionalidad semanal
# ==============================================================================
fig, ax = plt.subplots(figsize=(7, 2.5))
datos['dia_semana'] = datos.index.day_of_week + 1
datos.boxplot(column='Demand', by='dia_semana', ax=ax)
datos.groupby('dia_semana')['Demand'].median().plot(style='o-', linewidth=0.8, ax=ax)
ax.set_ylabel('Demand')
ax.set_title('Distribución demanda por día de la semana')
fig.suptitle('');
Se aprecia una estacionalidad semanal, con valores de demanda inferiores durante el fin de semana.
# Gráfico boxplot para estacionalidad diaria
# ==============================================================================
fig, ax = plt.subplots(figsize=(7, 2.5))
datos['hora_dia'] = datos.index.hour + 1
datos.boxplot(column='Demand', by='hora_dia', ax=ax)
datos.groupby('hora_dia')['Demand'].median().plot(style='o-', linewidth=0.8, ax=ax)
ax.set_ylabel('Demand')
ax.set_title('Distribución demanda por hora del día')
fig.suptitle('');
También existe una estacionalidad diaria, la demanda se reduce entre las 16 y las 21 horas.
Días festivos y no festivos
# Grafico violinplot
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 3))
sns.violinplot(
x = 'Demand',
y = 'Holiday',
data = datos.assign(Holiday = datos.Holiday.astype(str)),
palette = 'tab10',
ax = ax
)
ax.set_title('Distribución de la demanda entre festivos y no festivos')
ax.set_xlabel('demanda')
ax.set_ylabel('festivo');
Los días festivos tienden a tener menor consumo.
Gráficos de autocorrelación
# Gráfico autocorrelación
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_acf(datos.Demand, ax=ax, lags=60)
plt.show()
# Gráfico autocorrelación parcial
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_pacf(datos.Demand, ax=ax, lags=60)
plt.show()
Los gráficos de autocorrelación y autocorrelación parcial muestran una clara asociación entre la demanda de una hora y las horas anteriores, así como entre la demanda de una hora y la demanda de esa misma hora los días anteriores. Este tipo de correlación, es un indicativo de que los modelos autorregresivos pueden funcionar bien.
Modelo autoregresivo recursivo¶
Se crea y entrena un modelo autorregresivo recursivo (ForecasterAutoreg) a partir de un modelo de regresión lineal con penalización Ridge y una ventana temporal de 24 lags. Esto último significa que, para cada predicción, se utilizan como predictores la demanda de las 24 horas anteriores.
Entrenamiento del Forecaster¶
Se utiliza un StandardScaler como transformer para el preprocesado.
# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
regressor = Ridge(),
lags = 24,
transformer_y = StandardScaler()
)
forecaster.fit(y=datos.loc[:fin_validacion, 'Demand'])
forecaster
Predicción (backtest)¶
Se evalúa el comportamiento que habría tenido el modelo si se hubiese entrenado con los datos desde 2012-01-01 00:00 al 2014-11-30 23:59 y, después, a las 23:59 de cada día, se predijesen las 24 horas siguientes. A este tipo de evaluación se le conoce como backtesting, y puede aplicarse fácilmente con la función backtesting_forecaster(). Esta función devuelve, además de las predicciones, una métrica de error.
Note
A partir de la versión skforecast 0.9.0, todas las funciones de backtesting y búsqueda de hiperparámetros del módulomodel_selection incluyen el parámetro n_jobs que permite su paralelización y, por tanto, mejora su rendimiento.
Los beneficios de la paralelización dependen de varios factores, incluyendo el regresor seleccionado, el número de entrenamientos a realizar y el volumen de datos involucrados. Cuando el parámetro n_jobs se establece en 'auto', el nivel de paralelización se selecciona automáticamente basándose en reglas heurísticas que pretenden elegir la mejor configuración posible para cada escenario.
# Backtest
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['Demand'],
steps = 24,
metric = 'mean_absolute_error',
initial_train_size = len(datos.loc[:fin_validacion]),
refit = False,
n_jobs = 'auto',
verbose = True,
show_progress = True
)
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 3.5))
datos.loc[predicciones.index, 'Demand'].plot(ax=ax, linewidth=2, label='test')
predicciones.plot(linewidth=2, label='predicción', ax=ax)
ax.set_title('Predicción vs demanda real')
ax.legend();
# Error backtest
# ==============================================================================
print(f'Error backtest: {metrica}')
Optimización de hiperparámetros (tuning)¶
En el objeto ForecasterAutoreg entrenado, se han utilizado los primeros 24 lags y un modelo Ridge con los hiperparámetros por defecto. Sin embargo, no hay ninguna razón por la que estos valores sean los más adecuados.
Con el objetivo de identificar la mejor combinación de lags e hiperparámetros, se recurre a un Grid Search. Este proceso consiste en entrenar un modelo con cada combinación de hiperparámetros y lags, y evaluar su capacidad predictiva mediante backtesting. En el proceso de búsqueda, es importante evaluar los modelos utilizando únicamente los datos de validación y no incluir los de test, estos se utilizan solo en último lugar para evaluar al modelo final.
# Grid search de hiperparámetros
# ==============================================================================
forecaster = ForecasterAutoreg(
regressor = Ridge(),
lags = 24, # Este valor será remplazado en el grid search
transformer_y = StandardScaler()
)
# Lags utilizados como predictores
lags_grid = [5, 24, [1, 2, 3, 23, 24, 25, 47, 48, 49]]
# Hiperparámetros del regresor
param_grid = {'alpha': np.logspace(-3, 5, 10)}
resultados_grid = grid_search_forecaster(
forecaster = forecaster,
y = datos.loc[:fin_validacion, 'Demand'],
steps = 24,
metric = 'mean_absolute_error',
param_grid = param_grid,
lags_grid = lags_grid,
initial_train_size = len(datos[:fin_train]),
refit = False,
n_jobs = 'auto',
return_best = True,
verbose = False
)
# Resultados Grid Search
# ==============================================================================
resultados_grid
Los mejores resultados se obtienen si se utilizan los lags [1, 2, 3, 23, 24, 25, 47, 48, 49] y una configuración de Ridge {'alpha': 215.44}. Al indicar return_best = True en la función grid_search_forecaster(), al final del proceso, se reentrena automáticamente el objeto forecaster con la mejor configuración encontrada y el set de datos completo (train + validation).
forecaster
Backtest con el conjunto de test¶
Una vez identificado y entrenado el mejor modelo, se calcula su error al predecir los datos de test.
# Backtest modelo final
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['Demand'],
steps = 24,
metric = mean_absolute_error,
initial_train_size = len(datos[:fin_validacion]),
refit = False,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
fig, ax = plt.subplots(figsize=(8, 3.5))
datos.loc[predicciones.index, 'Demand'].plot(linewidth=2, label='test', ax=ax)
predicciones.plot(linewidth=2, label='predicción', ax=ax)
ax.set_title('Predicción vs demanda real')
ax.legend();
# Error backtest
# ==============================================================================
print(f'Error backtest: {metrica}')
Tras la optimización de lags e hiperparámetros, se ha conseguido reducir el error de predicción de 289.5 a 251.9.
Intervalos de predicción¶
Un intervalo de predicción define el intervalo dentro del cual es de esperar que se encuentre el verdadero valor de $y$ con una determinada probabilidad. Por ejemplo, es de esperar que el intervalo de predicción (1, 99) contenga el verdadero valor de la predicción con un 98% de probabilidad.
Rob J Hyndman y George Athanasopoulos, listan en su libro Forecasting: Principles and Practice mútiples formas de estimar intervalos de predicción, la mayoría los cuales requieren que los resudios (errores) del modelo se distribuyan de forma normal. Cuando no se puede asumir esta propiedad, se puede recurrir a bootstrapping, que solo asume que los residuos no están correlacionados. Este es el método utilizado en la librería skforecast.
# Backtest del conjunto de test con intervalos de predicción
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['Demand'],
steps = 24,
metric = 'mean_absolute_error',
initial_train_size = len(datos.loc[:fin_validacion]),
refit = False,
interval = [10, 90],
n_boot = 500,
in_sample_residuals = True,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
print('Métrica backtesting:', metrica)
predicciones.head(5)
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 3.5))
datos.loc[predicciones.index, 'Demand'].plot(linewidth=2, label='test', ax=ax)
predicciones['pred'].plot(linewidth=2, label='predicción', ax=ax)
ax.set_title('Predicción vs demanda real')
ax.fill_between(
predicciones.index,
predicciones['lower_bound'],
predicciones['upper_bound'],
alpha = 0.3,
color = 'red',
label = 'Intervalo predicción'
)
ax.legend();