はじめに

これまで電力需要予測という言葉を使って何度も扱ってきたテーマですが、今年の５月のデータがとても特徴的だったので再度記事を書いてみます。

［過去の記事］

動作環境

Colaboratory
Python 3.6.9
Pandas 1.0.4
Numpy 1.18.5
scikit-learn 0.22.2
TensorFlow 2.2.0
Keras 2.3.1

作成したファイルは以下のURLにて参照可能です。
https://colab.research.google.com/drive/1j_XFPLm7ucsRi63AzPbM0UaZ5Jjk4625?usp=sharing

使用データ

東京電力電気使用量実績
気象庁気温

ライブラリ読込

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from sklearn import model_selection
from sklearn.preprocessing import StandardScaler

import keras
from keras.models import Sequential
from keras.layers import Dense
from keras import optimizers

学習データ準備

電力使用量データの取得と加工

2019年のデータを取得し、学習用に使用します。

# URLを指定して電力データをダウンロード
url = "https://www.tepco.co.jp/forecast/html/images/juyo-2019.csv"
df_MW = pd.read_csv(url, encoding="shift_jis", skiprows=2)
df_MW.shape

# 日付と時刻のデータを合わせて日時のデータとする
df_MW["timestamp"] = df_MW.DATE + " " + df_MW.TIME

# 日時のデータをシリアル値に変換する
df_MW.timestamp = df_MW.timestamp.map(lambda _: pd.to_datetime(_))

# 使用する列を取得
cols = df_MW.columns[2:]
df_MW = df_MW[cols]

# 列名を使用しやすくなるよう変更
df_MW.columns = ["MW", "timestamp"]

# 1万kW = 10,000kW = 10MW なので、実績値を10倍して列名をMWとする
df_MW["MW"] = df_MW["MW"] * 10

# 日時データをインデックスに指定
df_MW.index = df_MW.pop("timestamp")

気温データの取得と加工

気象庁 | 過去の気象データ・ダウンロード

上記URLより、東京の2019年1月1日から12月31日の時別の気温データをダウンロードします。
ダウンロード後、data.csvをColaboratoryにアップロードします。

df_temp = pd.read_csv("data.csv", encoding="shift_jis", skiprows=4)

# 日時と気温のデータのみに変換
cols = df_temp.columns[:2]
df_temp = df_temp[cols]

# 列名を変更
df_temp.columns = ["timestamp", "temp"]

# 日時のデータをシリアル値に変換する
df_temp.timestamp = df_temp.timestamp.map(lambda _: pd.to_datetime(_))

# 日時データをインデックスに指定
df_temp.index = df_temp.pop("timestamp")

データの結合

電力使用量と気温のデータを日時をそろえて結合する。

# インデックスを基準にデータフレームを結合
df = pd.concat([df_MW, df_temp], axis=1)

# 無効データ(NaN)の削除
df = df.dropna()

学習データの作成

学習データとするために以下の処理を行う

月、時、曜日のデータを追加し、One-Hotエンコーディングを施す
説明変数と目的変数に分割
学習データと検証データに分割
正規化

# 月、時、曜日のデータを追加
df["month"] = df.index.month
df["hour"] = df.index.hour
df["weekday"] = df.index.weekday

# One-Hotエンコーディング
cols = ["month","hour","weekday"]
for col in cols:
    df = df.join(pd.get_dummies(df[col], prefix=col))
    df.pop(col)

# 説明変数と目的変数に分割
X_cols = df.columns[1:]
y_cols = ['MW']

X = df[X_cols]
y = df[y_cols]

# 学習データと検証データに分割
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size=.1, random_state=42)

# 正規化
scaler = StandardScaler()
scaler.fit(X_train)

X_train = scaler.transform(X_train)
X_test = scaler.transform(X_test)

学習

モデルの作成

# モデルの作成
model = Sequential()

model.add(Dense(64, input_dim=len(X_cols), activation='relu'))
model.add(Dense(64, activation='relu'))
model.add(Dense(1))

# 損失関数
# 参考URL https://keras.io/ja/losses/
loss = 'mean_squared_error'

# 最適化
# 参考URL https://keras.io/ja/optimizers/
optimizer = optimizers.RMSprop()

# 評価関数
# 参考URL https://keras.io/ja/metrics/
metrics = ['mean_absolute_error']

