H∞制御

H∞制御(H-infinity control)の基礎を解説。H∞ノルムの定義、外乱抑制のための制御器設計、ミニマックス問題としての定式化を紹介します。

H∞制御(H-infinity control)は、ロバスト制御理論の一つであり、システムに加わる外乱モデルの不確かさに対して、制御性能を保証することを目的とした制御器設計手法です。特に、最悪ケースの外乱に対するシステムの応答を最小化することを目指します。

古典制御( PID制御ゲイン余裕・位相余裕による安定余裕の評価 )は、あらかじめ決められた1つのプラントモデルに対して制御器を設計し、事後的に安定余裕でロバスト性を確認します。これに対しH∞制御は、プラントの不確かさそのものを設計変数に組み込み、想定される不確かさの集合すべてに対して性能を保証する制御器を、最適化問題として直接求める点が本質的に異なります。本記事では、その理論的な柱である小ゲイン定理H∞ノルム標準H∞制御問題リカッチ方程式による制御器の存在条件を、可能な限り導出を追いながら解説し、最後に python-control(SLICOTベースの hinfsyn)を使った数値実験で、PID制御との頑健性の違いを実測します。

H∞ノルム

定義

外乱 \(w\) から評価出力 \(z\) への閉ループ伝達関数を \(T_{zw}(s)\) とすると、そのH∞ノルムは周波数応答の最大特異値の上限として定義されます。

\[ \|T_{zw}\|_\infty = \sup_{\omega \in \mathbb{R}} \bar\sigma\big(T_{zw}(j\omega)\big) \]

ここで \(\bar\sigma(\cdot)\) は行列の最大特異値です。SISO(1入力1出力)系では \(\bar\sigma(T_{zw}(j\omega)) = |T_{zw}(j\omega)|\) となるので、H∞ノルムは単純にボード線図のゲイン曲線のピーク値(最悪周波数でのゲイン)に一致します。

なぜ「最悪ケースのゲイン」なのか:誘導L2ノルムとしての意味

H∞ノルムのもう一つの顔は、信号のエネルギー比の上限です。

\[ \|z\|_2^2 = \int_0^\infty z^T(t) z(t)\, dt, \qquad \|w\|_2^2 = \int_0^\infty w^T(t) w(t)\, dt \]

とL2ノルム(信号エネルギー)を定義すると、安定なシステム \(T_{zw}\) について次の等式が成り立ちます。

\[ \|T_{zw}\|_\infty = \sup_{w \in L_2,\, w \neq 0} \frac{\|z\|_2}{\|w\|_2} \]

これは「H∞ノルムは有限エネルギー入力に対する誘導ノルム(induced norm)である」という主張で、天下り的に見えるので証明の要点を追っておきます。

上界の証明(\(\leq\) ): パーセバルの定理より、時間領域のエネルギーは周波数領域のエネルギーに等しく保存されます。

\[ \|z\|_2^2 = \frac{1}{2\pi}\int_{-\infty}^{\infty} \|\hat Z(j\omega)\|^2 d\omega = \frac{1}{2\pi}\int_{-\infty}^{\infty} \|T_{zw}(j\omega)\hat W(j\omega)\|^2 d\omega \]

最大特異値の定義より、任意のベクトル \(v\) に対し \(\|T_{zw}(j\omega) v\| \leq \bar\sigma(T_{zw}(j\omega))\|v\|\) が成り立つので、

\[ \|z\|_2^2 \leq \frac{1}{2\pi}\int_{-\infty}^{\infty} \bar\sigma(T_{zw}(j\omega))^2 \|\hat W(j\omega)\|^2 d\omega \leq \|T_{zw}\|_\infty^2 \cdot \frac{1}{2\pi}\int_{-\infty}^{\infty} \|\hat W(j\omega)\|^2 d\omega = \|T_{zw}\|_\infty^2 \|w\|_2^2 \]

したがって \(\|z\|_2/\|w\|_2 \leq \|T_{zw}\|_\infty\) が任意の \(w\) について成り立ちます。

下界の証明(\(\geq\) 、達成可能性): ピーク周波数を \(\omega_0 = \arg\sup_\omega \bar\sigma(T_{zw}(j\omega))\) とし、\(T_{zw}(j\omega_0)\) の最大特異値に対応する右特異ベクトル方向 \(v_0\) の入力を考えます。\(w(t)\) として、周波数 \(\omega_0\) を中心とする帯域幅 \(\varepsilon\) の狭帯域信号(方向は \(v_0\) )を作ると、\(T_{zw}\) が連続なので帯域内では \(\bar\sigma(T_{zw}(j\omega)) \approx \bar\sigma(T_{zw}(j\omega_0))\) とみなせ、出力エネルギーと入力エネルギーの比は \(\varepsilon \to 0\) の極限で \(\bar\sigma(T_{zw}(j\omega_0)) = \|T_{zw}\|_\infty\) に収束します。

以上より上界と下界が一致し、\(\sup_w \|z\|_2/\|w\|_2 = \|T_{zw}\|_\infty\) が成立します。つまりH∞ノルムは、あらゆる有限エネルギー外乱の中で最も出力を増幅させる方向・周波数を選んだときのゲインであり、これが「最悪ケースのゲイン」と呼ばれる所以です。

