大学の電力消費実績データの解析を行うための外れ値・欠損値処理に関する研究

全文

(1)A-44 大学の電力消費実績データの解析を行うための外れ値・欠損値処理に関する研究 A study on pre-processing of outliers and missing values for energy consumption analysis at university 学生会員学生会員. ○佐橋. 寛也（中部大学）. 技術フェロー. 山羽. 基（中部大学）. Nepal Bishnu（中部大学） Hiroya SAHASHI*1. Motoi YAMAHA*1. Bishnu NEPAL*1. *1 Chubu University In the paper, we describe pre-processing methods based on the actual power consumption data of a comprehensive university. Data obtained from buildings is huge and its measurement items are diverse. Those obtained from buildings are expected to be used in the advanced data science field in the future. Therefore, the pre-processing methods described in the paper submitted earlier have been summarized, and the target data have newly coped with missing data by methods such as linear interpolation and time-series processes. はじめに建築物において、様々な事項（e.g. 消費電力, 太陽光発電量, 室内の温湿度 etc.）の計測を行う機器が取り付けられるケースが増加し、一建物から取得できるデータが多くなりつつあり、これらのデータが高度に利用される（e.g. 計測データのモデリングによる将来予測 etc.）ことが期待される。既報では、VAR モデルによる時系列分析 1)、k 平均法によるクラスタリング分析 2)を建築設備のデータに適用する場合、欠損値及び外れ値が原因で分析過程の計算が実行できないことがある。そのため、分析の対象とするデータを整える必要がある。そこで、本報での前処理は、時系列分析などの分析に繋げるためのデータの整備を目的とし、外れ値の検出精度や欠測部の予測精度については、特に議論を行わない。なお、ここで述べる前処理とは具体的には、外れ値・欠損値への対処法を指す。 1.. 外れ値・欠損値の処理方法研究活動では、主に中部大学に設置されているスマート BEMS のサーバーから取得できるデータを用いている。これは、中部大学キャンパス内に多数設置されている積算電力計（単位：kWh）からの通信によって送られるデータである。そのデータ内には少なからず異常値が存在する場合が多い。このような異常値は積算電力計の積算値が前の値より低い場合、積算電力計から積算値（外れ値）が送信される、あるいは通信エラーにより値が送信されない（欠損値）のような積算電力計の特性によるものだと考えられる。3) 2.. 空気調和・衛生工学会大会学術講演論文集｛2019.9.18 〜 20（札幌）｝第9巻. 2.1. 外れ値の検出外れ値とは、他の観測値と比べて非常に大きいまたは小さい観測値のことであり、既報 3)では標準偏差により Z スコアに標準化したデータに信頼区間を適用し、外れ値を検出する方法を報告した。しかし、Z スコアは外れ値を含む標本の平均を用いるので、外れ値が検出されなくなるまで検出処理を繰り返す必要があった。その後の報告 4)では、Tukey’s fences という手法を用いた外れ値の検出を行った。この手法は、四分位範囲の尺度に基づき、式（1）に示す下限閾と上限閾による区間を定め、区間外に該当するデータを外れ値とする方法である。 [ 𝑄1 − 𝑘(𝑄3 − 𝑄1 ) , 𝑄3 + 𝑘(𝑄3 − 𝑄1 ) ] …（1） 𝑄1 ：標本の第 1 四分位 𝑄3 ：標本の第 3 四分位 𝑘：負でない定数この報告では、下限閾の定数𝑘を 1.5 かつ上限閾の定数𝑘 を 3 として外れ値を検出することが扱うデータに則しているのではないかと考えられた。さらに通常、電力消費が 0 以下になることは考えられないため、0 以下の値は欠損値に置き換え、欠損値処理の段階でまとめて対処を図ってよいと考えられる。結果的に後者の Tukey’s fences 手法の方が、外れ値を含むおそれのある標本平均を利用しないため、一度の処理ですべての外れ値が検出されやすいと考えられ、今回扱ったデータにより則した外れ値検出法であるといえる。また、契約電力等の情報が既知である場合は、その値を超える観測値を外れ値としてしまってもよいであろう。そのような情報が手に入らない場合は、この外れ値検出法は有効になるであろう。. -181-.

