ロバスト制御の基礎 — H∞制御入門

工場の制御システムは、設計時に想定した条件通りに常に動くとは限りません。温度変化で材料の特性が変わったり、経年劣化でパラメータがずれたり、想定外の外乱が加わったりします。LQR制御やLQG制御は「ノミナルモデル（理想的なモデル）」が正確である前提で最適化しますが、モデルと実際のシステムにずれがあると性能が大幅に劣化する場合があります。

たとえば、ロボットアームの制御器を設計するとき、アームの質量が設計値の0.5 kgではなく実際には0.7 kgだったら何が起きるでしょうか。LQR制御器はこの40%の誤差に対して大きなオーバーシュートや振動を起こすかもしれません。あるいは、航空機の飛行制御では、高度や速度が変わるたびに空力特性が大きく変化します。特定の飛行条件でチューニングした制御器が、他の条件で不安定になっては困ります。

ロバスト制御（robust control）は、こうしたモデルの不確実性を明示的に考慮し、「最悪のケースでも一定の性能を保証する」制御を設計する手法です。その代表がH∞制御（H-infinity control）であり、外乱から出力へのゲインの最大値（$\mathcal{H}_\infty$ ノルム）を最小化します。

H∞制御を理解すると、以下のような応用が開けます。

パラメータ変動への耐性: プラントのパラメータが変化しても安定かつ一定の性能を保証
未モデル化ダイナミクスへの対応: 高周波域の不確実性に対するロバスト安定性
航空宇宙: 飛行条件が大きく変化する航空機の飛行制御
混合感度設計: 追従性能、外乱抑制、ロバスト安定性のバランス設計

本記事の内容

不確実性のモデリング（加法的、乗法的）
$\mathcal{H}_\infty$ ノルムの定義と意味
ロバスト安定性条件（スモールゲイン定理）
混合感度問題の定式化
一般化プラント表現
Pythonでの設計と検証

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

フィードバック制御の基礎 — 感度関数・相補感度関数
ボード線図の描き方と読み方 — 周波数応答の基礎
LQG制御 — H2最適制御との比較

不確実性のモデリング

なぜ不確実性を明示的にモデル化するのか

全ての数学モデルは「近似」です。実際のシステムとモデルの間には必ず差異があります。この差異をモデル不確実性（model uncertainty）と呼びます。

不確実性の発生源は大きく分けて以下の3種類です。

パラメータ不確実性: 物理パラメータ（質量、ばね定数、抵抗値など）の公称値からのずれ。製造公差や経年劣化により発生する
未モデル化ダイナミクス: 設計モデルには含めなかった高次のダイナミクス（構造物の高次モード、アクチュエータの遅延、非線形性など）
外部外乱: 風、振動、電磁干渉、測定ノイズなど外部から加わる信号

古典制御では、ゲイン余裕や位相余裕といった「余裕」で不確実性に対する耐性を暗黙的に評価していました。しかし、MIMO系ではゲイン余裕と位相余裕だけでは不十分です。H∞制御は、不確実性を周波数領域の重み付きノルムで定量的にモデル化し、それに対する最悪ケースの性能を最適化します。

加法的不確実性

実際のプラント $G_{\text{true}}(s)$ がノミナルモデル $G(s)$ と不確実性 $\Delta_a(s)$ の和で表される場合を加法的不確実性と呼びます。

$$ G_{\text{true}}(s) = G(s) + \Delta_a(s) $$

加法的不確実性の大きさは周波数に依存します。低周波域ではモデルが比較的正確で不確実性が小さく、高周波域ではモデル化が難しく不確実性が大きくなるのが一般的です。

この周波数依存性を表現するために、不確実性を重み関数 $W_a(s)$ と正規化された不確実性 $\Delta(s)$（$\|\Delta\|_\infty \leq 1$）の積で表します。

$$ \Delta_a(s) = W_a(s)\Delta(s), \quad \|\Delta\|_\infty \leq 1 $$

$W_a(s)$ は不確実性の「周波数ごとの上界」を表す安定な伝達関数です。$|W_a(j\omega)|$ が各周波数における不確実性の最大の大きさを示します。

乗法的不確実性

不確実性をノミナルプラントからの相対的なずれとして表現するのが乗法的不確実性です。

$$ G_{\text{true}}(s) = (1 + W_m(s)\Delta(s))G(s), \quad \|\Delta\|_\infty \leq 1 $$

ここで $W_m(s)$ は乗法的不確実性の重みです。乗法的不確実性の物理的な意味は

$$ \frac{G_{\text{true}}(j\omega) – G(j\omega)}{G(j\omega)} \leq |W_m(j\omega)| $$

であり、$|W_m(j\omega)|$ は周波数 $\omega$ における相対誤差の上界を表します。

たとえば、$|W_m(j\omega)|$ が低周波で0.1（10%の誤差）、高周波で10（1000%の誤差）であるとき、低周波域ではモデルが正確だが、高周波域ではモデルの信頼性が低いことを表現しています。

不確実性の重み関数の設定例

実際の設計では、不確実性の重み関数を具体的にどう設定するかが重要です。典型的な乗法的不確実性の重みは次のような1次の伝達関数で表されます。

$$ W_m(s) = \frac{\tau s + r_0}{(\tau/r_\infty)s + 1} $$

