周波数領域のOctave

Octaveを用いた周波数領域の制御工学

主な命令は以下のとおりです。

 1s=tf('s');
 2G=1/(s^2+s+1);%伝達関数の設定その1
 3G=tf([分子の係数],[分母の係数]);%伝達関数の設定2
 4G=zpk([零点の値],[極の値],ゲイン);伝達関数の設定3
 5[Mag Phi]=bode(G);戻り値を取得すると描画しません
 6pzmap(G);%極と零点を描画
 7rlocus(G);%根軌跡
 8nyquist(G);%ナイキスト線図
 9impulse(G);%インパルス応答
10step(G);%ステップ応答
11lsim(G,u,t);%任意の入力uに対する応答

細かい使い方の例を以下に示します。

サンプルソースと結果を示します。

ボード線図
ボード線図
極位置
極位置(minrealで極零整理)
根軌跡
ナイキスト線図
ランプ応答
ステップ応答
ステップ応答(重ね書き)
ステップ応答(window分割)
ステップ応答(文字挿入)
ステップ応答(重ね書き,拡大図)

github/b2_control_step1.m

  1clear all%変数クリア
  2close all%図のwindowをすべて閉じる
  3clc%コマンドウィンドゥをクリア
  4
  5pkg load control%tfなど制御関係の命令を使うためにcontrol packageを読み込む。
  6
  7%伝達関数の設定その1
  8s=tf('s');
  9G=1/(s^2+s+1);
 10%伝達関数の設定その2
 11G1=tf(1,[1 2 1]);%係数だけ設定1/s^2+2s+1
 12%伝達関数の設定その3
 13%零点、極、システムゲインを設定。符号に注意
 14G2=zpk([-3],[-1 -2],4);%(s+3)/[(s+1)(s+2)] * 4
 15G3=zpk([],[-1 -2],4);%零点がないとき:1/(s+1)(s+2) * 4
 16
 17%分母分子の係数取り出し
 18[num den]=tfdata(G1);
 19[z p k]=zpkdata(G2);
 20
 21n=den{1,1};%極を手動で計算する時,sの高次式の係数を取りだす
 22roots(n);%高次式の解を求める
 23
 24
 25%%以下の描画において、関数の結果を引数で戻した場合、描画されないので
 26%%目的用途に応じて使い分けすること。
 27
 28%bode線図
 29[Mag Phi]=bode(G1);%描画または数値を取得(自分で描きたいとき)
 30w=0.6;%指定周波数(rad/s)
 31[Ma Ph]=bode(G1,w);%指定した周波数wの周波数伝達関数の大きさと位相角を取得
 32%いずれの場合もMagは大きさなので、dBにするときは20log10(Mag)とすること
 33bode(G1)%bode線図自動描画
 34
 35%自分で周波数範囲を指定するとき
 36w2=logspace(-2,3,100);%(start,end,個数),(10^start)から(10^end)まで 個数分のw
 37[Mag2 phi2]=bode(G2,w2);
 38
 39[gm pm wg wp]=margin(G1);%ゲイン交差周波数と位相余有、位相交差周波数とゲイン余有
 40
 41figure
 42%自分でボード線図を描くとき(片対数グラフの書き方)
 43semilogx(w2,20*log10(Mag2))
 44
 45
 46
 47%%極と零点
 48%分母分子を整理(約分):高次-->低次にする。
 49G4=G1*G2/(1+G1*G2);
 50G4m=minreal(G4);
 51
 52pole(G4);%極位置を取得
 53zero(G4);%零点だけ求まる。
 54[p z]=pzmap(G4);%極と零点を取得
 55figure(10)%fig window番号指定
 56pzmap(G4)%極と零点を描画
 57figure(11)
 58pzmap(G4m)%minrealの前後を比較
 59
 60%根軌跡
 61[RLDATA K]=rlocus(G);%数値だけ取得
 62%Kは実軸と交差するときの値???要確認
 63figure%fig window番号を指定しないと前の続きになる
 64rlocus(G1);%図を描画
 65sgrid on
 66%または
 67zeta=0.5;
 68sgrid(zeta,[])
 69
 70daspect([1 1])%図の縦横比を1:1にするとき
 71
 72%ナイキスト線図
 73[Re Im W]=nyquist(G);%数値だけ取得
 74figure
 75nyquist(G);%自動描画
 76
 77
 78
 79%%時間応答
 80dt=0.01;%離散間隔(刻み時間)
 81ti=0:dt:10;%時間指定
 82
 83%impules応答
 84[y t x]=impulse(G);
 85%step応答
 86[y t x]=step(G,ti);%時間tiを指定
 87[y]=step(G);%だけでもok
 88%ramp応答
 89figure
 90ramp(G);%傾き(du/dt)1
 91
 92%任意の傾きのランプ入力を作る場合
 93du=0.1;%du/dt
 94u=zeros(1,length(ti));
 95for i=2:length(ti);
 96    u(1,i)=u(1,i-1)+du;
 97end                 %
 98
 99%任意の入力uの応答
