Прогнозирование нагрузки на энергосистему с визуализацией

В этом руководстве мы покажем, как использовать NeuralProphet для прогнозирования нагрузки энергосистемы и визуализации результатов в интерактивном инструменте. Прогноз основан на наборе данных RE-Europe, который моделирует континентальную электроэнергетическую систему Европы, включая потребности в энергии и притоки возобновляемой энергии за период с 2012 по 2014 год. Общий набор данных содержит 1494 узла, которые мы сократим до 5 узлов для целей данного примера.

# Импорт библиотек для работы с данными:
import pandas as pd  # Библиотека pandas для анализа и обработки данных
import numpy as np  # Библиотека numpy для численных вычислений

# Импорт библиотек NeuralProphet для прогнозирования временных рядов:
from neuralprophet import NeuralProphet, set_log_level  # Класс NeuralProphet и функция set_log_level

# Импорт библиотек Plotly для создания визуализаций:
import plotly.graph_objects as go  # Модуль go из Plotly для создания статических элементов визуализации
import plotly.express as px  # Модуль px из Plotly для создания выразительных диаграмм с минимальным кодом
from plotly.subplots import make_subplots  # Функция make_subplots из Plotly для создания подсюжетов

# Импорт библиотек ipywidgets для создания интерактивных элементов:
import ipywidgets as widgets  # Модуль widgets из ipywidgets для создания интерактивных элементов в Jupyter Notebook
from ipywidgets import interact_manual  # Функция interact_manual из ipywidgets для создания интерактивных элементов управления

# Установка режима визуализации по умолчанию:
plotting_backend = "plotly-static"  # Установить режим визуализации по умолчанию на "plotly-static" для статических изображений

Набор данных включает для каждого автобуса и каждого часа года следующую информацию:

• нагрузка (МВт)

• фактическое производство солнечной энергии (МВт)

• прогнозируемое производство солнечной энергии (МВт)

The forecasts are at hour 00 of every day for the next 24 hours.

# Загрузка данных
df = pd.read_csv(

  "https://raw.githubusercontent.com/ourownstory/neuralprophet-data/main/datasets/multivariate/ER_Europe_subset_10nodes.csv",
  # Чтение данных из CSV файла по ссылке
)

df["ID"] = df["ID"].astype(str)  # Преобразование столбца "ID" в строковый тип
IDs = df["ID"].unique()  # Получение уникальных значений из столбца "ID"
df["ds"] = pd.to_datetime(df["ds"])  # Преобразование столбца "ds" в формат даты и времени

# Использование данных за один год для ускорения обучения
df = df[df["ds"] > "2014-01-01"]  # Фильтрация данных с датой позже 01.01.2014
# Вывод первых 5 строк измененного dataframe
df.head()

	ds	ID	y	solar	solar_fcs
17545	2014-01-01 01:00:00	1	72.3142	0.0	0.0
17546	2014-01-01 02:00:00	1	67.5617	0.0	0.0
17547	2014-01-01 03:00:00	1	63.0677	0.0	0.0
17548	2014-01-01 04:00:00	1	59.9401	0.0	0.0
17549	2014-01-01 05:00:00	1	58.9296	0.0	0.0

# Построение графика первого месяца данных

fig = px.line(  # Создать линейный график с помощью plotly.express
    df[df["ds"] < "2014-02-01"],  # Фильтрация данных за январь 2014
    x="ds",  # Значения по оси X - дата
    y="y",  # Значения по оси Y - энергопотребление
    color="ID",  # Цвет линии определяется по столбцу "ID"
)

# Настройка оформления графика
fig.update_layout(
    xaxis_title="Дата",  # Заголовок оси X - Дата
    yaxis_title="Потребление энергии",  # Заголовок оси Y - Потребление энергии
    legend_title="ID",  # Заголовок легенды - ID
    height=600,  # Высота графика
    width=1000,  # Ширина графика
)