model.compile(loss=loss, optimizer=optimizer, metrics=metrics)

model.summary()

学習の実行

class PrintDot(keras.callbacks.Callback):
  def on_epoch_end(self, epoch, logs):
    if epoch % 10 == 0:
      #print(logs)
      print("\r{:4d} {:5.2f}, {:5.2f}".format(epoch, logs['mean_absolute_error'], logs['val_mean_absolute_error']))
    print('.', end='')
 
EPOCHS = 200
 
history = model.fit(X_train, y_train, 
                    epochs=EPOCHS, validation_split=0.2, 
                    batch_size=16, verbose=0,
                    callbacks=[PrintDot()])

学習の経過を可視化

plt.figure()
plt.xlabel('Epoch')
plt.ylabel('Mean Abs Error [MW]')
plt.plot(history.epoch, np.array(history.history['mean_absolute_error']),
          label='Train')
plt.plot(history.epoch, np.array(history.history['val_mean_absolute_error']),
          label = 'Val')
plt.legend()
plt.ylim([1000,1500])

plt.figure()
plt.xlabel('Epoch')
plt.ylabel('Loss [MW]')
plt.plot(history.epoch, np.array(history.history['loss']),
          label='Train')
plt.plot(history.epoch, np.array(history.history['val_loss']),
          label = 'Val')
plt.legend()
plt.ylim([2500000,4000000])

学習結果の確認

決定係数により学習結果を確認します。

from sklearn.metrics import r2_score

# 予測値の取得
result = model.predict(X_test)
print(r2_score(y_test, result))

実行すると以下の値が表示される。

0.9004605188281153

なかなかいい感じ。

検証

新規データによる予測

2020年のデータを使って学習結果の検証を行います。
2020年の電力使用量実績と気温のデータを使用するため、気温のデータをあらかじめダウンロードしておきます。
なお、ここではダウンロードした気温のデータはdata-2.csvとなっている。

# 電力使用量実績

url = "https://www.tepco.co.jp/forecast/html/images/juyo-2020.csv"
df_MW = pd.read_csv(url, encoding="shift_jis", skiprows=2)

df_MW["timestamp"] = df_MW.DATE + " " + df_MW.TIME
df_MW.timestamp = df_MW.timestamp.map(lambda _: pd.to_datetime(_))

cols = df_MW.columns[2:]
df_MW = df_MW[cols]
df_MW.columns = ["MW", "timestamp"]

df_MW.MW = df_MW.MW * 10

df_MW.index = df_MW.pop("timestamp")

# 気温

df_temp = pd.read_csv("data-2.csv", encoding="shift_jis", skiprows=4)

cols = df_temp.columns[:2]
df_temp = df_temp[cols]

df_temp.columns = ["timestamp", "temp"]

df_temp.timestamp = df_temp.timestamp.map(lambda _: pd.to_datetime(_))
df_temp.index = df_temp.pop("timestamp")

# データの結合
df = pd.concat([df_MW, df_temp], axis=1)
df = df.dropna()

df["month"] = df.index.month
df["hour"] = df.index.hour
df["weekday"] = df.index.weekday

cols = ["month","hour","weekday"]
for col in cols:
    df = df.join(pd.get_dummies(df[col], prefix=col))
    df.pop(col)

# 不足する列を追加
for col in X_cols:
  if col not in df.columns:
    df[col] = 0

X_valid = df[X_cols]
y_valid = df[y_cols]

# 正規化
X_valid = scaler.transform(X_valid)

# 予測値の取得
result = model.predict(X_valid)

# 決定係数の確認
print(r2_score(y_valid, result))

実行すると以下の数値が表示され、学習時より若干悪い値となっている。

0.8207052892915271

可視化と考察

# 可視化

df_result = y_valid.copy()
df_result["predicted"] = result

df_result[-60*24:].plot(figsize=(15,5))

この結果を見て、5月の新型コロナウイルス感染症予防やオンライン帰省の影響からか前年より電力使用量が少なくなっているのかなと考えられますね。
全体的に夜間の電気使用量が少ないのも飲食店などが早めに閉店するなどの影響なのかもしれません。

考察すると楽しいですね。

解説動画

本記事の内容をYoutubeで解説しています。あわせてご覧下さい。

20200616 電力使用量予測 with Keras