ナイキスト線図・根軌跡・安定余裕の読み方とPython実装：ゲイン余裕・位相余裕で閉ループ安定性を評価する

ナイキスト線図・根軌跡・ゲイン余裕/位相余裕による制御系の安定性評価を、伝達関数の数学的基礎からscipy.signal.freqresp・numpy.rootsによるPython実装まで解説。開ループ L(s)=K/(s(s+1)(s+2)) を題材に、位相余裕53.4°・ゲイン余裕15.6dB・安定限界ゲインK=6を数値実験で確認します。

4 min read

はじめに

ボード線図は、周波数応答（振幅・位相）を対数スケールで読み解くためのツールでした。しかし「このフィードバック制御系は安定か？」「ゲインをどこまで上げてよいか？」という設計上の中心的な問いに答えるには、開ループ特性から閉ループ安定性を判定する枠組みが必要です。

本記事では、その 3 本柱である ナイキスト線図・根軌跡・安定余裕（ゲイン余裕・位相余裕） を、伝達関数の数学的基礎から Python 実装まで体系的に解説します。題材として、制御工学の教科書で頻出する開ループ伝達関数

\[ L(s) = \frac{K}{s(s+1)(s+2)} \]

を使い、位相余裕・ゲイン余裕・安定限界ゲインを数値実験で確認します。control パッケージに依存せず、scipy.signal と numpy のみで実装します。

フィードバック系と閉ループ安定性

開ループ伝達関数 \(L(s) = C(s)P(s)\) （コントローラ \(C\) とプラント \(P\) の直列）に対し、単位負帰還を施した閉ループ伝達関数は

\[ T(s) = \frac{L(s)}{1 + L(s)} \]

です。閉ループ系が安定であるための必要十分条件は、特性方程式

\[ 1 + L(s) = 0 \]

のすべての根（＝閉ループ極）が左半平面（実部が負）にあることです。問題は、\(L(s)\) の開ループ特性という「測りやすい量」から、閉ループの安定性という「知りたい量」をどう導くかにあります。ナイキスト線図・根軌跡・安定余裕は、いずれもこの橋渡しを行う手法です。

安定余裕：ゲイン余裕と位相余裕

安定余裕は「安定な系が不安定になるまでにどれだけ余裕があるか」を定量化します。開ループ周波数応答 \(L(j\omega)\) を使って 2 つの指標を定義します。

ゲイン余裕（Gain Margin, GM）

位相が \(-180^\circ\) になる周波数を位相交差周波数 \(\omega_{pc}\) と呼びます。この周波数で開ループゲインをあと何倍まで上げると \(|L(j\omega_{pc})| = 1\) に達する（＝不安定になる）かを表すのがゲイン余裕です。

\[ \mathrm{GM} = \frac{1}{|L(j\omega_{pc})|}, \quad \mathrm{GM}[\mathrm{dB}] = -20\log_{10}|L(j\omega_{pc})| \]

位相余裕（Phase Margin, PM）

ゲインが \(1\) （\(0\,\mathrm{dB}\) ）になる周波数をゲイン交差周波数 \(\omega_{gc}\) と呼びます。この周波数で位相があと何度回れば \(-180^\circ\) に達するかを表すのが位相余裕です。

\[ \mathrm{PM} = 180^\circ + \angle L(j\omega_{gc}) \]

実務では 位相余裕 \(45^\circ \sim 60^\circ\) が安定性と応答速度のバランスが取れた目安とされます。位相余裕が小さいほど、閉ループ応答は振動的（オーバーシュート大）になります。

ナイキスト線図

定義

ナイキスト線図は、複素平面上に \(L(j\omega)\) を \(\omega : -\infty \to +\infty\) にわたってプロットした軌跡です。ボード線図が振幅と位相を別々の 2 枚に分けるのに対し、ナイキスト線図は 1 枚の複素平面上に軌跡として描きます。

ナイキストの安定判別法

閉ループ安定性は、点 \(-1 + 0j\) （臨界点）まわりの軌跡の回転数で判定できます。ナイキストの安定判別法は次のように述べられます。

