CT53 分離定理

Home Work

[M] 数学的には、

\displaystyle{(5.43)\quad \frac{d}{dt} \left[\begin{array}{c} x(t)\\ \hat{x}(t)-x(t) \end{array}\right] = \left[\begin{array}{cc} A-BF & -BF \\ 0  & A-HC \end{array}\right] \left[\begin{array}{c} x(t)\\ \hat{x}(t)-x(t) \end{array}\right] }

において、行列A-BFA-HCが対角ブロックに位置し、上三角ブロック行列になるので、2つの行列の固有値がそのまま分離することを言っているのかな。

[C] 上式は、状態オブザーバの推定誤差の振舞いは、状態フィードバックによる制御動作に影響を及ぼさないことを意味しているよ。でも状態推定がもたもたしていると、望ましい状態フィードバックができなくなるので、行列A-BFA-HCの固有値の相対的位置関係に注意が必要となるのかな。Flipped Classroomで考えるようだね。

[P] 制御と観測の話が分離してしまうというのは、なんだか理論がよくできていると思うけど、偶然なのかな、それとも必然なのかな?

Flipped Classroom
[1] 1次系

\displaystyle{(1)\quad \left\{\begin{array}{ll} \dot{x}(t)=ax(t)+bu(t)&(x(t),u(t)\in{\bf R})\\ y(t)=cx(t)&(y(t)\in{\bf R}) \end{array}\right. }

に対するオブザーバベース・コントローラ

\displaystyle{(2)\quad \left\{\begin{array}{ll} \dot{\hat x}(t)=(a-hc)\hat{x}(t)+hy(t)+bu(t)&(\hat{x}(t)\in{\bf R})\\ u(t)=-f\hat{x}(t)& \end{array}\right. }

による閉ループ系の状態方程式は

\displaystyle{(3)\quad \left\{\begin{array}{l} \dot{x}(t)=ax(t)+bu(t)=ax(t)-bf\hat{x}(t) \\ \dot{\hat x}(t)=(a-hc)\hat{x}(t)+hy(t)+bu(t)=(a-hc-bf)\hat{x}(t)+hcx(t) \end{array}\right. }

すなわち

\displaystyle{(4)\quad \left[\begin{array}{c} \dot{x}(t)\\ \dot{\hat x}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} a & -bf \\ hc & a-hc-bf \end{array}\right] }_{A_F} \left[\begin{array}{c} x(t)\\ \hat{x}(t) \end{array}\right] }

[2] 次式が成り立ちます。

\displaystyle{(5)\quad \begin{array}{l} {\rm det}(\lambda I_2-A_F)=(\lambda-a)(\lambda-a+hc+bf)+bfhc\\ =\lambda^2-(2a-hc-bf)\lambda+a(a-hc-bf)+bfhc\\ =\lambda^2-(2a-hc-bf)\lambda+(a-bf)(a-hc)\\ =(\lambda-a+bf)(\lambda-a+hc)\\ =(\lambda+\underbrace{\frac{1}{T_1}}_{\lambda_1})(\lambda+\underbrace{\frac{1}{T_2}}_{\lambda_2}) \end{array} }

ここで、2通りの対応付けが考えられます。

\displaystyle{ \begin{array}{ll} 1^\circ & a-bf=-\frac{1}{T_1}, a-hc=-\frac{1}{T_2}\\ 2^\circ & a-bf=-\frac{1}{T_2}, a-hc=-\frac{1}{T_1} \end{array} }

一方、T_2>T_1の関係より、\frac{1}{T_2}<\frac{1}{T_1}-\frac{1}{T_2}>-\frac{1}{T_1}だから、上の対応付けでは、次が成り立ちます。

\displaystyle{ \begin{array}{ll} 1^\circ & a-hc=-\frac{1}{T_2}>-\frac{1}{T_1}=a-bf\\ 2^\circ & a-bf=-\frac{1}{T_2}>-\frac{1}{T_1}=a-hc \end{array} }

ここで、1^\circは、状態推定の速さ(T_2)が制御動作(T_1)より遅いこと意味し、望ましくありません。しかし、2^\circは、状態推定の速さ(T_1)が制御動作(T_2)より速いこと意味し、こちらの方が望ましいと言えます。

[3] n次系

\displaystyle{(6)\quad \left\{\begin{array}{ll} \dot{x}(t)=Ax(t)+Bu(t)&(x(t)\in{\bf R}^n,u(t)\in{\bf R}^m)\\ y(t)=Cx(t)&(y(t)\in{\bf R}^p) \end{array}\right. }

に対する状態オブザーバ

\displaystyle{(7)\quad \left\{\begin{array}{ll} \dot{z}(t)=\hat{A}z(t)+\hat{B}y(t)+\hat{J}u(t)&(z(t)\in{\bf R}^{n-p})\\ \hat{x}(t)=\hat{C}z(t)+\hat{D}y(t)& \end{array}\right. }

\displaystyle{(8)\quad \exist U\in{\bf R}^{n-p\times n}: \left\{\begin{array}{ll} UA=\hat{A}U+\hat{B}C\\ UB=\hat{J}\\ I_n=\hat{C}U+\hat{D}C\\ \end{array}\right. }

の出力\hat{x}を用いて、制御則

\displaystyle{(7)\quad u(t)=-F\hat{x}(t)=-F(\hat{C}z(t)+\hat{D}y(t)) }

を実施するとします。このときの閉ループ系の状態方程式は

\displaystyle{(9)\quad \left\{\begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)=Ax(t)-BF(\hat{C}z(t)+\hat{D}y(t)) \\ \dot{z}(t)=\hat{A}z(t)+\hat{B}y(t)+\hat{J}u(t) =\hat{A}z(t)+\hat{B}y(t)-\hat{J}F(\hat{C}z(t)+\hat{D}y(t)) \end{array}\right. }

すなわち

\displaystyle{(10)\quad \left[\begin{array}{c} \dot{x}(t)\\ \dot{z}(t) \end{array}\right]= \left[\begin{array}{cc} A-BF\hat{D}C & -BF\hat{C} \\ (\hat{B}-\hat{J}F\hat{D})C  & \hat{A}-\hat{J}F\hat{C} \end{array}\right] \left[\begin{array}{c} x(t)\\ z(t) \end{array}\right] }

ここで、座標変換

\displaystyle{(11)\quad \begin{array}{l} \left[\begin{array}{c} x(t)\\ z(t)-Ux(t) \end{array}\right]= \left[\begin{array}{cc} I_n & 0 \\ -U  & I_{n-p} \end{array}\right] \left[\begin{array}{c} x(t)\\ z(t) \end{array}\right]\\ \Leftrightarrow \left[\begin{array}{c} x(t)\\ z(t) \end{array}\right]= \left[\begin{array}{cc} I_n & 0 \\ U  & I_{n-p} \end{array}\right] \left[\begin{array}{c} x(t)\\ z(t)-Ux(t) \end{array}\right] \end{array} }

を行うと

\displaystyle{(12)\quad \begin{array}{l} \left[\begin{array}{cc} I_n & 0 \\ U  & I_{n-p} \end{array}\right] \frac{d}{dt} \left[\begin{array}{c} x(t)\\ z(t)-Ux(t) \end{array}\right]\\ =\left[\begin{array}{cc} A-BF\hat{D}C & -BF\hat{C} \\ (\hat{B}-\hat{J}F\hat{D})C  & \hat{A}-\hat{J}F\hat{C} \end{array}\right] \left[\begin{array}{cc} I_n & 0 \\ U  & I_{n-p} \end{array}\right] \left[\begin{array}{c} x(t)\\ z(t)-Ux(t) \end{array}\right]\\ =\left[\begin{array}{cc} A-BF\hat{D}C-BF\hat{C}U & -BF\hat{C} \\ (\hat{B}-\hat{J}F\hat{D})C+(\hat{A}-\hat{J}F\hat{C})U  & \hat{A}-\hat{J}F\hat{C} \end{array}\right] \left[\begin{array}{c} x(t)\\ z(t)-Ux(t) \end{array}\right]\\ =\left[\begin{array}{cc} A-BF & -BF\hat{C} \\ \hat{A}U+\hat{B}C-\hat{J}F  & \hat{A}-\hat{J}F\hat{C} \end{array}\right] \left[\begin{array}{c} x(t)\\ z(t)-Ux(t) \end{array}\right] \end{array} }

したがって

\displaystyle{(12)\quad \begin{array}{l} \frac{d}{dt} \left[\begin{array}{c} x(t)\\ z(t)-Ux(t) \end{array}\right]\\ =\left[\begin{array}{cc} I_n & 0 \\ -U  & I_{n-p} \end{array}\right] \left[\begin{array}{cc} A-BF & -BF\hat{C} \\ \hat{A}U+\hat{B}C-\hat{J}F  & \hat{A}-\hat{J}F\hat{C} \end{array}\right] \left[\begin{array}{c} x(t)\\ z(t)-Ux(t) \end{array}\right]\\ = \left[\begin{array}{cc} A-BF & -BF\hat{C} \\ 0  & \hat{A} \end{array}\right] \left[\begin{array}{c} x(t)\\ z(t)-Ux(t) \end{array}\right] \end{array} }

CT52 低次元化

Home Work

[P] もし状態変数の一部が観測変数そのものであれば、状態オブザーバは既知の観測変数も推定しており、無駄なことをしているというのが、低次元化の動機かな?

[M] 出力方程式y=Cxの他に、別の補完的かつ仮想的な観測方程式z=Uxを考えて、これを推定する仕組みを導入したわけだよね。

[C] 状態フィードバックu=-Fxの実施のために、状態オブザーバの出力を用いるのであれば、一気にu=-Fxそのものを推定することも考えられているよ。

Flipped Classroom 5.2
[1]  図5.4の左図は観測変数がそのままでており、右図は観測できない状態変数を推定する様子を示しています。ここで、低次元化されたオブザーバは観測出力からの直達項の影響を受けていることに注意してください。

訂正 図5.4を出力するProgram 52は、最新のものではありませんでした。最新のp52a.m、p52a.sceを入手してください。図5.4はOKです。
MATLAB
%MATLAB: p52a.m (2024/5/1修正)
Tm=0.01; KE=0.05;
x0=[1;10];
t=0:0.001:0.05;
[t1,x1]=ode45(@(t,x) [0 1; 0 -1/Tm]*x+[0;1/Tm/KE],t,x0);
%-----Part 1
alpha=1.2
L=-(1-alpha)/Tm;
%[t2,x2]=ode45(@(t,x) (-1/Tm-L)*x,t,-x0);
[t2,x2]=ode45(@(t,x) (-1/Tm-L)*x,t,-[-L 1]*x0);
%-----Part 2
alpha=1.8
L=-(1-alpha)/Tm;
%[t3,x3]=ode45(@(t,x) (-1/Tm-L)*x,t,-x0);
[t3,x3]=ode45(@(t,x) (-1/Tm-L)*x,t,-[-L 1]*x0);
%-----
subplot(121), plot(t,x1(:,1)), hold on
subplot(121), plot(t,x1(:,1),'-.')
subplot(121), plot(t,x1(:,1),'--')
grid, axis([0 0.05 0 2])
legend('theta','alpha=1.2','alpha=1.8');
xlabel('t'), ylabel('theta(t)')
subplot(122), plot(t,x1(:,2)), hold on
%subplot(122), plot(t,x1(:,2)+x2(:,2),'-.')
%subplot(122), plot(t,x1(:,2)+x3(:,2),'--')
subplot(122), plot(t,x1(:,2)+x2,'-.')
subplot(122), plot(t,x1(:,2)+x3,'--')
%grid, axis([0 0.05 0 30])
grid, axis([0 0.05 0 80])
legend('omega','alpha=1.2','alpha=1.8');
xlabel('t'), ylabel('omega(t)')
SCILAB
//SCILAB: p52a.sce (2024/5/1修正)
Tm=0.01; KE=0.05;
function dx=f1(t,x),dx=[0 1; 0 -1/Tm]*x+[0;1/Tm/KE], endfunction
x0=[1;10];
t0=0; t=0:0.001:0.05;
x1=ode(x0,t0,t,f1);
//-----Part 1
alpha=1.2
L=-(1-alpha)/Tm;
function dx=f2(t,x),dx=(-1/Tm-L)*x, endfunction
//x2=ode(-x0,t0,t,f2);
x2=ode(-[-L 1]*x0,t0,t,f2);
//-----Part 2
alpha=1.8
L=-(1-alpha)/Tm;
//x3=ode(-x0,t0,t,f2);
x3=ode(-[-L 1]*x0,t0,t,f2);
//-----
subplot(121),plot(t,x1(1,:))
subplot(121),plot(t,x1(1,:),'-.')
subplot(121),plot(t,x1(1,:),'--')
mtlb_grid, mtlb_axis([0 0.05 0 2])
legend(['theta';'alpha=1.2';'alpha=1.8']);
xlabel('t'), ylabel('theta(t)')
subplot(122),plot(t,x1(2,:))
//subplot(122),plot(t,x1(2,:)+x2(2,:),'-.')
//subplot(122),plot(t,x1(2,:)+x3(2,:),'--')
subplot(122),plot(t,x1(2,:)+x2,'-.')
subplot(122),plot(t,x1(2,:)+x3,'--')
//mtlb_grid, mtlb_axis([0 0.05 0 30])
mtlb_grid, mtlb_axis([0 0.05 0 80])
xlabel('t'), ylabel('omega(t)')
legend(['omega';'alpha=1.2';'alpha=1.8']);

[2] (5.23)と(5.24)はそのままで、(5.26)が次式となります。

\displaystyle{(1)\quad \hat{C}U+\hat{D}C=K }

CT51 状態オブザーバ

Home Work 5.1

[C] 状態オブザーバは、当時学生だったルーエンバーガーが企業研修相当の機会で発明したものだと聞いたことがあるよ。実際、wikipediaには、次の記述があるね。

In his dissertation Luenberger introduced new methods for construction of state observers. The celebrated Luenberger observer is named after him.

[P] えー、すごいね。どうやって、見つけたのだろうか?やはり、推定誤差をネガティブ・フィードバックすれば、その誤差は小さくなるはずだと思ったのだろうね?

[M] 状態フィードバックによる安定化問題は、行列ABを用いて、可安定性や可制御性として論じられるけど、状態オブザーバの構成問題は、行列ACを用いて、可検出性や可観測性として論じられているね。何せセンサの不足を補えるのだから、とてつもなく便利だよね。

Flipped Classroom 5.1
[1]  (5.18)の左辺を計算すると

\displaystyle{(1)\quad A-HC=\left[\begin{array}{cc}0 & 1\\0 & -\frac{1}{T_m}\end{array}\right] -\left[\begin{array}{c}h_1\\h_2\end{array}\right] \left[\begin{array}{cc}1 & 0\end{array}\right] =\left[\begin{array}{cc}-h_1 & 1\\-h_2 & -\frac{1}{T_m}\end{array}\right] }

だから

\displaystyle{(2)\quad \begin{array}{ll} \det(\lambda I_n-A+HC)=\det\left[\begin{array}{cc}\lambda+h_1 & -1\\h_2 & \lambda+\frac{1}{T_m}\end{array}\right]\\ =(\lambda+h_1)(\lambda+\frac{1}{T_m})+h_2 =\lambda^2+(h_1+\frac{1}{T_m})\lambda+\frac{h_1}{T_m}+h_2 \end{array} }

これが右辺

\displaystyle{(3)\quad (\lambda+\alpha\frac{1}{T_m})^2=\lambda^2+2\alpha\frac{1}{T_m}\lambda+\frac{\alpha^2}{T_m^2} }

に等しいので

\displaystyle{(2)\quad \begin{array}{lll} h_1+\frac{1}{T_m}=2\alpha\frac{1}{T_m} &\Rightarrow & h_1=\frac{2\alpha-1}{T_m}\\ \frac{h_1}{T_m}+h_2=\frac{\alpha^2}{T_m^2} &\Rightarrow & h_2=\frac{\alpha^2-(2\alpha-1)}{T_m^2}=\frac{(\alpha-1)^2}{T_m^2} \end{array} }

[2] \alphaの値によって状態推定の様子がどう違うを観察し、どちらの値が望ましいかをどう判定すればよいかを考えてみてください。

(図5.3左図は縦軸の単位はradですが、Program 51ではdegに変換しています。配布するプログラムは、図5.3をそのまま出力します。)

CT43 補遺4

[C] 状態フィードバックに関する条件として次を学んだことになるね。

【可安定性の定義とその等価な条件】

定義DS: 状態フィードバックにより安定化可能

条件S1: {\rm rank}\, \left[\begin{array}{cc} B & A-\lambda I_n \end{array}\right] =n (\lambdaAのすべての不安定固有値)

条件S2: B^Tw=0,\ A^Tw=\lambda w \Rightarrow w=0 (\lambdaAのすべての不安定固有値)

【可制御性の定義とその等価な条件】

定義DC: 任意初期状態を,任意有限時間内に,任意状態に移動可能

条件C1: \displaystyle{\int_0^t \exp(A\tau)BB^T\exp(A^T\tau)\,d\tau>0 \quad (\forall t>0)}

条件C2: \boxed{{\rm rank}\, \left[\begin{array}{cccc} B & AB & \cdots & A^{n-1}B \end{array}\right]=n}

条件C3: Fを選んで,A-BFの固有値を任意に設定可能

条件C4: \boxed{{\rm rank}\, \left[\begin{array}{cc} B & A-\lambda I_n \end{array}\right] =n} (\lambdaAのすべての固有値)

条件C5: B^Tw=0,\ A^Tw=\lambda w \Rightarrow w=0 (\lambdaAのすべての固有値)

[P] 可制御性は、敢えて言えば、瞬間移動が可能であると言えるのかな?実際にはありえないから、線形理論の限界でもあるね。

[M] これらの証明については、次を参照するらしいよ。

可制御性

(詳しい説明は、あとで行う予定です。)

CT42 いつ安定化できるか

Home Work 4.2

[P] 手のひらに2本の棒を倒立させる制御は考えたことがなかったけど、とても興味深いな。全く同じ長さではできなけど、長さが違えばできるらしいよ。

[M] 安定化できるかどうかのチェックを、行列ABを用いてできるなんて、なんとすばらしい!

[C] いわゆる可安定性とか可制御性が成り立つかどうかをチェックすれば良いわけだ。これらがどう違うかを理解したり、様々な条件の等価性については、補遺で説明するらしいけど、まだかな?

Flipped Classroom 4.2
[1] 可安定性を数値計算で調べる場合は、条件S1を用います。

[2] 可制御性を数値計算で調べる場合は、条件C4を用います。

[3] 条件S1や条件C4が成り立たないことを確認してください。

[4] 初期値をいろいろ変えて、シミュレーションしてみてください。

CT41 状態フィードバック

Home Work 4.1

[M] 数学的には、行列Aが行列A-BFに変わるので、不安定固有値を安定固有値に移動させる可能性さえあるということかな。

[P] 物理的には、ばね定数や摩擦係数を変えることになって、デバイスを変えたわけではないのにダイナミックスが変わってしまい、ちょっと不思議な気がするね。

[C] ある高名な研究者の解説には「状態フィードバックは最強である」という表現があるらしいよ。これはすべての状態変数の動きを監視していて、変化があれば即座にアクチュエータにフィードバックできれば、これ以上の制御動作はないという意味らしいよ。

Flipped Classroom 4.1
[1]  プログラムを解析して、どのような固有値を設定しようとしているか、また正しい状態フィードバックが求まっているか調べてください。また、自分で指定した固有値に設定する状態フィードバックゲインを求めてみてください。

[2]  多入力の場合は、閉ループ系の固有値を指定しただけでは、状態フォードバックは一意の定まらなことを、しっかり認識してください。たとえば、適当な初期値を与えて応答が異なることを確認してみてください。

CT34 補遺3

安定判別

[C] 漸近安定性の等価な条件として次を学んだことになるね。

【漸近安定性の定義とその等価な条件】

定義DA: \forall x(0)\ne 0: x(t)=\exp(At)x(0)\rightarrow 0\quad(t\rightarrow\infty)

条件A0: \exp(At)\rightarrow 0_{n\times n}\quad(t\rightarrow\infty)

条件A1: {\rm Re}(\lambda_i(A))<0\quad(i=1,\cdots,n)

[P] なぜ行列Aだけを用いるのか、その理由を物理的に説明できま~す。

[M] 2次系の場合は分かったからは、あとは一般の場合の証明ができるかどうかだ。そのためには、まず実Jordan標準形を理解する必要があり、次が参考になるらしいよ。

実Jordan標準形

その上で、実Jordan標準形の行列指数関数を考えるのだそうだ。次が参考になるらしいよ。

行列指数関数

(以上、とりあえず手がかりとなる材料を示しておきます。)

CT33 パッシブ制御

Home Work 3.3

[P] 日産の技術紹介でスカイフック制御のビデオを見たことがあるよ。油圧で減衰係数を200通りも変えられると説明されていたね。これにより、高速運転のときはしっかりした足回りにしたり、低速で段差を乗り越えるような場合は柔らかく受け止めるようにできるらしいよ。高度なアクティブ制御の一例かな。

[M] 状況に応じて制御則が変わるのはすごいと思うが、閉ループ系はいつも漸近安定と言えるのだろうか?

[C] ゲインスケジューリング制御の授業では、減衰係数を可変にする制御系設計がレポート課題だったな。確か閉ループ系の安定性が議論されていたと思うけど。

Flipped Classroom 3.3
[1]  行列A=\left[\begin{array}{cc} 0 & 1\\ -\frac{K}{M} & -\frac{D}{M} \end{array}\right] の特性方程式は

\displaystyle{(1)\quad {\rm det}(\lambda I_2-A)={\rm det}\left[\begin{array}{cc} \lambda & -1\\ \frac{K}{M} & \lambda+\frac{D}{M} \end{array}\right]=\lambda^2+\frac{D}{M}\lambda+\frac{K}{M}=0 }

行列Aの固有値は

\displaystyle{(2)\quad \lambda=\frac{1}{2}(-\frac{D}{M}\pm\sqrt{\frac{D^2}{M^2}-\frac{K}{M}}) }

となって、必ず複素左半平面にあることがわかります。

[2] 応答の評価は次の通りです。速応性と静粛性はトレードオフの関係にあります。妥協点はD=\sqrt{2}でしょうか?

D=0.02のとき、速応性は〇、静粛性は×
D=\sqrt{2}のとき、速応性は△、静粛性は△
D=2のとき、速応性はx、静粛性は〇

CT32 安定判別

Home Work 3.2

[C] 安定判別といえば、古典制御ではナイキストの安定定理、ラウス・フルビッツの判定法を思い出すね。それらの手続きについては習ったけど、特にラウス・フルビッツの判定法の証明は習ったかな?

[M] 数学的には行列Aの固有値がすべて複素左半平面にあるかどうかだけどね。

[P] 固有値とか、固有ベクトルとか、線形代数で習ったけど、Jordan標準形の導出などちょっとついていけなかった。ましてや、これらの物理的な意味など説明できるのかな?

Flipped Classroom 3.2
[1],[2],[3],[4],[5]  次のサイトを参照してください。

行列指数関数

[6] 不安定な固有値には赤色をつけてみてください。

CT31 漸近安定性

Home Work 3.1

[P] 漸近安定性の意味は、平衡状態が乱されたときに復帰できるかどうかなので、物理的には明らかだね。

[C] でも状態空間表現におけ3つの行列A,B,Cのうち、どうして行列Aだけしか関係しないのかな?

[M] u\ne0の場合は、状態方程式は微分方程式なので、その解は複雑なはずだよね。なぜu=0とするのかな?

Flipped Classroom 3.1
[1] 漸近安定性は、平衡入力のもとで平常状態の振舞い(状態方程式の解)に関係しています。線形状態方程式において平衡入力はu=0で表されるから、行列Aだけしか関係しません。この説明を行なえる人は少ないようです。

[2] 平衡状態はx=0で表されるから、零に漸近しているのはa=-5の場合だけです。テキストはモノクロ印刷なので、必ずプログラムを実行してください。