この記事は「【リアルタイム電力予測】需要・価格・電源最適化ダッシュボード構築記」シリーズの十六日目です
Day1はこちら | 全体構成を見る

前々回は変数・目的関数・制約を定義、前回は上限やコストといった定数をどのように決めるかについて記載していきました。今日は全体のモデルの完成を目指していきます。

以下のように内容を分類していきます。結果分析と可視化は明日行っていきます。

数理最適化モデルについて
変数定義
目的関数（コスト最小化）
制約の追加（需給、上限、ランプ、エネルギー収支）
データ読み込み（timeseries.parquet, params_resolved.yaml）
求解
結果を DataFrame に落とす

1. 数理最適化モデルについて

今回のモデルは「経済負荷配分（Economic Dispatch）」で、基本的には連続変数のみで構成されるため、問題は 線形計画問題（LP）になります。

ライブラリ

問題を「数式モデル」として記述するためのインターフェース
例:
PuLP: 線形計画問題を Python で書けるようにする
Pyomo: 線形・非線形・整数計画など幅広く記述可能
cvxpy: 凸最適化を宣言的に書ける

ソルバー

与えられた数理モデルを実際に解くアルゴリズム
例：
CBC（オープンソースのMILPソルバー）
GLPK（線形計画ソルバー）
Gurobi, CPLEX（商用ソルバー、非常に高速）
OSQP（凸二次計画ソルバー）

今回 Pyomo でモデルを作り、HiGHSソルバー で解いていきます。

2. 変数定義

2.1 時間インデックス

変数の定義には決まった型があるので順番に設定していきます。
まず、モデルオブジェクトを作り、インデックス集合（Set）を定義します。

df = df_ts.copy().reset_index(drop=True)
T = list(df.index) # 30分刻みの時間集合

m = ConcreteModel()
m.T = Set(initialize=T)

その集合を使って変数や制約を「集合に属する要素ごと」に定義していきます。
m.Tは時間集合を表すインデックスで、後でVar(m.T)のようにＴの各要素に対応する変数を書きます。

2.1 発電（供給側）の変数

電源別に発電量（MW）を定義しています（すべて非負）。

太陽光：$P^{pv}_t$
風力：$P^{wind}_t$
水力：$P^{hydro}_t$
火力：$P^{lng}_t, P^{coal}_t, P^{oil}_t, P^{th_other}_t$
バイオマス：$P^{biomass}_t$
その他：$P^{misc}_t$

変数の領域にNonNegativeRealsを設定することで0以上の連続値しか取れないように制約します。

Reals → 制約なしの実数
NonNegativeIntegers → 0以上の整数
Binary → 0か1のみ
PositiveReals → 0より大きい実数

m.P_pv   = Var(m.T, within=NonNegativeReals)
m.P_wind = Var(m.T, within=NonNegativeReals)
m.P_hy   = Var(m.T, within=NonNegativeReals)
m.P_lng  = Var(m.T, within=NonNegativeReals)
m.P_coal = Var(m.T, within=NonNegativeReals)
m.P_oil  = Var(m.T, within=NonNegativeReals)
m.P_tho  = Var(m.T, within=NonNegativeReals)
m.P_biom = Var(m.T, within=NonNegativeReals)
m.P_misc = Var(m.T, within=NonNegativeReals)

2.2 出力抑制

「潜在的には出せるが、系統都合で捨てる分」を変数で持ちます。

太陽光出力抑制：$C^{pv}_t$
風力出力抑制：$C^{wind}_t$

m.C_pv   = Var(m.T, within=NonNegativeReals)
m.C_wind = Var(m.T, within=NonNegativeReals)

2.3 予備力

予備力制約を明示的に扱うために「予備力（余力）」も変数として定義します。

m.R_hy   = Var(m.T, domain=NonNegativeReals)
m.R_coal = Var(m.T, domain=NonNegativeReals)
m.R_oil  = Var(m.T, domain=NonNegativeReals)
m.R_lng  = Var(m.T, domain=NonNegativeReals)
m.R_th_other  = Var(m.T, domain=NonNegativeReals)
m.R_biomass = Var(m.T, domain=NonNegativeReals)
m.R_misc   = Var(m.T, domain=NonNegativeReals)

