制御系の過渡応答・周波数特性解析
目次
1.アナログコンピュータによる周波数特性解析
電子式アナログコンピュータによって、伝達関数 \(G(j\omega)\) で表されるシステムの周波数応答解析を行うことができる。その方式として第1に挙げられるのは、straightforwardに開ループでプログラムし、振幅と位相についてボード線図を描出する方法である。振幅については、絶対値要素と積分器（LPF）による包絡線検波ブロック、対数関数演算回路など、基本的な演算要素で構成可能であり、位相特性については \(\arctan{(y/x)}\) を出力する演算回路を乗算器で構成したり、比較器と積分器を用いる参照信号との位相比較回路をプログラムするなどして実現可能であるが、積分時間など、様々なパラメータの調整が難点である。
上記方法による周波数応答解析は参考書籍 [1] に譲るとして、今回は文献 [2], [3], [4]を参考に、再急降下法を利用して時不変システムの周波数応答を出力するアナログコンピュータプログラムを作成し、XYレコーダをもってそのベクトル軌跡（つまり\(G(j\omega)\)の複素数平面上における軌跡を \(\omega=0\rightarrow\infty\) の範囲で描出したもの ）を作成することを目標とする。
2.アナログコンピュータプログラム
2.1 演算の原理
解析対象のシステムの伝達関数を \(G(j\omega)\) とし、入力参照信号を
\[x(t)=a\,\sin{(\omega\,t)}\tag{1}\]
とし、観測される出力信号を
\[y(t)=b(\omega)\,\sin{(\omega\,t)}+c(\omega)\,\cos{(\omega\,t)}\tag{2}\]
する。\(a>0\) は入力の振幅、\(b>0\) は出力信号の入力に対して同相な成分の振幅、\(c>0\) は出力信号の入力に対し直交する成分の振幅である。通常の方法では、入出力信号の振幅比および位相シフトを取得する演算回路を組んで、ゲインと位相についてボード線図などを描出するわけだが、ここでは直交形式で得ることを目標とし、その際の成分分離には再急降下法を用いる。すなわち、目標値 \(\hat{y}(t)\) と実測値 \(y(t)\) の差を示す誤差関数
\[\varepsilon(t)=y(t)-\hat{y}(t)=y(t)-\hat{b}(\omega)\,\sin{(\omega\,t)}-\hat{c}(\omega)\,\cos{(\omega\,t)}\tag{3}\]
を定め、\(f=\frac{1}{2}\varepsilon^2\) を評価関数として設定し、正のゲイン \(K_b,\,K_c\) に対して
\begin{eqnarray}\begin{cases}\displaystyle{\frac{\mathrm{d}\hat{b}}{\mathrm{d}t}=-K_b\,\frac{\partial{f}}{\partial{\hat{b}}}}=K_b\,\varepsilon(t)\sin{(\omega\,t)}\\\\\displaystyle{\frac{\mathrm{d}\hat{c}}{\mathrm{d}t}=-K_c\,\frac{\partial{f}}{\partial{\hat{c}}}}=K_c\,\varepsilon(t)\cos{(\omega\,t)}\end{cases}\tag{4}\end{eqnarray}
をプログラムする。これにより \(\hat{y}(t)\) がすみやかに \(y(t)\) の直交2成分に吸着するように動き、\(x(t)\) の角周波数 \(\omega\) を変えていけば、 \(\hat{b}\)-\(\hat{c}\)平面上でベクトル軌跡を描く。
ここで、\(u(t)=a\sin{\omega\,t}\) と \(v(t)=a\cos{\omega\,t}\) を生成する回路は別に用意しなくてはならないのだが、通常の角周波数可変型の調和振動型微分方程式ブロック
\begin{eqnarray}\begin{cases}\displaystyle{\frac{\mathrm{d}u(t)}{\mathrm{d}t}=\omega\,v(t)}\\\\\displaystyle{\frac{\mathrm{d}v(t)}{\mathrm{d}t}=-\omega\,u(t)}\end{cases}\tag{5}\end{eqnarray}
を \(v(0)=a\,,u(0)=0\) の初期条件のもと解いた解を利用するだけでは、時間がたつとともに振幅変動誤差が生じうる（サークルテストの回路と同じ）。ゆえに、\(a^2\sin^2{\omega\,t}+a^2\cos^2{\omega\,t}\) の和を一定に保つようなサーボ系を付す工夫を行う。 具体的には、符号付き誤差
\[\epsilon(t)=a^2-\bigl\{u(t)^2+v(t)^2\bigr\}=a^2-r(t)^2\tag{6}\]
について、(5)式に
\begin{eqnarray}\begin{cases}\displaystyle{\frac{\mathrm{d}u(t)}{\mathrm{d}t}=\omega\,v(t)+K_1\,\epsilon(t)\,u(t)}\\\\\displaystyle{\frac{\mathrm{d}v(t)}{\mathrm{d}t}=-\omega\,u(t)+K_1\,\epsilon(t)\,v(t)}\end{cases}\tag{7}\end{eqnarray}
というように比例制御項を導入することになる（\(K_1\)は正の定数）。これが振幅の安定化に貢献する理由は以下のように考えれば理解される。
\(r^2(t)=u(t)^2+v(t)^2\) であるから、
\[\frac{\mathrm{d}r^2}{\mathrm{d}t}=2\,u\,\frac{\mathrm{d}u}{\mathrm{d}t}+2\,v\,\frac{\mathrm{d}v}{\mathrm{d}t}\tag{8}\]
(8)式に微分方程式(7)を代入すると
\[\frac{\mathrm{d}r^2}{\mathrm{d}t}=2\,u\,\biggr(\omega\,v+K_1\,\epsilon\,u\biggr)+2\,v\,\biggr(-\omega\,v+K_1\,\epsilon\,v\biggr)\tag{9}\]
回転項が打ち消されて
\begin{eqnarray}\frac{\mathrm{d}r^2}{\mathrm{d}t}&=&2\,u^2\,K_1\,\epsilon+2\,v^2\,K_1\,\epsilon\\\\&=&2\,K_1\,\epsilon\,r^2\tag{10}\end{eqnarray}
と整理される。\(r\neq{0}\) の場合、この式は、現状の振幅 \(r(t)\) が目標値 \(a\) より小さい場合（すなわち \(\epsilon>0\) のとき）は \(u-v\)平面において \(\frac{\mathrm{d}r^2}{\mathrm{d}t}>0\) の方向に軌道が移動し、発振振幅が大きくなることを示す。逆に現状の振幅 \(r(t)\) が目標値 \(a\) より大きい場合（すなわち \(\epsilon<0\) のとき）は振幅が小さくなることがわかる。以上から、(7)式の相平面上解軌道は半径 \(a\) の円に吸着し、振幅安定化の効果を有する線形発振方程式となっている。
ただし、ここでは(7)式の2式のうち一方だけに振幅安定化項を付加するだけとし、演算器の節約を図るものとする。たとえば
\begin{eqnarray}\begin{cases}\displaystyle{\frac{\mathrm{d}u(t)}{\mathrm{d}t}=\omega\,v(t)}\\\\\displaystyle{\frac{\mathrm{d}v(t)}{\mathrm{d}t}=-\omega\,u(t)+K_1\,\epsilon(t)\,v(t)}\end{cases}\tag{11}\end{eqnarray}
とすれば、
\begin{eqnarray}\frac{\mathrm{d}r^2}{\mathrm{d}t}=2\,K_1\,\epsilon\,v^2\tag{12}\end{eqnarray}
となり、\(v\neq{0}\) の場合において \(v\) 軸方向の吸着効果がはたらく。実際、初期状態が発振状態であれば \(v=0\) となるのはごく一瞬であるので、積分器の誤差による振幅減衰・振幅過大の累積誤差を修正する程度であれば(11)式のみで十分実用的である。
※ 演算器を節約したい場合は、より簡易的に \(V_z=\) 10 V~11 V 程度のツェナーダイオードによる振幅制限回路を使用してもよいだろう。
2.2 演算方程式・プログラム
(3), (4), (6), (11)を原方程式とする。\(K_b\)、\(K_c\) はともに共通値 \(K_0\) に設定し、振幅目標値を \(a=1.0\) とすれば、演算方程式群は以下の通りとなる。
\begin{eqnarray}\begin{cases}\displaystyle{a_\tau\,\frac{\mathrm{d}\hat{B}}{\mathrm{d}\tau}=-K_0\,E' \,U}\\\\\displaystyle{a_\tau\,\frac{\mathrm{d}\hat{C}}{\mathrm{d}\tau}=K_0\,E'\,V}\\\\\displaystyle{a_\tau\,\frac{\mathrm{d}U}{\mathrm{d}\tau}=\omega\,V}\\\\\displaystyle{a_\tau\frac{\mathrm{d}V}{\mathrm{d}\tau}=\omega\,U+K_1\,E\,V}\\\\\displaystyle{E'=Y-\hat{B}\,U-\hat{C}\,V}\\\\\displaystyle{E=1-U^2-V^2}\end{cases}\tag{13}\end{eqnarray}
ただし、\(\omega\) は係数器で変化させる可変係数とし、0から1までの値を取るものとする。また、各演算変数 \(U\)、\(V\)、\(Y\)、\(\hat{B}\)、\(\hat{C}\) の電圧換算係数はすべて1.0とした。\(E\)、\(E'\) はそれぞれ原変数の \(\epsilon\)、\(\varepsilon\) に対応する演算変数であり、\(Y\) は測定対象系の出力信号の演算変数である。測定対象系へ入力する信号は \(U\) である。 XYレコーダには x軸（実軸）に \(\hat{B}\)を、y軸（虚軸）には\(\hat{C}\) を入力する。
ここで、\(\hat{C}\) と \(\hat{B}\)を生成する再急降下ブロックのゲイン \(K_0\) と、 \(U\)、\(V\) を生成する部分の振幅修正ブロックのゲイン \(K_1\) は、それぞれの吸着速度が実質的時間換算係数 \(a_\tau/K_0\)、 \(a_\tau/K_1\)、 に支配されることに着目して値を選定する必要がある。
周波数を係数器のつまみで変更しながら、XYレコーダにベクトル軌跡を描きだすことを目標とする以上、\(\hat{B}\)、\(\hat{C}\) の周波数成分は \(U\)、\(V\)、\(Y\)のそれよりも何桁か低くなければならない。\(Y\) 波形の振動に合わせてペンレコーダが動揺すれば、周波数特性を描画できなくなるからである。ゆえに、今回は \(U\)、\(V\)、\(Y\) の各変数を生成する積分器（つまり\(U\)、\(V\) 発振系の2つの積分器と、被測定対象の伝達関数を模擬する部分のすべての積分器群）を高速モード（R=1 MΩ, C=1 nF）にし、他の積分器（ \(\hat{B}\)、\(\hat{C}\) を生成する積分器）は低速モード（R=1 MΩ, C=1 µF）にすることにした。式の上でいえば、 時間換算係数を \(a_\tau=10^{-3}\) としたうえで、再急降下ゲインを \(K_0=10^3\) と設定することにほかならない。これにより \(a_\tau/K_0=1.0\) となり、 \(U\)、\(V\)、\(Y\) の系の時間換算係数 \(a_\tau=10^{-3}\) と比較すると1000倍となっている。これは、時間換算係数の定義式 \(\tau=a_\tau\,t\) から、再急降下ブロックは \(U\)、\(V\) 生成部と測定対象系より 1000倍 ゆっくりな時定数をもつことを示している。これでXYレコーダのペンは、ある角周波数値において \(\hat{B}\)-\(\hat{C}\)平面上の1点に位置するようになる。
\(U\)、\(V\) 生成部の振幅修正ブロックのゲイン \(K_1\) については、\(K_1=10\) を選んだ。 \(a_\tau/K_1=10^{-2}\) となり、\(U\)、\(V\) の周波数成分の10倍の時定数で誤差が修正されることを意味する。
以上を勘案すると、最終的な演算方程式は以下のようになる。
\begin{eqnarray}\begin{cases}\displaystyle{\frac{\mathrm{d}\hat{B}}{\mathrm{d}\tau}=\,E' \,U}\\\\\displaystyle{\frac{\mathrm{d}\hat{C}}{\mathrm{d}\tau}=\,E'\,V}\\\\\displaystyle{10^{-3}\,\frac{\mathrm{d}U}{\mathrm{d}\tau}=\omega\,V}\\\\\displaystyle{10^{-3}\,\frac{\mathrm{d}V}{\mathrm{d}\tau}=-\omega\,U+10\,E\,V}\\\\\displaystyle{E'=Y-\hat{B}\,U-\hat{C}\,V}\\\\\displaystyle{E=1-U^2-V^2}\end{cases}\tag{14}\end{eqnarray}
測定対象系を表す式は書かれていないが、そのブロック内の全積分器が高速モードになっていればよい。ブロックダイヤグラムは以下のようになる。
図1 ベクトル軌跡描画ブロックダイヤグラム
3.結果
3.1 n次遅れ系
測定対象系が
\[G(j\omega)=\frac{1}{(1+j\omega)^n}\tag{15}\]
であるとする（\(n\in{\mathbb{N}}\)）。たとえば3次遅れ系は1次遅れ要素の3縦続接続と解釈できるので、その演算回路は以下のようになる
図2 3次遅れ系のブロックダイヤグラム
\(n=2\)~\(5\) の場合のベクトル軌跡の描出結果（数値計算解とアナログ演算解）を以下に示す。
図 3 n次遅れ要素のベクトル軌跡
\(0.01<\omega<10.0\) [rad/s] 程度の範囲で描画可能で、ある程度信頼のおける解がえられたのは \(0.10<\omega<2.00\) [rad/s] 程度の範囲であった。次数が大きくなるほど、位相遅れが大きくなっていることがわかる。
3.2 不安定零点
3次遅れ系
\[G_1(j\omega)=\frac{1}{(1+j\omega)^3}\tag{16}\]
に右半平面零点を導入した、単位時定数・単位ゲインの非最小位相系
\[G_2(j\omega)=\frac{1-j\omega}{(1+j\omega)^3}\tag{17}\]
\[G_3(j\omega)=\frac{(1-j\omega)^2}{(1+j\omega)^3}\tag{18}\]
を調べる。例として \(G_3(j\omega)\) のブロックダイヤグラムを作成するため、伝達関数を展開していく。
\begin{eqnarray}G_3(j\omega)&=&\frac{(1-j\omega)^2}{(1+j\omega)^3}\\\\&=&\frac{1}{1+j\omega}\biggl(\frac{1-j\omega}{1+j\omega}\biggr)^2\\\\&=&\frac{1}{1+j\omega}\biggl(\frac{2}{1+j\omega}-1\biggr)^2\\\\&=&\frac{1}{1+j\omega}\Biggl\{\biggl(\frac{2}{1+j\omega}\biggr)^2-2\,\biggl(\frac{2}{1+j\omega}\biggr)+1\Biggr\}\tag{19}\end{eqnarray}
1次遅れ要素と1次オールパス要素×2の縦続接続とみても間違いではないが、上述展開式の結果に基づけば、演算器の数を最小に抑えられる。
ここでの注意点は振幅スケーリングである。伝達関数をみると、中段・後段の1次遅れはゲイン 2 のものを用意したほうが見通しが良いと思われるかもしれないが、ゲインが1以上の積分器は低周波側で出力飽和が生じるため、積分器のゲインは1のままにしておいて、最終段の加算器側とその直前の係数器でつじつまを合わせたほうがいい。すなわち、実装すべき式は以下の通りになる。
\begin{eqnarray}G_3(j\omega)&=&\frac{1}{1+j\omega}\Biggl\{4\,\biggl(\frac{1}{1+j\omega}\biggr)^2-4\,\biggl(\frac{1}{1+j\omega}\biggr)+1\Biggr\}\tag{20}\end{eqnarray}
被測定系のブロックダイヤグラムおよびベクトル軌跡の1例を図に示す。
図 4 被測定系 \(G_3\) のブロックダイヤグラム
図 5 \(G_1\)、\(G_2\)、\(G_3\) のベクトル軌跡
右半平面零点（不安定零点）が増えるほど、同一周波数点での位相の遅れが目立っているうえ、高周波になってもゲインは落ちない。これこそが伝達関数に付加されたオールパス要素の効果にほかならない。負の実軸との交点（虚部=0）はそれぞれ \(G_1\,:-0.125, \;G_2:\,-0.5,\;G_3\,:-0.809...\) にあり、これらの系を開ループ伝達関数とみなした場合、 \(G_1(j\omega)=1/(1+j^\omega)^3\) が最も安定余裕があることが読み取れる。このように、開ループ伝達関数の不安定零点は位相余裕を悪化させ、閉ループの安定性制御を困難にする。
つぎに、伝達関数 \(G_3\) のステップ応答を低速モードで描出したものを以下に示す。これは被測定対象系にステップ波形を入力した場合の出力時系列波形をプロットしたもので、図1のようなベクトル軌跡描出演算回路は利用していない。
図 6 \(G_3\) と2次遅れ、3次遅れのステップ応答
\(G_3\) の過渡応答波形において、初期時刻近辺において逆応答が存在する。
このようにあらためてグラフ化されたものを目にすると、\(G_3(j\omega)=\frac{(1-j\omega)^2}{(1+j\omega)^3}\) のステップ応答の解析的表現は複雑だろうと予想されるかもしれないが、ブロックダイヤグラムからわかるように、実のところは1次遅れ～3次遅れのステップ応答の重ね合わせにすぎない（線形系だから当然であるが。）
4.参考文献
[1] 藤田広一. 「アナログ計算機のプログラム」昭晃堂, 第12版, pp.202-206, 1978.
[2] R. H. Kohr, C. H. Sprague. "A Multipurpose Frequency-Response Analyzer." ASME. J. Basic Eng. vol 92,No 2, pp. 223-232, 1970. https://doi.org/10.1115/1.3424976
[3] R. H. Kohr. "An automatic frequency response analyzer." SIMULATION, vol.8, No.1, pp.13-15, 1967.
[4] C. H. Sprague. "A flexible sinusoidal response analyzer." SIMULATION, vol.11, No.2, pp.60-64, 1968.