# Отображение графика в зависимости от режима визуализации
if plotting_backend == "plotly-static":
    fig.show("svg")  # Показать график в формате SVG (статическое изображение)

1. Подготовка данных

Чтобы лучше улавливать разные времена года и дни недели, мы добавляем условную сезонность. Для каждого времени года (лето, зима, осень и весна) и будни/выходные моделируется отдельная сезонность.

# Условная сезонность: 4 сезона для годовой и недельной сезонности

df["summer"] = 0  # Создать столбец "summer" (лето) с начальными нулями
df.loc[df["ds"].dt.month.isin([6, 7, 8]), "summer"] = 1  # Лето - июнь, июль, август

df["winter"] = 0  # Создать столбец "winter" (зима) с начальными нулями
df.loc[df["ds"].dt.month.isin([12, 1, 2]), "winter"] = 1  # Зима - декабрь, январь, февраль

df["fall"] = 0  # Создать столбец "fall" (осень) с начальными нулями
df.loc[df["ds"].dt.month.isin([9, 10, 11]), "fall"] = 1  # Осень - сентябрь, октябрь, ноябрь

df["spring"] = 0  # Создать столбец "spring" (весна) с начальными нулями
df.loc[df["ds"].dt.month.isin([3, 4, 5]), "spring"] = 1  # Весна - март, апрель, май


# Условная сезонность: 4 сезона, различие будний день/выходные для каждого сезона, для суточной сезонности

df["summer_weekday"] = 0  # Создать столбец "summer_weekday" (будний день летом) с начальными нулями
df.loc[(df["ds"].dt.month.isin([6, 7, 8])) & (df["ds"].dt.dayofweek.isin([0, 1, 2, 3, 4])), "summer_weekday"] = 1
# Лето будний день - пн, вт, ср, чт, пт (dt.dayofweek - 0 - понедельник)

df["summer_weekend"] = 0  # Создать столбец "summer_weekend" (выходные летом) с начальными нулями
df.loc[(df["ds"].dt.month.isin([6, 7, 8])) & (df["ds"].dt.dayofweek.isin([5, 6])), "summer_weekend"] = 1
# Лето выходные - сб, вс (dt.dayofweek - 5 - суббота, 6 - воскресенье)


df["winter_weekday"] = 0  # Аналогично для зимы (winter*)
df.loc[(df["ds"].dt.month.isin([12, 1, 2])) & (df["ds"].dt.dayofweek.isin([0, 1, 2, 3, 4])), "winter_weekday"] = 1
df["winter_weekend"] = 0
df.loc[(df["ds"].dt.month.isin([12, 1, 2])) & (df["ds"].dt.dayofweek.isin([5, 6])), "winter_weekend"] = 1


df["spring_weekday"] = 0  # Аналогично для весны (spring*)
df.loc[(df["ds"].dt.month.isin([3, 4, 5])) & (df["ds"].dt.dayofweek.isin([0, 1, 2, 3, 4])), "spring_weekday"] = 1
df["spring_weekend"] = 0
df.loc[(df["ds"].dt.month.isin([3, 4, 5])) & (df["ds"].dt.dayofweek.isin([5, 6])), "spring_weekend"] = 1


df["fall_weekday"] = 0  # Аналогично для осени (fall*)
df.loc[(df["ds"].dt.month.isin([9, 10, 11])) & (df["ds"].dt.dayofweek.isin([0, 1, 2, 3, 4])), "fall_weekday"] = 1
df["fall_weekend"] = 0
df.loc[(df["ds"].dt.month.isin([9, 10, 11])) & (df["ds"].dt.dayofweek.isin([5, 6])), "fall_weekend"] = 1


# Вывод первых 5 строк измененного dataframe
df.head()

	ds	ID	y	solar	solar_fcs	summer	winter	fall	spring	summer_weekday	summer_weekend	winter_weekday	winter_weekend	spring_weekday	spring_weekend	fall_weekday	fall_weekend
