2023年の夏は一番暑い夏といわれ、電気料金の節約が叫ばれていました。
そこで電気料金を分析したら面白いかもと思って予想してました。
ただ、やってみてわかったことは電気の構造って難しい！！

１．自己紹介

IT業界でプロジェクトのPMO/PLを担当してる、プログラムのかけない人間です。ふと、データ分析勉強しようと思って、Aidemy受講しました。

２．実施環境

・Windows11
・Jupyter notebook
・Python 3

３．データ準備

①データ調査
・経済産業省資源エネルルギーサイトの電力調査統計の統計表一覧から、
　火力が現時点で約80%を占めることが判明しました。
    統計表一覧｜資源エネルギー庁 (meti.go.jp)
　火力エネルギーの中でも、ほぼ石炭とLNGが占めていることが判明した
　ので、石炭とLNGの価格を取得することに決めました。
　（下記、参考）
　　⽯炭２５４．０億 kWh（３５．３％）、
　　ＬＮＧ ２６７．１億 kWh（３７．２％）
　・天気によっても電力の利用料は変化しており、実際の単価にも影響を
　　与えていることが、自然エネルギー財団からの指標で分かったので
　　気温も取得することに決めました。
　　（下記、参考URL）
　　https://www.renewable-ei.org/statistics/electricitymarket/?cat=retail
　②では、データ準備です。
　　東京の2019年４月～2022年3月で、下記のサイトからデータを
　　集めました。
　　・東京の電気データ　：スポット価格・需給価格
　　　https://www.renewable-ei.org/statistics/electricity/#demand
　   ・石炭：https://pps-net.org/statistics/coal3
　　・LNG：https://pps-net.org/statistics/gas
　　・天気：https://www.data.jma.go.jp/gmd/risk/obsdl/index.php　

４．実装

　(1)データ準備
　　　①データ集め
　　　　各サイトのデータを金額フォルダに格納して、データ読み込み
　　　・電気データの読み込み

import pandas as pd
import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt
from tensorflow.keras.datasets import mnist
from tensorflow.keras.layers import Activation, Dense, Dropout
from tensorflow.keras.models import Sequential, load_model
from tensorflow.keras import optimizers
import matplotlib.pyplot as plt
import seaborn as sns

# 電気データの読み込み
kw_df1 = pd.read_csv("./金額/download_201904.csv",header=1)
kw_df2 = pd.read_csv("./金額/download_202004.csv",header=1)
kw_df3 = pd.read_csv("./金額/download_202104.csv",header=1)

kw_df=pd.concat([kw_df1,kw_df2,kw_df3])
kw_df["日付"] = pd.to_datetime(kw_df["日付"])
kw_df.head()

　・時刻のデータも読込みましたが、石炭が1か月ごと・LNGは１日の
      ごとの単価だったこともあり、日ベースで平均を取得することに。

kw_df=kw_df.groupby("日付").mean()
kw_df.head()

　・天気の読み込み

# 気温データの読み込み
temp_df1 = pd.read_csv("./金額/data_201904.csv",header=3,encoding="shift-jis")
temp_df2 = pd.read_csv("./金額/data_202004.csv",header=3,encoding="shift-jis")
temp_df3 = pd.read_csv("./金額/data_202104.csv",header=3,encoding="shift-jis")
temp_df=pd.concat([temp_df1,temp_df2,temp_df3])
temp_df["年月日"] = pd.to_datetime(temp_df["年月日"])
temp_df=temp_df.rename(columns={"年月日":"日付"})
temp_df=temp_df.groupby(by=["日付"]).mean()
temp_df.describe()

・LNG価格の読み込み

# LNG価格データの読み込み
gas_df = pd.read_csv("./金額/gas.csv",encoding="shift-jis")
gas_df["日付"] = pd.to_datetime(gas_df["日付"])
gas_df.head()

・石炭価格の読み込み

# 石炭価格データの読み込み
coal_df = pd.read_csv("./金額/coal.csv",header=4,encoding="shift-jis")
coal_df.head()

      ②データ結合
　　・石炭価格が月単位なので、石炭以外で結合します。

df=pd.merge(temp_df,gas_df,on="日付",how='outer')
df=pd.merge(df,kw_df,on="日付",how='outer')

　　・日付のフォーマットを”年”,“月”,”日“,”曜”に分けます。

hour = pd.to_datetime(df['日付'], format='%Y/%m/%d')

df["日付"]=hour
year_int=[]
month_int=[]
day_int=[]
hour_int=[]
Weekday_int=[]
for i in hour:
    year_int.append(int(i.year))
    month_int.append(int(i.month))
    day_int.append(int(i.day))
    Weekday_int.append(i.weekday())
df["年"],df["月"],df["日"],df["曜日"]=year_int,month_int,day_int,Weekday_int
df.head()

・年・月・日をもとに石炭を結合します。

df_columns=["年","月","日"]
df=pd.merge(df,coal_df,on=df_columns,how="outer")
df.head()

