ゲリラ豪雨が発生する日の太陽光発電量を予測する(2回目)

先日、ゲリラ豪雨が発生する日の太陽光発電電力量を予測するの第一回目の記事を書きました。

一回目の記事では、ゲリラ豪雨が発生した日の特定と、発生時刻に太陽光発電量が落ち込むことが特定できました。今回は、その第二回目、「直前にゲリラ豪雨を検知して、太陽光発電量予測を修正することができるのか！？」を検証してみたいと思います。

ここでは、降水の予測として、「過去気象予報データ」と「速報版降⽔短時間予報データ」の履歴を使いたいと思います。また、わかりやすいように、前回の記事でもご紹介した雨量の実績値「速報版解析雨量データ」も使いたいと思います。
速報版降水短時間予報とは、作成頻度が10分間隔に⾼頻度化された降⽔短時間予報となります。

カテゴリ	データ名	備考
予報	過去気象予報	予報①
予報	速報版降⽔短時間予報	予報②
実績	速報版解析雨量	参照データ

過去気象予報データ（予報値①）の加工

・あらかじめ、過去気象予報データを「WEATHER_FORCAST_RESULT.csv」として配置する。
・過去気象予報データの降雨量については、1時間ごとの後1時間の積算値であるため、9:00～18:00のデータのみとするためには、8:00～17:00のデータに絞る必要がある。
・予測対象日時の重複を排除するため、今回は、最新の予報のみに絞る。

#過去気象予報データの読み込み
df_forecast = pd.read_csv("WEATHER_FORCAST_RESULT.csv")
df_forecast = df_forecast.sort_values(["日時","予報演算実施日時"])
df_forecast["予報演算実施日時"] = pd.to_datetime(df_forecast["予報演算実施日時"])
df_forecast["日時"] = pd.to_datetime(df_forecast["日時"])

df_forecast['予報演算実施日'] = df_forecast['予報演算実施日時'].apply(lambda x : x.strftime('%Y/%m/%d'))
df_forecast['予報演算実施時'] = df_forecast['予報演算実施日時'].apply(lambda x : x.strftime('%H:%M:%S'))
df_forecast['日'] = df_forecast['日時'].apply(lambda x : x.strftime('%Y/%m/%d'))
df_forecast['時'] = df_forecast['日時'].apply(lambda x : x.strftime('%H:%M:%S'))

df_forecast = df_forecast.dropna().reset_index(drop=True)

#最新の予測演算実施時のデータに絞る
df_forecast_latest = df_forecast[~df_forecast.duplicated(subset='日時', keep = 'last')].reset_index(drop=True)

df_forecast_latest = df_forecast_latest[(df_forecast_latest['時']>"08:00:00") & (df_forecast_latest['時']<"18:00:00")].reset_index(drop=True)
df_forecast_latest = df_forecast_latest[["日時", "降雨量（mm/h）", "日射量（W/m^2）","中層雲量（%）"]].copy()
df_forecast_latest = df_forecast_latest.rename(columns={'日時':'timestamp', \
                                    '降雨量（mm/h）':'rain_forecast_latest','日射量（W/m^2）':'isolation_forecast_latest','中層雲量（%）':'cloudy_forecast_latest'})

速報版降水短時間予報データ（予報値②）の加工

・あらかじめ、データフレーム「df_rain_forecast_q_short」に、「速報版降⽔短時間予報データ」の履歴を格納する。
・速報版降⽔短時間予報については、10分ごとの後1時間の積算値であるため、8:00～18:00のデータのみとするためには、8:00～17:00のデータに絞る必要がある。
→ただ、8時台のデータを蓄積していなかったので、8時台のデータがありませんでした。
・10分ごとにデータが出力されるため、各時間の0分のデータのみとする。
・最新の予報のみとする。

import pandas as pd
import json
import csv
from pandas import json_normalize

df_rain_forecast_q_short["timestamp"] = pd.to_datetime(df_rain_forecast_q_short["timestamp"])

#日時が9時から18時のみのデータにしてみる
df_rain_forecast_q_short['date'] = df_rain_forecast_q_short["timestamp"].apply(lambda x : x.strftime('%Y/%m/%d'))
df_rain_forecast_q_short['time'] = df_rain_forecast_q_short["timestamp"].apply(lambda x : x.strftime('%H:%M:%S'))

df_rain_forecast_q_short = df_rain_forecast_q_short\
    [(df_rain_forecast_q_short['time']>="08:00:00") & (df_rain_forecast_q_short['time']<"18:00:00")].reset_index(drop=True)

#0分のデータのみに抽出する
df_rain_forecast_q_short['minites'] = df_rain_forecast_q_short["timestamp"].apply(lambda x : x.strftime('%M:%S'))
df_rain_forecast_q_short = df_rain_forecast_q_short\
    [df_rain_forecast_q_short['minites'] == "00:00"].reset_index(drop=True)

df_rain_forecast_q_short['timestamp'] = df_rain_forecast_q_short['timestamp'].astype(str)
df_rain_forecast_q_short['timestamp'] = df_rain_forecast_q_short['timestamp'].str[:-6]

#予測時間と予測対象時間が等しいデータのみとする
df_rain_forecast_q_short = df_rain_forecast_q_short\
    [df_rain_forecast_q_short['datetime'] == df_rain_forecast_q_short['timestamp']]

予報値①、予測値②、実績値と発電量データの加工

前回の記事でご紹介した手法で結合し、下記のようなデータ（dataframe:df_pv_）に整形します。

