実際に使える制御工学：基礎編（３）「周波数特性とボード線図。SIN波でシステムの素性を暴く。」
＜この記事で伝えたい事＞
周波数特性とは、入力の周波数を変えると、出力の振幅や位相がどうなるかを表すものです
システムへの入力はSIN波（正弦波）なことが前提です。ボード線図とは、周波数特性をグラフ化する方法の1つです
ナイキスト線図というのもありますが、個人的にはボード線図の方が見やすいと思います。
手描きするのは大変ですが、コンピュータの力を借りれば楽勝です。
（使うソフトウェアはMATLABやScilab）
1. はじめに
自己紹介はこちらです。
この記事は、主に制御系の学生や若手エンジニアに向けた基礎知識紹介・解説記事シリーズの第３弾です。
今回は、制御工学でよく使う概念・ツールである周波数特性とボード線図を紹介・解説します。
2. 周波数特性とは？
冒頭で述べたとおりですが、文章で述べたことを図で表現してみます。
SIN波の周波数を変えると、振幅の小さくなる度合い・タイミングの遅れ度合いも変わります。
こういう性質・特徴をまとめて周波数特性と呼ぶのだと理解してもらえば、だいだい合っていると思います。
これだけ聞くと「その話、SIN波でしか成り立たないんじゃない？」と思う方が多いでしょう。
しかし、安心してください。
なぜなら、フーリエ変換によって、あらゆる波形は異なる周波数・振幅・位相をもったSIN波の組み合わせであることが数学的に証明されています。
（極一部の例外はありますが実用上は問題ないです、興味があればフーリエ変換について調べてみてください）
こうしたわけで、周波数特性の話をするときはSIN波を軸にするのが合理的・効率的なのです。
3. MATLABを使ったボード線図の描き方
これも冒頭で述べたとおりですが、ボード線図とは周波数特性をグラフ化する方法の１つです。
システム（伝達関数）の周波数特性は数式だけでも表現できますが、グラフ化した方が理解しやすいです。
ちょうど、物体の運動を数式だけで表現するよりグラフ化した方が理解しやすいのと似ています。
しかし、このボード線図を手描きするには大変な手間がかかります。
加えて、パラメータを1つ変えるだけで最初からやり直しです。
とてもやってられません。
そこで、MATLABの出番です。
コンピュータの力を借りましょう。
基礎編（２）で求めたDCモータ全体の伝達関数を具体例に、MATLABを使ったボード線図の描き方を紹介します。
DCモータの伝達関数
$$
\frac{\omega_{(s)}}{e_{(s)}}=\frac{K_t}{J\cdot L\cdot s^2+J\cdot R\cdot s+K_t\cdot K_e}
$$
$$
\begin{array}{cl}
e & ：\text{モータにかける電圧} \\
\omega & ：\text{モータシャフトの角速度}
\end{array}
$$
MATLABで以下のスクリプトを実行すればボード線図を描いてくれます。
コンピュータのやることなので具体的な数値が必要です。
% パラメータの定義
Kt=0.01; Ke=0.01; R=1.0; L=0.5e-3; J=3e-4;
% Kt:トルク定数[Nm], Ke:逆起電圧定数[V/(rad/s)], R:電気抵抗[Ω]
% L:コイルインダクタンス[H], J:モータ慣性モーメント[kgm^2]
% 伝達関数の定義
s=tf(‘s’); % MATLABに変数sをラプラス変数として認識させる
Gs=Kt/(J*L*s^2 + J*R*s + Kt*Ke);
% ボード線図の描画
opts=bodeoptions; % ボード線図のオプション設定のためのおまじない
opts.FreqUnits=’Hz’; % ボード線図の横軸単位をHzに指定(標準の横軸単位はrad/s)
bode(Gs,opts);
grid on;おそらくコピー&ペーストで問題なく動作します。
ただ、’（シングルクォーテーション）が半角でないとエラーが出ます。
そこだけ注意してください。
4. ボード線図の読み方
3章のスクリプトを実行すると下図のボード線図が描かれます。
上側はゲイン線図、下側は位相線図と呼ばれます。
この２つがセットでボード線図です。
この図を見て最初に戸惑うことは
「ゲインの単位dB(デシベル)って何だろう？」
という点ではないでしょうか？
私の推測ですが、ゲインの単位に[dB]を使うのはボード線図の考案者であるヘンドリック・W・ボーディ博士が電子・電気通信分野の人でもあったことの名残りと思います。
電子・電気通信分野では増幅率の単位に好んで[dB]を使う文化があります。
（ボーディ博士はあの有名なアメリカのベル研究所で働いていました）
さて、本題に戻ります。
まず、ゲインというのは増幅率のことです。
例えば、ゲインが5のシステム（伝達関数）は入力を５倍して出力します。
この5倍というゲインが[dB]換算でいくつかは次の変換式で分かります。
$$
20\cdot log_{10}(5)=13.9794_{[dB]}
$$
逆に、[dB]で表されたゲインを定数に変換するには次の変換式を使います。
$$
10^{\frac{13.9794[dB]}{20}}=5
$$
対数という普段あまり使わないだろう概念に抵抗感あるかもしれませんが、ここでしか使わないのでボード線図を見るときだけガマンしてやってください。
[dB]の意味が分かったところでボード線図の読み方の説明に移ります。
ボード線図を読むときは
「このシステムに何[Hz]のSIN波を入力したとき、
何倍増幅されて出力されるのか？
何[deg]位相が遅れて出力されるのか？」
という目線で見ると情報を読み取りやすいです。
（degは角度の°と同じ意味です）
区切りが良い10[Hz]でこの読み方をしてみましょう。
まず、ゲイン線図に目をやると周波数10[Hz]でのゲインは-5.55[dB]です。
$${10^{\frac{-5.55[dB]}{20}}=0.5278}$$なので、周波数10[Hz]のSIN波を入力すると振幅が0.5278倍されて出力されると分かります。
（ちなみに0[dB]はゲイン1倍=等倍出力を意味します、デシベル値がプラスのときは増幅・マイナスのときは減衰を意味します）
次に、位相線図に目をやると周波数10[Hz]での位相は-91.5[dB]です。
周波数10[Hz]のSIN波を入力すると、入力波形に対して91.5[deg]遅れたSIN波が出力されると分かります。
（位相のマイナスは遅れ・プラスは進みを意味します）
総合的に見ると具体例のDCモータは10[Hz]のSIN波を入力したとき、
入力波形の振幅を0.5278倍して出力する
入力波形に対して出力波形の位相が91.5[deg]遅れる
という動きをすると言えます。
別の周波数でも同じ要領で情報を読み取れます。
視覚的に分かりやすいようDCモータに10[Hz]・振幅1のSIN波を入力したときのシミュレーション波形を下図に示します。
（システムの動きが十分落ち着いた10[s]以降の波形です）
ボード線図どおり、出力波形は振幅0.5278で入力波形より91.5[deg]位相が遅れる（時間換算で0.0254[s]遅れる）様子が分かります。
ただし、動きが落ち着く前（過渡状態）だと少し話が違ってきます。
そのあたりは別の機会に話せると思います。
5. 周波数特性とモノの動きとの対応付け
純粋な制御技術の観点からは蛇足な話かもしれませんが、せっかく具体例にDCモータを引き合いに出したので、この話もしようと思います。
（興味がなければ読み飛ばしてOKです）
まず、皆さんの人差し指を立てて、ゆっくり左右に動かしてみてください。
（メトロノーム風に指先を振るというより、反復横跳びのように手首より先の全体を振るイメージ）
そこから、意識的にだんだん指を振る速さを上げてみてください。
やってくれましたか？
実際やると指を振る速さを上げるほど思いどおりに動かなくなり、最終的にほとんど動かせなかったのではないでしょうか？
そりゃそうでしょう。
物体というのは質量や力が同じでも速い動きはしづらいですから。
これは人間の体に限らず、万物に共通することです。（私の知る限り）
これを踏まえ、DCモータのボード線図をもう一度見てみましょう。
周波数が上がるほど、
ゲインは右肩下がりに下がっていく
＝かける電圧の振幅が同じでも、周波数が上がればシステム出力である角速度は小さくなっていく位相遅れがだんだん大きくなっていく
＝電圧と角速度とのピークのタイミングのズレが大きくなっていく
という様子が分かると思います。
DCモータにかける電圧の周波数を上げることは、人間が指を振る速さを上げることと同じです。
DCモータも物体である以上、速すぎる動かし方にはついてこれません。
その性質が周波数上昇に伴うゲイン低下・位相遅れ増大として、ボード線図にも表れるのです。
このように、ボード線図を通して様々なモノの動きを日常生活の中で体験できる動きと重ね合わせて理解できます。
そうすれば、モノの動きをより直感的に・深く理解できるのではないかと思います。
さて、これまではちょっと捻った話をしました。
次は、わりと正統派の話をしようと思います。
それは、パラメータを変えたときボード線図はどう変わるか？の話です。
ここまで、DCモータのパラメータは3章のMATLABサンプルコードの数値を使っています。
例えば、モータの慣性モーメントは$${J=3.0×10^{-4}[kgm^2]}$$です。
この$${J}$$だけ10倍（$${30.0×10^{-4}[kgm^2]}$$）にするとボード線図にどのような変化が出るかを見てみましょう。
$${J}$$を10倍にすると全体的にゲインが低下する様子が確認できます。
今回の場合、システムのゲイン低下は同じ電圧に対して角速度が小さくなる（＝動きにくくなる）ことを意味します。
これは「慣性モーメント$${J}$$を大きくするほど、そのモノは動きにくくなる」という直感的理解とマッチする結果です。
このように、ボード線図はモノやシステムの物理的特性（動特性）を見える化する道具でもあります。
例えば、慣性モーメント$${J}$$だけでなく他のパラメータも変えて最適な組み合わせ（またはそれに近い落としどころ）を探り当てる使い方もあります。
6. おわりに
この記事では、周波数特性とは何なのか、それをグラフ化したものであるボード線図の描き方・読み方を紹介・解説しました。
周波数特性やボード線図というのは、制御技術を学ぶ上で誰もがつまずきやすいところです。
かくいう私自身もそうでした。
（一部の天才と呼ばれる人々はどうか知りませんが）
ただし、これらを理解できれば制御技術への解像度がグッと上がります。
一人前の制御技術者への第一歩といったところでしょうか。
なお、今回の記事はかつて某企業で若手・新人教育をしたときの内容をベースにブラッシュアップしたものです。
（そのときの受講者の反応はかなり良かったです）
基礎編のメインターゲット層である学生・若手エンジニアの皆さんの助けになれば幸いです。
前回：基礎編（２）「伝達関数とブロック線図」はこちら
次回：基礎編（４）「PID制御」はこちら
付録1. 手計算でのボード線図の描き方
本編中では、ボード線図を描くのはMATLABにお任せしました。
ここでは、MATLABが裏でどんな計算をしているのかを紹介します。
我々に見えないところでコンピュータがどれだけ大量の計算をやっているのかを知れば、その有難みが分かることでしょう。
まず、下式の伝達関数のボード線図を描きたいと仮定します。
$$
G_{(s)}=\frac{\alpha}{s+\alpha}
$$
ボード線図を描くのは$${G_{(s)}}$$に$${s=j\cdot \omega}$$を代入することから始まります。
$${j}$$は虚数です。
あの2乗すると$${-1}$$になるという不可解な数です。
$${\omega}$$は角周波数（単位は[rad/s]）と呼ばれるもので、ボード線図の横軸になっていた周波数$${f}$$の親戚のようなものです。
$${\omega}$$と$${f}$$との間には次の関係性があります。（$${\pi}$$は円周率）
$$
\omega_{[rad/s]}=2\cdot\pi\cdot f_{[Hz]}
$$
$${s=j\cdot \omega}$$を代入すると$${G_{(s)}}$$は次のように変わります。
$$
G_{(j\cdot\omega)}=\frac{\alpha}{j\cdot\omega+\alpha}
$$
パッと見で分かるように分母が複素数（＝実数＋虚数）になっています。
このままでは何が何だか分かりません。
せめて分母だけでも実数化（分母の有理化）して、少しでも理解に近づくのを試みます。
$$
\begin{array}{ccl}
G_{(j\cdot\omega)} & = & \cfrac{\alpha}{j\cdot\omega+\alpha}\cdot\cfrac{-j\cdot\omega+\alpha}{-j\cdot\omega+\alpha}\\
& = &\cfrac{\alpha^2-j\cdot\alpha\omega}{\alpha^2+\omega^2}\\
& = &\cfrac{\alpha^2}{\alpha^2+\omega^2}-j\cdot\cfrac{\alpha\cdot\omega}{\alpha^2+\omega^2}
\end{array}
$$
分子には虚数が残っているものの、分母は実数だけにできました。
ここで、ちょっとした発想の転換が必要です。
複素数は平面グラフ上に表現できることを思い出してください。
例えば、複素数$${2+j\cdot 3}$$は下図で表現できます。
これを$${G_{(j\cdot\omega)}}$$にあてはめてグラフを描いたものが下図です。
なんとなく$${G_{(j\cdot\omega)}}$$をイメージできてきたのではないでしょうか？
実は、このグラフにおいて、
原点$${O}$$と点$${\bm{G_{(j\cdot\omega)}}}$$との距離＝伝達関数のゲイン
実数軸（横軸）と線分$${\bm{O\text{-}G_{(j\cdot\omega)}}}$$のなす角＝伝達関数の位相（グラフ中の角度はマイナス）
なのです。
伝達関数のゲインを$${|G_{j\cdot\omega)}|}$$、位相を$${\angle G_{(j\cdot\omega)}}$$とします。
三平方の定理から$${|G_{j\cdot\omega)}|}$$は次のように求められます。
$$
|G_{j\cdot\omega)}|=\sqrt{
\begin{pmatrix}
\cfrac{\alpha^2}{\alpha^2+\omega^2}\end{pmatrix}^2
+
\begin{pmatrix}
\cfrac{-\alpha\cdot\omega}{\alpha^2+\omega^2}\end{pmatrix}^2}
$$
これは[dB]換算前の数値です。
本文中の説明どおり、dB換算するには$${log_{10}}$$をとって20倍する必要があります。
また、逆三角関数計算によって$${\angle G_{(j\cdot\omega)}}$$は次のように求められます。
$$
\begin{array}{ccl}
\angle G_{(j\cdot\omega)} & = & tan^{-1}
\begin{pmatrix}
\cfrac{(-\alpha\cdot\omega)/(\alpha^2+\omega^2)}{\alpha^2/(\alpha^2+\omega^2)}
\end{pmatrix}\\
& = &tan^{-1}
\begin{pmatrix}
\cfrac{-\alpha\cdot\omega}{\alpha^2}
\end{pmatrix}
\end{array}
$$
$$
\begin{pmatrix}
tan\angle G_{(j\cdot\omega)}=\cfrac{(-\alpha\cdot\omega)/(\alpha^2+\omega^2)}{\alpha^2/(\alpha^2+\omega^2)}
\end{pmatrix}
$$
そして、$${|G_{j\cdot\omega)}|}$$と$${\angle G_{(j\cdot\omega)}}$$は共に$${\omega}$$によって変わるのを忘れてはいけません。
ボード線図を手描きするには$${\omega}$$を変えて$${20\cdot log_{10}(|G_{j\cdot\omega)}|)}$$と$${\angle G_{(j\cdot\omega)}}$$を計算し、点をいくつもプロットする必要があります。
最終的に、プロットした点を線でつなぎ合わせてボード線図が完成します。
点1つをプロットするだけでそこそこの量の計算が必要であり、ボード線図を描き終わるには膨大とも言える量の計算が必要だと想像できたのではないでしょうか？
通常、コンピュータは数千～数万とおりの$${\omega}$$について計算を実行してボード線図を描画します。
人間がこれをやれば、1つのボード線図を描き終えるのに数日～数週間かかってもおかしくはありません。
そこをコンピュータはものの数秒でやってくれます。
得意なことはコンピュータに任せてしまいましょう。
そうすれば、我々人間はコンピュータには不可能な高度な思考（設計判断など）に力を注ぎ込めるので合理的です。
付録2. ボード線図とナイキスト線図との対比
記事冒頭で「周波数特性のグラフ化方法にはナイキスト線図というのもある」と事を言いました。
このナイキスト線図が何なのかを説明します。
ナイキスト線図とは、付録1で説明した複素数平面グラフ上にを変えたを何点もプロットして1つの線でつなぎ合わせたものです。
本質的にナイキスト線図とボード線図は同じものです。
2つの違いは見せ方・まとめ方だけです。
具体例を図示します。
本質的に同じでもパッと見の印象はだいぶ違って見えるでしょう。
ボード線図のゲインがdB換算値なのに対し、ナイキスト線図はdB換算前のゲインなのもそう見える原因の1つと思います。
しかし、2つが違って見える根本的な原因は別のところにあります。
それは周波数特性というものが、
周波数
ゲイン
位相
という3つの要素から成り立っている点です。
要素が2つなら、2次元平面のグラフ上に問題なく情報を表現できます。
しかし、3つの要素を1つのグラフに無理なく情報をまとめるのは困難です。
（2次元グラフは平面＝2Dなので）
ボード線図はグラフをゲイン線図と位相線図の2つに分けることで、この問題をクリアしています。
一方、ナイキスト線図は3つの要素を1つのグラフにまとめています。
このため、注釈がなければどの周波数でゲイン・位相がいくつかを読み取れません。
注釈が不要なボード線図の方がよりスマートな方法と言えます。
周波数特性を表すのにボード線図の方がよく使われる理由はここにあると思います。
しかし、ナイキスト線図には特有のメリットがあります。
それはナイキストの安定判別法との相性が非常に良いことです。
制御システムの安定/不安定は制御設計における最重要項目です。
これが一目で分かるだけでもナイキスト線図は価値ある作図方法と言えるでしょう。
（ナイキストの安定判別法は別記事で紹介したいと思います）
付録3. 制御工学とボード線図の相性の良さ
基礎編（２）の2章で「ブロック線図の直列接続は伝達関数のかけ算」という話をしました。
制御工学をやると頻繁に伝達関数のかけ算をすることになります。
仮に、定数ゲインの伝達関数$${G_{1(s)},G_{2(s)}}$$があるとしましょう。
これらのかけ算は$${G_{1(s)}\cdot G_{2(s)}=2 \times 3=6}$$です。
これだけで済めばとても簡単です。
しかし、だいたいの場合、伝達関数ゲインは周波数によって変わります。
周波数を変えながらゲインを計算してかけ算するのもできなくはないですが、なかなか手間がかかります。
そこで出てくるのが対数（ゲインのdB換算時に使う計算方法）です。
試しに、2と3のdB換算値を計算してみましょう。
$$
20\cdot log_{10}(2)=6.0206_{[dB]}\\
20\cdot log_{10}(3)=9.5424_{[dB]}\\
$$
ここで、2と3のかけ算結果である6のdB換算値を計算して上2式と見比べてみましょう。
$$
\begin{array}{cl}
20\cdot log_{10}(6)&=15.5630_{[dB]}\\
&=6.0206_{[dB]}+9.5424_{[dB]}
\end{array}
$$
こうしてみると分かるようにdB換算後の数値を足すと、かけ算結果のdB換算値と一致します。
この性質は数学公式の積の対数として知られています。
$$
log_a(M\cdot N)=log_a(M)+log_a(N)
$$
これは利用できる性質です。
$${G_{1(s)},G_{2(s)}}$$のゲイン線図を描いた後、それらを足し算するだけで$${G_{1(s)}\cdot G_{2(s)}}$$のゲイン線図が求められるからです。
位相線図については対数を持ち出すまでもありません。
仮に、$${G_{1(s)},G_{2(s)}}$$がそれぞれ-90[deg]の位相を持つとします。
まず、$${G_{1(s)}}$$を通ると-90[deg]位相が遅れます。
それが$${G_{2(s)}}$$を通ると更に-90[deg]遅れます。
このため、全体ではそれらを足し算した-180[deg]の位相遅れになります。
$${G_{1(s)}=\cfrac{6}{s+4},G_{2(s)}=\cfrac{8}{2s+3}}$$としたときの$${G_{1(s)},G_{2(s)}}$$と、それらのかけ算$${G_{1(s)}\cdot G_{2(s)}}$$のボード線図を下図に示します。
コンピュータの普及した現在では意識することは少ないですが、昔はこういう知識も大事だったのだろうと思います。
