発電所はオペレ，ションズ，リサ，チで支えられている

发布的井上渡，2022年10月24日

452次观看(过去30天)|0喜欢|0评论

           今日はアプリケ，ションエンジニアの岩本さんによる投稿です。

           昨今,温暖化問題やロシアのウクライナ侵攻に端を発するエネルギー情勢の変動により,国内外の電力供給に不安が生じています。

           かねてから電力会社は,発電から送電に関わる運営に対する一括した責任によって,電源開発や系統設備の導入を継続的に行って,安定的な電力の供給を担ってました。しかし,3.11の震災を機に,それ以前から議論のあった電力自由化や総括原価方式の見直しが急速に進み,2016年からオープンな電力市場が形成されました。

           さらに一部の電力会社は，発送電分離化され，発電，送電の会社へと分社化されました当初,これらは電気料金を下げる目的で政策が施行されましたが,再生可能エネルギーの大量導入と固定価格買取制度(适合)による需要者負担増(いわゆる再エネ賦課金),液化天然气などの燃料高騰等によってここ数年の電気料金は増加しています。

引用:資源エネルギ，庁日本のエネルギ，2021年度版"エネルギ，の今を知る10の質問"

            また電力需給の面でもかなり厳しい状態になっています。要因の一は系統への再エネ導入によって，周波数が不安定化するためです。

引用:関西電力送配電でんき予報より

            上の図は，関西地域におけるとある日の電力の実績です。電力のトレンドは各社がでんき予報として公開しています。

            図中の黄色の曲線は太陽光の発電実績になりますが,見ての通り,太陽が真上にくる正午をピークとして,小高い丘のような形状をしています。日が落ちる夕方から夜間の発電実績は0(千瓦)です。このように再エネは気象現象によって発電できない時間帯が発生します。(注記:再エネの中でも水力や地熱は比較して安定的な電源です)

            我々が消費する電気は，今消費している電力とピッタシ同じ量だけ発電することで正常に回っています。少しでも需要と供給に過不足が発生すると,系統の周波数(東は50 Hz,西は60 (Hz))がだんだんと乱れ,やがては発電機が耐え切れずに系統から続々と解列していき,ブラックアウト(大規模停電)します。ブラックアウトと言えば，北海道の例が記憶に新しいですね。（参考：資源エネルギ，庁hpより日本初の“ブラックアウト”，その時一体何が起きたのか）

            この需要と供給が完全に一致することを同時同量と呼びますが,再エネが出力をバンバン振ってしまうと系統が不安定化しやすくなります。そのため,不測の事態に備えるために,調整用の火力発電が出力を下げて,常時運転をキープしたり,細かな起動停止を行う必要が出てきます。ただし,火力発電は定格負荷で熱効率が最大になるように設計されているので,部分負荷帯で運転をすると,効率が悪くなります。

            また燃料の在庫制約を考慮したり，発電の経済性を考えたりと，各発電所の運営を考えるのは至難の業です。実際,電力会社によっては100を超える発電所の運営を行っているケースもあるので,その一つ一つの計画を人間が考えるとなると作業の膨大さが容易に想像できるかと思います。

            こうした発電所の運営にまつわる問題背景に対し,従来から或(オペレーションズ・リサーチ)の技術が活用されてきました。ここで，OR とは数学や統計モデル、最適化手法を活用し、経営や運用など実社会における様々な問題に対し、ある指標に対し最も効率的になるような意思決定を行う科学的アプローチの総称です。（参考：電力システム改革に対応した火力発電計画作成システムの最適化計算手法の開発）

            実際に上記論文では,最新の或システムの導入によって,燃料費ベースで年間0.48%もの経済性を改善したとの報告が上がっています。仮に年間数兆円規模の燃料調達費がかかると想定すると，数十億円の規模でコストカットしている計算です。さらに環境問題となる二氧化碳排出量もモデル化し,或で解く問題に取り込むことで,二氧化碳排出量が最も少ない発電所の運営を立案できます。これは，数学(数学)の力を活用することで実社会に対し大きな貢献ができることを示す良い例かと思います。

