非線形システムの最適化

発電プラントの最適運転

発電プラントの最適運転とは、電力需要を賄いながら、システム全体の損失が最も小さくなる発電プラントの稼働スケジュールを求めることになります。

最適化問題の定式化

発電プラントの特性は非線形です。この最適化問題を逐次2次計画法で解くため、問題を定式化していきます。

まずは、

• 発電プラントの発電電力　\(Pg\)
• 揚水式発電プラントの発電電力  \(Ph\)
• 揚水式発電プラントが水をポンプアップ（揚水）する時の電力  \(Pp\)

とします。

供給される電力\(Ps\)は   \(Ps ＝ Pg + Ph - Pp\)  です。

発電プラントの最適化問題は、電力需要 \(Ps\) を賄いながら、システム全体の損失が最も小さくなる\(Pg,Ph,Pp\)の稼働スケジュールを求めることになります。

稼働時間を分割して離散化

この最適化問題を連続時間で考えるのは難しいので、稼働時間をK分割して離散化して扱うことにします。ある時点 \(i\) ( \(i\) :1～K)での供給される電力は \[Ps_i ＝ Pg_i + Ph_i - Pp_i\]で表します。

既に、最初の電力需要で示しているように、24時間を3時間ごとに8分割して考えることにします。(K=8)

3時間フェーズ毎の電力需要（供給電力[Mkw]）は、時刻0時から、3時、6時、9時、12時、15時、18時、21時の順に、

　Ps[K] = { 40, 40, 60, 80,100, 90, 70, 60 }

と定義します。

目的変数のベクトルX[N]

発電プラントの最適化問題は、システム全体の損失が最も小さくなる\(Pg,Ph,Pp\)の稼働スケジュールを求めることなので、目的変数のベクトルX[N]は、3時間毎に8フェーズに分割した各発電電力 \(Pg_i,Ph_i,Pp_i\)になります。\(i\) は1～8のフェーズを表します。

目的変数ベクトルX[N]の並びは、\(Pg_i,Ph_i,Pp_i\) の順に、

　X[N] = { Pg1,・・・, Pg8,  Ph1,・・・Ph8,  Pp1,・・・, Pp8 } 

とします。 (N=24)

尚、目的変数のベクトルX[N]の初期値は、電力需要に応じた発電として、

   X0[N] = {　40, 40, 60, 80,100, 90, 70, 60,
               0,  0,  0,  0,  0,  0,  0,  0,
               0,  0,  0,  0,  0,  0,  0,  0　}
としています。

目的関数

ある時刻フェーズ \(i\) の発電プラントの損失は、\[LossPg_i = 0.000213599*Pg_i^3 - 0.030697784*Pg_i^2 + 1.766374137*Pg_i + 25.35128198\]

揚水発電と揚水蓄電の損失は、\[LossPh_i = 0.10*Ph_i\] \[LossPp_i = 0.10*Pp_i\]
目的関数は、発電プラントの損失と揚水発電の損失、揚水蓄電の損失の積算値になります。
\[OBJ = \sum_{i=1}^{8}(LossPg_i +LossPh_i+LossPp_i)\]
逐次2次計画法では、損失の積算値 \(OBJ\) を最小にする目的変数ベクトルX[N]を求めます。

制約関数c(x)

制約関数 \(c(x)\) は下記の64式。( \(i\) :1～8) 

各時刻フェーズごとの電力収支
(1)～(8)      \(Pg_i+Ph_i-Pp_i-Ps_i = 0\)

発電プラントの出力上限
(9)～(16)   \(-Pg_i+100[Mkw] ≧0\)

揚水発電の出力上限
(17)～(24)   \(-Ph_i+20[Mkw] ≧0\)

揚水蓄電の電力上限
(25)～(32)   \(-Pp_i+20[Mkw] ≧0\)

各電力は正の値
(33)～(40)   \(Pg_i ≧0\)
(41)～(48)   \(Ph_i ≧0\)
(49)～(56)   \(Pp_i ≧0\)