数値検証

この等式を実際に確認します。軽く減衰した2次系 \(G(s) = \dfrac{1}{s^2+0.4s+4}\) (固有振動数 \(\omega_n=2\,\mathrm{rad/s}\) 、減衰比 \(\zeta=0.1\) )を例に、H∞ノルムを2通りの方法で計算しました。

import numpy as np
import control as ct

G = ct.tf([1.0], [1, 0.4, 4])
norm_inf, wpeak = ct.linfnorm(G)
print(f"control.linfnorm: ||G||inf = {norm_inf:.6f} at w = {wpeak:.4f} rad/s")

w = np.logspace(-2, 3, 200000)
mag = np.abs(G(1j*w))
print(f"brute force sup: max|G(jw)| = {mag.max():.6f} at w = {w[np.argmax(mag)]:.4f} rad/s")

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

control.linfnorm: ||G||inf = 1.256297 at w = 1.9799 rad/s
brute force sup: max|G(jw)| = 1.256297 at w = 1.9799 rad/s

python-control の内部アルゴリズム(周波数グリッドの最大値ではなく、状態空間の固有値問題を反復して求める厳密解法)とブルートフォースな周波数走査が、小数点以下6桁まで完全に一致しました。ピーク周波数 \(\omega=1.98\,\mathrm{rad/s}\) が共振周波数 \(\omega_n=2\) にほぼ一致することも、軽減衰系の直感と整合します。

H∞制御問題の目的

H∞制御問題の目的は、このH∞ノルム \(\|T_{zw}\|_\infty\) を、ある正の定数 \(\gamma\) 以下に抑えるような制御器 \(K\) を設計することです。

\[ \|T_{zw}\|_\infty < \gamma \]

を満たす制御器を見つけることが目標であり、\(\gamma\) は性能指標と呼ばれ、小さいほど外乱抑制性能が高いことを意味します。ゲーム理論的な視点からは、以下のミニマックス問題として表現されます。

\[ \min_u \max_w \int_0^\infty \big(z^T(t)z(t) - \gamma^2 w^T(t)w(t)\big)\, dt \]

制御入力 \(u\) が \(z^Tz\) を最小化しようとする一方、最悪ケースの外乱 \(w\) は \(w^Tw\) を最大化しようとする——この \(u\) と \(w\) のゼロサムゲームの鞍点を求める問題が、後述するリカッチ方程式の導出の出発点になります。

小ゲイン定理

H∞制御がロバスト性を保証できる理論的な根拠が**小ゲイン定理(Small Gain Theorem)**です。

設定

安定な公称システム \(M(s)\) と、安定だが大きさ以外は未知の不確かさ \(\Delta(s)\) が、次のようにフィードバックループを構成しているとします。

\[ y = Mu, \qquad u = w + \Delta y \]

これは、外部入力 \(w\) に不確かさ由来の摂動 \(\Delta y\) が加算的に加わる、標準的な \(M\) –\(\Delta\) 構造です。\(u\) について解くと、

\[ u = w + \Delta M u \;\; \Longrightarrow \;\; (I - \Delta M) u = w \;\; \Longrightarrow \;\; u = (I - \Delta M)^{-1} w \]

この閉ループが内部安定であるための条件を考えます。

定理と証明

小ゲイン定理: \(M, \Delta\) がともに安定(右半平面に極を持たない)で、

\[ \|M\|_\infty \, \|\Delta\|_\infty < 1 \]

が成り立つならば、任意のそのような \(\Delta\) に対して閉ループ \((I-\Delta M)^{-1}\) は内部安定である。

証明1(縮小写像・ノイマン級数による証明): 誘導ノルムの劣乗法性(submultiplicativity)より、各周波数 \(\omega\) で

\[ \bar\sigma(\Delta(j\omega)M(j\omega)) \leq \bar\sigma(\Delta(j\omega))\,\bar\sigma(M(j\omega)) \leq \|\Delta\|_\infty \|M\|_\infty < 1 \]

が成り立ちます。よって作用素 \(\Delta M\) のスペクトル半径(各周波数での最大特異値)が一様に1未満となり、\((I-\Delta M)^{-1}\) はノイマン級数

\[ (I - \Delta M)^{-1} = \sum_{k=0}^{\infty} (\Delta M)^k \]

として収束します(等比級数の収束条件 \(\|\Delta M\|<1\) と同じ理屈です)。この級数は安定な有理関数の積・和の可算和として構成されているため、右半平面で有界・解析的であり、\((I-\Delta M)^{-1}\) 自体が安定な伝達関数になります。つまり、ループを構成する信号は常に有界となり、内部安定性が保証されます。

