はじめに

これまで以下の記事で、電力系統における発電機の動揺方程式や、ガバナフリー制御・負荷周波数制御（LFC）による周波数維持のしくみを扱ってきた。

発電機の周波数と電気的出力の関係
単純な電力系統における過渡解析について
電力系統におけるGFとLFCの役割

今回はもっとも単純な「1台の発電機と1つの負荷だけからなる系統」を対象に、ダンパ定数を無視した動揺方程式を差分法（オイラー法）で数値的に解き、発電機出力 $P_S$ を変化させたときに負荷 $P_L$ がどのように追従するか、また系統周波数 $f$ がどのように応答するかを調べた。あわせて、$f$ と $P_L$ の関係から負荷特性（負荷の周波数依存性）を確認する。

系統モデル

発電機1台と負荷1つが、系統（母線）を介して直結しているだけの単純なモデルを考える。発電機は電力 $P_S$ を系統に供給し、負荷は系統から電力 $P_L$ を消費する。両者の電力に過不足が生じると、その差分が系統の周波数 $f$（角周波数 $\omega=2\pi f$）を変化させる。

発電機Gが系統に $P_S$ を供給し、負荷が系統から $P_L$ を消費する。$P_S \neq P_L$ のとき系統周波数 $f$ が変化する。

動揺方程式（ダンパ定数を無視）

発電機の動揺方程式は本来

2H\frac{d\omega}{dt} + D\omega = P_S - P_L

の形で表される。$H$ は慣性定数、$D$ はダンパ定数（制動係数）である。今回はダンパ項を無視し（$D=0$）、さらに簡略化のため慣性の項を定数 $2H$ ではなく瞬時角速度 $\omega$ そのもので置き換え、

\frac{d\omega}{dt} = \frac{P_S - P_L}{\omega}

という式で系統の慣性による周波数変化を表すことにした。$P_S>P_L$（供給過多）なら周波数は上昇し、$P_S<P_L$（供給不足）なら周波数は低下する。

ガバナ特性と負荷特性

発電機側には、周波数が下がれば出力を増やす（上がれば減らす）というガバナの応答特性があり、負荷側にも、周波数が上がれば消費電力が増える（下がれば減る）という負荷の周波数特性がある。これらを周波数偏差 $\Delta f$ に比例する形で

\Delta P_S = -k_1 \Delta f, \qquad \Delta P_L = k_2 \Delta f

と表す。$k_1$ をガバナ定数、$k_2$ を負荷の周波数特性定数とする。この2つの関係が、$P_S$ と $P_L$ の差によって生じた周波数変化を打ち消す方向に働き、系統全体としての慣性・安定性を再現する役割を果たす。

プログラム上は、角周波数 $\omega$ の変化量を使って

\frac{dP_S}{dt} = -2\pi k_1 \frac{d\omega}{dt}, \qquad \frac{dP_L}{dt} = 2\pi k_2 \frac{d\omega}{dt}

という形で実装した。$f=\omega/2\pi$ なので、これは $\Delta f$ に比例して $P_S$ が減り、$P_L$ が増える関係を表している。

差分法による数値計算

上記の微分方程式を、時間刻み $\Delta t$ のオイラー法（前進差分）で数値的に解いた。時刻 $t_n$ における値を添字 $n$ で表すと、更新式は次のようになる。

\omega_{n+1} = \omega_n + \frac{P_{S,n}-P_{L,n}}{\omega_n}\Delta t

P_{S,n+1} = P_{S,n} - 2\pi k_1\frac{P_{S,n}-P_{L,n}}{\omega_n}\Delta t

P_{L,n+1} = P_{L,n} + 2\pi k_2\frac{P_{S,n}-P_{L,n}}{\omega_n}\Delta t

f_n = \frac{\omega_n}{2\pi}

3つの状態量 $P_S, P_L, \omega$ をまとめて1つのベクトルとして扱い、共通の係数ベクトル $(-2\pi k_1,\ 2\pi k_2,\ 1)$ を掛けて一括で更新するという実装にした。

シミュレーション条件

初期状態は $P_L=1,\ f=1,\mathrm{Hz}$（$\omega=2\pi$）としたうえで、発電機出力 $P_S$ を時間に応じて次のように階段状に変化させた。

$t<0$：$P_S=1$（定常状態）
$0\le t<1$：$P_S=0.5$（発電機出力が急減）
$1\le t<9$：$P_S=1$（元の出力に復帰）
$t\ge9$：$P_S=1.5$（発電機出力が急増）

その他のパラメータは、ガバナ定数 $k_1=2$、負荷の周波数特性定数 $k_2=1$、時間刻み $\Delta t=10^{-5},\mathrm{s}$、計算範囲 $-1\le t\le10,\mathrm{s}$ とした。

プログラム

power_omega.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math

#P_S=1
P_L=1
f=1
omega=2*math.pi*f

t=-1
delta_t=1e-5

