【第12回】

<内容>

1.  周波数応答

 入力として振幅1で角周波数w の正弦波u(t) = sinwt を加えた時の出力の定常応答は

  図面1.jpg

と表される.ここで,M(w) をゲイン,f(w) を位相と呼び,これらの値により周波数応答を表現する.すなわち,ゲインM(w) = 1,位相f(w) = 0 は入力波形がそのまま出力に現れることを意味する.

 

周波数伝達関数G(jw)は一般に以下のような複素数となり,

  図面1.jpg

ゲインM(w) と位相f(w) は次式により求められる.

  図面1.jpg

 周波数応答は線図によって分かりやすく表現される.ここでは,よく利用されるボード線図Bode diagram)の描画方法について学ぶ.

 

 

ボード線図】(教科書p.137の図3.63参照)

 横軸に角周波数w をとり,ゲインM(w) をプロットしたゲイン線図と位相f(w) をプロットした位相線図の両者から成る.通常,ゲイン線図はdB(デシベル)単位で,位相線図はdeg(ディグリー:度)単位で描かれる.

 

M(w) [dB]は以下のように計算される.

  M(w) [dB] = 20 log10 M(w)

すなわち,M(w) = 1の場合は,0[dB]となる.

 

 

2.  ラウス・フルビッツの安定判別法

【安定性の定義】

  eq1

 系の伝達関数G(s) が上式のように与えられるとき,「分母多項式=0」とおいた特性多項式の根特性根),すなわち,G(s) pi, i=1,2,,nの実部が全て負であることが,系が安定であるための必要十分条件となる

  図面1.jpg

 

 

2.1  ラウスの安定判別法

 特性多項式の係数,すなわち,伝達関数G(s) の分母多項式の係数ai, i=0,1,,n-1 を用い,以下のラウス表を作る.

  eq2

ここで,

  eq3

である.ただし,存在しない係数は0とおく.

 

 系が安定であるための必要十分条件は,特性多項式の係数が全て正で,かつ,ラウス表の第1列の係数が全て正,すなわち,

  A1 > 0, B1 > 0, …, Z1 > 0

である.

 

 

<演習>

1周波数応答

11.ボード線図bode関数」

 

1次遅れ要素

図面1.jpg

 

 

練習12-1

 例題3.39(教科書p.136)を行い,結果を確認せよ.

¨         ゲイン特性は,w = 1/T まで0dBで,その点から-20dB/decの傾きで下がる折れ線に近似できる.w = 1/T 折れ点周波数と呼ぶ.

 

 

現在のシステムでは,bode関数は非プロパーな伝達関数も扱うことができるが,Mファイルの作成方法を学習するため,bode1.mを作成する.

 

 

練習12-2

 ゲイン線図を描く関数gbodeMファイル「gbode.m」,位相線図を描く関数pbodeMファイル「pbode.m」,さらにこれらを用いてボード線図を描く関数bode1Mファイル「bode1.m」を作成した後,以下の伝達関数のボード線図をbodebode1の両方で描き,同じ結果になっていることを確認せよ.

  図面1.jpg

  bode関数で描画する伝達関数はtf関数で,bode1関数で描画する伝達関数はtf1関数で入力すること.

  現在のシステムでは,bodeで描いても,位相特性が不連続になることはない(教科書p.136下部参照).

 

(注)gbode.mpbode.mbode1.mは,教科書ではなく,下記の通りに入力せよ.

 

gbode.m

function [gain,w]=gbode(sys,w,c)

% [gain,w]=gbode(sys,w,c)

if nargin<3

   c=0;

endif

x=G(sys,j*w);

gain=20*log10(abs(x));

if nargout==0

   semilogx(w,gain,c)

   xlabel('Frequency [rad/s]')

   ylabel('Gain [dB]')

   title('Gain characteristic')

   grid('on')

   clear gain w

endif

endfunction

 

gbode.mの解説>

¨         absは絶対値(ゲイン)を求める関数である.

¨         semilogxx軸を対数とする片対数グラフを描画する関数である.

¨         gbode3番目の引数cは,グラフの線の色を指定するものである.3番目の引数を指定しない場合は,デフォルトの色で描画される.

 

pbode.m

function [phase,w]=pbode(sys,w,c)

% [phase,w]=pbode(sys,w,c)

if nargin<3

   c=0;

endif

x=G(sys,j*w);

ph=angle(x)*180/pi;

n=length(ph);

p1=ph(1);

base=0;

for k=2:n

   p2=p1;

   p1=ph(k);

   if abs(p1-p2)>180

      base=base-sign(p1-p2)*360;

   endif

   ph1(k)=ph(k)+base;

endfor

phase=ph1;

if nargout==0

   semilogx(w,phase,c)

   xlabel('Frequency [rad/s]')

   ylabel('Phase [deg]')

   title('Phase characteristic')

   grid('on')

   clear phase w

endif

endfunction

 

pbode.mの解説>

¨         angleは位相[rad]を求める関数である.phの単位は[deg]である.pip である.

¨         abs(p1-p2)>180の条件で不連続点を検出している.不連続点では,隣接するデータの位相が突然180[deg]以上飛ぶことになる.

 

