業界別最適化Example集: マイクログリッドの最適バッテリーサイジング (エネルギー)
nAG数値計算ライブラリ > 業界別最適化Example集 > マイクログリッドの最適バッテリーサイジング (エネルギー)

Keyword: エネルギー, 最適化, マイクログリッド, バッテリーサイジング, MILP

エネルギー分野の最適化問題:マイクログリッドの最適バッテリーサイジング

問題の概要

microgrid

バッテリーエネルギー貯蔵システム(BESS)は、電力供給の管理、再生可能エネルギー源の信頼性向上、電力系統の安定化において重要な役割を果たしています。効率的なエネルギー貯蔵ソリューションに対する需要が高まるにつれ、洗練された最適化手法の重要性も高まっています。混合整数線形計画法(MILP)は、BESSの運用管理など、一部の変数が整数値である必要がある最適化問題を解くために使用される数学的手法です。

本問題では、外部グリッドに接続され、BESS、複数の発電機、負荷で構成されるマイクログリッドを対象に、運用コストを最小化するバッテリーの最適容量とシステムの最適運用を求めます。BESSの最適化のための数理モデルには、発電機とバッテリーの総運用コストを最小化する目的関数、バッテリーと発電機のスケジュールや仕様、購入電力などの変数、需給バランスや出力制約などの制約条件が含まれます。

マイクログリッドの最適バッテリーサイジングと運用(具体例)

本問題では、外部グリッドに接続された以下のようなマイクログリッドを考えます。需要は発電機、バッテリー、外部グリッドから賄われます。発電機は電力会社が所有しているものとし、エネルギー管理システムは、マイクログリッドを構成するバッテリーの運用制約や需要と供給の一致などを考慮しつつ、発電機とバッテリーの運用コストを最小化するように各構成要素の最適スケジュールを求めます。計画期間は24時間とします。

グリッド接続型マイクログリッドモデル

この問題を解くことで、マイクログリッドに導入すべき最適なバッテリー容量と、24時間の計画期間におけるバッテリーの最適な充放電スケジュールや外部グリッドからの最適な電力供給量が得られます。これにより、発電機の燃料費、バッテリーの導入コスト、外部グリッドからの電力購入費を最小化することができます。

問題の定式化

パラメータ

パラメータ 説明
ngen 発電機の数 3
ai 発電機iの燃料費用関数の1次項係数 a1=0.010694,a2=0.018761,a3=0.007612
bi 発電機iの燃料費用関数の定数項 b1=142.7348,b2=168.9075,b3=313.9102
pi 発電機iの最小出力 [kW] p1=30,p2=50,p3=30
pi 発電機iの最大出力 [kW] p1=70,p2=100,p3=120
UTi 発電機iの最小稼働時間 [h] UTi=8(i=1,2,3)
DTi 発電機iの最小停止時間 [h] DTi=6(i=1,2,3)
cpbat バッテリー出力の単位コスト [$/kW] 0.20
cebat バッテリー容量の単位コスト [$/kWh] 0.25
cjim 時刻jの電力供給単価 [$/kWh] 下図参照
Dj 時刻jの需要 [kW] 下図参照
nhours 計画期間の時間数 24

以下に各時刻における電力単価と需要 を示します。

時刻 電力供給単価 (USD/kWh) 需要 (kW)
0 0.09 200
1 0.08 180
2 0.08 170
3 0.08 160
4 0.08 150
5 0.08 150
6 0.09 170
7 0.11 250
8 0.13 320
9 0.12 300
10 0.11 280
11 0.10 260
12 0.10 270
13 0.10 280
14 0.11 290
15 0.12 300
16 0.13 320
17 0.14 350
18 0.13 340
19 0.12 330
20 0.11 320
21 0.10 280
22 0.09 240
23 0.09 220

グラフにすると、以下のようになります。

電力価格 需要データ

決定変数

変数 説明 範囲
pi,jG 時刻jにおける発電機iの出力 [kW] pipi,jGpi
si,j 時刻jにおける発電機iの起動状態(起動中1、停止中0) si,j{0,1}
pjb 時刻jにおけるバッテリーの充放電電力 [kW] (放電正、充電負) cbpjbcb
cb バッテリー出力の定格 [kW] cb0
eb バッテリー容量の定格 [kWh] eb0
pjim 時刻jにおける外部グリッドからの電力供給 [kW] 0pjim15

