Octave Control Systems Toolbox を使う

Octave ver.2.1 以降では Contol Systems Toolbox が強化され，制御系のシミュ
レーションや，制御系の表現同士の変換が容易になっている．ここでは Debian Etch 
で配布されている ver.2.9.9 版をベースにしているが，Debian Sarge の ver.2.1.69 
でもほぼ同様に動作する．

□ ２次系の過渡応答

   ２次系

   y''(t) + 2ζω y'(t) + ω^2 y(t) = x(t)

   で ζ = 0.2, ω = 1, x(t) = 1 (単位ステップ入力)の場合の過渡応答を求める．

   このとき，伝達関数は
                   1
     G(s) = --------------
            s^2 + 0.4s + 1
   で与えられる．

    octave: > sys = tf([1],[1, 0.4, 1], 0);
    octave: > sysout(sys)
    octave: > step(sys)

    あるいは，次のようにして応答と時間を変数に保存する事ができる．

    octave: > [y, t] = step(sys);
    octave: > grid "on"
    octave: > plot(t, y)

    計算の精度をあげるには終了時間と分割数を明示すると良い．

    octave: > [y, t] = step(sys, 1, 30, 300);
    octave: > plot(t, y)
    octave: > closeplot

□ ベクトル軌跡（ナイキスト線図）

   伝達関数が
                   1
     G(s) = -----------------
            (1+s)(1+2s)(1+3s)
   で与えられる制御系のベクトル軌跡を描く．

    octave: > sys = zp([],[-1, -1/2, -1/3], 1/6);
    octave: > sysout(sys)
    octave: > nyquist(sys)
    octave: > closeplot

□ ボード線図

   伝達関数が
               (1+3s)
     G(s) = ------------
            (1+s)(1+2s)
   で与えられる制御系のボード線図を描く．

    octave: > sys = zp([-1/3], [-1, -1/2], 3/2);
    octave: > sysout(sys)
    octave: > bode(sys)
    octave: > closeplot

□ 根軌跡

   前向き伝達関数が
                 K
     G(s) = -------------
             s(s+1)(s+2)
   で与えられるユニティフィードバック制御系の根軌跡を描く．

    octave: > sys = zp([], [0, -1, -2], 1);
    octave: > sysout(sys)
    octave: > rlocus(sys, 1/60, 0, 30)
    octave: > closeplot

    ver.2.1.69 では途中でエラーが出て計算が止まるが，ver.2.9.9 では表示される．

    尚，ver.2.1.69 では OCST のデモ(rldemo)でもエラーが発生している． 

□ 伝達関数で与えられるシステムの状態変数表示

   伝達関数が
                 s+2
     G(s) = ------------
            s^2 + 5 s + 2
   で与えられる制御系の状態変数表示を求める．

    octave: > sys = tf([1, 2], [1, 5, 2]);
    octave: > sysout(sys)
    octave: > [A, b, c, d] = sys2ss(sys)

□ 状態変数表示されたシステムの取り扱い．

   システム

     x'(t) = A x(t) + b u(t)
     y(t) = c x(t)

     A = [-1 , -2; -1, -4], b = [1; 1], c = [1, 0]

    に対応する伝達関数と，ステップ応答を求める

    octave: > A = [-1 , -2; -1, -4];
    octave: > b = [1; 1];
    octave: > c = [1, 0];
    octave: > sys = ss(A, b, c);
    octave: > sysout(sys)

    伝達関数の分子多項式と分母多項式の係数を求める．
    octave: > [n, d] = sys2tf(sys)
    octave: > tfout(n, d)

    ステップ応答を求める．
    octave: > [y, t] = step(sys, 1);
    octave: > grid("on");
    octave: > plot(t, y)
    octave: > closeplot

□ ユニティフィードバック制御系の伝達関数の計算

    前向き伝達関数が
                 s+1
     G(s) = ------------
            s^2 + 2 s + 2
    で与えられるユニティフィードバック制御系の伝達関数を求める．

    octave: > G = tf([1, 1], [1, 2, 2]);
    octave: > sys = buildssic([1, -1],[],[1],[1],G);
    octave: > sysout(sys)
    octave: > [n, d] = sys2tf(sys)
    octave: > tfout(n, d)

    buildssic の代わりに，sysdup, sysscale, sysconnect, sysprune を使
    う方法もある．

    octave: > G0 = sysdup(G, [], [1]);
    octave: > G1 = sysscale(G0, [], diag([1, -1]));
    octave: > G2 = sysconnect(G1, [1], [2]);
    octave: > sys = sysprune(G2, 1, 1);
    octave: > sysout(sys)
    octave: > [n, d] = sys2tf(sys)
    octave: > tfout(n, d)

□ フィードバック制御系の伝達関数の計算

    前向き伝達関数が
                 1
     G(s) = ------------
            s^2 + 3 s + 2
    で，フィードバック要素の伝達関数が
                1
     H(s) = ---------
             1 + 2 s
    で与えられるフィードバック制御系の伝達関数を求める．

    octave: > G = tf([1], [1, 3, 2]);
    octave: > H = tf([1], [2, 1]);
    octave: > sys = buildssic([1, -2; 2, 1],[],[1],[1], G, H);
    octave: > sysout(sys)
    octave: > [n, d] = sys2tf(sys)
    octave: > tfout(n, d)

    buildssic の代わりに，sysgroup, sysdup, 他を使う方法もある．

    octave: > GH = sysgroup(G, H);
    octave: > GH0 = sysdup(GH, [], [1]);
    octave: > GH1 = sysscale(GH0,[],diag([1, 1, -1]));
    octave: > GH2 = sysconnect(GH1,[1, 2],[2, 3]);
    octave: > sys = sysprune(GH2, 1, 1);
    octave: > sysout(sys)
    octave: > [n, d] = sys2tf(sys)
    octave: > tfout(n, d)