Pronosticando la producción de energía photo voltaic con collection de tiempo y el modelo XGBoost
La energía photo voltaic es por mucho la fuente de energía más abundante y prácticamente inagotable. Se estima que la tierra recibe anualmente 44 cuatrillones (4.4 x 10¹⁶) de watts de energía photo voltaic. En comparación, una planta eléctrica de gran escala produce aproximadamente 1 billón de watts al año. Se necesitarían 44 millones de estas para igualar al sol. [1]
Además de los retos tecnológicos y económicos, otro issue que influye en la adopción de fuentes renovables para producir energía eléctrica, es la dificultad para predecir su comportamiento. Contar con pronósticos confiables para la producción de energía photo voltaic es esencial para aprovechar más eficientemente este recurso.
En este artículo, usaremos datos de las instalaciones fotovoltaicas en la ciudad de Calgary, Canadá. Haremos un análisis de collection de tiempo, y usaremos el algoritmo XGBoost en Python para pronosticar la producción de energía photo voltaic.
Puedes encontrar el proyecto completo en mi repositorio de GitHub.
Los datos se pueden descargar manualmente del portal de datos abiertos de la ciudad de Calgary.
Alternativamente puedes trabajar con el dataset desde Kaggle aquí.
Explorando los datos
Iniciamos por importar librerías y leer el dataset:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import xgboost as xgbdf = pd.read_csv('./../knowledge/Solar_Energy_Production.csv')
Existen 11 instalaciones fotovoltaicas en la ciudad y observamos que la capacidad instalada fue en aumento desde 2015 hasta finales de 2017:
site_install = df.groupby('title')['installationDate'].distinctive()
site_install = pd.DataFrame(site_install).sort_values(by='installationDate')site_install
Visualizamos la producción complete en GigaWatt-hr de cada instalación:
counts = df.groupby('title')['kWh'].sum()
site_totals = pd.DataFrame(counts).sort_values(by='kWh', ascending=True)site_totals.plot(figsize=(15, 6), sort='barh', legend=False)
plt.title('Calgary Photo voltaic Energy Technology per Website [september 2015 to march 2023]')
plt.xlabel('GigaWatt-hours')
plt.ylabel('Website')
plt.present()
Para el análisis de collection de tiempo necesitamos solo la cantidad de energía (en kWh) y la fecha de producción photo voltaic, en este caso las columnas ’kWh’ y ‘date’ respectivamente. Asignaremos la columna ‘date’ al índice del dataframe, pero antes debemos convertirla de tipo object a datetime utilizando el método to_datetime de la librería Pandas. Además, nos resulta mas útil tener la producción de energía en una fecuencia de acumulado diario:
#convertir a datetime
df['date'] = pd.to_datetime(df['date'])
df_pw = df.drop(columns= ['name', 'id', 'address', 'public_url', 'installationDate', 'uid'])
df_pw = df_pw.set_index('date')#crear acumulado diario
count_date = df_pw.groupby(df_pw.index.date)['kWh'].sum()
pw_clean = pd.DataFrame(count_date)
pw_clean['date'] = pd.to_datetime(pw_clean.index)
pw_clean = pw_clean.set_index('date')
El primer paso en cualquier análisis de collection de tiempo es graficar los datos. La visualización nos permite detectar patrones, lecturas atípicas, variaciones en el tiempo y la relación entre variables:
pw_clean.plot(fashion='-', figsize=(20, 7), lw=1,
title='Calgary Photo voltaic Energy Manufacturing in kWh')
plt.present()
La gráfica nos confirma que la producción de energía photo voltaic se estabilizó a partir del año 2018. Como ya habíamos observado, fue a finales de 2017 cuando se completó la ultima instalación fotovoltaica en la ciudad.
Tomando esto en cuenta, utilizaremos solo los datos a partir de 2018 para el entrenamiento de nuestro modelo.
Breve reseña del modelo XGBoost
XGBoost (Excessive Gradient Boosting) es una implementación del algoritmo Gradient Boosting (GBM) para tareas de aprendizaje supervisado (regresión y clasificación). El algoritmo XGBoost utiliza ensambles de árboles de decisión (Tree ensembles), combinando multiples modelos “débiles” para generar un modelo predictivo con mejor desempeño. La técnica de boosting es usada como método de ensamble y gradient descent para optimizar el proceso. XGBoost nos ofrece gran flexibilidad en cuanto al ajuste de hiperparámentros para controlar la complejidad y evitar el sobreentrenamiento (overfitting). Además, proporciona soporte para trabajar con valores faltantes y variables categóricas.
Para nuestro análisis solo estaremos ajustando los siguientes hiperparámetros:
booster: El tipo de modelo a utilizar, por defectogbtreeusa árboles de decisión.n_estimators: El número de modelos que utilizará el ensamble. Un número muy bajo puede causar underfitting, mientras que un valor muy alto puede llevar a un modelo sobreentrenado (overfitting).early_stopping_rounds: El número de iteraciones para detener el entrenamiento automáticamente una vez que el rating de validación deja de mejorar. Es común combinar un valor alto den_estimatorsy dejar queearly_stopping_roundsencuentre el numero óptimo de iteraciones.goal: Especifica la tarea de aprendizaje y el objetivo de optimización. Para nuestro análisis de regresión y error cuadrático, especificamos'reg:squarederror'.reg_lambda: Es el parámetro de regularización equivalente a Ridge (L2) en regresiones múltiples. En normal este parámetro penaliza uniformemente todos los pesos (weights) del modelo. Incrementar el valor hace un modelo más conservador.max_depth: La profundidad o número máximo de splits de los árboles de decisión que usará el algoritmo. Con mayor profundidad podemos obtener mejores resultados, pero puede resultar en overfitting.eta: El tamaño de paso en cada iteración o taza de entrenamiento. Un valor mayor llegará más rápido al mínimo de la función objetivo, pero puede llegar a “saltarse” el valor óptimo. Un valor pequeño tomará mas tiempo o incluso nunca llegar al mínimo de la función objetivo.
Puedes consultar la documentación completa de hiperparámetros aquí
Creación de Atributos (characteristic engineering)
De momento nuestro dataset solo contiene la variable de salida (producción de energía en kWh):
El siguiente paso es agregar características que ayuden a explicar el comportamiento de esta variable. Para ello, aprovechamos los métodos de la librería Pandas para datos tipo datetime. Definimos una función para crear variables basadas en el índice de la serie de tiempo:
def create_attributes(df):
df = df.copy()
df['day'] = df.index.day
df['dayofweek'] = df.index.dayofweek
df['month'] = df.index.month
df['quarter'] = df.index.quarter
df['year'] = df.index.12 months
df['dayofyear'] = df.index.dayofyear
return dfpw_clean = create_attributes(pw_clean)
Entrenando el modelo
Antes del entrenamiento dividimos el set de prueba y validación. Recordamos que para collection de tiempo no es adecuado dividir aleatoriamente el set para validación cruzada (cross validation), debido a la correlación entre observaciones cercanas en el tiempo. Para ello, usamos la clase TimeSeriesSplit de la librería Scikitlearn. Dividimos los datos en 4 partidas para validación cruzada y la métrica de evaluación será el error cuadrático medio (RMSE):
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_squared_error #dividir en 4 folds
ts_cv = TimeSeriesSplit(n_splits=4)
fold = 0
preds = []
scores = []
for train_idx, val_idx in ts_cv.cut up(pw_clean):
prepare = pw_clean.iloc[train_idx]
take a look at = pw_clean.iloc[val_idx]
#agregar atributos al set de entrenamiento y prueba
prepare = create_attributes(prepare)
take a look at = create_attributes(take a look at)
options = ['day','dayofweek','month','quarter','year','dayofyear']
goal = ['kWh']
#dividir atributos y variable de salida
X_train = prepare[features]
y_train = prepare[target]
X_test = take a look at[features]
y_test = take a look at[target]
#crear instancia del regresor
xgb_reg = xgb.XGBRegressor(booster='gbtree',
seed=42,
n_estimators=1000,
early_stopping_rounds=50,
goal='reg:squarederror',
reg_lambda=0.001,
max_depth=5,
eta=0.01)
#entrenar modelo
xgb_reg.match(X_train, y_train,
eval_set=[(X_train, y_train), (X_test, y_test)],
verbose=100)
#predicciones y evaluacion
y_pred = xgb_reg.predict(X_test)
preds.append(y_pred)
rating = np.sqrt(mean_squared_error(y_test, y_pred))
scores.append(rating)
Observamos el rating obtenido en cada partida y el promedio ultimate:
print('Fold scores:', scores)
print('Avg. Rating:', np.imply(scores))
Fold scores: [2598.514863235386, 2435.9595990180537, 2117.653168183066, 1949.2793308876999]
Avg. Rating: 2275.3517403310516
Graficar el proceso de validación cruzada nos da una perspectiva visible de como se entrenó nuestro modelo. Generamos la visualización con el siguiente código:
fig, axs = plt.subplots(4, 1, figsize=(15, 15), sharex=True)fold = 0
for train_idx, val_idx in ts_cv.cut up(pw_clean):
prepare = pw_clean.iloc[train_idx]
take a look at = pw_clean.iloc[val_idx]
prepare['kWh'].plot(ax=axs[fold], lw=1,
label='Coaching Set',
title=f'Prepare/Check Break up Fold {fold}')
take a look at['kWh'].plot(ax=axs[fold], lw=1,
label='Check Set')
axs[fold].axvline(take a look at.index.min(), coloration='black', ls='--')
fold += 1
plt.present()
Prediciendo el futuro
Primero re-entrenamos el modelo con todos los datos. Deseamos aprovechar toda la información disponible ya que estaremos haciendo predicciones a partir de la última fecha del dataset:
#crear dataframe completo
pw_clean = create_attributes(pw_clean)options = ['day','dayofweek','month','quarter','year','dayofyear']
goal = ['kWh']
X_full = pw_clean[features]
y_full = pw_clean[target]
#crear instancia del regresor
xgb_regf = xgb.XGBRegressor(booster='gbtree',
seed=42,
n_estimators=1000,
early_stopping_rounds=50,
goal='reg:squarederror',
reg_lambda=0.001,
max_depth=5,
eta=0.01)
#entrenar modelo
xgb_regf.match(X_full, y_full,
eval_set=[(X_full, y_full)],
verbose=100)
A continuación, creamos un nuevo dataframe en el cual agregamos las fechas que deseamos pronosticar. Posteriormente, concatenamos el nuevo dataframe con el dataframe authentic, agregando una columna de tipo booleano para identificar las fechas en el futuro.
Probamos con 6 meses en el futuro. El método date_range de Pandas nos devuelve un rango de fechas en una frecuencia específica, para nuestro análisis será por día:
#rango de marzo a septiembre 2023
pred_dates = pd.date_range('2023-03-16','2023-09-17', freq='D')
preds_df = pd.DataFrame(index=pred_dates)#Crear columna
preds_df['Future'] = True
pw_clean['Future'] = False
#Concatenar dataframes
pred_pw = pd.concat([pw_clean.loc[pw_clean.index >= '01-01-2018'], preds_df])
Agregamos atributos con nuestra función create_attributes y finalmente producimos un dataframe con las fechas futuras para hacer el pronóstico (usaremos el dataframe “futuro” como enter del modelo y el dataframe “completo” para visualizar la serie de tiempo ultimate):
#agregar atributos
pred_pw = pred_pw.copy()
pred_pw = create_attributes(pred_pw)#seleccionar fechas
future_pred_pw = pred_pw.question('Future').copy()
future_pred_pw.head()
Pronóstico
Asignamos los valores pronosticados a una nueva columna ‘prediction’ y finalmente creamos una gráfica con los datos históricos y el pronóstico obtenido:
#predicciones del modelo
future_pred_pw['prediction'] = xgb_regf.predict(future_pred_pw[features])#graficar
ax = pw_clean['kWh'].loc[pw_clean.index >= '01-01-2018']
.plot(figsize=(20, 6), lw=1.5, title='Calgary Photo voltaic Energy Manufacturing in kWh')
future_pred_pw['prediction'].plot(fashion='-', lw=1.5)
ax.axvline('2023-03-16', coloration='black', ls='--', lw=1.5)
plt.legend(['Historic Data','Predictions'], fontsize=14)
plt.present()
La gráfica nos muestra que el modelo logra capturar el comportamiento de la producción de energía photo voltaic, siguiendo la variación estacional observada en los datos históricos.
Siguientes pasos
Para obtener mejores resultados podemos considerar lo siguiente:
- Obtener y agregar datos del estado del tiempo como un issue para la producción de energía photo voltaic.
- Crear atributos “lag” al igual que en modelos auto-regresivos.
- Ajuste de hiperparámetros mas robusto (experimentar con Grid search).