目的関数

目的
運用コストの最小化 min j=1nhoursi=1ngen(aipi,jG+bisi,j)+cpbatcb+cebateb+j=1nhourscjimpjim

ユーティリティの総運用コストは、発電機の燃料費、バッテリー出力・容量の設備コスト、電力供給費用の合計で表されます。これを最小化することで、マイクログリッドの経済運用を目指します。

制約条件

制約 説明
需要と供給の一致 i=1ngenpi,jG+pjb+pjim=Dj,j 各時刻において、発電機出力、バッテリー充放電電力、系統からの電力供給の和と需要が一致する
発電機出力上下限 pisi,jpi,jGpisi,j,i,j 発電機の出力は、起動中は上下限の範囲内、停止中は0となる
発電機の最小稼働時間 τ=jmin(j+UTi1,nhours)si,τUTi(si,jsi,j1),i,j2 発電機は、起動すると最小稼働時間を満たす期間は連続して運転する
発電機の最小停止時間 τ=jmin(j+DTi1,nhours)(1si,τ)DTi(si,j1si,j),i,j2 発電機は、停止すると最小停止時間を満たす期間は連続して停止する
バッテリー出力上下限 cbpjbcb,j バッテリーの充放電電力は出力定格の範囲内
バッテリー容量上下限 ebτ=1jpτb0,j バッテリーの充放電量の累積値(充電量)は容量定格の範囲内
系統電力供給上限 0pjim15,j 外部グリッドからの電力輸入は15kW以下に制限

これらの制約条件は、需給バランス、発電機やバッテリーの運用制約、系統連系の制約を表しています。これにより、システムの安定運用とコスト最小化を両立します。

コード例

以下に、この混合整数線形計画問題を nAG Library for Python の MILP ソルバー handle_solve_milp を用いて解くコード例を示します。

今回の問題は、連続変数と離散変数(バイナリ変数)が混在し、目的関数と制約条件がすべて線形であるため、混合整数線形計画法のソルバーである、handle_solve_milpを選択しました。

# 必要なライブラリをインポート
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 発電機のパラメータを定義
n_gen = 3  # 発電機の数
# 発電機の仕様をデータフレームに格納
df = pd.DataFrame(
    {
    'Gen1':[0.010694, 142.7348, 30, 70, 8, 6],
    'Gen2':[0.018761, 168.9075, 50, 100, 8, 6],
    'Gen3':[0.0076121, 313.9102, 30, 120, 8, 6],
    }
)
df.index=['a','b','pl','pu','ut','dt']
print(df.T.to_string())
spec_gen = df.T.to_numpy() 

# バッテリーのパラメータを設定
cost_bat_power = 0.2
cost_bat_energy = 0.25
# 外部グリッドからの電力供給価格
cost_imp_energy = np.array([
    0.09, 0.08, 0.08, 0.08, 0.08, 0.08, 0.09,  # 12:00 AM - 6:00 AM
    0.11, 0.13, 0.12, 0.11, 0.10, 0.10, 0.10,  # 7:00 AM - 1:00 PM
    0.11, 0.12, 0.13, 0.14, 0.13, 0.12, 0.11,  # 2:00 PM - 7:00 PM
    0.10, 0.09, 0.09                         # 8:00 PM - 11:00 PM
])

# 計画期間は24時間
n_hours = 24

# 電力供給コストをプロット
plt.figure(figsize=(7.5, 2.5))
plt.plot(np.arange(n_hours), cost_imp_energy, marker='o')
plt.title('Hourly Electricity Prices')
plt.xlabel('Time (Hour)')
plt.ylabel('Electricity Price (USD/kWh)')
plt.xticks(np.arange(0, 24, step=1))
plt.grid(True)
plt.show()

# 負荷需要を定義
load_demand = np.array([
    200, 180, 170, 160, 150, 150, 170, 250, 320, 300, 280, 260, 270, 280, 290, 
    300, 320, 350, 340, 330, 320, 280, 240, 220
])

# 負荷需要をプロット
plt.figure(figsize=(7.5, 2.5))
plt.plot(np.arange(n_hours), load_demand, marker='o')
plt.title('Load Demand')
plt.xlabel('Time (Hour)')
plt.ylabel('Load Demand (kW)')
plt.xticks(np.arange(0, 24, step=1))
plt.grid(True)
plt.show()