発電機の起動停止計画問題(单位承诺)

            では,実際に火力発電機群の起動停止問題を解くといった場合にどうするかというと,最適化問題に定式化してソルバーによって解かせます。MATLABの製品ファミリ，である优化工具箱™には，発電機のスケジュ，リング問題の例題があるので,これをベースにちょこっと改良して,簡単な発電所のスケジューリング問題を解かせてみることにしました。

            例題として6台の火力ユニット(石2基炭,LNG3基,石油1基)の構成で,1日の需要予測に対し,燃料費が最も安くなる構成パターンを組み合わせ最適化問題によって解かせます。なお,問題では1日24 hを30分毎に刻んだ時断面(メッシュ)について,計48メッシュにおけるユニット構成を決めていきます。

            上の図は,火力ユニットに対する電力需要(MW)を想定しており,昼間に太陽光がバンバン発電して,火力の需要が一時減ってしまうものの,夕方にかけてから太陽光の出力が減少し,再び火力の需要が一気に増加,ピークで使用率100%となるような厳しいパターンを想定しています。(通称ダックカ，ブ，参考:太陽光発電の増加に伴う課題“ダックカ，ブ現象とは”）

            問題では,この需要に対し,各種運転制約(後述)を考慮しながら,同時同量の原則で最も経済的な運転パターンを決めます。

            続いて，各発電機の特性は下表のとおりとします.(数値はダミ)

            項目概要：

a, b, c:兆瓦(発電出力)に対する燃料費特性(￥)※2次、関数(￥= ax ^ 2 + bx + c)で近似
P_max, p_min:発電機が出力できるmwの上下限(mw)
启动费用:発電機が停止から起動する際にかかるコスト(￥)
燃料:燃料のタ▪プ※石炭，LNG，石油の3種LNGはGTCCを想定
Minimum_Operation:一度起動したら運転を継続しなければならない最小時間(メッシュ)
最小化停止:一度停止したら停止を継続しなければならない最小時間(メッシュ)

            さて,発電機の起動停止計画問題は,どの発電機を選択するかを0 - 1の2値による表現,すなわち整数で表すため,整数と発電出力兆瓦(実数)の両方を最適化の決定変数に含む,混合整数計画問題(混合整数规划:MIP)と呼ばれる，解くことが難しい問題(np困難)に帰着されます。

            优化工具箱™では，このMIPに対応するソルバ，として混合整数線形計画問題(混合整数线性规划:MILP)を解くintlinprogがあるので，これを使って問題を解かせてみようと思います。

            詳細には,MILPで解くと時間ががかることが想定されるので,グリーディ法によってユニットの仮決め(初期推定解)を行ってから,MILPによって最適なユニットの構成を決め,最後にQP(2次計画問題)によって具体的なMWを決めるという段階的なアプローチで時間の短縮を図ろうと思います。

            Qpの解法ではquadprogというソルバ，を使います。

milpによる並解列の計算

            MILPでは1メッシュ毎(30分毎)の電力需要に対する経済的なユニットの運転状態(離散値)を決めていきます。

            運転状態としては,大まかに最低出力(P_MIN)中間出力((P_MAX-P_MIN) / 2),最大出力(P_MAX)の3つを用意します。各運転状態は0-1の設計変数として用意し，ソルバ，によって選択させます。なお，优化工具箱™では，最適化問題のプログラミングにおいて問題ベ，スのアプロ，チを行うことが可能となっています。

            問題ベースのプログラミングとは,要するにシンボリックな変数名(とかxとかの記述)で最適化の定式をそのままプログラム化できるということです。たとえば,運転状態を(低,中,高)yという0-1の決定変数に当てはめる場合には，以下のように定義します。

              nPeriods = 48;总目数百分比
             

              nGens = 6;总单位数%
             

              Y = optimvar(“y”nPeriods ngen, {“低”，“中间”，“高”}，“类型”，“整数”，下界的0,.．.
             

              “UpperBound”1);
             

            各ユニットが一つのメッシュで取り得る運転状態は1つであるため,以下のようにすべての総和が1以下になるという制約を考慮します。

            なお，場合によっては運転状態のどれも選択しない，即，停止もこの制約によって加味できます。

              Powercons = y(:，:，“低”) + y(:，:，“中间”) + y(:，:，“高”) <= 1;
             

            続いて，各メッシュにおけるmwの需要に対し，必要な供給量が満足されるように不等式制約を設けます。ここで，我がメッシュ，jがユニット，kが運転状態の添え字(低から順に1 ~ 3),创が美元(i, j, k)美元において取り得る出力(MW)を表します。

            美元MW_i \ leq \ sum_ {k = 1} ^ {3} \ sum_ {j = 1} ^ {ngen}创(i, j, k) + s_i美元，$ i = 1、2、3…,48美元，$ nGens = 6 $
           

            上式において厳密な等号としないのは,発電機の稼働を運転状態という離散化された有限のパターンの内から選択するようにしているため,厳密な出力のマッチングが保証されないためです。また不等式右辺に年代という変数がありますが，これは不等式を緩和するために設けたスラック変数(非負の実数)です。

            不等式を緩和する意図は,積み上げた供給力が足りないと制約が満たせなくなり,実行可能解がないという状態に陥ってしまいます。実行可能解が見つからないとしてもどこが厳しいのかを把握するためにも解を出力できる状態にしておきたいと思います。なお，スラック変数もソルバ，によって選択させる決定変数なので以下のように定義します。

              S = optimvar(“年代”nPeriods,“类型”，“连续”，下界的, 0);松弛变量%
             

              loadcons = optimineq(nPeriods,1);
             

              为i = 1:nPeriods
             

              Loadcons (i) = sum(gen(:，1) ' .*y(i，:，“低”) + gen(:，2) ' .*y(i，:，“中间”) + gen(:，3) ' .*y(i，:，“高”)) + s(i) >= MW(i);
             

              结束
             

            またソルバ，が発電機の停止と起動をむらみやたらに選択させないように考慮します。ここで，運転状態の変数である $y(i, j, k)$ のオン１とオフ０の期間に応じて決まる補助変数zを導入します。

              Z = optimvar(“z”nPeriods ngen,“类型”，“整数”，下界的0,.．.
             

              “UpperBound”1);
             

            ソルバ，が変数zを設定できるようにするには，次のように定式化します。

            ·条件$z(i,j) = 1$が厳密に$\sum_k y(i,j,k) = 0と$\sum_k y(i+1,j,k) = 1という条件下において満たされる必要がある

            ・$z(i,j) = 1$が厳密に満たされる場合に，$\sum_k (y(i+1,j,k) - y(i,j,k)) > 0$となる

            したがって，問題の定式化として次の線形不等式制約を含めます。

            $ \sum_k (y(i+1,j,k) - y(i,j,k)) - z(i,j) < = 0 $
           

            また，zは最適化の評価関数(後述)に含めるため，最小化問題を前提とするとソルバ，はzに対し，その値を低くしようとします。つまり,ソルバーはすべてのz(停止から起動への変化)を0にしようとします。

            しかし，需要増加によって発電機の起動が迫られるタaaplミングでは，線形不等式によって$ $z(i,j)$ $が1に等しくなるよう強制されるという具合になります。ここで利便性のために，$ $y(i+1,j,k) - y(i,j,k)$ $を表す補助変数wを作成しておきます。

              w = optimexpr(nPeriods,nGens);%分配w
             

              idx = 1:(nPeriods-1);% nPeriods = 48
             

              W (idx，:) = y(idx+1，:，“低”) - y(idx，:，“低”) + y(idx+1，:，“中间”) - y(idx，:，“中间”) + y(idx+1，:，“高”) - y(idx，:，“高”）;
             

              w(nPeriods，:) = y(1，:，“低”) - y(nPeriods，:，“低”) + y(1，:，“中间”) - y(nPeriods，:，“中间”) + y(1，:，“高”) - y(nPeriods，:，“高”）;
             

              开关= w - z <= 0;
             

            あとは最小起動時間と最小停止時間の制約を設定します。

            まず最小起動時間ですが,前述のとおり一度発電機が起動したら継続しなければならない最小の運転時間になります。定式化すると次式の通りです。

            $ \ sum_ {k = 1} ^ {3} (y_ {ijk} -y_ {i-1jk}) \ leq \ sum_ {k = 1} ^ {3} (y_{\τjk})美元，$ \tau = i+1，.....min(i+ l_j -1,48) $，$ I = 2,3，…，48 $，$ j = 1,2,3，…，nGens(=6) $
           

            上式においてlは各ユニットの最小起動時間(メッシュ)です。

              minOpecons = optimineq(nPeriods-1,nGens,max(minOperation));% minOperation是L的向量
             

              为i = 2:nPeriods% nperots = 48
             

              为j = 1:nGens% nGens = 6
             

              Tau = i;
             

              为l = 1:min(i+minOperation(j)-1,nPeriods)-(i+1)+1
             

              Tau = Tau + 1;
             

              minOpecons(i-1,j,l) = y(i,j，“低”- y(i-1,j，“低”) + y(i,j，“中间”- y(i-1,j，“中间”) + y(i,j，“高”- y(i-1,j，“高”) <= y(tau,j，“低”) + y(，j，“中间”) + y(，j，“高”）;
             

              结束
             

              结束
             

              结束
             

            同じ原理で，最小停止時間の制約も定式化します。最小停止時間は一度停止したら停止を継続しなければならない最小の時間です。

            $ \ sum_ {k = 1} ^ {3} (y_ {i-1jk} -y_ {ijk}) \ leq \ sum_ {k = 1} ^ {3} (1-y_{\τjk})美元，$ \tau = i+1，.....min(i+l_j-1,48) $，$ I = 2,3，…，48 $，$ j = 1,2,3，…，nGens(=6) $
           

              minStopcons = optimineq(nPeriods-1,nGens,max(minStop));% minStop是l的向量
             

              为i = 2:nPeriods
             

              为j = 1:nGens
             

              Tau = i;
             

              为l = 1:min(i+minOperation(j)-1,nPeriods)-(i+1)+1
             

              Tau = Tau + 1;
             

              minStopcons(i-1,j,l) = y(i-1,j，“低”) - y(i,j，“低”+ y(i-1,j，“中间”) - y(i,j，“中间”+ y(i-1,j，“高”) - y(i,j，“高”) <= (1-y(tau,j，“低”) + (1-y(tau,j，“中间”) + (1-y(tau,j，“高”));
             

              结束
             

              结束
             

              结束
             

            制約は以上で全てです。

            最後に最適化の評価関数を定式化します。

            評価関数は,発電機の運転時の燃料費,停止中の発電機を起動する際に生じる起動損失,スラック変数から構成します。

            $ J = \ sum_ {k = 1} ^ {3} \ sum_ {J = 1} ^ {ngen} \ sum_ {i = 1} ^ f {jk} {nPeriods} y (i, J, k) + \ sum_ {J = 1} ^ {ngen} \ sum_ {i = 1} ^ {nPeriods} h_jz (i, J) + \ sum_ {i = 1} ^ {nPeriods} s_i美元
           

            上式においてfは運転状態に応じた燃料費の特性，hは起動損失($ $z=1$ $の時に掛かってくる)です。

            これをコ，ドで書くと以下の通りです。なお，コ，ド上では，コストの正規化処理とスラック変数を小さくするための重みをそれぞれ掛けています。

              fuel = fuelArray.*y;
             

              startingCost = z*startCost;
             

              scaleFactor = 1e+8;
             

              Weight_slack = 1e+6;
             

              totalCost = (sum(sum(sum(sum(fuel))) + sum(sum(startingCost)))/scaleFactor + weight_slack*sum(s);
             

            以上でmilpの定式化が整ったので，これをソルバに食わせて問題を解いてくれるのを待ます。

            最適化問題を設定するために専用のapiであるoptimproblemを呼び出してオブジェクトを作成し，これまで作った制約条件と評価関数を指定していきます。

              unitCommit = optimproblem(“ObjectiveSense”，“最小化”）;
             

              unitCommit。目标=总成本;
             

              unitCommit.Constraints.switchcons = switchcons;
             

              unitCommit.Constraints.powercons = powercons;
             

              unitCommit.Constraints.loadcons = loadcons;
             

              unitCommit.Constraints.minOpecons = minOpecons;
             

              unitCommit.Constraints.minStopcons = minStopcons;
             

            最後にintlinprogのオプションを諸々設定して，解决コマンドで解きます。ここで，年代olve コマンドの引数に与えている unitCommit0 は初期推定解で、事前にグリーディ法で求めています。（グリーディ法は割愛）

            なみに初期推定解を空欄([])にしても解かせることは可能です。

              选项= optimoptions(“intlinprog”，“显示”，“通路”，“CutMaxIterations”, 50岁,“CutGeneration”，“高级”，“IntegerPreprocess”，“高级”，“HeuristicsMaxNodes”, 1800,“启发式”，“高级”，“BranchRule”，“strongpscost”，“MaxTime”, 7200,“ObjectiveImprovementThreshold”1的军医,“RootLPAlgorithm”，“primal-simplex”）;
             

              [unitCommit,fval,exitflag,output] = solve(unitCommit,unitCommit0，“选项”、选择);
             

            実際に実行すると以下のような結果になりました。

            上図は上からユニット1~6と並んでいます。また各運転状態に対する燃料費の比較図とユニットごとの起動損失の比較図を以下に示します。

            結果から出力が大きく,燃料費の安い石炭火力(ユニット1,2)がベースロードでずっと運転しています。

            ただ,ユニット2に至っては,正午で需要が下がるタイミングで出力を下げて調整を行い,需要が上がる18時を前に最大出力に戻します。また液化天然气で出力が大きく,比較して燃料費の抑えられるユニット6も負荷を調整しながら運転しています。lNGのユニット 3，4と石油のユニット5が残りの分担を行いますが、特に一日の中で起動と停止を行う、いわゆるdss(日間起動停止)運転を行っていることが分かります。

            需要に対する供給力の曲線を表示すると次の通りです。

            赤線で示す需要に対し，各メッシュで積み上げの供給量が同じか，上回っていることが分かります。なみにこのパタンでの運転費の概算は1億6063万4339円です。

            MILPによる並解列の計算で運転するユニットを決められたので,あとは完全に需要=供給力となるように各プラントの出力上下限(P_MAX P_MIN)を考慮しながら燃料費が最小となる出力をQPによって決めます。詳細は割愛しますが，結果のみ掲載します。

            結果から需要に対し，ピッタシ供給できていることが分かります。最終的な燃料費は1億7222万2969円となります。

最後に

            このように発電所の運営は最適化問題として解かれます。

            実際のシステムはもっと複雑で様々な運転制約や燃料在庫問題を考慮し,また作成する計画の長さも1年分・カ1月分・1週間分・1日単位など様々です。普段何気なく利用している電気ですが，一部ではこのように計画し運営されています。

            この冬は電力需給がかなり厳しいと予測されていますが，限りある電気を大事に使っていきたいですね。

余談アプリ化の話

            最後に作ったMATLABのコードは应用设计师を使って共有アプリ化し,MATLAB编译器™を使ってMATLABをもっていない第三者に配布することもできます。

类别:: オペレ，ションズリサ，チ

如欲留言，请点击在这里登录您的MathWorks帐户或创建一个新帐户。

Matlabユザコミュニティ
MATLAB & Simulinkユ，ザ，コミュニティ，向け日本語ブログ

MATLAB & Simulinkユ，ザ，コミュニティ，向け日本語ブログ

発電所はオペレ，ションズ，リサ，チで支えられている

発電機の起動停止計画問題(单位承诺)

milpによる並解列の計算

最後に

余談アプリ化の話

评论

発電機の起動停止計画問題(单位承诺)

milpによる並解列の計算

最後に

余談アプリ化の話

评论

选择网站

美洲

欧洲

亚太地区