\[ Z = N + P \]
  • \(Z\) ：閉ループ極のうち右半平面にある個数（\(Z = 0\) なら安定）
  • \(P\) ：開ループ極のうち右半平面にある個数
  • \(N\) ：ナイキスト軌跡が臨界点 \(-1\) を時計回りに回る回数

\(L(s)\) が開ループ安定（\(P = 0\) ）なら、軌跡が \(-1\) を囲まなければ（\(N = 0\) ）閉ループも安定です。安定余裕は、この臨界点に軌跡がどれだけ近いかを別の角度から数値化したものと解釈できます。

根軌跡

定義

根軌跡（root locus）は、ゲイン \(K\) を \(0 \to \infty\) と変化させたときに、閉ループ極（特性方程式 \(1 + K L_0(s) = 0\) の根）が複素平面上をどう動くかを描いた軌跡です。\(L(s) = K L_0(s)\) と書くと、特性方程式は

\[ \underbrace{\mathrm{den}(s)}_{L_0 \text{の分母}} + K \cdot \underbrace{\mathrm{num}(s)}_{L_0 \text{の分子}} = 0 \]

となり、各 \(K\) についてこの多項式の根を求めれば軌跡が得られます。\(K = 0\) では根は開ループ極に一致し、\(K\) を増やすと極が移動します。極が虚軸を横切る瞬間のゲインが安定限界であり、そこで系は持続振動（マージナル安定）になります。

Python 実装

周波数応答から安定余裕を計算

scipy.signal.freqresp で開ループ周波数応答を求め、ゲイン交差・位相交差を数値的に探して安定余裕を算出します。

import numpy as np
from scipy import signal

# 開ループ伝達関数 L(s) = 1 / (s(s+1)(s+2))   (K=1)
num = [1.0]
den = np.polymul(np.polymul([1, 0], [1, 1]), [1, 2])  # s(s+1)(s+2) = s^3+3s^2+2s
sys = signal.TransferFunction(num, den)

# 周波数応答
w = np.logspace(-2, 2, 2000)
w, H = signal.freqresp(sys, w)
mag = np.abs(H)
phase = np.unwrap(np.angle(H))

# ゲイン交差周波数（|L|=1）→ 位相余裕
idx_gc = np.argmin(np.abs(mag - 1.0))
wgc = w[idx_gc]
pm = 180 + np.degrees(phase[idx_gc])
print(f"ゲイン交差周波数 wgc = {wgc:.4f} rad/s,  位相余裕 PM = {pm:.2f} deg")

# 位相交差周波数（∠L=-180°）→ ゲイン余裕
idx_pc = np.argmin(np.abs(phase - (-np.pi)))
wpc = w[idx_pc]
gm_db = -20 * np.log10(mag[idx_pc])
print(f"位相交差周波数 wpc = {wpc:.4f} rad/s,  ゲイン余裕 GM = {gm_db:.2f} dB")

実行結果は次の通りです。

ゲイン交差周波数 wgc = 0.4455 rad/s,  位相余裕 PM = 53.43 deg
位相交差周波数 wpc = 1.4160 rad/s,  ゲイン余裕 GM = 15.59 dB

理論値と照合してみましょう。\(L(j\omega) = \dfrac{1}{j\omega(j\omega+1)(j\omega+2)}\) の位相が \(-180^\circ\) になるのは、分母の虚部が \(0\) になる \(\omega_{pc} = \sqrt{2} \approx 1.414\) のときです。このとき \(|L(j\omega_{pc})| = \dfrac{1}{\sqrt{2}\cdot\sqrt{3}\cdot\sqrt{6}} = \dfrac{1}{6}\) なので、ゲイン余裕は \(20\log_{10} 6 \approx 15.56\,\mathrm{dB}\) 。数値計算の \(15.59\,\mathrm{dB}\) とよく一致します。ゲイン余裕が \(6\) 倍ということは、\(K\) を \(6\) 倍まで上げると安定限界に達することを意味します。

ナイキスト線図の描画

import matplotlib.pyplot as plt

# 正の周波数と負の周波数（複素共役）で一周ぶんを描く
w_pos = np.logspace(-2, 2, 2000)
_, H_pos = signal.freqresp(sys, w_pos)

fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(H_pos.real, H_pos.imag, "b", label="ω > 0")
ax.plot(H_pos.real, -H_pos.imag, "b--", label="ω < 0")  # 実軸対称
ax.plot(-1, 0, "rx", markersize=12, label="臨界点 -1")   # 臨界点
ax.set_xlabel("Real")
ax.set_ylabel("Imaginary")
ax.set_title("Nyquist plot: L(s)=1/(s(s+1)(s+2))")
ax.grid(True)
ax.axhline(0, color="k", lw=0.5)
ax.axvline(0, color="k", lw=0.5)
ax.legend()
ax.set_xlim(-1.5, 0.5)
ax.set_ylim(-2, 2)
plt.tight_layout()
plt.show()

\(K = 1\) では軌跡が臨界点 \(-1\) を囲まないため（\(N = 0\) 、\(P = 0\) より \(Z = 0\) ）、閉ループは安定です。軌跡が実軸を横切る点が \(-1/6 \approx -0.167\) であり、ここからゲインを \(6\) 倍すればちょうど \(-1\) に到達する——ナイキスト線図とゲイン余裕が同じ事実を別表現していることが読み取れます。

根軌跡の描画

特性方程式 \(\mathrm{den}(s) + K \cdot \mathrm{num}(s) = 0\) の根を、\(K\) を掃引しながら numpy.roots で求めます。

Ks = np.linspace(0, 12, 400)
roots_all = []
for K in Ks:
    charpoly = den.astype(float).copy()
    charpoly[-len(num):] += K * np.array(num, dtype=float)
    roots_all.append(np.sort_complex(np.roots(charpoly)))
roots_all = np.array(roots_all)

fig, ax = plt.subplots(figsize=(6, 6))
for i in range(roots_all.shape[1]):
    ax.plot(roots_all[:, i].real, roots_all[:, i].imag, ".", ms=2)
ax.axvline(0, color="r", lw=1, ls="--", label="虚軸（安定限界）")
ax.set_xlabel("Real")
ax.set_ylabel("Imaginary")
ax.set_title("Root locus")
ax.grid(True)
ax.legend()
plt.tight_layout()
plt.show()

# 安定限界ゲインの確認
for K in [0, 2, 4, 6, 8]:
    charpoly = den.astype(float).copy()
    charpoly[-len(num):] += K * np.array(num, dtype=float)
    r = np.roots(charpoly)
    print(f"K={K:2d}: max Re = {r.real.max():+.3f}  極={np.round(r,3)}")

実行結果は次の通りです。

K= 0: max Re = +0.000  極=[-2. -1.  0.]
K= 2: max Re = -0.239  極=[-2.521+0.j    -0.239+0.858j -0.239-0.858j]
K= 4: max Re = -0.102  極=[-2.796+0.j    -0.102+1.192j -0.102-1.192j]
K= 6: max Re = +0.000  極=[-3.+0.j     0.+1.414j  0.-1.414j]
K= 8: max Re = +0.083  極=[-3.166+0.j     0.083+1.587j  0.083-1.587j]

\(K = 6\) でちょうど一対の極が虚軸上 \(\pm j\sqrt{2}\) に乗り（マージナル安定）、\(K > 6\) で右半平面に入って不安定化します。これは周波数応答から求めたゲイン余裕 \(6\) 倍・位相交差周波数 \(\sqrt{2}\) と完全に一致します。3 つの手法（安定余裕・ナイキスト・根軌跡）が同じ安定限界を指し示すことが数値実験で確認できました。

scipy.signal のみで安定余裕を関数化

再利用しやすいよう、周波数応答から安定余裕を返す関数にまとめておきます。

def stability_margins(num, den, w=None):
    """開ループ (num/den) の GM[dB], PM[deg], wpc, wgc を返す。"""
    if w is None:
        w = np.logspace(-3, 3, 5000)
    sys = signal.TransferFunction(num, den)
    w, H = signal.freqresp(sys, w)
    mag, phase = np.abs(H), np.unwrap(np.angle(H))
    i_gc = np.argmin(np.abs(mag - 1.0))
    i_pc = np.argmin(np.abs(phase - (-np.pi)))
    gm_db = -20 * np.log10(mag[i_pc])
    pm = 180 + np.degrees(phase[i_gc])
    return dict(GM_dB=gm_db, PM_deg=pm, wpc=w[i_pc], wgc=w[i_gc])

