Машинное обучение в энергетике, или не только лишь все могут смотреть в завтрашний день

Точное предсказание будущих событий — перспективная и интересная задача во многих сферах: от прогноза погоды до финтеха (котировки акций, курсы валют). Машинное обучение уже сегодня позволяет значительно сократить время и трудозатраты на принятие управленческих решений. 

Наша Data Science команда в НОРБИТ около полугода экспериментировала с использованием различных моделей машинного обучения для решения задач по классификации и регрессии, и по оптимизации бизнес-процессов в сфере b2b. Но когда появилась задача по предсказанию временных рядов, оказалось, что доступных материалов на эту тема в сети недостаточно для разработки быстрого решения. Суть задачи заключалась в том, чтобы с максимальной точностью (средняя абсолютная ошибка в процентах, или MAPE < 3%) предсказывать почасовое потребление целого города на следующие трое суток суток (72 часа) для нужд одной крупной электрогенерирующей компании. Максимальная допустимая ошибка в 3% — это очень высокий показатель по точности прогнозов. До применения машинного обучения ошибка находилась в пределах 3-23%, что напрямую приводило к финансовым убыткам, так как при недостаточной выработке электрической энергии дополнительные мощности необходимо было покупать на оптовом рынке электроэнергии (что сильно дороже собственной выработки), при перепроизводстве — продавать на рынке дешевле, чем могли бы продать городу. Также негативным эффектом отклонения плановой нагрузки от фактической является изменение частоты переменного тока в сети.
Достижению высокой точности прогнозирования мешают много факторов, например, при внезапном повышении температуры люди сразу тянутся за пультами от кондиционеров, при понижении — достают с антресолей обогреватели; во время праздников или больших футбольных матчей включают телевизоры и т.д. Если некоторые события цикличны и потенциально предсказуемы, то другие — совершенно нет. Допустим, внезапно по причине аварии перестал работать сталелитейный цех и все прогнозы превратились в тыкву (в нашем целевом городе обошлось без крупных производственных предприятий). Был довольно интересный пример проблемы при построении прогнозов в энергетике, когда один из всплесков изменений в объеме вырабатываемой электроэнергии был вызван тем, что корабль сел на мель на озере, и капитан договорился с владельцами гидроэлектростанций, чтобы они уменьшили сброс воды. Естественно, что модель машинного обучения не смогла предугадать это событие на основании накопленных исторических данных.

Хорошая новость для нас, что заказчик предусмотрел в техническом задании такие «особенные дни», в которые была допустима средняя ошибка прогноза в 8%.

Все, что предоставил заказчик, — это исторические данные почасового потребления электроэнергии за 3 года и название города.

Первая задача для подготовки построения модели машинного обучения — визуализация и поиск дополнительных факторов, которые могут оказывать влияние на прогноз. Степень такого влияния удобно визуально оценивать с помощью тепловой карты корреляции признаков. Темный квадрат в данном случае означает безусловную обратную зависимость величин (чем больше одна величина, тем меньше другая, и наоборот), белый — прямую:

import pandas as pd import seaborn as sns import matplotlib.pyplot as plt   # df - pandas DataFrame с признаками sns.heatmap(df.corr());

Насколько я знаю, в моделях, которые используют для прогнозов потребления электроэнергии в РФ, используют 2 фактора — история потребления и прогноз температуры.  В зимний и осенний периоды потребление энергии выше — сказывается низкая температура и сравнительно короткий световой день. Фактор «день недели» также оказывает влияние на количество расходуемой энергии — по выходным оно резко падает. Пики потребления приходятся на среду или четверг. Внутри одних суток по будням пиковые значения энергопотребления приходятся на 11-12 часов и постепенно падают с резким обрывом после девяти часов вечера. 

Все как в учебниках: Суточные зимний и летний графики потребления промышленного и культурного центров. Источник

Моделирование

Prophet

Первая идея, как сделать задачу максимально быстро — это взять Prophet от Facebook. У меня уже был опыт его использования, и я его запомнил как «штука, которая просто и быстро работает из коробки». Prophet для работы хочет видеть в pandas-датафрейме колонки «ds» и «y», т.е. дату в формате YYYY-MM-DD HH:MM:SS, и таргет — потребление в час, соответственно (примеры работы есть в документации).

df = pd.read_pickle('full.pkl') pjme = df.loc[:'2018-06-01 23:00'].rename(columns={'hour_value': 'y'})  pjme['ds'] = pjme.index split_date = '2018-03-01 23:00' pjme_train = pjme.loc[pjme.index <= split_date] pjme_test = pjme.loc[pjme.index > split_date]   playoffs = pd.DataFrame({   'holiday': 'playoff',   'ds': df[(df.rest_day == 1)|(df.is_weekend == 1)].index,   'lower_window': 0,   'upper_window': 1, })   model = Prophet(holidays=playoffs) model.fit(pjme_train.reset_index()) forecast = model.predict(df=pjme_test.reset_index())

