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

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

22 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\) ）閉ループも安定です。安定余裕は、この臨界点に軌跡がどれだけ近いかを別の角度から数値化したものと解釈できます。

偏角の原理による厳密な証明

\(Z = N + P\) は経験則ではなく、複素解析の**偏角の原理（Argument Principle）**から厳密に導かれます。天下り的に使う前に、なぜこの式が成り立つのかを示します。

偏角の原理：複素関数 \(f(s)\) が単純閉曲線 \(C\) の内部および周上で有理型（極を除いて正則）であり、\(C\) 上に零点・極を持たないとき、

\[ \frac{1}{2\pi i}\oint_C \frac{f'(s)}{f(s)}\,ds = Z_f - P_f \]

が成り立ちます。ここで \(Z_f, P_f\) はそれぞれ \(C\) の内部にある \(f\) の零点・極の個数（重複度込み）です。

証明：\(C\) の内部で \(f\) が位数 \(m\) の零点 \(s=a\) を持つとすると、\(f(s) = (s-a)^m g(s)\) （\(g(a) \neq 0\) で正則）と書けるので、

\[ \frac{f'(s)}{f(s)} = \frac{m}{s-a} + \frac{g'(s)}{g(s)} \]

であり、右辺第2項は \(s=a\) で正則なので、\(f'/f\) は \(s=a\) に留数 \(m\) の単純極を持ちます。同様に位数 \(n\) の極 \(s=b\) （\(f(s) = (s-b)^{-n} h(s)\) 、\(h(b) \neq 0\) ）では、

\[ \frac{f'(s)}{f(s)} = \frac{-n}{s-b} + \frac{h'(s)}{h(s)} \]

となり、留数は \(-n\) です。留数定理より、\(C\) に沿った周回積分は内部のすべての特異点での留数の和の \(2\pi i\) 倍に等しいので、

\[ \oint_C \frac{f'(s)}{f(s)}\,ds = 2\pi i \sum_k m_k - 2\pi i \sum_l n_l = 2\pi i (Z_f - P_f) \]

が得られ、主張が示されます。さらに \(\dfrac{f'(s)}{f(s)} = \dfrac{d}{ds}\log f(s)\) に注意すると、この積分は \(\log f(s) = \ln|f(s)| + i\arg f(s)\) の \(C\) に沿った変化量に等しく、\(C\) が閉曲線で \(|f(s)|\) は一周して元の値に戻るため実部の変化はゼロ、残るのは偏角の変化だけです。したがって

\[ Z_f - P_f = \frac{1}{2\pi}\Delta_C \arg f(s) \]

——右辺は、\(s\) が \(C\) を一周する間に像 \(f(s)\) が原点のまわりを反時計回りに何周するか（巻き数）に他なりません。

ナイキスト判別式への適用：\(f(s) = 1 + L(s)\) とし、\(C\) を右半平面全体を囲むナイキストDコンター（虚軸 \(-jR \to +jR\) を、\(L(s)\) が虚軸上に極を持つ場合はそこだけ小さな半円で迂回しながら通り、半径 \(R \to \infty\) の右半円で閉じる経路）に取ります。\(1+L(s)\) の \(C\) 内部の零点は特性方程式 \(1+L(s)=0\) の根、すなわち右半平面にある閉ループ極（\(Z\) 個）に一致し、\(1+L(s)\) の極は \(L(s)\) の極と同じなので、右半平面にある開ループ極（\(P\) 個）に一致します。教科書的なナイキスト経路は虚軸を下から上に進んでから右半円で閉じる向き——すなわち右半平面を右手に見ながら進む、複素解析の標準（反時計回り）とは逆向き——に走らせるため、偏角の原理の符号が反転し、

\[ Z - P = -\frac{1}{2\pi}\Delta_C \arg\bigl(1+L(s)\bigr) \]

となります。右辺は「\(1+L(s)\) が原点を反時計回りに回る回数の符号を反転させたもの」であり、これはちょうど「\(L(s)\) が点 \(-1\) を時計回りに回る回数」に等しくなります。これを \(N\) と定義すれば、

\[ Z = N + P \]

が導かれます。\(N\) が符号込みで数えられる回転数であり、\(Z, P\) は個数（非負整数）であるという非対称性は、この符号反転に由来します。

エッジケース：開ループ極が虚軸上にある場合の迂回

本記事で扱ってきた \(L(s) = K/(s(s+1)(s+2))\) は、実は \(s=0\) に極を持ちます（積分器 1 個）。これは虚軸上（右半平面の境界）にあるため、上記のナイキストDコンターがそのまま \(s=0\) を通ってしまい、\(L(s)\) が定義できません。この種の開ループ極が虚軸上にある系（積分器を含む系はすべてこれに該当）では、コンターに半径 \(\epsilon \to 0\) の小さな半円を追加し、極を迂回させる必要があります。

慣習として、この小半円は右半平面側に膨らませます（\(s = \epsilon e^{i\theta}\) 、\(\theta: -90^\circ \to +90^\circ\) ）。これにより \(s=0\) の極はコンターの内部（右半平面）に含まれない——つまり「不安定極として数えない」——という取り扱いになり、\(P\) の勘定と整合します（安定な積分制御はこの意味で「境界上」を安定側として扱う、というのが標準的な工学的規約です）。

\(s=0\) 近傍で \(L(s) \approx K/(2s)\) （\((s+1)(s+2) \approx 2\) より）なので、小半円上では

\[ L(\epsilon e^{i\theta}) \approx \frac{K}{2\epsilon} e^{-i\theta} \]

となり、\(\theta\) が \(-90^\circ \to +90^\circ\) と動く間、像は半径 \(K/(2\epsilon)\) （\(\epsilon \to 0\) で無限大に発散）の円弧を、角度 \(+90^\circ \to -90^\circ\) と時計回りに \(180^\circ\) だけ掃引します。これが、通常のナイキスト線図で \(\omega \to 0^+\) 側と \(\omega \to 0^-\) 側の軌跡の間に見える「隙間」を無限遠で閉じる弧です。この弧が \(-1\) を回り込まない限り、\(N\) の勘定に影響しません。

ナイキストDコンターの原点迂回とその像

上図 (a) は \(s\) 平面上のDコンター（原点の極を右に迂回する模式図）、(b) はその \(L(s)\) による像です。正負の周波数側の軌跡（青の実線・破線）の間を、原点迂回の像である大きな弧（赤）が無限遠側で滑らかに接続しており、コンター全体が閉曲線になっていることが分かります。

数値実験：偏角の原理と迂回処理の検証

偏角の原理の主張——「\(1+L(s)\) が \(C\) を一周する像の巻き数が、\(C\) 内部の零点数から極数を引いた値に一致する」——を、実際に上記の迂回コンターに沿って \(L(s)\) を数値評価し、巻き数を計算することで検証します。比較対象として、根軌跡の節で使った特性方程式の根を直接 numpy.roots で数えた \(Z\) （右半平面の閉ループ極数）を使います。

import numpy as np

den0 = np.polymul(np.polymul([1, 0], [1, 1]), [1, 2])  # s(s+1)(s+2)

def charpoly(K, num=np.array([1.0]), den=den0.astype(float)):
    cp = den.copy()
    cp[-len(num):] += K * num
    return cp

def L_of_s(s, K):
    return K / (s * (s + 1) * (s + 2))

def nyquist_D_contour(R=1e4, eps=1e-3, n=200000):
    # -jR→-jeps(虚軸) → 原点迂回(小半円、右へ) → +jeps→+jR(虚軸) → 大半円(右半平面)
    pts = []
    pts.append(1j * np.linspace(-R, -eps, n))
    theta = np.linspace(-np.pi / 2, np.pi / 2, n)
    pts.append(eps * np.exp(1j * theta))
    pts.append(1j * np.linspace(eps, R, n))
    theta2 = np.linspace(np.pi / 2, -np.pi / 2, n)
    pts.append(R * np.exp(1j * theta2))
    return np.concatenate(pts)

def winding_number_raw(vals):
    # 反時計回りを正として原点まわりの巻き数を計算
    angles = np.unwrap(np.angle(vals))
    return (angles[-1] - angles[0]) / (2 * np.pi)

def count_RHP_roots(K):
    r = np.roots(charpoly(K))
    return int(np.sum(r.real > 1e-9))

s_contour = nyquist_D_contour()
print(f"{'K':>6} | {'N (時計回り)':>12} | {'Z (直接計算)':>12} | {'P':>3} | {'N+P':>6}")
for K in [1, 4, 6, 8, 12]:
    ret_diff = 1 + L_of_s(s_contour, K)
    N = -winding_number_raw(ret_diff)  # コンターの向き(教科書的向き)に合わせ符号反転
    Z = count_RHP_roots(K)
    print(f"{K:6.2f} | {N:12.4f} | {Z:12d} | {0:3d} | {N:6.4f}")

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

     K | N (時計回り) | Z (直接計算) |   P |    N+P
  1.00 |      -0.0000 |            0 |   0 | -0.0000
  4.00 |      -0.0000 |            0 |   0 | -0.0000
  6.00 |      -0.0000 |            0 |   0 | -0.0000
  8.00 |       2.0000 |            2 |   0 | 2.0000
 12.00 |       2.0000 |            2 |   0 | 2.0000

\(K \le 6\) （安定・限界安定）では巻き数 \(N=0\) で \(Z=0\) と一致し、\(K=8, 12\) （不安定、極が2個右半平面へ移動）では \(N=2\) で \(Z=2\) と完全に一致しています。迂回処理を正しく実装したナイキストDコンターに沿った偏角の原理の数値積分が、根から直接数えた不安定極数と厳密に一致することが確認できました。\(K=6\) ちょうど（マージナル安定、閉ループ極が虚軸上）でも \(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]))

むだ時間系の安定余裕

遅れがゲイン交差周波数を変えないこと

輸送遅れ・通信遅延・センサのサンプリング遅延などは、伝達関数にむだ時間 \(e^{-sT}\) （\(T>0\) ）を乗じる形でモデル化されます。

\[ L_T(s) = L(s)\,e^{-sT} \]

\(s = j\omega\) 上では \(|e^{-j\omega T}| = 1\) なので、

\[ |L_T(j\omega)| = |L(j\omega)|, \qquad \angle L_T(j\omega) = \angle L(j\omega) - \omega T \]

——むだ時間はゲインを一切変えず、位相だけを周波数に比例して遅らせます。したがって \(|L(j\omega)|=1\) となるゲイン交差周波数 \(\omega_{gc}\) は、\(T\) の値によらず不変です。これは近似ではなく厳密な性質で、位相余裕へのむだ時間の影響を厳密に計算できる理由になっています。

許容むだ時間の導出

ゲイン交差周波数 \(\omega_{gc}\) での位相余裕は、むだ時間がないとき \(\mathrm{PM} = 180^\circ + \angle L(j\omega_{gc})\) でした。むだ時間を加えると、同じ \(\omega_{gc}\) において

\[ \mathrm{PM}'(T) = \mathrm{PM} - \omega_{gc} T \cdot \frac{180^\circ}{\pi} \]

だけ位相余裕が目減りします（\(\mathrm{PM}\) は度、\(\omega_{gc}T\) はラジアン）。位相余裕がちょうどゼロになる——ゲイン交差周波数において位相がちょうど \(-180^\circ\) に達し、安定限界に達する——むだ時間の上限は、\(\mathrm{PM}'(T_{max})=0\) を解いて

\[ T_{max} = \frac{\mathrm{PM}[\mathrm{rad}]}{\omega_{gc}} \]

と求まります。\(\omega_{gc}\) が \(T\) に依存しないため、この式は近似式ではなく厳密な安定限界を与えます（ただし、むだ時間の増加によってゲイン交差周波数より低い別の周波数で新たなゲイン交差が生じない、という条件のもとで成り立ちます。今回の題材ではこれを数値的に確認します）。

数値実験：許容むだ時間と安定限界の比較検証

まず、既存のゲイン交差周波数・位相余裕から \(T_{max}\) を予測し、むだ時間込みの位相余裕が実際にその点でゼロになることを確認します。次に、\(e^{-sT}\) を独立した数値的手法（6次パデ近似による有理関数近似）で置き換えた特性方程式の根を掃引し、根が虚軸を横切る臨界むだ時間を二分法で求めて、解析的な \(T_{max}\) と比較します。

import numpy as np
from scipy import signal
import math
from math import comb, factorial

num0 = [1.0]
den0 = np.polymul(np.polymul([1, 0], [1, 1]), [1, 2])
sys0 = signal.TransferFunction(num0, den0)

w = np.logspace(-2, 2, 2000)
w, H = signal.freqresp(sys0, w)
mag, phase = np.abs(H), np.unwrap(np.angle(H))
idx_gc = np.argmin(np.abs(mag - 1.0))
wgc = w[idx_gc]
PM_deg = 180 + np.degrees(phase[idx_gc])
PM_rad = np.radians(PM_deg)
T_max = PM_rad / wgc
print(f"wgc={wgc:.4f} rad/s, PM={PM_deg:.2f} deg, T_max=PM[rad]/wgc={T_max:.4f} s")

for T in [0.5, 1.0, T_max, 2.0, 3.0]:
    pm_with_delay = 180 + np.degrees(phase[idx_gc] - wgc * T)
    print(f"  T={T:.4f}s: 遅れ込み位相余裕 PM'={pm_with_delay:+.3f} deg")

# --- 独立検証: N次パデ近似による特性方程式の根で臨界むだ時間を求める ---
def D_coeffs(N):
    coeffs = np.zeros(N + 1)
    for k in range(N + 1):
        coeffs[k] = factorial(2 * N - k) / factorial(2 * N) * comb(N, k)
    return coeffs  # 昇べき (x^0..x^N)

def pade_delay_num_den(T, N=6):
    d = D_coeffs(N)
    num_asc = d * np.array([(-1) ** k for k in range(N + 1)]) * np.array([T ** k for k in range(N + 1)])
    den_asc = d * np.array([T ** k for k in range(N + 1)])
    return num_asc[::-1], den_asc[::-1]  # 降べきに反転

def max_real_part_with_delay(T, N=6):
    pn, pd = pade_delay_num_den(T, N)
    charpoly = np.polyadd(np.polymul(den0.astype(float), pd), np.polymul(np.array(num0), pn))
    return np.roots(charpoly).real.max()

lo, hi = 0.5, 3.5
for _ in range(60):
    mid = (lo + hi) / 2
    if max_real_part_with_delay(lo) * max_real_part_with_delay(mid) <= 0:
        hi = mid
    else:
        lo = mid
T_crit_numeric = (lo + hi) / 2
print(f"Pade(6次)近似による臨界むだ時間: T_crit = {T_crit_numeric:.4f} s")
print(f"解析的予測 T_max = {T_max:.4f} s, 相対誤差 = {abs(T_crit_numeric-T_max)/T_max*100:.4f} %")
print(f"T=0.98*T_crit: max Re = {max_real_part_with_delay(T_crit_numeric*0.98):+.4f}")
print(f"T=1.02*T_crit: max Re = {max_real_part_with_delay(T_crit_numeric*1.02):+.4f}")

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

wgc=0.4455 rad/s, PM=53.43 deg, T_max=PM[rad]/wgc=2.0934 s
  T=0.5000s: 遅れ込み位相余裕 PM'=+40.669 deg
  T=1.0000s: 遅れ込み位相余裕 PM'=+27.907 deg
  T=2.0934s: 遅れ込み位相余裕 PM'=+0.000 deg
  T=2.0000s: 遅れ込み位相余裕 PM'=+2.383 deg
  T=3.0000s: 遅れ込み位相余裕 PM'=-23.141 deg
Pade(6次)近似による臨界むだ時間: T_crit = 2.0913 s
解析的予測 T_max = 2.0934 s, 相対誤差 = 0.0992 %
T=0.98*T_crit: max Re = -0.0027
T=1.02*T_crit: max Re = +0.0026

位相余裕から逆算した許容むだ時間 \(T_{max} = 2.0934\,\mathrm{s}\) に対し、パデ近似で構成した特性方程式の根を独立に追跡して求めた臨界むだ時間は \(T_{crit} = 2.0913\,\mathrm{s}\) で、相対誤差はわずか \(0.0992\%\) でした。\(T\) が \(T_{crit}\) をわずかに下回れば最大実部は負（安定）、わずかに上回れば正（不安定）に転じることも確認でき、「位相余裕からむだ時間の許容量を逆算する」という設計則が、根の位置に基づく安定判定と定量的に一致することが検証できました。むだ時間 \(2.09\,\mathrm{s}\) は、この系の時定数（最も遅い極が \(s=-1\) 付近）と比べて決して無視できる量ではなく、通信遅延やセンサ処理遅延が大きい制御対象では、位相余裕に十分な余裕を見込んでおく必要があることが分かります。

設計指針

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

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 つは同じ安定性を異なる視点で可視化した相補的な道具です。

MIMO系への拡張：一般化ナイキスト判別法

ここまでは入力・出力とも1つのSISO（Single-Input Single-Output）系を扱ってきましたが、開ループが \(n \times n\) の伝達関数行列 \(L(s)\) になる多入力多出力（MIMO）系でも、偏角の原理に基づく安定判別は一般化できます。

SISO の場合と異なり、MIMO では単純にスカラー \(L(j\omega)\) を複素平面に描くことはできません。その代わりに、各周波数 \(\omega\) で行列 \(L(j\omega)\) の固有値 \(\lambda_1(j\omega), \dots, \lambda_n(j\omega)\) を計算し、それぞれを \(\omega\) について複素平面上に軌跡として描いたものを特性軌跡（characteristic loci）と呼びます。MacFarlane と Postlethwaite（1977）が確立した一般化ナイキスト判別法（generalized Nyquist criterion）は、\(Z=N+P\) と同じ形の判別式が特性軌跡についても成り立つことを示します——ここで \(N\) は「\(n\) 本の特性軌跡が臨界点 \(-1\) を時計回りに囲む回数の合計」です。

証明の考え方はSISOと同じ偏角の原理ですが、固有値 \(\lambda_i(s)\) は一般に \(s\) の代数関数であり、固有値どうしが交差する分岐点（branch point）を持つため、単一の複素平面ではなくリーマン面上で偏角の原理を適用する必要があります。この技術的な複雑さのため、各対角要素だけを個別にSISO的にナイキスト判別するのは（ループ間の干渉を無視することになり）一般には誤りで、特性軌跡による全体としての判別が必要になります。本記事のスコープを超えるため深入りしませんが、「MIMO系の安定判別は個々のループの安定余裕の単純な組み合わせでは代用できない」という点は、多変数制御を扱う際の重要な注意点です。

離散時間系での注意

ディジタル制御では、\(s\) 平面の左半平面が \(z\) 平面の単位円内に対応します。安定判別の臨界点や虚軸の役割が単位円に置き換わる点に注意が必要です（ Z変換の定義から安定性判定まで で扱った通り、全極が単位円内 \(|z_k|<1\) であることが離散系の安定条件でした）。

連続系設計を 双一次変換（Tustin変換） \(s = \frac{2}{T_s}\cdot\frac{1-z^{-1}}{1+z^{-1}}\) でディジタル化する場合、実周波数 \(\omega_a\) （アナログ）と \(z=e^{j\omega_d T_s}\) に対応する \(\omega_d\) （ディジタル）の関係は、\(s=j\omega_a\) を代入して解くと

\[ \omega_d = \frac{2}{T_s}\arctan\!\left(\frac{\omega_a T_s}{2}\right) \]

という非線形な対応になります。これを**周波数ワーピング（frequency warping）**と呼びます。\(\omega_a T_s \ll 1\) （サンプリング周波数に対して十分低い周波数）では \(\arctan x \approx x\) なので \(\omega_d \approx \omega_a\) とほぼ線形に対応しますが、ナイキスト周波数 \(\pi/T_s\) に近づくにつれて \(\arctan\) が飽和し、対応が大きく歪みます。ゲイン交差周波数・位相交差周波数がサンプリング周波数に対して十分低ければ安定余裕の値はほぼそのまま引き継がれますが、高速な応答を狙って交差周波数がナイキスト周波数に近づくと、双一次変換によって位相余裕が設計値より目減りすることがあるため、 標本化定理 に基づくサンプリング周波数の選定と合わせて注意が必要です。

近年の研究動向：データ駆動制御・学習ベース制御との接続

ナイキスト線図に代表される古典的な周波数領域解析は100年近く前の枠組みですが、モデルを陽に用いないデータ駆動制御・学習ベース制御が発展する中でも、その考え方は形を変えて生き続けています。

  • 円判別法（Circle Criterion）・Popov判別法は、ナイキストの安定判別法をセクタ条件を満たす非線形要素を含む系に一般化したものです。Richardson, Turner, and Gunn (2023, IEEE Control Systems Letters) は、フィードバックループにReLU活性化関数を持つニューラルネットワークを含む系（Lur’e系）に対して、繰り返しReLU特有の性質を利用した新しい二次制約に基づき円判別法・Popov判別法を強化し、保守性を抑えた安定判別条件を導出しています。ナイキスト平面上の幾何学的な判別条件が、学習済みニューラルネットワークをコントローラやプラントモデルの一部に含む系の検証にも応用されている例です。
  • Umar, Ryu, Back, and Kim (2024, Mathematics) は、観測データのみからLQR型コントローラを設計する「データ駆動LQR」に対して、入力チャネルの不確かさに対する安定余裕（ディスクマージン）を定量化し、モデルベース設計における古典的な安定余裕の概念が、データ駆動設計でも有効な頑健性の指標として引き継がれることを示しています。
  • Berberich et al. (2023, American Control Conference) など、データ駆動制御の理論的基盤である Willems の基本補題まわりの構成的な証明や頑健性拡張も活発に進められており、モデルベースの周波数領域解析とデータ駆動アプローチの橋渡しは現在進行形の研究領域です。

これらはいずれも「古典的な安定余裕・周波数領域判別の考え方そのものが不要になった」のではなく、「学習・データ駆動という新しい設計手法に対しても、同種の頑健性の物差しを再定義し直す」方向の研究である点が共通しています。

まとめ

  • 安定余裕（ゲイン余裕・位相余裕）は、開ループ周波数応答 \(L(j\omega)\) から閉ループ安定性の「余裕」を定量化する。位相余裕 \(45^\circ \sim 60^\circ\) が実用目安。
  • ナイキスト線図は \(L(j\omega)\) を複素平面に描き、臨界点 \(-1\) まわりの回転数で安定判別する（\(Z = N + P\) ）。この式は複素解析の偏角の原理から厳密に導出でき、実際にナイキストDコンターに沿った周回積分（巻き数）を数値計算し、根から直接数えた不安定極数と厳密に一致することを確認した。
  • 開ループ極が虚軸上にある（積分器を含む）系では、ナイキストDコンターに小さな半円の迂回が必要になる。今回の題材 \(L(s)=K/(s(s+1)(s+2))\) で実際にこの迂回を実装し、迂回の像が無限遠を回る弧として現れることを可視化・数値検証した。
  • 根軌跡はゲイン \(K\) を掃引したときの閉ループ極の軌跡で、極が虚軸を横切るゲインが安定限界。
  • 題材 \(L(s) = K/(s(s+1)(s+2))\) では、位相余裕 \(53.4^\circ\) ・ゲイン余裕 \(15.6\,\mathrm{dB}\) ・安定限界 \(K = 6\) の 3 者が完全に整合することを scipy.signal と numpy.roots で確認した。
  • むだ時間 \(e^{-sT}\) はゲイン交差周波数を変えず位相のみを遅らせるため、許容むだ時間 \(T_{max}=\mathrm{PM}[\mathrm{rad}]/\omega_{gc}\) が厳密に求まる。今回の題材では \(T_{max}=2.093\,\mathrm{s}\) で、パデ近似で構成した特性方程式の根から求めた臨界むだ時間 \(2.091\,\mathrm{s}\) （相対誤差 \(0.099\%\) ）とよく一致した。
  • MIMO系では特性軌跡（characteristic loci）を用いた一般化ナイキスト判別法、離散系では単位円と周波数ワーピングへの考慮が必要になる。
  • control パッケージがなくても、周波数応答（freqresp）と多項式の根（roots）だけで安定余裕・ナイキスト・根軌跡・むだ時間・迂回処理はすべて実装できる。

関連記事

参考文献

  • Ogata, K. (2010). Modern Control Engineering (5th ed.). Prentice Hall.
  • Åström, K. J., & Murray, R. M. (2008). Feedback Systems: An Introduction for Scientists and Engineers. Princeton University Press.
  • Nyquist, H. (1932). “Regeneration Theory.” Bell System Technical Journal, 11(1), 126–147.
  • Ahlfors, L. V. (1979). Complex Analysis (3rd ed.). McGraw-Hill.（偏角の原理の標準的な証明）
  • MacFarlane, A. G. J., & Postlethwaite, I. (1977). “The generalized Nyquist stability criterion and multivariable root loci.” International Journal of Control, 25(1), 81–127.
  • Richardson, C., Turner, M., & Gunn, S. (2023). “Strengthened Circle and Popov Criteria for the Stability Analysis of Feedback Systems With ReLU Neural Networks.” IEEE Control Systems Letters, 7, 2635–2640.
  • Umar, A. A., Ryu, K., Back, J., & Kim, J.-S. (2024). “Stability Margin of Data-Driven LQR and Its Application to Consensus Problem.” Mathematics, 12(2), 199.
  • scipy.signal.freqresp — SciPy documentation
  • numpy.roots — NumPy documentation

関連ツール

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

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

関連記事