Pembelajaran mesin di sektor energi, atau tidak hanya semua orang bisa menonton besok

Prediksi akurat tentang peristiwa di masa depan adalah tugas yang menjanjikan dan menarik di banyak bidang: mulai dari prakiraan cuaca hingga fintech (harga saham, nilai tukar). Pembelajaran mesin hari ini dapat secara signifikan mengurangi waktu dan tenaga yang terlibat dalam pengambilan keputusan manajemen.  Selama sekitar setengah tahun,

tim Ilmu Data kami di NORBIT bereksperimen menggunakan berbagai model pembelajaran mesin untuk menyelesaikan masalah klasifikasi dan regresi, dan untuk mengoptimalkan proses bisnis b2b. Tetapi ketika tugas memprediksi deret waktu muncul, ternyata materi yang tersedia tentang topik ini di jaringan tidak cukup untuk mengembangkan solusi cepat.


Inti dari tugas ini adalah memperkirakan dengan akurasi maksimum (kesalahan absolut rata-rata dalam persen, atau MAPE <3%) konsumsi per jam seluruh kota selama tiga hari ke depan (72 jam) untuk kebutuhan satu perusahaan pembangkit listrik besar.

Kesalahan maksimum yang diijinkan sebesar 3% adalah indikator yang sangat tinggi untuk akurasi perkiraan. Sebelum menggunakan pembelajaran mesin, kesalahannya berada di kisaran 3-23%, yang secara langsung menyebabkan kerugian finansial, karena dengan pembangkit energi listrik yang tidak mencukupi, kapasitas tambahan harus dibeli di pasar grosir listrik (yang jauh lebih mahal daripada yang dihasilkan sendiri), dengan kelebihan produksi - dijual di pasar lebih murah daripada yang bisa mereka jual ke kota. Juga, efek negatif dari penyimpangan beban yang direncanakan dari yang sebenarnya adalah perubahan frekuensi arus bolak-balik dalam jaringan.

Banyak faktor yang mengganggu pencapaian akurasi perkiraan yang tinggi, misalnya, dengan kenaikan suhu yang tiba-tiba, orang-orang segera meraih konsol dari pendingin udara, sementara menurunkan, mereka mendapatkan pemanas dari mezzanine; selama liburan atau pertandingan sepak bola besar termasuk TV, dll. Sementara beberapa peristiwa bersifat siklis dan berpotensi dapat diprediksi, yang lain benar-benar tidak ada. Misalkan, tiba-tiba, karena kecelakaan, toko baja berhenti bekerja dan semua ramalan berubah menjadi labu (di kota target kami tidak ada perusahaan produksi besar). Ada contoh yang agak menarik dari masalah dalam membuat perkiraan di sektor energi, ketika salah satu lonjakan perubahan dalam jumlah listrik yang dihasilkan disebabkan oleh kapal yang terdampar di danau, dan kapten setuju dengan pemilik pembangkit listrik tenaga air untuk mengurangi debit air.Secara alami, model pembelajaran mesin tidak dapat memprediksi peristiwa ini berdasarkan data historis.

Kabar baiknya bagi kami adalah bahwa pelanggan memberikan "hari istimewa" seperti itu dalam kerangka referensi yang memungkinkan kesalahan perkiraan rata-rata 8%.

Semua yang diberikan pelanggan adalah data historis konsumsi listrik per jam selama 3 tahun dan nama kota.


Tugas pertama untuk mempersiapkan konstruksi model pembelajaran mesin adalah visualisasi dan mencari faktor-faktor tambahan yang dapat mempengaruhi perkiraan. Tingkat pengaruh tersebut dinilai secara visual dengan menggunakan peta panas korelasi atribut. Kotak gelap dalam hal ini berarti hubungan terbalik tak bersyarat dari jumlah (semakin besar satu nilai, semakin kecil yang lain, dan sebaliknya), putih - langsung:

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt


# df - pandas DataFrame  
sns.heatmap(df.corr());