17545	2014-01-01 01:00:00	1	72.3142	0.0	0.0	0	1	0	0	0	0	1	0	0	0	0	0
17546	2014-01-01 02:00:00	1	67.5617	0.0	0.0	0	1	0	0	0	0	1	0	0	0	0	0
17547	2014-01-01 03:00:00	1	63.0677	0.0	0.0	0	1	0	0	0	0	1	0	0	0	0	0
17548	2014-01-01 04:00:00	1	59.9401	0.0	0.0	0	1	0	0	0	0	1	0	0	0	0	0
17549	2014-01-01 05:00:00	1	58.9296	0.0	0.0	0	1	0	0	0	0	1	0	0	0	0	0

2. Определение модели

Лучшие параметры для модели можно получить с помощью настройки гиперпараметров. Установив квантили, моделируется интервал прогнозирования, а не точный прогноз.

# Модель и прогнозирование

quantiles = [0.015, 0.985]  # Установить квантили для прогнозных интервалов (1.5% и 98.5%)

# Параметры модели
params = {
    "n_lags": 7 * 24,  # Количество прошлых наблюдений для прогноза (7 дней * 24 часа)
    "n_forecasts": 24,  # Количество прогнозируемых шагов (24 часа)
    "n_changepoints": 20,  # Количество точек изменения тренда
    "learning_rate": 0.01,  # Скорость обучения
    "ar_layers": [32, 16, 16, 32],  # Слои авторегрессии с шириной 32, 16, 16, 32 единиц
    "yearly_seasonality": 10,  # Учитывать сезонность с периодом 10 (вероятно, подразумевается год)
    "weekly_seasonality": False,  # Игнорировать недельную сезонность
    "daily_seasonality": False,  # Игнорировать суточную сезонность
    "epochs": 40,  # Количество эпох обучения
    "batch_size": 1024,  # Размер пакета данных для обучения
    "quantiles": quantiles,  # Квантили для прогнозных интервалов
}

# Создать модель NeuralProphet с указанными параметрами
m = NeuralProphet(**params)

# Установить режим визуализации для модели
m.set_plotting_backend(plotting_backend)

# Отключить вывод сообщений ниже уровня ERROR
set_log_level("ERROR")

2.1 Запаздывающие и будущие регрессоры

С дополнительной информацией наша модель может стать лучше. Мы используем данные о фактическом производстве солнечной энергии за последние 10 дней в качестве отстающего регрессора и прогноз солнечной активности в качестве будущего регрессора.

# Добавление запаздывающего регрессора
m.add_lagged_regressor(names="solar", n_lags=10 * 24)  # Добавить солнечную активность за 10 дней (24 часа/день)
# Добавление регрессора для будущих значений
m.add_future_regressor(name="solar_fcs")  # Добавить прогноз солнечной активности ("solar_fcs")

<neuralprophet.forecaster.NeuralProphet at 0x12cc57b20>

2.2 Условные сезонности

m.add_seasonality(name="summer_weekly", period=7, fourier_order=14, condition_name="summer")
m.add_seasonality(name="winter_weekly", period=7, fourier_order=14, condition_name="winter")
m.add_seasonality(name="spring_weekly", period=7, fourier_order=14, condition_name="spring")
m.add_seasonality(name="fall_weekly", period=7, fourier_order=14, condition_name="fall")

m.add_seasonality(name="summer_weekday", period=1, fourier_order=6, condition_name="summer_weekday")
m.add_seasonality(name="winter_weekday", period=1, fourier_order=6, condition_name="winter_weekday")
m.add_seasonality(name="spring_weekday", period=1, fourier_order=6, condition_name="spring_weekday")
m.add_seasonality(name="fall_weekday", period=1, fourier_order=6, condition_name="fall_weekday")

