「制御技術」カテゴリーアーカイブ

最小次元状態オブザーバ

投稿日時: 投稿者:

最小次元状態オブザーバ…Homework

[1]　制御対象のモデルが

\displaystyle{(1)\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, {\rm rank}C=p)\right. \end{array} }

で与えられるとします。これに対する状態オブザーバは、微分方程式

\displaystyle{(2)\quad \begin{array}{l} \dot{\hat{x}}(t)=(A-HC)\hat{x}(t)+Bu(t)+Hy(t),\ \hat{x}(0)=0\\ A-HC: Stable Matrix \end{array} }

で表され、解ベクトル\hat{x}は、制御対象の状態ベクトルxと同じ次元数nを持ちます。

いま、出力方程式y=Cx（Cのサイズはp\times n、行フルランク）を通して、状態に関する情報が部分的に得られていることに注目します。そこで、もう一つの補完的な出力方程式z=Ux（Uのサイズはn-p\times n、行フルランク）を\left[\begin{array}{c} C\\ U \end{array}\right]が正則となるように適切に選んで

\displaystyle{(3)\quad \left[\begin{array}{c} y(t)\\ z(t) \end{array}\right] = \left[\begin{array}{c} C\\ U \end{array}\right]x(t) \Rightarrow x(t)= \left[\begin{array}{c} C\\ U \end{array}\right]^{-1} \left[\begin{array}{c} y(t)\\ z(t) \end{array}\right] }

とすれば、状態に関する全情報が得られそうです。もちろん

\displaystyle{(4)\quad \hat{z}(t)-\underbrace{Ux(t)}_{z(t)}\rightarrow 0\quad (t\rightarrow\infty) }

を達成する仕組みが必要です。いま、その仕組みを

\displaystyle{(5)\quad \boxed{\left\{\begin{array}{lll} \dot{\hat{z}}(t)=\hat{A}\hat{z}(t)+\hat{B}y(t)+\hat{J}u(t),\ \hat{z}(0)=0\\ \hat{x}(t)=\hat{C}\hat{z}(t)+\hat{D}y(t) \end{array}\right.} }

とします。(5)の第１式から、Uを掛けた(1)の第１式を辺々引くと

\displaystyle{(6)\quad \begin{array}{l} \frac{d}{dt}(\hat{z}(t)-Ux(t))\\ =\hat{A}\hat{z}(t)+\hat{B}y(t)+\hat{J}u(t)-(UAx(t)+UBu(t))\\ =\hat{A}\hat{z}(t)-(UA-\hat{B}C)x(t)+(\hat{J}-UB)u(t) \end{array} }

を得ます。ここで

\displaystyle{(7)\quad \boxed{\hat{A}U+\hat{B}C=UA} }

を仮定し、\hat{J}を

\displaystyle{(8)\quad \boxed{\hat{J}=UB} }

とおけば

\displaystyle{(9)\quad \frac{d}{dt}(\hat{z}(t)-Ux(t))=\hat{A}(\hat{z}(t)-Ux(t)) }

を得ます。もし\hat{A}が安定行列であれば、(4)を達成できます。また

\displaystyle{(10)\quad \boxed{\hat{C}U+\hat{D}C=I_n} }

を仮定すると

\displaystyle{(11)\quad \hat{x}(t)=\hat{C}\hat{z}(t)+\hat{D}y(t) \rightarrow \underbrace{(\hat{C}U+\hat{D}C)}_{I_n}x(t)=x(t)\quad (t\rightarrow\infty) }

を得て、状態推定が行われていることがわかります。状態オブザーバ(2)と比べると、計算機内で解くべき微分方程式の解ベクトルの次元数がnからn-pに低次元化されていることになります。\left[\begin{array}{c} C\\ U \end{array}\right]が正則となるためには、Uの行数はn-p以下にすることはできませんから、最小次元状態オブザーバと呼ばれています。

[2] 以上から、適当なUを選んで、(7)、(8)、(10)を満足し、安定行列\hat{A}をもつ(3)を設計することを考えます。たとえばC=[I_p\ 0]のときはU=[-L\ I_{n-p}]のように選びます（Lのサイズはn-p\times p）。このとき、(7)と(10)をまとめた

\displaystyle{(12)\quad \left[\begin{array}{cc} \hat{A}&\hat{B}\\ \hat{C}&\hat{D} \end{array}\right] = \left[\begin{array}{cc} UA\\ I_n \end{array}\right] \left[\begin{array}{cc} U\\ C \end{array}\right]^{-1} }

において、A=\left[\begin{array}{cc} A_{11}& A_{12}\\ A_{21}& A_{22} \end{array}\right]（A_{11}のサイズはp\times p）のように分割すると

\displaystyle{(13)\quad \left[\begin{array}{cc} -L& I_{n-p} \end{array}\right] \left[\begin{array}{cc} A_{11}& A_{12}\\ A_{21}& A_{22} \end{array}\right] = \left[\begin{array}{cc} -LA_{11}+A_{21}& -LA_{12}+A_{22} \end{array}\right] }

\displaystyle{(14)\quad \left[\begin{array}{cc} -L& I_{n-p}\\ I_p & 0 \end{array}\right]^{-1} = \left[\begin{array}{cc} 0& I_{p}\\ I_{n-p} & L \end{array}\right] }

より、\hat{A}、\hat{B}、\hat{C}、\hat{D}は次式から求められます。

\displaystyle{(15)\quad \boxed{\left\{\begin{array}{lll} \hat{A}=A_{22}-LA_{12}\\ \hat{B}=A_{21}-LA_{11}+\hat{A}L\\ \hat{C}= \left[\begin{array}{cc} 0\\ I_{n-p} \end{array}\right]\\ \hat{D}= \left[\begin{array}{cc} I_{p}\\ L \end{array}\right] \end{array}\right.} }

ここで、Lは\hat{A}が安定行列となるように選ぶ必要があります。そのためには(A,C)が可観測対であればよいのですが、これについてはあとで触れます。また、\hat{J}は(10)から定めます。

● モータの状態方程式

\displaystyle{(16)\quad \left[\begin{array}{c} \dot{\theta}(t)\\ \dot{\omega}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 0 & -\frac{1}{T_m} \end{array}\right] }_{A} \left[\begin{array}{c} \theta(t)\\ \omega(t) \end{array}\right]+ \underbrace{ \left[\begin{array}{c} 0\\ \frac{1}{K_ET_m} \end{array}\right] }_{B} u(t) }

と出力方程式

\displaystyle{(17)\quad y(t)= \underbrace{ \left[\begin{array}{cc} 1 & 0 \end{array}\right] }_{C} \left[\begin{array}{c} \theta(t)\\ \omega(t) \end{array}\right] }

に対する状態オブザーバの低次元化を行ってみます。ちょうどC=[1\ 0]となっているので、U=[-L\ 1]とおきます。設計式(15)にA_{11}=0、A_{12}=1、A_{21}=0、A_{22}=-\frac{1}{T_m}を代入して

\displaystyle{(18)\quad \left\{\begin{array}{lll} \hat{A}=-\frac{1}{T_m}-L\\ \hat{B}=\hat{A}L\\ \hat{C}= \left[\begin{array}{cc} 0\\ 1 \end{array}\right]\\ \hat{D}= \left[\begin{array}{cc} 1\\ L \end{array}\right]\\ \hat{J}=\frac{1}{K_ET_m} \end{array}\right. }

ここで、Lは\hat{A}=-\alpha\frac{1}{T_m}とするのであればL=-(1-\alpha)\frac{1}{T_m}と決めます。以上から、次の低次元化された状態オブザーバが得られました。

\displaystyle{(19)\quad \left\{\begin{array}{lll} \dot{\hat{z}}(t)=\underbrace{-\frac{\alpha}{T_m}}_{\hat{A}}\hat{z}(t) +\underbrace{\frac{\alpha(1-\alpha)}{T_m^2}}_{\hat{B}}y(t) +\underbrace{\frac{1}{K_ET_m}}_{\hat{J}}u(t),\ \hat{z}(0)=0\\ \hat{x}(t)= \underbrace{ \left[\begin{array}{cc} 0\\ 1 \end{array}\right]}_{\hat{C}}\hat{z}(t) + \underbrace{ \left[\begin{array}{cc} 1\\ L \end{array}\right] }_{\hat{D}}y(t) \end{array}\right. }

●次図はT_m=0.01、K_E=0.05の下で、\alpha=1.2,1,8と変えたとき、状態オブザーバが真の状態を推定する様子をシミュレーションしています。

図1 状態推定の例

演習 A42…Flipped Classroom

次のコードを用いて図1のグラフを描け。

MATLAB 
%a42.m
%-----
 clear all, close all
 Tm=0.01; KE=0.05;
 A=[0 1; 0 -1/Tm]; B=[0;1/Tm/KE]; C=[1 0];
 sys=ss(A,B,eye(2,2),[]); 
 t=0:0.001:0.05;
 u=ones(1,length(t));
 x0=[1;10];
 x=lsim(sys,u,t,x0);
%-----;
 alpha=1.2
 L=-(1-alpha)/Tm;
 U=[-L 1]
 Aob=-1/Tm*alpha;
 Bob=alpha*(1-alpha)/Tm^2;
 Job=1/Tm/KE;
 Cob=[0;1];
 Dob=[1;-(1-alpha)/Tm];
 err1=ss(Aob,[],1,[]);
 xob0=-U*x0;
 e1=initial(err1,xob0,t);  
 w1=[Cob Dob]*[e1'+U*x';x(:,1)']; 
%----- 
 alpha=1.8
 L=-(1-alpha)/Tm;
 U=[-L 1]
 Aob=-1/Tm*alpha;
 Bob=alpha*(1-alpha)/Tm^2;
 Job=1/Tm/KE;
 Cob=[0;1];
 Dob=[1;-(1-alpha)/Tm];
 err2=ss(Aob,[],1,[]);
 xob0=-U*x0;
 e2=initial(err2,xob0,t);
 w2=[Cob Dob]*[e2'+U*x';x(:,1)']; 
%-----  
 subplot(121),plot(t,x(:,1),t,w1(1,:),t,w2(1,:)), 
 legend('theta','apha=1.2','alpha=1.8'), grid
 subplot(122),plot(t,x(:,2),t,w1(2,:),t,w2(2,:)), 
 legend('omega','apha=1.2','alpha=1.8'), grid 
%-----
%eof
SCILAB 
coming soon

Note A42　線形関数オブザーバ

制御対象の状態空間表現

\displaystyle{(1)\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)\right. \end{array} }

に対して、もう一つの状態空間表現

\displaystyle{(2)\quad \boxed{\left\{\begin{array}{lll} \dot{\hat{z}}(t)=\hat{A}\hat{z}(t)+\hat{B}y(t)+\hat{J}u(t),\ \hat{z}(0)=0&(x(t)\in{\bf R}^q)\\ \hat{y}(t)=\hat{C}\hat{z}(t)+\hat{D}y(t)&(\hat{y}(t)\in{\bf R}^{\hat{p}}) \end{array}\right.} }

を考えます。ただし、サイズq\times nの行列Uとサイズm\times nの行列Kに対して

\displaystyle{(3)\quad \boxed{\left\{\begin{array}{lll} \hat{A}: Stable\ Matrix\\ \hat{A}U+\hat{B}C=UA\\ \hat{J}=UB\\ \hat{C}U+\hat{D}C=K \end{array}\right.} }

が満足されるものとします。このとき

\displaystyle{(6)\quad \begin{array}{l} \frac{d}{dt}(\hat{z}(t)-Ux(t))\\ =\hat{A}\hat{z}(t)+\hat{B}y(t)+\hat{J}u(t)-(UAx(t)+UBu(t))\\ =\hat{A}\hat{z}(t)-(UA-\hat{B}C)x(t)+(\hat{J}-UB)u(t)\\ =\hat{A}(\hat{z}(t)-Ux(t))) \end{array} }

から

\displaystyle{(7)\quad \hat{z}(t)-Ux(t)\rightarrow 0\quad (t\rightarrow\infty) }

を得ます。また

\displaystyle{(8)\quad \hat{y}(t)=\hat{C}\hat{z}(t)+\hat{D}y(t) \rightarrow \underbrace{(\hat{C}U+\hat{D}C)}_{K}x(t)=Kx(t)\quad (t\rightarrow\infty) }

を得て、状態の線形関数Kx(t)の推定が行われていることがわかります。

したがって、条件(3)を満足する状態空間表現(2)は、線形関数オブザーバと呼ばれます。特にK=I_nのとき、は最小次元状態オブザーバとなります。

●興味深いのは、状態フィードバック

\displaystyle{(9)\quad u(t)=\underbrace{-F}_{K}x(t) }

を推定させる場合で、線形関数オブザーバの次元数qが最小次元状態オブザーバの次数n-pより小さくできる可能性があります。

たとえば機械系でよく見られる、次の１入力p出力2p次系

\displaystyle{(10)\quad \left\{\begin{array}{l} \left[\begin{array}{c} \dot{y}(t)\\ \ddot{y}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} 0 & I_p \\ A_{21} & A_{22} \end{array}\right] }_{A} \left[\begin{array}{c} y(t)\\ \dot{y}(t) \end{array}\right]+ \underbrace{ \left[\begin{array}{c} 0\\ B_2 \end{array}\right] }_{B} u(t)\\ y(t)= \underbrace{ \left[\begin{array}{cc} I_n & 0 \end{array}\right] }_{C} \left[\begin{array}{c} y(t)\\ \dot{y}(t) \end{array}\right] \end{array}\right. }

および、これに対する状態フィードバック

\displaystyle{(11)\quad u(t)=\underbrace{ \left[\begin{array}{cc} K_1 & K_2 \end{array}\right] }_{K}x(t) }

を考えます。この状態フィードバックを推定する線形関数オブザーバは、(p+1)入力１出力１次系

\displaystyle{(12)\quad \left\{\begin{array}{l} \dot{\hat{z}}(t)= \underbrace{ \lambda }_{\hat{A}} \hat{z}(t)+ \underbrace{ K_2(A_{21}+\lambda L) }_{\hat{B}} y(t)+ \underbrace{ K_2B_2 }_{\hat{J}} u(t)\quad(\lambda<0)\\ \hat{y}(t)= \underbrace{1}_{\hat{C}}\hat{z}(t)+ \underbrace{ (K_1+K_2L) }_{\hat{D}} y(t) \end{array}\right. }

として得られます。実際、LとUを

\displaystyle{(13)\quad \begin{array}{l} L=A_{22}-\lambda I_p,\ U=K_2 \left[\begin{array}{cc} -L & I_p \end{array}\right] \end{array}}

と選べば、(3)の第２、３、４式を満足していることを次のように確かめることができます。

\displaystyle{(14)\quad \begin{array}{l} \underbrace{ \lambda K_2\left[\begin{array}{cc} -L & I_p \end{array}\right] + K_2(A_{21}+\lambda L) \left[\begin{array}{cc} I_n & 0 \end{array}\right] }_{\hat{A}U+\hat{B}C} = \underbrace{ K_2\left[\begin{array}{cc} -L & I_p \end{array}\right] \left[\begin{array}{cc} 0 & I_p \\ A_{21} & A_{22} \end{array}\right] }_{UA}\\ \underbrace{ K_2B_2 }_{\hat{J}} = \underbrace{ K_2\left[\begin{array}{cc} -L & I_p \end{array}\right] \left[\begin{array}{c} 0\\ B_2 \end{array}\right] }_{UB}\\ \underbrace{ 1\cdot K_2\left[\begin{array}{cc} -L & I_p \end{array}\right] + (K_1+K_2L) \left[\begin{array}{cc} I_n & 0 \end{array}\right] }_{\hat{C}U+\hat{D}C} = \underbrace{ \left[\begin{array}{cc} K_1 & K_2 \end{array}\right] }_{K} \end{array}}

●状態フィードバックを推定する線形関数オブザーバ（したがって最小次元状態オブザーバ）に対しても分離定理を示すことができます。すなわち制御対象のモデル(1)に対して、K=-Fの場合の(3)を満たす線形関数オブザーバ(2)による閉ループ系は

\displaystyle{(15)\quad e(t)=z(t)-Ux(t) }

を定義すると、次式となります。

\displaystyle{(16)\quad \boxed{\left[\begin{array}{c} \dot{x}(t) \\ \dot{e}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A-BF & -B\hat{C} \\ 0 & \hat{A} \end{array}\right] }_{A_{CL}} \left[\begin{array}{c} x(t) \\ e(t) \end{array}\right]} }

この第１式は(1)の状態方程式のu(t)に、(2)の出力方程式の\hat{y}(t)を代入して

\displaystyle{(17)\quad \begin{array}{l} \dot{x}(t)=Ax(t)+B(\hat{C}\hat{z}(t)+\hat{D}y(t))\\ =Ax(t)+B(\hat{C}\hat{z}(t)+\underbrace{\hat{D}C}_{-F-\hat{C}U}x(t))\\ =(A-BF)x(t)+B\hat{C}\underbrace{(\hat{z}(t)-Ux(t))}_{e(t)} \end{array} }

のように得られます。また、第２式は(6)そのものです。

状態オブザーバ

投稿日時: 投稿者:

状態オブザーバ…Homework

[1]　状態フィードバックは、すべての状態変数の変化を捉えてフィードバックするので、この意味で最強の手段ですが、一般には絵に描いた餅ですと言ったら驚くでしょうか？これは実際にはすべての状態変数にセンサを配置できるわけではないからです。でもちょっと待ってください。状態方程式が分かっているのですから、これを計算機内で解いてやれば良いのではないでしょうか？そのアイデアを図1に示します。

図1 状態オブザーバの考え方（１）

●いま制御対象の状態空間表現を

\displaystyle{(1)\quad \left\{\begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t),\ x(0)\ne0\\ y(t)=Cx(t) \end{array}\right. }

とします。これを計算機で解きますが、その状態と出力を区別して

\displaystyle{(2)\quad \left\{\begin{array}{lll} \dot{\hat{x}}(t)=A\hat{x}(t)+Bu(t),\ \hat{x}(0)=0\\ \hat{y}(t)=C\hat{x}(t) \end{array}\right. }

で表します。(2)から(1)の第１式どうしを辺々引くと

\displaystyle{(3)\quad \frac{d}{dt}(\hat{x}(t)-x(t))=A(\hat{x}(t)-x(t)),\ \hat{x}(0)-x(0)\ne0 }

を得ます。ここでAが安定行列でない限り

\displaystyle{(4)\quad \hat{x}(t)-x(t)\rightarrow 0\quad (t\rightarrow\infty) }

とはなりません。

●そこで、図2に示すような工夫が行われました。

図2 状態オブザーバの考え方（２）

(2)の状態方程式に次のような補正項を加えます。

\displaystyle{(5)\quad \dot{\hat{x}}(t)=A\hat{x}(t)+Bu(t)-H(C\hat{x}(t)-y(t)),\ \hat{x}(0)=0 }

すなわち

\displaystyle{(6)\quad \boxed{\dot{\hat{x}}(t)=(A-HC)\hat{x}(t)+Bu(t)+Hy(t),\ \hat{x}(0)=0} }

ここで、Hは適合するサイズを持つ行列です。(6)から(1)の第１式どうしを辺々引くと

\displaystyle{(7)\quad \frac{d}{dt}(\hat{x}(t)-x(t))=(A-HC)(\hat{x}(t)-x(t)),\ \hat{x}(0)\ne0 }

を得ます。ここでA-HCが安定行列であるとすると

\displaystyle{(8)\quad \hat{x}(t)-x(t)\rightarrow 0\quad (t\rightarrow\infty) }

が成り立ちます。すなわち(6)の状態\hat{x}(t)は、制御対象の状態x(t)を漸近的に推定することができます。A-HCが安定行列である(6)を状態オブザーバ、行列Hをオブザーバゲインと呼びます。

さて、任意の正方行列Xに対して\det X=\det X^Tであることに注意して

\displaystyle{(9)\quad \left.\begin{array}{lll} \det(\lambda I_n-A+HC)&=&\det(\lambda I_n-A+HC)^T\\ &=&\det(\lambda I_n-A^T+C^TH^T) \end{array}\right. }

が成り立ちます。これからA-HCを安定行列とするためには、状態方程式

\displaystyle{(10)\quad \boxed{\dot{w}(t)=A^Tw(t)+C^Tv(t)} }

に対する状態フィードバック

\displaystyle{(11)\quad \boxed{v(t)=-H^Tw(t)} }

による閉ループ系が

\displaystyle{(12)\quad \dot{w}(t)=(A^T-C^TH^T)w(t) }

となるので、これを安定化する(12)を求めて、オブザーバゲインHを決定すればよいことが分かります。

状態方程式(10)をもつ物理系は存在せず、双対システムと呼ばれる仮想的なシステムであり、このシステムに対する状態フィードバック(11)の設計を通して、状態オブザーバの設計ができることになります。

[2]　モータの状態方程式

\displaystyle{(13)\quad \left[\begin{array}{c} \dot{\theta}(t)\\ \dot{\omega}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 0 & -\frac{1}{T_m} \end{array}\right] }_{A} \left[\begin{array}{c} \theta(t)\\ \omega(t) \end{array}\right]+ \underbrace{ \left[\begin{array}{c} 0\\ \frac{1}{K_ET_m} \end{array}\right] }_{B} u(t) }

と出力方程式

\displaystyle{(14)\quad y(t)= \underbrace{ \left[\begin{array}{cc} 1 & 0 \end{array}\right] }_{C} \left[\begin{array}{c} \theta(t)\\ \omega(t) \end{array}\right] }

に対する状態オブザーバ

\displaystyle{(15)\quad \left[\begin{array}{c} \dot{\hat{\theta}}(t)\\ \dot{\hat{\omega}}(t) \end{array}\right] =(A-HC) \left[\begin{array}{c} \hat{\theta}(t)\\ \hat{\omega}(t) \end{array}\right] +Bu(t)+Hy(t) }

を、A-HCの固有値が

\displaystyle{(16)\quad \lambda'_1=\lambda'_2=-\alpha\frac{1}{T_m}\ (\alpha>1) }

となるようにオブザーバゲインHを求めます。ここでは

\displaystyle{(17)\quad \det(\lambda I_2-A+HC)=(\lambda+\alpha\frac{1}{T_m})^2 }

における係数比較から直接求めてみます。

\displaystyle{(18)\quad A-HC= \left[\begin{array}{cc} 0 & 1 \\ 0 & -\frac{1}{T_m} \end{array}\right] - \left[\begin{array}{cc} 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] }

を(17)に代入して

\displaystyle{(19)\quad \begin{array}{l} \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+h_1\frac{1}{T_m}+h_2\\ =\lambda^2+2\alpha\frac{1}{T_m}\lambda+(\alpha\frac{1}{T_m})^2\\ \end{array} }

ここで係数比較により次式を得ます。

\displaystyle{(20)\quad \begin{array}{l} h_1=\frac{2\alpha-1}{T_m}\\ h_2=(\alpha\frac{1}{T_m})^2-h_1\frac{1}{T_m}=(\alpha\frac{1}{T_m})^2-\frac{2\alpha-1}{T_m}\frac{1}{T_m}=\frac{(\alpha-1)^2}{T_m^2} \end{array} }

これらから定まるオブザーバゲインHを用いて、状態オブザーバは次式のように構成されます。

\displaystyle{(21)\quad \begin{array}{l} \left[\begin{array}{c} \dot{\hat{\theta}}(t)\\ \dot{\hat{\omega}}(t) \end{array}\right] =\underbrace{ \left[\begin{array}{cc} -\frac{2\alpha-1}{T_m} & 1 \\ -\frac{(\alpha-1)^2}{T_m^2} & -\frac{1}{T_m} \end{array}\right] }_{A-HC} \left[\begin{array}{c} \hat{\theta}(t)\\ \hat{\omega}(t) \end{array}\right] +\underbrace{ \left[\begin{array}{c} 0\\ \frac{1}{K_ET_m} \end{array}\right] }_{B}u(t)\\ + \underbrace{\left[\begin{array}{c} \frac{2\alpha-1}{T_m} \\ \frac{(\alpha-1)^2}{T_m^2} \end{array}\right] }_{H}y(t) \end{array} }

●次図はT_m=0.01、K_E=0.05の下で、\alpha=1.2,1,8と変えたとき、状態オブザーバが真の状態を推定する様子をシミュレーションしています。

図3 状態推定の例

演習 A41…Flipped Classroom

次のコードを用いて図３のグラフを描け。

MATLAB 
%a41.m
%-----
 clear all, close all
 Tm=0.01; KE=0.05;
 A=[0 1; 0 -1/Tm]; B=[0;1/Tm/KE]; C=[1 0];
 sys=ss(A,B,eye(2,2),[]);
 t=0:0.001:0.05;
 u=ones(1,length(t));
 x0=[1;10];
 x=lsim(sys,u,t,x0);
%-----;
 alpha=1.2
 p=[-1/Tm*alpha -1/Tm*alpha];
 H=sf(A',C',p); H=H'
%  h1=-(1-2*alpha)/Tm; h2=(1-alpha)^2/Tm^2;
%  H=[h1;h2]
 Aob=A-H*C; 
 err1=ss(Aob,[],eye(2,2),[]);
 xob0=[0;0];
 e1=initial(err1,-x0,t);
%----- 
 alpha=1.8
 p=[-1/Tm*alpha -1/Tm*alpha];
 H=sf(A',C',p); H=H' 
%  h1=-(1-2*alpha)/Tm; h2=(1-alpha)^2/Tm^2;
%  H=[h1;h2]
 Aob=A-H*C; 
 err2=ss(Aob,[],eye(2,2),[]);
 xob0=[0;0];
 e2=initial(err2,-x0,t);
%-----  
 subplot(121),plot(t,x(:,1),t,x(:,1)+e1(:,1),t,x(:,1)+e2(:,1)), 
 legend('theta','apha=1.2','alpha=1.8'), grid
 subplot(122),plot(t,x(:,2),t,x(:,2)+e1(:,2),t,x(:,2)+e2(:,2)), 
 legend('omega','apha=1.2','alpha=1.8'), grid
%----- 
 function F=sf(A,B,p)
   n=size(A,1); r=eig(A);
   a=[1 -r(1)]; for i=2:n, a=conv(a,[1 -r(i)]); end, a=a(n+1:-1:2);
   b=[1 -p(1)]; for i=2:n, b=conv(b,[1 -p(i)]); end, b=b(n+1:-1:2);
   X=zeros(n,n); for i=1:n, X(i,1:n-i+1)=[a(i+1:n) 1]; end
   Y=B; for i=2:n, Y=[B A*Y]; end
   F=real((b-a)/X/Y);
 end
%-----
%eof
SCILAB 
coming soon

[3] n次系

\displaystyle{(1)\quad \left\{\begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t),\ x(0)\ne0\\ y(t)=Cx(t) \end{array}\right. }

に対する状態フィードバックu=-Fxは、状態オブザーバを用いて、次のように実施されます。

\displaystyle{(22)\quad \left\{\begin{array}{l} \dot{\hat{x}}(t)=(A-HC)\hat{x}(t)+Hy(t)+Bu(t)\\ u(t)=-F\hat{x}(t) \end{array}\right. }

すなわち

\displaystyle{(23)\quad \boxed{\left\{\begin{array}{l} \dot{\hat{x}}(t)=(A-HC-BF)\hat{x}(t)+Hy(t)\\ u(t)=-F\hat{x}(t) \end{array}\right.} }

これをオブザーバベースコントローラと呼びます。ここで、A-BFとA-HCは安定行列ですが、A-HC-BFが安定行列と限らないことに注意してください。また、

\displaystyle{(24)\quad e(t)=\hat{x}(t)-x(t) }

とおくとき、オブザーバベースコントローラは次式でも表せることに注意してください。

\displaystyle{(25)\quad \boxed{\left\{\begin{array}{l} \dot{e}(t)=(A-HC)e(t)\\ u(t)=-F(x(t)+e(t)) \end{array}} }

さて、オブザーバベースコントローラ(23)による閉ループ系は次式となります。

\displaystyle{(26)\quad \left[\begin{array}{c} \dot{x}(t) \\ \dot{\hat{x}}(t) \end{array}\right]= \left[\begin{array}{cc} A & -BF \\ HC & A-HC-BF \end{array}\right] \left[\begin{array}{c} x(t) \\ \hat{x}(t) \end{array}\right] }

ここで、座標変換

\displaystyle{(27)\quad \left[\begin{array}{c} x(t) \\ e(t) \end{array}\right]= \left[\begin{array}{cc} I & 0 \\ -I & I \end{array}\right] \left[\begin{array}{c} x(t) \\ \hat{x}(t) \end{array}\right] }

を行うと、

\displaystyle{(28)\quad \left[\begin{array}{cc} I & 0 \\ -I & I \end{array}\right] \left[\begin{array}{cc} A & -BF \\ HC & A-HC-BF \end{array}\right]\\ = \left[\begin{array}{cc} A-BF & -BF \\ 0 & A-HC \end{array}\right] \left[\begin{array}{cc} I & 0 \\ -I & I \end{array}\right] }

に注意して、閉ループ系の状態方程式として、次式を得ます。

\displaystyle{(29)\quad \boxed{\left[\begin{array}{c} \dot{x}(t) \\ \dot{e}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A-BF & -BF \\ 0 & A-HC \end{array}\right] }_{A_{CL}} \left[\begin{array}{c} x(t) \\ e(t) \end{array}\right]} }

したがって、A_{CL}の固有値は、A-BFとA-HCの固有値からなり、オブザーバベースコントローラによる閉ループ系は漸近安定であることが分かります。このことは状態フィードバックゲインFの設計は状態オブザーバゲインHの設計とは独立して行ってよいことを意味しており、大変重要な性質（分離定理）として知られています。

可制御正準形

投稿日時: 投稿者:

可制御正準形…Homework

[0]　n次系 \dot{x}=Ax+Bu\ (x\in{\rm\bf R}^n,u\in{\rm\bf R}^m)に対して、次を仮定します。

\displaystyle{(1)\quad {\rm rank}\left[\begin{array}{cccc} B & AB &\cdots & A^{n-1}B \end{array}\right]=n }

このとき、以下では、２つの座標変換を続けて行います。

\displaystyle{(2)\quad x'=T_1x\ \Rightarrow\ \dot{x}'=A_1x'+B_1u }

\displaystyle{(3)\quad x''=T_2x'\ \Rightarrow\ \dot{x}''=A_2x''+B_2u }

そして、A_2-B_2F_2の特性多項式を任意に設定するためのF_2の１つを示します。

[1]　１入力系の場合を考えます。Aの特性多項式を

\displaystyle{(4)\quad \det(\lambda I_n-A)=\lambda^n+a_1\lambda^{n-1}+\cdots+a_n }

としますと、ケーリ―・ハミルトンの定理より

\displaystyle{(5)\quad A^n+a_1A^{n-1}+\cdots+a_nI_n=0 }

が成り立ちます。これから

\displaystyle{(6)\quad A^nB=-a_nB-\cdots-a_1A^{n-1}B }

を得ます。これに基づいて

\displaystyle{(7)\quad \begin{array}{l} A \underbrace{ \left[\begin{array}{cccc} B & AB &\cdots & A^{n-1}B \end{array}\right] }_{T_1^{-1}}\\ = \underbrace{ \left[\begin{array}{cccc} B & AB &\cdots & A^{n-1}B \end{array}\right] }_{T_1^{-1}} \underbrace{ \left[\begin{array}{ccccc} 0 & \cdots & 0 & -a_n \\ 1 & \cdots & 0 & -a_{n-1} \\ \vdots & \ddots & \vdots &\vdots\\ 0 & \cdots & 1 & -a_1 \end{array}\right] }_{A_1} \end{array} }

を得ます。これから次の第１番目の正準形が定義されます。

\displaystyle{(8)\quad \boxed{ \begin{array}{l} A_1=T_1AT_1^{-1}= \left[\begin{array}{ccccc} 0 & \cdots & 0 & -a_n \\ 1 & \cdots & 0 & -a_{n-1} \\ \vdots & \ddots & \vdots &\vdots\\ 0 & \cdots & 1 & -a_1 \end{array}\right]\\ B_1=T_1B= \left[\begin{array}{cc} 1 \\ 0\\ \vdots\\ 0 \end{array}\right] \end{array}} }

さて、次式が成り立ちます。

\displaystyle{(9)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{ccccc} 0 & \cdots & 0 & -a_n \\ 1 & \cdots & 0 & -a_{n-1} \\ \vdots & \ddots & \vdots &\vdots\\ 0 & \cdots & 1 & -a_1 \end{array}\right] }_{A_1} \underbrace{ \left[\begin{array}{ccccc} a_{n-1} & a_{n-2} & \cdots & a_1 & 1 \\ a_{n-2} & a_{n-3} & \cdots & 1 & 0 \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ a_1 & 1 & \cdots & 0 & 0 \\ 1 & 0 & \cdots & 0 & 0 \end{array}\right] }_{T_2^{-1}}\\ =\underbrace{ \left[\begin{array}{ccccc} a_{n-1} & a_{n-2} & \cdots & a_1 & 1 \\ a_{n-2} & a_{n-3} & \cdots & 1 & 0 \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ a_1 & 1 & \cdots & 0 & 0 \\ 1 & 0 & \cdots & 0 & 0 \end{array}\right] }_{T_2^{-1}} \underbrace{ \left[\begin{array}{ccccc} 0 & 1 & \cdots & 0\\ \vdots & \vdots & \cdots & \vdots \\ 0 & 0 & \cdots & 1\\ -a_n & -a_{n-1} & \cdots & -a_1 \end{array}\right] }_{A_2} \end{array} }

実際、n=2のときは、次のように示されます。

\displaystyle{(10)\quad {\rm LHS}=\left[\begin{array}{cc} 0 & -a_2 \\ 1 & -a_1 \end{array}\right] \left[\begin{array}{cc} a_1 & 1 \\ 1 & 0 \end{array}\right] = \left[\begin{array}{cc} -a_2 & 0 \\ 0 & 1 \end{array}\right] }

\displaystyle{(11)\quad {\rm RHS}= \left[\begin{array}{cc} a_1 & 1 \\ 1 & 0 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ -a_2 & -a_1 \end{array}\right] = \left[\begin{array}{cc} -a_2 & 0 \\ 0 & 1 \end{array}\right] }

一般には、次のように示されます。

\displaystyle{(12)\quad \begin{array}{l} {\rm LHS}=\\ \left[\begin{array}{cc} 0_{1\times n-1} & -a_n \\ I_{n-1\times n-1} & \left[\begin{array}{c} -a_{n-1} \\ \vdots\\ -a_1 \end{array}\right] \end{array}\right] \left[\begin{array}{cccc|c} a_{n-1} & a_{n-2} & \cdots & a_1 & 1 \\ a_{n-2} & a_{n-3} & \cdots & 1 & 0 \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ a_1 & 1 & \cdots & 0 & 0 \\\hline 1 & 0 & \cdots & 0 & 0 \end{array}\right]\\ = \left[\begin{array}{cccc|c} -a_n & 0 & \cdots & 0 & 0 \\\hline 0 & a_{n-2} & \cdots & a_1 & 1 \\ 0 & a_{n-3} & \cdots & 1 & 0 \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ 0 & 1 & \cdots & 0 & 0 \end{array}\right] \end{array} }

\displaystyle{(13)\quad \begin{array}{l} {\rm RHS}=\\ \left[\begin{array}{cccc|c} a_{n-1} & a_{n-2} & \cdots & a_1 & 1 \\ a_{n-2} & a_{n-3} & \cdots & 1 & 0 \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ a_1 & 1 & \cdots & 0 & 0 \\\hline 1 & 0 & \cdots & 0 & 0 \end{array}\right] \left[\begin{array}{cc} 0_{n-1\times 1} & I_{n-1\times n-1} \\ -a_n & \left[\begin{array}{ccc} -a_{n-1} & \cdots & -a_1 \end{array}\right] \end{array}\right]\\ = \left[\begin{array}{c|cccc} -a_n & 0 & \cdots & 0 & 0 \\ 0 & a_{n-2} & \cdots & a_1 & 1 \\ 0 & a_{n-3} & \cdots & 1 & 0 \\ \vdots & \vdots & \cdots & \vdots & \vdots \\\hline 0 & 1 & \cdots & 0 & 0 \end{array}\right] \end{array} }

上の関係式(9)に基づいて、次の第２番目の正準形が定義されます。

\displaystyle{(14)\quad \boxed{ \begin{array}{l} A_2=T_2A_1T_2^{-1}= \left[\begin{array}{ccccc} 0 & 1 & \cdots & 0\\ \vdots & \vdots & \cdots & \vdots \\ 0 & 0 & \cdots & 1\\ -a_n & -a_{n-1} & \cdots & -a_1 \end{array}\right]\\ B_2=T_2B_1= \left[\begin{array}{cc} 0 \\ \vdots\\ 0\\ 1 \end{array}\right] \end{array}} }

したがって、A_2-B_2F_2の特性多項式を

\displaystyle{(15)\quad \det(\lambda I_n-A_2+B_2F_2)=\lambda^n+a_1'\lambda^{n-1}+\cdots+a_n' }

とするF_2は、次式で与えられます。

\displaystyle{(16)\quad F_2=\left[\begin{array}{ccccc} a_n'-a_n & a_{n-1}'-a_{n-1} & \cdots & a_1'-a_1 \end{array}\right] }

[2]　多入力系の場合を数値例を用いて説明します。

\displaystyle{(17)\quad \begin{array}{l} A= \left[\begin{array}{ccc} 0 & 7 & 4 \\ 1 & -1 & -2\\ 0 & 3 & 1 \end{array}\right],\ B= \left[\begin{array}{cc} 1 & 0 \\ 0 & 0 \\ 0 & 1 \end{array}\right] \end{array} }

に対して、可制御性は

\displaystyle{(18)\quad \begin{array}{l} {\rm rank} \left[\begin{array}{ccc} B & AB & A^2B \end{array}\right]= {\rm rank}\ \left[\begin{array}{cc|cc|cc} 1 & 0 & 0 & 4 & 7 & -10\\ 0 & 0 & 1 & -2 & -1 & 4\\ 0 & 1 & 0 & 1 & 3 & 5 \end{array}\right]=3 \end{array} }

となって成立しています。そこでBの第１列ベクトルと第２ベクトルをそれぞれB_1とB_2で表すとき、

\displaystyle{(19)\quad \begin{array}{l} B_1,B_2,AB_1,AB_2,A^2B_1,A^2B_2 \end{array} }

の順にベクトルの１次独立性を調べて、それらを取り出すと

\displaystyle{(20)\quad \begin{array}{l} B_1,B_2,AB_1 \end{array} }

を得ます。これを

\displaystyle{(21)\quad \begin{array}{l} B_1,AB_1,B_2 \end{array} }

のように並べ替えておきます。i=1,2に対して、A^\nu B_iが１次従属となる最小の\nuを可制御性指数と呼びます。まずi=1については、次式よりn_1=2となります。

\displaystyle{(22)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} 7 \\ -1 \\ 3 \end{array}\right] }_{A^2B_1} = \underbrace{7}_{\alpha_{110}} \underbrace{ \left[\begin{array}{c} 1 \\ 0 \\ 0 \end{array}\right] }_{B_1} +\underbrace{(-1)}_{\alpha_{111}} \underbrace{ \left[\begin{array}{c} 0 \\ 1 \\ 0 \end{array}\right] }_{AB_1} +\underbrace{3}_{\alpha_{120}} \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ 1 \end{array}\right] }_{B_2} \end{array} }

次にi=2については、次式よりn_2=1となります。

\displaystyle{(23)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c} 4 \\ -2 \\ 1 \end{array}\right] }_{AB_2} = \underbrace{4}_{\alpha_{210}} \underbrace{ \left[\begin{array}{c} 1 \\ 0 \\ 0 \end{array}\right] }_{B_1} +\underbrace{(-2)}_{\beta_{21}} \underbrace{ \left[\begin{array}{c} 0 \\ 1 \\ 0 \end{array}\right] }_{AB_1} +\underbrace{1}_{\alpha_{220}} \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ 1 \end{array}\right] }_{B_2} \end{array} }