ここで、$r_0$ は低周波域での相対誤差、$r_\infty$ は高周波域での相対誤差、$\tau$ はクロスオーバー周波数を決めるパラメータです。

たとえば、低周波域で10%、高周波域で200%の不確実性を想定する場合は $r_0 = 0.1$、$r_\infty = 2$ と設定します。

不確実性のモデリングができたところで、次にこの不確実性に対する安定性と性能を評価するための数学的道具である $\mathcal{H}_\infty$ ノルムを導入します。

$\mathcal{H}_\infty$ ノルムの定義と性質

ノルムの定義

伝達関数 $G(s)$ の $\mathcal{H}_\infty$ ノルムは、周波数応答のゲインの最大値として定義されます。

$$ \begin{equation} \|G\|_\infty = \sup_\omega |G(j\omega)| \end{equation} $$

MIMO系（多入力多出力系）では、スカラーのゲイン $|G(j\omega)|$ の代わりに最大特異値 $\bar{\sigma}[G(j\omega)]$ を使います。

$$ \begin{equation} \|G\|_\infty = \sup_\omega \bar{\sigma}[G(j\omega)] \end{equation} $$

ここで最大特異値 $\bar{\sigma}[G(j\omega)]$ とは、行列 $G(j\omega)$ の特異値分解（SVD）における最大の特異値です。

物理的意味 — 最悪ケースの増幅率

$\mathcal{H}_\infty$ ノルムの物理的意味を2つの観点から理解できます。

周波数領域の解釈: $\|G\|_\infty$ は、全ての周波数にわたるゲインの最大値です。つまり、「最も増幅される周波数でどれだけ増幅されるか」を示します。ボード線図のゲイン曲線のピーク値に対応します。

時間領域の解釈: $\|G\|_\infty$ は、エネルギー有限の入力信号に対する出力のエネルギーゲインの最大値でもあります。

$$ \|G\|_\infty = \sup_{u \neq 0} \frac{\|y\|_2}{\|u\|_2} $$

ここで $\|u\|_2 = \sqrt{\int_0^\infty |u(t)|^2 dt}$ は信号の $L_2$ ノルム（エネルギーの平方根）です。

この等式はパーセバルの定理から導かれます。時間領域の $L_2$ ノルムと周波数領域の $L_2$ ノルムが等しいことを利用して

$$ \|y\|_2^2 = \frac{1}{2\pi}\int_{-\infty}^{\infty} |G(j\omega)U(j\omega)|^2 d\omega \leq \|G\|_\infty^2 \cdot \frac{1}{2\pi}\int_{-\infty}^{\infty} |U(j\omega)|^2 d\omega = \|G\|_\infty^2 \|u\|_2^2 $$

と上界が得られ、等号は $|G(j\omega)|$ が最大となる周波数に集中した入力で達成されます。

$\mathcal{H}_2$ ノルムとの比較

H∞制御を理解する上で、LQG制御で使われる $\mathcal{H}_2$ ノルムとの違いを知っておくことが重要です。

$$ \|G\|_2 = \sqrt{\frac{1}{2\pi}\int_{-\infty}^{\infty} |G(j\omega)|^2 d\omega} $$

性質	$\mathcal{H}_2$ ノルム	$\mathcal{H}_\infty$ ノルム
定義	周波数応答の2乗積分の平方根	周波数応答のゲインの最大値
時間領域	インパルス応答のエネルギー	最悪ケースのエネルギーゲイン
最適化の意味	平均的な性能の最適化	最悪ケースの性能の最適化
対応する制御	LQG制御	H∞制御
ノイズ仮定	ホワイトノイズ（確率的）	有界エネルギー（確定的）

$\mathcal{H}_2$ 最適化は「平均的に良い」制御器を設計し、$\mathcal{H}_\infty$ 最適化は「最悪でも悪くない」制御器を設計します。外乱の統計的性質が既知のときは $\mathcal{H}_2$ が有効で、不確実性が大きく最悪ケースを保証したいときは $\mathcal{H}_\infty$ が適しています。

$\mathcal{H}_\infty$ ノルムの定義が明確になったところで、このノルムを使ってロバスト安定性の条件を定式化していきます。

スモールゲイン定理

定理の主張

ロバスト安定性の最も基本的な結果がスモールゲイン定理（Small Gain Theorem）です。

定理（スモールゲイン定理）: 安定なシステム $M(s)$ と安定な不確実性 $\Delta(s)$（$\|\Delta\|_\infty \leq 1$）のフィードバック結合を考えます。このフィードバック系が、$\|\Delta\|_\infty \leq 1$ を満たす全ての $\Delta$ に対して内部安定であるための十分条件は

$$ \begin{equation} \|M\|_\infty < 1 \end{equation} $$

です。

直感的な理解

スモールゲイン定理は、フィードバックループの「一周ゲイン」が1未満であれば信号が一周するたびに減衰し、系が安定に保たれるという直感に基づいています。

マイクとスピーカーの間でハウリング（フィードバック発振）が起きる状況を思い浮かべてください。スピーカーの音がマイクに戻り、増幅されてまたスピーカーから出る。このループゲインが1を超えると発振しますが、1未満であれば発振しません。スモールゲイン定理はこの直感を厳密化したものです。

