Octave で見る古典制御

□ ２次系の過渡応答

   ２次系

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

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

   y1' = y2
   y2' = -2*0.2*y2 - y1 + 1

    octave: > function ydot = f(y, t)
              ydot(1) = y(2);
              ydot(2) = -2*0.2*y(2) - y(1) + 1;
              endfunction
    octave: > y0 = [0; 0]
    octave: > t = linspace(0, 14, 1000);
    octave: > y = lsode("f", y0, t);
    octave: > grid "on"
    octave: > plot(t, y(:, 1))
    octave: > closeplot

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

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

    octave: > function g = G(s)
              g = 1 ./ ((1 .+ s).*(1 .+ 2 .* s).*(1 .+ 3 .* s));
              endfunction
    octave: > w = logspace(-2, 4, 300);
    octave: > s = G(j.*w);
    octave: > grid "on"
    octave: > plot(real(s), imag(s))
    octave: > closeplot

□ ボード線図

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

    octave: > function g = G(s)
              g = (1 .+ 3.* s) ./ ((1 .+ s).*(1 .+ 2 .* s));
              endfunction
    octave: > w = logspace(-2, 1, 100);
    octave: > s = G(j.*w);
    octave: > grid "on"

    ゲイン線図
    octave: > semilogx(w, 20.*log10(abs(s)))

    位相線図
    octave: > semilogx(w, 180.*angle(s)/pi)

     位相が 180(-180)度を越えるとグラフが不連続になることがあるので注意が
     必要である．

    octave: > closeplot

□ 安定性の判別

   特性方程式が
     4s^3 + 6s^2 + 6s + 2 = 0
   で与えられる制御系の安定判別を行う．

    octave: > p = [4 6 6 2]
    octave: > r = roots(p)

    r の実部が全て負なので，制御系は安定である．

   特性方程式が
     s^4 + s^3 + 6s^2 + 6s + 2 = 0
   で与えられる制御系の安定判別を行う．

    octave: > p = [1 1 6 6 2]
    octave: > r = roots(p)

    r の実部に正の値を取るものがあるので，制御系は不安定である．

□ 根軌跡

   一巡伝達関数が
                     K
     G(s)H(s) = ------------
                 s(s+1)(s+2)
   で与えられる制御系の根軌跡を描く．

   特性方程式は
     s^3 + 3s^2 + 2s + K = 0
   で与えられる．

    octave: > function r = f(k)
              p = [1 3 2 k];
              r = roots(p);
              endfunction
    octave: > k = linspace(0, 30, 60);
    octave: > [nr, nc] = size(k);
    octave: > r = f(k(1));
    octave: > for i = 2:nc
                   r = [r;f(k(i))];
                 endfor
    octave: > plot(real(r), imag(r), 'x')
    octave: > closeplot