# Python用nAGライブラリをインポート
from naginterfaces.base import utils
from naginterfaces.library import opt, mip

# モデルの変数の総数を計算
nvar = 3*24 + 24  + 1  + 1  + 24 + 3*24 

# モデルデータを保持するための問題ハンドルを作成
handle = opt.handle_init(nvar=nvar)

# 目的関数の係数を設定
idxc = list(range(1, 3*24 + 24  + 1  + 1 + 1))
c = np.tile(spec_gen[:,0], 24)
c = np.concatenate((c, cost_imp_energy, [cost_bat_power, cost_bat_energy]))

# モデルの目的関数を設定
opt.handle_set_quadobj(handle=handle, idxc=idxc, c=c)

# 時間ごとの電力バランス制約を追加
for hour in range(n_hours):
    irowa = np.full(5, 1, dtype=int)
    icola = np.array([3*hour+1, 3*hour+2, 3*hour+3, 3*24+hour+1, 3*24+24+1+1+hour+1], dtype=int)
    a = np.full(5, 1.0, dtype=float)
    bl = np.full(1, load_demand[hour], dtype=float)
    bu = np.full(1, load_demand[hour], dtype=float)
    opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)
    
# 発電機の出力上下限制約を追加
for generator in range(n_gen):
    for hour in range(n_hours):
        irowa = [1, 1, 2, 2]
        icola = np.tile([3*24+24+1+1+24+hour*3+generator+1, hour*3+generator+1],2)
        a = np.array([spec_gen[generator, 2], -1.0, spec_gen[generator, 3], -1.0], dtype=float)
        bl = np.array([-1.e20, 0.])
        bu = np.array([0., 1.e20])
        opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)
        
# 最小稼働時間と最小停止時間制約を追加
ut = 8.0
dt = 6.0
s_start = 3*24 + 24  + 1  + 1  + 24
for generator in range(n_gen):
    for hour in range(1, n_hours):
        irowa = np.array([],dtype=int)
        icola = np.array([],dtype=int)
        a = np.array([],dtype=float)
        for i in range(hour-1, min(hour+ 8, n_hours)):
            icola = np.append(icola, s_start+i*3+generator+1)
            irowa = np.append(irowa, 1)
            a = np.append(a, 1.0)
        a[0] = ut 
        a[1] = 1 - ut
        bl = 0.0
        bu = 1.e20
        opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)
for generator in range(n_gen):
    for hour in range(1, n_hours):
        irowa = np.array([],dtype=int)
        icola = np.array([],dtype=int)
        a = np.array([],dtype=float)
        bu = -1.0
        for i in range(hour-1, min(hour+ 6, n_hours)):
            icola = np.append(icola, s_start+i*3+generator+1)
            irowa = np.append(irowa, 1)
            a = np.append(a, 1.0)
            bu += 1.0
        a[0] = dt 
        a[1] = 1 - dt
        bl = -1.e20
        opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)
        
# バッテリーの出力上下限制約を追加
for hour in range(n_hours):
    irowa = [1, 1, 2, 2]
    icola = [3*24+24+1+1+hour+1, 3*24+24+1, 3*24+24+1+1+hour+1, 3*24+24+1]
    a = [1.0, -1.0, 1.0, 1.0]
    bl = [-1.e20, 0.]
    bu = [0., 1.e20]
    opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)
    
# バッテリーのエネルギー上下限制約を追加
for hour in range(n_hours):
    irowa = np.array([],dtype=int)
    icola = np.array([],dtype=int)
    a = np.array([],dtype=float)
    for i in range(hour+1):
        icola = np.append(icola, 3*24+24+1+1+i+1)
        irowa = np.append(irowa, 1)
        a = np.append(a, 1.0)
    bl = -1.e20
    bu = 0.
    opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)
    irowa = np.array([1],dtype=int)
    icola = np.array([3*24+24+1+1],dtype=int)
    a = np.array([1.],dtype=float)
    for i in range(hour+1):
        icola = np.append(icola, 3*24+24+1+1+i+1)
        irowa = np.append(irowa, 1)
        a = np.append(a, 1.0)
    bl = 0.
    bu = 1.e20
    opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)
    
# 外部グリッドからの電力供給の上限を設定
for i in range(n_hours):
    opt.handle_set_bound(handle=handle, comp='Var', idx=3*24+i+1, bli=0., bui=15.)