証明2(ナイキストの安定判別法との関係): ナイキスト線図・根軌跡の記事 で見たように、開ループ伝達関数の特性方程式 \(\det(I - L(s)) = 0\) の根の位置は、ナイキスト軌跡の臨界点まわりの回転数 \(Z=N+P\) で判定できます。ここでは \(L=\Delta M\) が開ループ伝達関数に相当します。

  • \(\Delta=0\) (不確かさがない公称ループ)のとき、\(M\) が安定なら明らかに閉ループも安定(\(Z=0\) )です。
  • \(\Delta\) を \(0\) から実際の値まで連続的にスケーリングするホモトピー \(\Delta_\tau = \tau\Delta\ (\tau \in [0,1])\) を考えると、小ゲイン条件より全ての \(\tau, \omega\) に対して \(\bar\sigma(\Delta_\tau(j\omega)M(j\omega)) \leq \tau\|\Delta\|_\infty\|M\|_\infty < 1\) なので、\(\det(I-\Delta_\tau(j\omega)M(j\omega))\) の固有値(特性利得)は複素平面上で常に半径1未満の円内にあり、臨界点 \(1+0j\) に到達も包囲もできません。
  • したがって \(\tau\) を連続的に動かす間、ナイキスト軌跡が臨界点をまたぐことはなく、回転数 \(N\) は変化しません。\(\tau=0\) で \(Z=0\) (安定)だったので、\(\tau=1\) (実際の \(\Delta\) )でも \(Z=0\) 、すなわち閉ループは安定であり続けます。

この2つの証明はそれぞれ「作用素論的な縮小写像」と「複素関数論的な連続変形(ホモトピー)」という異なる立場からの説明ですが、同じ結論——ループゲインの積が1未満なら不安定化しない——を与えます。小ゲイン定理は十分条件であり必要条件ではない(つまり\(\|M\|_\infty\|\Delta\|_\infty \geq 1\) でも安定な場合はある)ことに注意が必要です。これは、\(\Delta\) の位相情報を一切使わず大きさのみで評価する、保守的な(安全側の)判定だからです。

数値検証:小ゲイン条件のタイトさ

小ゲイン定理が「十分条件としてどれだけタイトか」を、後述の数値実験で構成する \(M(s) = W_3(s)T(s)\) (ロバスト性重み × 相補感度関数)を使って確認します。実定数ゲイン \(\delta\) による摂動 \(\Delta = \delta\) を、\(M\) を含むループ \(1 - \delta M(s) = 0\) に対してスイープし、特性方程式の根が虚軸を横切る \(\delta\) を探索しました。

w = np.logspace(-2, 3, 20000)
Mpeak = np.abs(M(1j*w)).max()
print(f"||M||_inf = {Mpeak:.5f}  ->  small-gain bound 1/||M||_inf = {1/Mpeak:.4f}")

for delta in np.linspace(0, 20, 4001):
    poles = ct.poles(ct.feedback(1, -delta*M))
    if not np.all(poles.real < 0):
        print(f"instability first occurs near delta = {delta:.3f}")
        break

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

||M||_inf = 0.13930  ->  small-gain bound 1/||M||_inf = 7.1786
instability first occurs near delta = 7.255

小ゲイン定理が保証する安定限界 \(1/\|M\|_\infty = 7.179\) に対し、実際に不安定化するのは \(\delta \approx 7.255\) でした。理論的な下限(保証値)が実際の不安定化点よりわずかに小さい(約1%の保守性)ことが確認でき、小ゲイン定理が「十分条件」として機能しつつも、この例ではほぼタイトであることが分かります。

標準H∞制御問題の定式化

一般化プラント

H∞制御では、制御対象そのものではなく、外乱・評価出力・重み関数まで含めた「一般化プラント」 \(P(s)\) を設計対象とします。\(P\) は外乱 \(w\) ・制御入力 \(u\) を入力とし、評価出力 \(z\) (性能を測る仮想信号)・観測 \(y\) (制御器に渡される信号)を出力とする、次のブロック構造を持ちます。

\[ \begin{bmatrix} z \\ y \end{bmatrix} = P(s) \begin{bmatrix} w \\ u \end{bmatrix} = \begin{bmatrix} P_{11}(s) & P_{12}(s) \\ P_{21}(s) & P_{22}(s) \end{bmatrix} \begin{bmatrix} w \\ u \end{bmatrix} \]

制御器 \(K(s)\) は観測 \(y\) から制御入力を生成します(\(u = Ky\) )。この閉ループ構造は、信号線で書けば

      +------------------+
  w ->|                  |-> z
      |     P(s)         |
  u ->|                  |-> y
      +------------------+
        ^                |
        |     K(s)       |
        +----------------+

に相当します。

下側線形分数変換(Lower LFT)の導出

\(w\) から \(z\) への閉ループ伝達関数を、\(u=Ky\) を代入して求めます。

\[ z = P_{11}w + P_{12}u, \qquad y = P_{21}w + P_{22}u, \qquad u = Ky \]

第2式に \(u=Ky\) を代入すると \(y = P_{21}w + P_{22}Ky\) となり、\(y\) について解くと

\[ (I - P_{22}K)y = P_{21}w \;\;\Longrightarrow\;\; y = (I-P_{22}K)^{-1}P_{21}w \]

これを \(u=Ky\) に代入し、さらに \(z\) の式に代入すると、

\[ z = P_{11}w + P_{12}K(I-P_{22}K)^{-1}P_{21}w \]