乗法的不確実性に対するロバスト安定条件

乗法的不確実性 $G_{\text{true}}(s) = (1 + W_m(s)\Delta(s))G(s)$ のもとで、制御器 $C(s)$ によるフィードバック系のロバスト安定条件を導出します。

ノミナルの閉ループ系の相補感度関数は

$$ T(s) = \frac{C(s)G(s)}{1 + C(s)G(s)} $$

です。不確実性を含む系の特性方程式を整理すると、ロバスト安定条件は

$$ \begin{equation} \|W_m T\|_\infty < 1 \end{equation} $$

と書けます。

この条件の導出を示します。実際のプラントでの閉ループ特性方程式は

$$ 1 + C(s)G_{\text{true}}(s) = 0 $$

$G_{\text{true}} = (1 + W_m \Delta)G$ を代入すると

$$ 1 + C(1 + W_m\Delta)G = 1 + CG + CG W_m\Delta = (1 + CG)(1 + \frac{CG}{1+CG}W_m\Delta) = 0 $$

ノミナル系が安定（$1 + CG \neq 0$）であれば、ロバスト安定性の条件は

$$ 1 + T W_m \Delta \neq 0, \quad \forall \|\Delta\|_\infty \leq 1 $$

スモールゲイン定理より、これは $\|W_m T\|_\infty < 1$ と等価です。

周波数ごとに書くと

$$ |W_m(j\omega)| \cdot |T(j\omega)| < 1, \quad \forall \omega $$

加法的不確実性に対するロバスト安定条件

同様に、加法的不確実性 $G_{\text{true}} = G + W_a\Delta$ に対するロバスト安定条件は

$$ \begin{equation} \left\|W_a \frac{C}{1+CG}\right\|_\infty < 1 \end{equation} $$

となります。$\frac{C}{1+CG} = CS$ を入力感度関数と呼びます。

ロバスト安定条件が明確になったところで、次に安定性だけでなく性能も同時に設計する枠組みである混合感度問題を見ていきます。

混合感度問題

設計のトレードオフ

実際の制御設計では、以下の複数の要求を同時に満たす必要があります。

追従性能: 目標値変化に対して出力が素早く正確に追従する → 感度関数 $S$ を小さくしたい
外乱抑制: 外乱の影響を小さくする → 感度関数 $S$ を小さくしたい
ノイズ抑制: 計測ノイズの影響を小さくする → 相補感度関数 $T$ を小さくしたい
ロバスト安定性: 不確実性に対して安定を保つ → 相補感度関数 $T$ を小さくしたい
制御入力の制限: アクチュエータの飽和を避ける → 入力感度関数 $CS$ を小さくしたい

しかし、$S + T = 1$ の関係（全パスバンドで成立）があるため、$S$ と $T$ を同時に全周波数で小さくすることは原理的に不可能です。

$$ S(s) + T(s) = \frac{1}{1+L} + \frac{L}{1+L} = 1 $$

このため、周波数帯域ごとに $S$ と $T$ のどちらを優先するかを設定するのが混合感度問題の考え方です。

混合感度問題の定式化

混合感度問題は、重み付き感度関数と相補感度関数の $\mathcal{H}_\infty$ ノルムを最小化する問題です。

$$ \begin{equation} \min_C \left\|\begin{bmatrix} W_1 S \\ W_2 T \end{bmatrix}\right\|_\infty \end{equation} $$

ここで

$W_1(s)$: 感度関数への重み — 低周波域での追従性能と外乱抑制を指定
$W_2(s)$: 相補感度関数への重み — 高周波域でのロバスト安定性とノイズ抑制を指定

この $\mathcal{H}_\infty$ ノルムを $\gamma$ 未満にするための条件は

$$ \left\|\begin{bmatrix} W_1 S \\ W_2 T \end{bmatrix}\right\|_\infty < \gamma $$

であり、これは全ての周波数 $\omega$ で

$$ |W_1(j\omega) S(j\omega)|^2 + |W_2(j\omega) T(j\omega)|^2 < \gamma^2 $$

を意味します。特に各項について

$$ |W_1(j\omega) S(j\omega)| < \gamma, \quad |W_2(j\omega) T(j\omega)| < \gamma $$

が成り立つので、$\gamma = 1$ を目指す場合

$$ |S(j\omega)| < \frac{1}{|W_1(j\omega)|}, \quad |T(j\omega)| < \frac{1}{|W_2(j\omega)|} $$

つまり、$1/|W_1|$ は感度関数の上界、$1/|W_2|$ は相補感度関数の上界を規定します。

重み関数の典型的な設定

$W_1(s)$ の設定方針

$W_1$ は低周波域で大きく（追従性を要求）、高周波域で小さく設定します。典型的な形は

$$ W_1(s) = \frac{s/M + \omega_b}{s + \omega_b \epsilon} $$

ここで – $\omega_b$: 所望の帯域幅 — この周波数以下で追従性を要求 – $M$: $|S|$ のピーク値の上界 — 通常2以下（$M > 2$ はロバスト性が低い） – $\epsilon$: 定常偏差の許容値 — 小さいほど定常偏差が小さい