(2) 2.2. 補定 2.1 で検出された外れ値を欠損値に置き換え、データに欠損値がある場合はこれと併せて適切な値を当てはめる。この操作を補定と呼ぶ。（1）多変量の扱い既報 3)ではオープンソースソフトウェア R のパッケージ mice を利用した代入法による欠損値対処について報告した。mice の代入モデルの選択肢は 24 種類あり、代入モデルの選択については分析者のデータ解析能力に依存する部分が大きい。そこで、代入モデルに選択の違いによって代入値にどの程度のばらつきがでるのかを確認した。代入モデルの訓練データとしては、2016 年 1 月 18 日（月）から対象日の 1 月 21 日（木）までの電力消費実績データ（単位： kWh，1 時間間隔）を用いた。データの変数列（5 変量）は以下に箇条書きに示すとおりである。 ⚫ 東キャンパス受電電力量（変数 1） ⚫ 西キャンパス受電電力量（変数 2） ⚫ 中地区方面電力量（変数 3） ⚫ 西地区方面電力量（変数 4） ⚫ 北地区方面電力量（変数 5）. 確認できるが、3 種類すべてで観測値から離れた代入値がみられる。したがって、この図-1 の結果のみで扱うデータに最も適した mice の代入モデルを定めてしまうことは芳しくない。また、mice は欠測を伴う変数列のモデルを構築するために他の変数列を説明変数として利用する。すなわち、単変量のデータには適用できないのである。さらに、代入モデルの選択を誤るとデータのトレンドから外れた値を欠損値に補ってしまうおそれもある。（2）単変量の扱いそこで、単変量時系列データの欠損値に対応できる欠損値の対処法についても報告 4)を行った。欠測を伴う日付の 3 日前のデータを訓練データとし、R の auto.arima 関数を用い、赤池情報量基準（AIC）を最小とする時系列過程を構築し、その過程の欠測に対応する予測値を代入するというものであった。しかし、時系列過程の選定にかなりの時間を要し、さらに日曜日のような大学が休日のデータが訓練データに含まれると観測値から外れた予測をしてしまうという問題点が見つかった。 3. 単変量データの補定方法の検討多変量データの補定は mice パッケージを用いることで適正に備える見通しができたが、多変量のエネルギー消. そして、対象日の 6, 9, 12, 15, 18, 21 時に欠測が生じていると仮定して単一代入法による検証を行った。単一代入法とは、一つの代入モデルから得られた予測値によって欠損値を置き換える手法である。その検証結果が図-1 である。図-1 の凡例において、value.norm がベイズ線形回帰モデルによる代入値、value.norm.nob が無作為な乱数による誤差項を追加した線形回帰モデルによる代入値、 value.pmm が最近隣ホットデック法（数式モデルを仮定せ. 費データが入手できる場合はあまり多くない。また、単変量に対して補定を行うために、、時系列データモデル及びパラメータを決めて行うことが適切と考えられる。2015 年度のデータにて、疑似欠損値を発生させ、前述の auto.arima 関数でモデルを特定した。 2.1（1）の 5 変量の 2015 年度の欠測パターンについて確認すると、変数 1 と 2 が同時刻に欠測が生じており、変数 3，4 ならびに 5 でも同時刻で欠測（変数 1 と 2 のパ. ず、観測値を補う手法）の特殊形態による代入値である。. ターンと異なる）がみられた。変数 1 と 2 の欠測パターンには連続するものは存在せず、もう一方の欠測パターンには連続する欠測がみられたが、続いても 2 時間分までであった。この結果から、扱うデータにおいては連続する欠測の発生は起こりにくいと考えられた。そこで、自己回帰（AR）過程、自己回帰移動平均（ARMA）過程を代入モデルとして24時間分のデータを用いてモデルを作成するという方針とした。想定する状況としては、時刻𝑡において利用可能な観測値からなる変数 1 の情報集合 Ω𝑡 = {𝑦24 , 𝑦23 , ⋯ , y1 }. 図 - 1 mice によるモデルの代入値の比較図 1 の結果をみると、代入モデルに違いによって代入値がバラついていることが分かる。観測値に近い代入値も空気調和・衛生工学会大会学術講演論文集｛2019.9.18 〜 20（札幌）｝. …（2）. を用いて、1 時間先の𝑦𝑡+1 を予測することを考える。5)この式（2）に示す情報集合から AIC を最小にする定常性のある時系列過程の検討を行った。具体的には、欠測があった 17 日分のデータを除いた 2015 年度の 349 日分（うるう年の 2 月 29 日を含む）の時系列過程の作成を行い、どのパラメータを伴う過程の頻出頻度が高いのかを調査した。2015 年度の日毎の情報集合を対象に 24 時間分、時刻. -182-.

