はじめに

同期発電機の動揺方程式は、同期安定性を議論する上で基本的な式である。一般に慣性定数が大きいほど安定性は高く、また電圧位相角が90度を超えると脱調が生じ、送電不能となる。
一方、動揺方程式を解析的に解くことは極めて困難であり、数値計算によるアプローチも無限大の時間を扱えないという本質的な限界を抱えている。そのため、従来の手法では安定性の厳密な評価が難しかった。
そこで本稿では、微分方程式の安定性をベクトル場として可視化・評価する手法を導入し、動揺方程式の安定性解析に応用する。

導入

同期発電機の動揺方程式は、参考文献1より以下のように与えられる。

{\frac{M}{\omega_n}\frac{d\delta^2}{dt^2}+\frac{D}{\omega_n}\frac{d\delta}{dt}=P
}

ただし、


{P=P_m-P_e=P_m-\frac{V_rV_s}{X}\sin \delta
}

ここで、参考文献2のように、横軸を$\delta$とし、縦軸を$\dot{\delta}$として、ベクトル場を作成する。

つまり、


\frac{d}{dt}

\begin{pmatrix}
\delta  \\
\dot{\delta} \\
\end{pmatrix}

=

\begin{pmatrix}
0 & 1 \\
-\frac{V_rV_s}{X}\frac{\omega_n}{M} & -\frac{D}{M} \\
\end{pmatrix}


\begin{pmatrix}
\delta  \\
\dot{\delta} \\
\end{pmatrix}

+

\begin{pmatrix}
0  \\
P_m \\
\end{pmatrix}

の左辺のベクトル場を作成する。

プログラム

python vector_P.py

"""
電力系統の動揺方程式（スイング方程式）の安定性解析 - ベクトル場による可視化
=============================================================================
参考: https://qiita.com/arairuca/items/3d698ce923a086118090

■ 動揺方程式（スイング方程式）とは
  電力系統において、発電機のロータ（回転子）の運動を表す方程式。
  機械的入力パワーと電気的出力パワーの不均衡によって位相角が変動する様子を記述する。

    M * d²δ/dt² + D * dδ/dt = Pm - Pe

  ただし:
    δ  : 位相角 [rad]（系統基準に対する発電機ロータの角度差）
    Pm : 機械的入力パワー [pu]（タービンから与えられる動力）
    Pe : 電気的出力パワー [pu]（Pe = Vr*Vs/(X) * sin(δ)）
    M  : 慣性定数 [s]（ロータの慣性を表す）
    D  : 制動係数（ダンパー巻線や空気摩擦などによる減衰）

■ 状態空間表現（1階連立ODEへの変換）
  δ の微分を d_delta = dδ/dt = ω（角速度）と置くことで
  2階ODEを1階連立ODEに変換する:

    dδ/dt       = ω                              ... (1)
    dω/dt       = -D/M * ω - ωn*Vr*Vs/(X*M) * sin(δ) + Pm/M  ... (2)

  これにより、状態変数 (δ, ω) の2次元位相面（phase plane）上で
  各点の「速度ベクトル」 (dδ/dt, dω/dt) を計算できる。

■ 平衡点（ベクトルが零になる点）
  dδ/dt = 0 かつ dω/dt = 0 となる点:
    - ω = 0
    - sin(δ*) = Pm * X / (ωn * Vr * Vs)

  これを満たす δ は2つ存在する:
    - 安定平衡点（SEP）: δ* = arcsin(Pm*X / (ωn*Vr*Vs))
    - 不安定平衡点（UEP）: δ* = π - arcsin(Pm*X / (ωn*Vr*Vs))
"""

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語フォントをmatplotlibで使用可能にするライブラリ
import math

# ---------------------------------------------------------------------------
# グリッドパラメータ
# ---------------------------------------------------------------------------
# ベクトル場を計算するグリッドの解像度（点の数）
# 大きいほど滑らかだが計算コストが増大する
n = 1000

# 横軸（位相角 δ）の描写範囲: -π ～ +π [rad]
# 位相角はこの範囲で周期的な挙動を示すため、1周期分を表示する
L = math.pi

# x軸方向のグリッド刻み幅（参考値: 実際のメッシュはlinspaceで生成）
dx = 2 * L / (n - 1)

