制御工学への応用

Copyright (C) 29Feb2004 coskx
last update Jan2012

【１】制御工学への応用

この文書はScilabの制御工学への応用を紹介をしています。

【２】一次遅れ要素

一次遅れ要素，あるいは一次遅れ系は

Ｔ d/dt ｙ(t) + ｙ(t) = ｋｘ(t)
　　（　ｘ(t)：入力　ｙ(t)：出力　Ｔ，ｋ：定数　）

で表されます。ラプラス変換により，伝達関数は

Ｇ(s)　＝　ｋ　/　（ Ｔs + 1 ）

となります。ここで Ｔ＝1，ｋ＝2 のもとで過渡応答を求めてみよう。

（１）単位ステップ応答

clf(); s=poly(0,'s'); G=2/(s+1); mySys1=syslin('c',G); t=[0:0.01:10]; y=csim('step',t,mySys1); plot2d(t,y,5,rect=[0,0,10,2.1]); a=get("current_axes"); //get the current axes a.grid=[2 2]; //make x,y-grid a.labels_font_size=4; //目盛の数字 a.title.text="unit step response G(s)=2/(s+1)"; a.title.font_size=4; a.x_label.text="time[sec]"; a.x_label.font_size=4; a.y_label.text="response"; a.y_label.font_size=4;

変数 s の定義（宣言としておこう） 伝達関数 G の定義 リニアシステム mySys1 の定義　「G」をもとに連続系（'c'）として定義 時間行ベクトル t の定義　[0 0.01 0.02 0.03 .... 9.98 9.99 10] 線形システム mySys1 のシミュレーション 　　ステップ入力とし，時刻データは時間行ベクトル t を用いる 　　シミュレーション結果は行ベクトル y に格納される 赤でｔ-ｙグラフを描く。ただし値の範囲はt：0〜10，y：0〜2.1 グラフの詳細設定

単位ステップ応答のグラフ

 

（２）周波数応答

リニアシステムの定義が済めば， bode 関数を用いてボード線図を描くことができます。

clf(); s=poly(0,'s'); G=2/(s+1); mySys1=syslin('c',G); bode(mySys1); f=get("current_figure"); f.children(1).grid=[2 2]; f.children(2).grid=[2 2]; f.children(1).labels_font_size=2; //目盛の数字 f.children(2).labels_font_size=2; //目盛の数字 f.children(1).x_label.font_size=2; f.children(2).x_label.font_size=2; f.children(1).y_label.font_size=2; f.children(2).y_label.font_size=2;

変数 s の定義（宣言としておこう 伝達関数 G の定義 リニアシステム mySys1 の定義　「G」をもとに連続系（'c'）として定義 ボード線図を描く 図の詳細設定 （オブジェクト"current_figure"のchildren(1)(2)は2つの図のAxes をそれぞれ指している。）

ボード線図　（周波数の単位はHz）

リニアシステムの定義が済めば， nyquist 関数を用いてナイキスト線図を描くことができます。

clf(); s=poly(0,'s'); G=2/(s+1); mySys1=syslin('c',G); nyquist(mySys1); a=get("current_axes"); //get the current axes a.grid=[2 2]; //make x,y-grid a.labels_font_size=2; //目盛の数字 a.title.font_size=2; //a.x_location="origin"; //a.y_location="origin"; a.x_label.font_size=2; a.y_label.font_size=2;

変数 s の定義（宣言としておこう） 伝達関数 G の定義 リニアシステム mySys1 の定義　「G」をもとに連続系（'c'）として定義 ナイキスト線図を描く 図の詳細設定

ナイキスト線図　（周波数の単位はHz）

 

（３）状態空間のアプローチ

先に与えた一次遅れ要素の微分方程式の記述は

d/dt ｙ(t) = -1/Ｔ ｙ(t) + ｋ/Ｔ ｘ(t)　

で表すこともできます。

状態空間表現ではＡ，ｂ，ｃはそれぞれ1行1列の特別な場合の行列となります。

