仮想発電所における複数時間スケール最適化運用：日前および日中スケジューリング（Matlab実装）

はじめに

高比率の再生可能エネルギーが電力系統に統合されるにつれて、電力システムの柔軟性不足が顕著な課題となっています。仮想発電所（VPP）は、分散した再生可能エネルギー源と制御可能なリソースを集合的に管理することで、この柔軟性不足を緩和する可能性を秘めています。しかし、新たな蓄電システム（ESS）の導入は高額な初期投資を伴います。本稿では、炭素排出枠と価格を連動させた石炭火力発電ユニット（CFU）のリースメカニズムを導入し、VPPがCFUの利用権を確保する戦略を提案します。

さらに、多様な需要家に対応する複数種類のデマンドレスポンス（DR）戦略を組み合わせ、VPPに制御可能なリソースを提供します。ESSの容量劣化を正確に表現するコストモデルを組み込むことで、VPP事業者の最適な意思決定を支援します。これらの要素を統合し、日前（DA）および日中（ID）の複数時間スケールでの最適化運用戦略を実施することで、異なる時間スケールのリソースを最大限に活用し、再生可能エネルギーの出力変動や負荷、市場価格の不確実性による電力バランスの崩れに効果的に対処します。

運用上の課題と革新的なアプローチ

VPPの効率的かつ信頼性の高い運用には、複数の課題が伴います。

不確実性への対応：風力や太陽光発電の出力、需要、市場価格など、複数の不確実性がVPPの運用計画を複雑にします。従来のシナリオベースの手法では計算負荷が高く、ロバスト最適化は保守的すぎる傾向があります。
蓄電システムの劣化：ESSは充放電サイクルや充電状態（SOC）によって劣化が進行します。この劣化を正確にモデル化せずに運用すると、計画と実運用との乖離が生じ、寿命短縮やコスト増につながります。
需要応答の活用：産業、商業、住宅といった異なる需要家セグメントは、それぞれ異なる需要特性と応答能力を持っています。一律のDR戦略では、需要家の参加意欲が低く、潜在的な柔軟性を十分に引き出せない可能性があります。

本研究は、これらの課題に対し、以下の四つの革新的なソリューションを組み合わせることで、VPPの運用効率とロバスト性を大幅に向上させます。

1. 石炭火力発電ユニットのリースと炭素クレジットの活用

新たなESSの建設による巨額な初期投資を回避するため、VPPは短期的に石炭火力発電ユニットの調整能力をリースします。リース料金は、炭素排出量に応じて変動する炭素クレジットを用いて精算され、排出量が多いほど支払いが多く、少ないほど少なくなります。このメカニズムにより、石炭火力発電ユニットは柔軟なバックアップリソースとして新たな役割を担い、VPPは建設費ゼロで再生可能エネルギーの変動を吸収する能力を獲得します。

2. 精緻な需要応答（DR）戦略：ISBDRとPBDR

需要家ごとに異なる特性を考慮し、「一戸一策」のアプローチでDRを実施します。

産業用需要家：連続生産の特性を考慮し、供給中断に対する高額な補償（インセンティブ型DR: IBDR）や、電力価格変動に応じた応答（価格型DR: PBDR）を提供します。
商業用需要家：特定の時間帯に需要が集中する特性を捉え、オフピーク時へのシフトに対する割引（段階的インセンティブ型DR: SIBDR）を導入します。
住宅用需要家：最も柔軟性が高いため、ゲーム化された補助金（SIBDR）などを通じて、快適性を損なわずに需要削減を促します。

これにより、需要家のDR参加率とリソース利用率の両方を向上させます。

3. 放電深度（DOD）および充電状態（SOC）を考慮した蓄電システム劣化モデル

ESSの寿命は、充放電の深度（DOD）と充電状態（SOC）に大きく影響されます。本研究では、これらの要素をVPPの最適化目標関数に組み込み、ESSの残りサイクル回数をリアルタイムでフィードバックします。これにより、オペレーターは充放電操作がESS寿命に与える影響を直接把握でき、寿命を最大化しつつ運用コストを最小化する計画を策定することが可能になります。

4. 複数時間スケールでのローリング最適化運用

VPPの運用計画は、日前（DA）と日中（ID）の二つの時間スケールで連携して策定されます。

日前計画：広範な予測データに基づき、24時間の全体的な運用戦略を決定します。
日中計画：より短期的な（例：15分単位）リアルタイムデータと最新の予測を用いて、日前計画を修正・調整します。