したがって、\(w\) から \(z\) への閉ループ伝達関数(下側線形分数変換 \(F_l(P,K)\) )は

\[ T_{zw}(s) = F_l(P,K) = P_{11} + P_{12}K(I-P_{22}K)^{-1}P_{21} \]

となります。標準H∞制御問題とは、この \(T_{zw}\) の内部安定性を保ちながら \(\|T_{zw}\|_\infty < \gamma\) を達成する \(K\) を求める問題です。

具体例:混合感度問題への当てはめ

後の数値実験で使う「混合感度(S/KS/T)」設計問題を、この一般化プラントの枠組みで具体的に書き下します。制御対象を \(G(s)\) 、性能重み \(W_1\) 、制御努力重み \(W_2\) 、ロバスト性重み \(W_3\) とし、\(w=r\) (目標値)、\(e=r-Gu\) (偏差)とすると、評価出力を

\[ z_1 = W_1 e = W_1 w - W_1 G u, \qquad z_2 = W_2 u, \qquad z_3 = W_3 G u, \qquad y = w - Gu \]

と定義すれば、

\[ P_{11} = \begin{bmatrix} W_1 \\ 0 \\ 0 \end{bmatrix}, \quad P_{12} = \begin{bmatrix} -W_1 G \\ W_2 \\ W_3 G \end{bmatrix}, \quad P_{21} = I, \quad P_{22} = -G \]

となります。感度関数 \(S=(I+GK)^{-1}\) 、相補感度関数 \(T=I-S\) を使うと、閉ループは \(z_1 = W_1 S w\) 、\(z_2 = W_2 KS w\) 、\(z_3 = W_3 T w\) と書け、\(\|T_{zw}\|_\infty<\gamma\) は3つの周波数重み付き伝達関数 \(W_1S,\ W_2KS,\ W_3T\) を同時に \(\gamma\) 未満に抑える問題に帰着します。python-controlcontrol.mixsyn(G, W1, W2, W3) は、まさにこの一般化プラントを内部で構成し hinfsyn に渡す関数です。

リカッチ方程式による存在条件

状態フィードバック(全情報)問題の導出

一般の出力フィードバックH∞制御問題の完全な導出は長くなるため、まず本質が見える状態フィードバック(full-information)問題を完全に導出し、一般の場合は結果と参考文献を示します。

プラントを

\[ \dot x = Ax + B_1 w + B_2 u, \qquad z = \begin{bmatrix} C_1 x \\ u \end{bmatrix}, \quad D_{12}^TD_{12}=I,\ C_1^TD_{12}=0 \]

とし、状態 \(x\) を直接観測できる(\(y=x\) )状況を考えます。目的は、コスト

\[ J_\gamma(u,w) = \int_0^\infty \big(x^TC_1^TC_1x + u^Tu - \gamma^2w^Tw\big)\,dt \]

の鞍点 \(\min_u\max_w J_\gamma\) を求めることです(\(J_\gamma<0\) を全ての \(w\neq0\) で達成できれば \(\|T_{zw}\|_\infty<\gamma\) が示せます)。価値関数を \(V(x)=x^TXx\) (\(X=X^T\geq0\) )と仮定し、ハミルトン・ヤコビ・アイザックス方程式

\[ \frac{d}{dt}(x^TXx) + x^TC_1^TC_1x + u^Tu - \gamma^2w^Tw = 0 \]

を \(\dot x = Ax+B_1w+B_2u\) に沿って評価すると、

\[ x^T(A^TX+XA+C_1^TC_1)x + 2w^TB_1^TXx + 2u^TB_2^TXx + u^Tu - \gamma^2w^Tw = 0 \]

\(u,w\) についてそれぞれ平方完成すると、

\[ u^Tu+2u^TB_2^TXx = (u+B_2^TXx)^T(u+B_2^TXx) - x^TXB_2B_2^TXx \] \[ -\gamma^2w^Tw+2w^TB_1^TXx = -\gamma^2(w-\gamma^{-2}B_1^TXx)^T(w-\gamma^{-2}B_1^TXx) + \gamma^{-2}x^TXB_1B_1^TXx \]

これらを代入すると、\(x\) に依存する残りの項がすべて消えるために \(X\) が満たすべき条件として、次の**H∞代数リカッチ方程式(ARE)**が導かれます。

\[ A^TX+XA+X\big(\gamma^{-2}B_1B_1^T-B_2B_2^T\big)X+C_1^TC_1=0 \]

このとき鞍点を与える制御則・最悪外乱は、平方完成の非負項をゼロにする選択

\[ u^\star = -B_2^TXx \quad (\text{最適制御}), \qquad w^\star = \gamma^{-2}B_1^TXx \quad (\text{最悪外乱}) \]

です。

ハミルトン行列の固有値条件とγ反復

このAREが安定化解 \(X\geq0\) を持つかどうかは、対応するハミルトン行列

\[ H_\gamma = \begin{bmatrix} A & -(B_2B_2^T-\gamma^{-2}B_1B_1^T) \\ -C_1^TC_1 & -A^T \end{bmatrix} \]

