読者です 読者をやめる 読者になる 読者になる

粗大メモ置き場

個人用,たまーに来訪者を意識する雑記メモ

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

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

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

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

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


中身としては

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

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


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

f:id:ossyaritoori:20160622025020p:plain
でもよく見たらこれどのブロックがどの変数か見えなくなってるから意味わからんなw


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

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

%% プラント
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 )