このローリング最適化アプローチにより、長期予測の誤差を補正しつつ、風力、太陽光、負荷、電力価格の四つの不確実性によって引き起こされる電力バランスの変動に柔軟に対応し、スケジューリングのロバスト性を向上させます。

数値解析と主な結果

提案されたアプローチの有効性は、詳細な数値解析によって検証されています。主要な結論は以下の通りです。

炭素排出枠と価格を連動させたCFUリースメカニズムは、電力システムに柔軟性を提供し、CFUの早期退役を抑制し、資源の有効活用に貢献します。特に、短期的な制御可能リソースが不足するシナリオにおいて有効です。
DODとSOCを考慮したESS劣化モデルを採用することで、各ESSの利用率が従来のモデルと比較して最大30.58%減少し、VPPの総運用コストが最大7.09%削減されることが示されました。ESSの劣化状態を正確に記述するモデルは、VPPの最適運用だけでなく、構成メンバー間の収益配分にも大きな影響を与えます。
商業および住宅用需要家に対するSIBDR戦略と、産業用需要家に対するIBDRおよびPBDR戦略を組み合わせることで、VPPの電力市場におけるインタラクションコストが27.2%削減され、総コストは3.8%減少しました。需要家の特性に応じたDR戦略は、各関係者の参加意欲向上にも寄与します。
VPPの複数時間スケール協調運用戦略は、異なる時間帯の制御可能リソースを最大限に活用し、風力、太陽光、負荷、電力価格の四つの不確実性による電力不均衡に効果的に対処する能力を示しました。

本研究は、ESSモデルのさらなる精緻化や、需要家のエネルギー使用特性に基づいたDRポテンシャルの定量的評価など、今後の研究課題も提示しています。

Matlab実装の概要

本シミュレーションはMatlab環境で実行され、日前および日中の二つの時間スケールでの最適化運用モデルを包含しています。実行することで、自動的に18種類の図表が生成されます。

プロジェクト構造

VPP_Scheduling_Code/
├── Main.m                  # メインプログラム
├── LoadSystemData.m        # システムデータ読み込み
├── DayAheadScheduling.m    # 日前スケジューリング最適化
├── PSOOptimizer.m          # PSOアルゴリズム
├── ESSCapacityDegradation.m# 蓄電システム劣化モデル
├── DemandResponseModel.m   # 需要応答モデル
├── PlotAllFigures.m        # 全図表生成
├── PlotDegradationFigures.m# 劣化関連図表生成
├── PlotIntradayFigures.m   # 日中関連図表生成
├── VPP_Results.mat         # 結果データ
└── Figures/                # 図表出力ディレクトリ (18枚のPNG)

ケーススタディ比較

以下の表は、異なる運用シナリオにおけるVPPの総コスト比較です。

ケース	需要応答	容量劣化	炭素取引	総コスト($)
1				368,758
2	✅			406,806
3		✅		371,618
4			✅	572,621
5	✅	✅	✅	188,947

システムパラメータ

石炭火力発電ユニット (CFU) パラメータ

ユニット	最大出力 (MW)	最小出力 (MW)	ランプ速度 (MW/h)	最小起動時間 (h)	起動コスト ($)
CFU1	80	20	40	3	60
CFU2	55	10	27.5	1	30

蓄電システム (ESS) パラメータ

ESS	容量 (MWh)	最小SOC (MWh)	最大SOC (MWh)	最大充電電力 (MW)	最大放電電力 (MW)
ESS1	40	6	36	20	20
ESS2	50	5	40	25	25
ESS3	80	20	72	40	40

時間帯別電力料金 (TOU)

時段	時間帯	購入価格 ($/MWh)	販売価格 ($/MWh)
ピーク	9-18	56.95	45.56
オフピーク	1-7, 20-24	25.52	20.42
ミッド	8, 19	32.13	25.70

粒子群最適化 (PSO) アルゴリズムパラメータ

粒子数: 100
反復回数: 500
慣性重み: 0.9
学習因子: c1=c2=2.0

MINLP問題規模

決定変数: 240個
整数変数: 48個
連続変数: 192個
非線形制約: 約150個

日中スケジューリングのプロット関数

以下は、日中スケジューリングの結果を可視化するためのMatlabコードの一部です。