Sejauh yang saya tahu, dalam model yang digunakan untuk memperkirakan konsumsi listrik di Federasi Rusia, 2 faktor digunakan - sejarah konsumsi dan perkiraan suhu. 


Di musim dingin dan musim gugur, konsumsi energi lebih tinggi - suhu rendah dan siang hari yang relatif singkat. Faktor "hari dalam seminggu" juga mempengaruhi jumlah energi yang dihabiskan - pada akhir pekan turun tajam. Puncak konsumsi terjadi pada hari Rabu atau Kamis.





Di dalam satu hari pada hari kerja, nilai puncak konsumsi energi terjadi pada 11-12 jam dan secara bertahap turun dengan penurunan tajam setelah jam sembilan malam. 


Semuanya seperti dalam buku teks:

Jadwal konsumsi harian musim dingin dan musim panas harian untuk pusat-pusat industri dan budaya. Sumber

Pemodelan


Nabi


Gagasan pertama bagaimana membuat tugas secepat mungkin adalah mengambil Nabi dari Facebook. Saya sudah memiliki pengalaman menggunakannya, dan saya mengingatnya sebagai "sesuatu yang sederhana dan cepat keluar dari kotak." Nabi ingin melihat kolom "ds" dan "y" di bingkai data panda, yaitu tanggal dalam format YYYY-MM-DD HH: MM: SS, dan targetnya adalah konsumsi per jam, masing-masing (contoh pekerjaan ada dalam dokumentasi ).

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())



Prediksi "nabi" tampak objektif, tetapi kesalahannya 7-17% lebih tinggi dari yang diizinkan, oleh karena itu, saya tidak melakukan percobaan lebih lanjut pada kerangka ini.

Sarimaax


Upaya kedua untuk mengatasi masalah ini adalah menggunakan SARIMAX. Metode ini cukup memakan waktu untuk memilih koefisien, tetapi dimungkinkan untuk mengurangi kesalahan perkiraan dibandingkan dengan Nabi menjadi 6-11%. 

Sayangnya, hanya hasil untuk prediksi mingguan yang tersisa, tetapi data prediksi inilah yang saya gunakan sebagai fitur tambahan dalam model peningkatan.

Untuk perkiraan, Anda perlu memilih parameter SARIMA (p, d, q) (P, D, Q, s):

  • p adalah urutan autoregresif (AR), yang memungkinkan Anda untuk menambahkan nilai masa lalu dari deret waktu (dapat dinyatakan sebagai "besok, mungkin akan turun salju jika salju turun dalam beberapa hari terakhir");
  • d adalah urutan integrasi data sumber (ini menentukan jumlah titik waktu sebelumnya yang akan dikurangkan dari nilai saat ini);
  • q adalah urutan moving average (MA), yang memungkinkan Anda untuk mengatur kesalahan model dalam bentuk kombinasi linear dari nilai-nilai kesalahan yang diamati sebelumnya;
  • P - p, tetapi musiman;
  • D - d, tetapi musiman;
  • Q - q, tapi musiman;
  • s adalah dimensi musiman (bulan, kuartal, dll.).

Pada ekspansi musiman dari sejarah mingguan, Anda dapat melihat tren dan musiman, yang memungkinkan musiman_dekomposisi ditampilkan:

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


Menurut grafik autokorelasi dan autokorelasi parsial, ia memilih perkiraan awal parameter untuk model 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)


Dan nilai akhir yang dipilih adalah (0, 2, 2) x (1, 2, 0, 52). Akibatnya, grafik perkiraan konsumsi tampak seperti ini:


Meningkatkan


Pada titik ini, menjadi jelas bahwa tanpa penggunaan faktor eksternal tambahan, akurasi yang diperlukan tidak dapat dicapai. Pertama-tama, ini adalah faktor cuaca: suhu, tekanan, kelembaban, kekuatan dan arah angin, curah hujan. 

