粗大メモ置き場

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

極配置の実装:アッカーマン法,EigenvalueStructure法のMATLAB実装

序論・参考文献

制御工学の授業で一度は計算させられる極配置ですが MATLABにはplaceなる関数があって,実務を始めてからは全く手計算をしなくなります。

よくある疑問とエラー

授業では忘却の彼方にありますがB行列のrankより高い極配置はできないというエラーは誰しも見たことがあるかと思います。 以下のような現象です。

A = [1 2;3 4];
B = [1; -2];
pole = [-3,-3];
place(A,B,pole)

Bのランクは1ですので,このコマンドを実行すると以下のエラーを起こします。

エラー: place (line 78)
"place" コマンドは、rank(B) より大きい多重度を持つ極を配置できません。

そういえば中身どうなっているんだ?Cとかで実装するとどうなるんや?という疑問がありませんか。 MATLABのは論文読んでないのでまだ実装しませんが, 本稿で紹介するアルゴリズムを用いた場合はこの場合でも,ほぼ極配置ができます。

  -3.0000 + 0.0000i
  -3.0000 - 0.0000i

ひょっとしたらちょっと誤差があるのかもしれない。

アッカーマン法

SISO系限定ですがアッカーマン法であればこのような状況でも対応できます。 なお、MATLABではackerという関数名で実装されています。(昔はReferenceがあった気がするのですが。。。?)

アッカーマン(Ackermann)はこの人?某ヒロインが検索汚染してきます。

アッカーマン法の詳細な記述についてはなんとwikipediaが一番スッキリしていて見やすくなってます。(英語)

原著の論文はこれっぽいです。

J. Ackermann. Der Entwurf linearer Regelungssysteme im Zustandsraum.
Regelungstechnik und Prozessdatenverarbeitung, 7:297–300,
1972.

参考資料

どこかの授業スライドですがこれが具体例込でわかりやすいです。 また,参考になりそうなスライドがありました。(要DL,MIMOや不可制御系の時にどうする的なのが丁寧に書いてあったスライド。わかりやすい5-10分で読める。)

Outline

  • Step1. 求めたい極を展開して特性方程式にする
  • Step2. 可制御性行列を求める
  • Step3. 以下の通りに計算する(雑)

f:id:ossyaritoori:20180516102843p:plain 添付は小林H23スライドより拝借。

MATLABコード

多項式展開が思ったより実装が面倒そうだったので一般化してないですが例えば2次元の場合は次の通り。

function F = placeAckermann2d(A,B,pole)
p = - pole;
d1 = p(1)+p(2);
d0 = p(1)*p(2);
X = A^2 + d1*A + d0*eye(2);

Uc = [B A*B];

F = [0 1]/Uc*X;
end

これで冒頭の結果を得ることができます。Bのrankが落ちていても重根極配置できるので試してみてください。 この程度の計算であれば行列と逆行列計算ライブラリさえあればCやPythonでも実装余裕ですね。

MIMOの場合:EigenvalueStructure法

アッカーマン法ではMIMOの時に行列のサイズがおかしくなります。 そこでこの手法です。1981年のこの文献が元

原理としては以下のような感じです。 f:id:ossyaritoori:20180516131828p:plain

Nullspaceの推定にsvgを使うので計算量には不安があるが,見事に動作します。 この場合自由度の足りない重根極配置をするとVの正則性が失われるので最初に出たようなエラーが出るということのようです。

function F = placeMIMO2D(A,B,pole)
[u1 w1 v1] = svd([pole(1)*eye(2)-A B])
[u2 w2 v2] = svd([pole(2)*eye(2)-A B])
en = size(w1,2)
VV = [v1(:,3:en) v2(:,3:en)];
F = VV(3:en,:)/VV(1:2,:)
end