計 測 自 動 制 御 学 会 論 文 集 Vol.xx, No.xx, 1/9（2019）
予測残差補正による太陽光発電量予測及び階層型 BEMS
山 内
賢
∗・滑
川
徹
∗,∗∗PV Power Prediction via Residual Correction and
Hierarchical Multiple Buildings Energy Management
Ken Yamauchi
∗and Toru Namerikawa
∗,∗∗This paper deals with photovoltaic (PV) power prediction and building energy management system (BEMS) in multiple buildings. First, we consider PV power prediction. We predict PV power by using k-means method and SVR. However, there is an error between prediction and actual generation. In order to solve this issue, we estimate prediction error in advance by markov process and revise the prediction. And, by applying copula to calculated point prediction, calculate interval prediction. Second, we consider 2 level BEMS. At first level, in each building, power consumption is determined considering comfort. At second level, in order to reduce total cost in all buildings, battery scheduling is calculated with robustness by applying scenario robust theory. And, at the end of each section, simulation results show the effectiveness of proposed prediction method and energy management method.
Key Words: PV power prediction, markov process, energy management, model predictive control
1.
は じ め に 近年,地球温暖化や化石燃料の枯渇化など環境問題が深刻 化しており,再生可能エネルギーの都市への大量導入が注目 されている. しかし,再生可能エネルギーは自然エネルギーで あるため,発電量や周波数が気候条件に大きく左右されてし まうという問題点がある. 電力ネットワークなどの大規模シ ステムに導入する場合,周波数や電圧が大幅に乱れてしまい 大規模停電につながる恐れもある1) . そのため,再生可能エネ ルギーの発電量を事前に予測することで発電計画を行い,シ ステム全体を通して電力供給量を安定させる必要がある. ま たその一方で,予測誤差が大きいと需給バランスの許容値を 上回ってしまい問題となること, 規模の小さな発電機では自 然環境の影響を受けやすいこと, 予測誤差を低減できれば不 安定性に対処するために設置する蓄電池の容量が抑えられ経 済的であることなどから,予測精度の向上が重要である. このように,再生可能エネルギーや他のエネルギーを組み合 わせ,これらを制御するようなシステムがエネルギーマネジメントシステム(EMS:Energy Management System)である. EMSの対象は様々であり,ビルを対象としたBEMS(Building
EMS),家庭を対象としたHEMS(Home EMS),地域社会を
∗ 慶應義塾大学大学院理工学研究科 横浜市港北区日吉 3-14-1 ∗∗ 独立行政法人科学技術振興機構 CREST 川口市本町 4-1-8
∗ Graduate School of Science and Engineering, Keio
Uni-versity, 3-14-1 Hiyoshi, Kohoku-ku, Yokohama
∗∗ JST, CREST, 4-1-8, Honcho, Kawaguchi
対象としたCEMS(Community EMS)などがある. 具体的 には,節電のための空調機器などの制御や,蓄電池の導入,電気 自動車との連携によりエネルギーを有効利用するための制御 が行われている. また,将来的には,より規模の大きなCEMS の運用が注目されている. 本論文では, CEMSを見据え,コ ミュニティに含まれるビル群を対象としたエネルギーマネジ メントについて考える.これらの研究背景から,本研究では発 電量予測と複数のビルを対象としたエネルギーマネジメント システムについて扱う. 再生可能エネルギーの発電予測の研究は盛んに行われてお り,発電予測を行う対象としては,主に風力発電と太陽光発電 の二つが存在する. 風力発電予測手法としては, JIT(Just-in-time)モデリングを用いた手法2) ,機械学習を用いた手法3) , 太陽光発電に関してはSVMを用いた手法4) ,ニューラルネッ トワークを用いた手法5)などがある. また,本研究でも扱う ように,発電量を確率的に予測し,幅を持たせた区間によって 評価を行う予測手法も考えられている6), 7) . そして,これら 発電予測を用いた研究としては,発電予測値に基づいたマイク ログリッドの制御8)などが行われている. エネルギーマネジ メントシステムの研究としては,再生可能エネルギーの発電量 の不確かさを考慮し,ロバスト性を考慮したEMSの研究9) , 熱資源に基づいた空調管理システムを対象とした研究10)が 挙げられる. また,近隣の建物との電力融通を考えることで電 力コストの低減を目指す研究11)などがある. 本論文では,より精度の高い太陽光発電予測手法を提案す ること,複数のビルを対象としたBEMSにおいて総コストを TR 00xx/19/xxxx–0001 c⃝ 2019 SICE
削減する運用計画を提案することを主な目的とする. 太陽光 発電予測については,まず従来手法を用いて予測値を算出し予 測値と実測値との予測残差データを収集することで,予測残差 の推定を行うことを考える. そして,推定された予測残差によ り予測値の修正を行う. BEMSにおいては,複数のビル間に おいて電力の融通を図ることによって全てのビルにおける総 コストを削減することを考える. その際,発電予測値の予測誤 差により運用計画と実際の運用に差が生まれてバッテリー制 約を犯す可能性があるため,シナリオロバスト理論に適用す ることでロバスト性を持ったBEMSの運用を行う. 文献12) では,ビルが1台の場合の問題を扱っているが,本論文では複 数台に拡張した問題について扱う. 本論文の構成としては以下の通りである. まず, 2章で太陽 光発電予測手法について説明し, 3章でエネルギーマネージメ ントについて説明する. 各章末において,それぞれシミュレー ション結果を示す. 最後に, 4章で結論を述べる.
2.
予測残差補正による太陽光発電予測 文献13)では ,予測残差データを集めた後,ある時刻での予 測残差とその次の時刻での予測残差にマルコフ性があると考 えている. つまり,ある時刻での予測残差が得られたとき,次 の時刻の予測残差をそのマルコフ性から推定し,発電量予測値 を修正することができる. 文献13)では,過去の予測誤差デー タに対してクラスタリングを行い,そのクラスの遷移確率から 予測残差を推定する. それに対して本論文では,過去の予測誤 差データに対してコピュラ14)を適用し ,誤差データの推移を 同時分布関数として表すことで予測残差を推定する. コピュ ラとは, 1次元周辺分布関数を接合して多次元分布関数を生成 する関数のことである. 離散的なクラスタリングを用いずに 連続的な確率分布を用いることで,予測残差の推定をより正確 に行うことができ,高精度な予測が可能となると考えられる. 2. 1 予測アルゴリズム 本研究の予測アルゴリズムの概要をFig. 1に示す. Step Ⅱ PV prediction by k-means method, SVRPV prediction after revision
+ + Learning stage Step Ⅴ Calculation copula parameter Step Ⅲ
Calculation prediction residual Step Ⅳ
Normalization Data
Step Ⅵ
Calculation joint distribution function of Step Ⅶ Creation probability density function of PV prediction algorithm Step Ⅰ Collection data Past weather Past PV generation Weather forecast
Fig. 1 Flow chart of forecast algorithm
次に,各stepの内容について説明する. step I. データの収集 過去の発電データ,天気データ,天気予報データを収集する. step II. k-平均法, SVRによる予測値の算出 文献12)を参考に , k-平均法, SVR(サポートベクター回帰) を用いて予測値を算出する. step III. Rb, Raの算出 StepIIにおける予測の結果から予測発電量Pを,実データ から過去の発電量T Gを得る. 予測残差Rは,相対残差 (relative residual)を用いて以下の式から算出される. R(t) =P (t)− T G(t) PST C （1） PST C は予測対象の定格発電量を表す.時刻tについては, 本研究では6時∼18時の12時間を予測時間とするため, t = 1, 2,· · · , 12とする. そして,予測残差の推移を表すた めに, Rb, Raについては以下のように表す. Rb(k) = R(t| t = 1, 2, · · · , 11) （2） Ra(k) = R(t| t = 2, 3, · · · , 12) （3） (k = 1, 2,· · · , 11) step IV. データの正規化 データを扱いやすくするために正規化を行う. Rb, Raの最 大値,最小値を各々Rb
max, Rbmin, Ramax, Rmina として正規 化されたデータRbr, Rarを以下のように計算する. Rbr = Rb− Rb min Rb max− Rbmin （4） Rar = Ra− Ramin Ra max− Ramin （5） step V. コピュラパラメータの算出 正規化された予測残差データRbr, Rar の同時分布を算出す るために,本研究ではtコピュラを用いる. そのため, tコ ピュラパラメータである自由度νと相関行列Σを以下のよ うに求める. ν = 4.0275, Σ = [ 1.0000 0.9188 0.9188 1.0000 ] （6） ここで,相関行列Σはデータの相関を表しており,値が1 に近いことから, tコピュラが適していることが分かる. ま た,本論文ではMATLABのcopulafitを用いて算出する. step VI. Rb, Raの同時分布hの算出 step Vで求めたコピュラパラメータν, Σを用いて,Rbr, Rar の同時分布hを生成する. このstepまでが過去のデータを 用いた学習期間となっており,次のstepが実際の時刻にお ける予測段階であることに注意していただきたい. step VII. 確率密度関数fRˆr(xRˆr)の算出と予測値の修正 同時分布hから,現時刻の予測残差Rr(t)に対する次の時刻 の予測残差Rrˆ (t)の確率密度関数fRˆr(xRˆr)を式(5)より求 める. xRb r, xR a rは各々R b r, Rarに対する確率変数, f (xRa r) はxRa r を確率変数とする確率密度関数である. fRˆr(t+1)(xRˆr(t+1)) = fRa r(xRar|xRbr= Rr(t)) = h(xRa r, xRb r = Rr(t)) （7） 得られた確率密度関数を用いて,時刻t + 1での推定予測残 差Rrˆ (t + 1)を以下のように求める.
計測自動制御学会論文集 第 xx 巻 第 xx 号 2019 年 xx 月 3 ˆ Rr(t + 1) = ∫ 1 0 fRˆr(xRˆr)xRˆrdxRˆr （8） ˆ Rr(t + 1)は正規化された値であるため,式(9)より実際の スケールに戻す. ˆ
R(t + 1) = Rmin+ ˆRr(t + 1)(Rmax− Rmin)（9） そして最後に,推定予測残差R(t + 1)ˆ を用いて点予測によ る予測値P (t + 1)を式(10)のように修正し,新しい予測 値Pp(t + 1)を算出する. Pp(t + 1) = P (t + 1) 1 + ˆR(t + 1) （10） 2. 2 シミュレーション 2. 2. 1 シミュレーション条件 太陽光発電データに関しては，本研究では慶應義塾大学矢 上キャンパス25棟屋上に設置されている発電システムから得 られる過去のデータを用いる. 予測を行う時刻は6時から18 時の12時間で, 1時間毎の予測を行う. また,提案手法による 残差の修正を行うためには基となる予測結果が必要であるた め,過去の発電データ取得期間,従来手法による予測及び予測 残差データ取得期間,提案手法による予測修正期間の3つに 分けて予測を行った. 予測期間についてFig. 2にまとめた. Data1~Data70
70 days Data71~Data9525 days Data96~Data13439 days Past data collection
period by conventional method PV prediction period PV prediction revision periodby proposed method
PV generation data
Weather data PV prediction residual dataPV prediction data PV prediction after revision
Fig. 2 prediction schedule
予測手法についてはk-平均法, SV Rのみを用いた従来手法
(Without revision),マルコフ連鎖に基づいた予測残差補正手
法(Method1)13),提案手法(Proposed)で比較を行った. ま
た,評価方法については以下のMRE(Mean Relative Error)
を用いる. T Gは発電実測値, P は発電予測値である. M RE = 1 PST C 1 h h ∑ t=1 ∥T G(t) − P (t)∥ × 100[%]（11） 2. 2. 2 シミュレーション結果
Fig. 3に39日分のMRE値を, Table 1にその平均値を
示す. 5 10 15 20 25 30 35 Day 0 5 10 15 MRE[%]
Mean Relative Error Without revision (Average: 7.057[%]) Method 1 (Average: 6.615[%]) proposed (Average: 6.383[%])
Fig. 3 Point Forecast result
Table 1 Point Forecast result
Without revision Method1 Proposed MRE[%] 7.057 6.615 6.383 Fig. 3, Table 1から,残差補正を行わない場合や手法1と 比較して,提案手法では予測誤差が小さくなっていることが 確認できた. 離散的なクラスタリングではなく連続的な分布 関数を用いたことで,予測誤差が小さくなったと考えられる.
3.
階層型BEMS
3. 1 問題設定 この章では,複数ビルにおける階層型エネルギーマネージ メントについて扱う. まず,本論文で想定する問題設定につい て説明する. 各ビルは太陽光発電機とバッテリーを備えてお り,太陽光発電によって得られた電力をバッテリーに蓄電し, 電力を使用する際はバッテリーから電力を放電する. また,各 ビルはバッテリーから放電する電力を互いに融通することが できる,という状況を想定する. そのため,バッテリーから放 電する電力,他のビルから融通される電力,グリッドから購入 する電力,太陽光発電から供給される電力によって消費電力 を賄う. 消費電力量に関しては,快適性(室温)を考慮しつつ 消費電力量を抑えるようなスケジューリングを行う. なお,こ の問題において以下の仮定が全て成り立つものとする. ［仮定1］ 1.現時刻までの発電量T Gは既知である． 2. その日における発電予測値Pと需要量T Lは全て既知で ある． 3. バッテリー残量の初期値H0は既知であり,各時刻におけ る残量はデータとして入手可能である. 3. 2 階層型エネルギーマネージメント この節では,本研究で扱う階層型エネルギーマネージメント のモデルについて説明する. 以下に,システムモデルを示す. なお,このモデルは文献11)を参考にした.: Information communication : Power flow L-BEMS L-BEMS L-BEMS G-BEMS Energy scheduling Global Local Building PV generator Battery Building PV generator Battery Building PV generator Battery
Grid
Fig. 4 Concenptional architecture of BEMS
Fig. 4に示すように, L-BEMS, G-BEMSの二つのシステ
ムが存在する階層型BEMSを考える. まず, L-BEMSでは 各ビルにおいて快適性を考慮した消費電力量を各々決定する. 次にG-BEMSでは,全てのビルにおいてバッテリー容量制約 を満たしつつ総購入価格を抑えるような最適購入量wgird,融 通量wtradeを決定する. 3. 2. 1 L-BEMS L-BEMSでは,各ビルにおいて快適性を考慮した消費電力 量を決定する.本論文では,各ビルの快適性として部屋の温度 のみを考え,室温が所望の温度制約内に収まるようにしつつ, 消費電力を削減することを考える. まず,部屋の温度を考慮す るにあたり,室温に関するダイナミクスは式(12),式(13)の ように表す. CdT dt = q c +T oa− T R （12） qc= cpm(Tla− T ) （13） また,各記号について以下のTable 2にまとめる. Table 2 Symbols C[kJ/K] thermal capacity qc[kW] coolability T [◦C] room temperature Toa[◦C] outside temperature Tla[◦C] supply air temperature
R[K/kW] thermal resistance cp[kJ/kg K] specific heat capacity
m[kg/s] mass flow rate
次に,後退差分法を用いて離散化することで,以下の式(14) が得られる. Tt+1 = ATt+ BTtla+ DT oa t （14） A = 1− Ts RC − Ts Ccpm, B = Ts Ccpm, D = Ts RC 冷房による冷却能力qcと冷房による消費電力Upowerとの 関係は,成績係数COP(Coefficient of Performance)によっ て表現される15) . COP = q c t Utpower （15） 業務用エアコンのCOPは一般的に3.0∼3.5であることから, 本論文ではCOP = 3.2とする. 次に, L-BEMSで用いる評価関数を以下に表す. min Tla t J = Hp ∑ t=1 {∥ Tt+1− TF∥2 Q1 +∥ Ttla− T oa∥2 R1} （16） s.t. Tt+1= ATt+ BTtla+ DTtoa Tlb≤ Tt≤ Tub Tla,lb≤ Ttla≤ T la,ub （17） ここで, Tlb, Tub, Tla,lb, Tla,ubは室温,冷房による供給空 気温度の下限値および上限値とし, TF は室温の目標値とす る. そして, Hpは予測ホライズン, Q1, R1は各々状態,入力 に関する重みである. 空調の供給温度Ttlaを評価関数Jの最 小化問題として求めている. 予測ホライゾンについては, 1日の運用時間をTとすれば， 時刻はt = 1,· · · , T となり，予測ホライゾンNpは運用終 了時間までにするためN p = T− t + 1とした．今回想定し ているサンプリング時間は1時間であり，6時から18時の 12時間の運用を考えているため，T = 12となる. そして, L-BEMSで算出された消費電力を,空調によって1部屋で消 費される電力としてG-BEMSに用いる. 3. 2. 2 G-BEMS G-BEMSでは, L-BEMSで決定された各ビルの消費電力 T L,太陽光発電予測値P ,バッテリー残量Hから1step先の 最適な電力購入量wgird,融通量wtradeを算出する. まずは, 1stepずつ解く問題(予測ホライゾンNp= 1)とすると,以下 のようになる． min
wgridi,t+1,wtrade
i,t+1 Jt+1= N ∑ i=1 {Ri,t+1wgrid i,t+1 2
+ Ai,t+1( ˆHi,t+1− Hi,t+1F )
2} （
18） s.t. T Li,t+1 = wi,t+1grid + wi,t+1trade
+ T Gi,t+1+ ∆ ˆHi,t+1 （19） ˆ
Hi,t+1 ∈ [Hmin, Hmax] （20）
wi,t+1grid ≥ 0 （21） n ∑ i=1 wtradei,t+1 = 0 （22） 式(18)において,電力購入量wgrid,電力融通量wtradeを 評価関数Jの最小化問題として求める. 時刻tで次の時刻 t + 1での電力購入量wt+1grid,電力融通量w trade t+1 を求めるた め,評価関数における時刻のインデックスはt + 1となってい る.ここで, ˆHをバッテリー残量予測値, HFをバッテリー残 量目標値とし, Rは購入量に関する重み, Aをバッテリーに 関する重みとする. 制約式については,式(19)は需要T Lの 制約,式(20)はバッテリー容量に関する制約である. また,式 (21)は電力購入量が常に正であることを表し,式(22)は各ビ
計測自動制御学会論文集 第 xx 巻 第 xx 号 2019 年 xx 月 5
ルが他のビルから融通を受ける電力量の合計が0,つまり電力
融通ネットワークが完結していることを表している. 考慮す
るビルの台数はnとし, iはビルのインデックスを表す. また,放電推定量∆ ˆHは以下のように表される.
∆ ˆHi,t+1= Hi,t− ˆHi,t+1 （23）
ここで,時刻t + 1でのバッテリー残量は未知であるため予 測値Hi,t+1ˆ を用いる. また,時刻tにおいてはt + 1での実 際の発電量についても未知であるため,発電量T Gi,t+1は点 予測値Pi,t+1を用いる. T Gi,t+1= Pi,t+1 （24） ここで，ビルiが得る電力量wtotal i,t+1を購入量,融通量の合 計値として
wi,t+1total = wi,t+1grid + w trade i,t+1 （25） とすると,式(19)，式(20)より wi,t+1total ∈ [w total,min i,t+1 , w total,max i,t+1 ] （26）
wtotal,mini,t+1 = T Li,t+1− Pi,t+1− Hi,t+ Hmin wi,t+1total,max = T Li,t+1− Pi,t+1− Hi,t+ Hmax となる. よって,式(26)の制約の中でこの問題を解くことで, 時刻t + 1での最適電力購入量wgridi,t+1,最適電力融通量w trade i,t+1 を算出し,これらの値で運用する. その後,実際に時刻がt + 1 になると発電量T Gt+1が既知となるので,式(19)に代入す ることで以下のようにして時刻t + 1での実際のバッテリー 残量Hi,t+1を式(27),式(28)から求めることができる. T Li,t+1= wi,t+1total+ T Gi,t+1+ ∆Hi,t+1 ∴ ∆Hi,t+1= T Li,t+1− T Gi,t+1− wtotali,t+1 Hi,t+1in = Hi,t− T Li,t+1+ T Gi,t+1+ wtotali,t+1 （27）
Hi,t+1= Hi,t+1in (H in i,t+1∈ [H min , Hmax]) Hmax (Hin i,t+1> Hmax) Hmin (Hin i,t+1< Hmin) （28） ここで,式(27)は,実際の値を用いて求められた時刻t + 1 でのバッテリー残量Hin i,t+1を表している. しかし,この問題 を解く際は実際の発電量T Gを発電量予測値P で計算して いるが,実測値T Gと予測値P には誤差があるため,以下の 式よりバッテリー残量推定値Hˆ と実際のバッテリー残量H に差が生じうることが分かる. ˆ
Hi,t+1= Hi,t− T Li,t+1+ Pi,t+1+ wi,t+1 Hi,t+1= Hi,t− T Li,t+1+ T Gi,t+1+ wi,t+1
ˆ
Hi,t+1− Hi,t+1= Pi,t+1− T Gi,t+1 （29） そのため, ˆHt+1に関しては制約条件Ht+1∈ [Hmin, Hmax] が満たされていることが保証されるが,実際のバッテリー残 量Ht+1 に関して制約が満たされることは保証されていな い. そこで, Hin i,t+1 > Hmax の場合は余剰電力を放電し, Hi,t+1in < Hminの場合は不足電力を追加購入するものとする. 次に,以上の問題を予測ホライズンNp,制御ホライズン Nc= NpとしたMPCを用いた問題に書き換える. min wgridt+1,wtrade t+1 Jt+1= wt+1grid T Rt+1wgridt+1 +( ˆHt+1− Ht+1F ) T At+1( ˆHt+1− Ht+1F )（30） s.t. T Lt+1 = wgridt+1 + wtradet+1 + T Gt+1+ ∆ ˆHt+1 （31） ˆ H ∈ [Hmin, Hmax] （32） wgridt+1 ≥ 0 （33） n ∑ i=1 wtradei,t+1 = 0 （34） この問題は先ほどと同様,時刻tにおいて次の時刻t + 1 での最適解wt+1を見つける最適化問題である．wt+1は wgridt+1, wtrade t+1 を含む以下のような行列となっている. wt+1= [ wgridt+1 wtradet+1 ] （35） 他の各行列については以下に示す.
wgridt+1 = [w1,t+1grid · · · wgrid1,t+N p· · · wn,t+1grid · · · wgridn,t+N p]T
wtradet+1 = [w trade 1,t+1· · · w trade 1,t+N p· · · w trade n,t+1· · · w trade n,t+N p] T ˆ Ht+1 = [H1,t+1· · · H1,t+N p· · · Hn,t+1· · · Hn,t+N p]T ˆ Ht+1F = [H F 1,t+1· · · H F 1,t+N p· · · H F n,t+1· · · H F n,t+N p] T T Lt+1 = [T L1,t+1· · · T L1,t+N p · · · T Ln,t+1· · · T Ln,t+N p]T T Gt+1 = [T G1,t+1· · · T G1,t+N p · · · T Gn,t+1· · · T Gn,t+N p]T Pt+1 = [P1,t+1· · · P1,t+N p· · · Pn,t+1· · · Pn,t+N p]T
Ri = diag[Ri,t+1· · · Ri,t+N p] Rt+1 = diag[R1· · · Rn]
Ai = diag[Ai,t+1· · · Ai,t+N p] At+1 = diag[A1· · · An] 予測ホライゾンについては, L-BEMSの場合と同様,すべ ての時間においてその日の最終運用時間18時までを考慮す る. 次に,制約について考える.時刻t + kにおける制約につ いては式(26)と同様に wi,t+ktotal ∈ [w total,min i,t+1 , w total,max i,t+1 ] （36）
wi,t+ktotal,min = T Li,t+k− Pi,t+k− Hi,t+k−1+ Hmin wtotal,maxi,t+k = T Li,t+k− Pi,t+k− Hi,t+k−1+ Hmax
となる. k = 1では現時刻のバッテリー残量Htが取得可能
な情報であるため，次の時刻の制約条件は容易に求まる．し かし，k = 2,· · · , Npでは,将来のバッテリー残量Ht+k−1が
そこで以下のように考える．
wi,t+ktotal = T Li,t+k− Pi,t+k+ Hi,t+k− Hi,t+k−1 wi,t+ktotal−1 = T Li,t+k−1− Pi,t+k−1+ Hi,t+k−1− Hi,t+k−2
. . .
wtotali,t+1 = T Li,t+1− Pi,t+1+ Hi,t− Hi,t+1
これらをすべて足し合わせると, k ∑ l=1 wi,t+ltotal= k ∑ l=1 T Li,t+l− k ∑ l=1
Pi,t+l+ Hi,t+k− Hi,t （37） となる．Hi,t+k∈ [Hmin, Hmax]より，
k ∑
l=1
wtotalt+l ∈ [Sw min
i,t+k, Swi,t+kmax] （38）
Swmini,t+k = k ∑ l=1 T Li,t+l− k ∑ l=1
Pi,t+l− Hi,t+ Hmin Swmaxi,t+k = k ∑ l=1 T Li,t+l− k ∑ l=1
Pi,t+l− Hi,t+ Hmax
となり,時刻tのみのバッテリー残量のみで制約を求められ る．これを行列で表現すれば， Ewt+1≤ bt+1 （39） E = 1 0 · · · 0 1 1 · · · 0 . . . . .. . .. ... 1 1 · · · 1 −1 0 · · · 0 −1 −1 · · · 0 . . . . .. . .. ... −1 −1 · · · −1 bt+1 = Swmax1,t+1 · · · Sw max N,t+1 . . . · · · ... Swmax1,t+N p · · · Sw max N,t+N p −Swmin 1,t+1 · · · −SwminN,t+1 . . . · · · ... −Swmin 1,t+N p · · · −SwN,t+N pmin となり,この制約の下で問題を解く. そして,時刻がtから t + 1に更新され,実際の発電量T Gが分かった後，求められ た最適解ベクトルwの内, wgrid, wtradeの時刻t + 1の成 分，つまりwgridi,t+1,wtradei,t+1 (i = 1,· · · , N)を用いることで,式
(27),式(28)からバッテリー残量Hi,t+1が更新される. 3. 3 シナリオロバスト理論の適用9), 16) 次に,シナリオロバスト理論を前節の問題に適用する. 時刻 tにおけるt + 1の実際の発電量について,シナリオjに依存 した予測値Pδ(j) i,t+1を用いて以下のように考える. T Gi,t+1 = Pδ(j) i,t+1 （40） この予測値が各シナリオδ(j)によって変化することが前節 の内容と異なる．まず,評価関数,制約についてMPCを用い ない場合について整理すると以下のようになる. Jt+1 = N ∑ i=1 {Ri,t+1wgrid i,t+1 2
+ Ai,t+1(Xi,t+1− wgridi,t+1− w trade
i,t+1)2} （41） (Xi,t+1= Hi,t+ Pi,t+1− T Li,t+1− Hi,t+1F )
wtotali,t+1∈ [w total,min
i,t+1 (δ(j)), wi,t+1total,max(δ(j))] （42）
wi,t+1total,min(δ(j)) = T Li,t+1− P δ(j)
i,t+1− Hi,t+ H min
wtotal,maxi,t+1 (δ(j)) = T Li,t+1− P δ(j) i,t+1− Hi,t+ H max しかし,評価関数の値がシナリオによって変化してしまうと 問題を解くことができないため,上記のように評価関数にお ける発電予測値Pについては, Pδ(j)ではなくPとして計算 を行う. これによって,制約(式(42))のみがシナリオによっ て変化する問題となり,問題を解くことができる. また,発電量予測値Pi,t+1δ(j) については,文献 12)の内容であ る区間予測手法によって予測値に対する信頼区間CIα(αは信 頼度)を得ることができるため,その範囲内で一様な確率でラン ダムにN個決定したものをPδ(j) i,t+1(j = 1, 2,· · · , N)とする. また,発電量予測値Pi,t+1δ(j) については,文献12)の内容である 区間予測手法によって予測値に対する信頼区間CIα(αは信頼 度)を得ることができるため,その範囲内で区間予測値の確率 分布に従ってN個決定したものをPδ(j) i,t+1(j = 1, 2,· · · , N) とする. 尚,コピュラを用いて算出した区間予測値の確率分布 は正規分布であった.そして,得られたN個のPδ(j) i,t+1を用い
て計算された制約wtotal,mini,t+1 (δ(j)), wi,t+1total,max(δ(j))の内,最
も制約の厳しいものをこの問題における制約とする. 以上を
踏まえて, MPCを用いた場合では,制約は以下のようになる．
k ∑ l=1
wi,t+l∈ [Swmini,t+k, Swmaxi,t+k] S wi,t+kmin(δ(j)) =
T Li,t+1− Pi,t+1δ(j) − H0+ Hmin (k = 1)
∑k l=1T Li,t+l− (P δ(j) i,t+1+ ∑k l=2Pi,t+l) − H0+ Hmin (k = 2,· · · , Np) S wi,t+kmax(δ(j)) =
T Li,t+1− Pi,t+1δ(j) − H0+ Hmax (k = 1)
∑k l=1T Li,t+l− (P δ(j) i,t+1+ ∑k l=2Pi,t+l) − H0+ Hmax (k = 2,· · · , Np) 導出された最適解ベクトルwt+1grid, wt+1tradeはN pステップ 先までの解を求めているが,実際に用いる値は第1要素のみ であるため第2要素以降の制約違反を考える必要はない．そ のため，1時間先の制約のみをシナリオ化させている． 3. 4 シミュレーション 3. 4. 1 シミュレーション条件 L-BEMSで想定する外気温度ToaをFig. 5に示す.
計測自動制御学会論文集 第 xx 巻 第 xx 号 2019 年 xx 月 7
Fig. 5 Outside Temperature(Day 1)
G-BEMSで想定するビルの台数N = 3とし,各ビルの部 屋数Rを各々50, 40, 30,ビルの外気温度Toaは上記のよう に既知であり,外気温度,太陽光発電量T Gは全てのビルにお いて等しいものとする. 太陽光発電予測値Pについては,前章 における提案手法によって算出された予測値を用い, 1時間 当たりの最大出力60 W の太陽光パネルが200枚，126 W h のバッテリーが12個置かれている状況を想定する．また,ビ ルを想定しているため空調による消費電力の他に15 kW の 基本電力を常に使用していると想定し,各ビルの重要量T Lは 以下の式で表すものとする. T Li = Utpower× Ri+ 15 （43） 購入電力に関する重みは電力料金を基に作成し,東京電力 の料金プランを参考に価格設定を行った17)．以下のFig. 6 が各ビルにおける電力料金であり,重みRである. 6 7 8 9 10 11 12 13 14 15 16 17 18 Time [h] 0 10 20 30 40 50 60
Power Price [yen/kW]
Building1 Building2 Building3
Fig. 6 Weight of Electrical Rate R
その他のパラメータについては以下のTable 3に示す.
Table 3 G-BEMS Simulation parameters
Parmameter Symbol Value sampling time[h] st 1 hour room temperature upper limit[℃] Tub 28
room temperature lower limit[℃] Tlb 24
room temperature reference value[℃] TF 26
room thermal capacity[kJ/K] C 1600 room thermal resistance[K/kW] R 4.5 room capacity[m3] V 350 room specific heat capacity[kJ/kgK] cp 1.007
mass flow rate[kg/s] m 0.15 battery capacity initial value[kWh] H0 7.5 battery capacity maximum value[kWh] Hmax 15
battery capacity minimum value[kWh] Hmin 0
battery capacity reference value[kWh] HF 7.5
simulation period[h] T 12 the number of scenarios N 200 degree of reliability α 80%
シミュレーションに用いる手法については以下のようにする.
Table 4 Methods
Method1 Method2 Method3(Proposed Method)
Power Trade × 〇 〇 Scinario Robust × × 〇 Method1では,シナリオロバスト理論は適用せず,電力の 融通も考慮していない. Method2では,電力の融通を考慮す るが,シナリオロバスト理論を適用していないため制約違反 に対するロバスト性は低い. それらに対してMethod3(提案 手法)では,その両方を考慮することで,コストを削減しつつ ロバスト性を有する解を算出することを考える. 3. 4. 2 シミュレーション結果 まず, L-BEMS におけるシミュレーション結果を以下の Fig. 7∼Fig. 9に示す.
Fig. 8 Supply Temperature (Day1)
Fig. 9 Electricity Consumption (Day1)
Fig. 7より,室温が制約温度範囲内に収められていることが 確認できた. 消費電力,供給温度については,外気温度に影響 を受けるため外気温度の推移と似たような推移になっている. 次に, G-BEMSに関して各手法によるシミュレーション結 果として,総コストと不足電力を以下のTable 5にまとめた. 不足電力量は, G-BEMSにおいてスケジューリングを行う際, 発電予測誤差に対応できず不足した電力量である. そのため, 不足電力量を比較することで制約違反について評価すること ができる.
Table 5 G-BEMS Result
Total Cost[thousand yen] Shortage Power[kWh] Method1 478.67 19.5777 Method2 467.74 11.6350 Method3 468.46 1.6587 まず, Method1とMethod3を比較すると,ビル間での電力 融通により総コスト(Total Cost)を削減でき,シナリオロバ スト理論を用いたことにより不足電力(Shortage Power)を 削減できたことが分かる. 次に, Method2とMethod3を比 較すると,シナリオロバスト理論を用いたことにより不足電 力を削減することができたが,電力コストに関しては多少で はあるが増大してしまっている.これはロバスト性を考慮し たことによって最適性が失われるためであると考えられるた め,妥当な結果であると考えられる.
4.
お わ り に 本論文では,まず予測残差補正を行う太陽光発電予測手法を 提案した. 従来手法(k-平均法, SVRを用いた手法)から予測 値を算出し,予測値と実測値との予測残差データを収集するこ とで,予測残差の推定を行うことを考えた. ある時刻の予測残 差が得られた時次の時刻の予測残差を推定し,推定された予測 残差により予測値の修正を行った. 最後に,有効性について実 データを用いたシミュレーションによって確認した. BEMS においては,複数のビル間において電力の融通を図ることに よって全てのビルにおける総コストを削減するような階層型 エネルギーマネージメントを提案した. ビル同士で電力を融 通することで,各時刻において価格の低いビルでの電力購入量 を増やし,全体的なコスト削減を目指した. 太陽光発電によっ て得られた電力を使用するが,発電予測値の予測誤差により運 用計画と実際の運用に差が生まれてバッテリー制約を犯す可 能性があるため,シナリオロバスト理論に適用することでロバ スト性を持ったBEMSの運用を行った. BEMSに関しても, 最後に有効性についてシミュレーションによって確認した. 今後の課題としては,複数のビルでの電力のやり取りにお いて送電にかかるコストや電力損失などを考慮することでよ り現実的なモデルを扱うことや,より大規模なスマートコミュ ニティでの問題に発展させることが考えられる. 参 考 文 献 1）滑川徹: スマートグリッドのための分散予測制御, 計測と制御, 51-1, 62/68 (2012) 2）石川友規, 小嶋昂明, 滑川徹: JIT モデリングに基づくカルマン フィルタによる短期風力発電量予測, 電気学会論文誌 C, 135-1, 81/89 (2015)3）P.Xiaosheng, W.Bo, Y.Fan, F.Gaofeng, W.Zheng, C.Kai: A Deep Learning Approach for Wind Power Prediction Based on Stacked Denoising Auto Encoders Optimized by Bat Algorithm, Proc. of China International Conference on Electricity Distribution (CICED), 945/948 (2018) 4）Han Seung Jang, Kuk Yeol Bae, Hong-Shik Park, Dan
Keun Sung: Solar Power Prediction Based on Satellite Im-ages and Support Vector Machine, IEEE Transactions on Sustainable Energy, 7-3, 1255/1263 (2016)
5）A.Asrari, T.X.Wu, B.Ramos: A Hybrid Algorithm for Short-Term Solar Power Prediction-Sunshine State Case Study, IEEE Transactions on Sustainable Energy, 8-2, 582/591 (2017)
6）Can Wan, Jianhui Wang, Jin Lin, Yonghua Song, Zhao Yang Dong: Nonparametric Prediction Intervals of Wind Power via Linear Programming, IEEE Transactions on Power Systems, 33-1, 1074/1076 (2018)
7）Mohammad Javad Sanjari, H. B. Gooi: Probabilistic Fore-cast of PV Power Generation Based on Higher Order Markov Chain, IEEE Transactions on Power Systems, 32-4, 2942/2952 (2017)
8）竹田皓貴, 鷹羽浄嗣: PV 発電予測に基づくマイクログリッ ドの確率的モデル予測制御, 計測自動制御学会論文誌, 54-2, 219/226 (2018)
9）H.V.Haghi, S.Lotfard, Z.Qu: Multivariate Predictive An-alytics of Wind Power Data for Robust Control of En-ergy Storage, IEEE Transactions on Industrial
Informat-計測自動制御学会論文集 第 xx 巻 第 xx 号 2019 年 xx 月 9 ics, 12-4, 1350/1360 (2016) 10）相澤亮太, 熊谷敏, 杵嶋修三: 仮想熱源ストレージへの熱資源 貯蓄と再配分による空調運転計画-居室の熱慣性応答推定と熱 資源配分アルゴリズム-, 電気学会論文誌 C, 136-2, 233/243 (2017)
11）Il-Young Joo, Dae-Hyun Choi: Distributed optimiza-tion framework for energy management of multiple smart homes with distributed energy resources, IEEE Access, 5, 15551/15560 (2017)
12）Shotaro Sato, Toru Namerikawa: Scenario-Based Robust MPC for Energy Management Systems with Renewable Generators, Proc. of 37th Chinese Control Conference (CCC), 2304/2309 (2018)
13）Kun Ding, Li Feng, Xiang Wang, Siyu Qin: Forecast of PV Power Generation Based on Residual Correction of Markov Chain, Proc. of International Conference on Control, Au-tomation and Information Sciences (ICCAIS), 355/359 (2015)
14）塚原英敦: 21 世紀の統計科学 III, 360, 東京大学出版会 (2008) 15）JSRAE 日本冷凍空調学会 成績係数:
https://www.jsrae.or.jp/annai/yougo/28.html, 2019/11/19 access.
16）M.C.Campi, S.Garatti, F.A.Ramponi: Non-convex sce-nario optimization with application to system identifica-tion, Proc.of 54th IEEE Conference on Decision and Con-trol (CDC), 4023/4028 (2015) 17）東京電力エナジーパートナー 料金プラン: http://www.tepco.co.jp/ep/private/plan/index-j.html, 2019/7/31 access. ［著 者 紹 介］ 山 内 賢 2019 年慶應義塾大学理工学部システムデザイン 工学科卒業．同年同大学大学院理工学研究科総合 デザイン工学専攻に入学. 滑 川 徹（正会員） 1994 年金沢大学大学院自然科学研究科システ ム科学専攻博士課程中退．同年金沢大学工学部電 気・情報工学科助手．同講師を経て 2002 年長岡 技術科学大学機械系助教授．2006 年金沢大学大学 院自然科学研究科電子情報科学専攻助教授．2009 年慶應義塾大学理工学部システムデザイン工学科 准教授，2014 年同教授となり現在に至る．ロバス ト制御理論，分散協調制御理論とその電力ネット ワークシステムへの応用に関する研究に従事．博 士 (工学)．2014 年計測自動制御学会制御部門パ イオニア技術賞，2017 年計測自動制御学会論文賞 を受賞．IEEE CSS，システム制御情報学会，電 気学会などの会員．