# 発電機のスイッチはバイナリ変数である
switches = np.arange(3*24+24+1+1+24+1, nvar+1)
opt.handle_set_property(handle=handle, ptype='Bin', idx=switches)

# オプションを設定
for option in [
        'Print Options = NO',
        'Print Level = 1',
]:
    opt.handle_opt_set(handle, option)

# 略式の反復出力用に明示的なI/Oマネージャを使用
iom = utils.FileObjManager(locus_in_output=False)

# 問題を解く
x, _, _ = mip.handle_solve_milp(handle, io_manager=iom)

# ハンドルを破棄
opt.handle_free(handle)

# 24時間の発電スケジュールを計算
switch1 = np.array([x for i, x in enumerate(x[-3*24:]) if i % 3 ==0])
switch2 = np.array([x for i, x in enumerate(x[-3*24:]) if i % 3 ==1])
switch3 = np.array([x for i, x in enumerate(x[-3*24:]) if i % 3 ==2])
# 24時間の発電量を計算
generator1 = np.array([x for i, x in enumerate(x[:3*24]) if i % 3 ==0])
generator2 = np.array([x for i, x in enumerate(x[:3*24]) if i % 3 ==1])
generator3 = np.array([x for i, x in enumerate(x[:3*24]) if i % 3 ==2])
# 電力供給量を計算
power_import = x[3*24 : 3*24+24]
battery_power_rate = x[3*24+24]
battery_energy_rate = x[3*24+24+1]
battery_power = x[3*24+24+1+1 : 3*24+24+1+1+24]

df = pd.DataFrame(
    {
    'Battery':[battery_power_rate, battery_energy_rate],
    }
)
df.index=['Power Rate','Energy Rate']
print(df.T.to_string())

# バッテリー貯蔵エネルギーのディスパッチをプロット
plt.figure(figsize=(7.5, 2.5))
plt.plot(np.arange(1, n_hours+1), battery_power, marker='o', linestyle='-')
plt.title('Battery Dispatching Power for 24 Hours')
plt.xlabel('Time (hours)')
plt.ylabel('Power (kW)')
plt.xticks(np.arange(0, 25, step=1))
plt.grid(True)
plt.show()

# 外部グリッドからの電力供給をプロット
plt.figure(figsize=(7.5, 2.5))
plt.plot(np.arange(1, n_hours+1), power_import, marker='o', linestyle='-')
plt.title('Import Power for 24 Hours')
plt.xlabel('Time (hours)')
plt.ylabel('Power (kW)')
plt.xticks(np.arange(0, 25, step=1))
plt.grid(True)
plt.show()

結果例

上記のコードを実行すると、以下のような出力結果とプロット図が得られます。

             a         b    pl     pu   ut   dt
Gen1  0.010694  142.7348  30.0   70.0  8.0  6.0
Gen2  0.018761  168.9075  50.0  100.0  8.0  6.0
Gen3  0.007612  313.9102  30.0  120.0  8.0  6.0
 H02BK, Solver for MILP problems
 Status: converged, an optimal solution found
 Final primal objective value  1.218345E+02
 Final dual objective bound    1.218345E+02
         Power Rate  Energy Rate
Battery        45.0        135.0

本問題の最適解として、バッテリーの最適出力率が45 kW、最適エネルギー率が135 kWhと求められました。以下のグラフは、24時間にわたるバッテリーの放電スケジュールと外部グリッドからの電力供給のスケジュールを示しています。これにより、ピーク時に外部からの電力供給を抑えることが可能となり、バッテリーはオフピーク時に充電し、ピーク時に放電して需要を満たしています。供給される電力は15 kW以下に抑えられ、特に第二のピーク期間にのみ発生します。

バッテリー放電スケジュールのプロット図 Battery Dispatching Power for 24 Hours

外部グリッドからの電力供給スケジュールのプロット図 Import Power for 24 Hours

コードの詳細説明

以下のセクションでは、今回のソースコードを上から順番に、部分ごとで何をやっているのかをもう少し詳細に説明します。 Jupyter Notebookを使用している場合は、以下の手順でコードを実行し、結果を確認しながら進めることができます:

  1. コードを部分ごとに新しいセルにコピー&ペーストします。
  2. セルを実行します(Shift + Enterキー、または「Run」ボタンをクリック)。
  3. 実行結果を確認します。
  4. 次の部分のコードに進み、手順を繰り返します。