これらに基づいて、次式を得ます。

\displaystyle{(24)\quad \begin{array}{l} A \underbrace{ \left[\begin{array}{ccc} B_1 & AB_1 & B_2 \end{array}\right] }_{T_1^{-1}} =\underbrace{ \left[\begin{array}{ccc} B_1 & AB_1 & B_2 \end{array}\right] }_{T_1^{-1}} \underbrace{ \left[\begin{array}{ccc} 0 & \alpha_{110} & \alpha_{210} \\ 1 & \alpha_{111} & \beta_{21}\\ 0 & \alpha_{120} & \alpha_{220} \end{array}\right] }_{A_1}\\ &&B= \underbrace{ \left[\begin{array}{ccc} B_1 & AB_1 & B_2 \end{array}\right] }_{T_1^{-1}} \underbrace{ \left[\begin{array}{cc} 1 & 0 \\ 0 & 0 \\ 0 & 1 \end{array}\right] }_{B_1} \end{array} }

これから次の第１番目の正準形が定義されます。

\displaystyle{(25)\quad \begin{array}{l} A_1=T_1AT_1^{-1}= \left[\begin{array}{ccc} 0 & \alpha_{110} & \alpha_{210} \\ 1 & \alpha_{111} & \beta_{21}\\ 0 & \alpha_{120} & \alpha_{220} \end{array}\right] =\left[\begin{array}{ccc} 0 & 7 & 4 \\ 1 & -1 & -2\\ 0 & 3 & 1 \end{array}\right]\\ & B_1=T_1B= \left[\begin{array}{cc} 1 & 0 \\ 0 & 0 \\ 0 & 1 \end{array}\right] \end{array} }

さて、次式が成り立ちます。

\displaystyle{(26)\quad \begin{array}{l} &&\underbrace{ \left[\begin{array}{cc|c} 0 & \alpha_{110} & \alpha_{210} \\ 1 & \alpha_{111} & \beta_{21}\\\hline 0 & \alpha_{120} & \alpha_{220} \end{array}\right] }_{A_1} \underbrace{ \left[\begin{array}{cc|c} -\alpha_{111} & 1 & -\beta_{21} \\ 1 & 0 & 0\\\hline 0 & 0 & 1 \end{array}\right] }_{T_2^{-1}}\\ &=\underbrace{ \left[\begin{array}{cc|c} -\alpha_{111} & 1 & -\beta_{21} \\ 1 & 0 & 0\\\hline 0 & 0 & 1 \end{array}\right] }_{T_2^{-1}} \underbrace{ \left[\begin{array}{cc|c} 0 & 1 & 0 \\ \alpha_{110}+\beta_{21}\alpha_{120} & \alpha_{111} & \alpha_{210}+\beta_{21}\alpha_{220}\\\hline \alpha_{120} & 0 & \alpha_{220} \end{array}\right] }_{A_2} \end{array} }

すなわち

\displaystyle{(27)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cc|c} 0 & 7 & 4 \\ 1 & -1 & -2\\\hline 0 & 3 & 1 \end{array}\right] }_{A_1} \underbrace{ \left[\begin{array}{cc|c} 1 & 1 & 2 \\ 1 & 0 & 0\\\hline 0 & 0 & 1 \end{array}\right] }_{T_2^{-1}} =\underbrace{ \left[\begin{array}{cc|c} 1 & 1 & 2 \\ 1 & 0 & 0\\\hline 0 & 0 & 1 \end{array}\right] }_{T_2^{-1}} \underbrace{ \left[\begin{array}{cc|c} 0 & 1 & 0 \\ 1 & -1 & 2\\\hline 3 & 0 & 1 \end{array}\right] }_{A_2} \end{array} }

実際、次のように示されます。

\displaystyle{(28)\quad \begin{array}{l} {\rm LHS}= \left[\begin{array}{ccc} 0 & \alpha_{110} & \alpha_{210} \\ 1 & \alpha_{111} & \beta_{21}\\ 0 & \alpha_{120} & \alpha_{220} \end{array}\right] \left[\begin{array}{ccc} -\alpha_{111} & 1 & -\beta_{21} \\ 1 & 0 & 0\\ 0 & 0 & 1 \end{array}\right]\\ &= \left[\begin{array}{ccc} \alpha_{110} & 0 & \alpha_{210} \\ 0 & 1 & 0 \\ \alpha_{120} & 0 & \alpha_{220} \end{array}\right] \end{array} }

\displaystyle{(29)\quad \begin{array}{l} {\rm RHS}=\\ & \left[\begin{array}{ccc} -\alpha_{111} & 1 & -\beta_{21} \\ 1 & 0 & 0\\ 0 & 0 & 1 \end{array}\right] \left[\begin{array}{ccc} 0 & 1 & 0 \\ \alpha_{110}+\beta_{21}\alpha_{120} & \alpha_{111} & \alpha_{210}+\beta_{21}\alpha_{220}\\ \alpha_{120} & 0 & \alpha_{220} \end{array}\right]\\ &= \left[\begin{array}{ccc} \alpha_{110} & 0 & \alpha_{210} \\ 0 & 1 & 0 \\ \alpha_{120} & 0 & \alpha_{220} \end{array}\right] \end{array} }

上の関係式(??)に基づいて、次の第２番目の正準形が定義されます。

\displaystyle{(30)\quad \begin{array}{l} &A_2=T_2A_1T_2^{-1}= \left[\begin{array}{ccc} 0 & 1 & 0 \\ 1 & -1 & 2\\ 3 & 0 & 1 \end{array}\right]\\ & B_2=T_2B_1= \left[\begin{array}{cc} 0 & 0 \\ 1 & -2 \\ 0 & 1 \end{array}\right] \end{array} }

これは次のような表現が可能です。

\displaystyle{(31)\quad \begin{array}{l} & A_2=A_0+B_0G_0^{-1}F_0\\ & B_2=B_0G_0^{-1} \end{array} }

ただし