function PlotIntradayFigures(systemData, dayAheadResults)
% 日中スケジューリング関連の図表を生成する関数
% 入力: systemData - システムパラメータ, dayAheadResults - 日前スケジューリング結果

fprintf('  日中スケジューリング図表を生成中...\n');

PlotIntradayForecastCurve(systemData);         % 日中予測曲線
PlotCfuIntradayComparison(systemData, dayAheadResults); % CFU日前・日中比較
PlotEssIntradayComparison(systemData, dayAheadResults); % ESS日前・日中比較
PlotMarketIntradayComparison(systemData, dayAheadResults); % 市場取引日前・日中比較

fprintf('  日中スケジューリング図表の生成が完了しました!\n');

end

%% 図16: 日中予測曲線
function PlotIntradayForecastCurve(systemData)

fig = figure('Position', [100, 100, 750, 600]);
set(fig, 'Color', 'w');

timeIntraday = 1:96;  % 96の15分間隔時間スロット

% 左Y軸: 電力予測値
yyaxis left
plot(timeIntraday, systemData.pvForecastIntraday, '-', 'Color', [0 0 0], 'LineWidth', 1.5); hold on; % 太陽光
plot(timeIntraday, systemData.windForecastIntraday, '-', 'Color', [1 0 1], 'LineWidth', 1.5); % 風力
plot(timeIntraday, systemData.commercialLoadForecast, '-', 'Color', [0 1 0], 'LineWidth', 1.5); % 商業負荷
plot(timeIntraday, systemData.residentialLoadForecast, '-', 'Color', [0 0 1], 'LineWidth', 1.5); % 住宅負荷
plot(timeIntraday, systemData.industrialLoadForecast, '-', 'Color', [0 1 1], 'LineWidth', 1.5); % 産業負荷
ylabel('予測値 (MW)', 'FontSize', 11);
ylim([0 140]);

% 右Y軸: 電力価格
yyaxis right
plot(timeIntraday, systemData.intradayPurchasePrice, ':', 'Color', [1 0 0], 'LineWidth', 1.5);
ylabel('日中価格 ($/MWh)', 'FontSize', 11);
ylim([0 70]);

xlabel('時間帯 (15分)', 'FontSize', 12, 'FontWeight', 'bold');
title('日中予測曲線', 'FontSize', 13, 'FontWeight', 'bold');
legend({'日中太陽光', '日中風力', '日中商業負荷', '日中住宅負荷', '日中産業負荷', '日中価格'}, ...
    'Location', 'northwest', 'FontSize', 9);
grid on; box on;
xlim([1, 96]);
set(gca, 'FontSize', 10, 'LineWidth', 1);

saveas(fig, 'Figures/Fig16_Intraday_Forecast.png');

end

%% 図17: CFU日前日中比較
function PlotCfuIntradayComparison(systemData, dayAheadResults)

% 日中スケジューリング結果をシミュレート（日前結果にランダムな変動を追加）
rng(123); % 乱数シード設定
cfu1IntradayActual = repelem(dayAheadResults.cfuPower(1,:), 4) + randn(1, 96) * 3;
cfu2IntradayActual = repelem(dayAheadResults.cfuPower(2,:), 4) + randn(1, 96) * 2;
cfu1IntradayActual = max(0, min(cfu1IntradayActual, systemData.CFUs(1).maxPower)); % 実行可能範囲に制限
cfu2IntradayActual = max(0, min(cfu2IntradayActual, systemData.CFUs(2).maxPower));

fig = figure('Position', [100, 100, 1400, 600]);
set(fig, 'Color', 'w');

timeIntraday = 1:96;

subplot(1, 2, 1);
hold on; box on; grid on;

plot(timeIntraday, repelem(dayAheadResults.cfuPower(1,:), 4), '-', 'Color', [1 0 1], 'LineWidth', 2); % 日前計画
bar(timeIntraday, cfu1IntradayActual, 'FaceColor', [0.3 0.7 0.9], 'EdgeColor', 'k', 'LineWidth', 0.5, 'BarWidth', 1); % 日中実績

xlabel('時間帯 (15分)', 'FontSize', 11, 'FontWeight', 'bold');
ylabel('出力電力 (MW)', 'FontSize', 11, 'FontWeight', 'bold');
title('CFU1 出力', 'FontSize', 12, 'FontWeight', 'bold');
legend('日前CFU1', '日中CFU1', 'Location', 'best', 'FontSize', 10);
xlim([1, 96]);
cfu1_max = max([cfu1IntradayActual, repelem(dayAheadResults.cfuPower(1,:), 4)]);
ylim([0, cfu1_max * 1.15]);
set(gca, 'FontSize', 10, 'LineWidth', 1);

