scilabで制御系設計

2013.4.23:周波数応答追記
2012.5.14:初版,6.8:r1
scilabを使って制御系を設計しましょう。

1.installと実行

scilabはじめましょうを参考に

helpを使おう!!

Scilabはhelpがとても充実しています。
Scilabコンソール上で
help 知りたい命令
とすると、使い方やオプション設定などの解説が現れます。
あちこちの webサイトにScilabの情報もありますが、古い情報も多いので、ぜひhelpを利用してください。

2.よく使う命令

//以降の1行はコメントになります。ブロックコメントはないようなので、 SciNotesの編集機能(行選択後、編集->コメント選択)を使用してください。
また、行末にセミコロン;をつけると、 Scilabコンソールでは結果出力が抑制されますが、Scinotesに書くスクリプトファイルでは効き目がないようです。 注意してください。
  1. (伝達)関数の設定

  2. 連続時間システムの定義 上記で定義したsの関数を連続時間システムの関数として定義します。
    G=syslin('c',Gs);//Gsは上記で定義したsの関数
    
    これをしないと、ボード線図を描いたり時間応答を計算したり、他で利用できません。忘れずに。

    状態方程式表現の行列を用いる場合は、
    A=[0,1;-1;-1];
    B=[0;1];
    C=[1,0];
    D=0;
    G=syslin('c',A,B,C,D)
    
    で良いようです。A,B,C,Dは定番どおりdx=Ax+Bu,y=Cx+Duの行列です。

  3. 出力フィードバック
  4. PID制御などの出力フィードバックは、コントローラの伝達関数をC として閉ループ伝達関数Gを構築して、シミュレーションを行います。
    C=Kp(1+1/Ti/s+Td*s/(1+0.1*Td*s))
    G=P*C/(1+P*C)
    y=csim('c',t,G)
    

  5. 状態フィードバック
  6. u=-fx+rなので、閉ループの状態方程式はdx=(A-Bf)x+rとなります(r:目標値)。状態方程式を書き換えて
    Af=A-B*f;
    Bf=[1; 0];//目標値r=x1
    G=syslin('c',Af,Bf,C,D)
    
    とします。
    fはフィードバック係数です。極配置で指定するときは、
    poles=[-2+%i, -2-%i];//%iは虚数単位
    f=ppol(A,B,poles);
    
    とすれば、fが求まります。スゴク楽に設計できます!! spec(Af)で、固有値を表示してくれます。

  7. 時間応答
    [y [,x]]=csim(入力,シミュレーション時間,連続時間システムの伝達関数)
    

  8. 各種図
  9. 伝達関数の評価

    2013.4.22 New
    2つめのhornerでもよいですが、周波数応答を直接評価する方法があります。
    repfreq(G,wg/2/%pi)                    //Gは連続時間システムの伝達関数。wg rad/sでの複素数が求まります
    [db phi]=dbphi(repfreq(G,wg/2/%pi));   //wg rad/sでのゲイン(db dB)と位相(phi deg)が求まります
    [phi db]=phasemag(repfreq(G,wg/2/%pi));//dbphiとは戻り値が逆順に求まります
    
    戻り値が2つあるので注意してください。

    ある周波数のゲインや位相を取り出したいときは、hornerも使えます。
    gfs=horner(G,%i*2*%pi*f_s); //G(f_s*j2pi)の計算
    g_fs=20*log10(abs(gfs));    //ゲイン
    p_fs=phasemag(gfs,'c');     //位相
    
    お好みでどうぞ

  10. ゲイン余裕、位相余裕を取り出す。 結果を見る(数値を知りたい)ときは、ワークスペースを見るか、コンソールで変数名(enter)とするか、 disp(変数名)としましょう。
    前述のとおり周波数の単位はHzなので注意してください。

  11. ズームしたい
    zoom_rect([xmin, ymin , xmax, ymax])
    
    または
    ax=gca();
    ax.data_bounds=[xmin, ymin ; xmax, ymax];
    
  12. 数値を出力したい。
  13. sceファイルでのスクリプトでは、変数や結果を画面などに出力するためにはおまじないが必要です。
    1. 画面に出力する
      disp(変数名)
      
      Scilabの dispはヒジョーにくせがあって後入れ先出し(LIFO)です。 つまり、
      disp(ho, ge)
      
      とすると、
      ge
      ho
      
      の順に現れます。
      disp(ho,"ho", ge,"ge")
      
      などとして、何を出力しているのかわかるようにしたほうが、間違いないと思います。
    2. ファイルに出力する
      print('filename',x1,....)
      
      help printとしてみてください。
      mprintfなどなど 他にもいろいろあるようです。

3.グラフィックスの技

複数の伝達関数の応答を重ねて比較するなど、重ね描きはぜひしたいもの。 しかし、上記の命令を単純に並べると、うまく描画できない場合があります。
そんな時は、計算結果だけを取得して、自分で描きましょう。
plot2d([x座標データ],[y座標データ])が基本です。
したがって、(x1,y1)から(x2,y2)に直線を引きたいときは
plot2d([x1 x2],[y1 y2])
とします。 その他 もっといろんなことをしたり、 線の種類を変えたい時は以下のようにします。
  1. グラフィックエディタを使う。
      グラフィックウィンドウで、軸や線の種類を設定できます。
    1. グラフィックウィンドウの編集のボタン押す。
    2. エンティティ選択を開始、を選択し、 変更したい線の上で左クリックを押します。
    3. PolylineEditorが現れますので、Styleのタブの中から、線種と太さ(Line)、色(Foreground)、などを選択します

  2. エンティティを直接設定する。
    gce()で直近に書いたエンティティが取得できます。そのchildrenで各種設定を行います。 たとえば、以下の用に遊べます。
    fh=gce();//
    fhc=fh.children;
    fhc.line_mode="off";//線は書かず、markだけ
    fhc.thickness=1;
    fhc.line_style=4;
    //[0 8] の範囲の整数を指定します.
    // 0 および 1は実線, 2:破線,3:1点鎖線,4:長い1点鎖線,5:非常に大きい1点鎖線, 6:非常に大きく長い1点鎖線,点,7:2点)の選択を意味します.
    fhc.mark_style=1;//[0--14]
    fhc.mark_foreground=color("red");//markの線の色
    fhc.mark_background=color("red");//markの塗りつぶし
    

4.その他プログラミングなこと