2.4 未供給

m.Shed   = Var(m.T, domain=NonNegativeReals)  # 未供給（停電）

2.5 蓄電池・揚水・連系線

# 蓄電池
m.P_ch   = Var(m.T, domain=NonNegativeReals)  # 充電
m.P_dis  = Var(m.T, domain=NonNegativeReals)  # 放電
m.E_bat  = Var(m.T, domain=NonNegativeReals)  # SoC
m.R_bat  = Var(m.T, domain=NonNegativeReals)  # 上方予備

# 揚水
m.P_pump = Var(m.T, domain=NonNegativeReals)  # 揚水（消費）
m.P_gen  = Var(m.T, domain=NonNegativeReals)  # 発電
m.E_ps   = Var(m.T, domain=NonNegativeReals)
m.R_ps   = Var(m.T, domain=NonNegativeReals)

# 連系線（インポート）
m.P_imp  = Var(m.T, domain=NonNegativeReals)  # 流入
m.R_imp  = Var(m.T, domain=NonNegativeReals)  # 予備力としての余力

3. 目的関数（コスト最小化）

目的は「1日の総コスト最小化」です。
最適化で扱う発電量は MW（電力）ですが、燃料単価は円/MWh（電力量あたり）なので、30分スロットを MWh に直すために$Δt$を掛けます。

def obj_rule(_):
    return sum(
        dt * (
            costs["coal"] * m.P_coal[t]
            + costs["oil"] * m.P_oil[t]
            + costs["lng"] * m.P_lng[t]
            + costs["th_other"] * m.P_th_other[t]
            + costs["biomass"] * m.P_biomass[t]
            + costs["misc"] * m.P_misc[t]
            + c_curt_pv * m.Curt_pv[t]
            + c_curt_wind * m.Curt_wind[t]
            + c_import * m.P_imp[t]
            + c_shed * m.Shed[t]
        )
        for t in T
    )
m.OBJ = Objective(rule=obj_rule, sense=minimize)

再エネ抑制・未供給（停電）に大きなペナルティを置くことで、
「できる限り供給する」「余った再エネはなるべく捨てない」方向に誘導します
連系線の輸入にもコストを置くことで、無制限に頼らないようにします

4. 制約の追加（需給・上限・ランプ・エネルギー収支）

「物理的にあり得ない解を排除する」「運用上のルールを反映する」役割を持ちます。

4.1 再エネの可用性（上限＋抑制）

太陽光・風力は「使える上限（avail）」が与えられ、その範囲で「発電するか／捨てるか（抑制）」を選びます。
Constraint は「制約を一つずつ名前付きで定義する方法です。

m.PV_av = Param(m.T, initialize=pv_av, within=Reals)
m.W_av  = Param(m.T, initialize=w_av, within=Reals)

m.PVAvail   = Constraint(m.T, rule=lambda _, t: m.P_pv[t]   + m.Curt_pv[t]   == m.PV_av[t])
m.WindAvail = Constraint(m.T, rule=lambda _, t: m.P_wind[t] + m.Curt_wind[t] == m.W_av[t])

4.2 出力上限制約（「出力」＋「予備力」が設備容量を超えない）

# 出力上限（P_max - 予備力）
m.CapHy = Constraint(m.T, rule=lambda _, t: m.P_hy[t] <= Pmax_hy - m.R_hy[t])
m.CapCoal = Constraint(m.T, rule=lambda _, t: m.P_coal[t] <= Pmax_th["coal"] - m.R_coal[t])
m.CapOil = Constraint(m.T, rule=lambda _, t: m.P_oil[t] <= Pmax_th["oil"] - m.R_oil[t])
m.CapLng = Constraint(m.T, rule=lambda _, t: m.P_lng[t] <= Pmax_th["lng"] - m.R_lng[t])
m.CapOth = Constraint(m.T, rule=lambda _, t: m.P_th_other[t] <= Pmax_th["th_other"] - m.R_th_other[t])
m.CapBio = Constraint(m.T, rule=lambda _, t: m.P_biomass[t] <= Pmax_th["biomass"] - m.R_biomass[t])
m.CapMis = Constraint(m.T, rule=lambda _, t: m.P_misc[t] <= Pmax_th["misc"] - m.R_misc[t])