print(stability_margins([1.0], [1, 3, 2, 0]))

設計指針

位相余裕からダンピングを見積もる

2 次系の近似では、位相余裕 \(\mathrm{PM}\) と減衰比 \(\zeta\) の間に \(\zeta \approx \mathrm{PM}[\deg] / 100\) という経験則があります。今回の \(\mathrm{PM} = 53.4^\circ\) は \(\zeta \approx 0.53\) に相当し、適度なオーバーシュート（約 \(14\%\) ）を持つ良好な応答が期待できます。位相余裕を大きく取りすぎると応答が鈍く、小さすぎると振動的になるため、\(45^\circ \sim 60^\circ\) が実用上の目安になります。

ボード線図との対応

安定余裕はボード線図上でも直読できます。ゲイン交差周波数での位相と \(-180^\circ\) の差が位相余裕、位相交差周波数でのゲインと \(0\,\mathrm{dB}\) の差がゲイン余裕です。ナイキスト線図は「臨界点への近さ」という 1 枚の絵で両者を統合的に見せ、根軌跡は「ゲインを動かしたときの極の運動」を見せる——3 つは同じ安定性を異なる視点で可視化した相補的な道具です。

離散時間系での注意

ディジタル制御では、\(s\) 平面の左半平面が \(z\) 平面の単位円内に対応します。安定判別の臨界点や虚軸の役割が単位円に置き換わる点に注意が必要です。連続系設計をサンプリングして実装する場合、標本化定理に基づくサンプリング周波数の選定と、離散化に伴う位相遅れの増加（位相余裕の目減り）を考慮します。

まとめ

  • 安定余裕（ゲイン余裕・位相余裕）は、開ループ周波数応答 \(L(j\omega)\) から閉ループ安定性の「余裕」を定量化する。位相余裕 \(45^\circ \sim 60^\circ\) が実用目安。
  • ナイキスト線図は \(L(j\omega)\) を複素平面に描き、臨界点 \(-1\) まわりの回転数で安定判別する（\(Z = N + P\) ）。
  • 根軌跡はゲイン \(K\) を掃引したときの閉ループ極の軌跡で、極が虚軸を横切るゲインが安定限界。
  • 題材 \(L(s) = K/(s(s+1)(s+2))\) では、位相余裕 \(53.4^\circ\) ・ゲイン余裕 \(15.6\,\mathrm{dB}\) ・安定限界 \(K = 6\) の 3 者が完全に整合することを scipy.signal と numpy.roots で確認した。
  • control パッケージがなくても、周波数応答（freqresp）と多項式の根（roots）だけで安定余裕・ナイキスト・根軌跡はすべて実装できる。

関連記事

  • ボード線図の読み方と作成 - 安定余裕を直読できる周波数応答の対数スケール表現。本記事の姉妹編で、ゲイン余裕・位相余裕はボード線図上でも読み取れます。
  • PID制御の基礎 - 安定余裕を確保しながらゲインを調整する代表的なコントローラ設計です。
  • PID制御のPython実装 - 本記事の安定余裕を実際のPIDゲイン調整に応用する実装編です。
  • H∞制御の基礎 - 安定余裕を周波数領域の重み関数として一般化する、ロバスト制御の枠組みです。
  • 標本化定理とエイリアシング - 連続系設計をディジタル実装する際のサンプリング周波数選定の基礎です。
  • デジタルフィルタ設計指針ハブ - 同じ伝達関数・周波数応答の枠組みを、制御ではなく信号フィルタリングに適用した姉妹ハブです。

参考文献

関連ツール

この記事に関連する無料ツール

DevToolBox Moving Average Calculator SMA/EMA/WMA をブラウザで即座に計算 DevToolBox Frequency Converter Hz/kHz/MHz/rad/s 単位を相互変換 DevToolBox JSON Formatter JSONを整形・圧縮・バリデーション

関連記事