はじめに
本記事はStormworks 第1 Advent Calendar 2024の24日目の記事ですが，扱うのは古典制御理論の概要であるため制御理論を始めたばかりの人でも読めるように作っているつもりです（ただしラプラス変換は前提知識）が，後半はかなり難しい話になるので半分備忘録になっています．
さて，本記事を読んでいるということは諸兄は「Stormworksでも実際の工業製品に使われているような高度な制御をやってみたい！」という希望を少なからず抱いていると思います．しかし，実際の産業界で使われている制御則のほとんどがPID制御の類です．ただ，少なくない数のコントローラがロバスト制御や適応制御，H∞制御といった発展的な制御則を用いて非常に良い性能を実現しているのも事実です．
そこで本記事はこれらの発展的な制御則を学ぶための準備の一環として，古典制御理論を大まかに紹介していきます．
目次
本記事は「はじめての制御工学」に倣い，次のような構成で大雑把に書いています．詳しく知りたい人は教科書を買った方がいいと思います．
1. 古典制御理論とは？
2. 伝達関数
3. 極と安定性
4. 制御系の構成とその安定性
5. 周波数特性の解析
6. ゲイン余裕と位相余裕
7. ループ整形法
1. 古典制御理論とは？
まず「古典制御理論とは何ぞや？」と思われる読者が多いかもしれないのでざっくりと説明しておくと，「１入力１出力で完結しているコントローラの設計理論」という認識で（少なくとも初学者は）いいと思います．このため基本的には線形なコントローラになるため解析が非常にしやすく，簡単にいい感じ（※要出典）のコントローラを設計することができます．しかし現実のシステムはパラメータが未知で時不変でなかったりいくつかの線形近似を行っていたりと，ゲイン余裕や位相余裕を考えないと痛い目を見る可能性があります．
2. 伝達関数
制御工学という分野では，コントローラを含むシステム全体を「ブロック線図」という図で表していることが多いです．今回はブロック線図については省略し，ブロック線図内にあるシステムそのものを表す「伝達関数」について説明します．これを用いることで，（制御屋さん視点で）直感的でわかりやすく表現することができます．
いま，入力が$U(s)$，出力が$Y(s)$で二次遅れ系（二階線形微分方程式の系をラプラス変換したもの）のシステムが
U(s)=(s^2+as+b)Y(s), a>0, b>0
のとき，伝達関数$G(s)$を
G(s)=\frac{Y(s)}{U(s)}=\frac{1}{s^2+as+b}
のように定義できます．これを式変形することで
Y(s)=\frac{1}{s^2+as+b}U(s)=G(s)U(s)
と表すことができます．この伝達関数で表されるシステムの定常応答（システムに定数入力）は最終値の定理より
\lim_{t\to \infty}{y(t)}=\lim_{s\to 0}{sG(s)}
で簡単に求めることができる（※収束が前提条件のため，発散する場合や振動入力では最終値の定理は使えません）ように，非常に便利なものとなっています．
3. 極と安定性
さて，前章で例として挙げた二次遅れ系の伝達関数$G(s)$に対して，その分母を特性多項式と言い，これと$0$を等式関係にしたものを特性方程式と言います．具体的には，
s^2+as+b=0
が特性方程式となりこれの$s$についての解
s=\frac{-a\pm\sqrt{a^2-4b}}{2}
をシステムの固有値あるいはシステムの極と言います．この極の実部が負であればシステムは安定（無限時間後に一定）であると言えます．ちなみに虚数部は絶対値が大きければ大きいほど大きな振動をします．
4. 制御系の構成とその安定性
ここからが本題となります．システムが既知で
\ddot{y}(t)+a\dot{y}(t)+by(t)=u(t),a>0,b>0
y(0)=1,\dot{y}(0)=0
で表現できるとします．これをラプラス変換すると，
Y(s)(s^2+as+b)=U(s)+s+a \tag{1}
となります．このシステムは静的安定性を持つため無限時間後には$y(t)$が$0$となります（最終値の定理より求まる）．このシステムに対し，制御目標を
\lim_{t\to \infty}{y(t)}=r(t)
となるようにしてみましょう．このとき，偏差を
e(t)=r(t)-y(t)
E(s)=R(s)-Y(s)
とおくことができます．これを用いて式(1)に$Y(s)=E(s)-R(s)$を代入すると
E(s)=\frac{1}{s^2+as+b}U(s)+\frac{s+a}{s^2+as+b}+\frac{r}{s} \tag{2}
と表現できます．さて，これに最終値の定理を適用すると$r$が出てくるため制御目標を達成できません．そこで，入力$U(s)$を
U(s)=-(s+a)-(k_p+\frac{k_i}{s}+k_ds)E(s) \tag{3}
k_p>0, k_i>0,k_d>0
とおきます．このとき，$-(s+a)$は式(2)の$\frac{s+a}{s^2+as+b}$を打ち消す項，その後ろはいわゆるPID制御のものと同じです．こうすることで式(3)は
E(s)=\frac{s^3+as^2+bs}{s^3+(a+k_d)s^2+(b+k_p)s+k_i}R(s) \tag{4}
と表すことができます．さて，一般に$r(t)$は定数であるため$R(s)=\frac{r}{s}$となるから，式(4)は
E(s)=\frac{r(s^2+as+b)}{s^3+(a+k_d)s^2+(b+k_p)s+k_i} \tag{5}
と表すことができますね．ここで，最終値の定理より，
\lim_{t\to \infty}{e(t)}=\lim_{s\to 0}{sE(s)}=0
となることから，制御目標を達成できています．しかし，Stormworksをやっている皆さんはお分かりだと思いますがゲインチューニングしないとこの入力はまともに使えません．
そこで，前章での特性方程式を使います．式(5)より，偏差$E(s)$の特性方程式は
s^3+(a+k_d)s^2+(b+k_p)s+k_i=0 \tag{6}
となります．これの解は一般に３つの複素数（その中の特殊なものとして実数解があることも）となり，PIDの各パラメータを変えることによってある程度極を任意のものにすることができます．
余談ですが，PIDの各ゲインを正に限定したのは，ラウスの安定判別法やフルビッツの安定判別法の第1条件として
・特性方程式の全ての係数の符号が同じ
を確実に満たすためです．なのでこれを満たすのであればゲインが多少負であっても一応大丈夫です（解析的には）．
5. 周波数特性の解析
前章で最終値の定理というすこぶる便利なものを何の前触れもなく使い倒しましたが，この定理は万能ではなくシステムが収束するときのみ使用可能という制限があります．しかし実際のシステムは定数入力だけでなく周波数入力も発生し，このときの安定性も非常に重要になってきます．
さて，以下のようなシステム
U(s)=(s^2+as+b)Y(s), a>0, b>0
G(s)=\frac{Y(s)}{U(s)}=\frac{1}{s^2+as+b}
に対して，正弦波入力を考えていきます．ラプラス変換でももちろん解くことができますが，いろいろと面倒くさいです．しかも制御として考える場合，入力の周波数はノイズであり大体未知なので，虱潰しで探していくことになります．いくらラプラス変換が便利でもやってられません．そこで，
g=20\log_{10}{|G(jω)|}=20\log_{10}{\frac{B}{A}} \tag{7}
φ=\angle G(jω) \tag{8}
を考えます．ここで$A$は入力の振幅，$B$は出力の振幅となるため$g$はゲイン，$φ$は位相の遅れを意味します．ちなみに$n$次遅れ系は最大で$180n$°まで位相が遅れると決まっています．
これを用いることで，周波数入力に対する安定解析が非常に楽になります．ここでは紹介しませんが，折れ線近似で結構簡単に手書きできる方法もあるくらいです．今回は具体例として
G(s)=\frac{Y(s)}{U(s)}=\frac{1}{s^2+s+4}
の周波数解析をしていきましょう．このとき，
G(jω)=\frac{4-ω^2-jω}{(4-ω^2)^2+ω^2} \tag{9}
これをMATLABを用いてボード線図にすると，次のようになります．
グラフの横軸が入力の周波数，上のグラフの縦軸がゲイン$g$，下のグラフの縦軸が位相$φ$を表しています．いま，システムは二次遅れ系なので位相が最大で180°遅れており，90°ではゲインが局所的に大きくなっています．これが共振現象ですね．このように，ボード線図を用いるとシステムの周波数特性を視覚的に分かりやすくできます．
6. ゲイン余裕と位相余裕
一般に，システムの伝達関数を$P(s)$，コントローラの伝達関数を$C(s)$とした場合，システム全体の伝達関数を閉ループ伝達関数と呼び，
G(s)=\frac{P(s)C(s)}{1+P(s)C(s)} \tag{10}
と表現されます．このとき，（分母にあたる）特性多項式を用いて
1+P(s)C(s)=1+L(s) \tag{11}
の開ループ伝達関数$L(s)$を定義します．この開ループ伝達関数を周波数解析することで，望ましい応答なのかどうかを判別することができます．逆を言えば，この伝達関数を設計することで望ましい応答のコントローラを実現できます．この手法はループ整形法といい，次章で説明します．
さて，以下のようなシステムの開ループ伝達関数を考えるとします．
P(s)=\frac{1}{(s+1)^2}, C(s)=\frac{0.3}{s+0.1} \tag{12}
このシステムのボード線図は下のようになります．
ここで，
・ゲイン交差周波数$ω_{gc}$：$|L(jω_{gc})|=1$となる$ω$
・位相交差周波数$ω_{pc}$：$\angle G(jω_{pc})=-180°$となる$ω$
が定義でき，このときの$PM≈80°$，$GM≈20$dBがそれぞれ位相余裕，ゲイン余裕となります．ここで，$PM$はシステムがどのくらいの余裕をもって安定かを意味しており，$GM$は開ループ伝達関数$L(s)$をどこまで大きくしても安定かを意味しています．
7. ループ整形法
本章では簡単のため具体例のみ示します．いま，制御対象が
P(s)=\frac{0.01}{s^2+0.04s+0.01} \tag{13}
であるとき，ボード線図は次のようになります．
ここで，位相遅れコントローラを
C(s)=\frac{100(s+0.1)}{s}\frac{16(s+0.5)}{s+8} \tag{14}
とおいたとき，開ループ伝達関数$L_1(s)=P(s)C(s)$のボード線図は下のようになります．
このボード線図で注目すべき点として，ゲイン交差周波数$ω_{gc}$が高いため即応性が高く，低周波領域で$-20$dB/decとなっていることから定常偏差は$0$となります．また，位相余裕PMが確保されていることから応答も振動的でないことがわかります．
さて，ここでロールオフ特性（観測値に含まれるノイズを効率よく低減したい）を考える．$ω_s=25$rad/sでゲインを速やかに減衰させるため，$ω_{s0}=10$rad/s以上の帯域でゲインを減衰することにします．そこで新たなコントローラを
C_r(s)=C(s)\frac{{ω_{s0}}^2}{s^2+ω_{s0}s+{ω_{s0}}^2} \tag{15}
とします．これを用いた開ループ伝達関数は
となります．このコントローラを用いることにより，$ω_{s0}=10$rad/s以上の観測ノイズの影響を大きく低減することができます．
おわりに
いかがでしたか？ こんなもので理解できたなら苦労はしないと思います．実際証明は全て省略していますし何より具体例がほとんどありません．しかし，この記事を最後まで読んである程度理解できた方なら「システムさえわかれば簡単じゃん」とお思いでしょう．ですが残念なことに，Stormworksだけでなく現実でも既知なシステムなんて一握りですし，そのパラメータも未知となるのがほとんどです．
しかし，近年の研究では既知なパラメータを一つも用いない自動運転コントローラが登場してきているように，制御理論も高度化してできることも増え始めています．そうしたコントローラの開発は研究レベルとなるため専門家でないと理解が難しいですが，この記事でより高度な制御則が存在する，ということだけでも知って制御の勉強のモチベーションにしていただければ幸いです．