4.3 連系線（輸入）の容量制約

「使っている分」と「予備力として残している分」の合計が上限を超えないようにします。

m.CapTie = Constraint(m.T, rule=lambda _, t: m.P_imp[t] + m.R_imp[t] <= Cap_tie)

4.4 予備力制約（最低限必要な余力を確保）

予備率$𝜌$と需要$𝐿_t$に対して、上方予備力の合計が必要量以上であることを要求します。

def reserve_rule(_, t):
    return (
        m.R_hy[t]
        + m.R_coal[t] + m.R_oil[t] + m.R_lng[t]
        + m.R_th_other[t] + m.R_biomass[t] + m.R_misc[t]
        + m.R_bat[t] + m.R_ps[t] + m.R_imp[t]
    ) >= rho * load[t]
m.Reserve = Constraint(m.T, rule=reserve_rule)

4.5 蓄電池のエネルギー収支

蓄電池は「時間をまたぐ」ので、残量（MWh）を持ちます。

ConstraintList を使って、t=0 と t>=1 を分けて追加しています。
ConstraintList は「制約をリストに入れて、必要なだけ追加できる方法です。

# 蓄電池のエネルギーバランス
eta_ch = float(bat_conf["eta_ch"])
eta_dis = float(bat_conf["eta_dis"])
Emax_b = float(bat_conf["E_max"])
Pchmax = float(bat_conf["P_ch_max"])
Pdismax = float(bat_conf["P_dis_max"])
E0_b = float(bat_conf["E0"])

m.BatSOC = ConstraintList()
m.BatCap = ConstraintList()

# t=0
m.BatSOC.add(m.E_bat[0] == E0_b + eta_ch*m.P_ch[0]*dt - (1/eta_dis)*m.P_dis[0]*dt)

# t>=1
for t in T[1:]:
    m.BatSOC.add(m.E_bat[t] == m.E_bat[t-1] + eta_ch*m.P_ch[t]*dt - (1/eta_dis)*m.P_dis[t]*dt)
    m.BatCap.add(m.E_bat[t] <= Emax_b)
    m.BatCap.add(m.P_ch[t]  <= Pchmax)
    m.BatCap.add(m.P_dis[t] <= Pdismax)

容量・出力の上限：

m.BatCap.add(m.E_bat[0] <= Emax_b)
m.BatCap.add(m.P_ch[0]  <= Pchmax)
m.BatCap.add(m.P_dis[0] <= Pdismax)

蓄電池の予備力（上方予備）:
「これから追加で放電できる余力」を予備力とみなします。

    # バッテリー予備力（上方：これから放電できる余力）
    m.BatRes = ConstraintList()
    for t in T:
        m.BatRes.add(m.R_bat[t] <= Pdismax - m.P_dis[t])
        m.BatRes.add(m.R_bat[t] <= m.E_bat[t] / dt)

4.6 揚水のエネルギー収支（上池エネルギー）

揚水も同様に「時間をまたぐ」ので、上池相当のエネルギー残量を持ちます。

# 揚水のエネルギーバランス
eta_pump = float(ps_conf["eta_pump"])
eta_gen = float(ps_conf["eta_gen"])
Emax_ps = float(ps_conf["E_max"])
Ppumpmax = float(ps_conf["P_pump_max"])
Pgenmax = float(ps_conf["P_gen_max"])
E0_ps = float(ps_conf["E0"])

m.PSSOC = ConstraintList()
m.PSCap = ConstraintList()
    
m.PSSOC.add(m.E_ps[0] == E0_ps + eta_pump*m.P_pump[0]*dt - (1/eta_gen)*m.P_gen[0]*dt)