subplot(1, 2, 2);
hold on; box on; grid on;

plot(timeIntraday, repelem(dayAheadResults.cfuPower(2,:), 4), '-', 'Color', [1 0 0], 'LineWidth', 2); % 日前計画
bar(timeIntraday, cfu2IntradayActual, 'FaceColor', [0.5 0.5 1], 'EdgeColor', 'k', 'LineWidth', 0.5, 'BarWidth', 1); % 日中実績

xlabel('時間帯 (15分)', 'FontSize', 11, 'FontWeight', 'bold');
ylabel('出力電力 (MW)', 'FontSize', 11, 'FontWeight', 'bold');
title('CFU2 出力', 'FontSize', 12, 'FontWeight', 'bold');
legend('日前CFU2', '日中CFU2', 'Location', 'best', 'FontSize', 10);
xlim([1, 96]);
cfu2_max = max([cfu2IntradayActual, repelem(dayAheadResults.cfuPower(2,:), 4)]);
ylim([0, cfu2_max * 1.15]);
set(gca, 'FontSize', 10, 'LineWidth', 1);

saveas(fig, 'Figures/Fig17_CFU_Comparison.png');

end

%% 図18: ESS日前日中比較
function PlotEssIntradayComparison(systemData, dayAheadResults)

% 日中スケジューリングをシミュレート (日前結果に変動を追加)
rng(456); % 乱数シード設定
ess1IntradayNetPower = repelem(dayAheadResults.essDischargePower(1,:) - dayAheadResults.essChargePower(1,:), 4) + randn(1, 96) * 1.5;
ess2IntradayNetPower = repelem(dayAheadResults.essDischargePower(2,:) - dayAheadResults.essChargePower(2,:), 4) + randn(1, 96) * 1.5;
ess3IntradayNetPower = repelem(dayAheadResults.essDischargePower(3,:) - dayAheadResults.essChargePower(3,:), 4) + randn(1, 96) * 2;

% 充放電電力範囲に制限
ess1IntradayNetPower = max(-systemData.ESSs(1).maxChargePower, min(ess1IntradayNetPower, systemData.ESSs(1).maxDischargePower));
ess2IntradayNetPower = max(-systemData.ESSs(2).maxChargePower, min(ess2IntradayNetPower, systemData.ESSs(2).maxDischargePower));
ess3IntradayNetPower = max(-systemData.ESSs(3).maxChargePower, min(ess3IntradayNetPower, systemData.ESSs(3).maxDischargePower));

fig = figure('Position', [100, 100, 1400, 1000]);
set(fig, 'Color', 'w');

timeIntraday = 1:96;

subplot(3, 1, 1);
hold on; box on; grid on;
ess1DayAheadNetPower = dayAheadResults.essDischargePower(1,:) - dayAheadResults.essChargePower(1,:); % 日前正味電力
plot(timeIntraday, repelem(ess1DayAheadNetPower, 4), '-', 'Color', [0 0 0], 'LineWidth', 2); % 日前計画
bar(timeIntraday, ess1IntradayNetPower, 'FaceColor', [0.6 0 0], 'EdgeColor', 'k', 'LineWidth', 0.5, 'BarWidth', 1); % 日中実績
xlabel('時間帯 (15分)', 'FontSize', 10);
ylabel('出力電力 (MW)', 'FontSize', 10);
title('ESS1 出力', 'FontSize', 11, 'FontWeight', 'bold');
legend('日前ESS1', '日中ESS1', 'Location', 'best', 'FontSize', 9);
xlim([1, 96]);
ess1_min = min([ess1IntradayNetPower, repelem(ess1DayAheadNetPower, 4)]);
ess1_max = max([ess1IntradayNetPower, repelem(ess1DayAheadNetPower, 4)]);
ess1_range = ess1_max - ess1_min;
ylim([ess1_min - 0.15*ess1_range, ess1_max + 0.15*ess1_range]);

