Octave で見る現代制御

□ 可制御性・可観測性

   システム

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

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

    の可制御性と，可観測性を調べる．

    octave: > A = [-2, -1; -1, -2]
    octave: > b = [1; 1]
    octave: > c = [1, 0]
    octave: > Mc = [b, A*b]
    octave: > det(Mc)

      Mc の行列式が 0 なので可制御ではない．

    octave: > Mo = [c; c*A]
    octave: > det(Mo)

      Mo の行列式が 0 でないので可観測である．

□ 対角化

   システム

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

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

    を対角化する．

    octave: > A = [-2, -1; -1, -2]
    octave: > b = [1; 1]
    octave: > c = [1, 0]
    octave: > [T, V] = eig(A)
    octave: > Abar = inv(T)*A*T
    octave: > bbar = inv(T)*b
    octave: > cbar = c*T

    この結果，システムが可制御ではないが可観測であることがわかる．

□ 可制御正準形

   システム

     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: > p = poly(A)

      特性多項式が s^2 + 5 s + 6 となることがわかる．

    octave: > T = [b, A*b]*[p(2), 1; 1, 0]
    octave: > Abar = inv(T)*A*T
    octave: > bbar = inv(T)*b
    octave: > cbar = c*T

□ 安定性判別

   システム

     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: > p = poly(A)
    octave: > r = roots(p)

    A の特性方程式の根 r の実部が全て負なので，システムは安定である．

□ 状態フィードバック制御の極配置

   システム

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

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

    の極を -1+2j, -1-2j に設定する状態フィードバックゲイン f を求める．

    octave: > A = [-1, -2; -1, -4]
    octave: > b = [1; 1]
    octave: > c = [1, 0]
    octave: > p1 = poly(A)
    octave: > T = [b, A*b]*[p1(2), 1; 1, 0]
    octave: > p2 = poly([-1+2j, -1-2j])
    octave: > f = [p2(3) - p1(3), p2(2) - p1(2)]*inv(T)

    初期値 x(0) = [1; -1] の場合の応答を確かめる．

    octave: > global A; global b; global f;
    octave: > function xdot = F(x, t)
              global A; global b; global f;
              xdot = (A-b*f)*x;
              endfunction
    octave: > x0 = [1; -1]
    octave: > t = linspace(0, 5, 1000);
    octave: > y = lsode("F", x0, t);
    octave: > grid "on"
    octave: > plot(t, y)
    octave: > closeplot

□ 直接フィードバック制御における根軌跡

   システム

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

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

    に対して u(t) = -k y(t) という直接フィードバックを行ったときの根軌跡を
    求める．

    octave: > A = [-1, 0, 0; -1, -1, -1; -2, 1, -2]
    octave: > b = [1; 1; 1]
    octave: > c = [0, 0, 1]
    octave: > global A; global b; global c;
    octave: > function r = F(k)
              global A; global b; global c;
              p = poly(A-b*k*c);
              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

□ 同一次元オブザーバの構成

   システム

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

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

    に対して，-4, -5 に極を持つ同一次元オブザーバ

    z'(t) = A z(t) -k[c z(t) - y(t)] + b u(t) = (A - k c)z(t) + k y(t) + b u(t)
  
    のゲイン k を決定する．

    octave: > A = [-1, -2; -1, -4]
    octave: > b = [1; 1]
    octave: > c = [1, 0]
    octave: > p1 = poly(A)
    octave: > T = [c.', A.'*c.']*[p1(2), 1; 1, 0]
    octave: > p2 = poly([-4, -5])
    octave: > k = ([p2(3) - p1(3), p2(2) - p1(2)]*inv(T)).'

    オブザーバの推定の様子を確かめる．初期値 x(0) = [0; 0], z(0) = [1; 1]
    u(t) = 1 (単位ステップ入力)とする．

    octave: > AA = [A, zeros(2, 2); k*c, A-k*c]
    octave: > bb = [b; b]
    octave: > global AA; global bb;
    octave: > function xdot = F(x, t)
              global AA; global bb;
              xdot = AA*x + bb;
              endfunction
    octave: > x0 = [0; 0; 1; 1]
    octave: > t = linspace(0, 5, 1000);
    octave: > y = lsode("F", x0, t);
    octave: > grid "on"
    octave: > plot(t, y)
    octave: > closeplot

□ 同一次元オブザーバとフィードバック制御の併合システム

   システム

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

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

    に対して，-4, -5 に極を持つ同一次元オブザーバ

    z'(t) = A z(t) -k[c z(t) - y(t)] + b u(t) = (A - k c)z(t) + k y(t) + b u(t)
  
    の推定値を用いて，極を -1+2j, -1-2j に配置したフィードバック制御を
    行う．

    octave: > A = [-1, -2; -1, -4]
    octave: > b = [1; 1]
    octave: > c = [1, 0]
    octave: > p1 = poly(A)
    octave: > T = [c.', A.'*c.']*[p1(2), 1; 1, 0]
    octave: > p2 = poly([-4, -5])
    octave: > k = ([p2(3) - p1(3), p2(2) - p1(2)]*inv(T)).'

    octave: > T = [b, A*b]*[p1(2), 1; 1, 0]
    octave: > p2 = poly([-1+2j, -1-2j])
    octave: > f = [p2(3) - p1(3), p2(2) - p1(2)]*inv(T)

    併合システムの応答を確かめる．
    初期値 x(0) = [1; -1], z(0) = [0; 0] とする．

    octave: > AA = [A, -b*f; k*c, A-b*f-k*c]
    octave: > global AA;
    octave: > function xdot = F(x, t)
              global AA;
              xdot = AA*x;
              endfunction
    octave: > x0 = [1; -1; 0; 0]
    octave: > t = linspace(0, 5, 1000);
    octave: > y = lsode("F", x0, t);
    octave: > grid "on"
    octave: > plot(t, y)
    octave: > closeplot