列名	説明
timestamp	日時
kwh	発電量
rain_value	速報版解析雨量
rain_forecast_latest	過去気象予報(雨量)
isolation_forecast_latest	過去気象予報(日射量)
cloudy_forecast_latest	過去気象予報(雲量)
rain_value_short_forecast	速報版降⽔短時間予報

発電量の予測

予報値①を使った予測

まずは、予報値①の過去気象予報(雨量)と過去気象予報(日射量)を説明変数とし、発電量の予測をするモデルを作成したいと思います。予測は重回帰分析を採用します。
重回帰分析については、こちらの記事が参考になりました。

なお、予測する日は前回の記事でゲリラ豪雨が発生したと想定した、2021/7/24とし、直前の1か月間のデータでモデルを作成します。

from sklearn.linear_model import LinearRegression

df_pv_model = df_pv_[(df_pv_["timestamp"] > "2021-06-23") & (df_pv_["timestamp"] < "2021-07-24")].reset_index(drop=True)

x = df_pv_model[['isolation_forecast_latest', 'rain_forecast_latest']]
y = df_pv_model[['kwh']]

model_lr = LinearRegression()
model_lr.fit(x, y)

print('偏回帰係数')
print(model_lr.coef_)
print('切片 w2: %.3f' %model_lr.intercept_)
print('決定係数 R^2： ', model_lr.score(x, y))

偏回帰係数
[[0.00451768 0.02813854]]
切片 w2: 0.120
決定係数 R^2： 0.8138712966779383

過去気象予報(雨量)(予報値①)の偏回帰係数がプラスの値なので、雨と発電量の関係性的に、何かおかしいなという印象です。

それでは、7/24を予測してみましょう。

import numpy as np
import matplotlib.pyplot as plt

fig, ax1 = plt.subplots()
ax2 = ax1.twinx()

df_pv_somedays = df_pv_[(df_pv_["timestamp"] >= "2021-07-24") & (df_pv_["timestamp"] <"2021-07-25")].reset_index(drop=True)

x_somedays = df_pv_somedays[['isolation_forecast_latest', 'rain_forecast_latest']]
timestamp_somedays = df_pv_somedays[['timestamp']]

#予測を実施
predict_y_somedays = model_lr.predict(x_somedays)

ax1.plot(df_pv_somedays["timestamp"], y_somedays, label="measured")
ax1.plot(df_pv_somedays["timestamp"], predict_y_somedays, label="predicted")
ax2.scatter(df_pv_somedays["timestamp"], df_pv_somedays["rain_forecast_latest"], c="y", label="rain_forecast_latest")
ax2.scatter(df_pv_somedays["timestamp"], df_pv_somedays["rain_value_short_forecast"], c="b", label="rain_value_short_forecast")
ax2.scatter(df_pv_somedays["timestamp"], df_pv_somedays["rain_value"], c="r", label="rain_value")

labels = ax1.get_xticklabels()
plt.setp(labels, rotation=45, fontsize=10)

ax1.set_xlabel("time")
ax1.set_ylabel("pv")

#凡例の表示
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1 + h2, l1 + l2, loc='upper left', bbox_to_anchor=(1, 1), facecolor='gray')

ax1.grid()

グラフを見ると、7/24 9時の発電量の落ち込み（青の折れ線）を、発電量の予測（オレンジの折れ線）が捉えられていないことがわかります。
そもそも、「過去気象予報データ（予報値①：黄色い点）」と「速報版降⽔短時間予報データ（予報値②：青い点）」の両方とも、降水量を予測できていないですね。

予報値②

次に、予報値②の速報版降⽔短時間予報データの履歴と過去気象予報(日射量)を説明変数とし、発電量の予測をするモデルを作成したいと思います。

from sklearn.linear_model import LinearRegression

df_pv_model = df_pv_[(df_pv_["timestamp"] > "2021-06-23") & (df_pv_["timestamp"] < "2021-07-24")].reset_index(drop=True)

x = df_pv_model[['isolation_forecast_latest', 'rain_value_short_forecast']]
y = df_pv_model[['kwh']]

model_lr = LinearRegression()
model_lr.fit(x, y)

print('偏回帰係数')
print(model_lr.coef_)
print('切片 w2: %.3f' %model_lr.intercept_)
# print('y= %.3fx + %.3f' % (model_lr.coef_ , model_lr.intercept_))
print('決定係数 R^2： ', model_lr.score(x, y))

偏回帰係数
[[ 0.00445267 -0.0313882 ]]
切片 w2: 0.162
決定係数 R^2： 0.8141869060509984

速報版降⽔短時間予報データ(予報値②)の偏回帰係数はマイナスであるので、予報値①より説明変数としての意味をなしているのでは、と考えられます。
こちらも、グラフを表示してみましょう。

速報版降⽔短時間予報データで降水量が捉えられていないので、発電量の落ち込みは予測できないですね。
ただ、偏回帰係数を見ると、予報値②の方が良いですが、決定係数はあまり変わらないので、雨量の予報を使って発電量を予測する、というのは今回は難しい結果となってしまいました。

まとめ

考察

・そもそも、今回の記事で扱っている説明変数（予測値）では、該当日のゲリラ豪雨を捉えられていなかった。
・雨量が多いほど、太陽光発電が下がる、というのは、感覚的に違う気もするので、重回帰分析が正しいかもわからない。

次回からのアクション

・他の日はどうなのか検証してみる。
・別データ「（速報ではない）降⽔短時間予報データ」を使ってみる。
・雨量をそのまま使った重回帰分析以外も試してみたい。