Хабрахабр

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

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

Наша Data Science команда в НОРБИТ около полугода экспериментировала с использованием различных моделей машинного обучения для решения задач по классификации и регрессии, и по оптимизации бизнес-процессов в сфере b2b. Но когда появилась задача по предсказанию временных рядов, оказалось, что доступных материалов на эту тема в сети недостаточно для разработки быстрого решения.

Суть задачи заключалась в том, чтобы с максимальной точностью (средняя абсолютная ошибка в процентах, или MAPE

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

Достижению высокой точности прогнозирования мешают много факторов, например, при внезапном повышении температуры люди сразу тянутся за пультами от кондиционеров, при понижении — достают с антресолей обогреватели; во время праздников или больших футбольных матчей включают телевизоры и т.д. Если некоторые события цикличны и потенциально предсказуемы, то другие — совершенно нет. Допустим, внезапно по причине аварии перестал работать сталелитейный цех и все прогнозы превратились в тыкву (в нашем целевом городе обошлось без крупных производственных предприятий). Был довольно интересный пример проблемы при построении прогнозов в энергетике, когда один из всплесков изменений в объеме вырабатываемой электроэнергии был вызван тем, что корабль сел на мель на озере, и капитан договорился с владельцами гидроэлектростанций, чтобы они уменьшили сброс воды. Естественно, что модель машинного обучения не смогла предугадать это событие на основании накопленных исторических данных.

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

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

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

import pandas as pdimport seaborn as snsimport 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.indexsplit_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 smsm.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, StandardScalerfrom 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 lgbfrom lightgbm import LGBMRegressorfrom 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, или среднеквадратичная логарифмическая ошибка, которая является, по сути, расстоянием между двумя точками на плоскости — реальным значением и предсказанным). 

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

В следующий раз поделюсь другой нашей собственной интересной историей по прогнозированию оттока сотрудников для задачи планирования подбора в нашей компании.

Кстати, наша команда растёт, и мы ищем людей!
Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть