Maxima で見る現代制御
□ 伝達関数の計算
システム
x'(t) = A x(t) + b u(t)
y(t) = c x(t)
A = [1 , 3; -1, 2], b = [0; 1], c = [2, 1]
の伝達関数を求める．
(%i ) A: matrix([1, 3],[-1, 2]);
(%i ) b: matrix([0],[1]);
(%i ) c: matrix([2, 1]);
(%i ) G: ratsimp(c.invert(s*ident(2)-A).b);
□ 最小実現
伝達関数が
2
s - 1
-----------------
3 2
s + s + 2 s + 1
で与えられる系の最小実現を求める．
(%i ) G: (s^2-1)/(s^3 + s^2 + 2*s + 1);
(%i ) hd: hipow(denom(G), s) - 1;
(%i ) A: addrow(addcol(zeromatrix(hd,1), ident(hd)), makelist(-coeff(denom(G),s,i),i,0,hd));
(%i ) b: addrow(zeromatrix(hd,1), [1]);
(%i ) hn: hipow(num(G), s);
(%i ) c: matrix(makelist(coeff(num(G),s,i),i,0,hd));
□ 可制御性・可観測性
システム
x'(t) = A x(t) + b u(t)
y(t) = c x(t)
A = [-2 , -1; -1, -2], b = [1; 1], c = [1, 0]
の可制御性と，可観測性を調べる．
(%i ) A: matrix([-2, -1],[-1, -2]);
(%i ) b: matrix([1],[1]);
(%i ) c: matrix([1, 0]);
(%i ) Mc: transpose(matrix(transpose(b)[1], transpose(A.b)[1]));
または
(%i ) Mc: addcol(b, A.b);
(%i ) determinant(Mc);
Mc の行列式が 0 なので可制御ではない．
(%i ) Mo: matrix(c[1], (c.A)[1]);
または
(%i ) Mo: addrow(c, c.A);
(%i ) determinant(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]
を対角化する．
(%i ) A: matrix([-2, -1], [-1, -2]);
(%i ) b: matrix([1],[1]);
(%i ) c: matrix([1, 0]);
(%i ) load(eigen);
(%i ) E: eigenvectors(A);
(%i ) T: transpose(matrix(E[2], E[3]));
(%i ) Abar: invert(T).A.T;
(%i ) bbar: invert(T).b;
(%i ) 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]
を可制御正準形に変換する．
(%i ) A: matrix([-1, 2], [-1, -4]);
(%i ) b: matrix([1],[1]);
(%i ) c: matrix([1, 0]);
(%i ) load(nchrpl);
(%i ) p: NCHARPOLY(A, s);
特性多項式が s^2 + 5 s + 6 となることがわかる．なお，NCHARPOLY は大文字
で認識されるバージョン(5.9.1)と小文字で認識されるバージョン(5.10.0)がある
ようなので，注意が必要である．
NCHARPOLY の代わりに charpoly も使えるが，符号を反転させる必要がある．
(%i ) p: expand(-charpoly(A,s));
(%i ) Mc: addcol(b, A.b);
(%i ) S: matrix([coeff(p, s, 1), 1],[1, 0]);
(%i ) T: Mc.S;
(%i ) Abar: invert(T).A.T;
(%i ) bbar: invert(T).b;
(%i ) 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]
の安定性を判別する．
(%i ) A: matrix([-1, 2],[-1, -4]);
(%i ) load(nchrpl);
(%i ) p: NCHARPOLY(A, s);
(%i ) solve(p = 0, s);
A の特性方程式の根 r の実部が全て負なので，システムは安定である．
あるいは，行列 A の固有値を次のように求めても良い．
(%i ) load(eigen);
(%i ) eigenvalues(A);
得られた式が複雑すぎて判定が難しい場合には float を用いて実部の数値解を求め
ると良い．
(%i ) v: solve(p = 0, s);
(%i ) float(realpart(v));
□ 状態フィードバック制御の極配置
システム
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 を求める．
(%i ) A: matrix([-1, -2],[-1, -4]);
(%i ) b: matrix([1],[1]);
(%i ) c: matrix([1, 0]);
(%i ) load(nchrpl);
(%i ) p1: NCHARPOLY(A, s);
(%i ) Mc: addcol(b, A.b);
(%i ) S: matrix([coeff(p1,s,1), 1],[1, 0]);
(%i ) T: Mc.S;
(%i ) p2: expand((s - (-1+2*%i))*(s - (-1-2*%i)));
(%i ) f: [coeff(p2,s,0) - coeff(p1,s,0), coeff(p2,s,1) - coeff(p1,s,1)].invert(T);
初期値 x(0) = [1; -1] の場合の応答を確かめる．
(%i ) x: [x1(t), x2(t)];
(%i ) atvalue(x1(t), t=0, 1);
(%i ) atvalue(x2(t), t=0, -1);
(%i ) sol: desolve(map("=", diff(x,t), transpose((A-b.f).x)[1]), x);
(%i ) plot2d([ev(x1(t),sol), ev(x2(t),sol)], [t,0,5], [GNUPLOT_PREAMBLE, "set grid"]);
□ 同一次元オブザーバの構成
システム
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 を決定する．
(%i ) A: matrix([-1, -2],[-1, -4]);
(%i ) b: matrix([1],[1]);
(%i ) c: matrix([1, 0]);
(%i ) load(nchrpl);
(%i ) p1: NCHARPOLY(A, s);
(%i ) Mo: addrow(c, c.A);
(%i ) S: matrix([coeff(p1,s,1),1],[1,0]);
(%i ) T: transpose(Mo).S;
(%i ) p2: expand((s - (-4))*(s - (-5)));
(%i ) k: transpose([coeff(p2,s,0) - coeff(p1,s,0), coeff(p2,s,1) - coeff(p1,s,1)].invert(T));
オブザーバの推定の様子を確かめる．初期値 x(0) = [0; 0], z(0) = [1; 1]
u(t) = 1 (単位ステップ入力)とする．
(%i ) AA: addrow(addcol(A, zeromatrix(2,2)), addcol(k.c, A-k.c));
(%i ) bb: addrow(b, b);
(%i ) x: [x1(t), x2(t), x3(t), x4(t)];
(%i ) atvalue(x1(t), t=0, 0);
(%i ) atvalue(x2(t), t=0, 0);
(%i ) atvalue(x3(t), t=0, 1);
(%i ) atvalue(x4(t), t=0, 1);
(%i ) sol: desolve(map("=", diff(x,t), transpose(AA.x+bb)[1]), x);
(%i ) plot2d([ev(x1(t),sol), ev(x2(t),sol), ev(x3(t), sol), ev(x4(t), sol)], [t,0,5], [GNUPLOT_PREAMBLE, "set grid"]);
□ 同一次元オブザーバとフィードバック制御の併合システム
システム
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 に配置したフィードバック制御を
行う．
(%i ) A: matrix([-1, -2],[-1, -4]);
(%i ) b: matrix([1],[1]);
(%i ) c: matrix([1, 0]);
(%i ) load(nchrpl);
(%i ) p1: NCHARPOLY(A, s);
(%i ) Mo: addrow(c, c.A);
(%i ) S: matrix([coeff(p1,s,1),1],[1,0]);
(%i ) T: transpose(Mo).S;
(%i ) p2: expand((s - (-4))*(s - (-5)));
(%i ) k: transpose(makelist(coeff(p2,s,i)-coeff(p1,s,i),i,0,1).invert(T));
(%i ) Mc: addcol(b, A.b);
(%i ) T: Mc.S;
(%i ) p2: expand((s - (-1+2*%i))*(s - (-1-2*%i)));
(%i ) f: makelist(coeff(p2,s,i)-coeff(p1,s,i),i,0,1).invert(T);
併合システムの応答を確かめる．
初期値 x(0) = [1; -1], z(0) = [0; 0] とする．
(%i ) AA: addrow(addcol(A, -b.f), addcol(k.c, A-b.f-k.c));
(%i ) x: [x1(t), x2(t), x3(t), x4(t)];
(%i ) atvalue(x1(t), t=0, 1);
(%i ) atvalue(x2(t), t=0, -1);
(%i ) atvalue(x3(t), t=0, 0);
(%i ) atvalue(x4(t), t=0, 0);
(%i ) sol: desolve(map("=", diff(x,t), transpose(AA.x)[1]), x);
(%i ) plot2d([ev(x1(t),sol), ev(x2(t),sol), ev(x3(t), sol), ev(x4(t), sol)], [t,0,5], [GNUPLOT_PREAMBLE, "set grid"]);
以前のページでは AA の計算に誤りがあったため、解が求まっていなかった。
Maxima
naniwa@rbt.his.u-fukui.ac.jp