m.add_seasonality(name="summer_weekend", period=1, fourier_order=6, condition_name="summer_weekend")
m.add_seasonality(name="winter_weekend", period=1, fourier_order=6, condition_name="winter_weekend")
m.add_seasonality(name="spring_weekend", period=1, fourier_order=6, condition_name="spring_weekend")
m.add_seasonality(name="fall_weekend", period=1, fourier_order=6, condition_name="fall_weekend")

<neuralprophet.forecaster.NeuralProphet at 0x12cc57b20>

3. Обучение модели

Для последующего тестирования нашей модели данные разделяются на обучающие и тестовые. Тестовые данные составляют последние 5% от всего массива данных. Модель обучается на обучающих данных.

# Разделение данных на обучающую и тестовую выборки
df_train, df_test = m.split_df(df, valid_p=0.05, local_split=True)  # Разделить данные на 95% обучающих и 5% тестовых

print(f"Размер обучающей выборки: {df_train.shape}")  # Вывод формы обучающей выборки
print(f"Размер тестовой выборки: {df_test.shape}")   # Вывод формы тестовой выборки

Train shape: (41545, 17)
Test shape: (3090, 17)

# Обучение модели
metrics_fit = m.fit(df_train, freq="H", metrics=True)  # Обучить модель на обучающей выборке (freq - частота данных - часы)

# Прогнозирование
forecast = m.predict(df)  # Сделать прогнозы на основе всего dataframe (df)

# Извлечение столбцов с сезонностью и трендом и преобразование в числовой формат
season_cols = [col for col in forecast.columns if "season" in col or "trend" in col]  # Поиск столбцов с "season" или "trend" в названии

# Преобразование значений в столбцах season_cols в числовой формат
forecast[season_cols] = forecast[season_cols].apply(pd.to_numeric)

4. Визуализация

4.1 Модельные участки

NeuralProphet включает функции построения графиков для изучения прогнозируемых данных. Для обзора функций построения графиков см. учебное пособие по визуализации

m.highlight_nth_step_ahead_of_each_forecast(1).plot(forecast.iloc[-10 * 24 :])

m.plot_parameters(components=["trend", "trend_rate_change", "lagged_regressors"])

4.2 Интерактивный инструмент

Давайте изучим данные более интерактивно. Выбрав дату, будет показан прогноз на следующие 24 часа, что соответствует стандартам индустрии.

def get_day_ahead_format(date, ID, column_name):
  """
  Функция извлекает прогноз на следующие сутки для указанной даты, ID узла и названия столбца.

  Args:
      date (pandas.Timestamp): Дата, для которой требуется прогноз.
      ID (str): ID узла, для которого требуется прогноз.
      column_name (str): Название столбца с прогнозируемой величиной (например, "y" или "yhat").

  Returns:
      pandas.DataFrame: Таблица с прогнозами на следующие сутки.

  Raises:
      ValueError: Ошибка, если запрашиваемая дата отсутствует в прогнозе.
  """

  # Создать пустой DataFrame для хранения значений прогноза
  values = pd.DataFrame(columns=["ds", column_name])

  # Отфильтровать прогнозы по ID узла
  forecast_ID = forecast[forecast["ID"] == ID]

  # Цикл по 24 часам для формирования прогноза на следующие сутки
  for i in range(0, 24):
    # Формирование названия столбца с прогнозом для i-го часа
    if "%" in column_name:
      step_name = f"yhat{i+1} {column_name}"  # Например, yhat10 1.5% (прогноз с 10-го часа, 1.5% квантиль)
    else:
      step_name = f"{column_name}{i+1}"  # Например, yhat10 (прогноз на 10-й час)

    # Дата следующего часа
    ds = date + pd.Timedelta(hours=i + 1)

    # Проверка наличия даты в прогнозах для данного ID
    if ds not in forecast_ID["ds"].values:
      raise ValueError(
          f'Дата {ds} отсутствует в прогнозе. Выберите дату в диапазоне от {forecast_ID["ds"].min()} до {forecast_ID["ds"].max()}'
      )

    # Получение значения прогноза для i-го часа
    step_value = forecast_ID[forecast_ID["ds"] == ds][step_name].values[0]

    # Добавление строки с прогнозом для i-го часа в DataFrame
    values = pd.concat([values, pd.DataFrame({"ds": ds, column_name: step_value}, index=[0])], ignore_index=True)

  return values