for t in T[1:]:
    m.PSSOC.add(m.E_ps[t] == m.E_ps[t-1] + eta_pump*m.P_pump[t]*dt - (1/eta_gen)*m.P_gen[t]*dt)
    m.PSCap.add(m.E_ps[t] <= Emax_ps)
    m.PSCap.add(m.P_pump[t] <= Ppumpmax)
    m.PSCap.add(m.P_gen[t]  <= Pgenmax)

容量制約：

m.PSCap.add(m.E_ps[0] <= Emax_ps)
m.PSCap.add(m.P_pump[0] <= Ppumpmax)
m.PSCap.add(m.P_gen[0]  <= Pgenmax)

予備力

# 揚水予備力
m.PSRes = ConstraintList()
for t in T:
    m.PSRes.add(m.R_ps[t] <= Pgenmax - m.P_gen[t])
    m.PSRes.add(m.R_ps[t] <= m.E_ps[t] / dt)

4.7 ランプ制約（火力・水力）

急に出力を変えられない制約です。

# ランプ制約（火力と水力）
init = init_state or {k: 0.0 for k in ["hy"] + thermal_keys}

def _get_init(key: str) -> float:
    return float(init.get(key, 0.0))

m.Ramp = ConstraintList()
# t=0
m.Ramp.add(m.P_hy[0] - _get_init("hy") <= RU_hy)
m.Ramp.add(_get_init("hy") - m.P_hy[0] <= RD_hy)
for g in thermal_keys:
    Pg = getattr(m, f"P_{g}")
    m.Ramp.add(Pg[0] - _get_init(g) <= RU_th[g])
    m.Ramp.add(_get_init(g) - Pg[0] <= RD_th[g])

# t>=1
for t in T[1:]:
    m.Ramp.add(m.P_hy[t] - m.P_hy[t - 1] <= RU_hy)
    m.Ramp.add(m.P_hy[t - 1] - m.P_hy[t] <= RD_hy)
    for g in thermal_keys:
        Pg = getattr(m, f"P_{g}")
        m.Ramp.add(Pg[t] - Pg[t - 1] <= RU_th[g])
        m.Ramp.add(Pg[t - 1] - Pg[t] <= RD_th[g])

4.8 水力の日次エネルギー上限（1日の総発電量）

水力を「燃料費0」としていると出し放題になりやすいので、1日あたりの上限（MWh）を置きます。

m.HydroDay = ConstraintList()
df["date"] = df["timestamp"].dt.date
for d, idx in df.groupby("date").groups.items():
    m.HydroDay.add(sum(m.P_hy[t]*dt for t in idx) <= E_day_hy)

4.9 需給バランス（同時同量）

「供給（発電＋放電＋揚水発電＋輸入）」が「需要＋充電＋揚水ポンプ＋未供給」に一致します。

# 需給バランス
# 左辺：各電源の実出力 + インポート + 揚水発電 + 蓄電池放電
# 右辺：需要 + 揚水の揚水 + 蓄電池充電 + 未供給
def balance_rule(_, t):
    return (
        m.P_pv[t] + m.P_wind[t] + m.P_hy[t]
        + m.P_coal[t] + m.P_oil[t] + m.P_lng[t]
        + m.P_th_other[t] + m.P_biomass[t] + m.P_misc[t]
        + m.P_gen[t] + m.P_dis[t] + m.P_imp[t]
    ) == (
        load[t] + m.P_pump[t] + m.P_ch[t] + m.Shed[t]
    )

m.Balance = Constraint(m.T, rule=balance_rule)

5. データ読み込み

df_ts: timestamp, predicted_demand, pv_avail_MW, wind_avail_MW など
params: params_resolved.yaml を読み込んで「今日の上限が埋まった状態」
costs: fuel_cost.parquet 等から「今日の円/MWh」を渡す
init_state: 前日23:30の各電源の発電量

上記のdataframeやyamlを以下のように読み込んで使います。