bode1.m

function bode1(sys,w)

% bode1(sys,w)

if (nargin<2)

   w=logspace(-2,2,200);

endif

subplot (2, 1, 1);

gbode(sys,w)

subplot (2, 1, 2);

pbode(sys,w)

endfunction

 

bode1.mの解説>

¨         nargin < 2の場合,すなわち,wを指定していない場合は,デフォルト設定でwを与えている.

¨         subplotは複数のグラフを描画する際に使用する.上の例では,キャンバスを上下に分割し,上にゲイン線図を,下に位相線図を描くように指定している.詳細はヘルプを参照のこと.

 

 

―――――――――――――――――――――――――――――――――――――――――――――――――――――――

 

 

2次遅れ要素

図面1.jpg

 

 

練習12-3

 例題3.40(教科書p.140)を行い,結果を確認せよ.

(注)実行は,教科書の通りではなく,下記の通りにせよ.

 

【図3.65の描画方法】

>> close

>> g1=tf1(1,[1 0.6 1]);

>> g2=tf1(1,[1 1 1]);

>> g3=tf1(1,[1 2 1]);

>> g4=tf1(1,[1 4 1]);

>> w=logspace(-2,2,100);

>> gbode(g1,w,'1');

>> hold on

>> gbode(g2,w,'2');

>> gbode(g3,w,'3');

>> gbode(g4,w,'4');

>> legend('line 1','line 2','line 3','line 4')

 

【図3.66の描画方法】

>> close

>> pbode(g1,w,'1');

>> hold on

>> pbode(g2,w,'2');

>> pbode(g3,w,'3');

>> pbode(g4,w,'4');

>> legend('line 1','line 2','line 3','line 4')

 

注)

¨         hold onは重ね描きのための命令で,描画した図を消さない設定にする.この設定を取り消すには,hold off命令を入力するか,close命令を入力すればよい.

¨         グラフの色指定は,以下の通りである.

'0''k'black

'1''r'red

'2''g'green

'3''b'blue

'4''m'magenta

'5''c'cyan

'6''w'white

¨         システム既存のbodeコマンドで実行するときは下記のように入力する.

>> close

>> g1 = tf(1,[1 0.6 1]);

>> g2 = tf(1,[1 1 1]);

>> g3 = tf(1,[1 2 1]);

>> g4 = tf(1,[1 4 1]);

>> w = logspace(-2,2,100);

>> bode(g1,g2,g3,g4,w)

 

 

2.ラウス・フルビッツの安定判別法

21.ラウスの安定判別法

ラウス表(教科書p.144の表3.5)に基づいて,安定判別を行う.

 

 

練習12-4

 ラウスの安定判別を行う関数routhを記述したMファイル「routh.m」を作成した後,例題3.42(教科書p.147)を行い,結果を確認せよ.

(3.55)が間違っている.ただしくは,D(s)=s4+5s3+2s2+2s+1 である.

 

routh.mの解説>

¨         pは,特性多項式の係数を持つベクトルである.

¨         ベクトルpの要素数はn+1個である.ここで,nは特性多項式におけるsの最高次数を意味する.

¨         p1 = [ p(1)/p(1)  p(2)/p(1)  … p(n1)/p(1)  0]である.これは,特性多項式の係数をsの最高次数の係数p(1)で割り,sの最高次数の係数を1としている.さらに,ラウス表を作る際に,n1が奇数の場合は最終列の第2行目の係数が存在しなくなる.係数が存在しない所には0を入れるので,そのためにp1の最後に0を加えている.

¨         while文を抜けたときには,jには「ラウス表の列数」+1の値が入っている.したがって,ラウス表の列数n2j-1としている.

¨         n1行×n2列のラウス表rを作成する.先に求めた2行分の配列r2行×n2列)の下の全要素(n1-2行×n2列)に0を入れておく.

¨         rem(x,y)x/yの剰余を計算する関数である.

¨         2行目以降の第1列要素が10-6 より小さい場合は,breakによりfor文の処理を抜ける.なぜなら,次の行の係数を計算する際,その要素が分母にあるので,値が十分に小さいと計算値が無限大になり意味を成さないからである.

¨         fix(x)xの小数点以下切り捨てによる整数化を行う関数である.

¨         ラウス表の列数は2行ごとに1列だけ減少していくので,関数fix(x)を利用して以下のように処理している.

for j = 1 : n2 - fix(k/2)

  r(i, j) = (r(i-1, 1) * r(i-2, j+1) - r(i-2, 1) * r(i-1, j+1)) / r(i-1,1);

endfor

k = k + 1;

注)プログラム行が複数行にわたる場合は,末尾に「...」を入力すればよい.教科書では,上記のプログラムの r(i, j) の部分を2行で入力している.

¨         ラウス表rの第1列要素に0(ここでは,10-6 より小さい値)があれば,s = 0としている.

 

注)関数routhを使う前に,特性多項式の全係数が全て正という条件が満たされているかを確認した後に,この関数を使用する必要がある.

 

 

【レポート課題(11)(授業終了後にアップロード)

  RePORT