@interact_manual  # с этой строкой кода прогноз становится интерактивным, позволяя вам выбирать идентификатор и дату
def show_forecast_uncertainty(ID=IDs, date=widgets.DatePicker()):
    start = pd.to_datetime(date)
    end = start + pd.DateOffset(days=1)

    # получаем прогноз
    np_yhats = get_day_ahead_format(start, ID, "yhat")
    np_yhats_lower = get_day_ahead_format(start, ID, "1.5%")
    np_yhats_upper = get_day_ahead_format(start, ID, "98.5%")
    np_ar = get_day_ahead_format(start, ID, "ar")

    # получаем фактические значения, температуру и прогноз температуры
    forecast_window = forecast[forecast["ID"] == ID]
    forecast_window = forecast_window[(forecast_window["ds"] >= start) & (forecast_window["ds"] < end)]
    df_window = df[df["ID"] == ID]
    df_window = df_window[(df_window["ds"] >= start) & (df_window["ds"] < end)]

    # получаем компоненты
    seasonalities_to_plot = []
    for season in season_cols:
        if np.abs(forecast_window[season].sum()) > 0:
            seasonalities_to_plot.append(season)

    fig = make_subplots(
        rows=3,
        cols=1,
        shared_xaxes=True,
        subplot_titles=("Фактические значения и прогноз", "Солнечная энергия", "Компоненты прогноза"),
    )

    # рисуем фактические значения
    fig.add_trace(
        go.Scatter(
            x=forecast_window["ds"], y=forecast_window["y"], name="Фактические значения", mode="markers", marker=dict(symbol="x")
        ),
        row=1,
        col=1,
    )

    # рисуем прогноз и неопределенность
    fig.add_trace(
        go.Scatter(
            x=np_yhats_lower["ds"],
            y=np_yhats_lower["1.5%"],
            fill=None,
            mode="lines",
            line_color="#A6B168",
            name="1.5% квантиль",
        ),
        row=1,
        col=1,
    )
    fig.add_trace(
        go.Scatter(
            x=np_yhats_upper["ds"],
            y=np_yhats_upper["98.5%"],
            fill="tonexty",
            mode="lines",
            line_color="#A6B168",
            name="98.5% квантиль",
        ),
        row=1,
        col=1,
    )
    fig.add_trace(go.Scatter(x=np_yhats["ds"], y=np_yhats["yhat"], name="Прогноз", line_color="#7A863B"), row=1, col=1)

    # рисуем солнечную энергию
    fig.add_trace(
        go.Scatter(x=df_window["ds"], y=df_window["solar"], mode="lines", name="Фактическое производство солнечной энергии"), row=2, col=1
    )
    fig.add_trace(
        go.Scatter(x=df_window["ds"], y=df_window["solar_fcs"], mode="lines", name="Прогноз производства солнечной энергии"),
        row=2,
        col=1,
    )

    # рисуем различные сезоны из seasonalities_to_plot
    for season in seasonalities_to_plot:
        fig.add_trace(go.Scatter(x=forecast_window["ds"], y=forecast_window[season], name=season), row=3, col=1)
    fig.add_trace(go.Scatter(x=np_ar["ds"], y=np_ar["ar"], name="Авторегрессия"), row=3, col=1)

    fig.update_layout(
        title=f"Прогноз для ID {ID} на {date}",
        yaxis=dict(title="Нагрузка (МВт)"),
        width=1200,
        height=800,
    )

    if plotting_backend == "plotly-static":
        fig.show("svg")
    else:
        fig.show()