dt = float(params.get("dt_hours", 0.5))
rho = float(params["reserve"]["rho"])
c_curt_pv = float(params["penalty"]["curtail_pv"])
c_curt_wind = float(params["penalty"]["curtail_wind"])
c_shed = float(params["penalty"]["shed"])

# 時間インデックス
df = df_ts.copy().reset_index(drop=True)
T = list(df.index)

m = ConcreteModel()
m.T = Set(initialize=T)

# パラメータ（負荷と再エネ上限）
load  = {t: float(df.loc[t, "predicted_demand"])       for t in T}
pv_av = {t: float(df.loc[t, "pv_avail_MW"])   for t in T}
w_av  = {t: float(df.loc[t, "wind_avail_MW"]) for t in T}

# P_max & ramp を params から取得
Pmax_hy = float(params["hydro"]["P_max"])
RU_hy = float(params["hydro"]["RU"])
RD_hy = float(params["hydro"]["RD"])
E_day_hy = float(params["hydro"]["E_day"])

thermal_keys = ["coal", "oil", "lng", "th_other", "biomass", "misc"]
Pmax_th = {g: float(params["thermal"][g]["P_max"]) for g in thermal_keys}
RU_th = {g: float(params["thermal"][g]["RU"]) for g in thermal_keys}
RD_th = {g: float(params["thermal"][g]["RD"]) for g in thermal_keys}

Cap_tie = float(params["tie"]["cap"])
c_import = float(params["tie"]["c_import"])

bat_conf = params["battery"]
ps_conf = params["pstorage"]

6. 求解（ソルバー実行）

モデルを作ったら、HiGHS で解きます。

solver = SolverFactory("highs")
solver.solve(m, tee=False)

tee=True にするとソルバーのログが出るので、デバッグ時には便利です。

7. 結果を DataFrame にする

Pyomo 変数はそのままでは pandas にできないので、value() で数値に戻して DataFrame にします。
可視化・検算用として「total_cost（円）」も列として作っています。

# 結果DataFrame
out = pd.DataFrame({
    "timestamp": df["timestamp"],
    "pv":    [value(m.P_pv[t])   for t in T],
    "wind":  [value(m.P_wind[t]) for t in T],
    "hydro": [value(m.P_hy[t])   for t in T],
    "coal":  [value(m.P_coal[t]) for t in T],
    "oil":   [value(m.P_oil[t])  for t in T],
    "lng":   [value(m.P_lng[t])  for t in T],
    "th_other":   [value(m.P_th_other[t])  for t in T],
    "biomass": [value(m.P_biomass[t]) for t in T],
    "misc":   [value(m.P_misc[t])   for t in T],
    "pstorage_gen": [value(m.P_gen[t])  for t in T],
    "pstorage_pump":[value(m.P_pump[t]) for t in T],
    "pstorage_soc":  [value(m.E_ps[t])  for t in T],
    "battery_dis":  [value(m.P_dis[t])  for t in T],
    "battery_ch":   [value(m.P_ch[t])   for t in T],
    "battery_ch":   [value(m.P_ch[t])   for t in T],
    "import":       [value(m.P_imp[t])  for t in T],
    "curtail_pv":   [value(m.Curt_pv[t])   for t in T],
    "curtail_wind": [value(m.Curt_wind[t]) for t in T],
    "reserve": [
        value(m.R_hy[t] + m.R_coal[t] + m.R_oil[t] + m.R_lng[t]
              + m.R_bat[t] + m.R_ps[t] + m.R_imp[t])
        for t in T
    ],
    "shed":    [value(m.Shed[t])    for t in T],
    "predicted_demand":[load[t] for t in T],
 })
out["total_cost"] = dt * (
    out["coal"]*costs["coal"] +
    out["oil"]*costs["oil"] +
    out["lng"]*costs["lng"] +
    out["th_other"]*costs["th_other"] +
    out["biomass"]*costs["biomass"] +
    out["misc"]*costs["misc"] +
    out["curtail_pv"]*c_curt_pv +
    out["curtail_wind"]*c_curt_wind +
    out["shed"]*c_shed +
    out["import"]*c_import
)

明日

結果を可視化してみます！