一機無限大母線系統における過渡安定度シミュレーション（3相短絡事故と等面積法）

はじめに

電力系統の過渡安定度は、送電線の3相短絡事故のような大擾乱が発生したときに、発電機が同期を保ったまま運転を継続できるかを評価する問題である。この評価には、発電機の回転運動を表す動揺方程式という非線形微分方程式を解く必要があるが、コンピュータが普及する以前は手計算でこれを解くのは現実的でなかった。そこで、等面積法と呼ばれる、P-δ平面上の面積比較だけで安定・不安定を判定するグラフィカルな手法が広く使われてきた。
　　本稿では、このような等面積法を用いることが可能な一機無限大母線系統(2回線送電線)を対象に、

1. 動揺方程式を数値積分(RK4)して、3相短絡事故に対する位相角δ・電気出力Peの時間応答を求める
2. 同じ系統について等面積法による安定判定(臨界遮断角・臨界遮断時間)を解析的に求める
3. 両者が一致することを確認する

というシミュレーションをPythonで実装する。あわせて、「事故中は系統インピーダンスが大きくなるのか小さくなるのか」という、初学者が混乱しやすい論点についても整理する。

参考: 電気の真髄「過渡安定度と等面積法」

モデル: 一機無限大母線系統(2回線送電線)

