発電機の最適な運転スケジュール: 問題ベース

この例では、収益からコストを引いた後の利益を最大化するよう、2 台のガス火力発電機の最適な運転スケジュールを設定する方法を示します。必ずしも現実的な例ではありませんが、判定のタイミングに応じて変動するコストを考慮する方法が説明されています。

この問題に対するソルバーベースのアプローチについては、発電機の最適な運転スケジュール: ソルバーベースを参照してください。

問題の定義

電力市場の価格は、一日の間にその時間帯によって変動します。発電機を運転する場合、価格が高い時間帯に稼働するようスケジュールを設定することで、この価格の変動を有利に利用することができます。たとえば、2 台の発電機を管理していると仮定します。各発電機にオフ、低、高の 3 つの稼働レベルがあり、そのレベルごとに各発電機の燃料消費と電力生産量が決まっています。発電機がオフのときの燃料消費量は 0 になります。

1 日につき 30 分ごとの間隔 (つまり 24 時間に 48 回) で、各発電機に稼働レベルを割り当てることができます。過去の記録から、この各時間枠におけるメガワット時 (MWh) あたりの収益がわかっているとします。この例のデータは Australian Energy Market Operator 提供による https://www.nemweb.com.au/REPORTS/CURRENT/ 2013 年中期レポートから引用したもので、同社の利用条件 https://www.aemo.com.au/privacy-and-legal-notices/copyright-permissions のもとに使用されています。

load dispatchPrice; % Get poolPrice, which is the revenue per MWh
bar(poolPrice,.5)
xlim([.5,48.5])
xlabel('Price per MWh at each period')

発電機をいったんオフにした後で運転を再開するには、コストがかかります。さらに、1 日の最大燃料消費量の制約も適用されます。この制約があるのは、燃料は前日に購入し、購入しただけの分しか使用できないためです。

問題の表記とパラメーター

このスケジュール設定問題は、0-1 整数計画問題として定式化できます。インデックス i、j および k と、2 値のスケジュール設定ベクトル y を、次のように定義します。

• nPeriods = 期間数 (ここでは 48)。

• i = 1 期間、1 <= i <= 48。

• j = 発電機インデックス、ここでは 1 <= j <= 2。

• 期間 i に発電機 j がレベル k で稼働している場合、y(i,j,k) = 1 となります。ここで低稼働レベルを k = 1 とし、高稼働レベルを k = 2 とします。sum_k y(i,j,k) = 0 の場合、発電機はオフになっています。

オフにした発電機について、運転をいつ再開するかを判断します。そのため、期間 i に発電機 j の運転開始にコストがかかるかどうかを示す補助 2 値変数 z(i,j) を定義します。

• 発電機 j が期間 i にオフで、期間 i + 1 にオンの場合は z(i,j) = 1。それ以外の場合は z(i,j) = 0。言い換えれば、sum_k y(i,j,k) = 0 かつ sum_k y(i+1,j,k) = 1 のときに z(i,j) = 1 となります。

ここで、y の設定に基づいて z を自動的に設定する方法が必要になります。この設定は、次の線形制約によって処理します。

上記のほかに、コスト、各発電機の電力生産レベル、消費レベルおよび使用できる燃料を表すパラメーターも必要です。

• poolPrice(i) -- 期間 i における MWh あたりの収益 (ドル単位)

• gen(j,k) -- 発電機 j によって稼働レベル k で生産される MW 数

• fuel(j,k) -- 発電機 j によって稼働レベル k で消費される燃料

• totalFuel -- 1 日に使用できる燃料

• startCost -- いったんオフにした発電機を再稼働させるためのコスト (ドル単位)

• fuelPrice -- 燃料の単位あたりコスト

poolPrice は load dispatchPrice; の実行時に計算済みです。その他のパラメーターは次のように設定します。

fuelPrice = 3;
totalFuel = 3.95e4;
nPeriods = length(poolPrice); % 48 periods
nGens = 2; % Two generators
gen = [61,152;50,150]; % Generator 1 low = 61 MW, high = 152 MW
fuel = [427,806;325,765]; % Fuel consumption for generator 2 is low = 325, high = 765
startCost = 1e4; % Cost to start a generator after it has been off

発電機の効率性

2 つの稼働点における、2 台の発電機の効率性を調べます。