\displaystyle{(32)\quad \begin{array}{l} &A_0= \left[\begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 0\\ 0 & 0 & 0 \end{array}\right],\ B_0= \left[\begin{array}{cc} 0 & 0 \\ 1 & 0 \\ 0 & 1 \end{array}\right]\\ &F_0= \left[\begin{array}{cc|c} \alpha_{110} & \alpha_{111} & \alpha_{210}\\\hline \alpha_{120} & 0 & \alpha_{220} \end{array}\right] = \left[\begin{array}{cc|c} 7 & -1 & 4\\\hline 3 & 0 & 1 \end{array}\right]\\ &&G_0= \left[\begin{array}{cc} 1 & -\beta_{21} \\ 0 & 1 \end{array}\right] = \left[\begin{array}{cc} 1 & 2 \\ 0 & 1 \end{array}\right] \end{array} }

したがって、A_2-B_2F_2の特性多項式を

\displaystyle{(33)\quad \begin{array}{l} \det(\lambda I_3-A_2+B_2F_2)=(\lambda^2+a_1'\lambda+a_2')(\lambda+a_3') \end{array} }

とするF_2の一つは、次式で与えられます。

\displaystyle{(34)\quad \begin{array}{l} F_2=F_0+ \left[\begin{array}{ccccc} a_2' & a_1' & 0\\ 0 & 0 & a_3' \end{array}\right] \end{array} }

実際

\displaystyle{(35)\quad \begin{array}{l} A_2-B_2F_2= \left[\begin{array}{ccc} 0 & 1 & 0 \\ -a_2' & -a_1' & 0\\ 0 & 0 & -a_3' \end{array}\right] \end{array} }

[3]　多入力系の場合を考えます。Bのi番目の列ベクトルをB_iと表すとき

\displaystyle{(36)\quad \begin{array}{l} B_1,\cdots,B_m,AB_1,\cdots,AB_m,A^2B_1,\cdots,A^2B_m,\cdots \end{array} }

の順にベクトルの１次独立性を調べて、A^\nu B_iが１次従属となる最小の\nuとして可制御性指数

\displaystyle{(37)\quad \begin{array}{l} n_1,\ n_2,\ \cdots,\ n_m \end{array} }

を定めます。得られた１次独立なベクトルを

\displaystyle{(38)\quad \begin{array}{l} B_1,\cdots,A^{n_1-1}B_1,\cdots,B_m,\cdots,A^{n_m-1}B_m \end{array} }

のように並べ替えておきます。このとき

\displaystyle{(39)\quad \begin{array}{l} A^{n_i}B_i=\sum_{j=1}^{m} \left[\begin{array}{ccc} B_j\ \cdots\ A^{n_j-1}B_j \end{array}\right]\alpha_{ij} \end{array} }

ただし

\displaystyle{(40)\quad \begin{array}{l} \alpha_{ij}= \left\{\begin{array}{l} \left[\begin{array}{c} \alpha_{ij0} \\ \vdots \\ \alpha_{ij,n_j-1} \end{array}\right]\ (n_j\le n_i)\\ \left[\begin{array}{c} \alpha_{ij0} \\ \vdots \\ \alpha_{ij,n_i-1}\\ \beta_{ij} \\ 0_{n_j-n_i-1\times 1} \\ \end{array}\right]\ (j<i, n_j>n_i)\\ \left[\begin{array}{c} \alpha_{ij0} \\ \vdots \\ \alpha_{ij,n_i-1}\\ 0 \\ 0_{n_j-n_i-1\times 1} \\ \end{array}\right]\ (j>i, n_j>n_i) \end{array}\right. \end{array} }

が成り立ちます。これから、次式を得ます。

\displaystyle{(41)\quad \begin{array}{l} A \underbrace{ \left[\begin{array}{c|c|c} B_1\ \cdots\ A^{n_1-1}B_1 & \cdots & B_m\ \cdots\ A^{n_m-1}B_m \end{array}\right] }_{T_1^{-1}}\\ =\underbrace{ \left[\begin{array}{c|c|c} B_1\ \cdots\ A^{n_1-1}B_1 & \cdots & B_m\ \cdots\ A^{n_m-1}B_m \end{array}\right] }_{T_1^{-1}}\\ \times \underbrace{ \left[\begin{array}{c|c|c} \begin{array}{c|c} \begin{array}{c} 0_{1\times n_1-1} \\ I_{n_1-1} \end{array} & \alpha_{11} \end{array} &\cdots & \begin{array}{c|c} 0_{n_1\times n_m-1} & \alpha_{m1} \end{array} \\\hline \vdots & \cdots & \vdots \\\hline \begin{array}{c|c} 0_{n_m\times n_1-1} & \alpha_{1m} \end{array} & \cdots & \begin{array}{c|c} \begin{array}{c} 0_{1\times n_m-1} \\ I_{n_m-1} \end{array} & \alpha_{mm} \end{array} \end{array}\right] }_{A_1}\\ B= \underbrace{ \left[\begin{array}{c|c|c} B_1\ \cdots\ A^{n_1-1}B_1 & \cdots & B_m\ \cdots\ A^{n_m-1}B_m \end{array}\right] }_{T_1^{-1}}\\ \times \underbrace{ \left[\begin{array}{c|c|c} \begin{array}{c} 1 \\ 0_{n_1-1\times1} \end{array} &\cdots & 0_{n_1\times1} \\\hline \vdots & \cdots & \vdots \\\hline 0_{n_m\times1} & \cdots & \begin{array}{c} 1 \\ 0_{n_m-1\times1} \end{array} \end{array}\right] }_{B_1} \end{array} }

これに基づいて、次の第１番目の正準形が定義されます。

\displaystyle{(42)\quad \boxed{ \begin{array}{l} A_1=T_1AT_1^{-1}= \left[\begin{array}{c|c|c} \begin{array}{c|c} \begin{array}{c} 0_{1\times n_1-1} \\ I_{n_1-1} \end{array} & \alpha_{11} \end{array} &\cdots & \begin{array}{c|c} 0_{n_1\times n_m-1} & \alpha_{m1} \end{array} \\\hline \vdots & \cdots & \vdots \\\hline \begin{array}{c|c} 0_{n_m\times n_1-1} & \alpha_{1m} \end{array} & \cdots & \begin{array}{c|c} \begin{array}{c} 0_{1\times n_m-1} \\ I_{n_m-1} \end{array} & \alpha_{mm} \end{array} \end{array}\right]\\ B_1=T_1B= \left[\begin{array}{c|c|c} \begin{array}{c} 1 \\ 0_{n_1-1\times1} \end{array} &\cdots & 0_{n_1\times1} \\\hline \vdots & \ddots & \vdots \\\hline 0_{n_m\times1} & \cdots & \begin{array}{c} 1 \\ 0_{n_m-1\times1} \end{array} \end{array}\right] \end{array}} }

さて、次式が成り立ちます。(この証明はまだ完成していません)

\displaystyle{(43)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{c|c|c} \begin{array}{c|c} \begin{array}{c} 0_{1\times n_1-1} \\ I_{n_1-1} \end{array} & \alpha_{11} \end{array} &\cdots & \begin{array}{c|c} 0_{n_1\times n_m-1} & \alpha_{m1} \end{array} \\\hline \vdots & \cdots & \vdots \\\hline \begin{array}{c|c} 0_{n_m\times n_1-1} & \alpha_{1m} \end{array} & \cdots & \begin{array}{c|c} \begin{array}{c} 0_{1\times n_m-1} \\ I_{n_m-1} \end{array} & \alpha_{mm} \end{array} \end{array}\right] }_{A_1}\\ \times \underbrace{ \left[\begin{array}{ccc} L_{11} & \cdots & L_{m1} \\ \vdots & \cdots & \vdots \\ L_{1m} & \cdots & L_{mm} \end{array}\right] }_{T_2^{-1}} =\underbrace{ \left[\begin{array}{ccc} L_{11} & \cdots & L_{m1} \\ \vdots & \cdots & \vdots \\ L_{1m} & \cdots & L_{mm} \end{array}\right] }_{T_2^{-1}} \underbrace{ (A_0+B_0G_0^{-1}F_0) }_{A_2} \end{array} }

ただし

\displaystyle{(44)\quad \begin{array}{l} L_{ii}=-\left[\begin{array}{ccc} S_{n_i}\alpha_{ii} & \cdots & S_{n_i}^{n_i}\alpha_{ii} \end{array}\right]+J_{n_1}\\ L_{ij}=-\left[\begin{array}{ccc} S_{n_j}\alpha_{ij} & \cdots & S_{n_j}^{n_i}\alpha_{ij} \end{array}\right]\\ S_\nu= \left[\begin{array}{cc} 0_{\nu-1\times 1} & I_{\nu-1}\\ 0 & 0_{1\times\nu-1} \end{array}\right] \\ J_\nu= \left[\begin{array}{ccc} 0 & \cdots & 1 \\ \vdots & \cdot & \vdots \\ 1 & \cdots & 0 \end{array}\right] \end{array} }

および

\displaystyle{(45)\quad \begin{array}{l} A_0= \left[\begin{array}{c|c|c} S_{n_1} & \cdots & 0_{n_1\times n_m} \\\hline \vdots & \ddots & \vdots \\\hline 0_{n_m\times n_1} & \cdots & S_{n_m} \end{array}\right]\\ B_0= \left[\begin{array}{c|c|c} \begin{array}{c} 0_{n_1-1\times 1} \\ 1 \end{array} & \cdots & 0_{n_1\times1} \\\hline \vdots & \ddots & \vdots \\\hline 0_{n_m\times1} & \cdots & \begin{array}{c} 0_{n_m-1\times 1} \\ 1 \end{array} \end{array}\right]\\ F_0= \left[\begin{array}{ccc|c|ccc} \bar\alpha_{110} & \cdots & \bar\alpha_{11,n_1-1} & \cdots & \bar\alpha_{m10} & \cdots & \bar\alpha_{m1,n_m-1}\\ \vdots & \cdots & \vdots & \cdots & \vdots & \cdots & \vdots \\ \bar\alpha_{1m0} & \cdots & \bar\alpha_{1m,n_1-1} & \cdots & \bar\alpha_{mm0} & \cdots & \bar\alpha_{mm,n_m-1}\\ \end{array}\right]\\ \bar\alpha_{ijk}= \left\{\begin{array}{ll} \alpha_{ijk} & (k\verb|<| n_j)\\ 0 & (k\ge n_j)\end{array}\right.\\ G_0=\left[\begin{array}{cccl} 1&\bar\beta_{21}&\cdots&\bar\beta_{m1}\\ 0&1&\ddots&\vdots\\ \vdots&\ddots&\ddots&\bar\beta_{m,m-1}\\ 0&0&\cdots&1\end{array}\right]\\ \bar\beta_{ij}=\left\{\begin{array}{ll} -\beta_{ij} &(n_j\ge n_i)\\ 0 &(n_j\verb|<|n_i) \end{array}\right. \end{array} }

次の第２番目の正準形が定義されます。

\displaystyle{(46)\quad \boxed{\begin{array}{l} A_2=A_0+B_0G_0^{-1}F_0\\ B_2=B_0G_0^{-1} \end{array}} }

可制御標準形

投稿日時: 投稿者:

可制御標準形…Homework

[1] ある平衡状態周りの挙動が線形状態方程式

\displaystyle{(1)\quad \dot{x}(t)=Ax(t)+Bu(t)\quad(x(t)\in{\rm\bf R}^n,u(t)\in{\rm\bf R}^m) }

で表される制御対象に対して、座標変換

\displaystyle{(2)\quad x'(t)=Tx(t)\Leftrightarrow x(t)=T^{-1}x'(t) }

を行うと（Tは正則行列）、状態方程式は次式となります。

\displaystyle{(3)\quad \dot{x}'(t)=\underbrace{TAT^{-1}}_{A'}x'(t)+\underbrace{TB}_{B'}u(t) }

このとき、漸近安定性、可安定性、可制御性は、座標変換によって不変であることを確かめることができます。実際、漸近安定性については

\displaystyle{(4)\quad {\rm det}(\lambda I_n-\underbrace{TAT^{-1}}_{A'}) ={\rm det}T(\lambda I_n-A)T^{-1} ={\rm det}(\lambda I_n-A) }

から明らかです。また、(1)に対する状態フィードバック

\displaystyle{(5)\quad u(t)=-Fx(t) }

による閉ループ系

\displaystyle{(6)\quad \dot{x}(t)=(A-BF)x(t) }

は座標変換後は

\displaystyle{(7)\quad u(t)=-\underbrace{FT^{-1}}_{F'}x'(t) }

を(3)に代入して

\displaystyle{(8)\quad \dot{x}'(t)=(TAT^{-1}-TBFT^{-1})x'(t) }

となることに注意します。このとき、可安定性については

\displaystyle{(9)\quad \begin{array}{l} {\rm det}(\lambda I_n-\underbrace{TAT^{-1}}_{A'}+\underbrace{TB}_{B'}\underbrace{FT^{-1}}_{F'}) ={\rm det}T(\lambda I_n-A+BF)T^{-1}\\ ={\rm det}(\lambda I_n-A+BF) \end{array} }

から明らかです。また、可制御性については

\displaystyle{(10)\quad \begin{array}{ll} {\rm rank}\, [\begin{array}{cccc} \underbrace{TB}_{B'} & \underbrace{TAT^{-1}}_{A'}\underbrace{TB}_{B'} & \cdots & \underbrace{(TAT^{-1})^{n-1}}_{A'^{n-1}}\underbrace{TB}_{B'} \end{array}] \\ ={\rm rank}\, T[\begin{array}{cccc} B & AB &\cdots & A^{n-1}B \end{array}]\\ ={\rm rank}\, [\begin{array}{cccc} B & AB &\cdots & A^{n-1}B \end{array}] \end{array} }

から明らかです。

●可安定性と可制御性は座標変換によって不変であることに基づいて，可安定性や可制御性の判別をしやすいように座標変換行列Tを選ぶことを考えます。すなわち，n次系(1)に対して，次の可制御標準形に座標変換できるかどうかを調べます。

\displaystyle{(11)\quad \boxed{ \left[\begin{array}{c} \dot{x}_1'(t) \\ \dot{x}_2'(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} A_{\#} & X \\ 0 & A_k \end{array}\right] }_{TAT^{-1}} \left[\begin{array}{c} x_1'(t) \\ x_2'(t) \end{array}\right] + \underbrace{ \left[\begin{array}{c} B_{\#} \\ 0 \end{array}\right] }_{TB} u(t)} }

ただし，(A_{\#},B_{\#})は可制御対とします。このブロック線図を描いてみればわかるように，状態x_2'は入力uの影響を受けることはありません。したがって，可制御性は成立しないこと，そしてA_kが安定行列かどうかが可安定性と関係していることは想像がつくと思います。以下では，次の判定法が導かれること示します。

「可制御標準形に座標変換したとき，A_kが安定行列のとき可安定性が成り立ち，そうでないときは可安定性は成り立たない。また，可制御標準形においてA_kが存在しなければ，可制御性が成り立っている。」

[2]　(11)を得るような座標変換行列Tの候補は，階段化アルゴリズム（Note A33参照）を用いて、次のように得られます。

定理　サイズn\times nの行列Aとサイズn\times mの行列Bに対して，つぎのように変換する直交行列Tが存在する。

\displaystyle{(12)\quad \left[\begin{array}{cc} T^TB & T^TAT \end{array}\right]= \left[\begin{array}{c|cccc|c} B_1 & A_1 & X & \cdots & X & X \\ 0 & B_2 & A_2 & \cdots & X & X \\ \vdots & & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & B_{k-1} & A_{k-1} & X \\ \hline 0 & 0 & \cdots & 0 & B_k & A_k \end{array}\right] }

ここで，kは正整数，またB_jのサイズをm_j\times m_{j-1}とするとき(m_0=m,\ m_1+\cdots+m_k=n)，つぎが成り立つ（B_1,\cdots,B_{k-1}は横長の形状となり，行フルランクをもつ。また，B_kは横長で，行フルランクをもつ行列または零行列のどちらかである。）。

\displaystyle{(13)\quad {\rm rank}\,B_j=m_j\quad (j=1,\cdots,k-1) }
\displaystyle{(14)\quad {\rm rank}\,B_k=m_k\ または\ 0 }
\displaystyle{(15)\quad m\ge m_1 \ge \cdots \ge m_{k-1} \ge m_k }

このとき，(A,B)が可制御対であるための必要十分条件は

条件C0:　(12)において，B_k\ne 0

が成り立つことである。また，(A,B)が可安定対であるための必要十分条件は、つぎが成り立つことである。

条件S0:　(12)において，B_k=0のとき，A_kは安定行列

B_k=0のとき，(11)のような形式が得られていることは，(12)のTを上のT^{-1}=T^Tと読み替えれば明らかであろう。

●可制御性の定義とその等価な条件は次のようにまとめられます。

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

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

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

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

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

条件C0:　(12)において，B_k\ne 0

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

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

すでに定義DC\Leftrightarrow条件C1\Leftrightarrow条件C2を示していますので、以下では条件C3\Rightarrow条件C0\Rightarrow条件C4\Leftrightarrow条件C5\Rightarrow条件C2を示します。

＜条件C3\Rightarrow条件C0＞ B_k=0とすると，(12)を用いて

\displaystyle{(16)\quad {\rm det}(\lambda I_n-T^TAT+T^TB\cdot FT) ={\rm det}\left(\lambda I_n- \left[\begin{array}{cc} X & X \\ 0 & A_k \end{array}\right] \right) }

となり，閉ループ系においてA_kの固有値はそのまま残るが，C3に矛盾しています。

＜条件C0\Rightarrow条件C4＞　\lambdaを任意の複素数とする。[\begin{array}{cc} T^TB & T^TAT-\lambda I_n \end{array}]に，(12)を適用すると

\displaystyle{(17)\quad \left[\begin{array}{c|cccc|c} B_1 & A_1-\lambda I_{m_1} & X & \cdots & X & X \\ 0 & B_2 & A_2-\lambda I_{m_2} & \cdots & X & X \\ \vdots & & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & B_{k-1} & A_{k-1}-\lambda I_{m_{k-1}} & X \\ \hline 0 & 0 & \cdots & 0 & B_k & A_k-\lambda I_{m_k} \end{array}\right] }

を得ます。ここで，B_1,\cdots,B_kは横長の形状となり，行フルランク(行数に等しい階数)をもつので

\displaystyle{(18)\quad \left[\begin{array}{cc} T^TB & T^TAT-\lambda I_n \end{array}\right] = T^T \left[\begin{array}{cc} B & A-\lambda I_n \end{array}\right] \left[\begin{array}{cc} I_n & 0 \\ 0 & T \end{array}\right] }

に，シルベスターの公式を適用して

\displaystyle{(19)\quad {\rm rank}\, [\begin{array}{cc} T^TB & T^TAT-\lambda I_n \end{array}] = {\rm rank}\, [\begin{array}{cc} B & A-\lambda I_n \end{array}] =n }

を得ます。この階数がnより小さくなるのは，(17)から明らかなように，\lambdaがA_kの固有値に等しい場合です。したがって，\lambdaはAの固有値について調べれば十分である。すなわち，条件C4を得たことになります。

＜条件C4\Leftrightarrow条件C5＞　条件C4は，wをn次元ベクトルとして

\displaystyle{(20)\quad w^T\, [\begin{array}{cc} B & A-\lambda I_n \end{array}] =0 \Rightarrow w =0 }

すなわち

\displaystyle{(21)\quad B^Tw=0,\ A^Tw=\lambda w \Rightarrow w=0 }

と等価です。

＜条件C5\Rightarrow条件C2＞　条件C5より

\displaystyle{(22)\quad w^TA^{i-1}B=0\quad (i=0,1,2,\cdots) }

を得て

\displaystyle{(23)\quad w^T\, [\begin{array}{ccccc} B & AB & \cdots & A^{n-1}B \end{array}] =0 \Rightarrow w =0 }

が成り立つ。これは条件C2を意味します。

[3]　可安定性の定義とその等価な条件は次のようにまとめられます。

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

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

条件S0:　(12)において，B_k=0のとき，A_kは安定行列

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

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

＜定義DS\Leftrightarrow条件S0＞　可制御であれば可安定であるので，不可制御で可安定の場合をどのように特徴付けるかを考えます。不可制御であれば条件C0が不成立なので，B_k=0となり，A_kの固有値だけは状態フィードバックによってどうすることもできません。このとき可安定であるためには，A_kが安定行列であることを前提とするしかありません。これは条件S0そのものです。また，可安定性と可制御性と相違は，前者がB_k\ne0ばかりでなくB_k=0を許すところにあることに注意します。

＜条件S0\Leftrightarrow条件S1\Leftrightarrow条件S2＞　条件S1\Leftrightarrow条件S2は，条件C4\Leftrightarrow条件C5の場合と同様にして明らかです。そこで、条件S0との等価性を調べます。B_k=0のとき，(17)の行列のランクが落ちるのは

\displaystyle{(24)\quad {\rm rank}\, \left. [\begin{array}{cc} B_k & A_k-\lambda I_{m_k} \end{array}] \right|_{B_k=0} = {\rm rank}\, [\begin{array}{cc} 0 & A_k-\lambda I_{m_k} \end{array}] \le m_k }

の部分で，\lambdaがA_kの固有値に一致した場合になります。したがって，条件S0を保証するためには，\lambdaにA_kのすべての不安定固有値を入れて，(17)を調べ，行フルランクとなればよいことになります。

Note A33　階段化アルゴリズム　

アルゴリズム＜階段化アルゴリズム(staircase algorithm)＞
入力パラメータ m,\,n,\,A(n\times n),\,B(n\times m)
出力パラメータ k,\,m_1,\cdots,m_k,\,T(n\times n)
ステップ1 初期化

j=0,\,s=0,\,m_0=m,\,T^{(1)}=I_n,\,B^{(1)}=B,\,A^{(1)}=Aとおく。

ステップ2 B^{(j)}の階数決定

jをj+1に更新し，B^{(j)}(n-s\times m_{j-1})に対して，つぎの特異値分解を行い，B^{(j)}の階数m_jとU^{(j)}(n-s\times n-s)を求める。

\displaystyle{(1)\quad B^{(j)}=U^{(j)}\Sigma^{(j)}V^{(j)}\,^T }

\displaystyle{(2)\quad \Sigma^{(j)}= \left[\begin{array}{cc} \Sigma_1^{(j)} & 0 \\ 0 & 0 \end{array}\right] }

もしm_j=n-sまたはm_j=0ならば，k=j,\,T=T^{(j)}と設定して終了する。そうでないときは，ステップ3へ行く。

ステップ3 座標変換

W^{(j)}(n\times n)を，U^{(j)}(n-s\times n-s)を用いて

\displaystyle{(3)\quad W^{(j)}= \left[\begin{array}{cc} I_s & 0 \\ 0 & U^{(j)} \end{array}\right] }

と設定し，A^{(j+1)}(n\times n)とT^{(j+1)}(n\times n)を，次式から求める。

\displaystyle{(4)\quad A^{(j+1)}=W^{(j)}\,^TA^{(j)}W^{(j)},\ T^{(j+1)}=T^{(j)}W^{(j)} }

ステップ4 B^{(j)}の切り出し

sをs+m_jに更新し，A^{(j+1)}(n\times n)に対してつぎの分割を行い，B^{(j+1)}(n-s\times m_j)を得て,ステップ2へ戻る。

\displaystyle{(5)\quad A^{(j+1)}= \left[\begin{array}{ccc} X(s\times s-m_j) & X(s\times m_j) & X(s\times n-s) \\ 0 & B^{(j+1)} & X \end{array}\right] }

●階段化アルゴリズムを適用した数値例を示します。

\displaystyle{(6)\quad A= \left[\begin{array}{ccc} 0 & 0 & 0 \\ 0 & -1 & 1 \\ 0 & 0 & 0 \end{array}\right] ,\ B= \left[\begin{array}{cc} 1 & 0 \\ 0 & 1 \\ 0 & 1 \end{array}\right] }

に対して、階段化アルゴリズムの各ステップは，つぎのように計算されます。

ステップ1：j=0,\,s=0,\,m_0=2,\,T^{(1)}=I_3,\,B^{(1)}=B,\,A^{(1)}=Aとおく。
ステップ2：j=1とし，B^{(1)}に対する特異値分解

\displaystyle{(7)\quad B^{(1)}= \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 1 \\ \frac{1}{\sqrt{2}} & 0 & -\frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} \end{array}\right] }_{U^{(1)}} \underbrace{ \left[\begin{array}{cc} \sqrt{2} & 0 \\ 0 & 1 \\ 0 & 0 \end{array}\right] }_{\Sigma^{(1)}} \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right]^T }_{V^{(1)T}} }

から，B^{(1)}の階数m_1=2とU^{(1)}が求まる。
ステップ3：m_1\ne3かつm_1\ne0なので，ステップ4へ行く。
ステップ4：W^{(1)}=U^{(1)}とし，A^{(2)}=W^{(1)}\,^TA^{(1)}W^{(1)}を求めると

\displaystyle{(8)\quad A^{(2)}= \left[\begin{array}{cc|c} 0 & 0 & 1 \\ 0 & 0 & 0 \\ \hline 0 & 0 & -1 \end{array}\right] }

となり，T^{(2)}=W^{(1)}とする。
ステップ5：s=2とする。A^{(2)}に対する分割から，B^{(2)}= [\begin{array}{cc} 0 & 0 \end{array}]を得る。
ステップ2′：j=2とする。B^{(2)}は零行列なので，その階数はm_2=0である。
ステップ3′：m_2=0なので，k=2,\,T=T^{(2)}のように設定する。

以上から，(12)は，つぎのように得られた。

\displaystyle{(9)\quad \begin{array}{l} \left[\begin{array}{cc} T^TB & T^TAT \end{array}\right]= \left[\begin{array}{ccc} B_1 & A_1 & X \\ 0 & B_2 & A_2 \end{array}\right]= \left[\begin{array}{cc|cc|c} 0 & \sqrt{2} & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & 0 & -1 \end{array}\right] \end{array} }

ここで，B_2は零行列で，A_2=-1<0より，条件S0が成り立つ。よって，(A,B)は可安定である。

演習 A33…Flipped Classroom

1^\circ 階段化アルゴリズムのコードは次のように書ける。上の数値例で確かめよ。

MATLAB 
%staircase.m
%----- 
 function [T,m]=staircase(A,B,tol)
   [n,r]=size(B);
   j=0; s=0; T=eye(n); B1=B; A1=A;
   while j＜n
     j=j+1; [U1,S1,V1]=svd(B1);
     m(j)=rank(B1,tol);
     if (m(j)==n-s)|(m(j)==0),k=j; break, end
     W=[eye(s)       zeros(s,n-s) ;
        zeros(n-s,s) U1          ];
     T=T*W; A1=W'*A1*W;
     s=s+m(j); B1=A1(s+1:n,s-m(j)+1:s);
 end
%-----
%eof

可安定性と可制御性

投稿日時: 投稿者:

可安定性と可制御性…Homework

[1]　ある平衡状態周りの挙動が線形状態方程式

\displaystyle{(1)\quad \dot{x}(t)=Ax(t)+Bu(t)\quad(x(t)\in{\rm\bf R}^n,u(t)\in{\rm\bf R}^m) }

で表される制御対象に対して、状態フィードバック

\displaystyle{(2)\quad u(t)=-Fx(t) }

による閉ループ系

\displaystyle{(3)\quad \dot{x}(t)=(A-BF)x(t) }

を安定化できる条件について調べます。

●状態フィーバックにより安定化できることを可安定性と言います。

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

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

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

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

証明はあとで述べますが、可安定性の判定は行列AとBを用いて行われるので、可安定性が成り立つとき、対(A,B)は可安定対という言い方をします。

●可安定性の十分条件として次の可制御性が知られています。

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

定義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}　(\lambdaはAのすべての固有値)

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

可制御性の判定も行列AとBを用いて行われるので、可制御性が成り立つとき、対(A,B)は可制御対という言い方をします。条件C1の積分式を可制御性グラミアン行列、条件C2の行列を可制御性行列と呼びます。また条件C4を満足するAの固有値を\lambdaは可制御固有値、満足しない固有値を不可制御固有値と呼びます。

[2]　以下では、可制御性の条件について、まず
　定義DC\Leftrightarrow条件C1\Leftrightarrow条件C2
を証明します。特に、命題P\Rightarrow Qを証明するのに、
　not(P\Rightarrow Q)\quad \Leftrightarrow\quad P\ and\ not(Q)
を用いて矛盾が出ることを示していることに注意してください。

●＜定義C0\Rightarrow条件C1＞　ある時刻tで可制御性グラミアン行列が正定でないとします。このとき，n次元ベクトルw\ne0が存在して

{\displaystyle{(4)\quad w^T\left(\int_0^t \exp(A\tau)BB^T\exp(A^T\tau)\,d\tau\right)w=0 }

が成り立ちます。これより

{\displaystyle{(5)\quad \int_0^t ||B^T\exp(A^T\tau)w||^2\,d\tau=0 \\ \Rightarrow w^T\exp(A\tau)B=0\ (0\le \tau\le t) }

を得ます。いま，初期状態x(0)をx(t)=0に移すことを考えると

{\displaystyle{(6)\quad 0=\exp(At)x(0)+\int_0^{t} \exp(A(t-\tau))Bu(\tau)\,d\tau }

と書けます。ここで，変数変換により

{\displaystyle{(7)\quad \int_0^{t} \exp(A(t-\tau))Bu(\tau)\,d\tau= \int_0^{t} \exp(A\tau)Bu(t-\tau)\,d\tau }

となることに注意して，(6)の左からw^Tをかけると

{\displaystyle{(8)\quad w^T\exp(At)x(0)=0 }

を得ます。ここで，x(0)=\exp(-At)wのときw=0となり，矛盾が生じます。よって条件C1が成り立ちます。

●＜定義DC\Leftarrow条件C1＞　(1)に対して，0\le \tau\le tにおいて，入力を

{\displaystyle{(9)\quad \begin{array}{l} u(\tau) =B^T\exp(A^T(t-\tau)) \\ \times \displaystyle{\left(\int_0^{t} \exp(A\tau)BB^T\exp(A^T\tau)\,d\tau\right)^{-1}} (x^*-\exp(At)x(0)) \end{array} }

と定めれば，初期状態x(0)が移される先の状態x(t)は、(9)を(1)の解

{\displaystyle{(10)\quad x(t)=\exp(At)x(0)+\int_0^{t} \exp(A(t-\tau))Bu(\tau)\,d\tau }

に代入して，次式に注意すれば，x^*と計算されます。

{\displaystyle{(11)\quad \begin{array}{l} \displaystyle{\int_0^{t} \exp(A(t-\tau))BB^T\exp(A^T(t-\tau))\,d\tau} \nonumber \\ =\displaystyle{\int_0^{t} \exp(A\tau)BB^T\exp(A^T\tau)\,d\tau} \end{array} }

以上で，t,x(0),x^*は任意であるので、定義DCを得ていることになります。

●＜条件C1\Rightarrow条件C2＞　可制御性行列は行フルランク(行数に等しい階数)とはならないとすします。このとき，あるn次元ベクトルw\ne0が存在して

{\displaystyle{(12)\quad w^T\, [\begin{array}{cccc} B & AB & \cdots & A^{n-1}B \end{array}]=0 }

が成り立ちます。ケーリー・ハミルトンの定理を用いて

{\displaystyle{(13)\quad w^TA^iB=0\quad (\forall i\ge 0) }

したがって

{\displaystyle{(14)\quad w^T\exp(A\tau)B=0 }

が成り立ちます。これより，ある時刻tに対して

{\displaystyle{(15)\quad w^T\left(\int_0^t \exp(A\tau)BB^T\exp(A^T\tau)\,d\tau \right)w=0 }

を得ますが，これは条件C1と矛盾です。よって条件C2が成り立ちます。

●＜条件C1\Leftarrow条件C2＞　ある時刻tで可制御性グラミアン行列が正定でないとします。このとき，w\ne0に対して(14)が成り立ちます。この第i回微分を求めて，t=0とおくと

{\displaystyle{(16)\quad w^TA^iB=0\quad (i=0,1,\cdots,n-1) }

これより，(12)を得ますが，これは条件C2と矛盾しています。よって条件C1が成り立ちます。

●あとで＜条件C3\Rightarrow条件C4\Leftrightarrow条件C5\Rightarrow条件C2＞を可制御標準形を求めて、＜条件C2\Rightarrow条件C3＞を可制御正準形を求めて示します。

●上の可制御性と可安定性のさまざまな条件のうち、理論展開では条件C5と条件S2が、数値計算では条件C4と条件S1がよく用いられます。特に、条件C5と条件S2による判定法は、PBH(Popov-Blevitch-Hautus)法と呼ばれています。

演習 A32-1…Flipped Classroom

次のコードは、PBH法の骨格となる部分を示している。これを実行して得られる変数contとstabの解釈を行え。

MATLAB 
%a32_1.m
%-----
 clear all, close all
 A=[0 0;0 -1]; B=[1;0]; 
 n=size(A,1);
 r=eig(A);
%----- 
 cont=[];
 for i=1:n
   w1=rank([B A-r(i)*eye(n)]);
   w2=min(svd([B A-r(i)*eye(n)])); 
   cont=[cont; i,r(i),w1,w2];
 end
 cont
%-----
stab=[];
 for i=1:n
   if real(r(i))>=0  
     w1=rank([B A-r(i)*eye(n)]);
     w2=min(svd([B A-r(i)*eye(n)])); 
     stab=[stab; i,r(i),w1,w2];
   end 
 end
 stab
%-----
%eof

演習 A32-2…Flipped Classroom

倒立振子CIP、CIP2、AIP、PIP、DIPの可制御性を調べよ。

MATLAB 
%lPIP.m
%-----
 clear all, close all
 global mc m1 m2 ell1 ell2 g th10 th20
 mc=1; m1=0.1; m2=0.2; ell1=0.15; ell2=0.2; g=9.8;
%----- 線形化
 A=[zeros(3,3) eye(3);zeros(3,6)];
 A(4,1)=0; 
 A(4,2)=-3*m1*g/(m1+m2+4*mc); 
 A(4,3)=-3*m2*g/(m1+m2+4*mc);  
 A(5,1)=0; 
 A(5,2)=(12*m1+3*m2+12*mc)*g/((4*m1+4*m2+16*mc)*ell1);  
 A(5,3)=9*m2*g/((4*m1+4*m2+16*mc)*ell1);   
 A(6,1)=0; 
 A(6,2)=9*m1*g/((4*m1+4*m2+16*mc)*ell2);  
 A(6,3)=(3*m1+12*m2+12*mc)*g/((4*m1+4*m2+16*mc)*ell2);   
 B=zeros(6,1);
 B(4)=4/(m1+m2+4*mc);
 B(5)=-3/((m1+m2+4*mc)*ell1);
 B(6)=-3/((m1+m2+4*mc)*ell2); 
%----- 可制御性
 n=size(A,1);
 r=eig(A);
 cont=[];
 for i=1:n
   w1=rank([B A-r(i)*eye(n)]);
   w2=min(svd([B A-r(i)*eye(n)])); 
   cont=[cont; i,r(i),w1,w2];
 end
 cont
%-----
%end

Note A32　可制御性は状態フィードバックにより不変

(A,B)が可制御対ならば(A-BF,B)も可制御対であることは次のように証明されます。

条件C5より，任意の\lambdaに対して

\displaystyle{(1)\quad B^Tw=0,\ A^Tw=\lambda w \Rightarrow\ w=0 }

成り立つので，これを仮定して

\displaystyle{(2)\quad B^Tw=0,\ (A-BF)^Tw=\lambda w \Rightarrow\ w=0 }

を示します。実際，B^Tw=0より

\displaystyle{(3)\quad (A-BF)^Tw=A^Tw-F^TB^Tw=A^Tw=\lambda w }

を得ますが，このようなwは，仮定(1)より，w=0でなければなりません。

状態フィードバック

投稿日時: 投稿者:

状態フィードバック…Homework

[1]　ある平衡状態周りの挙動が線形状態方程式

\displaystyle{(1)\quad \dot{x}(t)=Ax(t)+Bu(t)\quad(x(t)\in{\rm\bf R}^n,u(t)\in{\rm\bf R}^m) }

で表される制御対象に対して、状態フィードバック

\displaystyle{(2)\quad \boxed{u(t)=-Fx(t)} }

を考えます。ここで、mを入力変数の個数、nを状態変数の個数とするとき、Fはサイズm\times nのゲイン行列です。(2)を(1)に代入して、次式を得ます。

\displaystyle{(3)\quad \boxed{\dot{x}(t)=(A-BF)x(t)} }

以上の状況は次図のように図示されます。

図1 状態フィードバックによる閉ループ系

これは、制御対象とコントローラ（状態フィードバック）が閉路を構成することから、閉ループ系と呼ばれます。ここで、興味深いのは、制御対象が不安定であっても（行列Aの固有値の一部が複素左半平面にあっても）、状態フィードバックを行って閉ループ系を構成すると、行列A-BFの固有値をすべて複素左半平面に移動させる可能性があることです。どのような条件の下で、どのようにゲイン行列Fを決定できるかを検討する問題を状態フィードバックによる安定化問題と言います。

●状態フィードバックゲイン行列Fを求めることは、適当な安定行列A_{F}を与えて、行列方程式

\displaystyle{(4)\quad A-BF=A_F\Leftrightarrow BF=A-A_F }

を解くことに他なりません。よく学生さんは次のように解きます。

\displaystyle{(5)\quad F=B^{-1}(A-A_F) }

でも一般にはBは正方行列ではないので、そもそも逆行列を考えることができません。A_Fを安定行列とする状態フィードバックが求まる条件の説明は次節で行いますが、ここではその準備として、状態フィードバックゲイン行列Fを求める公式を示します。

[2]　いま行列A_F=A-BFの特性多項式を次式で表します。

\displaystyle{(6)\quad {\rm det}(\lambda I_n-A_F)=\lambda^n+a_1'\lambda^{n-1}+\cdots+a_n' %=(\lambda-\lambda_1')\times\cdots\times(\lambda-\lambda_n') }

ケーリーハミルトンの定理から次式が成り立ちます。

\displaystyle{(7)\quad A_F^n+a_1'A_F^{n-1}+\cdots+a_{n-1}'A_F+a_n'I_n=0_{n\times n} }

これはn=2のとき、次のように書けます。

\displaystyle{(8)\quad \begin{array}{l} A_F^2+a_1'A_F+a_2'I_2=0\\ \quad\Downarrow\\ (A-BF)A_F+a_1'(A-BF)+a_2'I_2=0\\ \quad\Downarrow\\ AA_F-BFA_F+a_1'A-a_1'BF+a_2'I_2=0\\ \quad\Downarrow\\ A(A-BF)+a_1'A+a_2'I_2=BFA_F+a_1'BF\\ \quad\Downarrow\\ A^2+a_1'A+a_2'I_2=BF(A_F+a_1'I_2)+ABF\\ \quad\Downarrow\\ A^2+a_1'A+a_2'I_2=\left[\begin{array}{ccc} B &AB \end{array}\right] \left[\begin{array}{ccc} F(A_F+a_1'I_2)\\ F \end{array}\right] \end{array} }

したがって、１入力系で\left[\begin{array}{ccc} B &AB \end{array}\right]が正則ならば、次の公式を得ます。

\displaystyle{(9)\quad \boxed{F= \left[\begin{array}{ccc} 0 &1 \end{array}\right] \left[\begin{array}{ccc} B &AB \end{array}\right]^{-1} (A^2+a_1'A+a_2'I_2)} }

同様に、n=3のとき、次のように書けます。

\displaystyle{(10)\quad \begin{array}{l} A_F^3+a_1'A_F^2+a_2'A_F+a_3'I_3=0\\ \quad\Downarrow\\ (A-BF)(A_F^2+a_1'A_F+a_2'I_3)+a_3'I_3=0\\ \quad\Downarrow\\ A(A_F^2+a_1'A_F+a_2'I_3)+a_3'I_3\\=BF(A_F^2+a_1'A_F+a_2'I_3)\\ \quad\Downarrow\\ A((A-BF)A_F+a_1'(A-BF)+a_2'I_3)+a_3'I_3\\=BF(A_F^2+a_1'A_F+a_2'I_3)\\ \quad\Downarrow\\ A^2(A-BF)+a_1'A^2+a_2'A+a_3'I_3\\ =BF(A_F^2+a_1'A_F+a_2'I_3)+ABF(A_F+a_1'I_3)\\ \quad\Downarrow\\ A^3+a_1'A^2+a_2'A+a_3'I_3\\ =BF(A_F^2+a_1'A_F+a_2'I_3)+ABF(A_F+a_1'I_3)+A^2BF\\ \quad\Downarrow\\ A^3+a_1'A^2+a_2'A+a_3'I_3\\ =\left[\begin{array}{ccc} B &AB&A^2B \end{array}\right] \left[\begin{array}{ccc} F(A_F^2+a_1'A_F+a_2'I_3)\\ F(A_F+a_1'I_3)\\ F \end{array}\right] \end{array} }

したがって、１入力系で\left[\begin{array}{ccc} B &AB&A^2B \end{array}\right]が正則ならば、次の公式を得ます。

\displaystyle{(11)\quad F= \left[\begin{array}{ccc} 0 &0 &1 \end{array}\right] \left[\begin{array}{ccc} B &AB&A^2B \end{array}\right]^{-1} (A^3+a_1'A^2+a_2'A+a_3'I_3) }

n>3の場合も同様にして、A_F=A-BFの特性多項式(7)を

\displaystyle{(12)\quad %\begin{array}{l} A^n+a'_1A^{n-1}+\cdots+a'_nI_n =\left[\begin{array}{cccc} B & AB &\cdots & A^{n-1}B \end{array}\right] \left[\begin{array}{ccc} X_{n-1}\\ \vdots\\ X_1\\ F \end{array}\right] %\end{array} }

のように書き換えることができます。ここで、X_1,\cdots,X_{n-1}はサイズm\times nの適当な行列です。したがって、１入力系で\left[\begin{array}{cccc} B & AB &\cdots & A^{n-1}B \end{array}\right]が正則ならば、(12)の右辺の第２番目の行列の最下段にあるFを取りだすことにより、次の公式が得られます。

\displaystyle{(13)\quad \boxed{\begin{array}{l} F= \left[\begin{array}{cccc} 0 & \cdots & 0 & 1 \end{array}\right] \left[\begin{array}{cccc} B & AB &\cdots & A^{n-1}B \end{array}\right]^{-1}\\ \times (A^n+a'_1A^{n-1}+\cdots+a'_nI_n) \end{array}} }

このように、１入力n次系の場合は、Fのn個の要素が、A-BFの特性多項式のn個の係数（すなわちn個の固有値）を指定することにより決定されます。ところがm入力系の場合、A_Fの固有値だけでは決まらず、対応する固有ベクトルまで指定する必要があります。このことを以下で調べてみます。

[3]　いま、ジョルダン分解A_F=V\Lambda V^{-1}を仮定して

\displaystyle{(14)\quad A-BF=\underbrace{V\Lambda V^{-1}}_{A_F} \Rightarrow AV-V\Lambda=B\underbrace{FV}_{G} }

ここで、\Lambdaが対角行列の場合、A_Fの固有値\lambda_1',\cdots,\lambda_n'に対応する固有ベクトルをv_1,\cdots,v_nとして

\displaystyle{(15)\quad A \underbrace{ \left[\begin{array}{ccc} v_1&\cdots&v_n \end{array}\right] }_{V} - \underbrace{ \left[\begin{array}{ccc} v_1&\cdots&v_n \end{array}\right] }_{V} \underbrace{ {\rm diag}\{\lambda_1',\cdots,\lambda_m'\} }_{\Lambda}=}

\displaystyle{ \left[\begin{array}{ccc} (A-\lambda_1'I_n)v_1&\cdots&(A-\lambda_n'I_n)v_n \end{array}\right] =B \underbrace{ \left[\begin{array}{ccc} g_1&\cdots&g_n \end{array}\right] }_{G} }

これより、もし\lambda_1',\cdots,\lambda_n'がAの固有値と一致しないならば、固有ベクトルの表現式として次式を得ます。

\displaystyle{(16)\quad \underbrace{ \left[\begin{array}{ccc} v_1&\cdots&v_n \end{array}\right] }_{V}= \left[\begin{array}{ccc} (A-\lambda_1'I_n)^{-1}Bg_1&\cdots&(A-\lambda_n'I_n)^{-1}Bg_n \end{array}\right] }

したがって、状態フィードバックゲイン行列Fを求める次の手順が考えられます。

1) A_Fの固有値\lambda_1',\cdots,\lambda_n'をAの固有値とは異なるように指定する。

2) A_Fの固有ベクトルv_1,\cdots,v_nを決める自由度g_1,\cdots,g_nを、(15)のVが正則になるように与える。

3) 次式で Fを求める。

\displaystyle{(17)\quad \boxed{F=GV^{-1}} }

この手順は簡潔ですが、１）の制約があること、２）の具体的な手法が不明であることなどの難点があります。状態フィードバックを定める手続きは、ある最適化問題の解として得ることが多く、これについてはあと説明します。

演習 A31…Flipped Classroom

1^\circ 次のコードを実行し、A-BFの固有値を計算し、指定されているものになっていることを確かめよ。

MATLAB 
%a31_1.m
%-----
 clear all, close all
 A=[0 1;0 0]; B=[0;1]; p=[-1 -1];
 F=sf(A,B,p)
 lambda=eig(A-B*F)
%----- 
 A=[1 2;3 4]; B=[5;6]; p=[-1+i -1-i];
 F=sf(A,B,p)
 lambda=eig(A-B*F)
%----- 
 function F=sf(A,B,p)
   n=size(A,1); r=eig(A);
   a=[1 -r(1)]; for i=2:n, a=conv(a,[1 -r(i)]); end, a=a(n+1:-1:2);
   b=[1 -p(1)]; for i=2:n, b=conv(b,[1 -p(i)]); end, b=b(n+1:-1:2);
   X=zeros(n,n); for i=1:n, X(i,1:n-i+1)=[a(i+1:n) 1]; end
   Y=B; for i=2:n, Y=[B A*Y]; end
   F=real((b-a)/X/Y);
 end
%-----
%eof

2^\circ 次のコードを実行し、A-BFの固有値と固有ベクトルを計算し、指定されているものになっていることを確かめよ。

MATLAB 
%a31_2.m
%-----
 clear all, close all
 A=[0 0;0 -1]; B=[1 1;1 -1]; p=[-2 -3]; V=[1 0;0 1];
 [n,m]=size(B);
 G=[];
 for i=1:n, G=[G ((A-p(i)*eye(n,n))\B)\V(:,i)]; end
 F=G/V
 [V,R]=eig(A-B*F)
%-----
%eof

Note A31 不可制御の例
どのような制御対象に対しても，閉ループ系を安定化する状態フィードバックが求まるわけではありません。たとえば、皆さんは小さいころ棒立て遊びをやった記憶があると思います。不安定なモノを安定化できた楽しい思い出だったと思います。それでは、２本の棒を同時に立てることはできたと思いますか？この状況を次図のように表します。


図 ２本の棒立て制御

ここで、２本の棒の長さと質量をそれぞれ\ell_1,\ell_2とm_1,m_2とします。また鉛直線からの傾きを\theta_1,\theta_2とします。２本の棒は１つの軸に取り付けられており、自由に回転できるとします。この軸を左右に動かして回転トルクを与えるのですが、その結果として働くトルクをここではh_1\tau, h_2\tauで表します。このとき運動方程式は次式となります。ただし、J_1=\frac{4}{3}m_1\ell_1^2, J_2=\frac{4}{3}m_2\ell_2^2。

\displaystyle{(1)\quad \left\{\begin{array}{l} J_1\ddot\theta_1(t)=m_1g\ell_1\theta_1(t)+h_1\tau(t) \\ J_2\ddot\theta_2(t)=m_1g\ell_2\theta_2(t)+h_2\tau(t) \end{array}\right. }

これから次の状態方程式を得ます。

\displaystyle{(2)\quad \begin{equation} \frac{d}{dt} \left[\begin{array}{c} \theta_1\\ \theta_2\\ \dot\theta_1\\ \dot\theta_2\\ \end{array}\right]= \left[\begin{array}{cccc} 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ a_1 & 0 & 0 & 0\\ 0 & a_2 & 0 & 0\\ \end{array}\right] \left[\begin{array}{c} \theta_1\\ \theta_2\\ \dot\theta_1\\ \dot\theta_2\\ \end{array}\right]+ \left[\begin{array}{c} 0\\ 0\\ b_1\\ b_2\\ \end{array}\right]\tau }

ただし、a_1=\frac{3g}{4\ell_1}, a_2=\frac{3g}{4\ell_2}、b_1=\frac{h_1}{J_1}, b_2=\frac{h_2}{J_2}。

それでは、これを安定化する状態フィードバックゲインを求めてみます。上述の公式

\displaystyle{(3)\quad F= \left[\begin{array}{cccc} 0 & \cdots & 0 & 1 \end{array}\right] %\underbrace{ \left[\begin{array}{cccc} B & AB &\cdots & A^{n-1}B \end{array}\right]^{-1} %}_{可制御性行列} }

を参照すれば、次の行列が正則である必要があります。

\displaystyle{(4)\quad \begin{equation} \left[\begin{array}{cccc} B & AB & A^2B & A^3B \end{array}\right] = \left[\begin{array}{cccc} 0 &b_1 &0 &a_1b_1\\ 0 &b_2 &0 &a_2b_2\\ b_1&0 &a_1b_1 &0\\ b_2&0 &a_2b_2 &0\\ \end{array}\right] }

これは２つの棒の長さ違うときは正則となり、閉ループ系を安定化する状態フィードバックが求まります。ところが２つの棒が同じで、a_1=a_2, b_1=b_2となる場合は、第１行と第２行が同じとなり、第３行と第４行も同じとなり、正則性は失われます。これから全く同じ２本の棒の棒立ては成功しないと考えられます。これは「二兎を追うものは一兎も得ず」のことわざ通りですね。

２次系の周波数応答

投稿日時: 投稿者:

２次系の周波数応答…Homework

[1] 次のようにパラメタライズされた漸近安定な１入力１出力２次系の状態空間表現を考えます。

\displaystyle{(1a)\quad \underbrace{ \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{x}_2(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{cc} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{cc} 0 \\ K\omega_n^2 \end{array}\right] }_{B} u(t) }

\displaystyle{(1b)\quad y(t)= \underbrace{ \left[\begin{array}{cc} 1 & 0 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} }

ここで、\zeta<1とします。また以下では簡単のためK=1とします。このときインパルス応答は、次式で与えられました。

\displaystyle{(2a)\quad G(t)=\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\sin\lambda_I t}
\displaystyle{(2b)\quad \lambda_1,\lambda_2=\underbrace{-\zeta\omega_n}_{\lambda_R}\pm j\underbrace{\omega_n\sqrt{1-\zeta^2}}_{\lambda_I}}

この２次系に対して、正弦波入力

\displaystyle{(3)\quad u(t)=\sin\omega t}

を考えます。２次振動系（\zeta<1）の場合、零状態応答は次式で与えられます（導出は Note A24-2を参照）。

\displaystyle{(4a)\quad y(t)=\frac{1}{\sqrt{(1-(\frac{\omega}{\omega_n})^2)^2+4\zeta^2(\frac{\omega}{\omega_n})^2}}(\sin(\omega t+\theta)-e^{\lambda_Rt}\sin(\lambda_I t+\theta'))}
\displaystyle{(4b)\quad \theta=\tan^{-1}\frac{-2\zeta\frac{\omega}{\omega_n}}{1-(\frac{\omega}{\omega_n})^2}}
\displaystyle{(4c)\quad \theta'=\tan^{-1}\frac{-2\zeta\sqrt{1-\zeta^2}}{1-2\zeta^2-(\frac{\omega}{\omega_n})^2}}

ここで、t\rightarrow\inftyとすると

\displaystyle{(5)\quad \boxed{y(t)\simeq\underbrace{\frac{1}{\sqrt{(1-\frac{\omega^2}{\omega_n^2})^2+4\zeta^2\frac{\omega^2}{\omega_n^2}}}}_{|\hat{G}(j\omega)|} \sin(\omega t+\underbrace{\tan^{-1}\frac{-2\zeta\frac{\omega}{\omega_n}}{1-\frac{\omega^2}{\omega_n^2}}}_{\angle\hat{G}(j\omega)})} }

ただし

\displaystyle{(6)\quad \hat{G}(s)=\frac{\omega_n^2}{s^2+2\zeta\omega_ns+\omega_n^2}}

これは、入力が正弦波のときは、時間が十分立てば、出力も正弦波となることを示しています。その振幅と位相はそれぞれ\hat{G}(j\omega)の絶対値と偏角となっています。

|\hat{G}(j\omega)|をゲイン、\angle\hat{G}(j\omega)を位相と呼びます。このゲイン線図と位相線図をペアにして片対数グラフにプロットしたものをボード線図と呼びます。ゲインはdb値(20{\log_{10}|\hat{G}(j\omega)|)、位相はdeg値(\frac{180}{\pi}\angle\hat{G}(j\omega))です。

実際、２次系において\omega_n=1,\zeta=0.01,0.707,1のとき、ボード線図を描いてみると次のグラフが得られます

図1 ２次系の周波数応答の比較

周波数応答に基づく２次系の同定…Homework

[2] ２次振動系の場合、ゲイン線図

\displaystyle{(7)\quad |\hat{G}(j\omega)|=\frac{1}{\sqrt{(1-\frac{\omega^2}{\omega_n^2})^2+4\zeta^2\frac{\omega^2}{\omega_n^2}}} }

の頂点は次式で与えられます。

\displaystyle{(8)\quad \boxed{(\frac{\omega_p}{\omega_n},M_p)=\left(\sqrt{1-2\zeta^2},{1\over 2\zeta\sqrt{1-\zeta^2}}\right)\quad(\zeta\le{1\over\sqrt{2}})} }

ここで、頂点が現れる条件、\zeta\le{1\over\sqrt{2}}がついていることに注意してください。

実際、{d\over d\omega}|\hat G(j\omega_p)|=0を満足する\omega_pは

(9)\quad \begin{array}{l} \displaystyle{{d\over d\omega}|\hat G(j\omega)| =\left.-\frac{1}{2}\frac{2(1-(\frac{\omega}{\omega_n})^2)(-\frac{2\omega}{\omega_n^2})) +4\zeta^2\frac{2\omega}{\omega_n^2}}{((1-(\frac{\omega}{\omega_n})^2)^2+4\zeta^2(\frac{\omega}{\omega_n})^2)^{3/2}}\,\right|_{\omega=\omega_p}=0}\\ \displaystyle{\Rightarrow\ -(1-(\frac{\omega_p}{\omega_n})^2)+2\zeta^2=0}\\ \displaystyle{\Rightarrow\ (\frac{\omega_p}{\omega_n})^2=1-2\zeta^2} \end{array}

となって、次式のように表されます。

\displaystyle{(10)\quad \boxed{\frac{\omega_p}{\omega_n}=\sqrt{1-2\zeta^2}}}

また、このときの最大値は

(11)\quad \begin{array}{l} \displaystyle{M_p=|\hat G(j\omega_p)|=\frac{1}{\sqrt{(2\zeta^2)^2+4\zeta^2(1-2\zeta^2)}}}\\ \displaystyle{=\frac{1}{\sqrt{4\zeta^4+4\zeta^2-8\zeta^4}}=\frac{1}{\sqrt{4\zeta^2(1-\zeta^2)}}} \end{array}

となって、次式のように表されます。。

\displaystyle{(12)\quad \boxed{M_p=\frac{1}{2\zeta\sqrt{1-\zeta^2}}}}

(10)と(12)から、２次振動系はゲイン線図が頂点を持つ場合、その座標(\frac{\omega_p}{\omega_n},M_p)を求めて、減衰係数と固有角周波数(\zeta,\omega_n)が得られることを示しています。

●ちなみに、|\hat{G}(j\omega_b)|=\frac{|\hat{G}(j0)|}{\sqrt{2}}を満足する帯域幅\omega_bは次のように計算されます。

\displaystyle{(13)\quad \frac{1}{\sqrt{(1-(\frac{\omega_b}{\omega_n})^2)^2+4\zeta^2(\frac{\omega_b}{\omega_n})^2}} =\frac{1}{\sqrt{2}} }

において、\omega'=\frac{\omega_b}{\omega_n}とおくと

(14)\quad \begin{array}{l} \displaystyle{(1-\omega'\,^2)^2+4\zeta^2\omega'\,^2=1-2\omega'\,^2+\omega'\,^4+4\zeta^2\omega'\,^2=2}\\ \displaystyle{\Rightarrow\ \omega'^4-2(1-2\zeta^2)\omega'\,^2-1=0}\\ \displaystyle{\Rightarrow\ \omega'^2=1-2\zeta^2\pm\sqrt{4\zeta^4-4\zeta^2+2}} \end{array}

となって、次式のように表されます。

\displaystyle{(15)\quad \boxed{\frac{\omega_b}{\omega_n}=\sqrt{1-2\zeta^2\pm\sqrt{4\zeta^4-4\zeta^2+2}}１}}

演習A24…Flipped Classroom

図1のゲイン曲線#1を描け。

MATLAB 
%a24.m
%-----
 clear all, close all
 A1=[0 1;-1 -2*0.01]; 
 A2=[0 1;-1 -2/sqrt(2)]; 
 A3=[0 1;-1 -2*1]; 
 B=[0;1]; C=[1 0]; D=[];
 sys1=ss(A1,B,C,D);
 sys2=ss(A2,B,C,D);
 sys3=ss(A3,B,C,D);
 w=logspace(-1,1,1000);
 bode(sys1,sys2,sys3,w)
 grid
 title('Bode Diagrams of 2nd-order System')
 xlabel('Freq')
%ylabel('Gain')
 ylabel('Phase')
 legend('zeta=0.01','zeta=0.707','zeta=1')
%-----
%eof
SCILAB 
coming soon

Note A24-1　n次系の周波数応答

次の漸近安定な１入力１出力n次系の状態空間表現を考えます。

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

このとき、正弦波入力

\displaystyle{(2)\quad u(t)=\sin\omega t}

に対する零状態応答を計算します。そのために

\displaystyle{(3)\quad u(t)=e^{j\omega t}=\cos(\omega t)+j\sin(\omega t)}

を、零状態応答の式

\displaystyle{(4)\quad y(t)=\int_0^tC\exp(A(t-\tau))Bu(\tau)d\tau}

に代入して

(5)\quad \begin{array}{l} \displaystyle{y(t)=\int_0^tC\exp(A(t-\tau))Be^{j\omega\tau}d\tau}\\ \displaystyle{=C\exp(At)\int_0^te^{j\omega\tau}\exp(-A\tau)Bd\tau}\\ \displaystyle{=C\exp(At)\int_0^t\exp(j\omega\tau I_n)\exp(-A\tau)Bd\tau}\\ \displaystyle{=C\exp(At)\int_0^t\exp((j\omega I_n-A)\tau)Bd\tau}\\ \displaystyle{=C\exp(At) \left[\frac{}{}\exp((j\omega I_n-A)\tau)\right]_0^t(j\omega I_n-A)^{-1}B}\\ \displaystyle{=C\exp(At) (\exp((j\omega I_n-A)t)-I_n)(j\omega I_n-A)^{-1}B}\\ \displaystyle{=-C\exp(At)(j\omega I_n-A)^{-1}B+C\exp(At)\exp(j\omega t I_n)\exp(-At)(j\omega I_n-A)^{-1}B}\\ \displaystyle{=-C\exp(At)(j\omega I_n-A)^{-1}B+C(j\omega I_n-A)^{-1}Be^{j\omega t}} \end{array}

ここでt\rightarrow\inftyとすると

(6)\quad \begin{array}{l} \displaystyle{y(t)=C\underbrace{\exp(At)}_{\rightarrow 0\ (t\rightarrow\infty)}(j\omega I_n-A)^{-1}B+\underbrace{C(j\omega I_n-A)^{-1}B}_{G(j\omega) =|G(j\omega)|e^{j\angle G(j\omega)}}e^{j\omega t}}\\ \displaystyle{\simeq |G(j\omega)|e^{j(\omega t+\angle G(j\omega))}}\\ \displaystyle{=|G(j\omega)|\cos(\omega t+\angle G(j\omega))+j|G(j\omega)|\sin(\omega t+\angle G(j\omega))} \end{array}

これから正弦波入力(3)に対する零状態応答は、t\rightarrow\inftyのとき次式で与えられます。

\displaystyle{(7)\quad y(t)\simeq|G(j\omega)|\sin(\omega t+\angle G(j\omega))}

これは、入力が正弦波のときは、時間が十分立てば、出力も正弦波となることを示しています。その振幅と位相はそれぞれ\hat{G}(j\omega)の絶対値と偏角となっています。

Note A24-2　２次系の周波数応答

次の２次系を考えます。

\displaystyle{(1a)\quad \underbrace{ \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{x}_2(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{cc} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{cc} 0 \\ \omega_n^2 \end{array}\right] }_{B} u(t) }

\displaystyle{(1b)\quad y(t)= \underbrace{ \left[\begin{array}{cc} 1 & 0 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} }

ここで、\zeta<1とします。このときインパルス応答は

\displaystyle{(2a)\quad G(t)=\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\sin\lambda_I t}
\displaystyle{(2b)\quad \lambda_1,\lambda_2=\underbrace{-\zeta\omega_n}_{\lambda_R}\pm j\underbrace{\omega_n\sqrt{1-\zeta^2}}_{\lambda_I}}

でした。いま、正弦波入力

\displaystyle{(3)\quad u(t)=\sin\omega t}

に対する時間応答を考えます。このとき、公式

\displaystyle{(4a)\quad \int e^{ax}\sin px \cos qx\,dx}
\displaystyle{=\frac{e^{ax}}{2} \left(\frac{a\sin(p+q)x-(p+q)\cos(p+q)x}{a^2+(p+q)^2}\right.+\left.\frac{a\sin(p-q)x-(p-q)\cos(p-q)x}{a^2+(p-q)^2}\right)}
\displaystyle{(4b)\quad \int e^{ax}\sin px \sin qx\,dx}
\displaystyle{=\frac{e^{ax}}{2} \left(-\frac{a\cos(p+q)x+(p+q)\sin(p+q)x}{a^2+(p+q)^2}\right.+\left.\frac{a\cos(p-q)x+(p-q)\sin(p-q)x}{a^2+(p-q)^2}\right)}

を用いて、周波数応答が次のように得られます。

(5)\quad \begin{array}{l} \displaystyle{y(t)=\int_0^t\frac{\omega_n^2}{\lambda_I}e^{\lambda_R(t-\tau)}\sin\lambda_I(t-\tau)\sin\omega\tau\,d\tau}\\ \displaystyle{=\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\sin\lambda_It\int_0^te^{-\lambda_R\tau}\sin\omega\tau\cos\lambda_I\tau\,d\tau}\\ \displaystyle{-\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\cos\lambda_It\int_0^te^{-\lambda_R\tau}\sin\omega\tau\sin\lambda_I\tau\,d\tau}\\ \displaystyle{=\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\sin\lambda_It\times}\\ \displaystyle{(\left[\frac{e^{-\lambda_R\tau}(-\lambda_R\sin(\omega+\lambda_I)\tau-(\omega+\lambda_I)\cos(\omega+\lambda_I)\tau)}{2(\lambda_R^2+(\omega+\lambda_I)^2)}\right.}\\ \displaystyle{+\left.\frac{e^{-\lambda_R\tau}(-\lambda_R\sin(\omega-\lambda_I)\tau-(\omega-\lambda_I)\cos(\omega-\lambda_I)\tau)}{2(\lambda_R^2+(\omega-\lambda_I)^2)}\right]_0^t)}\\ \displaystyle{-\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\cos\lambda_It\times}\\ \displaystyle{(\left[-\frac{e^{-\lambda_R\tau}(-\lambda_R\cos(\omega+\lambda_I)\tau+(\omega+\lambda_I)\sin(\omega+\lambda_I)\tau)}{2(\lambda_R^2+(\omega+\lambda_I)^2)}\right.}\\ \displaystyle{+\left.\frac{e^{-\lambda_R\tau}(-\lambda_R\cos(\omega-\lambda_I)\tau+(\omega-\lambda_I)\sin(\omega-\lambda_I)\tau)}{2(\lambda_R^2+(\omega-\lambda_I)^2)^2}\right]_0^t)}\\ \displaystyle{=\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\sin\lambda_It\times}\\ \displaystyle{(\frac{e^{-\lambda_Rt}(-\lambda_R\sin(\omega+\lambda_I)t-(\omega+\lambda_I)\cos(\omega+\lambda_I)t)+\omega+\lambda_I}{2(\lambda_R^2+(\omega+\lambda_I)^2)}}\\ \displaystyle{+\frac{e^{-\lambda_Rt}(-\lambda_R\sin(\omega-\lambda_I)t-(\omega-\lambda_I)\cos(\omega-\lambda_I)t)+\omega-\lambda_I}{2(\lambda_R^2+(\omega-\lambda_I)^2)})}\\ \displaystyle{-\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\cos\lambda_It\times}\\ \displaystyle{(\frac{-e^{-\lambda_Rt}(-\lambda_R\cos(\omega+\lambda_I)t+(\omega+\lambda_I)\sin(\omega+\lambda_I)t)-\lambda_R}{2(\lambda_R^2+(\omega+\lambda_I)^2)}}\\ \displaystyle{+\frac{e^{-\lambda_Rt}(-\lambda_R\cos(\omega-\lambda_I)t+(\omega-\lambda_I)\sin(\omega-\lambda_I)t)+\lambda_R}{2(\lambda_R^2+(\omega-\lambda_I)^2)})} \end{array}

さらに変形すると

(6)\quad \begin{array}{l} \displaystyle{y(t)=\frac{\omega_n^2}{2\lambda_I}\sin\lambda_It\times\frac{-\lambda_R\sin(\omega+\lambda_I)t-(\omega+\lambda_I)\cos(\omega+\lambda_I)t}{\lambda_R^2+(\omega+\lambda_I)^2}}\\ \displaystyle{-\frac{\omega_n^2}{2\lambda_I}\cos\lambda_It\times\frac{-(-\lambda_R\cos(\omega+\lambda_I)t+(\omega+\lambda_I)\sin(\omega+\lambda_I)t)}{\lambda_R^2+(\omega+\lambda_I)^2}}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}\sin\lambda_It\times\frac{-\lambda_R\sin(\omega-\lambda_I)t-(\omega-\lambda_I)\cos(\omega-\lambda_I)t}{\lambda_R^2+(\omega-\lambda_I)^2}}\\ \displaystyle{-\frac{\omega_n^2}{2\lambda_I}\cos\lambda_It\times\frac{-\lambda_R\cos(\omega-\lambda_I)t+(\omega-\lambda_I)\sin(\omega-\lambda_I)t}{\lambda_R^2+(\omega-\lambda_I)^2}}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}e^{\lambda_Rt}\sin\lambda_It\times(\frac{\omega+\lambda_I}{\lambda_R^2+(\omega+\lambda_I)^2}+\frac{\omega-\lambda_I}{\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{-\frac{\omega_n^2}{2\lambda_I}e^{\lambda_Rt}\cos\lambda_It\times(\frac{-\lambda_R}{\lambda_R^2+(\omega+\lambda_I)^2}+\frac{\lambda_R}{\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{=\frac{\omega_n^2}{2\lambda_I}\times\frac{-\lambda_R\sin(\omega+\lambda_I)t\sin\lambda_It-(\omega+\lambda_I)\cos(\omega+\lambda_I)t\sin\lambda_It}{\lambda_R^2+(\omega+\lambda_I)^2}}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}\times\frac{-\lambda_R\cos(\omega+\lambda_I)t\cos\lambda_It+(\omega+\lambda_I)\sin(\omega+\lambda_I)t\cos\lambda_It}{\lambda_R^2+(\omega+\lambda_I)^2}}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}\times\frac{-\lambda_R\sin(\omega-\lambda_I)t\sin\lambda_It-(\omega-\lambda_I)\cos(\omega-\lambda_I)t\sin\lambda_It}{\lambda_R^2+(\omega-\lambda_I)^2}}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}\times\frac{+\lambda_R\cos(\omega-\lambda_I)t\cos\lambda_It-(\omega-\lambda_I)\sin(\omega-\lambda_I)t\cos\lambda_It}{\lambda_R^2+(\omega-\lambda_I)^2}}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}e^{\lambda_Rt}\sin\lambda_It\times(\frac{\omega+\lambda_I}{\lambda_R^2+(\omega+\lambda_I)^2}+\frac{\omega-\lambda_I}{\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}e^{\lambda_Rt}\cos\lambda_It\times(\frac{\lambda_R}{\lambda_R^2+(\omega+\lambda_I)^2}-\frac{\lambda_R}{\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{=\frac{\omega_n^2}{2\lambda_I}\times\frac{-\lambda_R\sin(\omega+\lambda_I)t\sin\lambda_It-\lambda_R\cos(\omega+\lambda_I)t\cos\lambda_It}{\lambda_R^2+(\omega+\lambda_I)^2}}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}\times\frac{-(\omega+\lambda_I)\cos(\omega+\lambda_I)t\sin\lambda_It+(\omega+\lambda_I)\sin(\omega+\lambda_I)t\cos\lambda_It}{\lambda_R^2+(\omega+\lambda_I)^2}}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}\times\frac{-\lambda_R\sin(\omega-\lambda_I)t\sin\lambda_It+\lambda_R\cos(\omega-\lambda_I)t\cos\lambda_It}{\lambda_R^2+(\omega-\lambda_I)^2}}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}\times\frac{-(\omega-\lambda_I)\cos(\omega-\lambda_I)t\sin\lambda_It-(\omega-\lambda_I)\sin(\omega-\lambda_I)t\cos\lambda_It}{\lambda_R^2+(\omega-\lambda_I)^2}}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}e^{\lambda_Rt}\sin\lambda_It\times(\frac{\omega+\lambda_I}{\lambda_R^2+(\omega+\lambda_I)^2}+\frac{\omega-\lambda_I}{\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}e^{\lambda_Rt}\cos\lambda_It\times(\frac{\lambda_R}{\lambda_R^2+(\omega+\lambda_I)^2}-\frac{\lambda_R}{\lambda_R^2+(\omega-\lambda_I)^2})} \end{array}

これを、公式

\displaystyle{(7a)\quad \sin(\alpha\pm\beta)=\sin\alpha\cos\beta\pm\cos\alpha\sin\beta}
\displaystyle{(7b)\quad \cos(\alpha\pm\beta)=\cos\alpha\cos\beta\mp\sin\alpha\sin\beta}

を用いて、以下のように変形します。

(8)\quad \begin{array}{l} \displaystyle{y(t)=\frac{\omega_n^2}{2\lambda_I}\times\frac{-\lambda_R\cos\omega t}{\lambda_R^2+(\omega+\lambda_I)^2}+\frac{\omega_n^2}{2\lambda_I}\times\frac{(\omega+\lambda_I)\sin\omega t}{\lambda_R^2+(\omega+\lambda_I)^2}}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}\times\frac{\lambda_R\cos\omega t}{\lambda_R^2+(\omega-\lambda_I)^2}+\frac{\omega_n^2}{2\lambda_I}\times\frac{-(\omega-\lambda_I)\sin\omega t}{\lambda_R^2+(\omega-\lambda_I)^2}}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}e^{\lambda_Rt}\sin\lambda_It\times(\frac{\omega+\lambda_I}{\lambda_R^2+(\omega+\lambda_I)^2}+\frac{\omega-\lambda_I}{\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{+\frac{\omega_n^2}{2\lambda_I}e^{\lambda_Rt}\cos\lambda_It\times(\frac{\lambda_R}{\lambda_R^2+(\omega+\lambda_I)^2}-\frac{\lambda_R}{(\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{=\sin\omega t\times\frac{\omega_n^2}{2\lambda_I}(\frac{\omega+\lambda_I}{\lambda_R^2+(\omega+\lambda_I)^2}-\frac{\omega-\lambda_I}{\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{-\cos\omega t\times\frac{\omega_n^2}{2\lambda_I}(\frac{\lambda_R}{\lambda_R^2+(\omega+\lambda_I)^2}-\frac{\lambda_R}{\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{+e^{\lambda_Rt}\sin\lambda_It\times\frac{\omega_n^2}{2\lambda_I}(\frac{\omega+\lambda_I}{\lambda_R^2+(\omega+\lambda_I)^2}+\frac{\omega-\lambda_I}{\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{+e^{\lambda_Rt}\cos\lambda_It\times\frac{\omega_n^2}{2\lambda_I}(\frac{\lambda_R}{\lambda_R^2+(\omega+\lambda_I)^2}-\frac{\lambda_R}{\lambda_R^2+(\omega-\lambda_I)^2})} \end{array}

ここで

(9)\quad \begin{array}{l} \displaystyle{\frac{\omega_n^2}{2\lambda_I}(\frac{\omega+\lambda_I}{\lambda_R^2+(\omega+\lambda_I)^2}-\frac{\omega-\lambda_I}{\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{=\frac{\omega_n^2}{2\lambda_I}(\frac{\omega+\lambda_I}{\omega_n^2+\omega^2+2\omega\lambda_I}-\frac{\omega-\lambda_I}{\omega_n^2+\omega^2-2\omega\lambda_I})}\\ \displaystyle{=\frac{\omega_n^2}{2\lambda_I} \frac{(\omega+\lambda_I)(\omega_n^2+\omega^2-2\omega\lambda_I)-(\omega-\lambda_I)(\omega_n^2+\omega^2+2\omega\lambda_I)} {(\omega_n^2+\omega^2)^2-4\omega^2\lambda_I^2}}\\ \displaystyle{=\frac{\omega_n^2}{2\lambda_I} \frac{-4\omega^2\lambda_I+2\lambda_I(\omega_n^2+\omega^2)} {(\omega_n^2+\omega^2)^2-4\omega^2\lambda_I^2}}\\ \displaystyle{=\frac{\omega_n^2}{1} \frac{-2\omega^2+(\omega_n^2+\omega^2)} {(\omega_n^2+\omega^2)^2-4\omega^2\omega_n^2(1-\zeta^2)}}\\ \displaystyle{=\frac{\omega_n^2(\omega_n^2-\omega^2)} {(\omega_n^2-\omega^2)^2+4\omega^2\omega_n^2\zeta^2}}\\ \displaystyle{=\frac{1-(\frac{\omega}{\omega_n})^2}{(1-(\frac{\omega}{\omega_n})^2)^2+4\zeta^2(\frac{\omega}{\omega_n})^2}} \end{array}

(10)\quad \begin{array}{l} \displaystyle{\frac{\omega_n^2}{2\lambda_I}(\frac{\omega+\lambda_I}{\lambda_R^2+(\omega+\lambda_I)^2}+\frac{\omega-\lambda_I}{\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{=\frac{\omega_n^2}{2\lambda_I}(\frac{\omega+\lambda_I}{\omega_n^2+\omega^2+2\omega\lambda_I}+\frac{\omega-\lambda_I}{\omega_n^2+\omega^2-2\omega\lambda_I})}\\ \displaystyle{=\frac{\omega_n^2}{2\lambda_I} \frac{(\omega+\lambda_I)(\omega_n^2+\omega^2-2\omega\lambda_I)+(\omega-\lambda_I)(\omega_n^2+\omega^2+2\omega\lambda_I)} {(\omega_n^2+\omega^2)^2-4\omega^2\lambda_I^2}}\\ \displaystyle{=\frac{\omega_n^2}{2\lambda_I} \frac{2\omega(\omega_n^2+\omega^2)-4\omega\lambda_I^2} {(\omega_n^2+\omega^2)^2-4\omega^2\lambda_I^2}}\\ \displaystyle{=\frac{\omega_n^2}{\lambda_I} \frac{\omega(\omega_n^2+\omega^2)-2\omega\omega_n^2(1-\zeta^2)} {(\omega_n^2+\omega^2)^2-4\omega^2\omega_n^2(1-\zeta^2)}}\\ \displaystyle{=\frac{\omega_n^2}{\lambda_I} \frac{\omega^3-\omega\omega_n^2(1-2\zeta^2)}{(\omega_n^2-\omega^2)^2+4\omega^2\omega_n^2\zeta^2}}\\ \displaystyle{=\frac{\omega}{\lambda_I}\frac{(\frac{\omega}{\omega_n})^2-(1-2\zeta^2)}{(1-(\frac{\omega}{\omega_n})^2)^2+4\zeta^2(\frac{\omega}{\omega_n})^2}} \end{array}

および

(11)\quad \begin{array}{l} \displaystyle{\frac{\omega_n^2}{2\lambda_I}(\frac{\lambda_R}{\lambda_R^2+(\omega+\lambda_I)^2}-\frac{\lambda_R}{\lambda_R^2+(\omega-\lambda_I)^2})}\\ \displaystyle{=\frac{\omega_n^2}{2\lambda_I}\lambda_R(\frac{1}{\omega_n^2+\omega^2+2\omega\lambda_I}-\frac{1}{\omega_n^2+\omega^2-2\omega\lambda_I})}\\ \displaystyle{=\frac{\omega_n^2}{2\lambda_I}\lambda_R \frac{(\omega_n^2+\omega^2-2\omega\lambda_I)-(\omega_n^2+\omega^2+2\omega\lambda_I)} {(\omega_n^2+\omega^2)^2-4\omega^2\lambda_I^2}}\\ \displaystyle{=\frac{\omega_n^2}{2\lambda_I}\lambda_R \frac{-4\omega\lambda_I} {(\omega_n^2+\omega^2)^2-4\omega^2\lambda_I^2}}\\ \displaystyle{=\frac{\omega_n^2}{1}(-\zeta\omega_n) \frac{-2\omega} {(\omega_n^2+\omega^2)^2-4\omega^2\omega_n^2(1-\zeta^2)}}\\ \displaystyle{=\frac{2\zeta\omega_n^3\omega} {(\omega_n^2-\omega^2)^2+4\omega^2\omega_n^2\zeta^2}}\\ \displaystyle{=\frac{2\zeta\frac{\omega}{\omega_n}}{(1-(\frac{\omega}{\omega_n})^2)^2+4\zeta^2(\frac{\omega}{\omega_n})^2}} \end{array}

だから

(12)\quad \begin{array}{l} \displaystyle{y(t)=\sin\omega t\times\frac{1-(\frac{\omega}{\omega_n})^2}{(1-(\frac{\omega}{\omega_n})^2)^2+4\zeta^2(\frac{\omega}{\omega_n})^2}}\\ \displaystyle{-\cos\omega t\times\frac{2\zeta\frac{\omega}{\omega_2}}{(1-(\frac{\omega}{\omega_n})^2)^2+4\zeta^2(\frac{\omega}{\omega_n})^2}}\\ \displaystyle{+e^{\lambda_Rt}\sin\lambda_It\times\frac{\omega}{\lambda_I}\frac{(\frac{\omega}{\omega_n})^2-(1-2\zeta^2)}{(1-(\frac{\omega}{\omega_n})^2)^2+4\zeta^2(\frac{\omega}{\omega_n})^2}}\\ \displaystyle{+e^{\lambda_Rt}\cos\lambda_It\times\frac{2\zeta\frac{\omega}{\omega_n}}{(1-(\frac{\omega}{\omega_n})^2)^2+4\zeta^2(\frac{\omega}{\omega_n})^2}} \end{array}

すなわち、周波数応答は

\displaystyle{(13a)\quad y(t)=\frac{1}{\sqrt{(1-(\frac{\omega}{\omega_n})^2)^2+4\zeta^2(\frac{\omega}{\omega_n})^2}}(\sin(\omega t+\theta)-e^{\lambda_Rt}\sin(\lambda_I t+\theta'))}
\displaystyle{(13b)\quad \theta=\tan^{-1}\frac{-2\zeta\frac{\omega}{\omega_n}}{1-(\frac{\omega}{\omega_n})^2}}
\displaystyle{(13c)\quad \theta'=\tan^{-1}\frac{2\zeta\frac{\omega}{\omega_n}}{\frac{\omega}{\lambda_I} ((\frac{\omega}{\omega_n})^2-(1-2\zeta^2))}=\tan^{-1}\frac{-2\zeta\sqrt{1-\zeta^2}}{1-2\zeta^2-(\frac{\omega}{\omega_n})^2}}

のように表されます。

２次系のステップ応答

投稿日時: 投稿者:

２次系のステップ応答…Homework

[0]　次のm入力p出力n次系の状態空間表現を考えます。

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

いま行列Bの第j列をB_j、行列Cの第i列をC_iで表すとき、第j番目の入力u_jから第i番目の出力y_iまでの時間応答は

\displaystyle{(2)\quad y_i(t)=\underbrace{C_i\exp(At)x(0)}_{zero-input\ response} +\underbrace{\int_0^tG_{ij}(t-\tau)u_j(\tau)d\tau}_{zero-state\ response}}

のように表されます。ここで、G_{ij}(t)は(1)のインパルス応答行列の(i,j)要素で、次のように与えられます。

\displaystyle{(3)\quad G_{ij}(t)=C_i\exp(At)B_j}

このときステップ入力

\displaystyle{(4)\quad u_j(t)= \left\{\begin{array}{ll} 0 & (t\le0)\\ 1 & (t>0) \end{array}\right.}

に対する零状態応答はステップ応答と呼ばれ、インパルス応答を積分して

\displaystyle{(5)\quad y_i(t)=\int_0^tG_{ij}(\tau)d\tau}

のように計算されます。

[1] さて、次のようにパラメタライズされた漸近安定な１入力１出力２次系の状態空間表現を考えます。

\displaystyle{(6a)\quad \underbrace{ \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{x}_2(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{cc} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{cc} 0 \\ K\omega_n^2 \end{array}\right] }_{B} u(t) }

\displaystyle{(6b)\quad y(t)= \underbrace{ \left[\begin{array}{cc} 1 & 0 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} }

以下では、簡単のためK=1とした２次系(7)に対するインパル応答を(6)から求めます。そのために、行列Aの固有値\lambda(A)=\{\lambda_1,\lambda_2\}を計算すると、\zeta>1、\zeta<1、\zeta=1に応じて、それぞれ次に示すようになります。

\displaystyle{(7)\quad \left\{\begin{array}{l} \lambda_1,\lambda_2=-\zeta\omega_n\pm\omega_n\sqrt{\zeta^2-1}\quad(\zeta>1)\\ \lambda_1,\lambda_2=\underbrace{-\zeta\omega_n}_{\lambda_R}\pm j\underbrace{\omega_n\sqrt{1-\zeta^2}}_{\lambda_I}\quad(\zeta<1)\\ \lambda_1=\lambda_2=\lambda=-\zeta\omega_n=-\omega_n\quad(\zeta=1) \end{array}\right. }

図1 A=\left[\begin{array}{cc} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right]の固有値分布

[3] ２次系のステップ応答(1)

\zeta>1のとき、\lambda_1,\lambda_2=-\zeta\omega_n\pm\omega_n\sqrt{\zeta^2-1}とおくと、インパルス応答とステップ応答はそれぞれ次式で与えられます。

\displaystyle{(8a)\quad G(t)=\frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}(e^{\lambda_2t}-e^{\lambda_1t})}
\displaystyle{(8b)\quad \boxed{S(t)=1+\frac{1}{\lambda_2-\lambda_1}(\lambda_1e^{\lambda_2t}-\lambda_2e^{\lambda_1t})}}

実際、インパルス応答を積分して、ステップ応答は次のように計算されます。

(9)\quad \begin{array}{l} \displaystyle{S(t)=\int_0^tG(\tau)d\tau=\int_0^t\frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}(e^{\lambda_2\tau}-e^{\lambda_1\tau})d\tau}\\ \displaystyle{= \frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}\left[\frac{e^{\lambda_2\tau}}{\lambda_2}-\frac{e^{\lambda_1\tau}}{\lambda_1}\right]_0^t}\\ \displaystyle{=\frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}(\frac{e^{\lambda_2t}}{\lambda_2}-\frac{e^{\lambda_1t}}{\lambda_1})-\frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}(\frac{1}{\lambda_2}-\frac{1}{\lambda_1})}\\ \displaystyle{=1+\frac{1}{\lambda_2-\lambda_1}(\lambda_1e^{\lambda_2t}-\lambda_2e^{\lambda_1t})} \end{array}

[4] ２次系のステップ応答(2)

\zeta<1のとき、\lambda_1,\lambda_2=\underbrace{-\zeta\omega_n}_{\lambda_R}\pm j\underbrace{\omega_n\sqrt{1-\zeta^2}}_{\lambda_I}とおくと、インパルス応答とステップ応答はそれぞれそれぞれ次式で与えられます。

\displaystyle{(10a)\quad G(t)=\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\sin\lambda_I t}
\displaystyle{(10b)\quad \boxed{S(t)=1-\frac{\omega_n}{\lambda_I}e^{\lambda_Rt}\sin(\lambda_I t+\tan^{-1}\frac{\lambda_I}{-\lambda_R})}}

実際、インパルス応答を積分して、公式

\displaystyle{(11)\quad \int e^{ax}\sin bx dx=\frac{e^{ax}}{a^2+b^2}(a\sin bx-b \cos bx)}

を用いて、ステップ応答は次のように計算されます。

(12)\quad \begin{array}{l} \displaystyle{S(t)=\int_0^tG(\tau)d\tau=\int_0^t\frac{\omega_n^2}{\lambda_I}e^{\lambda_R\tau}\sin\lambda_I \tau d\tau}\\ \displaystyle{=\frac{\omega_n^2}{\lambda_I}\left[ \frac{e^{\lambda_R\tau}}{\lambda_R^2+\lambda_I^2}(\lambda_R\sin\lambda_I\tau-\lambda_I\cos\lambda_I\tau)\right]_0^t}\\ \displaystyle{=\frac{\omega_n^2}{\lambda_I}\frac{1}{\omega_n^2} (e^{\lambda_Rt}(\lambda_R\sin\lambda_It-\lambda_I\cos\lambda_It)+\lambda_I)}\\ \displaystyle{=1-\frac{\omega_n}{\lambda_I}e^{\lambda_Rt} (\sin\lambda_It\times\frac{-\lambda_R}{\omega_n}+\cos\lambda_It\times\frac{\lambda_I}{\omega_n})}\\ \displaystyle{=1-\frac{\omega_n}{\lambda_I}e^{\lambda_Rt} (\sin\lambda_It\cos\phi+\cos\lambda_It\sin\phi)}\\ \displaystyle{=1-\frac{\omega_n}{\lambda_I}e^{\lambda_Rt}\sin(\lambda_It+\phi) \quad(\phi=\tan^{-1}\frac{\lambda_I}{-\lambda_R})} \end{array}

[5] ２次系のステップ応答(3)

\zeta=1のとき、\lambda_1=\lambda_2=\lambda=-\zeta\omega_n=-\omega_nとおくと、インパルス応答とステップ応答はそれぞれそれぞれ次式で与えられます。

\displaystyle{(13a)\quad G(t)=\lambda^2te^{\lambda t}}
\displaystyle{(13b)\quad \boxed{S(t)=1+(\lambda t-1)e^{\lambda t}}}

実際、インパルス応答を積分して、ステップ応答は次のように計算されます。

(14)\quad \begin{array}{l} \displaystyle{S(t)=\int_0^tG(\tau)d\tau=\int_0^t\lambda^2\tau e^{\lambda\tau}d\tau}\\ \displaystyle{=\lambda^2\left[\tau \frac{1}{\lambda}e^{\lambda\tau}\right]_0^t -\lambda^2\int_0^t \frac{1}{\lambda}e^{\lambda\tau}d\tau}\\ \displaystyle{=\lambda te^{\lambda t}-\lambda\left[\frac{1}{\lambda}e^{\lambda\tau}\right]_0^t}\\ \displaystyle{=1+(\lambda t-1)e^{\lambda t}} \end{array}

２次振動系の同定…Homework

[6]　２次系(6)において、\zeta=0.01,1/\sqrt{2},1、\omega_n=1、K=1の場合、ステップ応答のシミュレーション例を示します。

図2 ２次系のステップ応答の比較

\zeta<1の場合の(7)で表される２次系を２次振動系と呼びます。そのステップ応答

\displaystyle{(15)\quad S(t)=1-\frac{\omega_n}{\lambda_I}e^{\lambda_Rt}\sin(\lambda_I t+\phi)\quad(\phi=\tan^{-1}\frac{\lambda_I}{-\lambda_R})}

の第１オーバーシュートについて、次が成り立ちます。

\displaystyle{(16)\quad \boxed{(T_p,1+p_0)=\left(\frac{\pi}{\omega_n\sqrt{1-\zeta^2}},1+\exp(-\frac{\zeta\pi}{\sqrt{1-\zeta^2}})\right)} }

実際、\dot{S}(t)=0を計算すると

(17)\quad \begin{array}{l} \displaystyle{\dot{S}(t)=-\frac{\omega_n}{\lambda_I}(\lambda_Re^{\lambda_Rt}\sin(\lambda_It+\phi)+e^{\lambda_Rt}\cos(\lambda_It+\phi)\times\lambda_I)=0}\\ \displaystyle{\Rightarrow \sin(\lambda_It+\phi)\times\frac{-\lambda_R}{\omega_n}-\cos(\lambda_It+\phi)\times\frac{\lambda_I}{\omega_n}=0}\\ \displaystyle{\Rightarrow \sin(\lambda_It+\phi)\cos\phi-\cos(\lambda_It+\phi)\sin\phi=0}\\ \displaystyle{\Rightarrow \sin(\lambda_It+\phi-\phi)=0}\\ \displaystyle{\Rightarrow \sin\lambda_It=0}\\ \displaystyle{\Rightarrow \lambda_IT_p=\pi} \end{array}

となって、次式を得ます。

\displaystyle{(18)\quad T_p=\frac{\pi}{\lambda_I}=\frac{\pi}{\omega_n\sqrt{1-\zeta^2}}}

この時刻T_pにおけるS(t)の値は

\displaystyle{(19)\quad S(T_p)=1-\frac{\omega_n}{\lambda_I}e^{\lambda_R\frac{\pi}{\lambda_I}}\sin(\lambda_I\frac{\pi}{\lambda_I}+\phi) =1-\underbrace{\frac{\omega_n}{\lambda_I}}_{\frac{1}{\sin\phi}}e^{\frac{\lambda_R}{\lambda_I}\pi}\underbrace{\sin(\pi+\phi)}_{-\sin\phi}}

となって、次式を得ます。

\displaystyle{(20)\quad S(T_p)=1+e^{\frac{\lambda_R}{\lambda_I}\pi}=1+\underbrace{\exp(-\frac{\zeta\pi}{\sqrt{1-\zeta^2}})}_{p_0} }

●したがって、次のように、(T_p,1+p_0)から(\zeta,\omega_n)を得ることができます。

(21)\quad \begin{array}{l} \displaystyle{\left\{\begin{array}{l} T_p=\frac{\pi}{\omega_n\sqrt{1-\zeta^2}}\\ p_0=\exp(-\frac{\zeta\pi}{\sqrt{1-\zeta^2}}) \end{array}\right.}\\ \displaystyle{\Leftrightarrow \left\{\begin{array}{l} T_p^2=\frac{\pi^2}{\omega_n^2(1-\zeta^2)}\\ \log p_0=-\frac{\zeta\pi}{\sqrt{1-\zeta^2}} \end{array}\right.}\\ \displaystyle{\Leftrightarrow \left\{\begin{array}{l} \omega_n^2=\frac{\pi^2}{T_p^2(1-\zeta^2)}\\ (\log p_0)^2=\frac{\zeta^2\pi^2}{1-\zeta^2} \end{array}\right.}\\ \displaystyle{\Leftrightarrow \left\{\begin{array}{l} \zeta^2=\frac{(\log p_0)^2}{(\log p_0)^2+\pi^2}\\ \omega_n^2=\frac{\pi^2}{T_p^2}\frac{(\log p_0)^2+\pi^2}{\pi^2} \end{array}\right.}\\ \displaystyle{\Leftrightarrow \boxed{\left\{\begin{array}{l} \zeta=\sqrt{\frac{(\log p_0)^2}{(\log p_0)^2+\pi^2}}\\ \omega_n=\frac{\sqrt{(\log p_0)^2+\pi^2}}{T_p} \end{array}\right.}} \end{array}

演習A23…Flipped Classroom

図2のステップ応答#1を描き、第１オーバーシュートの頂点の座標(T_p,1+p_0)をマウスをクリックして取得し、減衰係数と固有角周波数(\zeta,\omega_n)を求めるプログラムを作成せよ。

MATLAB 
%a23.m
%-----
 clear all, close all
 A1=[0 1;-1 -2*0.01]; 
 A2=[0 1;-1 -2/sqrt(2)]; 
 A3=[0 1;-1 -2*1]; 
 B=[0;1]; C=[1 0]; D=[];
 sys1=ss(A1,B,C,D);
 sys2=ss(A2,B,C,D);
 sys3=ss(A3,B,C,D);
 t=0:0.1:10;
 step(sys1,sys2,sys3,t) 
 grid
 title('Step Resp of 2nd-order System')
 xlabel('time')
 ylabel('y(t)')
 legend('zeta=0.01','zeta=0.707','zeta=1')
%-----
 [Tp,p0]=ginput(1);
 p0=p0-1;
 zeta=sqrt(log(p0)^2/(log(p0)^2+pi^2))
 wn=sqrt(log(p0)^2+pi^2)/Tp
 hold on
 plot(Tp,1+p0,'o')
 A1a=[0 1;-wn^2 -2*zeta]; 
 sys1a=ss(A1a,B,C,D);
 step(sys1a,t) 
%-----
%eof
SCILAB

ｎ次系の時間応答

投稿日時: 投稿者:

ｎ次系の時間応答…Homework

[1] 次のm入力p出力n次系の状態空間表現を考えます。

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

これから、n次系の時間応答の表現式

\displaystyle{(2)\quad y(t)=C\exp(At)x(0)+\int_0^tC\exp(A(t-\tau))Bu(\tau)d\tau}

を得ます。ここで、インパルス応答行列

\displaystyle{(3)\quad \boxed{G(t)=C\exp(At)B}}

を定義しますと（G(t)\in{\rm\bf R}^{p\times m}）、上式は

\displaystyle{(4)\quad \boxed{y(t)=\underbrace{C\exp(At)x(0)}_{zero-input\ response} +\underbrace{\int_0^tG(t-\tau)u(\tau)d\tau}_{zero-state\ response}}}

すなわち、n次系の時間応答は、零入力応答（第１項）と零状態応答（第２項）の和となります。

いま、(1)をラプラス変換して、x(0)=0とすると（零状態応答だけを考えます）

\displaystyle{(5a)\quad \hat{x}(s)-\underbrace{x(0)}_0=A\hat{x}(s)+B\hat{u}(s) \Rightarrow \hat{x}(s)=(sI_n-A)^{-1}B\hat{u}(s) }

\displaystyle{(5b)\quad \hat{y}(s)=C\hat{x}(s) \Rightarrow \hat{y}(s)=\underbrace{C(sI_n-A)^{-1}B}_{\hat{G}(s)}\hat{u}(s) }

ここで現れた

\displaystyle{(6)\quad \boxed{\hat{G}(s)=C(sI_n-A)^{-1}B}}

を伝達関数行列と呼びます。

●n次系の時間応答についての相互関係は次図のように表されます。多変数系の場合、i番目の入力からj番目の出力までのステップ応答と周波数応答を考えます。

図1 多変数系時間等応答の相互関係

２次系の時間応答…Homework

[2] さて、次のようにパラメタライズされた漸近安定な１入力１出力２次系の状態空間表現を考えます。\zeta>0、\omega_n、K はそれぞれ減衰係数、固有角周波数、定常ケインと呼ばれます。

\displaystyle{(7a)\quad \boxed{\underbrace{ \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{x}_2(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{cc} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{cc} 0 \\ K\omega_n^2 \end{array}\right] }_{B} u(t)} }

\displaystyle{(7b)\quad y(t)= \underbrace{ \left[\begin{array}{cc} 1 & 0 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] }_{x(t)} }

以下では、簡単のためK=1とした２次系に対するインパルス応答(3)を求めます。そのために、行列Aの固有値\lambda(A)=\{\lambda_1,\lambda_2\}を計算すると、\zeta>1、\zeta<1、\zeta=1に応じて、それぞれ次に示すようになります。

\displaystyle{(8)\quad \boxed{\left\{\begin{array}{l} \lambda_1,\lambda_2=-\zeta\omega_n\pm\omega_n\sqrt{\zeta^2-1}\quad(\zeta>1)\\ \lambda_1,\lambda_2=\underbrace{-\zeta\omega_n}_{\lambda_R}\pm j\underbrace{\omega_n\sqrt{1-\zeta^2}}_{\lambda_I}\quad(\zeta<1)\\ \lambda_1=\lambda_2=\lambda=-\zeta\omega_n=-\omega_n\quad(\zeta=1) \end{array}\right.} }

図2 A=\left[\begin{array}{cc} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right]の固有値分布

[3] ２次系のインパルス応答(1)

\zeta>1のとき、\lambda_1,\lambda_2=-\zeta\omega_n\pm\omega_n\sqrt{\zeta^2-1}とおくと、インパルス応答は次式で与えられます。

\displaystyle{(9)\quad \boxed{G(t)=\frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}(e^{\lambda_2t}-e^{\lambda_1t})}}

実際

\displaystyle{(10a)\quad \underbrace{ \left[\begin{array}{cc} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] }_{A} = \underbrace{ \left[\begin{array}{cc} 1 & 1\\ \lambda_1 & \lambda_2 \end{array}\right] }_{V} \underbrace{ \left[\begin{array}{cc} \lambda_1& 0\\ 0 & \lambda_2 \end{array}\right] }_{\Lambda} \underbrace{ \left[\begin{array}{cc} 1 & 1\\ \lambda_1 & \lambda_2 \end{array}\right]^{-1} }_{V^{-1}} }

\displaystyle{(10b)\quad \exp(\Lambda t)= \left[\begin{array}{cc} e^{\lambda_1t}& 0\\ 0 & e^{\lambda_2t} \end{array}\right] }

となることから、インパルス応答は次のように計算されます。

(11)\quad \begin{array}{l} \displaystyle{G(t)=C\exp(At)B=CV\exp(\Lambda t)V^{-1}B}\\ \displaystyle{ =\left[\begin{array}{cc} 1 & 0 \end{array}\right] \left[\begin{array}{cc} 1 & 1\\ \lambda_1 & \lambda_2 \end{array}\right] \left[\begin{array}{cc} e^{\lambda_1t}& 0\\ 0 & e^{\lambda_2t} \end{array}\right] \frac{1}{\lambda_2-\lambda_1} \left[\begin{array}{cc} \lambda_2 & -1\\ -\lambda_1 & 1 \end{array}\right] \left[\begin{array}{cc} 0 \\ \omega_n^2 \end{array}\right]}\\ \displaystyle{ =\frac{\omega_n^2}{\lambda_2-\lambda_1} \left[\begin{array}{cc} e^{\lambda_1t}& e^{\lambda_2t} \end{array}\right] \left[\begin{array}{cc} -1\\ 1 \end{array}\right]}\\ \displaystyle{=\frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}(e^{\lambda_2t}-e^{\lambda_1t})} \end{array}

[4] ２次系のインパルス応答(2)

\zeta<1のとき、\lambda_1,\lambda_2=\underbrace{-\zeta\omega_n}_{\lambda_R}\pm j\underbrace{\omega_n\sqrt{1-\zeta^2}}_{\lambda_I}とおくと、インパルス応答とステップ応答はそれぞれそれぞれ次式で与えられます。

\displaystyle{(12)\quad \boxed{G(t)=\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\sin\lambda_I t}}

実際

\displaystyle{(13a)\quad \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] }_{A} = \underbrace{ \left[\begin{array}{cc} 1 & 0 \\ \lambda_R & \lambda_I \end{array}\right] }_{V} \underbrace{ \left[\begin{array}{cc} \lambda_R & \lambda_I \\ -\lambda_I & \lambda_R \end{array}\right] }_{\Lambda} \underbrace{ \frac{1}{\lambda_I} \left[\begin{array}{cc} \lambda_I & 0 \\ -\lambda_R & 1 \end{array}\right] }_{V^{-1}}}
\displaystyle{(13b)\quad \exp(\Lambda t)= e^{\lambda_R t} \left[\begin{array}{cc} \cos\lambda_I t & \sin\lambda_I t\\ \sin\lambda_I t & \cos\lambda_I t \end{array}\right]}

となることから、インパルス応答は次のように計算されます。

(14)\quad \begin{array}{l} \displaystyle{G(t)=C\exp(At)B=CV\exp(\Lambda t)V^{-1}B}\\ \displaystyle{=\left[\begin{array}{cc} 1 & 0 \end{array}\right] \left[\begin{array}{cc} 1 & 0\\ \lambda_R & \lambda_I \end{array}\right] e^{\lambda_R t} \left[\begin{array}{cc} \cos\lambda_I t & \sin\lambda_I t\\ \sin\lambda_I t & \cos\lambda_I t \end{array}\right] \frac{1}{\lambda_I} \left[\begin{array}{cc} \lambda_I & 0\\ -\lambda_R & 1 \end{array}\right] \left[\begin{array}{cc} 0 \\ \omega_n^2 \end{array}\right]}\\ \displaystyle{=\frac{\omega_n^2}{\lambda_I} e^{\lambda_R t} \left[\begin{array}{cc} 1 & 0 \end{array}\right] \left[\begin{array}{cc} \cos\lambda_I t & \sin\lambda_I t\\ \sin\lambda_I t & \cos\lambda_I t \end{array}\right] \left[\begin{array}{cc} 0\\ 1 \end{array}\right]}\\ \displaystyle{=\frac{\omega_n^2}{\lambda_I} e^{\lambda_R t}\sin\lambda_I t} \end{array}

[5] ２次系のインパルス応答(3)

\zeta=1のとき、\lambda_1=\lambda_2=\lambda=-\zeta\omega_n=-\omega_nとおくと、インパルス応答とステップ応答はそれぞれそれぞれ次式で与えられます。

\displaystyle{(15)\quad \boxed{G(t)=\lambda^2te^{\lambda t}}}

実際

\displaystyle{(16a)\quad \underbrace{ \left[\begin{array}{cc} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] }_{A} = \underbrace{ \left[\begin{array}{cc} 1 & 1\\ \lambda & \lambda+1 \end{array}\right] }_{V} \underbrace{ \left[\begin{array}{cc} \lambda & 1\\ 0 & \lambda \end{array}\right] }_{\Lambda} \underbrace{ \left[\begin{array}{cc} 1 & 1\\ \lambda & \lambda+1 \end{array}\right]^{-1} }_{V^{-1}}}
\displaystyle{(16b)\quad \exp(\Lambda t)= e^{\lambda t} \left[\begin{array}{cc} 1 & t\\ 0 & 1 \end{array}\right] }

となることから、インパルス応答は次のように計算されます。

(17)\quad \begin{array}{l} \displaystyle{G(t)=C\exp(At)B=CV\exp(\Lambda t)V^{-1}B}\\ \displaystyle{= \left[\begin{array}{cc} 1 & 0 \end{array}\right] \left[\begin{array}{cc} 1 & 1\\ \lambda & \lambda+1 \end{array}\right] e^{\lambda t} \left[\begin{array}{cc} 1 & t\\ 0 & 1 \end{array}\right] \frac{1}{\lambda+1-\lambda} \left[\begin{array}{cc} \lambda+1 & -1\\ -\lambda & 1 \end{array}\right] \left[\begin{array}{cc} 0 \\ \omega_n^2 \end{array}\right]}\\ \displaystyle{= \omega_n^2 e^{\lambda t} \left[\begin{array}{cc} 1 & t+1 \end{array}\right] \left[\begin{array}{cc} -1\\ 1 \end{array}\right]}\\ \displaystyle{=\lambda^2te^{\lambda t}} \end{array}

２次振動系の同定…Homework

[7] ２次系(7)において、\zeta=0.01,1/\sqrt{2},1、\omega_n=1、K=1の場合、インパルス応答のシミュレーション例を示します。

図3 ２次系のインパルス応答の比較

\zeta<1の場合の(7)で表される２次系を２次振動系と呼びます。そのインパルス応答

\displaystyle{(18)\quad G(t)=\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\sin\lambda_I t}

の第k+1番目(k=0,1,\cdots)の頂点の座標について、次が成り立ちます。

\displaystyle{(19)\quad \boxed{(t_k,G(t_k))=\left(\frac{k\pi+\phi}{\omega_n\sqrt{1-\zeta^2}},\omega_n\exp(-\frac{\zeta(k\pi+\phi)}{\sqrt{1-\zeta^2}})\right)} }

実際、\dot{G}(t)=0を計算すると

(20)\quad \begin{array}{l} \displaystyle{\dot{G}(t)=\frac{\omega_n^2}{\lambda_I}(\lambda_Re^{\lambda_Rt}\sin\lambda_It+e^{\lambda_Rt}\cos\lambda_It\times\lambda_I)=0}\\ \displaystyle{\Rightarrow \sin\lambda_It\times\frac{-\lambda_R}{\omega_n}-\cos\lambda_It\times\frac{\lambda_I}{\omega_n}=0}\\ \displaystyle{\Rightarrow \sin\lambda_It\cos\phi-\cos\lambda_It\sin\phi=0,\quad\phi=\tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta}}\\ \displaystyle{\Rightarrow \sin(\lambda_It-\phi)=0}\\ \displaystyle{\Rightarrow \lambda_It_k-\phi=k\pi\quad(k=0,1,\cdots)} \end{array}

となって、次式を得ます。

\displaystyle{(21)\quad t_k=\frac{k\pi+\phi}{\lambda_I}=\frac{k\pi+\phi}{\omega_n\sqrt{1-\zeta^2}}

この時刻t_kにおけるG(t)の値は

\displaystyle{(22)\quad \begin{array}{l} G(t_k)=\frac{\omega_n^2}{\lambda_I}e^{\lambda_R\frac{k\pi+\phi}{\lambda_I}}\sin(\lambda_I\frac{\phi}{\lambda_I})\\ =\omega_n\underbrace{\frac{\omega_n}{\lambda_I}}_{\frac{1}{\sin\phi}}e^{\frac{\lambda_R}{\lambda_I}(k\pi+\phi)}\underbrace{\sin(k\pi+\phi)}_{(-1)^k\sin\phi}\\=(-1)^k\omega_ne^{\frac{\lambda_R}{\lambda_I}(k\pi+\phi)} \end{array} }

となって、次式を得ます。

\displaystyle{(23)\quad G(t_k)=(-1)^k\omega_ne^{\frac{\lambda_R}{\lambda_I}(k\pi+\phi)} =(-1)^k\omega_n\exp(-\frac{\zeta}{\sqrt{1-\zeta^2}}(k\pi+\phi)) }}

●さて、(19)からそれぞれ次の関係式を得ます。

\displaystyle{(24)\quad t_{k+2}-t_k=\frac{(k+2)\pi+\phi}{\omega_n\sqrt{1-\zeta^2}}-\frac{k\pi+\phi}{\omega_n\sqrt{1-\zeta^2}} =\frac{2\pi}{\omega_n\sqrt{1-\zeta^2}}=T}

\displaystyle{(25)\quad \frac{G(t_{k+2})}{G(t_k)} =\exp(-\frac{\zeta((k+2)\pi+\phi)}{\sqrt{1-\zeta^2}} +\frac{\zeta(k\pi+\phi)}{\sqrt{1-\zeta^2}}) =\exp(-\frac{2\pi\zeta}{\sqrt{1-\zeta^2}})=\Lambda}

(24)はインパルス応答の隣合う山と山の時間間隔（または谷と谷の時間間隔）Tは一定であることを示しています。(25)はインパルス応答の隣合う山と山の振幅比（または谷と谷の振幅比）\Lambdaは一定であることを示しています。

したがって、次のように、(T,\Lambda)から(\zeta,\omega_n)を得ることができます。

(26)\quad \begin{array}{l} \displaystyle{ \left\{\begin{array}{l} T=\frac{2\pi}{\omega_n\sqrt{1-\zeta^2}}\\ \Lambda=\exp(-\frac{2\pi\zeta}{\sqrt{1-\zeta^2}}) \end{array}\right.}\\ \displaystyle{\Leftrightarrow \left\{\begin{array}{l} T^2=\frac{4\pi^2}{\omega_n^2(1-\zeta^2)}\\ \log\Lambda=-\frac{2\pi\zeta}{\sqrt{1-\zeta^2}} \end{array}\right.}\\ \displaystyle{\Leftrightarrow \left\{\begin{array}{l} \omega_n^2=\frac{4\pi^2}{T^2(1-\zeta^2)}\\ (\log\Lambda)^2=\frac{4\pi^2\zeta^2}{1-\zeta^2} \end{array}\right.}\\ \displaystyle{\Leftrightarrow \left\{\begin{array}{l} \zeta^2=\frac{(\log\Lambda)^2}{(\log\Lambda)^2+4\pi^2}\\ \omega_n^2=\frac{4\pi^2}{T^2}\frac{(\log\Lambda)^2+4\pi^2}{4\pi^2} \end{array}\right.}\\ \displaystyle{ \Leftrightarrow \boxed{\left\{\begin{array}{l} \zeta=\sqrt{\frac{(\log\Lambda)^2}{(\log\Lambda)^2+4\pi^2}}\\ \omega_n=\frac{\sqrt{(\log p_0)^2+4\pi^2}}{T} \end{array}\right.}} \end{array}

演習A22…Flipped Classroom

図2の応答#1を描き、２つの山の頂点の座標をマウスをクリックして取得し、減衰係数と固有角周波数を求めるプログラムを作成せよ。

MATLAB 
%a22.m
%-----
 clear all, close all
 A1=[0 1;-1 -2*0.01]; 
 A2=[0 1;-1 -2/sqrt(2)]; 
 A3=[0 1;-1 -2*1]; 
 B=[0;1]; C=[1 0]; D=[];
 sys1=ss(A1,B,C,D);
 sys2=ss(A2,B,C,D);
 sys3=ss(A3,B,C,D);
 t=0:0.1:10;
 x0=B;
 initial(sys1,sys2,sys3,x0,t) 
 grid
 title('Impulse Resp of 2nd-order System')
 xlabel('time')
 ylabel('y(t)')
 legend('zeta=0.01','zeta=0.707','zeta=1')
%-----
 [Tp,p0]=ginput(2);
 T=Tp(2)-Tp(1);
 Lambda=p0(2)/p0(1),
 zeta=sqrt(log(Lambda)^2/(log(Lambda)^2+4*pi^2))
 wn=sqrt(log(Lambda)^2+4*pi^2)/T
 hold on
 plot(Tp(1),p0(1),'o',Tp(2),p0(2),'o')
 A1a=[0 1;-wn^2 -2*zeta]; 
 sys1a=ss(A1a,B,C,D);
 initial(sys1a,x0,t) 
%-----
%eof
SCILAB

Note A22-1　状態方程式の解

n次系の状態方程式

\displaystyle{(1)\quad\dot{x}(t)=Ax(t)+Bu(t)}

の解を求めてみます。Note A21と同様にして

\displaystyle{(2)\quad\frac{d}{dt}(\exp(-At)x(t))=\exp(-At)Bu(t)}

までは問題ないと思います。これを0からt'まで積分して

\displaystyle{(3)\quad\left[\exp(-At)x(t)\right]_0^{t'}=\int_0^{t'}\exp(-At)Bu(t)dt}

すなわち

\displaystyle{(4)\quad\exp(-At')x(t')-x(0)=\int_0^{t'}\exp(-At)Bu(t)dt}

したがって

\displaystyle{(5)\quad x(t')=\exp(At')x(0)+\exp(At')\int_0^{t'}\exp(-At)Bu(t)dt}

ここで、tを\tauに、t'をtと置き換えれば

\displaystyle{(6)\quad\boxed{x(t)=\exp(At)x(0)+\int_0^t\exp(A(t-\tau))Bu(\tau)d\tau}}

が得られます。

Note A22-2　重ね合わせの原理

次のm入力p出力n次系の状態空間表現を考えます。

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

２つの入力u_1とu_2に対するそれぞれの零状態応答は

\displaystyle{(2)\quad y_1(t)=\int_0^tC\exp(A(t-\tau))Bu_1(\tau)d\tau}

\displaystyle{(3)\quad y_2(t)=\int_0^tC\exp(A(t-\tau))Bu_2(\tau)d\tau}

となります。このとき入力

\displaystyle{(4)\quad \alpha u_1(t)+\beta u_2(t)}

に対する零状態応答は

\displaystyle{(5)\quad \alpha y_1(t)+\beta y_2(t)}

となります。実際、

(6)\quad \begin{array}{l} \displaystyle{\int_0^tC\exp(A(t-\tau))B(\alpha u_1(\tau)+\beta u_2(\tau))d\tau}\\ \displaystyle{=\alpha \int_0^tC\exp(A(t-\tau))Bu_1(\tau)d\tau +\beta \int_0^tC\exp(A(t-\tau))Bu_2(\tau))d\tau}\\ =\alpha y_1(t)+\beta y_2(t)} \end{array}

このとき、重ね合わせの原理が成り立つ、n次系(1)は線形系(Linear System)であると言います。重ね合わせの原理が成り立たないときは非線形系(Nonlinear System)であると言います。たとえば入力を２倍したときの時間応答が元の時間応答の２倍にならなければ非線形系であると判断できます。

漸近安定性

投稿日時: 投稿者:

漸近安定性…Homework

[0]　次図に示すように軸支された剛体振り子は２つの静止した平衡状態(\alpha)と(\beta)を持ちます。

図1 剛体振り子

いまこの平衡状態から棒を少し傾けて手を離すと、明らかに違った振舞いをします。実際には軸回りの摩擦や空気抵抗による抗力があるので、平衡状態(\alpha)の方は元に戻りますが、平衡状態(\beta)の方は元に戻ることなく平衡状態(\alpha)に落ち着くでしょう。前者を漸近安定、後者を不安定とよびます。このような違いは平衡状態回りでの振舞いを表す線形状態方程式にどのように反映されているのか気になるところです。

●この場合の運動方程式は次式となります。

\displaystyle{(1)\quad J\ddot{\theta}(t)=-mg\ell\sin\theta(t)-c\dot{\theta}(t)\quad(J=\frac{4}{3}m\ell^2) }

ここで、右辺の第２項が角速度に比例する抗力を表しています（c>0は定数）。これから平衡状態(\alpha)と(\beta)に応じて次の状態方程式を得ます。

\displaystyle{(2)\quad \underbrace{ \left[\begin{array}{l} \dot{\theta}(t) \\ \dot{\omega}(t) \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -\frac{3g}{4\ell}&-\frac{c}{J} \end{array}\right] }_{A=A_1} \underbrace{ \left[\begin{array}{l} \theta(t) \\ \omega(t) \end{array}\right] }_{x(t)} }

\displaystyle{(3)\quad \underbrace{ \frac{d}{dt} \left[\begin{array}{l} \theta(t)-\pi \\ \omega(t)-0 \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ \frac{3g}{4\ell}&-\frac{c}{J} \end{array}\right] }_{A=A_2} \underbrace{ \left[\begin{array}{l} \theta(t)-\pi \\ \omega(t)-0 \end{array}\right] }_{x(t)} }

以下では、行列Aの固有値を調べて、漸近安定かどうかを判定できることを説明します。

[1]　制御対象が平衡状態にあることは線形状態方程式\dot{x}(t)=Ax(t)+Bu(t)において、x=0, u=0を意味します。そこで平衡状態が乱されてx\ne0となる時刻をt=0にとると、線形状態方程式は次式となります。

\displaystyle{(4)\quad \dot{x}(t)=Ax(t)\qquad(x(0)\ne0) }

漸近安定であることは平衡状態x=0に戻ることでしたから、次のように表すことができます。

\displaystyle{(5)\quad \boxed{x(t)\rightarrow 0\qquad(t\rightarrow\infty)} }

そのための条件を検討するため、微分方程式\dot{x}(t)=Ax(t)の解を求める必要があります。これは次式で表されます（Note A21-1参照）。

\displaystyle{(6)\quad \boxed{x(t)=\exp(At)x(0)} }

ここで、xはn次元ベクトル、Aはn次正方行列です。任意のn次正方行列Xに対して、行列指数関数\exp(X)は

\displaystyle{(7)\quad \boxed{\exp(X)=I_n+X+\frac{1}{2}X^2+\cdots+\frac{1}{k!}X^k+\cdots} }

で定義されます。そこで、解の表現式(6)をn=2の場合について詳しく説明します。

[2]　いま２次系に限定して、微分方程式\dot{x}(t)=Ax(t)を要素を明示して書くと

\displaystyle{(8)\quad \underbrace{ \left[\begin{array}{l} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} a_{11} & a_{12} \\ a_{21} & a_{22} \\ \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{l} x_1(t) \\ x_2(t) \end{array}\right] }_{x(t)} }

すなわち

\displaystyle{(9)\quad \left\{\begin{array}{l} \dot{x}_1(t)=a_{11}x_1(t)+a_{12}x_2(t) \\ \dot{x}_2(t)=a_{21}x_1(t)+a_{22}x_2(t) \end{array}\right. }

となります。もしAが対角行列で、a_{12}=a_{21}=0であれば

\displaystyle{(10)\quad \left\{\begin{array}{l} \dot{x}_1(t)=a_{11}x_1(t) \\ \dot{x}_2(t)=a_{22}x_2(t) \end{array}\right. }

となって、解は

\displaystyle{(11)\quad \left\{\begin{array}{l} x_1(t)=e^{a_{11}t}x_1(0) \\ x_2(t)=e^{a_{22}t}x_2(0) \\ \end{array}\right. }

すなわち

\displaystyle{(12)\quad \underbrace{ \left[\begin{array}{l} {x}_1(t) \\ {x}_2(t) \end{array}\right] }_{{x}(t)}= \underbrace{ \left[\begin{array}{cc} e^{a_{11}t} & 0 \\ 0 & e^{a_{22}t} \\ \end{array}\right] }_{\exp(At)} \underbrace{ \left[\begin{array}{l} x_1(0) \\ x_2(0) \end{array}\right] }_{x(0)} }

ですから、解の表現式(6)は納得できると思います。

[3]　それではAが対角行列でない場合はどうするかですが、一対一対応の変数変換

\displaystyle{(13)\quad \left\{\begin{array}{l} {y}_1(t)=t_{11}x_1(t)+t_{12}x_2(t) \\ {y}_2(t)=t_{21}x_1(t)+t_{22}x_2(t) \end{array}\right. }

すなわち

\displaystyle{(14)\quad \underbrace{ \left[\begin{array}{l} {y}_1(t) \\ {y}_2(t) \end{array}\right] }_{y(t)}= \underbrace{ \left[\begin{array}{cc} t_{11} & t_{12} \\ t_{21} & t_{22} \\ \end{array}\right] }_{T} \underbrace{ \left[\begin{array}{l} x_1(t) \\ x_2(t) \end{array}\right] }_{x(t)} }

を考えます。またこの逆変換を次式で表しておきます。

\displaystyle{(15)\quad \left\{\begin{array}{l} {x}_1(t)=v_{11}y_1(t)+v_{12}y_2(t) \\ {x}_2(t)=v_{21}y_1(t)+v_{22}y_2(t) \end{array}\right. }

すなわち

\displaystyle{(16)\quad \underbrace{ \left[\begin{array}{l} {x}_1(t) \\ {x}_2(t) \end{array}\right] }_{x(t)}= \underbrace{ \left[\begin{array}{cc} v_{11} & v_{12} \\ v_{21} & v_{22} \\ \end{array}\right] }_{V} \underbrace{ \left[\begin{array}{l} y_1(t) \\ y_2(t) \end{array}\right] }_{y(t)} }

ここで、TV=VT=I_2が成り立つので、VはTの逆行列T^{-1}（V=T^{-1}）となっています。逆行列をもつ行列は正則行列と呼ばれます。

すなわち、微分方程式\dot{x}(t)=Ax(t)に対して、正則行列Tによる変数変換y(t)=Tx(t)を行ないます。x(t)=Vy(t)（V=T^{-1}）と両辺を微分した\dot{x}(t)=V\dot{y}(t)を代入して次を得ます。

\displaystyle{(17)\quad V\dot{y}(t)=AVy(t)\Rightarrow\dot{y}(t)=\underbrace{V^{-1}AV}_{\Lambda}y(t) }

この変換された行列\Lambdaをできるだけ対角化します。その議論は、線形代数において\Lambda=V^{-1}AVを標準形とする方法として学びました。もし\exp(\Lambda t)を求めることができると

\displaystyle{(18)\quad x(t)=Vy(t)=V\exp(\Lambda t)y(0)=V\exp(\Lambda t)V^{-1}x(0) }

となって、漸近安定性(5)の判定は、

\displaystyle{(19)\quad \boxed{\exp(\Lambda t)\rightarrow 0\qquad(t\rightarrow\infty)} }

の判定に帰着されます。ここで、収束先の0はn次の零行列です。

[4]　さて、２次行列の標準形は次の３つに分類されます。

\displaystyle{(20)\quad \Lambda_1= \left[\begin{array}{cc} \lambda_1 & 0 \\ 0 & \lambda_2 \end{array}\right] \quad(\lambda(A)=\{\lambda_1,\lambda_2\}) }

\displaystyle{(21)\quad \Lambda_2= \left[\begin{array}{cc} \lambda_R & \lambda_I \\ -\lambda_I & \lambda_R \end{array}\right] \quad(\lambda(A)=\{\lambda_R\pm j\lambda_I\}) }

\displaystyle{(22)\quad \Lambda_3= \left[\begin{array}{cc} \lambda & 1 \\ 0 & \lambda \end{array}\right] \quad(\lambda(A)=\{\lambda,\lambda\}) }

ここで、jは虚数単位\sqrt{-1}、\lambda_1,\lambda_2、\lambda_R,\lambda_I、\lambdaはすべて実数です。また任意のn次正方行列Xに対する行列式を{\rm det} Xで表すと

\displaystyle{(23)\quad \lambda(A)=\{\lambda:{\rm det}(\lambda I_n-A)=0\} }

すなわち、\lambda(A)は行列Aの固有値の集合です。{\rm det}(\lambda I_n-A)は特性多項式と呼ばれます。

上の３種類の\Lambdaについて\exp(\Lambda t)を具体的に計算してみます。その結果は次式となります。

\displaystyle{(24)\quad \exp(\Lambda_1 t)= \left[\begin{array}{cc} e^{\lambda_1t} & 0 \\ 0 & e^{\lambda_2 t} \end{array}\right] }

\displaystyle{(25)\quad \exp(\Lambda_2 t)=e^{\lambda_R t} \left[\begin{array}{cc} \cos(\lambda_It) & \sin(\lambda_It) \\ -\sin(\lambda_It) & \cos(\lambda_It) \end{array}\right] }

\displaystyle{(26)\quad \exp(\Lambda_3 t)=e^{\lambda t} \left[\begin{array}{cc} 1 & t \\ 0 & 1 \end{array}\right] }

したがって、漸近安定性(5)すなわち(19)が成り立つのは、(24)のとき\lambda_1<0かつ\lambda_2<0、(25)のとき\lambda_R<0、(26)のとき\lambda<0であることがわかると思います。したがって、

漸近安定性のための条件は行列Aのすべての固有値の実部が負であること

です。対偶をとれば、一つでも実部が負の固有値があれば不安定と判定します。実部が負の固有値を安定固有値、実部が非負の固有値を不安定固有値と呼びます。

●この漸近安定性の条件は、高次系についても同様に成り立ちます（Note A21-2を参照）。

●上の剛体振り子の場合の漸近安定性を判定するために行列Aの固有値を調べてみます。まず、平衡状態(\alpha)の場合、次式を得ます。

(27)\quad \begin{array}{l} \displaystyle{{\rm det}(\lambda I_2-\left[\begin{array}{cc} 0 & 1 \\ -\frac{3g}{4\ell}&-\frac{c}{J} \end{array}\right]) =\lambda^2+\frac{c}{J}\lambda+\frac{3g}{4\ell}=0}\\ \displaystyle{\Rightarrow \lambda=\frac{1}{2}(-\frac{c}{J}\pm\sqrt{\frac{c^2}{J^2}-4\frac{3g}{4\ell}})} \end{array} }

これから行列Aの固有値の実部は必ず負となります（根号内が負の場合は明らか、正の場合は根号を開いても\frac{c}{J}より小）。したがって、平衡状態(\alpha)は漸近安定と判定できます。

次に、平衡状態(\beta)の場合、次式を得ます。

(28)\quad \begin{array}{l} \displaystyle{{\rm det}(\lambda I_2-\left[\begin{array}{cc} 0 & 1 \\ \frac{3g}{4\ell}&-\frac{c}{J} \end{array}\right]) =\lambda^2+\frac{c}{J}\lambda-\frac{3g}{4\ell}=0}\\ \displaystyle{\Rightarrow \lambda=\frac{1}{2}(-\frac{c}{J}\pm\sqrt{\frac{c^2}{J^2}+4\frac{3g}{4\ell}})} \end{array} }

これから行列Aの固有値は実数となり、一つは正、他は負となります。したがって、平衡状態(\beta)は漸近安定ではないと判定できます。

実は２次系の場合特性多項式の係数がすべて正であることが、漸近安定であるための条件であることがラウスフルビッツの判定法として知られています。

演習A21-1…Flipped Classroom

1^\circ 行列指数関数について次の指数法則が成り立つことを示せ。

\displaystyle{(29)\quad XY=YX\Rightarrow\exp(X+Y)=\exp(X)\exp(Y) }

2^\circ 行列指数関数の定義から、(24)を示せ。

3^\circ 行列指数関数の定義から、次式を用いて、(25)を示せ。

\displaystyle{(30)\quad \left[\begin{array}{cc} \lambda_R & \lambda_I \\ -\lambda_I & \lambda_R \end{array}\right]= \underbrace{\lambda_RI_2}_{X}+ \underbrace{\lambda_I \left[\begin{array}{cc} 0 & 1 \\ -1 & 0 \end{array}\right]}_{Y} }

4^\circ 行列指数関数の定義から、次式を用いて、(26)を示せ。

\displaystyle{(31)\quad \left[\begin{array}{cc} \lambda & 1 \\ 0 & \lambda \end{array}\right]= \underbrace{\lambda I_2}_{X}+ \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right]}_{Y} }

演習A21-2…Flipped Classroom

与えられた行列Aに基づいて漸近安定性を調べるプログラムを次に示す。行列Pを実部が大きな順に表示されるように改良せよ。

MATLAB 
%a21.m
%-----
 clear all, close all
 A=rand(10);
 p=eig(A)
%-----
 sys=ss(A,[],[],[]);
 damp(sys)
%----- 
 P=[real(p) imag(p) -cos(angle(p)) abs(p) -1./real(p)]
%-----
%eof
SCILAB 
coming soon

Note A21-1　自由系の時間応答

n次自由系を表す微分方程式

\displaystyle{(1)\quad\dot{x}(t)=Ax(t)}

の解は

\displaystyle{(2)\quad x(t)=\exp(At)x(0)}

と表されます。これが解であることは元の微分方程式に代入すればすぐに確かめられ、また次のようにして導くことができます。

元の微分方程式を

\displaystyle{(3)\quad\dot{x}(t)-Ax(t)=0}

と書いて、左から積分因数と呼ばれる\exp(-At)をかけると

\displaystyle{(4)\quad\exp(-At)\dot{x}(t)-\exp(-At)Ax(t)=0}

すなわち

\displaystyle{(5)\quad\frac{d}{dt}(\exp(-At)x(t))=0}

したがって、cを定数ベクトルとして

\displaystyle{(6)\quad\exp(-At)x(t)=c\ \Rightarrow\ x(t)=\exp(At)c}

ここで、t=0とおくとcは初期値x(0)に等しいので

\displaystyle{(7)\quad x(t)=\exp(At)x(0)}

と表されます。

Note A21-2　高次系の漸近安定性

●行列指数関数を用いると、微分方程式

\displaystyle{(1)\quad \dot{x}(t)=Ax(t)}

の解は次のように表されます。

\displaystyle{(2)\quad x(t)=\exp(At)x(0)=V\exp(\Lambda t)V^{-1}x(0)}

ここで

\displaystyle{(3)\quad \begin{array}{ll} &\exp(At)=V\exp(\Lambda t)V^{-1}={\displaystyle \sum_{i=1}^p\sum_{k=1}^{m_i}t^{k-1}}e^{\lambda_it}\,R_{ik} \\ &+{\displaystyle \sum_{i=1}^q\sum_{k=1}^{\ell_i}t^{k-1}}e^{\lambda_{Ri}t}\cos\lambda_{Ii}t\,C_{ik} +{\displaystyle \sum_{i=1}^q\sum_{k=1}^{\ell_i}t^{k-1}}e^{\lambda_{Ri}t}\sin\lambda_{Ii}t\,S_{ik} \end{array}

と書けることに注意します（R_{ik},C_{ik},S_{ik}は適当なn次実正方行列）。

したがって、任意のx(0)\ne0に対して、t\rightarrow\inftyのときx(t)\rightarrow 0となるための条件は

\displaystyle{(4)\quad \exp(\Lambda t)\rightarrow 0\quad(t\rightarrow\infty)}

となります。これはAのすべての固有値の実部が負を意味します。

●\dot{x}(t)=ax(t)（a=1,0,0.5）の解のグラフを見ると、a=0の場合は、漸近安定ではないが、発散はしないので、不安定とまではいえないのではないかと思うかもしれません。したがって零の固有値を不安定とみなすのか、安定とみなすか迷うところです。しかし、\dot{x}(t)=Ax(t)において、A=\left[\begin{array}{cc} 0& 1\\ 0& 0 \end{array}\right]の場合、解はx(t)=\left[\begin{array}{cc} 1& t\\ 0& 1 \end{array}\right]x(0)となって、x(0)の第２要素が零でない場合は発散します。

いま、\dot{x}(t)=Ax(t)において、Aのすべての固有値の実部は負または零の場合を考え、実部が零の固有値\lambdaの重複度をq\ge2とし、次が成り立つとします。

\displaystyle{(5)\quad n-{\rm rank}(A-\lambda I_n)=q}

左辺は幾何学的重複度、左辺は代数学的重複度と呼ばれています。すべての実部が零の固有値について、幾何学的重複度と代数学的重複度が等しいことが、解が発散しないために追加的に要請される条件であることが知られています。上の例の場合、幾何学的重複度は2-1=1、代数学的重複度は2ですから、この条件は成立していないことがわかります。

Note A21-3　自由系の時間応答は固有値と固有ベクトルからどのように構成されるか

●２次自由系

\displaystyle{(1)\quad \dot{x}(t)=Ax(t)}

の解は次のように表されます。

\displaystyle{(2)\quad x(t)=\exp(At)x(0)}

２次行列Aの標準形は次の３つに分類されます。

\displaystyle{(3a)\quad \Lambda_1= \left[\begin{array}{cc} \lambda_1 & 0 \\ 0 & \lambda_2 \end{array}\right] \quad(\lambda(A)=\{\lambda_1,\lambda_2\}) }

\displaystyle{(3b)\quad \Lambda_2= \left[\begin{array}{cc} \lambda_R & \lambda_I \\ -\lambda_I & \lambda_R \end{array}\right] \quad(\lambda(A)=\{\lambda_R\pm j\lambda_I\}) }

\displaystyle{(3c)\quad \Lambda_3= \left[\begin{array}{cc} \lambda & 1 \\ 0 & \lambda \end{array}\right] \quad(\lambda(A)=\{\lambda,\lambda\}) }

ここで、jは虚数単位\sqrt{-1}、\lambda_1,\lambda_2、\lambda_R,\lambda_I、\lambdaはすべて実数です。

このとき２次行列は、対応する固有ベクトルを用いて、次のように表されます。

\displaystyle{(4a)\quad A=\underbrace{\left[\begin{array}{cc} v_1 & v_2 \end{array}\right]}_{V} \underbrace{\left[\begin{array}{cc} \lambda_1 & 0 \\ 0 & \lambda_2 \end{array}\right]}_{\Lambda_1} \underbrace{\left[\begin{array}{cc} v_1 & v_2 \end{array}\right]^{-1}}_{V^{-1}} }

\displaystyle{(4b)\quad \begin{array}{l} A=\underbrace{\left[\begin{array}{cc} v_R+jv_I & v_R-jv_I \end{array}\right]}_{V} \left[\begin{array}{cc} \lambda_R+j\lambda_I & 0 \\ 0 & \lambda_R-j\lambda_I \end{array}\right]\\ \times\underbrace{\left[\begin{array}{cc} v_R+jv_I & v_R-jv_I \end{array}\right]^{-1}}_{V^{-1}}\\ =\underbrace{\left[\begin{array}{cc} v_R & v_I \end{array}\right]}_{V'} \underbrace{\left[\begin{array}{cc} \lambda_R & \lambda_I \\ -\lambda_I & \lambda_R \end{array}\right]}_{\Lambda_2} \underbrace{\left[\begin{array}{cc} v_R & v_I \end{array}\right]^{-1}}_{V'^{-1}} \end{array} }

\displaystyle{(4c)\quad A=\underbrace{\left[\begin{array}{cc} v & v' \end{array}\right]}_{V} \underbrace{\left[\begin{array}{cc} \lambda & 1 \\ 0 & \lambda \end{array}\right]}_{\Lambda_3} \underbrace{\left[\begin{array}{cc} v & v' \end{array}\right]^{-1}}_{V^{-1}} }

ここで、v_1,v_2、v_R,v_I、v、v'はすべて２次元の実ベクトルです（v'は一般化固有ベクトル）。

これらの行列指数関数は次式となります。

\displaystyle{(5a)\quad \exp(A t)=\underbrace{\left[\begin{array}{cc} v_1 & v_2 \end{array}\right]}_{V} \underbrace{\left[\begin{array}{cc} e^{\lambda_1t} & 0 \\ 0 & e^{\lambda_2 t} \end{array}\right]}_{\exp(\Lambda_1 t)} \underbrace{\left[\begin{array}{cc} v_1 & v_2 \end{array}\right]^{-1}}_{V^{-1}} }

\displaystyle{(5b)\quad \exp(A t)=\underbrace{\left[\begin{array}{cc} v_R & v_I \end{array}\right]}_{V'} \underbrace{e^{\lambda_R t} \left[\begin{array}{cc} \cos(\lambda_It) & \sin(\lambda_It) \\ -\sin(\lambda_It) & \cos(\lambda_It) \end{array}\right]}_{\exp(\Lambda_2 t)} \underbrace{\left[\begin{array}{cc} v_R & v_I \end{array}\right]^{-1}}_{V'^{-1}} }

\displaystyle{(5c)\quad \exp(A t)=\underbrace{\left[\begin{array}{cc} v & v' \end{array}\right]}_{V} \underbrace{e^{\lambda t} \left[\begin{array}{cc} 1 & t \\ 0 & 1 \end{array}\right]}_{\exp(\Lambda_3 t)} \underbrace{\left[\begin{array}{cc} v & v' \end{array}\right]^{-1}}_{V^{-1}} }

これらを(2)に代入して

\displaystyle{(6a)\quad \begin{array}{l} x(t)=\underbrace{\left[\begin{array}{cc} v_1 & v_2 \end{array}\right]}_{V} \underbrace{\left[\begin{array}{cc} e^{\lambda_1t} & 0 \\ 0 & e^{\lambda_2 t} \end{array}\right]}_{\exp(\Lambda_1 t)} \underbrace{\left[\begin{array}{cc} v_1 & v_2 \end{array}\right]^{-1}x(0)}_{\left[\begin{array}{c} c_1\\c_2 \end{array}\right]}\\ =c_1e^{\lambda_1t}v_1+c_2e^{\lambda_2t}v_2 \end{array} }

\displaystyle{(6b)\quad \begin{array}{l} x(t)=\underbrace{\left[\begin{array}{cc} v_R & v_I \end{array}\right]}_{V'} \underbrace{e^{\lambda_R t} \left[\begin{array}{cc} \cos(\lambda_It) & \sin(\lambda_It) \\ -\sin(\lambda_It) & \cos(\lambda_It) \end{array}\right]}_{\exp(\Lambda_2 t)} \underbrace{\left[\begin{array}{cc} v_R & v_I \end{array}\right]^{-1}x(0)}_{\left[\begin{array}{c} c_1\\c_2 \end{array}\right]}\\ =c_1e^{\lambda_Rt}(\cos(\lambda_It)v_R-\sin(\lambda_It)v_I)+c_2e^{\lambda_Rt}(\sin(\lambda_It)v_R+\cos(\lambda_It)v_I) \end{array} }

\displaystyle{(6c)\quad \begin{array}{l} x(t)=\underbrace{\left[\begin{array}{cc} v & v' \end{array}\right]}_{V} \underbrace{e^{\lambda t} \left[\begin{array}{cc} 1 & t \\ 0 & 1 \end{array}\right]}_{\exp(\Lambda_3 t)} \underbrace{\left[\begin{array}{cc} v & v' \end{array}\right]^{-1}x(0)}_{\left[\begin{array}{c} c_1\\c_2 \end{array}\right]}\\ =c_1e^{\lambda t}v+c_2e^{\lambda t}(tv+v') \end{array} }

●一般の自由系についても、同様の式が成り立ち、Aの固有値と固有ベクトルがどのように時間応答に寄与しているかの分析ができ、モード間の連成について考察することができます。