モデルパラメータと係数

モデルパラメータは、後でモデルを構築するために、以下のコード部分で直接定義されています。各パラメータの説明についてはコメントを参照してください。発電機の燃料コストは通常、二次関数を使用して計算されます。二次関数の特性に基づいて、実際には区分関数を使用して線形化し、強力なMILPソルバーを利用するのが合理的です。ここでは、説明のために二次関数を単純に線形近似して採用します。

# Import necessary libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Generator parameters
# Number of generators
n_gen = 3
# Data frame for the specifications of the generators
df = pd.DataFrame(
    {
    'Gen1':[0.010694, 142.7348, 30, 70, 8, 6],
    'Gen2':[0.018761, 168.9075, 50, 100, 8, 6],
    'Gen3':[0.0076121, 313.9102, 30, 120, 8, 6],
    }
)
df.index=['a','b','pl','pu','ut','dt']
print(df.T.to_string())
spec_gen = df.T.to_numpy() 
             a         b    pl     pu   ut   dt
Gen1  0.010694  142.7348  30.0   70.0  8.0  6.0
Gen2  0.018761  168.9075  50.0  100.0  8.0  6.0
Gen3  0.007612  313.9102  30.0  120.0  8.0  6.0

上記のDataFrameには、発電機の仕様がすべて含まれています。燃料コストは、次の線形関数として定義されます。 costi(pijg)=aipijg+bi ここで、iは発電機、jは時間(時)を表します。plpuは、それぞれi番目の発電機の最小出力と最大出力です。残りのパラメータは次のとおりです。 - ut: 発電機の最小稼働時間 - dt: 発電機の最小停止時間

以下のコード部分を実行すると、時間ごとの電力価格のプロット図が出力されます。

# Battery parameters
cost_bat_power = 0.2
cost_bat_energy = 0.25
# External grid import price
cost_imp_energy = np.array([
    0.09, 0.08, 0.08, 0.08, 0.08, 0.08, 0.09,  # 12:00 AM - 6:00 AM
    0.11, 0.13, 0.12, 0.11, 0.10, 0.10, 0.10,  # 7:00 AM - 1:00 PM
    0.11, 0.12, 0.13, 0.14, 0.13, 0.12, 0.11,  # 2:00 PM - 7:00 PM
    0.10, 0.09, 0.09                           # 8:00 PM - 11:00 PM
])

# Number of hours in the planning horizon
n_hours = 24

# Plot import energy cost
plt.figure(figsize=(7.5, 2.5))
plt.plot(np.arange(n_hours), cost_imp_energy, marker='o')
plt.title('Hourly Electricity Prices')
plt.xlabel('Time (Hour)')
plt.ylabel('Electricity Price (USD/kWh)')
plt.xticks(np.arange(0, 24, step=1))
plt.grid(True)
plt.show()

以下の部分を実行すると負荷需要のプロット図が表示されます。

# Load demand to be fulfilled by the energy management system
load_demand = np.array([
    200, 180, 170, 160, 150, 150, 170, 250, 320, 300, 280, 260, 270, 280, 290, 
    300, 320, 350, 340, 330, 320, 280, 240, 220
])

# Plot load demand
plt.figure(figsize=(7.5, 2.5))
plt.plot(np.arange(n_hours), load_demand, marker='o')
plt.title('Load Demand')
plt.xlabel('Time (Hour)')
plt.ylabel('Load Demand (kW)')
plt.xticks(np.arange(0, 24, step=1))
plt.grid(True)
plt.show()

nAG Libraryを用いて混合整数計画を解く

# Import the nAG Library for Python
from naginterfaces.base import utils
from naginterfaces.library import opt, mip

# Total number of variables in the model
# generators: 3 * 24 * 2
# battery: 24 + 2
# power import: 24
#       pg  + pim + cb + eb + pb +  s   
nvar = 3*24 + 24  + 1  + 1  + 24 + 3*24 

# Create a problem handle to hold model data
handle = opt.handle_init(nvar=nvar)

電力会社の総運転コストである目的関数は、次のように定義されます。 minj=1n_hours(i=1n_gencosti(pijg)+cost_imp_energyjpjim)+cost_bat_powercb+cost_bat_energyeb

# Set objective coefficient
idxc = list(range(1, 3*24

 + 24  + 1  + 1 + 1))