Ａ＝-1/Ｔ，　ｂ＝ｋ/Ｔ，　ｃ＝１

次の記述でもリニアシステムの記述が可能です。

-->A=[-1] -->b=[2]; -->c=[1]; -->mySys1=syslin('c',A,b,c);

補足　一般的な状態空間表現の扱いとsyslin関数

sl=syslin(dom,A,B,C [,D [,x0] ])
   
represents the system : 
      s x  = A*x + B*u
        y  = C*x + D*u
      x(0) = x0

 

【３】二次遅れ要素

二次遅れ要素，あるいは二次遅れ系は

d²/dt² ｙ(t)  ＋  2 ζ ω_n d/dt ｙ(t)  ＋  ω_n²ｙ(t)  ＝  ｋｘ(t)
　　（　ｘ(t)：入力　ｙ(t)：出力　ζ，ω_n，ｋ：定数　）

で表されます。ラプラス変換により，伝達関数は

Ｇ(s)　＝ （　ｋ　ω_n² ）／（　s ²　+ 2 ζ ω_n s + ω_n²　）

となります。ここで ζ＝0.2, ω_n＝10，ｋ＝2 のもとで過渡応答を求めてみよう。

（１）単位ステップ応答

clf(); s=poly(0,'s'); G=200/(s^2+2*0.2*10*s+100); mySys2=syslin('c',G); t=[0:0.002:5]; y=csim('step',t,mySys2); plot2d(t,y,5,rect=[0,0,5,3.5]); a=get("current_axes"); //get the current axes a.grid=[2 2]; //make x,y-grid a.labels_font_size=4; //目盛の数字 a.title.text="unit step response G(s)=200/(s^2+2*0.2*10*s+100)"; a.title.font_size=4; a.x_label.text="time[sec]"; a.x_label.font_size=4; a.y_label.text="response"; a.y_label.font_size=4;

変数 s の定義（宣言としておこう） 伝達関数 G の定義 リニアシステム mySys2 の定義　「G」をもとに連続系（'c'）として定義 時間行ベクトル t の定義　[0 0.002 0.004 0.006 .... 4.996 4.998 5] 線形システム mySys2 のシミュレーション ステップ入力とし，時刻データは時間行ベクトル t を用いる シミュレーション結果は行ベクトル y に格納される 赤でｔ-ｙグラフを描く。ただし値の範囲はt：0〜5，y：0〜3.5 グラフの詳細設定

単位ステップ応答のグラフ

 

（２）周波数応答

リニアシステムの定義が済めば， bode 関数を用いてボード線図を描くことができます。

clf(); s=poly(0,'s'); G=200/(s^2+2*0.2*10*s+100); mySys2=syslin('c',G); bode(mySys2); f=get("current_figure"); f.children(1).grid=[2 2]; f.children(2).grid=[2 2]; f.children(1).labels_font_size=2; //目盛の数字 f.children(2).labels_font_size=2; //目盛の数字 f.children(1).x_label.font_size=2; f.children(2).x_label.font_size=2; f.children(1).y_label.font_size=2; f.children(2).y_label.font_size=2;

変数 s の定義（宣言としておこう） 伝達関数 G の定義 リニアシステム mySys2 の定義　「G」をもとに連続系（'c'）として定義 ボード線図を描く 図の詳細設定 （オブジェクト"current_figure"のchildren(1)(2)は2つの図のAxes をそれぞれ指している。）

ボード線図　（周波数の単位はHz）

周波数応答における，ピーク周波数とピークを求めます。

-->s=poly(0,'s'); -->G=200/(s^2+2*0.2*10*s+100); -->mySys2=syslin('c',G); -->peakfreq=freson(mySys2)  peakfreq  =     1.5265606 -->peak=20*log(abs(repfreq(mySys2,peakfreq)))/log(10)  peak  =     14.156688