の固有値配置で判定できます。\(H_\gamma\) が虚軸上に固有値を持たない(\(\mathrm{dom}(\mathrm{Ric})\) に属する)ことが、安定化解 \(X=X^T\geq0\) が存在するための必要十分条件です。\(\gamma\) を大きくすると条件は緩くなり(\(\gamma\to\infty\) でH2制御に帰着)、\(\gamma\) を小さくしていくとどこかで \(H_\gamma\) の固有値が虚軸に達し、それ以下では解が存在しなくなります。この境界値 \(\gamma_{\min}\) を求めるために、\(\gamma\) を外側のループで二分探索しながら「その \(\gamma\) で安定化解が存在するか」を内側で判定する——これが**γ反復(gamma iteration)**の考え方です。実務で使われる hinfsyn(SLICOTのSB10AD)も、出力フィードバック版(\(X_\infty,Y_\infty\) の2本のリカッチ方程式とスペクトル半径条件 \(\rho(X_\infty Y_\infty)<\gamma^2\) )に対して内部的にこのγ反復を実行しています(Doyle, Glover, Khargonekar & Francis, 1989)。

数値検証

質量・ばね・ダンパ系 \(\ddot x + 0.4\dot x+4x=u+w\) を状態空間表現し、\(B_1=[1,0]^T\) (外乱が位置に直接加わる)、\(B_2=[0,1]^T\) (制御力)、\(C_1=\mathrm{diag}(1,0)\) (位置のみを評価)として、\(\gamma\) を二分探索でスイープし、ハミルトン行列の固有値が虚軸に最も近づく点を追跡しました。

import numpy as np
from scipy.linalg import eig, solve_continuous_are

A = np.array([[0.0, 1.0], [-4.0, -0.4]])
B1 = np.array([[1.0], [0.0]])
B2 = np.array([[0.0], [1.0]])
C1 = np.array([[1.0, 0.0], [0.0, 0.0]])
Q = C1.T @ C1

def min_abs_real(gamma):
    R = B2 @ B2.T - (1.0/gamma**2) * (B1 @ B1.T)
    H = np.block([[A, -R], [-Q, -A.T]])
    return np.min(np.abs(eig(H, right=False).real))

lo, hi = 1.5, 2.0  # lo: 虚軸上(解なし), hi: 虚軸外(解あり)
for _ in range(60):
    mid = (lo + hi) / 2
    if min_abs_real(mid) < 1e-9:
        lo = mid
    else:
        hi = mid
print(f"gamma_min (bisection) = {hi:.4f}")

gamma = hi + 0.01
R = B2 @ B2.T - (1.0/gamma**2) * (B1 @ B1.T)
X = solve_continuous_are(A, np.eye(2), Q, np.linalg.inv(R))
print("X =\n", X, "\neig(X) =", np.linalg.eigvals(X))
print("closed-loop eig(A - B2 B2^T X) =", np.linalg.eigvals(A - B2 @ B2.T @ X))

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

gamma_min (bisection) = 1.5992
X =
 [[2.36716484 0.37766372]
 [0.37766372 0.58509291]]
eig(X) = [2.44389691+0.j 0.50836083+0.j]
closed-loop eig(A - B2 B2^T X) = [-0.49254645+2.03348511j -0.49254645-2.03348511j]

境界値は \(\gamma_{\min}\approx1.599\) と求まりました。その少し上の \(\gamma=1.609\) で解いたリカッチ方程式の解 \(X\) は固有値 \(2.444,\ 0.508\) を持ち、確かに正定値(\(X\geq0\) )です。さらに、この \(X\) による状態フィードバック \(u=-B_2^TXx\) を施した閉ループ \(A-B_2B_2^TX\) の固有値は \(-0.493\pm2.033j\) となり、実部が負(安定)であることも確認できました。理論通り、ハミルトン行列が虚軸に固有値を持たなくなる境界 \(\gamma_{\min}\) を境に、安定化解の存在と閉ループ安定性が同時に成立/消滅することが数値的に確認できます。

有界実補題(Bounded Real Lemma)

制御器を介さず、固定された安定システム \(G(s)\) のH∞ノルムそのものをリカッチ方程式で特徴づける定理が有界実補題です。

定理: 状態空間実現を \(G(s)=C(sI-A)^{-1}B\) (\(A\) は安定, \(D=0\) の場合)とすると、次の2つは同値である。

  1. \(\|G\|_\infty < \gamma\)
  2. 次のARE(有界実リカッチ方程式)が安定化解 \(X=X^T\geq0\) を持つ。
\[ A^TX + XA + \gamma^{-2}XBB^TX + C^TC = 0 \]

これは前節の状態フィードバックのAREで \(B_2=0\) (制御入力がない)とした特別な場合に対応しており、上界・下界の証明も同様の完全平方の議論で行えます(詳細は Zhou, Doyle & Glover, Robust and Optimal Control, 1996, Ch.4 を参照)。直感的には、「\(G\) を外部から見たときの最悪ゲイン \(\gamma\) を、\(\gamma\) を境にリカッチ方程式が解けるか解けないかという代数的な問いに変換した」ものです。

数値検証