(3) 𝑡から 1 時間ずつラグ次数を増やしていった。その時刻𝑡 ごとの最頻出過程の表を表-1 に示す。結果としては、すべての時刻で AR(2)過程の頻出頻度が最も高かった。これはラグ次数を 2 とした場合自己相関に有意性がみられるためと考えられる。全体比が最も高いのが時刻 1 の 0.691、対して最も低かったのが時刻 10 の 0.312 であった。念のため、時刻 10 の選定モデル一覧を表-2 に示す。表 - 1 時刻𝑡ごとの最頻出時系列過程時刻𝑡. 最頻出過程. 頻度. 年間出現率. 0. AR(2). 225. 0.645. 1. AR(2). 241. 0.691. 2. AR(2). 228. 0.653. 3. AR(2). 227. 0.65. 4. AR(2). 235. 0.673. 5. AR(2). 221. 0.633. 6. AR(2). 222. 0.636. 7. AR(2). 141. 0.404. 8. AR(2). 161. 0.461. 9. AR(2). 161. 0.461. 10. AR(2). 109. 0.312. 11. AR(2). 194. 0.556. 12. AR(2). 214. 0.613. 13. AR(2). 220. 0.630. 14. AR(2). 194. 0.556. 15. AR(2). 205. 0.587. 16. AR(2). 203. 0.582. 17. AR(2). 196. 0.562. 18. AR(2). 217. 0.622. 19. AR(2). 214. 0.613. 20. AR(2). 217. 0.622. 21. AR(2). 222. 0.636. 22. AR(2). 216. 0.619. 23. AR(2). 216. 0.619. 表 - 2 時刻 10 の選定モデル一覧選定過程. 頻度. 年間出現率. AR(2). 109. 0.312. ARMA(2,1). 49. 0.14. AR(1). 34. 0.097. ARMA(1,1). 32. 0.092. ARMA(2,2). 30. 0.086. ARMA(1,2). 24. 0.069. AR(0). 23. 0.066. AR(3). 16. 0.046. ARMA(0,1). 10. 0.029. AR(4). 5. 0.014. 空気調和・衛生工学会大会学術講演論文集｛2019.9.18 〜 20（札幌）｝. ARMA(2,3). 4. 0.011. ARMA(0,2). 3. 0.009. ARMA(3,2). 3. 0.009. ARMA(1,3). 2. 0.006. ARMA(4,1). 2. 0.006. ARMA(4,2). 2. 0.006. ARMA(0,3). 1. 0.003. 表-2 の結果によると、時刻 10 において選定された時系列過程は 17 種類に及んでいるが、今回の点予測による補定は定常 AR(2)過程をすべての時刻で用いることとする。 AR（2）過程を式（3）に示す。 𝑦𝑡 = 𝑐 + 𝜙1 𝑦𝑡−1 + 𝜙2 𝑦𝑡−2 + 𝜖𝑡 , 𝜖𝑡 ∼ W. N. (𝜎 2 ) …（3） 𝑦𝑡 ：時刻𝑡の変数𝑦の値 𝑐：定数 𝜙1 , 𝜙2 ：係数 𝜖𝑡 ：確率変動を表現する擾乱項 W. N. (𝜎 2 )：分散𝜎 2 のホワイトノイズまた、定常とは時系列過程の弱定常性を満たすことを意味し、任意の𝑡と𝑘に対して、式（4）と（5）が成立することを指す。 𝐸 (𝑦𝑡 ) = 𝜇 …（4） Cov(𝑦𝑡 , 𝑦𝑡−𝑘 ) = 𝐸 [(𝑦𝑡 − 𝜇)(𝑦𝑡−𝑘 − 𝜇)] = 𝛾𝑘 …（5） 𝐸 (𝑦𝑡 )：時系列過程𝑦𝑡 の期待値 𝜇：時刻𝑡によらず一定な時系列過程𝑦𝑡 の平均 Cov(𝑦𝑡 , 𝑦𝑡−𝑘 )：𝑘次の自己共分散 𝛾𝑘 ：時刻𝑡によらず一定な𝑘次の自己共分散そして、欠損値の対処法として頻繁に用いられる線形補間（1 次補間）と点予測の比較を行うこととする。補間係数を𝛼（式（6））とすると、その線形多項式は式（7）として表される。 𝛼=. 𝑥−𝑥0 𝑥1 −𝑥0. …（6）. 𝑦 = (1 − 𝛼 )𝑦0 + 𝛼𝑦1 …（7） [ ] 𝑥：区間 𝑥0 , 𝑥1 にある変数𝑥 𝑦：𝑥によって定まる代入値ベクトル 𝑦0 , 𝑦1 ：欠測の前後の観測値これらの式（6）と（7）により、点予測の場合は観測値 𝑦0 , 𝑦1 の算術平均を代入することになることが分かる。図-2 と図-3 に 2016 年 1 月 3 日（日）と 1 月 21 日（木）の AR(2)過程の予測値及び線形補間の代入値の比較結果を示す。どちらの結果においても、線形補間の補定結果の方が観測値との精度面で優っているようにみえる。しかしながら、21 日の時刻 8 から 9 に 910kWh 増加に対して時刻 9 の観測値は 2422kWh であったのだが、線形補間による代入値が 2035kWh、AR（2）過程による予測値が 1900kWh であった。さらに、時刻 10 の AR(2)モデルの予測値が3269kWh とデータのトレンドから大きく外れた補. -183-.