低周波極限（$\omega \to 0$）では $|W_1| \to 1/\epsilon$ となり、$|S| < \epsilon$ が要求されます。これは定常偏差が $\epsilon$ 以下であることを意味します。

$W_2(s)$ の設定方針

$W_2$ は乗法的不確実性の重み $W_m$ と関連付けて設定します。ロバスト安定条件 $\|W_m T\|_\infty < 1$ を保証するために、$W_2 \geq W_m$ とすることが多いです。

$$ W_2(s) = \frac{s + \omega_t/M_t}{\epsilon_t s + \omega_t} $$

高周波域（$\omega \to \infty$）では $|W_2| \to 1/\epsilon_t$ となり、$|T| < \epsilon_t$ が要求されます。これはセンサノイズの減衰を保証します。

3重み混合感度問題

制御入力の大きさも制限したい場合は、入力感度関数 $CS$ も含めた3重みの混合感度問題を考えます。

$$ \begin{equation} \min_C \left\|\begin{bmatrix} W_1 S \\ W_2 T \\ W_3 CS \end{bmatrix}\right\|_\infty \end{equation} $$

ここで $W_3(s)$ は制御入力への重みで、アクチュエータの帯域や飽和制限を反映します。

混合感度問題を解くための一般的な枠組みが、次に述べる一般化プラント表現です。

一般化プラント表現

$P$-$K$ 構造

H∞制御問題を統一的に扱うために、一般化プラント（generalized plant）$P(s)$ と制御器 $K(s)$ の構造に帰着させます。

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

ここで – $w$: 外生入力（外乱、目標値、ノイズ） – $u$: 制御入力 – $z$: 評価出力（最小化したい信号） – $y$: 観測出力（制御器が利用できる信号）

制御器 $u = K(s)y$ で閉ループを構成すると、$w$ から $z$ への閉ループ伝達関数は下側線形分数変換（lower linear fractional transformation; lower LFT）で表されます。

$$ \begin{equation} T_{zw}(s) = \mathcal{F}_l(P, K) = P_{11} + P_{12}K(I – P_{22}K)^{-1}P_{21} \end{equation} $$

H∞制御問題は、このLFTの $\mathcal{H}_\infty$ ノルムを最小化する $K$ を求める問題に帰着します。

$$ \min_K \|\mathcal{F}_l(P, K)\|_\infty $$

混合感度問題の一般化プラント

混合感度問題の2重み版を一般化プラント形式に書き下します。

フィードバック構成で $L = CG$ とすると、$S = 1/(1+L)$、$T = L/(1+L)$ です。外生入力 $w = r$（目標値）、制御入力 $u$、評価出力 $z = [W_1 e, W_2 y]^T$、観測出力 $y_{\text{obs}} = e = r – y$ とすると、一般化プラントは

$$ P(s) = \begin{bmatrix} W_1 & -W_1 G \\ 0 & W_2 G \\ I & -G \end{bmatrix} $$

と書けます。ここで上の2行が $z$、最後の行が $y_{\text{obs}}$ に対応します。

リカッチ方程式による解法

H∞制御問題の解は、2つのリカッチ方程式の解から構成されます。$\gamma > \gamma^*$（最適値）のとき、以下の2つのリカッチ方程式の安定化解 $X \geq 0$ と $Y \geq 0$ が存在します。

状態フィードバック用リカッチ方程式:

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

オブザーバ用リカッチ方程式:

$$ AY + YA^T + B_1B_1^T + Y(\gamma^{-2}C_1^TC_1 – C_2^TC_2)Y = 0 $$

これらの解が存在し、スペクトル半径条件 $\rho(XY) < \gamma^2$ を満たすとき、達成可能な最小の $\gamma$ が $\gamma^*$ です。

制御器は以下の形で与えられます。

$$ K(s) = -F(sI – A_K)^{-1}L $$

ここで $F$、$L$、$A_K$ は $X$、$Y$、$\gamma$ から構成されます。数式の詳細は省略しますが、重要な点はLQG制御と同様に状態フィードバックとオブザーバの組み合わせ構造を持つことです。ただし、H∞制御ではLQGの分離定理は一般には成り立たず、$X$ と $Y$ が結合しています（スペクトル半径条件）。

実際の設計ではリカッチ方程式を手で解くことはなく、数値ソルバを使います。それでは、Pythonで実際に混合感度設計を行ってみましょう。

Pythonでの実装と検証

パラメータ変動に対するロバスト性の比較

まず、パラメータ変動を含む複数のプラントに対して、単純な制御器とロバスト補償器の性能を比較します。

import numpy as np
import matplotlib.pyplot as plt
import control as ctrl

# ノミナルプラント: G = 10/[s(s+1)(s+10)]
G_nom = ctrl.TransferFunction([10],
    np.polymul([1, 0], np.polymul([1, 1], [1, 10])))

# パラメータ変動を含む複数のプラント
plants = {"Nominal": G_nom}
for delta in [-0.3, -0.15, 0.15, 0.3]:
    G_pert = ctrl.TransferFunction(
        [10 * (1 + delta)],
        np.polymul([1, 0],
                   np.polymul([1, 1 + delta * 0.5],
                              [1, 10 + delta * 3])))
    plants[f"delta={delta:+.2f}"] = G_pert