先ほどの \(G(s)=1/(s^2+0.4s+4)\) に対し、有界実リカッチ方程式に対応するハミルトン行列の固有値が虚軸に最も近づく \(\gamma\) を二分探索で求め、control.linfnorm が計算した実際のH∞ノルムと突き合わせます。

ss = ct.tf2ss(G)
A, B, C, D = np.array(ss.A), np.array(ss.B), np.array(ss.C), np.array(ss.D)

def brl_min_re(gamma):
    R = -(gamma**2) * np.eye(1)
    H = np.block([[A, -B @ np.linalg.inv(R) @ B.T], [-C.T @ C, -A.T]])
    return np.min(np.abs(eig(H, right=False).real))

lo, hi = 1.0, 2.0
for _ in range(60):
    mid = (lo + hi) / 2
    if brl_min_re(mid) < 1e-9:
        lo = mid
    else:
        hi = mid
print(f"BRL gamma_min (Hamiltonian) = {hi:.12f}")
print(f"control.linfnorm ||G||_inf   = {norm_inf:.12f}")

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

BRL gamma_min (Hamiltonian) = 1.256297269074
control.linfnorm ||G||_inf   = 1.256297269074

小数点以下12桁まで、ハミルトン行列の固有値配置から求めた境界値と、linfnorm が別アルゴリズム(状態空間の固有値問題を直接解く手法)で求めたH∞ノルムが完全に一致しました。これは有界実補題「\(\|G\|_\infty<\gamma \iff\) リカッチ方程式が解を持つ」を、理論と独立に実装された2つの数値計算で相互検証したことになります。

数値例:質量・ばね・ダンパ系のH∞制御とPID制御の比較

以上の理論を踏まえ、モデル化誤差を持つ具体的なプラントに対してH∞制御器を設計し、PID制御と頑健性を比較します。実装には python-control(v0.10.2)と slycot(SLICOTのPythonラッパー、hinfsyn が内部で使用)を用いました。

プラントと不確かさ

質量・ばね・ダンパ系 \(m\ddot x + c\dot x + kx = u\) を公称値 \(m=1\,\mathrm{kg}\) , \(c=0.4\,\mathrm{Ns/m}\) , \(k=4\,\mathrm{N/m}\) (固有振動数 \(\omega_n=2\,\mathrm{rad/s}\) 、減衰比 \(\zeta=0.1\) の軽減衰系)でモデル化します。

\[ G(s) = \frac{1}{ms^2+cs+k} = \frac{1}{s^2+0.4s+4} \]

これに対し、設計モデルに含めないアクチュエータの遅れ(未モデル化動特性、遮断周波数 \(p=15\,\mathrm{rad/s}\) )を不確かさとして導入します。

\[ G_{\mathrm{true}}(s) = G(s)\cdot\frac{p}{s+p} \]

この乗法的不確かさの厳密な大きさは \(\Delta_{\mathrm{mult}}(s) = G_{\mathrm{true}}/G-1 = -s/(s+p)\) で、\(|\Delta_{\mathrm{mult}}(j\omega)|=\omega/\sqrt{\omega^2+p^2}\) は低周波でほぼ0、高周波で1に漸近する典型的な未モデル化動特性の形をしています。

重み関数の設計

混合感度問題として、以下の3つの重み関数を設計しました。

  • 性能重み \(W_1(s) = \dfrac{s/1.5+2}{s+0.002}\) :直流ゲインが極めて大きく(積分作用に相当し定常偏差をほぼ0にする)、ピーク感度上限 \(M_1=1.5\) 、帯域幅目安 \(2\,\mathrm{rad/s}\) 。
  • 制御努力重み \(W_2(s)=0.1\) :制御信号の過大な高周波ゲインを抑制。
  • ロバスト性重み \(W_3(s) = \dfrac{1.15s}{s+15}\) :未モデル化遅れの厳密な誤差形状 \(s/(s+p)\) に対し15%の安全マージンを載せたもの。\(\|\Delta_{\mathrm{mult}}/W_3\|_\infty = 1/1.15\approx0.87<1\) となるよう正規化しています。
import control as ct

m, c, k_nom = 1.0, 0.4, 4.0
p_act = 15.0

def plant(k, lag=False):
    G = ct.tf([1.0], [m, c, k])
    if lag:
        G = G * ct.tf([p_act], [1, p_act])
    return G

G_nom = plant(k_nom)
W1 = ct.tf([1/1.5, 2.0], [1, 0.002])
W2 = ct.tf([0.1], [1])
W3 = ct.tf([1.15, 0], [1, p_act])

K, CL, info = ct.mixsyn(G_nom, w1=W1, w2=W2, w3=W3)
gamma = info[0]
print(f"achieved gamma = {gamma:.4f}")

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

achieved gamma = 1.1493

達成された \(\gamma=1.149\) は、性能重み・制御努力重み・ロバスト性重み3つを同時に満たす閉ループが構成できたことを意味します(\(\gamma\) が1に近いほど、要求した重みの設計仕様がほぼそのまま達成されたことになります)。

感度関数の確認とロバスト安定性の検証

L_nom = G_nom * K
S_nom = ct.feedback(1, L_nom)
T_nom = ct.minreal(1 - S_nom, verbose=False)