Предсказания «пророка» выглядели объективными, но ошибка в 7–17% превышала допустимую, поэтому дальнейшие эксперименты я на этом фреймворке не проводил.

SARIMAX

Вторая попытка решения этой задачи была с помощью SARIMAX. Метод довольно затратный по времени подбора коэффициентов, но зато удалось снизить ошибку прогноза по сравнению с Prophet до 6-11%. 

К сожалению, остался только результат для недельных предсказаний, но именно эти прогнозные данные я использовал в качестве дополнительных фич в бустинг-моделях.

Для прогноза нужно подобрать параметры SARIMA(p,d,q)(P,D,Q,s):

• p — порядок авторегрессии (AR), который позволяет добавить прошлые значения временного ряда (можно выразить как «завтра, возможно, пойдет снег, если в последние несколько дней шел снег»);
• d — порядок интегрирования исходных данных (он определяет число предыдущих временных точек, которые нужно вычесть из текущего значения);
• q — порядок скользящего среднего (MA), который позволяет установить погрешность модели в виде линейной комбинации наблюдавшихся ранее значений ошибок;
• P — p, но сезонное;
• D — d, но сезонное;
• Q — q, но сезонное;
• s — размерность сезонности (месяц, квартал и т.д.).

На сезонном разложении по недельной истории можно увидеть тренд и сезонности, которые позволяет отобразить seasonal_decompose:

import statsmodels.api as sm sm.tsa.seasonal_decompose(df_week).plot()

По графикам автокорреляции и частичной автокорреляции выбрал начальные приближения параметров для модели SARIMAX: p=0, P=1, q=2, Q=2.

sm.graphics.tsa.plot_acf(df_week.diff_week[52:].values.squeeze(), lags=60, ax=ax) sm.graphics.tsa.plot_pacf(df_week.diff_week[52:].values.squeeze(), lags=60, ax=ax)

И итоговые подобранные значения (0, 2, 2)x(1, 2, 0, 52). В итоге график прогноза потребления выглядел так:

Бустинг

В этот момент стало очевидно, что без использования дополнительных внешних факторов не удастся достичь необходимой точности. В первую очередь, это погодные факторы: температура, давление, влажность, сила и направление ветра, осадки. 

Первой попыткой была попытка купить архив с погодой в Яндекс. У ребят отличное API и техподдержка, но конкретно по нашему городу в данных было довольно много пропусков. В итоге я купил архив в сервисе openweathermap.org за 10$. Важно отметить, что, несмотря на то, что погода — очень полезный фактор, ошибки в прогнозе очевидно будут сильно портить итоговую MAPE модели. Я встречал информацию, что правильным решением этой проблемы было бы в обучении использовать не фактические исторические значения погоды, а прогнозы погоды за трое суток, но у меня, к сожалению, такой возможности не было.  Также я добавил признаки времени суток, дня недели, праздники, порядковый номер дня в году, а также большое количество временных лагов и средних значений за предыдущие периоды.
Все эти данные после масштабирования (MinMaxScale) и One Hot Encoding (значения каждого категориального признака становятся отдельными столбцами с 0 и 1, таким образом, признак с тремя значениями становится тремя разными столбцами, и единица будет лишь в одном из них) я использовал в небольшом соревновании трех популярных моделей бустинга XGBoost, LightGBM и CatBoost. 

XGBoost

Об XGBoost’е написано много хорошего. Представленный на конференции SIGKDD в 2016 году, он произвел настоящий фурор и до сих пор держится на лидирующих позициях. Интересный материал для тех, кто про него не слышал. Подготовим данные для модели: 

from sklearn.preprocessing import MinMaxScaler, StandardScaler from sklearn.model_selection import train_test_split   def gen_data(     stop_day = '2018-06-01 23:00',      drop_columns = [],     test_size=0.15,     is_dummies=True,     is_target_scale=False):      df = pd.read_pickle('full.pkl')     df = df.loc[:stop_day]     scaler = MinMaxScaler()     target_scaler = StandardScaler()     y = df.hour_total          if is_target_scale:                 y = target_scaler.fit_transform(y.values.reshape(-1, 1)).reshape(1,-1)[0]              X = df.drop(['hour_value', *drop_columns], axis=1)     num_columns = X.select_dtypes(include=['float64']).columns     X[num_columns] = scaler.fit_transform(X[num_columns])          if is_dummies:         X = pd.get_dummies(X)      train_count_hours = len(X) - 72     valid_count_hours = len(X) - int(len(X) * 0.2)     X_test  = X[train_count_hours:]     y_test  = y[train_count_hours:]     X = X[:train_count_hours]     y = y[:train_count_hours]      # Тут можно пробовать разные варианты разбиения на обучающую и тестовую выборки, от этого сильно зависит итоговая точность модели     X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=test_size, random_state=42)     train_dates = X_train.index     val_dates = X_test.index     return X, y, X_train, X_val, y_train, y_val, X_test, y_test, target_scaler