efficiency = gen./fuel; % Calculate electricity per unit fuel use
rr = efficiency'; % for plotting
h = bar(rr);
h(1).FaceColor = 'g';
h(2).FaceColor = 'c';
legend(h,'Generator 1','Generator 2','Location','NorthEastOutside')
ax = gca;
ax.XTick = [1,2];
ax.XTickLabel = {'Low','High'};
ylim([.1,.2])
ylabel('Efficiency')

ここで、それぞれ対応する稼働点 (低および高) においては発電機 2 の方が発電機 1 よりも若干効率が高くなっていますが、高稼働点における発電機 1 は、低稼働点における発電機 2 よりも効率が良いことに注意してください。

解の変数

問題を設定するには、すべての問題データと制約を、問題形式にエンコードする必要があります。問題の解を表す変数 y(i,j,k) と、発電機の運転開始にコストがかかるかどうかを示す補助変数 z(i,j) を使用します。ここで、y は nPeriods-by-nGens-by-2 の配列、z は nPeriods-by-nGens の配列です。変数はすべて 2 値です。

y = optimvar('y',nPeriods,nGens,{'Low','High'},'Type','integer','LowerBound',0,...
    'UpperBound',1);
z = optimvar('z',nPeriods,nGens,'Type','integer','LowerBound',0,...
    'UpperBound',1);

線形制約

稼働レベルに、1 に等しい成分が複数指定されないように、線形不等式制約を設定します。

powercons = y(:,:,'Low') + y(:,:,'High') <= 1;

1 期間あたりの運転費が、その期間の燃料コストです。発電機 j がレベル k で稼働している場合、コストは fuelPrice * fuel(j,k) となります。

消費されるすべての燃料を考慮した式 fuelUsed を作成します。

yFuel = zeros(nPeriods,nGens,2);
yFuel(:,1,1) = fuel(1,1); % Fuel use of generator 1 in low setting
yFuel(:,1,2) = fuel(1,2); % Fuel use of generator 1 in high setting
yFuel(:,2,1) = fuel(2,1); % Fuel use of generator 2 in low setting
yFuel(:,2,2) = fuel(2,2); % Fuel use of generator 2 in high setting

fuelUsed = sum(sum(sum(y.*yFuel)));

制約は、消費される燃料が使用できる燃料以下であることです。

fuelcons = fuelUsed <= totalFuel;

発電機の稼働開始指標変数の設定

変数 y のオンとオフの期間に応じてソルバーが変数 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,j,k) + y(i+1,j,k) ) > 0 となることがわかります。

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

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

また、目的関数のコストに変数 z を含めます。目的関数内の変数 z に対し、ソルバーはその値を低くしようとします。つまり、ソルバーはすべての値を 0 に等しく設定しようとします。しかし、発電機がオンになっている期間では、線形不等式によって z(i,j) が 1 に等しくなるよう強制されます。

y(i+1,j,k) - y(i,j,k) を表す補助変数 w を作成します。w に関して発電機の運転開始時の不等式を表します。

w = optimexpr(nPeriods,nGens); % Allocate w
idx = 1:(nPeriods-1);
w(idx,:) = y(idx+1,:,'Low') - y(idx,:,'Low') + y(idx+1,:,'High') - y(idx,:,'High');
w(nPeriods,:) = y(1,:,'Low') - y(nPeriods,:,'Low') + y(1,:,'High') - y(nPeriods,:,'High');
switchcons = w - z <= 0;

目的の定義

目的関数には、発電機の運転時の燃料コスト、発電機の運転から得られる収益および発電機を運転開始するためのコストが含まれます。

generatorlevel  = zeros(size(yFuel));
generatorlevel(:,1,1) = gen(1,1); % Fill in the levels
generatorlevel(:,1,2) = gen(1,2);
generatorlevel(:,2,1) = gen(2,1);
generatorlevel(:,2,2) = gen(2,2); 

収益 = y.*generatorlevel.*poolPrice です。

revenue = optimexpr(size(y));
for ii = 1:nPeriods
    revenue(ii,:,:) = poolPrice(ii)*y(ii,:,:).*generatorlevel(ii,:,:);
end

合計燃料コスト = fuelUsed*fuelPrice です。

fuelCost = fuelUsed*fuelPrice;

発電機の運転開始コスト = z*startCost です。

startingCost = z*startCost;

利益 = 収益 - 合計燃料コスト - 運転開始コスト です。

profit = sum(sum(sum(revenue))) - fuelCost - sum(sum(startingCost));

