ラウス・フルヴィッツの安定判別法は、特性方程式の係数から直接安定性を判定できる代数的な方法です。特性方程式の根(極)を実際に求めることなく、すべての根が左半面にあるかどうかを判定できます。
本記事では、ラウス表の構成方法から安定判別の手順、特殊ケースの処理、Pythonでの実装までを解説します。
本記事の内容
- 安定性の代数的条件
- ラウス表の作り方と安定判別の手順
- 特殊ケース(ゼロ要素、ゼロ行)の処理
- Pythonによる実装と検証
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
フルヴィッツの安定条件
特性方程式が次のように与えられたとします。
$$ D(s) = a_n s^n + a_{n-1} s^{n-1} + \cdots + a_1 s + a_0 = 0 $$
フルヴィッツの必要条件: すべての根の実部が負であるためには、すべての係数 $a_i$ が同符号(通常正)でなければなりません。
ただし、これは必要条件であり十分条件ではありません。十分条件を得るためにラウスの方法を用います。
ラウス表の構成
特性方程式 $a_n s^n + a_{n-1} s^{n-1} + \cdots + a_0 = 0$ に対して、次のようにラウス表を構成します。
$$ \begin{array}{c|ccc} s^n & a_n & a_{n-2} & a_{n-4} & \cdots \\ s^{n-1} & a_{n-1} & a_{n-3} & a_{n-5} & \cdots \\ s^{n-2} & b_1 & b_2 & b_3 & \cdots \\ s^{n-3} & c_1 & c_2 & c_3 & \cdots \\ \vdots & \vdots & & & \\ s^0 & & & & \end{array} $$
ここで各要素は次のように計算します。
$$ b_1 = \frac{a_{n-1} \cdot a_{n-2} – a_n \cdot a_{n-3}}{a_{n-1}} $$
$$ b_2 = \frac{a_{n-1} \cdot a_{n-4} – a_n \cdot a_{n-5}}{a_{n-1}} $$
一般に、
$$ b_k = \frac{a_{n-1} \cdot a_{n-2k} – a_n \cdot a_{n-2k-1}}{a_{n-1}} $$
同様に $c$ 行は $a$ の行と $b$ の行から計算します。
$$ c_1 = \frac{b_1 \cdot a_{n-3} – a_{n-1} \cdot b_2}{b_1} $$
ラウスの安定判別定理
定理: 特性方程式のすべての根が左半面にある(システムが安定である)ための必要十分条件は、ラウス表の第1列のすべての要素が正であることです。
さらに、第1列の符号変化の回数が、右半面にある根の数に等しくなります。
具体例
例1: 安定なシステム
$$ D(s) = s^3 + 6s^2 + 11s + 6 = 0 $$
ラウス表を構成します。
$$ \begin{array}{c|cc} s^3 & 1 & 11 \\ s^2 & 6 & 6 \\ s^1 & b_1 & \\ s^0 & c_1 & \end{array} $$
$$ b_1 = \frac{6 \cdot 11 – 1 \cdot 6}{6} = \frac{66 – 6}{6} = 10 $$
$$ c_1 = \frac{10 \cdot 6 – 6 \cdot 0}{10} = 6 $$
第1列: $1, 6, 10, 6$ — すべて正なので安定です。
実際の根は $s = -1, -2, -3$ で、すべて左半面にあります。
例2: 不安定なシステム
$$ D(s) = s^4 + 2s^3 + 3s^2 + 4s + 5 = 0 $$
$$ \begin{array}{c|cc} s^4 & 1 & 3 & 5 \\ s^3 & 2 & 4 & \\ s^2 & b_1 & b_2 & \\ s^1 & c_1 & & \\ s^0 & d_1 & & \end{array} $$
$$ b_1 = \frac{2 \cdot 3 – 1 \cdot 4}{2} = 1, \quad b_2 = \frac{2 \cdot 5 – 1 \cdot 0}{2} = 5 $$
$$ c_1 = \frac{1 \cdot 4 – 2 \cdot 5}{1} = -6 $$
$$ d_1 = \frac{-6 \cdot 5 – 1 \cdot 0}{-6} = 5 $$
第1列: $1, 2, 1, -6, 5$ — 符号変化が2回($1 \to -6$ と $-6 \to 5$)なので、右半面に2つの根があり不安定です。
特殊ケース
ケース1: 第1列にゼロが現れる場合
第1列の要素がゼロになると、次の行の計算でゼロ除算が発生します。この場合、ゼロを微小な正の数 $\epsilon > 0$ で置き換えて計算を続け、最後に $\epsilon \to 0$ の極限を取ります。
ケース2: 行全体がゼロになる場合
これはシステムに虚軸上の根(純虚根)の対称なペアが存在することを意味します。ゼロ行の一つ上の行を用いて補助多項式を構成し、その微分で置き換えて計算を続けます。
Pythonでの実装
ラウス表の生成
import numpy as np
def routh_table(coeffs):
"""ラウス表を生成する"""
n = len(coeffs)
# ラウス表の行数はn、列数はceil(n/2)
rows = n
cols = (n + 1) // 2
table = np.zeros((rows, cols))
# 最初の2行を埋める
table[0, :len(coeffs[0::2])] = coeffs[0::2]
table[1, :len(coeffs[1::2])] = coeffs[1::2]
# 残りの行を計算
for i in range(2, rows):
for j in range(cols - 1):
if abs(table[i-1, 0]) < 1e-10:
# 第1列がゼロの場合、微小値で置換
table[i-1, 0] = 1e-10
table[i, j] = (table[i-1, 0] * table[i-2, j+1]
- table[i-2, 0] * table[i-1, j+1]) / table[i-1, 0]
return table
def check_stability(coeffs):
"""ラウス・フルヴィッツ安定判別"""
table = routh_table(coeffs)
first_col = table[:, 0]
# 符号変化の回数を数える
sign_changes = 0
for i in range(1, len(first_col)):
if first_col[i-1] * first_col[i] < 0:
sign_changes += 1
print("ラウス表:")
print(table)
print(f"\n第1列: {first_col}")
print(f"符号変化回数: {sign_changes}")
if sign_changes == 0:
print("判定: 安定(すべての極が左半面)")
else:
print(f"判定: 不安定(右半面に{sign_changes}個の極)")
return sign_changes == 0
# 例1: 安定なシステム
print("=== 例1: s^3 + 6s^2 + 11s + 6 ===")
stable1 = check_stability([1, 6, 11, 6])
print("\n=== 例2: s^4 + 2s^3 + 3s^2 + 4s + 5 ===")
stable2 = check_stability([1, 2, 3, 4, 5])
# numpyのrootsで検証
print("\n=== 検証(根の計算) ===")
roots1 = np.roots([1, 6, 11, 6])
print(f"例1の根: {roots1}")
roots2 = np.roots([1, 2, 3, 4, 5])
print(f"例2の根: {roots2}")
print(f"例2の右半面の根の数: {np.sum(np.real(roots2) > 0)}")
フィードバック系の安定範囲の解析
import numpy as np
import matplotlib.pyplot as plt
def count_unstable_poles(coeffs):
"""不安定な極の数を返す"""
n = len(coeffs)
cols = (n + 1) // 2
table = np.zeros((n, cols))
table[0, :len(coeffs[0::2])] = coeffs[0::2]
table[1, :len(coeffs[1::2])] = coeffs[1::2]
for i in range(2, n):
for j in range(cols - 1):
if abs(table[i-1, 0]) < 1e-10:
table[i-1, 0] = 1e-10
table[i, j] = (table[i-1, 0] * table[i-2, j+1]
- table[i-2, 0] * table[i-1, j+1]) / table[i-1, 0]
first_col = table[:, 0]
sign_changes = sum(1 for i in range(1, len(first_col))
if first_col[i-1] * first_col[i] < 0)
return sign_changes
# プラント: P(s) = 1/(s(s+1)(s+2)), コントローラ: C(s) = K
# 特性方程式: s^3 + 3s^2 + 2s + K = 0
K_values = np.linspace(0.01, 10, 500)
unstable_counts = []
for K in K_values:
coeffs = [1, 3, 2, K]
unstable_counts.append(count_unstable_poles(coeffs))
# 安定範囲の可視化
plt.figure(figsize=(10, 5))
plt.plot(K_values, unstable_counts, 'b-', linewidth=2)
plt.xlabel('Gain K')
plt.ylabel('Number of unstable poles')
plt.title('Routh-Hurwitz Stability Analysis: $s^3 + 3s^2 + 2s + K = 0$')
plt.grid(True, alpha=0.3)
plt.axhline(y=0, color='g', linestyle='--', alpha=0.5, label='Stable region')
plt.legend()
plt.tight_layout()
plt.show()
# 臨界ゲインを求める(ラウスの条件: 3*2 - K > 0 → K < 6)
print(f"安定条件: 0 < K < 6")
print(f"臨界ゲイン: K_cr = 6")
ラウス表の $s^1$ 行の第1列は $\frac{3 \cdot 2 – 1 \cdot K}{3} = \frac{6 – K}{3}$ となるため、安定条件は $K < 6$ です。
まとめ
本記事では、ラウス・フルヴィッツの安定判別法について解説しました。
- ラウス表は特性方程式の係数から機械的に構成できる
- 第1列のすべての要素が正であれば、システムは安定
- 第1列の符号変化回数が、右半面にある極の数に等しい
- 特殊ケース(ゼロ要素、ゼロ行)にも対処法がある
- パラメータ $K$ に関する安定範囲を代数的に求められる
次のステップとして、以下の記事も参考にしてください。