c = np.tile(spec_gen[:,0], 24)
c = np.concatenate((c, cost_imp_energy, [cost_bat_power, cost_bat_energy]))

# Set model objective
opt.handle_set_quadobj(handle=handle, idxc=idxc, c=c)

毎時の負荷バランスを満たす必要があるため、j時の電力収支は次のようになります。 i=1n_genpijg+pjb+pjim=load_demandj.

# Add power balance constraint for each hour
for hour in range(n_hours):
    irowa = np.full(5, 1, dtype=int)
    icola = np.array([3*hour+1, 3*hour+2, 3*hour+3, 3*24+hour+1, 3*24+24+1+1+hour+1], dtype=int)
    a = np.full(5, 1.0, dtype=float)
    bl = np.full(1, load_demand[hour], dtype=float)
    bu = np.full(1, load_demand[hour], dtype=float)
    opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)

発電機の出力は、発電機の仕様でplpuとして与えられた制限内にある必要があります。したがって、i番目の発電機のj時の発電機制限は、次のように定義されます。 plisijpijgpuisij.

# Add generator limits constraint
for generator in range(n_gen):
    for hour in range(n_hours):
        irowa = [1, 1, 2, 2]
        icola = np.tile([3*24+24+1+1+24+hour*3+generator+1, hour*3+generator+1],2)
        a = np.array([spec_gen[generator, 2], -1.0, spec_gen[generator, 3], -1.0], dtype=float)
        bl = np.array([-1.e20, 0.])
        bu = np.array([0., 1.e20])
        opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)

各発電機には、最小稼働時間と最小停止時間の制約があります。 t=jj+ut1situt(sijsij1), t=jj+dt1(1sit)dt(sij1sij). 簡単のため、すべての発電機の最小稼働時間と最小停止時間は同じであると仮定していることに注意してください。

# Add minimum up and down time constraints
ut = 8.0
dt = 6.0
s_start = 3*24 + 24  + 1  + 1  + 24
# minimum up time constraints
for generator in range(n_gen):
    for hour in range(1, n_hours):
        irowa = np.array([],dtype=int)
        icola = np.array([],dtype=int)
        a = np.array([],dtype=float)
        for i in range(hour-1, min(hour+ 8, n_hours)):
            icola = np.append(icola, s_start+i*3+generator+1)
            irowa = np.append(irowa, 1)
            a = np.append(a, 1.0)
        a[0] = ut 
        a[1] = 1 - ut
        bl = 0.0
        bu = 1.e20
        opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)
# minimum down time constraints
for generator in range(n_gen):
    for hour in range(1, n_hours):
        irowa = np.array([],dtype=int)
        icola = np.array([],dtype=int)
        a = np.array([],dtype=float)
        bu = -1.0
        for i in range(hour-1, min(hour+ 6, n_hours)):
            icola = np.append(icola, s_start+i*3+generator+1)
            irowa = np.append(irowa, 1)
            a = np.append(a, 1.0)
            bu += 1.0
        a[0] = dt 
        a[1] = 1 - dt
        bl = -1.e20
        opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)

蓄電池の出力はいかなる時も制限を外れることはできません。そこで、次の出力定格制限を定義します。 cbpjbcb.

# Add battery power limits constraints
for hour in range(n_hours):
    irowa = [1, 1, 2, 2]
    icola = [3*24+24+1+1+hour+1, 3*24+24+1, 3*24+24+1+1+hour+1, 3*24+24+1]
    a = [1.0, -1.0, 1.0, 1.0]
    bl = [-1.e20, 0.]
    bu = [0., 1.e20]
    opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)

蓄電池システムのエネルギー定格は次のようになります。 ebj=1tpjb0,for t=1,,n_hours.

# Add battery energy limits constraints
for hour in range(n_hours):
    # rhs
    irowa = np.array([],dtype=int)
    icola = np.array([],dtype=int)
    a = np.array([],dtype=float)
    for i in range(hour+1):
        icola = np.append(icola, 3*24+24+1+1+i+1)
        irowa = np.append(irowa, 1)
        a = np.append(a, 1.0)
    bl = -1.e20
    bu = 0.
    opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)
    # lhs
    irowa = np.array([1],dtype=int)
    icola = np.array([3*24+24+1+1],dtype=int)
    a = np.array([1.],dtype=float)
    for i in range(hour+1):
        icola = np.append(icola, 3*24+24+1+1+i+1)
        irowa = np.append(irowa, 1)
        a = np.append(a, 1.0)
    bl = 0.
    bu = 1.e20
    opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)