# ---------------------------------------------------------------------------
# 電力系統パラメータ（per unit 系）
# ---------------------------------------------------------------------------
Vr = 1        # 受電端電圧 [pu]（受電側の電圧大きさ）
Vs = 1        # 送電端電圧 [pu]（発電機端の電圧大きさ）
omega_n = 2 * math.pi * 50  # 系統基準角速度 ωn = 2π × 50 [rad/s]（50Hz系統）
M = 0.01       # 慣性定数 [s]（小さいほどロータが素早く加速する）
D = 0.1       # 制動係数（正値が安定化に寄与; 0にすると純粋な保守系）
P_m = 1     # 機械的入力パワー [pu]（タービン出力; 一定と仮定）
XX = 1        # 送電線のリアクタンス X [pu]（Vr と Vs を結ぶ系統インピーダンス）

# ---------------------------------------------------------------------------
# 位相面（Phase Plane）グリッドの生成
# ---------------------------------------------------------------------------
# x軸: 位相角 δ ∈ [-π, π] [rad]
x = np.linspace(-math.pi, math.pi, n)

# y軸: 角速度 ω = dδ/dt ∈ [-2ωn, +2ωn] [rad/s]
#   ωn を基準に ±2倍の範囲を見ることで、安定領域と不安定領域の両方を捉える
y = np.linspace(-2 * omega_n, 2 * omega_n, n)

# 2次元メッシュグリッドを生成 → 各格子点の (δ, ω) 座標を行列で保持
X, Y = np.meshgrid(x, y)

# ベクトル場の各成分（U: δ方向、V: ω方向）を格納する配列を初期化
U = np.zeros((n, n))  # U[p][q] = dδ/dt = ω（位相角の時間変化率）
V = np.zeros((n, n))  # V[p][q] = dω/dt（角速度の時間変化率）

# ---------------------------------------------------------------------------
# 各グリッド点でのベクトル（dδ/dt, dω/dt）を計算
# ---------------------------------------------------------------------------
# 全格子点をループし、状態空間表現の右辺を評価する
for p in range(n):
    for q in range(n):
        delta = X[p][q]    # 現在の格子点の位相角 δ [rad]
        d_delta = Y[p][q]  # 現在の格子点の角速度偏差 ω = dδ/dt [rad/s]

        # 式(1): dδ/dt = ω（定義そのもの）
        U[p][q] = d_delta

        # 式(2): dω/dt = -D/M * ω - ωn*Vr*Vs/(X*M) * sin(δ) + Pm/M
        #   第1項 (-D/M * ω)              : 制動トルク（ω に比例した減衰）
        #   第2項 (-ωn*Vr*Vs/(X*M)*sin(δ)): 電気的同期トルク（復元力）
        #   第3項 (Pm/M)                  : 機械的駆動トルク（一定入力）
        V[p][q] = (
            -D / M * d_delta
            - omega_n * Vr * Vs / (XX * M) * math.sin(delta)
            + P_m / M
        )

# ---------------------------------------------------------------------------
# ベクトル場（流線）のプロット
# ---------------------------------------------------------------------------
# streamplot: 各点のベクトルを流線として描画
#   color    : 流線の色
#   density  : 流線の密度（大きいほど線が多い）
#   linewidth: 流線の太さ
#   arrowsize: 矢印の大きさ
plt.streamplot(X, Y, U, V, color='blue', density=1.5, linewidth=0.5, arrowsize=1.5)

# 平衡点のプロット
# ここでは δ = π/2, ω = ωn の点を安定限界として示す
# （実際の安定平衡点は ω=0, δ=arcsin(Pm*X/(ωn*Vr*Vs)) の点）
plt.plot(math.pi / 2, omega_n, 'ro', label='安定限界', color='red')

# 軸ラベル・凡例・保存・表示
plt.xlabel('位相角 δ [rad]')
plt.ylabel('角速度 ω [rad/s]')
plt.title('動揺方程式の位相面（ベクトル場）')
plt.legend()
plt.tight_layout()
plt.savefig('動揺方程式_安定性.png', dpi=150)
plt.show()

結果

これを実行すると、以下のようになる。

赤点の部分（通常周波数かつ脱調ぎりぎりの位相角）周辺で、安定か不安定に分岐する。

まとめ

本稿では、ベクトル場を用いた動揺方程式の安定性解析を行った。この手法は方程式を直接解くのではなく位相空間上のベクトル場として振る舞いを把握するため、計算コストが低く、数値解法のように解の誤差や妥当性を別途検証する必要もない。安定性の定性的な評価を目的とする場面において、特に有効なアプローチといえる。

参考文献

1.同期発電機の動揺方程式と数値計算

2.微分方程式とベクトル場