100[y t x]=lsim(G,u,ti);%時間刻みtを指定
101
102
103%%図を描く
104figure%新しいfigure windowを開ける
105step(G1)%図は自動で描かれる
106
107%自分で描きたいときは、数値を取得する。
108[y1 t1]=step(G1,ti);%刻み時間指定:不安定になることもあるので注意
109
110figure(20)%番号を指定して、新しいfigure windowを開ける
111plot(ti,y1);
112[y2 t2]=step(G2,ti);%刻み時間指定:不安定になることもあるので注意
113hold on%図を重ね描きするおまじない
114lh=plot(ti,y2);%オプション指定時は、ハンドル(lh)を取得して設定すると複数設定時は省力でき便利。
115set(lh,'linewidth',2,'linestyle','-','color','r','linewidth',4);
116%線のスタイルの設定,plot(t,y,'r-','linewidth',4)と同じ,入力がラクな方を選んでください
117set(gca,'xlim',[0 5],'ylim',[0 1]);%x,y軸の最小最大値の設定
118set(gca,'xtick',0:0.5:10,'ytick',0:0.1:1);%x,y軸の目盛り間隔を指定するとき
119%set(gca,'xticklabel',0:1:10,'yticklabel',0:0.2:1);%x,y軸の目盛り数値表示を指定するとき
120set(gca,'xlabel','time','ylabel','y out');%x,y軸の名前の設定, xlabel('time')でもok
121
122%set(gca,...)は一行にまとめても可
123legend('G1','G2','location','southeast');%注釈を図中に記すとき。詳細はhelp legend参照
124title('time step')%タイトルをつけたいとき:レポートでは図の下に個別に明記することが望ましいので非オススメ
125
126%%1つのwindowに複数の図を描く
127%重ね描きするときはhold onが必要
128figure%番号を指定しないと、以前のfig windowの続き番号になる
129plot(t1,y1);
130hold on%図を重ね書きするときのおまじない
131plot(t2,y2);
132%任意の2点(x0,y0)(x1,y1)を結びたいときは
133%plot([x0 x1],[y0 y1]);
134
135%画面を分割するときは、subplotで分ける
136figure
137subplot(2,1,1)%2行1列に分割し、1番目の場所に描く
138plot(t1,y1);
139subplot(2,1,2)%2行1列に分割し、2番目の場所に描く
140plot(t2,y2);
141
142figure
143step(G1,G2);%自動で重ね描きするとき。線オプションの手動加工はできない。
144hold on
145plot([0 1.2],[3.79 3.79])
146text(0.8,3.9,'63.2%')%図中に文字を置く,text(x座標,y座標,'文字')
147
148%%図をファイルに保存する
149fname='fig_';%図のファイル名の共通部.
150ext1='.svg';%scalable vector graph:拡大/縮小しても図はきれいで品質良い.wordでは張り込めないかも
151%svgはinkscapeやtext editorで編集できる,xml(HTML5)形式で書かれている
152
153ext2='.png';%portable network graphics. bitmap形式。 pbm,jpgもあり
154for i=[1:11 20]%保存する図番を指定:そのためfigure(n)として描画したほうがわかりやすいかも
155    print(i,[fname,num2str(i),ext2],'-dpng','-S640,480')
156    %print(i,[fname,num2str(i),ext1],'-dsvg','-S640,480')
157end
158
159%保存した図が小さい場合やファイルが生成できないなどの不具合がある場合は、
160%PC環境によるものなので、手動で図を保存すること。
161%拡大表示して、FILE->SAVE_AS> file type :pngを選択して名前をつけて保存

このシリーズの投稿