# 単純な比例制御器
C_simple = ctrl.TransferFunction([5], [1])

# ロバスト補償器（位相進み補償）
C_robust = ctrl.TransferFunction([3, 6], [1, 20])

omega = np.logspace(-2, 2, 500)
fig, axes = plt.subplots(2, 2, figsize=(14, 9))

controllers = [(C_simple, "Simple (K=5)"),
               (C_robust, "Robust compensator")]

# 感度関数の比較
for idx, (C, C_name) in enumerate(controllers):
    ax = axes[0, idx]
    for label, G in plants.items():
        L = C * G
        S = ctrl.feedback(1, L)
        mag_S, _, w_S = ctrl.frequency_response(S, omega)
        ax.semilogx(w_S, 20 * np.log10(np.abs(mag_S)),
                    linewidth=1.5, label=label)
    ax.axhline(0, color="black", linestyle="--", alpha=0.4)
    ax.set_xlabel(r"$\omega$ [rad/s]", fontsize=11)
    ax.set_ylabel("|S| [dB]", fontsize=11)
    ax.set_title(f"Sensitivity: {C_name}", fontsize=12)
    ax.legend(fontsize=8, loc="lower right")
    ax.grid(True, which="both", alpha=0.3)
    ax.set_ylim([-30, 15])

# ステップ応答の比較
t = np.linspace(0, 8, 500)
for idx, (C, C_name) in enumerate(controllers):
    ax = axes[1, idx]
    for label, G in plants.items():
        L = C * G
        G_cl = ctrl.feedback(L)
        poles = ctrl.poles(G_cl)
        if all(p.real < 0 for p in poles):
            t_out, y_out = ctrl.step_response(G_cl, T=t)
            ax.plot(t_out, y_out, linewidth=1.5, label=label)
        else:
            ax.text(0.3, 0.5, f"{label}: Unstable!",
                    transform=ax.transAxes, color="red",
                    fontsize=10, ha="center")
    ax.axhline(1, color="black", linestyle="--", alpha=0.4)
    ax.set_xlabel("Time [s]", fontsize=11)
    ax.set_ylabel("Output", fontsize=11)
    ax.set_title(f"Step Response: {C_name}", fontsize=12)
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("h_infinity_robust_comparison.png",
            dpi=150, bbox_inches="tight")
plt.show()

この結果から、ロバスト性の違いが明確に確認できます。

感度関数（上段）: 単純な制御器（左）ではパラメータ変動に対して感度関数のピークが大きく変動し、一部のケースで0 dBを大きく超えています。ピークが大きい（$\|S\|_\infty$ が大きい）ということは、特定の周波数の外乱が増幅されることを意味します。ロバスト補償器（右）では変動が抑制され、全てのケースで感度関数のピークが制限されています
ステップ応答（下段）: 単純な制御器はノミナルモデルでは良好な応答を示しますが、パラメータが $\pm 30\%$ 変動すると大きなオーバーシュートや振動が発生し、不安定になるケースもあります。ロバスト補償器は全てのパラメータ変動に対して安定な応答を維持し、オーバーシュートも制限されています

混合感度設計の可視化

次に、重み関数と設計結果の関係を可視化します。

import numpy as np
import matplotlib.pyplot as plt
import control as ctrl

# ノミナルプラント
G = ctrl.TransferFunction([10],
    np.polymul([1, 0], np.polymul([1, 1], [1, 10])))

# 重み関数の設定
wb = 1.0    # 所望の帯域幅
M = 2.0     # 感度のピーク上界
eps = 0.01  # 定常偏差の許容値

# W1: 低周波で大きく（追従性要求）
W1 = ctrl.TransferFunction([1/M, wb], [1, wb * eps])

# W2: 高周波で大きく（ロバスト安定性要求）
wt = 50.0
Mt = 2.0
eps_t = 0.01
W2 = ctrl.TransferFunction([1, wt/Mt], [eps_t, wt])

# 手動設計のロバスト補償器
C = ctrl.TransferFunction(
    np.polymul([0.8, 1.6], [1, 3]),
    np.polymul([1, 20], [1, 0.01]))

# 閉ループ伝達関数
L = C * G
S = ctrl.feedback(1, L)
T = ctrl.feedback(L)

omega = np.logspace(-3, 3, 1000)

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# (a) 重み関数の逆数とS, T
ax = axes[0]
mag_S, _, w_S = ctrl.frequency_response(S, omega)
mag_T, _, w_T = ctrl.frequency_response(T, omega)
mag_1W1, _, w_1 = ctrl.frequency_response(1/W1, omega)
mag_1W2, _, w_2 = ctrl.frequency_response(1/W2, omega)

ax.semilogx(w_S, 20*np.log10(np.abs(mag_S)),
            'b-', linewidth=2, label='$|S|$')
ax.semilogx(w_T, 20*np.log10(np.abs(mag_T)),
            'r-', linewidth=2, label='$|T|$')
ax.semilogx(w_1, 20*np.log10(np.abs(mag_1W1)),
            'b--', linewidth=1.5, label='$1/|W_1|$')
ax.semilogx(w_2, 20*np.log10(np.abs(mag_1W2)),
            'r--', linewidth=1.5, label='$1/|W_2|$')