　　　③データ下準備
　　　　a. 石炭価格のデータを前データで穴埋め作業

 #穴埋め 
df=df.fillna(method='ffill')
df.describe()

　　　　b.重複の列・列の順番を整理。併せて、Nanのデータも削除

datecolumns=["年","月","日","曜日","気温(℃)","終値","政策コスト除く","火力","水力","バイオマス","太陽光","風力","需要","東京"]
df=df[datecolumns]
df.dropna()

　　　　c.データを眺めます。

　　　d.外れ値の除去
データを眺めた結果、下記のことがわかりました。
・東京（電気単価）は、15以上はほぼなし
・政策コスト（石炭単価）は20以上ほぼなし
・終値（LNG単価）は６以上ほぼなし
・日に大きな差はなかったので、曜日だけ残して日は削除してよさそう。

 #外れ値除去 
df = df.loc[df["東京"] <= 15]
df = df.loc[df["政策コスト除く"] <= 20]
df = df.loc[df["終値"] <= 6]
cols = ["年","月","曜日","気温(℃)","終値","政策コスト除く","火力","水力","太陽光","需要","東京"] 
df=df[cols]
df.describe()

　　(2)各モデルで分析
　　①trainデータとtestデータに分けます。

X=df.iloc[:,:-1]
y=df.iloc[:,-1]
X.shape

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test =train_test_split(X, y, test_size=0.3, random_stat＃e=0) #正規 
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
scaler.fit(X_train)

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

　②モデル分析・スコア確認
　　# RandomForestRegressor

from sklearn.ensemble import RandomForestRegressor

model = RandomForestRegressor()
model.fit(X_train, y_train)

from sklearn.metrics import r2_score
# 予測結果
y_pred = model.predict(X_test)

r2_score(y_test, y_pred)

　　　　　出力結果：0.7507445789670499
　 #リッジ回帰

from sklearn.linear_model import Ridge
# リッジ回帰
model = Ridge()                         # リッジ回帰モデル
model.fit(X_train, y_train)
result = model.predict(X_test) #mse = mean_squared_error(y_test, result)# 評価

print(r2_score(y_test, result))

　　　　　　出力結果：0.5625766847064508
　　# ディープラーニング（keras）

# ニューラルネットワークのモデルを学習させる
# エポック数は100、バッチサイズは32
model.fit (X_train,y_train , batch_size=32, epochs=100, verbose=0)

# ニューラルネットワークのモデルの性能を評価する
# テストデータを使って予測値と実測値の誤差を計算する
test_loss, test_mae = model.evaluate(X_test, y_test)
print('Test loss:', test_loss)
print('Test MAE:', test_mae)

# 予測値の取得
y_pred = model.predict(X_test)

from sklearn.metrics import mean_squared_error

 #二乗平方根で誤差を算出 
mse = mean_squared_error(y_test, y_pred)
print("KERAS REG RMSE : %.2f" % (mse ** 0.5))
print("KERAS REG SCORE: %3.2f" % (1.0 - mse ** 0.5 / y_test.mean()))

出力結果

# LightGBM

import lightgbm as lgb


# early_stopping用の評価データをさらに分割
X_train_v, X_valid, y_train_v, y_valid = train_test_split(X_train, y_train, test_size=0.25, random_state=42)

# データをDatasetクラスに格納
dtrain = lgb.Dataset(X_train_v, label=y_train_v)  # 学習用
dvalid = lgb.Dataset(X_valid, label=y_valid)  # early_stopping用
# 使用するパラメータ
params = {'objective': 'regression',  # 最小化させるべき損失関数
         'metric': 'rmse',  # 学習時に使用する評価指標(early_stoppingの評価指標にも同じ値が使用される)
         'random_state': 42,  # 乱数シード
         'boosting_type': 'gbdt',  # boosting_type
         'verbose': -1  # これを指定しないと`No further splits with positive gain, best gain: -inf`というWarningが表示される
         }
verbose_eval = 0  # この数字を1にすると学習時のスコア推移がコマンドライン表示される
# early_stoppingを指定してLightGBM学習
gbm = lgb.train(params, dtrain,
                valid_sets=[dvalid],  # early_stoppingの評価用データ
                num_boost_round=10000,  # 最大学習サイクル数。
                callbacks=[lgb.early_stopping(stopping_rounds=10, 
                                verbose=True), # early_stopping用コールバック関数
                           lgb.log_evaluation(verbose_eval)] # コマンドライン出力用コールバック関数
                )

# スコア(RMSE)算出
y_pred = gbm.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
print("LightGBM REG RMSE : %.2f" % (mse ** 0.5))
print("LightGBM SCORE: %3.2f" % (1.0 - mse ** 0.5 / y_test.mean()))

出力結果

５．結論

・分析の結果、ディープラーニング（keras）の分析で一番いい指標が
　出ました！　パラメータをもうちょっと調整すればもっといい指標が
　出ると思いましたが、体力が及ばず。。
・2019年～2022年度指標をもとにしましたが、自然エネルギー（風力等）の　
　調整指標も今後見ていく必要ありますね。