対象とするモデルは、同期発電機(過渡リアクタンス$X_d'$)→変圧器($X_t$)→2回線並列の送電線($X_l$)→無限大母線、という構成である。

\text{事故前} \quad X_0 = X_d' + X_t + \frac{X_l}{2} \qquad (\text{2回線並列})

\text{事故除去後} \quad X_2 = X_d' + X_t + X_l \qquad (\text{事故回線を遮断し1回線運用})

事故が発生すると、片方の回線上で3相短絡が生じ、一定時間後(遮断時間)に保護リレーがその回線を系統から切り離す、という一連の流れを想定する。事故中のリアクタンス(発電機内部起電力Eと無限大母線電圧Vを結ぶ等価伝達リアクタンス)については後述する。

発電機の内部起電力$E$・無限大母線電圧$V$は事故中も一定とみなし、電気出力は

P_e(\delta) = \frac{EV}{X}\sin\delta

で与えられる。$X$は上記3つの段階(事故前・事故中・事故除去後)でそれぞれ異なる値をとる。

動揺方程式と数値シミュレーション(RK4)

動揺方程式は

\frac{2H}{\omega_s}\frac{d^2\delta}{dt^2} = P_m - P_e(\delta)

であり、$H$は慣性定数[s]、$\omega_s$は同期角周波数[rad/s]である。これを1階の連立微分方程式に書き直すと

\frac{d\delta}{dt} = \omega, \qquad \frac{d\omega}{dt} = \frac{\omega_s}{2H}\bigl(P_m - P_e(\delta)\bigr)

となる。ここで$\omega$は同期速度からの速度偏差であり、$\delta$と$\omega$の2つの状態変数を持つ2階のODEであることに注意する。実装の初期段階で、この速度偏差$\omega$を状態として持たず、加速度$d\omega/dt$を$\delta$へ直接積分してしまう(=慣性による「勢い」が消えてしまう)という誤りをしやすい。角運動量を持つ振り子のシミュレーションと同様、必ず位置(δ)と速度(ω)の2状態をペアで時間発展させる必要がある。

本実装ではRK4(4次のルンゲ=クッタ法)で数値積分する。

keitou_toumenseki.py(抜粋)

def deriv(state, t):
    delta, omega = state
    ddelta = omega
    domega = (omega_s / (2 * H)) * (Pm - Pe_of(delta, t))
    return np.array([ddelta, domega])

state = np.array([delta0, 0.0])
t = 0.0
for i in range(n_steps):
    k1 = deriv(state, t)
    k2 = deriv(state + 0.5 * dt * k1, t + 0.5 * dt)
    k3 = deriv(state + 0.5 * dt * k2, t + 0.5 * dt)
    k4 = deriv(state + dt * k3, t + dt)
    state = state + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
    t += dt

パラメータは以下のとおり(系統周波数60Hz、$H=5,\mathrm{s}$、$P_m=0.8,\mathrm{pu}$、$E=1.1,\mathrm{pu}$、$V=1.0,\mathrm{pu}$、$X_d'=0.3$、$X_t=0.1$、$X_l=0.5$)。

keitou_toumenseki.py(パラメータ)

f0 = 60.0
omega_s = 2 * np.pi * f0
H = 5.0
Pm = 0.8

E = 1.1
V = 1.0
Xd = 0.3
Xt = 0.1
Xl = 0.5

シミュレーション結果: δ・Pe の時間応答

事故発生時刻を$t=1.0,\mathrm{s}$、遮断時間を0.15s(=$t=1.15,\mathrm{s}$で事故除去)とした場合の応答を示す。

位相角δの時間応答。事故発生と同時にδが急上昇し始め、遮断後も慣性でしばらく増加を続けた後、ピークを迎えて減少に転じ、ダンパ項(制動係数D)を含めていないため減衰せずに一定振幅で振動を続ける。

電気出力Peの時間応答。事故発生と同時にPeがほぼ0まで落ち込み(後述)、遮断の瞬間にPe2=Pmax2 sinδの曲線へ不連続にジャンプし、その後はδの振動に応じてPeも振動する。

速度偏差ωの時間応答。事故中はPm-Pe>0となる一方的な加速により単調に増加し、遮断後はPe2>Pmが優勢な間だけ減速に転じる。

この3枚のグラフから、事故に対する系統の応答は次の3段階で理解できる。

• 事故前: $\delta=\delta_0$、$P_e=P_m$で釣り合った定常運転点にある。
• 事故中(0.10〜0.15秒間): 発電機端子近傍の3相ボルト短絡を仮定しているため$P_e\approx0$となり、機械入力$P_m$がそのままロータを加速するトルクになる。この間、$\delta$は加速度一定の放物線的な増加(等加速度運動)を示し、$\omega$は単調増加する。
• 事故除去後: 事故回線が遮断され1回線運用の$P_{e2}=P_{max2}\sin\delta$曲線に切り替わる。遮断時の$\delta$はすでに$\delta_0$より進んでいるため、$P_{e2}>P_m$となり減速トルクが働く。このシミュレーションでは十分早く遮断できているため$\delta$はピークを打って減速に転じ、その後は(制動項を入れていないので)一定振幅の等時的な振動を続ける。

等面積法とは

上記の時間応答を得るには非線形微分方程式を数値積分する必要があるが、現代のようにコンピュータで数値計算するのではなく、かつては次の等面積法(グラフィカルな手法)で安定性を評価していた。

P-δ平面上で、事故中は$P_m>P_e$となり加速エネルギー

S_1 = \int_{\delta_0}^{\delta_c} (P_m - P_e)\, d\delta

が$\delta_0$から遮断角$\delta_c$までロータに蓄積される。事故除去後は$P_e>P_m$となる区間があり、減速エネルギー

S_2 = \int_{\delta_c}^{\delta_{max}} (P_e - P_m)\, d\delta

がこれを打ち消す方向に働く。$S_1 \le S_2$を満たせる限り(=$\delta$が事故除去後曲線の不安定平衡点$\delta_{max}=\pi-\arcsin(P_m/P_{max2})$に達する前に減速エネルギーを使い切れる限り)ロータは脱調せず、$S_1=S_2$となる$\delta_c$が臨界遮断角である。

\cos\delta_c = \frac{P_m(\delta_{max}-\delta_0) + P_{max2}\cos\delta_{max} - P_{max1}\cos\delta_0}{P_{max2}-P_{max1}}

等面積法の意義は、動揺方程式そのものを解かなくても、$P_{max0}, P_{max1}, P_{max2}$という3本のsinカーブの形状(=事故前・事故中・事故除去後のリアクタンス)さえ分かれば、代数的に安定限界(臨界遮断角、さらに事故中は$P_e\approx0$の等加速度運動とみなせるので臨界遮断時間$t_c$も)を求められる点にある。手回し計算機の時代に多入力・多出力の非線形微分方程式を直接解くことは非現実的であったため、この「エネルギーの釣り合い」という物理的な洞察に基づく近似手法が実用上大きな価値を持っていた。

現代では、本稿のようにRK4等で動揺方程式を直接数値積分することが容易であり、時々刻々のδ, ωの波形もそのまま得られる。そのため実務上、等面積法だけに頼る場面は少なくなったが、①手計算でも安定余裕をすぐ見積もれる、②数値シミュレーション結果の妥当性をチェックする物理的な指標になる、という点で今なお有用な考え方である。本稿でも、両者(数値積分による時間応答と、等面積法による解析的な臨界値)を突き合わせて一致することを確認する。

P-δ図と等面積法による安定判定

上記の系統パラメータから、

P_{max0}=\frac{EV}{X_0}=1.692,\quad P_{max1}=0,\quad P_{max2}=\frac{EV}{X_2}=1.222 \ [\mathrm{pu}]

\delta_0=28.2^\circ,\quad \delta_{max}=139.1^\circ,\quad \delta_c=59.3^\circ,\quad t_c=0.190\,\mathrm{s}

が得られる。遮断時間0.15s(< $t_c$)としたケースのP-δ図を以下に示す。

橙色がS1(加速エネルギー)、青色がS2(減速可能量)。遮断時の$\delta$(≈48.0°)における$S_1=0.276$は、$\delta_{max}$(≈139.1°)までに得られる$S_2=0.470$より小さいため安定と判定できる。

実際に遮断時間を変えてシミュレーションすると、この判定は数値積分の結果とも一致する。

青(遮断0.15s)は臨界遮断角(破線)を一瞬超えるが振動しながら減速し脱調しない。赤(遮断0.30s)は$t_c=0.190\mathrm{s}$を超えて遮断が遅れたため$\delta$が単調に増加し続け、脱調(不安定)に至る。

遮断時間	遮断時δ	S1	S2	判定
0.15 s	48.0°	0.276	0.470	安定
0.30 s	106.6°	1.094	0.121	不安定(脱調)

等面積法による解析値($t_c\approx0.190\mathrm{s}$)と、数値積分で実際に脱調が起きる境界がほぼ一致しており、両アプローチの整合性が確認できる。

事故中の系統インピーダンスは大きくなるのか、小さくなるのか

ここまでのシミュレーションでは「事故中は$P_e\approx0$」、言い換えれば「事故中はPeを計算する式のXが非常に大きくなる」という仮定を置いた。この点は直感に反すると感じられやすいので整理しておく。「インピーダンス」という言葉が指している対象が2つあることに注意すると分かりやすい。

① 短絡電流を計算するための(事故点から見た)インピーダンス

これはボルト短絡であればほぼ0Ωである。だからこそ短絡電流は定格電流の何倍〜何十倍にも達する。このインピーダンスは事故前・事故除去後と比べて小さくなる。短絡容量計算やリレー整定で使う「系統インピーダンス」は基本的にこちらを指す。

② 動揺方程式・等面積法で使う、発電機内部起電力Eと無限大母線電圧Vを結ぶ等価伝達(移相)リアクタンスX

$P_e=EV\sin\delta/X$の$X$はこちらであり、事故中は事故前より大きくなる(発電機端子近傍の完全な3相ボルト短絡であれば理論上は無限大、$P_{max}=0$とみなせる)。本稿のシミュレーションが仮定しているのはこちらである。

同じ「X」という文字を使うため混同しやすいが、①は「事故点をグラウンドさせたときに、事故点から電源側を見た合成インピーダンス」であり、②は「電源(E)から受電端(V)まで有効電力を運ぶ経路の実効的なインピーダンス」であって、指している回路上の区間がそもそも異なる。

これを踏まえると、筆者や読者が抱く質問の各項目には次のように答えられる。

• 「短絡時は母線-対地間の抵抗が極小になるから、受電側の有効電力は大きくなる?」 → 逆である。対地インピーダンスがほぼ0になるということは、その母線の電圧もほぼ0まで潰れるということであり、有効電力$P=VI\cos\varphi$は電圧がほぼ0では電流がいくら大きくてもほとんど伝わらない。したがって受電側(無限大母線側)に届く有効電力はむしろ激減する。
• 「短絡事故時は無効電力のみが事故点に流れるので、発電機から有効電力はほとんど供給されないの?」 → 概ねその通り。送電系統の主要素(発電機の過渡リアクタンス・変圧器・送電線)は$R\ll X$の純リアクタンス性とみなせるため、事故電流は位相が約90°遅れた無効性の電流が支配的になる。発電機は変わらず機械入力$P_m$相当のトルクで駆動され続けているのに、外部へ取り出せる有効電力$P_e$がほぼ0になるため、その差分$P_m-P_e$がロータを加速するトルクとなる。これがまさに等面積法の加速エネルギー$S_1$の物理的な起源である。
• 「したがって事故時の系統インピーダンスXは極大となる?」 → 動揺方程式で使う②の意味のXについてはその通りである。発電機端子ごく近傍の完全な3相ボルト短絡であれば$X\to\infty$($P_{max1}=0$)とみなせる。事故点が送電線の途中(発電機側から離れた地点)にある場合は、生成される事故点・発電機・無限大母線の3端子網をY-Δ変換で消去することで、事故前より大きいが有限の等価Xが得られる。いずれの場合も、①の意味の短絡インピーダンス(小さい)とは逆に、②の意味のXは事故前より必ず大きくなる。

まとめ

一機無限大母線系統(2回線送電線)を対象に、動揺方程式をRK4で数値積分することで、3相短絡事故に対するδ・ω・Peの時間応答を再現した。事故中は電気出力Peがほぼ0まで落ち込み、その間ロータが機械入力Pmでほぼ一定加速度で加速される一方、事故除去後は1回線運用の新しいP-δ曲線に切り替わり、遮断が十分早ければ減速に転じて脱調を免れることを確認した。

さらに、この応答は等面積法(S1=S2となる臨界遮断角・臨界遮断時間)による解析的な安定判定とも整合しており、遮断時間0.15s(< 臨界遮断時間0.190s)では安定、0.30s(> 臨界遮断時間)では不安定という数値積分の結果を、等面積法の面積比較だけからも事前に予測できることを示した。等面積法は、コンピュータで直接微分方程式を解ける現代においても、事故除去時間に対する安定余裕を素早く見積もったり、数値シミュレーション結果を検算したりするための物理的な指標として引き続き有用である。

また、事故中に動揺方程式で使う実効的な系統リアクタンスXが(短絡電流計算上のインピーダンスとは逆に)大きくなる理由についても整理した。3相ボルト短絡では事故点電圧がほぼ0まで潰れるため、事故電流自体は非常に大きくなる一方、それはほぼ無効性の電流であり、発電機から系統の反対側(無限大母線側)へ運ばれる有効電力はむしろ激減する。この「有効電力が運べなくなる」という現象を、動揺方程式の枠組みでは「実効的なXが大きくなる(=Pmaxが小さくなる)」という形で表現している。

参考文献

• 電気の真髄「過渡安定度と等面積法」 https://denki-no-shinzui.com/transient-stability/

付録: プログラム全文

keitou_toumenseki.py

"""
等面積法(equal area criterion)による過渡安定度シミュレーション

参考: https://denki-no-shinzui.com/transient-stability/

モデル: 同期発電機 - 変圧器 - 二回線送電線 - 無限大母線
  事故前  : X0 = Xd + Xt + Xl/2         (二回線並列)
  事故中  : 発電機端子近傍の三相短絡を仮定 -> Pe1 = 0 (r1=0)
  事故除去後: X2 = Xd + Xt + Xl           (事故回線を遮断し一回線運用)

動揺方程式:
  d(delta)/dt = omega
  d(omega)/dt = (omega_s / (2H)) * (Pm - Pe(delta, t))
"""

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

# ---------------- パラメータ ----------------
f0 = 60.0                    # 系統周波数 [Hz]
omega_s = 2 * np.pi * f0     # 同期角周波数 [rad/s]
H = 5.0                      # 慣性定数 [s]
Pm = 0.8                     # 機械入力(一定) [pu]

E = 1.1                      # 発電機内部起電力 [pu]
V = 1.0                      # 無限大母線電圧 [pu]

Xd = 0.3                     # 発電機過渡リアクタンス Xd'
Xt = 0.1                     # 変圧器リアクタンス
Xl = 0.5                     # 送電線1回線分のリアクタンス

t_fault_on = 1.0             # 事故発生時刻 [s]
t_clear = t_fault_on + 0.15  # 事故除去(遮断)時刻 [s] ここを動かして安定/不安定を比較する

dt = 0.001
t_max = 6.0

# ---------------- 各状態のリアクタンス・最大電力 ----------------
X0 = Xd + Xt + Xl / 2         # 事故前
X2 = Xd + Xt + Xl             # 事故除去後(1回線運用)

Pmax0 = E * V / X0
Pmax1 = 0.0                   # 事故中(発電機端子近傍三相短絡を仮定)
Pmax2 = E * V / X2

def Pe_of(delta, t):
    if t < t_fault_on:
        return Pmax0 * np.sin(delta)
    elif t < t_clear:
        return Pmax1 * np.sin(delta)
    else:
        return Pmax2 * np.sin(delta)

# ---------------- 等面積法(解析計算) ----------------
delta0 = np.arcsin(Pm / Pmax0)                 # 事故前運転点
delta_max = np.pi - np.arcsin(Pm / Pmax2)      # 事故除去後曲線の不安定平衡点(上限角)

r1 = Pmax1 / Pmax0
r2 = Pmax2 / Pmax0

# delta_c = cos^-1[ (Pm(delta_max-delta0) + Pmax2*cos(delta_max) - Pmax1*cos(delta0)) / (Pmax2-Pmax1) ]
cos_delta_c = (Pm * (delta_max - delta0) + Pmax2 * np.cos(delta_max) - Pmax1 * np.cos(delta0)) / (Pmax2 - Pmax1)

if -1.0 <= cos_delta_c <= 1.0:
    delta_c = np.arccos(cos_delta_c)
    # 事故中は Pe1=0 なので d2delta/dt2 = omega_s*Pm/(2H) (一定加速度) から臨界遮断時間を解析的に求める
    t_c = np.sqrt(4 * H * (delta_c - delta0) / (omega_s * Pm))
else:
    delta_c = None
    t_c = None
    print("警告: この系統では等面積条件を満たす臨界遮断角が存在しません(常に安定 or 常に不安定)。")

def S1_area(delta_a, delta_b):
    """加速エネルギー: delta_a(事故前運転点) -> delta_b(遮断時角度)"""
    return Pm * (delta_b - delta_a) + Pmax1 * (np.cos(delta_b) - np.cos(delta_a))

def S2_area(delta_a, delta_b):
    """減速エネルギー: delta_a(遮断時角度) -> delta_b(到達可能な最大角)"""
    return Pmax2 * (np.cos(delta_a) - np.cos(delta_b)) - Pm * (delta_b - delta_a)

# ---------------- 時間domainシミュレーション(RK4) ----------------
def deriv(state, t):
    delta, omega = state
    ddelta = omega
    domega = (omega_s / (2 * H)) * (Pm - Pe_of(delta, t))
    return np.array([ddelta, domega])

n_steps = int(t_max / dt)
t_ary = np.zeros(n_steps)
delta_ary = np.zeros(n_steps)
omega_ary = np.zeros(n_steps)
Pe_ary = np.zeros(n_steps)

state = np.array([delta0, 0.0])
t = 0.0
for i in range(n_steps):
    t_ary[i] = t
    delta_ary[i] = state[0]
    omega_ary[i] = state[1]
    Pe_ary[i] = Pe_of(state[0], t)

    k1 = deriv(state, t)
    k2 = deriv(state + 0.5 * dt * k1, t + 0.5 * dt)
    k3 = deriv(state + 0.5 * dt * k2, t + 0.5 * dt)
    k4 = deriv(state + dt * k3, t + dt)
    state = state + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
    t += dt

# 遮断時点の実際の位相角(シミュレーション値)
idx_clear = np.searchsorted(t_ary, t_clear)
delta_at_clear = delta_ary[idx_clear]

S1_actual = S1_area(delta0, delta_at_clear)
S2_available = S2_area(delta_at_clear, delta_max)
stable = S1_actual <= S2_available

# ---------------- 結果表示 ----------------
print(f"Pmax0={Pmax0:.3f}, Pmax1={Pmax1:.3f}, Pmax2={Pmax2:.3f} [pu]")
print(f"delta0    = {np.degrees(delta0):.2f} deg")
print(f"delta_max = {np.degrees(delta_max):.2f} deg (事故除去後曲線の不安定平衡点)")
if delta_c is not None:
    print(f"delta_c(臨界遮断角) = {np.degrees(delta_c):.2f} deg")
    print(f"t_c(臨界遮断時間)   = {t_c:.4f} s (事故発生から)")
print(f"設定した遮断時間      = {t_clear - t_fault_on:.4f} s -> delta_at_clear = {np.degrees(delta_at_clear):.2f} deg")
print(f"S1(加速エネルギー) = {S1_actual:.4f},  S2(利用可能な減速エネルギー) = {S2_available:.4f}")
print("=> 安定" if stable else "=> 不安定(脱調)")

# ---------------- プロット ----------------
# (1) delta - t
plt.figure()
plt.plot(t_ary, np.degrees(delta_ary))
plt.axvspan(t_fault_on, t_clear, color="red", alpha=0.15, label="事故継続中")
if delta_c is not None:
    plt.axhline(np.degrees(delta_c), color="gray", linestyle="--", linewidth=1, label="臨界遮断角 δc")
plt.xlabel("t [s]")
plt.ylabel("delta [deg]")
plt.title("位相角の時間応答")
plt.legend()
plt.grid()

# (2) omega - t
plt.figure()
plt.plot(t_ary, omega_ary)
plt.axvspan(t_fault_on, t_clear, color="red", alpha=0.15)
plt.xlabel("t [s]")
plt.ylabel("omega (d(delta)/dt) [rad/s]")
plt.title("速度偏差の時間応答")
plt.grid()

# (3) P-delta 図 + 等面積(S1, S2)
delta_plot = np.linspace(0, np.pi, 500)
plt.figure()
plt.plot(np.degrees(delta_plot), Pmax0 * np.sin(delta_plot), label="事故前 Pe0")
plt.plot(np.degrees(delta_plot), Pmax1 * np.sin(delta_plot), label="事故中 Pe1")
plt.plot(np.degrees(delta_plot), Pmax2 * np.sin(delta_plot), label="事故除去後 Pe2")
plt.axhline(Pm, color="black", linewidth=1, label="Pm")

# S1: delta0 -> delta_at_clear, Pm と Pe1 の間
d_S1 = np.linspace(delta0, delta_at_clear, 100)
plt.fill_between(np.degrees(d_S1), Pmax1 * np.sin(d_S1), Pm, color="tab:orange", alpha=0.4, label="S1(加速)")

# S2: delta_at_clear -> delta_max, Pe2 と Pm の間
d_S2 = np.linspace(delta_at_clear, delta_max, 100)
plt.fill_between(np.degrees(d_S2), Pm, Pmax2 * np.sin(d_S2), color="tab:blue", alpha=0.4, label="S2(減速可能量)")

for d, label in [(delta0, "δ0"), (delta_at_clear, "δ(遮断時)"), (delta_max, "δmax")]:
    plt.axvline(np.degrees(d), color="gray", linestyle=":", linewidth=1)
    plt.text(np.degrees(d), Pmax0 * 1.02, label, ha="center")

plt.xlabel("delta [deg]")
plt.ylabel("P [pu]")
plt.title("P-δ 曲線と等面積法")
plt.legend()
plt.grid()

plt.show()