MATLABシンボリックでサーボ系設計

シンボリックでいろいろ作ったコードを保存するコーナーです。

＊ コードの内容
プラントのABCD行列が数値で与えられた時に
状態フィードバックのF行列と状態推定器のK行列，1型積分器の係数を決定するところまで
半自動化しました。
キモはこれです。

• 特性方程式det|sI-(A-BF)|の極を自由に配置するFを自動で決定する(係数比較

Fの計算方法は僕の知る限り
係数比較による力技か、変換行列を使うものかの2つに分かれますが
折角なので他の所でも使えるような直接的な手法にして実装しました。
6/23　追記：placeコマンドを使えばこういうことに頭を悩ませなくてもいいですね。ハイ。

中身としては

• solveにおいて作成するeqnsはベクトルor行列形式でもいい

• detの特性方程式を計算するcharpoly関数を確認した

って所の2点を新たに勉強したということになります。
後輩のレポートを手伝っていたら出来上がったので残しておくことにします。

これがsimulinkのブロック線図。
良い復習になりました。

でもよく見たらこれどのブロックがどの変数か見えなくなってるから意味わからんなｗ

重根極配置ではオーバーシュートが避けられず
極の周波数をかなり落として応答性が低い結果になってしまいました。

それを吟味するのは今回の趣旨じゃないのでコードがキモイのも含めてご勘弁。

%% プラント
A = [0 1 0; 0 0 1; 9 7 -1]
B = [0 ; 0 ; 1]
C = [ 1 1 0]

%% 1型サーボ系にするためのフィードバック設計
Abar = [A  [0; 0; 0];-C [0]]
Bbar = [B ;0]

% syms f1;
% syms f2;
% syms f3;
% syms f4;
%Fbar = [f1 f2 f3 f4]

Fbar = sym('f',[1 4])


abf = sym(Abar)-sym(Bbar)*Fbar

% snum = simplify(eig(abf))
Seq = charpoly(abf);

%% pole placement
freq = 0.5;
omega = 2 * pi * freq;
% size(Seq)
pole = [1,4*omega, 6 *omega^2,4*omega^3, omega^4];

%% solve equation
eqns = pole == Seq;
Ans = solve(eqns);

%% outputs
K1 = double(- Ans.f4)
f = double([Ans.f1,Ans.f2,Ans.f3])

%% 同一次元オブザーバー

Kmat = sym('K',[3; 1])

Seq2 = charpoly(A - Kmat*C);
 
%% pole placement2
freq2 = 25;
omega2 = 2 * pi * freq2;

syms s;
pole2 = fliplr(double(coeffs((s+omega2)^size(A,1),s)));
% pole2 = [1,3*omega2, 3 *omega2^2,1*omega^3];

%% solve
eqns2 = pole2 == Seq2;
Ans2 = solve(eqns2);
K = double([Ans2.K1;Ans2.K2;Ans2.K3])

AKC =double( A - K*C )