w = np.logspace(-2, 2, 4000)
print(f"max|S| = {np.abs(S_nom(1j*w)).max():.4f} at w={w[np.argmax(np.abs(S_nom(1j*w)))]:.2f}")
print(f"max|T| = {np.abs(T_nom(1j*w)).max():.4f}")

W3T = ct.minreal(W3 * T_nom, verbose=False)
w3t_mag = np.abs(W3T(1j*w))
print(f"max|W3*T| = {w3t_mag.max():.4f}  -> robust stability {'GUARANTEED' if w3t_mag.max()<1 else 'NOT guaranteed'}")
max|S| = 1.2118 at w=4.08
max|T| = 0.9989
max|W3*T| = 0.1393  -> robust stability GUARANTEED

ピーク感度 \(\max|S|=1.212\) は良好な値(一般に2未満が目安)です。小ゲイン定理の核心である \(\max|W_3T|=0.139<1\) が成立しており、これは「アクチュエータ遅れとして想定した不確かさの正規化された大きさ \(\|\Delta\|_\infty\leq1\) 」の範囲内であれば、どんな不確かさが実現してもこの閉ループは安定であり続けることを理論的に保証しています(先述の小ゲイン定理の数値検証と同じ \(M=W_3T\) です)。

下図は、感度関数 \(S\) ・相補感度関数 \(T\) と、それぞれの設計仕様 \(1/W_1\) (性能上限)・\(1/W_3\) (ロバスト性上限)を重ねたものです。\(S\) が低周波で \(1/W_1\) の下に収まり(良好な追従性能)、\(T\) が高周波で \(1/W_3\) の下に収まっている(ロバスト性の確保)ことがトレードオフとして視覚化されています。

Mixed-sensitivity H-infinity design: S, T, and weight bounds

PID制御器との比較

比較対象として、ITAE最適規範(3次系)でチューニングしたPID制御器を用意しました。閉ループ特性多項式を \(s^3+1.75\omega_ns^2+2.15\omega_n^2s+\omega_n^3\) (\(\omega_n=3\,\mathrm{rad/s}\) )に一致させると、

\[ K_d = 1.75\omega_n - c = 4.850,\quad K_p = 2.15\omega_n^2-k_{\mathrm{nom}}=15.350,\quad K_i=\omega_n^3=27.000 \]

というゲインが得られます。

wn_c = 3.0
Kd = 1.75*wn_c - c
Kp = 2.15*wn_c**2 - k_nom
Ki = wn_c**3
PID = ct.tf([Kd, Kp, Ki], [1, 0])
print(f"PID gains: Kp={Kp:.3f} Ki={Ki:.3f} Kd={Kd:.3f}")
PID gains: Kp=15.350 Ki=27.000 Kd=4.850

この2つの制御器を、次の4つのシナリオでステップ応答比較しました。

シナリオ制御器オーバーシュート2%整定時間閉ループ最大 \(\mathrm{Re}(\text{極})\)
公称(\(k=4.0\) 、遅れなし)H∞0.90%1.300 s−0.2000
公称PID17.08%1.984 s−1.5629
設計想定の不確かさ(\(k=4.0\) +遅れ)H∞4.03%2.040 s−0.2000
設計想定の不確かさPID31.98%2.322 s−1.7792
ストレステスト(\(k=3.2\) +遅れ)H∞20.74%97.19 s−0.0230
ストレステスト(\(k=3.2\) +遅れ)PID35.33%2.454 s−1.6922
ストレステスト(\(k=4.8\) +遅れ)H∞8.15%6.369 s−0.4106
ストレステスト(\(k=4.8\) +遅れ)PID28.83%2.178 s−1.8475

「設計想定の不確かさ」(W3で明示的に考慮したアクチュエータ遅れのみ)の範囲では、H∞制御はオーバーシュート4.03%とPIDの31.98%を大きく下回り、小ゲイン定理の保証通りロバスト性の恩恵が明確に表れています。下図左(nominal)・右(ストレステスト)でこの違いを視覚的に確認できます。

Step response comparison: H-infinity vs ITAE-PID under nominal and stress-test plants

注意点:小ゲイン定理の保証は「想定した不確かさの範囲内」に限られる

一方で「ストレステスト」(アクチュエータ遅れに加え、\(W_3\) の設計に含めていないばね定数 \(k\) の変動を同時に加えたケース)では、興味深い逆転が起こりました。\(k=3.2\) (公称より20%低下)のケースでは、H∞制御の閉ループ最大極の実部はわずか \(-0.023\) まで虚軸に接近し、2%整定時間は97.19秒(公称時の1.30秒の約75倍)にまで悪化しています。対照的に、PID制御はこのシナリオでも2.45秒で収束し、オーバーシュートこそ大きいものの実用上はむしろ頑健に振る舞いました。

さらに \(k\) を下げていくと、H∞制御は \(k\approx3.086\) (公称より約22.9%低下、アクチュエータ遅れとの複合)で閉ループが不安定化することが二分探索で確認できました。

lo, hi = 2.8, 4.0
for _ in range(40):
    mid = (lo + hi) / 2
    G_true = plant(mid, lag=True)
    stable = np.all(ct.poles(ct.feedback(G_true*K, 1)).real < 0)
    if stable:
        hi = mid
    else:
        lo = mid