揚水発電と揚水蓄電の電力収支。ポンプアップ（揚水）した分しか発電できない。
(57)～(64) \[0.90*\sum_{i=1}^{1}(Pp_i)-1.10*\sum_{i=1}^{1}(Ph_i)≧0\]\[0.90*\sum_{i=1}^{2}(Pp_i)-1.10*\sum_{i=1}^{2}(Ph_i)≧0\]\[0.90*\sum_{i=1}^{3}(Pp_i)-1.10*\sum_{i=1}^{3}(Ph_i)≧0\]\[0.90*\sum_{i=1}^{4}(Pp_i)-1.10*\sum_{i=1}^{4}(Ph_i)≧0\]\[0.90*\sum_{i=1}^{5}(Pp_i)-1.10*\sum_{i=1}^{5}(Ph_i)≧0\]\[0.90*\sum_{i=1}^{6}(Pp_i)-1.10*\sum_{i=1}^{6}(Ph_i)≧0\]\[0.90*\sum_{i=1}^{7}(Pp_i)-1.10*\sum_{i=1}^{7}(Ph_i)≧0\]\[0.90*\sum_{i=1}^{8}(Pp_i)-1.10*\sum_{i=1}^{8}(Ph_i)≧0\]

目的関数f(x)の勾配ベクトル∇f(x)

目的関数 \(f(x)\) の勾配ベクトル \(∇f(x)\) は

発電プラントの損失

　\( f(x) = 0.000213599*Pg^3 - 0.030697784*Pg^2 + 1.766374137*Pg + 25.35128198\)

より、勾配ベクトル∇f(x)は

   \(∇f(x) = 0.000640797*Pg^2 - 0.061395568*Pg + 1.766374137\)

揚水発電の損失勾配ベクトル = 0.10 (10%)
揚水蓄電の損失勾配ベクトル = 0.10 (10%)

制約関数c(x)の勾配ベクトル∇c(x)

(1)～(64)の制約関数 \(c(x)\) の勾配 \(∇c(x)\)（偏微分係数 \(∂Cm/∂Xn\) ）をCONGRAT[M][N]に格納しています。

ダウンロード

上記の発電プラントの最適化問題をプログラムSQPに実装し、Windows VisualStudio(x64)で計算しています。

発電プラント最適化問題のCソースコードはこちらからダウンロードできます→リンク

計算結果

VisualStudio(x64)で計算した結果、演算回数25回で、こうなりました。
制約式条件も全て満たしています。

やはり夜間電力で揚水発電の水を汲み上げています。
この時、発電プラントの出力を40[Mkw]から60[Mkw]に上げており、より効率の良い運転ポイントで発電プラントを稼働させています。

さらに、電力需要のピークで揚水発電を実施することで、電力需要ピーク時の発電プラントの出力を最大効率点の80[Mkw]とし、ここで大幅に損失を低減しています。

揚水発電を稼働させることで、発電プラントの運転を全区間において最大効率点の80[Mkw]に近づける稼働スケジュールになっています。

発電システム全体の損失エネルギーの累積比較の時間推移を見ると、夜間電力でポンプアップ（揚水）している時間帯に発生する損失エネルギーの追加分を、電力需要ピーク時間帯の揚水発電アシストでカバーしつつ、さらに損失エネルギーを削減できていることから、高負荷運転時の発電効率の改善がポイントであることが分かります。

稼働スケジュールの最適化により、電力の需要を満たしながら、発電プラントシステムの損失を 71.4[Mkw/h] (3.92%) 低減しています。

この損失エネルギーの削減量は、発電プラントを1時間程度休止するのと同等のエネルギー量にあたります。

終わりに

逐次2次計画法のソフトウエアSQPに、目的変数24、制約条件64式の非線形最適化問題を実装してみました。

実際の電力インフラの運転は、特性の異なる複数の発電プラントを電力需要に応じて運転するわけで、今回前提とした単純なモデルのようには扱えないと思います。また、風力発電や太陽光発電など、予測が難しく変動も大きいエネルギー源もあり、揚水式水力発電所の運用目的もより幅広く考える必要があると思います。

最後までお読みいただき、ありがとうございました。
何かしらのお役に立てば嬉しいかぎりです。

次は、このSQPプログラムをDLL化してEXCELのソルバーにしてみます → リンク