自動制御系の計算では線形システムを扱うので、微分方程式をラプラス変換して得られる伝達関数の振舞を調べることが多いようです。x(t) という入力を G というシステムに入力して y(t) という出力が得られたとします。そのとき x(t) と y(t) の間の関係を記述した微分方程式をラプラス変換すると Y(s) = G(s) * X(s) という形になり、入力 X(s) と伝達関数 G(s) の積が出力 Y(s) になるという単純な関係になります。この G(s) = Y(s) / X(s) で得られる伝達関数の性質を調べればシステムの制御のためのいろいろな性質を知ることができます。

ボード線図

とくに入力が調和的な周期関数であるときは、出力も調和的な周期関数になるので、伝達関数 G(s) の s の変わりに j*w (j は虚数単位、w は角周波数ωのつもり) を代入した G(j*w) は周波数伝達関数になります。G(j*w) の振舞を半対数グラフにプロットしたものがボード線図です。Octave を使って伝達関数 G(s) のボード線図を作ってみましょう。G(s) = 1 / (1 + s) の場合の位相とゲインのボード線図を作成します。

octave:> function g = G (s)
> g = 1 ./ ( 1 .+ s );
> endfunction
octave:> w = logspace (-2, 2, 100);
octave:> s = G (j*w);
octave:> subplot (2, 1, 1);
octave:> semilogx (w, 180*angle(s)/pi)
octave:> subplot (2, 1, 2);
octave:> semilogx (w, 20*log10(abs(s)))

G(s) を定義し、G(j*w) の値を求め(jは虚数単位を表す定数)、半対数グラフに G(j*w) の位相とゲインをプロットしただけです。ゲインをデシベルで表すため絶対値の 20*log10() を取っています。subplot (2, 1, 1) は見馴れないコマンドですが、gnuplot の画面を２行１列に分割し上の方の画面を指定しています。指定された分割画面に semilogx のグラフが表示されます。semilogx は x 軸が対数の半対数グラフを描きます。

ベクトル軌跡

G(j*w) のベクトル軌跡は G(j*w) の実部と虚部をプロットするだけなので次のように簡単に描くことができます。G(s) = 1 /(1 + s)(1 + 2s) のベクトル軌跡は次のようになります。

octave:> function g = G(s)
> g = 1 ./ ((1 .+ s) .* (1 .+ 2 .* s));
> endfunction
octave:> w = logspace(-2, 4, 100);
octave:> s = G (j*w);
octave:> oneplot
octave:> grid "on"
octave:> plot(real(s), imag(s))

上のコマンドの oneplot は上の操作で一つのウィンドウに複数のグラフをプロットしていたためそれをリセットするためのコマンドです。

逆ラプラス変換

G(s) の逆ラプラス変換は G(s) = Y(s)/X(s) を微分方程式に復原することによって行います。例えば G(s) = 1 / (1 + s) の場合は次のようになります。

1 / (1 + s) = Y(s) / X(s)
(1 + s) * Y(s) = X(s)

y'(t) + y(t) = x(t)

従って x(t) にデルタ関数を代入して上の微分方程式を解いたときの y(t) が伝達関数 G(s) の時間領域でのインパルス応答関数 g(t) になります。x(t) はデルタ関数なので離散化すると x(0) = 1 で他の t では x(t) = 0 になります。従って上の微分方程式は y'(t) + y(t) = 0 を y0 = 1 (x0 = 1 より)の初期条件で解くことになります。それでは Octave を使って計算してみましょう。

octave:> function ydot = g (y, t) 
> ydot = -y;
> endfunction
octave:> y0 = 1;

y0 = 1 としたのは 入力の x(t) がデルタ関数なので、y(t) を離散化したとき y'0 + y0 = x0 = 1 となるからです( y'0 = 0 と考える )。

octave:> t = linspace (0, 10, 100);
octave:> y = lsode ("g", y0, t);
octave:> plot (y)

従って G(s) = 1 / (1 + s) の逆フーリエ変換は g(t) = exp( -t )となります。