Upaya pertama adalah upaya untuk membeli arsip cuaca di Yandex. Orang-orang memiliki API dan dukungan teknis yang sangat baik, tetapi secara khusus di kota kami ada cukup banyak kesenjangan dalam data. Hasilnya, saya membeli arsip di layanan openweathermap.org dengan harga $ 10. Penting untuk dicatat bahwa, meskipun fakta bahwa cuaca merupakan faktor yang sangat berguna, kesalahan perkiraan jelas akan merusak model MAPE akhir. Saya bertemu informasi bahwa solusi yang tepat untuk masalah ini adalah dalam pelatihan untuk menggunakan bukan nilai historis cuaca yang sebenarnya, tetapi ramalan cuaca selama tiga hari, tetapi, sayangnya, saya tidak memiliki kesempatan seperti itu. 

Saya juga menambahkan tanda-tanda waktu, hari dalam seminggu, liburan, nomor seri hari dalam setahun, serta sejumlah besar jeda waktu dan nilai rata-rata untuk periode sebelumnya.

Semua data setelah penskalaan (MinMaxScale) dan One Hot Encoding (nilai dari setiap item kategorikal menjadi kolom terpisah dengan 0 dan 1, sehingga item dengan tiga nilai menjadi tiga kolom yang berbeda, dan unit hanya akan menjadi salah satu dari mereka) yang saya gunakan dalam kecil kompetisi tiga model pendongkrak populer XGBoost, LightGBM dan CatBoost. 

XGBoost


Banyak hal baik telah ditulis tentang XGBoost. Dipresentasikan pada konferensi SIGKDD pada tahun 2016, ia membuat percikan dan masih memegang posisi terdepan. Materi yang menarik bagi mereka yang belum pernah mendengarnya.

Mari kita siapkan data untuk model: 

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


Saya mengajar model menggunakan hyperparameter terbaik yang dipilih sebelumnya:

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 adalah kerangka kerja dari Microsoft yang keunggulan utamanya adalah kecepatan belajar pada set data yang sangat besar. Dan juga, tidak seperti XGBoost, LightGBM dapat bekerja dengan kategori, menggunakan lebih sedikit memori. Di sini kemudian ada lagi.

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 adalah pustaka gradien canggih yang meningkatkan pohon keputusan sumber terbuka "dari Yandex . Perbedaan yang menguntungkan dari model ini adalah kenyamanan bekerja dengan kumpulan data yang berisi atribut kategorikal - Anda dapat mentransfer data yang berisi kategori ke model tanpa konversi ke angka, dan ia mengungkapkan pola sendiri, dengan tangan tidak ada yang perlu diputar, dan kualitas prediksi tetap tinggi.

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


Agar tidak mengulangi banyak artikel, saya akan memberikan tabel perbandingan dari sini .


Untuk menghitung kesalahan MAPE, saya mengambil bulan terakhir yang diketahui (28 hari) dan, memindahkan jendela dengan resolusi satu jam, saya membuat perkiraan untuk 72 jam berikutnya.

Dan dalam kompetisi dadakan saya hadiahnya adalah:

Juara 3: favorit saya - XGBoost - dengan kualitas perkiraan terendah;
Juara 2: CatBoost sebagai yang paling nyaman dan "semuanya di luar kotak";
Posisi Pertama: LightGBM sebagai yang tercepat dan paling akurat.

Untuk perbandingan yang lebih akurat dari model, saya menggunakan R2 (R-squared atau koefisien determinasi, menunjukkan berapa banyak varian bersyarat dari model berbeda dari varian nilai nyata) dan RMSLE (Root Mean Squared Logarithmic Error, atau root mean square logarithmic error, yang, pada kenyataannya, jarak antara dua titik di pesawat - nilai riil dan prediksi). 

MetrikLightgbmCatboostXGBoostNabiSarimaax
r20.941370,939840,929090.814350.73778
Nyonya0,024680,024770,012190,008290,00658

Tapi semua ini, tentu saja, terkait dengan tugas saya. Pada data lain di tangan yang berbeda, semuanya bisa sangat berbeda. 

Lain kali, saya akan berbagi cerita lain tentang prediksi arus keluar karyawan untuk tugas merencanakan rekrutmen di perusahaan kami.

Source: https://habr.com/ru/post/undefined/


All Articles