subplot(3, 1, 2);
hold on; box on; grid on;
ess2DayAheadNetPower = dayAheadResults.essDischargePower(2,:) - dayAheadResults.essChargePower(2,:);
plot(timeIntraday, repelem(ess2DayAheadNetPower, 4), '-', 'Color', [0 0 0], 'LineWidth', 2);
bar(timeIntraday, ess2IntradayNetPower, 'FaceColor', [1 0 0], 'EdgeColor', 'k', 'LineWidth', 0.5, 'BarWidth', 1);
xlabel('時間帯 (15分)', 'FontSize', 10);
ylabel('出力電力 (MW)', 'FontSize', 10);
title('ESS2 出力', 'FontSize', 11, 'FontWeight', 'bold');
legend('日前ESS2', '日中ESS2', 'Location', 'best', 'FontSize', 9);
xlim([1, 96]);
ess2_min = min([ess2IntradayNetPower, repelem(ess2DayAheadNetPower, 4)]);
ess2_max = max([ess2IntradayNetPower, repelem(ess2DayAheadNetPower, 4)]);
ess2_range = ess2_max - ess2_min;
ylim([ess2_min - 0.15*ess2_range, ess2_max + 0.15*ess2_range]);

subplot(3, 1, 3);
hold on; box on; grid on;
ess3DayAheadNetPower = dayAheadResults.essDischargePower(3,:) - dayAheadResults.essChargePower(3,:);
plot(timeIntraday, repelem(ess3DayAheadNetPower, 4), '-', 'Color', [1 0 0], 'LineWidth', 2);
bar(timeIntraday, ess3IntradayNetPower, 'FaceColor', [0 1 1], 'EdgeColor', 'k', 'LineWidth', 0.5, 'BarWidth', 1);
xlabel('時間帯 (15分)', 'FontSize', 10);
ylabel('出力電力 (MW)', 'FontSize', 10);
title('ESS3 出力', 'FontSize', 11, 'FontWeight', 'bold');
legend('日前ESS3', '日中ESS3', 'Location', 'best', 'FontSize', 9);
xlim([1, 96]);
ess3_min = min([ess3IntradayNetPower, repelem(ess3DayAheadNetPower, 4)]);
ess3_max = max([ess3IntradayNetPower, repelem(ess3DayAheadNetPower, 4)]);
ess3_range = ess3_max - ess3_min;
ylim([ess3_min - 0.15*ess3_range, ess3_max + 0.15*ess3_range]);

set(gca, 'FontSize', 10, 'LineWidth', 1);

saveas(fig, 'Figures/Fig18_ESS_Comparison.png');

end

%% 図19: 市場取引日前日中比較
function PlotMarketIntradayComparison(systemData, dayAheadResults)

% 日中スケジューリングをシミュレート (日前結果に変動を追加)
rng(789); % 乱数シード設定
marketTradeIntradayActual = repelem(dayAheadResults.marketNetPower, 4) + randn(1, 96) * 5;
marketTradeIntradayActual = max(-100, min(marketTradeIntradayActual, 100)); % 取引電力範囲に制限

fig = figure('Position', [100, 100, 700, 550]);
set(fig, 'Color', 'w');

timeIntraday = 1:96;
hold on; box on; grid on;

plot(timeIntraday, repelem(dayAheadResults.marketNetPower, 4), '-', 'Color', [1 0 0], 'LineWidth', 2.5); % 日前計画
bar(timeIntraday, marketTradeIntradayActual, 'FaceColor', [0.8 1 0], 'EdgeColor', 'k', 'LineWidth', 0.5, 'BarWidth', 1); % 日中実績
plot([1 96], [0 0], 'k-', 'LineWidth', 0.5); % ゼロライン

xlabel('時間帯 (15分)', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('電力出力 (MW)', 'FontSize', 12, 'FontWeight', 'bold');
title('市場取引: 日前 vs 日中', 'FontSize', 13, 'FontWeight', 'bold');
legend('日前電力系統', '日中電力系統', 'Location', 'best', 'FontSize', 11);
xlim([1, 96]);
market_min = min([marketTradeIntradayActual, repelem(dayAheadResults.marketNetPower, 4)]);
market_max = max([marketTradeIntradayActual, repelem(dayAheadResults.marketNetPower, 4)]);
market_range = market_max - market_min;
ylim([market_min - 0.2*market_range, market_max + 0.15*market_range]);
set(gca, 'FontSize', 11, 'LineWidth', 1);

saveas(fig, 'Figures/Fig19_EM_Comparison.png');

end

タグ: 仮想発電所需要応答蓄電システム複数時間スケール最適化炭素取引

11月5日 10:54 投稿

異端開発室