(4) 定となってしまっている。AR モデルでは過去のデータの影響を受けるので、増加する傾向が継続するような結果になっている。. 連続した 2 時間分の欠測の場合については、そのような欠測がみられた変数 3 における 2015 年 8 月 24 日（月）の補定結果を図-4 に示す。電力消費トレンドが AR(2)過程のような推移となることは考えにくく、それまでの増加の傾向の影響を受けている。総括本報では、大学の電力消費実績データに則した前処理について述べた。はじめに、建築物から得られるデータは膨大でかつ、その計測項目が多岐にわたることに触れ、建築設備に関するデータの高度利用が期待されると論じた。そして、データの分析・活用には扱うデータの整備（前処理）が重要であることを示唆した。そこで、大学の電力消費実績データの傾向に悪影響を与えることなく、そのデータ中の異常値を対処する手法について取り纏め、新たに AR(2)モデルと線形補間による欠損値対処の比較を試みた。本報では、AIC による判定により AR(2)モデルが選択される割合が高く、これを用いた補定を行ったが、1 点の欠損値に対しては、線形補間よりも精度が適切であるとは判断できなかった。AIC 4.. 図 - 2 AR(2)過程と線形補間による欠測対処（1）. ではモデルの単純さが判定の基準になるので、次数の大きいモデルを適用するほうが良い場合も考えられる。本報の検討の範囲における中部大学の電力消費実績データに則した前処理として、外れ値に関しては Tukey’s fences 手法、欠損値に対しては線形補間が最適な方法であると現時点では結論付けられる。. 図 -3. 参考文献 1）佐橋寛也, 山羽基, 横江彩, Nepal Bishnu : 多棟建物のピークカットに寄与する節電目標算出方法の検討, 平成 30 年度空気調和・衛生工学会大会（愛知）学術講演論文 2） Nepal Bishnu, 山羽基, 河村貢, 横江彩, 佐橋寛也 : 多棟の既存建物の低炭素化に向けたエネルギーマネジメントに関する研究（第 3 報）クラスタリング手法を用いた要因別エネルギー消費の分析とエネルギーマネジメント, 2018.9, 空気調和・衛生工学会論文集, 258 号, 21-28 3）佐橋寛也, 山羽基, 横江彩, Nepal Bishnu : BEMS データの欠損値・外れ値の扱いに関する研究, 平成 30 年度日本建築学会大会（東北） 4）佐橋寛也, 山羽基, Nepal Bishnu : BEMS データの外れ値・欠損値の前処理方法について, 空気調和・衛生工学会中部支部平成 30 年度学術研究発表会 5）沖本竜義, 『経済・ファイナンスデータの軽量時系列分析』第 13 刷, 朝倉書店, p.56-72, 2018.6. AR(2)過程と線形補間による欠測対処（2）. 図 - 4 2 時間分の連続欠測に対する補定の比較空気調和・衛生工学会大会学術講演論文集｛2019.9.18 〜 20（札幌）｝. -184-.

(5)