ax.set_xlabel(r"$\omega$ [rad/s]", fontsize=11)
ax.set_ylabel("Magnitude [dB]", fontsize=11)
ax.set_title("S, T vs Weight Bounds", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, which="both", alpha=0.3)
ax.set_ylim([-60, 20])

# (b) 重み付き感度関数
ax = axes[1]
mag_W1S, _, w_ws = ctrl.frequency_response(W1*S, omega)
mag_W2T, _, w_wt = ctrl.frequency_response(W2*T, omega)

ax.semilogx(w_ws, 20*np.log10(np.abs(mag_W1S)),
            'b-', linewidth=2, label='$|W_1 S|$')
ax.semilogx(w_wt, 20*np.log10(np.abs(mag_W2T)),
            'r-', linewidth=2, label='$|W_2 T|$')
ax.axhline(0, color="black", linestyle="--",
           linewidth=1.5, alpha=0.7, label='$\\gamma = 1$ (0 dB)')
ax.set_xlabel(r"$\omega$ [rad/s]", fontsize=11)
ax.set_ylabel("Magnitude [dB]", fontsize=11)
ax.set_title("Weighted Sensitivity", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, which="both", alpha=0.3)
ax.set_ylim([-40, 10])

# (c) ナイキスト線図
ax = axes[2]
mag_L, phase_L, w_L = ctrl.frequency_response(L, omega)
re_L = np.real(mag_L * np.exp(1j * phase_L))
im_L = np.imag(mag_L * np.exp(1j * phase_L))
ax.plot(re_L, im_L, 'b-', linewidth=1.5, label='$L(j\\omega)$')
ax.plot(-1, 0, 'rx', markersize=10, markeredgewidth=2,
        label='Critical point')
theta = np.linspace(0, 2*np.pi, 100)
ax.plot(np.cos(theta) - 1, np.sin(theta),
        'g--', linewidth=1, alpha=0.5,
        label='Unit circle at -1')