print(f"H-infinity controller unstable below k = {hi:.4f}")
H-infinity controller unstable below k = 3.0858

これは決してH∞制御が劣っているという意味ではなく、小ゲイン定理・BRLが与える頑健安定性の保証は、設計時に \(\Delta\) としてモデル化した不確かさの構造とサイズにのみ及ぶという、極めて重要な限界を示しています。今回の \(W_3\) はアクチュエータ遅れという特定の不確かさだけを表現しており、ばね定数という別種の(構造化された・パラメトリックな)不確かさは一切考慮されていません。\(\max|W_3T|=0.139\) という余裕のある小ゲイン条件は、あくまで「アクチュエータ遅れ方向の不確かさ」に対してのみ意味を持ち、それ以外の方向の摂動に対しては何の保証も与えません。

このような複数種類・構造化された不確かさを同時に扱うには、多入力多出力の不確かさブロックに対して個別に正規化を行う構造化特異値 \(\mu\) (structured singular value)とD-K反復によるμ-synthesisが本来必要です(Doyle, 1982; Skogestad & Postlethwaite, Multivariable Feedback Control, 2005, Ch.8)。本記事のスコープを超えるため深入りしませんが、「H∞制御=万能にロバスト」ではなく「設計時に明示した不確かさに対してのみ、数学的に証明された保証を持つ」という点は、実務で不確かさをモデル化する際の重要な注意点です。

外乱抑制性能

最後に、公称プラントに対する出力外乱の抑制性能を比較します。

Sdist_hinf = ct.feedback(1, G_nom*K)
Sdist_pid = ct.feedback(1, G_nom*PID)

積分作用を持つ点は両者共通のため定常偏差はいずれも0に収束しますが、過渡応答の収束速度はステップ応答と同じ傾向(H∞: 1.30秒、PID: 1.98秒、公称プラントにおいて)を示しました。

まとめ

  • H∞ノルム \(\|T_{zw}\|_\infty=\sup_\omega\bar\sigma(T_{zw}(j\omega))\) は、有限エネルギー(\(L_2\) )入力に対する誘導ノルムであり、最悪ケースのゲインを表す。数値実験で control.linfnorm とブルートフォース探索が小数点以下6桁まで一致することを確認した。
  • 小ゲイン定理は \(\|M\|_\infty\|\Delta\|_\infty<1\) を、ノイマン級数の収束(作用素論)とナイキストの安定判別法の連続変形(複素関数論)の2通りで証明できる十分条件であり、数値実験では理論的な安定限界(\(1/\|M\|_\infty=7.179\) )と実際の不安定化点(\(\delta\approx7.255\) )が約1%の保守性で一致した。
  • 標準H∞制御問題は一般化プラント \(P\) の下側線形分数変換 \(F_l(P,K)=P_{11}+P_{12}K(I-P_{22}K)^{-1}P_{21}\) として定式化され、混合感度問題はその特別な場合として具体的に書き下せる。
  • H∞代数リカッチ方程式は、状態フィードバック問題の鞍点条件(HJI方程式)の平方完成から導出でき、対応するハミルトン行列が虚軸上に固有値を持たないことが解の存在条件(γ反復の判定基準)となる。数値例では \(\gamma_{\min}\approx1.599\) の境界で解の存在と閉ループ安定性が同時に消滅することを確認した。
  • 有界実補題はこの考え方を制御器なしの固定システムに適用したもので、\(\|G\|_\infty<\gamma\) とリカッチ方程式の可解性が同値であることを、独立な2つのアルゴリズム(ハミルトン行列の固有値判定と linfnorm)で小数点以下12桁一致という形で数値的に検証した。
  • 数値例では、質量・ばね・ダンパ系にアクチュエータ遅れの不確かさを持たせ、control.mixsyn(\(\gamma=1.149\) )で設計したH∞制御器が、想定した不確かさの範囲内ではPIDより大幅にロバスト(オーバーシュート4.03% vs 31.98%)である一方、設計時に想定していない種類の不確かさ(ばね定数の変動)に対しては保証が及ばず、\(k\approx3.086\) で不安定化することを確認した。ロバスト制御の保証は、常に「何を不確かさとしてモデル化したか」に依存することを、実測をもって示した。

関連記事

参考文献

  • Doyle, J. C., Glover, K., Khargonekar, P. P., & Francis, B. A. (1989). “State-space solutions to standard H2 and H∞ control problems.” IEEE Transactions on Automatic Control, 34(8), 831–847.
  • Zhou, K., Doyle, J. C., & Glover, K. (1996). Robust and Optimal Control. Prentice Hall.
  • Skogestad, S., & Postlethwaite, I. (2005). Multivariable Feedback Control: Analysis and Design (2nd ed.). Wiley.
  • Doyle, J. C. (1982). “Analysis of feedback systems with structured uncertainties.” IEE Proceedings D, 129(6), 242–250.
  • python-control: hinfsyn documentation
  • SLICOT: Subroutine Library in Systems and Control Theory

関連ツール