変数 s の定義（宣言としておこう） 伝達関数 G の定義 リニアシステム mySys2 の定義　「G」をもとに連続系（'c'）として定義 ピークゲインを与える周波数を求めます。 1.52Hzでした。 ピークゲインを求めます。 14dBでした。

リニアシステムの定義が済めば， nyquist 関数を用いてナイキスト線図を描くことができます。

clf(); s=poly(0,'s'); G=200/(s^2+2*0.2*10*s+100); mySys2=syslin('c',G); nyquist(mySys2); a=get("current_axes"); //get the current axes a.grid=[2 2]; //make x,y-grid a.labels_font_size=2; //目盛の数字 a.title.font_size=2; //a.x_location="origin"; //a.y_location="origin"; a.x_label.font_size=2; a.y_label.font_size=2;

変数 s の定義（宣言としておこう） 伝達関数 G の定義 リニアシステム mySys2 の定義　「G」をもとに連続系（'c'）として定義 ナイキスト線図を描く 図の詳細設定

ナイキスト線図　（周波数の単位はHz）

 

（３）状態空間のアプローチ

先に与えた二次遅れ要素の微分方程式の記述は

ｙ(t)　＝　ｙ₁，　d/dt ｙ(t)　＝　ｙ₂　とおくことにより

| d/dt ｙ₁(t) |   |   0           1  | | ｙ₁ |   |  0 |
|             | = |                  | |     | + |    | ｘ(t)
| d/dt ｙ₂(t) |   | −ω_n²   −2 ζ ω_n | | ｙ₂ |   | ｋ |

                | ｙ₁ |
ｙ＝ [  1  0  ] |     |
                | ｙ₂ |

状態空間表現ではＡ，ｂ，ｃはそれぞれ2行2列，2行1列，1行2列の行列となります。

     |   0        1     |         |  0 |
Ａ＝ |                  |，　ｂ＝ |    |，　ｃ＝ [  1  0  ] 
     | −ω_n²   −2 ζ ω_n |         | ｋ |

次の記述でもリニアシステムの記述が可能です。
ここで ζ＝0.2, ω_n＝10，ｋ＝2 のもとで記述すると次のようになります。

-->A=[0 1; -100 -2*0.2*10]; -->b=[0 2]'; -->c=[1 0]; -->mySys2=syslin('c',A,b,c);

 

（４）ζを変化させた時の単位ステップ応答を同時に描画

先に与えた二次遅れ要素の単位ステップ応答で，複数のζについてグラフを描いてみます。

clf(); s=poly(0,'s'); G1=200/(s^2+2*0.1*10*s+100); G2=200/(s^2+2*0.3*10*s+100); G3=200/(s^2+2*0.5*10*s+100); G4=200/(s^2+2*0.7*10*s+100); G5=200/(s^2+2*0.9*10*s+100); mySys21=syslin('c',G1); mySys22=syslin('c',G2); mySys23=syslin('c',G3); mySys24=syslin('c',G4); mySys25=syslin('c',G5); t=[0:0.01:5]'; //転置して列ベクトルにする y1=csim('step',t,mySys21)'; //転置して列ベクトルにする y2=csim('step',t,mySys22)'; //転置して列ベクトルにする y3=csim('step',t,mySys23)'; //転置して列ベクトルにする y4=csim('step',t,mySys24)'; //転置して列ベクトルにする y5=csim('step',t,mySys25)'; //転置して列ベクトルにする plot2d(t,[y1 y2 y3 y4 y5],[1 3 4 5 6],rect=[0,0,5,3.5]); legend(["zeta=0.1","zeta=0.3","zeta=0.5","zeta=0.7","zeta=0.9"]); a=get("current_axes"); //get the current axes a.grid=[2 2]; //make x,y-grid a.labels_font_size=3; //目盛の数字 a.title.text="unit step response"; a.title.font_size=3; a.x_label.text="time[sec]"; a.x_label.font_size=3; a.y_label.text="response"; a.y_label.font_size=3;

ζによる単位ステップ応答の変化