Учу модель, используя лучшие подобранные ранее гиперпараметры:

X, y, X_train, X_val, y_train, y_val, X_test, y_test, _ = gen_data(     stop_day = stop_day,      drop_columns = drop_columns,     test_size=0.1,     is_dump=True,     is_target_scale=False)   params_last = {     'base_score': 0.5,      'booster': 'gbtree',      'colsample_bylevel': 1,      'colsample_bynode': 1,      'colsample_bytree': 0.4,      'gamma': 0,     'max_depth': 2,      'min_child_weight': 5,      'reg_alpha': 0,      'reg_lambda': 1,      'seed': 38,     'subsample': 0.7,      'verbosity': 1,     'learning_rate':0.01 }   reg = xgb.XGBRegressor(**params_last, n_estimators=2000) print(X_train.columns) reg.fit(X_train, y_train,         eval_set=[(X_val, y_val)],         early_stopping_rounds=20,         verbose=False)   y_pred = reg.predict(X_test)

LightGBM

LightGBM — фреймворк от Microsoft, основное преимущество которого состоит в скорости обучения на действительно больших массивах данных. А также в отличие от XGBoost’а LightGBM умеет работать с категориями, использует меньше памяти. Вот тут есть подробнее.

import lightgbm as lgb from lightgbm import LGBMRegressor from sklearn.metrics import mean_squared_error   stop_day = '2018-06-01 23:00' start_day_for_iterate = datetime.strptime(stop_day, '%Y-%m-%d %H:%M') test_size = 0.2  X, y, X_train, X_val, y_train, y_val, X_test, y_test, _ = gen_data(     stop_day = stop_day,      drop_columns = drop_columns,     test_size=test_size,     is_dump=True,     drop_weekend=False,     is_target_scale=False)   model = LGBMRegressor(boosting_type='gbdt',                  num_leaves=83,                  max_depth=-1,                  learning_rate=0.008,                  n_estimators=15000,                  max_bin=255,                  subsample_for_bin=50000,                  min_split_gain=0,                  min_child_weight=3,                 min_child_samples=10,                  subsample=0.3,                  subsample_freq=1,                  colsample_bytree=0.5,                  reg_alpha=0.1,                  reg_lambda=0,                  seed=38,                 silent=False,                  nthread=-1)   history = model.fit(X_train, y_train,              eval_metric='rmse',             eval_set=[(X_val, y_val)],             early_stopping_rounds=40,             verbose = 0)   y_pred = model.predict(X_test, num_iteration=model.best_iteration_)

CatBoost

CatBoost — продвинутая библиотека градиентного бустинга на деревьях решений с открытым исходным кодом" от Яндекса. Выгодное отличие модели состоит в удобстве работы с датасетами, содержащими категориальные признаки — передавать в модель содержащие категории данные можно без преобразования в числа, причем она выявляет закономерности самостоятельно, руками ничего подкручивать не требуется, и качество предсказания остается на высоте.

cat_features=np.where(X.dtypes == 'category')[0] eval_dataset = Pool(X_test, y_test)  model = CatBoostRegressor(learning_rate=0.5,                           eval_metric='RMSE',                            leaf_estimation_iterations=3,                           depth=3)  model.fit(X_train, y_train,           eval_set=(X_val, y_val),           cat_features=cat_features,           plot=True,           verbose=True)  y_pred = model.predict(X_test)

XGBoost vs. LightGBM vs. CatBoost

Чтобы не повторять многочисленные статьи, приведу сравнительную таблицу отсюда.

Для вычисления ошибки по MAPE я брал последний известный месяц (28 дней) и, двигая окно с дискретностью в один час, делал прогноз на следующие 72 часа. А в моем импровизированном соревновании призовые места заняли:

3-е место: мой любимчик — XGBoost — с самым низким качеством прогноза;
2-е место: CatBoost как самый удобный и «все из коробки»;
1-е место: LightGBM как самый быстрый и точный.

Для более точного сравнения моделей я использовал R2 (R-квадрат или коэффициент детерминации, показывающий, как сильно условная дисперсия модели отличается от дисперсии реальных значений) и RMSLE (Root Mean Squared Logarithmic Error, или среднеквадратичная логарифмическая ошибка, которая является, по сути, расстоянием между двумя точками на плоскости — реальным значением и предсказанным). 

Но это все, разумеется, применительно к моей задаче. На других данных в других руках все может быть совсем по-другому.  В следующий раз поделюсь другой нашей собственной интересной историей по прогнозированию оттока сотрудников для задачи планирования подбора в нашей компании.