fc2ブログ

制御工学その23

6.4 MATLABを利用した演習



6.4.1 ボード線図

MATLAB使えばボード線図もかける。例えば

P(s) = 1/(s^2+s+1)

のボード線図は


>> sysP = tf([1],[1 1 1]);
>> bode(sysP);
>>


んで結果が



また、周波数領域を指定して描図することもできる。

>> bode(sysP,{10^(-1),10^1});

結果が




defaultでは対数スケールで表示される。
この辺dBと親和してるからか。

ゲイン特性と位相特性を別々にも描ける。





問題6.9

それぞれの伝達関数に関して、ボード線図を描く。

(1)

P(s) = 1/(s+1)(10s+1)

コマンド

>> sysP = 1/((s+1)*(10*s+1))

伝達関数:
1
-----------------
10 s^2 + 11 s + 1

>> bode(sysP);






(2)

P(s) = s(s+1)/10(10s+1)

コマンド略




(3)

P(s) = 1/(s+1)^4

コマンド略





6.4.2 周波数応答

y(s) = P(s)U(s)において、

P(s) = 1/(s^2+s+1)

u(t) = sin2t

とした時のy(t)の定常状態

y_app(t) = |P(2i)|*sin(2t+∠P(2i))

のグラフを描くことを考える。

コマンドは以下


>> sysP = tf([1],[1 1 1]);
>> w = 2; t = 0:0.01:20;
>> [Gg,Gp] = bode(sysP,w);
>> Gp_rad = Gp*pi/180;
>> yapp = Gg*sin(w*t + Gp_rad);
>> plot(t,yapp);hold on
>> u = sin(w*t);
>> lsim(sysP,u,t);hold off

ただし、このままだと余分な入力が描画されているのでちょちょっといじって削除。
ちょっとインチキくさいけど、グラフ描く目的は達成されてるのでいいか。




これを見る限り、5秒を過ぎたあたりからy(t)とy_app(t)が一致していることがわかる。

6.4.3 共振ピーク、共振周波数

共振ピークと共振周波数はmax関数を利用すれば求まる。



>> lsim(sysP,u,t);
>> w = logspace(-1,1,10000);
>> [Gg,Gp] = bode(sysP,w);
>> [Mp,i] = max(Gg)

Mp =

1.1547


i =

4248

>> wp = w(i)

wp =

0.7071

>>


つまり、共振ピークの式

Mp = max(|P(iω)|) = max(G_g)

をそのままmax関数で求め、さらにその時のiの値=周波数ωの値(離散値)でもって
ピーク周波数を求めたことになる。


6.4.4 ベクトル軌跡とナイキスト軌跡

ベクトル軌跡、ナイキスト軌跡も関数一発で出せる。すごい。




問題6.10

それぞれのベクトル軌跡を描く。

それぞれの伝達関数に関して、ボード線図を描く。

(1)

P(s) = 1/(s+1)(10s+1)

コマンド

>> sysP = 1/((s+1)*(10*s+1))

伝達関数:
1
-----------------
10 s^2 + 11 s + 1

>> nyquist(sysP);





(2)

P(s) = s(s+1)/10(10s+1)

コマンド略




本当に合ってるのだろうか。ちょっと調べてみる。

軌跡の最初のほうを拡大すると、(ここからscilabにバトンタッチ)




実軸付近の細かい部分を更に拡大してみる。





ω=0.3あたりで急激に曲がっているのがわかる。

P(iω)に対して、

P(iω) = (-ω^2+iω)/(100ωi+10) = (-ω^2+iω)(10-100ωi)/(100 + 10000ω^2)
= -10ω^2 + 10ωi + 100ω^3i + 100ω^2 /10000ω^2+100 = 90ω^2 + (100ω^3+10ω)i /(10000ω^2+100)
= 9ω^2 /1000ω^2+10 + (10ω^3+ω)i/1000ω^2+10

実軸部分の式

Re[P(iω)] = 9ω^2/(1000ω^2+10)

これをグラフ化すると、




大体0~8.9の範囲となっている。
さらに角度に関して、

∠P(iω) = tan-1(10ω^2+1/9ω)

調べてみようとしたけどうまくいかない。
仕方が無いので別のグラフ描画ツールを使ってみる。

y = arctan(10ω/9 + 1/9ω)*180/πのグラフ




ω=0.3を境に角度の減少がとまり、増加し始めている。
このへんから、そんなに間違ってないのではないだろうか。



(3)

P(s) = 1/(s+1)^4

コマンド略





まとめ
MATLABだんだんわかってきたような

わからないところ
配列演算がうまくいってないのか角度のグラフがうまく描けない。要検討。
スポンサーサイト



コメントの投稿

非公開コメント

検索フォーム
RSSリンクの表示
リンク
ブロとも申請フォーム

この人とブロともになる

QRコード
QR