系統からの電力購入は最大15 kWに制限されています。 0pjim15.

# Add limits to the imported power
for i in range(n_hours):
    opt.handle_set_bound(handle=handle, comp='Var', idx=3*24+i+1, bli=0., bui=15.)
# Switches for the generators are binary
switches = np.arange(3*24+24+1+1+24+1, nvar+1)
opt.handle_set_property(handle=handle, ptype='Bin', idx=switches)
# Set options
for option in [
        'Print Options = NO',
        'Print Level = 1',
]:
    opt.handle_opt_set(handle, option)

# Use an explicit I/O manager for abbreviated iteration output:
iom = utils.FileObjManager(locus_in_output=False)

# Solve the problem.
x, _, _ = mip.handle_solve_milp(handle, io_manager=iom)

# Destroy the handle:
opt.handle_free(handle)
 H02BK, Solver for MILP problems
 Status: converged, an optimal solution found
 Final primal objective value  1.218345E+02
 Final dual objective bound    1.218345E+02

これで最適な構成を解から抽出できます。

# Generator schedule for 24 hours
switch1 = np.array([x for i, x in enumerate(x[-3*24:]) if i % 3 ==0])
switch2 = np.array([x for i, x in enumerate(x[-3*24:]) if i % 3 ==1])
switch3 = np.array([x for i, x in enumerate(x[-3*24:]) if i % 3 ==2])
# Generator energy for 24 hours
generator1 = np.array([x for i, x in enumerate(x[:3*24]) if i % 3 ==0])
generator2 = np.array([x for i, x in enumerate(x[:3*24]) if i % 3 ==1])
generator3 = np.array([x for i, x in enumerate(x[:3*24]) if i % 3 ==2])
# Power import
power_import = x[3*24 : 3*24+24]
battery_power_rate = x[3*24+24]
battery_energy_rate = x[3*24+24+1]
battery_power = x[3*24+24+1+1 : 3*24+24+1+1+24]

df = pd.DataFrame(
    {
    'Battery':[battery_power_rate, battery_energy_rate],
    }
)
df.index=['Power Rate','Energy Rate']
print(df.T.to_string())
         Power Rate  Energy Rate
Battery        45.0        135.0

このモデルでは、最適な蓄電池の出力定格は45 kW、エネルギー定格は135 kWhとなります。

蓄電池の放電スケジュールと電力輸入スケジュールを以下のようにプロットできます。出力されるプロット図を見ると、蓄電池はオフピーク時に充電し、ピーク時に放電して外部からのエネルギー輸入の必要性を軽減していることがわかります。輸入電力は15 kW以下に抑えられ、第2ピーク時のみ発生しています。

# Plot dispatching from battery energy storage system
plt.figure(figsize=(7.5, 2.5))
plt.plot(np.arange(1, n_hours+1), battery_power, marker='o', linestyle='-')
plt.title('Battery Dispatching Power for 24 Hours')
plt.xlabel('Time (hours)')
plt.ylabel('Power (kW)')
plt.xticks(np.arange(0, 25, step=1))
plt.grid(True)
plt.show()

次の部分では、24時間の輸入電力をプロットします。

# Plot power import from external grid
plt.figure(figsize=(7.5, 2.5))
plt.plot(np.arange(1, n_hours+1), power_import, marker='o', linestyle='-')
plt.title('Import Power for 24 Hours')
plt.xlabel('Time (hours)')
plt.ylabel('Power (kW)')
plt.xticks(np.arange(0, 25, step=1))
plt.grid(True)
plt.show()

まとめ

本例題では、nAG Library の MILP ソルバーを用いて、マイクログリッドのモデリングとバッテリー容量の最適設計を行いました。また、マイクログリッド内の各構成要素の運用スケジュールも同時に最適化できました。

提案モデルは、太陽光発電や風力発電などの追加電源、その他の制御可能な負荷などを容易に拡張できます。

参考文献


関連情報
© 譌・譛ャ繝九Η繝シ繝。繝ェ繧ォ繝ォ繧「繝ォ繧エ繝ェ繧コ繝�繧コ繧ー繝ォ繝シ繝玲�ェ蠑丈シ夂、セ 2024
Privacy Policy  /  Trademarks