問題を解く

最適化問題を作成し、目的と制約を含めます。

dispatch = optimproblem('ObjectiveSense','maximize');
dispatch.Objective = profit;
dispatch.Constraints.switchcons = switchcons;
dispatch.Constraints.fuelcons = fuelcons;
dispatch.Constraints.powercons = powercons;

スペースを節約するため、反復表示を抑制します。

options = optimoptions('intlinprog','Display','final');

問題を解きます。

[dispatchsol,fval,exitflag,output] = solve(dispatch,'options',options);

Solving problem using intlinprog.

Optimal solution found.

Intlinprog stopped because the objective value is within a gap tolerance of the optimal value, options.RelativeGapTolerance = 0.0001. The intcon variables are integer within tolerance, options.ConstraintTolerance = 1e-06.

解の検証

解を時間の関数としてプロットします。

subplot(3,1,1)
bar(dispatchsol.y(:,1,1)*gen(1,1)+dispatchsol.y(:,1,2)*gen(1,2),.5,'g')
xlim([.5,48.5])
ylabel('MWh')
title('Generator 1 Optimal Schedule','FontWeight','bold')
subplot(3,1,2)
bar(dispatchsol.y(:,2,1)*gen(2,1)+dispatchsol.y(:,2,2)*gen(2,2),.5,'c')
title('Generator 2 Optimal Schedule','FontWeight','bold')
xlim([.5,48.5])
ylabel('MWh')
subplot(3,1,3)
bar(poolPrice,.5)
xlim([.5,48.5])
title('Energy Price','FontWeight','bold')
xlabel('Period')
ylabel('$ / MWh')

発電機 2 は発電機 1 より長い時間稼働します。これは発電機 2 の方が効率性が高いことから容易に予想できます。発電機 2 が稼働するときは、常に最高のレベルで稼働します。発電機 1 は主に高レベルで稼働しますが、1 期間に限り低レベルに落として稼働します。各発電機は毎日連続した期間にわたって稼働するため、運転開始コストが発生するのは毎日 1 度だけです。

発電機が運転を開始する期間の z 変数が 1 であることを確認してください。

starttimes = find(round(dispatchsol.z) == 1); % Use round for noninteger results
[theperiod,thegenerator] = ind2sub(size(dispatchsol.z),starttimes)

theperiod = 2×1

    23
    16

thegenerator = 2×1

     1
     2

発電機がオンになる期間は、プロットに一致しています。

運転開始のペナルティを下げた場合との比較

startCost の値を小さくすると、解で複数の電力生産期間が使用されます。

startCost = 500; % Choose a lower penalty for starting the generators
startingCost = z*startCost;
profit = sum(sum(sum(revenue))) - fuelCost - sum(sum(startingCost));
dispatch.Objective = profit;
[dispatchsolnew,fvalnew,exitflagnew,outputnew] = solve(dispatch,'options',options);

Solving problem using intlinprog.

Optimal solution found.

Intlinprog stopped because the objective value is within a gap tolerance of the optimal value, options.RelativeGapTolerance = 0.0001. The intcon variables are integer within tolerance, options.ConstraintTolerance = 1e-06.

subplot(3,1,1)
bar(dispatchsolnew.y(:,1,1)*gen(1,1)+dispatchsolnew.y(:,1,2)*gen(1,2),.5,'g')
xlim([.5,48.5])
ylabel('MWh')
title('Generator 1 Optimal Schedule','FontWeight','bold')
subplot(3,1,2)
bar(dispatchsolnew.y(:,2,1)*gen(2,1)+dispatchsolnew.y(:,2,2)*gen(2,2),.5,'c')
title('Generator 2 Optimal Schedule','FontWeight','bold')
xlim([.5,48.5])
ylabel('MWh')
subplot(3,1,3)
bar(poolPrice,.5)
xlim([.5,48.5])
title('Energy Price','FontWeight','bold')
xlabel('Period')
ylabel('$ / MWh')

starttimes = find(round(dispatchsolnew.z) == 1); % Use round for noninteger results
[theperiod,thegenerator] = ind2sub(size(dispatchsolnew.z),starttimes)

theperiod = 3×1

    22
    16
    45

thegenerator = 3×1

     1
     2
     2

関連するトピック

• 発電機の最適な運転スケジュール: ソルバーベース
• 整数モデリングと論理モデリング
• 問題ベースの最適化ワークフロー