t_max=10

P_S_ary=[]
P_L_ary=[]
omega_ary=[]
f_ary=[]
t_ary=[]
k_1=2
k_2=1

while(t<t_max):
    
    if t<0:
        P_S=1
    elif t<1:
        P_S=0.5
    elif t<9:
        P_S=1
    else:
        P_S=1.5
    
    A=np.array([P_S,P_L,omega]).T
    A+=(P_S-P_L)/omega*np.array([-k_1*2*math.pi,k_2*2*math.pi,1]).T*delta_t     
    P_S=A[0]
    P_L=A[1]
    omega=A[2]
    P_S_ary.append(P_S)
    P_L_ary.append(P_L)
    omega_ary.append(omega)
    f_ary.append(omega/(2*math.pi))
    t_ary.append(t)
    t+=delta_t
## グラフ1    
plt.plot(t_ary,P_S_ary,label="P_S")
plt.plot(t_ary,P_L_ary,label="P_L")
plt.xlabel("時間[s]")
plt.ylabel("電力")
plt.legend()    
plt.savefig("power_omega.png")
plt.show()
## グラフ2
plt.plot(t_ary,f_ary)
plt.xlabel("時間[s]")
plt.ylabel("周波数")
plt.savefig("f_P_L.png")
plt.show()
## グラフ3
plt.plot(f_ary,P_L_ary,label="P_L")
plt.xlabel("周波数")
plt.ylabel("電力")
plt.legend()
plt.savefig("P_L_f.png")
plt.show()

結果

グラフ1

$P_S$と$P_L$の時間変化

青線が発電機出力 $P_S$、オレンジ線が負荷 $P_L$。

$P_S$ が階段状に変化するのに対し、$P_L$ は瞬時には追従せず、周波数の変化を介して徐々に $P_S$ に近づいていく様子が確認できた。$t=0$ で $P_S$ が $1\to0.5$ に急減すると、$P_L$ は一旦 $P_S$ を下回るまで下がったのち、$t=1$ で $P_S$ が $1$ に戻るまでの間に緩やかに $1$ へ収束する。$t=9$ で $P_S$ が $1.5$ に急増した際も同様に、$P_L$ は遅れながら追いかける挙動を示した。これは、$P_S$ と $P_L$ の差が直接 $P_L$ を動かすのではなく、一旦周波数の変化を引き起こし、その周波数偏差が負荷特性 $\Delta P_L=k_2\Delta f$ を通じて $P_L$ を動かす、という間接的な経路になっているためである。

グラフ2

周波数の時間応答

系統周波数 $f$ の時間変化。

$P_S<P_L$ となる $0\le t<1$ の区間では供給不足のため周波数が低下し、最も低下したところで $f\approx0.992$ まで下がった。$t=1$ で $P_S$ が元に戻ると、供給と消費の差が縮小するにつれて周波数はゆるやかに $f=1$ へ回復する。逆に $t=9$ 以降は $P_S$ が $P_L$ を上回るため周波数は上昇し、$t=10$ の時点で $f\approx1.008$ に達した。周波数の変化の向きが $P_S-P_L$ の符号と一致しており、動揺方程式が意図通りに機能していることが分かる。

グラフ3

負荷$P_L$と周波数$f$の関係（負荷特性）

横軸が周波数 $f$、縦軸が負荷 $P_L$。

$P_L$ と $f$ を散布図的にプロットすると、シミュレーション区間全体にわたってほぼ完全な直線関係になった。これは、モデルに組み込んだ負荷の周波数特性 $\Delta P_L=k_2\Delta f$ がそのまま現れたものであり、系統周波数が上昇するほど負荷の消費電力が線形に増加するという負荷特性を、このシンプルな差分モデルでも再現できていることを確認できた。

まとめ

1発電機-1負荷という最小構成の系統モデルに対し、ダンパ定数を無視した動揺方程式と、ガバナ特性 $\Delta P_S=-k_1\Delta f$・負荷特性 $\Delta P_L=k_2\Delta f$ を組み合わせ、オイラー法による差分計算でシミュレーションを行った。発電機出力 $P_S$ をステップ状に変化させると、負荷 $P_L$ は周波数変化を介して遅れながら追従すること、また周波数と負荷の間には想定通りの線形関係（負荷特性）が成り立つことを確認できた。今回は慣性定数 $H$ を定数として扱わず瞬時角速度 $\omega$ で代用する簡略化を行ったが、この点の妥当性や、$k_1,k_2$ を変化させたときの応答の違いについては今後の記事で扱いたい。

参考文献

発電機の周波数と電気的出力の関係 https://qiita.com/arairuca/items/fa088166c4072cc6ab01
単純な電力系統における過渡解析について https://qiita.com/arairuca/items/93c3fe6f49f3630d8183
電力系統におけるGFとLFCの役割 https://qiita.com/arairuca/items/5ef5fa34768ffb69b5a6