ax.set_xlabel("Real", fontsize=11)
ax.set_ylabel("Imag", fontsize=11)
ax.set_title("Nyquist Diagram", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
ax.set_xlim([-3, 1])
ax.set_ylim([-2, 2])

plt.tight_layout()
plt.savefig("h_infinity_mixed_sensitivity.png",
            dpi=150, bbox_inches="tight")
plt.show()

このグラフから、混合感度設計の本質が読み取れます。

左図（S, T と重み関数の逆数）: 感度関数 $|S|$（青実線）が低周波域で $1/|W_1|$（青破線）の下にあり、追従性の要求を満たしています。相補感度関数 $|T|$（赤実線）が高周波域で $1/|W_2|$（赤破線）の下にあり、ロバスト安定性の要求を満たしています。中間の周波数帯が $S$ から $T$ への遷移が起こる領域であり、制御帯域幅に対応します
中央図（重み付き感度関数）: $|W_1 S|$ と $|W_2 T|$ の両方が0 dB以下であれば、$\gamma < 1$ の条件を満たしています。0 dBに近づいている部分が設計のボトルネック（最も厳しい制約）です
右図（ナイキスト線図）: 一巡伝達関数 $L(j\omega)$ が臨界点 $-1$ から十分離れていることが確認でき、ゲイン余裕と位相余裕が十分であることがわかります。$-1$ を中心とする単位円との距離が $1/\|S\|_\infty$ に対応します

感度のピーク値と安定余裕の関係

H∞制御の文脈では、感度関数のピーク値 $M_s = \|S\|_\infty$ が重要な設計指標です。$M_s$ とゲイン余裕・位相余裕の間には以下の関係があります。

import numpy as np
import matplotlib.pyplot as plt

# M_sとゲイン余裕・位相余裕の関係
Ms_values = np.linspace(1.01, 5, 200)

# ゲイン余裕の下界: GM >= Ms/(Ms-1)
GM_lower = Ms_values / (Ms_values - 1)

# 位相余裕の下界: PM >= 2*arcsin(1/(2*Ms))
PM_lower = 2 * np.arcsin(1 / (2 * Ms_values)) * 180 / np.pi

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# (a) M_s vs ゲイン余裕
ax = axes[0]
ax.plot(Ms_values, 20*np.log10(GM_lower), 'b-', linewidth=2)
ax.axhline(6, color='red', linestyle='--', alpha=0.7,
           label='GM = 6 dB (guideline)')
ax.fill_between(Ms_values, 20*np.log10(GM_lower), 0,
                alpha=0.1, color='blue')
ax.set_xlabel(r"$M_s = \|S\|_\infty$", fontsize=12)
ax.set_ylabel("Gain Margin lower bound [dB]", fontsize=12)
ax.set_title(r"$M_s$ vs Gain Margin", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim([1, 5])
ax.set_ylim([0, 30])

# (b) M_s vs 位相余裕
ax = axes[1]
ax.plot(Ms_values, PM_lower, 'r-', linewidth=2)
ax.axhline(30, color='blue', linestyle='--', alpha=0.7,
           label='PM = 30 deg (guideline)')
ax.fill_between(Ms_values, PM_lower, 0, alpha=0.1, color='red')
ax.set_xlabel(r"$M_s = \|S\|_\infty$", fontsize=12)
ax.set_ylabel("Phase Margin lower bound [deg]", fontsize=12)
ax.set_title(r"$M_s$ vs Phase Margin", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim([1, 5])
ax.set_ylim([0, 90])

plt.tight_layout()
plt.savefig("ms_stability_margins.png",
            dpi=150, bbox_inches="tight")
plt.show()

このグラフから、$M_s$ と安定余裕の関係が読み取れます。

左図（$M_s$ vs ゲイン余裕）: $M_s$ が小さいほどゲイン余裕の下界が大きくなります。$M_s = 2$（実用的な上限）のとき、ゲイン余裕は6 dB以上が保証されます。$M_s = 1.5$ であればゲイン余裕は約9.5 dB以上です
右図（$M_s$ vs 位相余裕）: 同様に、$M_s = 2$ のとき位相余裕は約29度以上が保証されます。$M_s \leq 2$ という設計指針が、古典制御の「ゲイン余裕6 dB以上、位相余裕30度以上」というガイドラインにほぼ対応していることがわかります

つまり、$M_s$ は古典制御のゲイン余裕・位相余裕を1つの指標で統一的に表現するものです。MIMO系ではゲイン余裕や位相余裕の定義が自明でないため、$M_s$ がより基本的な安定余裕の指標として使われます。

ロバスト性能条件

ロバスト安定性とロバスト性能の区別

ここまでに2つの条件を見てきました。

ノミナル性能（NP）: ノミナルモデルで性能仕様を満たす → $\|W_1 S\|_\infty < 1$
ロバスト安定性（RS）: 全ての不確実プラントで安定 → $\|W_m T\|_\infty < 1$

実際に必要なのは、これらの組み合わせであるロバスト性能（RP）です。

ロバスト性能: 不確実性の範囲内の全てのプラントに対して、性能仕様を満たす。

乗法的不確実性に対するロバスト性能条件は

$$ \begin{equation} |W_1(j\omega) S_\Delta(j\omega)| < 1, \quad \forall \|\Delta\|_\infty \leq 1, \forall \omega \end{equation} $$

ここで $S_\Delta$ は不確実プラントに対する感度関数です。

構造化特異値 $\mu$

ロバスト性能の正確な条件は構造化特異値 $\mu$（structured singular value; SSV）で記述されます。

$$ \mu_\Delta[M(j\omega)] < 1, \quad \forall \omega $$

$\mu$ はスモールゲイン定理の保守性を低減した指標で、不確実性の「構造」（どのパラメータが独立に変動するか）を考慮します。

$\mu$ の計算は一般にNP困難ですが、上界と下界の計算は効率的に行えます。上界は $D$ スケーリングと呼ばれる手法で計算されます。

$$ \mu_\Delta[M] \leq \inf_{D \in \mathcal{D}} \bar{\sigma}[DMD^{-1}] $$

$D$-$K$ イテレーションは、$\mu$ の上界を最小化する制御器を求める反復アルゴリズムです。

$D$ を固定して $K$ を $\mathcal{H}_\infty$ 最適化で求める
$K$ を固定して $D$ を周波数ごとに最適化し、有理関数で近似する
収束するまで1と2を繰り返す

この反復は $\mu$ の上界を単調に減少させますが、大域的最適解への収束は保証されません。

ロバスト安定条件の実用的なまとめ

不確実性の種類	ロバスト安定条件
乗法的不確実性	$\\|W_m T\\|_\infty < 1$
加法的不確実性	$\\|W_a CS\\|_\infty < 1$
逆乗法的不確実性	$\\|W_{im} S\\|_\infty < 1$
フィードバック不確実性	$\\|W_f CS\\|_\infty < 1$

いずれも、「不確実性が大きい周波数帯域で対応する閉ループ伝達関数を小さくする」という同じ構造を持っています。

ロバスト性能条件の理解が深まったところで、最後に設計手順の全体像をまとめます。

H∞制御設計の手順まとめ

H∞混合感度設計の実践的な手順は以下の通りです。

Step 1: プラントのモデリング

物理法則やシステム同定から、ノミナルプラント $G(s)$ を求めます。

Step 2: 不確実性の評価

パラメータ変動の範囲や未モデル化ダイナミクスの大きさを評価し、不確実性の重み $W_m(s)$（または $W_a(s)$）を設定します。

Step 3: 性能仕様の設定

所望の帯域幅、定常偏差、外乱抑制性能から、感度関数の重み $W_1(s)$ を設定します。

Step 4: 一般化プラントの構成

$G(s)$、$W_1(s)$、$W_2(s)$（$\geq W_m(s)$）から一般化プラント $P(s)$ を構成します。

Step 5: H∞最適化

$\|\mathcal{F}_l(P, K)\|_\infty$ を最小化する $K(s)$ を求めます。数値的にはリカッチ方程式のソルバやLMI（線形行列不等式）ソルバを使用します。

Step 6: 検証

ノミナル性能: $\|W_1 S\|_\infty < 1$ を確認
ロバスト安定性: $\|W_2 T\|_\infty < 1$ を確認
ステップ応答、感度関数のボード線図で時間領域・周波数領域の性能を確認
パラメータ変動を含むプラントでモンテカルロシミュレーション

Step 7: 制御器の低次元化（必要に応じて）

H∞最適制御器はプラントの次数と同程度の次数を持つため、実装のために平衡化打ち切りなどで低次元化が必要になることがあります。

import numpy as np
import matplotlib.pyplot as plt

# H∞設計手順のフローチャートを可視化
fig, ax = plt.subplots(1, 1, figsize=(10, 8))
ax.set_xlim([0, 10])
ax.set_ylim([0, 10])
ax.axis('off')

steps = [
    (5, 9.2, "Step 1: Plant Modeling\n$G(s)$"),
    (5, 8.0, "Step 2: Uncertainty Evaluation\n$W_m(s)$"),
    (5, 6.8, "Step 3: Performance Spec\n$W_1(s), W_2(s)$"),
    (5, 5.6, "Step 4: Generalized Plant\n$P(s)$"),
    (5, 4.4, "Step 5: $H_\\infty$ Optimization\n$\\min_K \\|\\mathcal{F}_l(P,K)\\|_\\infty$"),
    (5, 3.2, "Step 6: Verification\nNP + RS + RP check"),
    (5, 2.0, "Step 7: Controller Reduction\n(if needed)")
]

for x, y, text in steps:
    bbox = dict(boxstyle="round,pad=0.4", facecolor="lightblue",
                edgecolor="navy", linewidth=1.5)
    ax.text(x, y, text, ha='center', va='center',
            fontsize=11, bbox=bbox)

# 矢印
for i in range(len(steps) - 1):
    ax.annotate('', xy=(5, steps[i+1][1] + 0.45),
                xytext=(5, steps[i][1] - 0.45),
                arrowprops=dict(arrowstyle='->', color='navy',
                               linewidth=1.5))

# フィードバック矢印（Step 6からStep 3へ）
ax.annotate('Redesign\nweights',
            xy=(2.5, 6.8), xytext=(2.5, 3.2),
            arrowprops=dict(arrowstyle='->',
                           color='red', linewidth=1.5,
                           connectionstyle='arc3,rad=0.3'),
            fontsize=9, color='red', ha='center')

ax.set_title("$H_\\infty$ Mixed Sensitivity Design Procedure",
             fontsize=14, fontweight='bold', pad=10)
plt.tight_layout()
plt.savefig("h_infinity_design_procedure.png",
            dpi=150, bbox_inches="tight")
plt.show()

このフローチャートから、H∞設計が反復的なプロセスであることがわかります。Step 6の検証で仕様を満たさない場合は、Step 3に戻って重み関数を調整します。この反復は通常3-5回程度で収束し、経験を積むと重み関数の初期設定が上手くなり、反復回数が減ります。

LQG制御との比較

H∞制御の位置づけを明確にするために、LQG制御との比較を行います。

項目	LQG制御	H∞制御
最適化基準	$\mathcal{H}_2$ ノルム（平均性能）	$\mathcal{H}_\infty$ ノルム（最悪ケース）
外乱仮定	ホワイトガウスノイズ	有界エネルギーの任意信号
ロバスト性保証	なし（SISO系では60度の位相余裕）	スモールゲイン定理による保証
設計パラメータ	$Q, R$（状態・入力の重み行列）	$W_1, W_2$（周波数重み関数）
制御器の次数	プラントと同じ	プラントと同じ
分離定理	成立	一般には不成立

LQG制御は「モデルが正確で外乱がガウスノイズ」という理想的な状況で最適ですが、モデルが不正確なときに性能が劣化するという弱点があります（特にMIMO系でのLQGのロバスト性の不保証はDoyleの反例として知られています）。H∞制御はモデル誤差を明示的に考慮するため、ロバスト性が保証されますが、ノミナル性能はLQGに劣る場合があります。

実用的には、LQGとH∞のいいとこ取りとして、ループ整形（loop shaping）を行った後にH∞最適化を適用するMcFarlane-Glover法や、$\mathcal{H}_2$ と $\mathcal{H}_\infty$ の制約を同時に扱うマルチオブジェクティブ制御が使われることもあります。

まとめ

本記事では、ロバスト制御の基礎としてH∞制御の概念を解説しました。

モデル不確実性は制御設計で避けられない問題であり、加法的・乗法的不確実性として周波数領域でモデル化する
$\mathcal{H}_\infty$ ノルムは周波数応答のゲインの最大値であり、「最悪ケースでの入出力エネルギーゲイン」を表す
スモールゲイン定理はロバスト安定性の基本条件 $\|M\|_\infty < 1$ を与え、乗法的不確実性に対しては $\|W_m T\|_\infty < 1$ となる
混合感度問題で感度関数 $S$ と相補感度関数 $T$ の周波数重み付きノルムを最小化し、追従性能とロバスト安定性のトレードオフを設計する
一般化プラント表現とリカッチ方程式により、H∞最適制御器が系統的に求められる
$M_s = \|S\|_\infty$ は古典制御のゲイン余裕・位相余裕を統一的に表す指標であり、$M_s \leq 2$ が実用的なガイドラインである

次のステップとして、以下の記事も参考にしてください。

モデル予測制御 — 制約付き最適制御
適応制御の基礎 — パラメータ変動への別のアプローチ
LQG制御 — H2最適制御との対比
ゲイン余裕と位相余裕 — 古典的な安定余裕