制御系CAD

制御を意識したモノつくりを目指して

状態オブザーバ

投稿日時: 2021年10月10日投稿者: cacsd

状態オブザーバ…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$ の設計とは独立して行ってよいことを意味しており、大変重要な性質（分離定理）として知られています。

可制御正準形

投稿日時: 2021年10月9日投稿者: cacsd

可制御正準形…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}} }$

可制御標準形

投稿日時: 2021年10月9日投稿者: cacsd

可制御標準形…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

可安定性と可制御性

投稿日時: 2021年10月9日投稿者: cacsd

可安定性と可制御性…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$ でなければなりません。

状態フィードバック

投稿日時: 2021年10月9日投稿者: cacsd

状態フィードバック…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$ となる場合は、第１行と第２行が同じとなり、第３行と第４行も同じとなり、正則性は失われます。これから全く同じ２本の棒の棒立ては成功しないと考えられます。これは「二兎を追うものは一兎も得ず」のことわざ通りですね。

２次系の周波数応答

投稿日時: 2021年10月9日投稿者: cacsd

２次系の周波数応答…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))}$

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

次の２次系を考えます。

$\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{(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}}$

のように表されます。

２次系のステップ応答

投稿日時: 2021年10月9日投稿者: cacsd

２次系のステップ応答…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

ｎ次系の時間応答

投稿日時: 2021年10月9日投稿者: cacsd

ｎ次系の時間応答…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)

$\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)であると言います。たとえば入力を２倍したときの時間応答が元の時間応答の２倍にならなければ非線形系であると判断できます。

漸近安定性

投稿日時: 2021年10月9日投稿者: cacsd

漸近安定性…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$ の固有値と固有ベクトルがどのように時間応答に寄与しているかの分析ができ、モード間の連成について考察することができます。

ドローン

投稿日時: 2021年10月9日投稿者: cacsd

上図のようなドローンについて、次の文献をフォローしていきます。

Randal Beard: Quadrotor Dynamics and Control, Rev0.1, Brigham Young University, 2008

以下では、 $\sin(\cdot)$ 、 $\cos(\cdot)$ をそれぞれ $S_{(\cdot)}$ 、 $C_{(\cdot)}$ と略記します。

準備

●上図のように、ベクトル ${\bf p}$ を、ベクトル $\hat{n}$ のまわりに $\mu$ だけ回転して得られるベクトル ${\bf q}$ の表現式を求めます。

$\displaystyle{{\bf q}=\vec{ON}+\vec{NW}+\vec{WQ}}$
ただし
$\displaystyle{\vec{ON}=({\bf p}\cdot \hat{n}) \hat{n} }$
$\displaystyle{\vec{NW}=\frac{ {\bf p}-({\bf p}\cdot \hat{n}) \hat{n} }{ || {\bf p}-({\bf p}\cdot \hat{n}) \hat{n} ||}\overline{NQ}\,C_\mu =({\bf p}-({\bf p}\cdot \hat{n}) \hat{n})\, C_\mu \quad(\overline{NQ}=|| {\bf p}-({\bf p}\cdot \hat{n}) \hat{n} ||)}$
$\displaystyle{\vec{WQ}=\frac{{\bf p}\times\hat{n}}{||{\bf p}||S_\phi}\overline{NQ}\,S_\mu=-(\hat{n}\times{\bf p})\,S_\mu \quad(\overline{NQ}=||{\bf p}||S_\phi)}$

から、次式を得ます。

$\displaystyle{{\bf q}=(1-C_\mu)({\bf p}\cdot \hat{n}) \hat{n}+C_\mu{\bf p} -S_\mu(\hat{n}\times{\bf p})}$

ここで

$\displaystyle{{\bf p}=p_x^0\hat{i}^0+p_y^0\hat{j}^0+p_z^0\hat{k}^0}$
$\displaystyle{{\bf q}=q_x^0\hat{i}^0+q_y^0\hat{j}^0+q_z^0\hat{k}^0}$
$\displaystyle{\hat{n}=\hat{n}_x \hat{i}^0+\hat{n}_y \hat{j}^0+\hat{n}_z \hat{k}^0}$

とおくと

$\displaystyle{{\bf p}\cdot \hat{n}=p_x^0\hat{n}_x+p_y^0\hat{n}_y+p_z^0\hat{n}_z }$
$\displaystyle{\hat{n}\times{\bf p}= {\rm det} \left[\begin{array}{ccc} \hat{i}^0 & \hat{j}^0 & \hat{k}^0\\ \hat{n}_x & \hat{n}_y & \hat{n}_z\\ p_x^0 & p_y^0 & p_z^0 \end{array}\right]}$
$\displaystyle{=(\hat{n}_yp_z^0-\hat{n}_zp_y^0)\hat{i}^0+(\hat{n}_zp_x^0-\hat{n}_xp_z^0)\hat{j}^0+(\hat{n}_xp_y^0-\hat{n}_yp_x^0)\hat{k}^0 }$

● $\hat{n}=\hat{i}$ の場合（left-handed rotation about z-axis）

$\displaystyle{{\bf q}=(1-C_\psi)({\bf p}\cdot \hat{n}) \hat{n}+C_\psi{\bf p} -S_\psi(\hat{n}\times{\bf p}) }$
$\displaystyle{ \underbrace{ \left[\begin{array}{c} q_x^0 \\ q_y^0 \\ q_z^0 \end{array}\right] }_{{\bf q}} =(1-C_\psi) p_z^0 \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ 1 \end{array}\right] }_{\hat{i}} \hat{n}+C_\psi \left[\begin{array}{c} p_x^0 \\ p_y^0 \\ p_z^0 \end{array}\right] -S_\psi \left[\begin{array}{c} -p_y^0 \\ p_x^0 \\ 0 \end{array}\right] }$
$\displaystyle{ = \underbrace{ \left[\begin{array}{ccc} C_\psi & S_\psi & 0 \\ -S_\psi & C_\psi & 0 \\ 0 & 0 & 1 \end{array}\right] }_{R(\phi)} \underbrace{ \left[\begin{array}{c} p_x^0 \\ p_y^0 \\ p_z^0 \end{array}\right] }_{{\bf p}} }$

● $\hat{n}=\hat{j}$ の場合（left-handed rotation about y-axis)

$\displaystyle{{\bf q}=(1-C_\theta)({\bf p}\cdot \hat{n}) \hat{n}+C_\theta{\bf p} -S_\theta(\hat{n}\times{\bf p}) }$
$\displaystyle{ \underbrace{ \left[\begin{array}{c} q_x^0 \\ q_y^0 \\ q_z^0 \end{array}\right] }_{{\bf q}} =(1-C_\theta) p_y^0 \underbrace{ \left[\begin{array}{c} 0 \\ 1 \\ 0 \end{array}\right] }_{\hat{j}} \hat{n}+C_\theta \left[\begin{array}{c} p_x^0 \\ p_y^0 \\ p_z^0 \end{array}\right] -S_\theta \left[\begin{array}{c} p_z^0 \\ 0\\ -p_x^0 \\ \end{array}\right] }$
$\displaystyle{ =\underbrace{ \left[\begin{array}{ccc} C_\theta & 0 & -S_\theta \\ 0 & 1 & 0 \\ S_\theta & 0 & C_\theta \end{array}\right] }_{R(\theta)} \underbrace{ \left[\begin{array}{c} p_x^0 \\ p_y^0 \\ p_z^0 \end{array}\right] }_{{\bf p}} }$

● $\hat{n}=\hat{k}$ の場合（left-handed rotation about x-axis）

$\displaystyle{{\bf q}=(1-C_\phi)({\bf p}\cdot \hat{n}) \hat{n}+C_\phi{\bf p} -S_\phi(\hat{n}\times{\bf p}) }$
$\displaystyle{ \underbrace{ \left[\begin{array}{c} q_x^0 \\ q_y^0 \\ q_z^0 \end{array}\right] }_{{\bf q}} =(1-C_\phi) p_x^0 \underbrace{ \left[\begin{array}{c} 1 \\ 0 \\ 0 \end{array}\right] }_{\hat{k}} \hat{n}+C_\phi \left[\begin{array}{c} p_x^0 \\ p_y^0 \\ p_z^0 \end{array}\right] -S_\phi \left[\begin{array}{c} 0\\ -p_z^0 \\ p_y^0 \\ \end{array}\right] }$
$\displaystyle{ =\underbrace{ \left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & C_\phi & S_\phi \\ 0 & -S_\phi & C_\phi \end{array}\right] }_{R(\psi)} \underbrace{ \left[\begin{array}{c} p_x^0 \\ p_y^0 \\ p_z^0 \end{array}\right] }_{{\bf p}} }$

以上のようなベクトルの回転は、相対的に座標軸を回転させても、把握できます。

●right-handed rotation about z-axis

$\displaystyle{{\bf p}^0=p_x^0\hat{i}^0+p_y^0\hat{j}^0+p_z^0\hat{k}^0}$
$\displaystyle{{\bf p}^1=p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1}$
$\qquad\Downarrow$
$\displaystyle{p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1=p_x^0\hat{i}^0+p_y^0\hat{j}^0+p_z^0\hat{k}^0}$
$\qquad\Downarrow$
$\displaystyle{\hat{i}^1\cdot (p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1)=p_x^1=p_x^0\hat{i}^1\cdot\hat{i}^0+p_y^0\hat{i}^1\cdot\hat{j}^0+p_z^0\hat{i}^1\cdot\hat{k}^0}$
$\displaystyle{\hat{j}^1\cdot (p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1)=p_y^1=p_x^0\hat{j}^1\cdot\hat{i}^0+p_y^0\hat{j}^1\cdot\hat{j}^0+p_z^0\hat{j}^1\cdot\hat{k}^0}$
$\displaystyle{\hat{k}^1\cdot (p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1)=p_z^1=p_x^0\hat{k}^1\cdot\hat{i}^0+p_y^0\hat{k}^1\cdot\hat{j}^0+p_z^0\hat{k}^1\cdot\hat{k}^0}$
$\qquad\Downarrow$
$\displaystyle{ \underbrace{ \left[\begin{array}{c} p_x^1 \\ p_y^1 \\ p_z^1 \end{array}\right] }_{{\bf p}^1} = \underbrace{ \left[\begin{array}{ccc} \hat{i}^1\cdot\hat{i}^0&\hat{i}^1\cdot\hat{j}^0&\hat{i}^1\cdot\hat{k}^0\\ \hat{j}^1\cdot\hat{i}^0&\hat{j}^1\cdot\hat{j}^0&\hat{j}^1\cdot\hat{k}^0\\ \hat{k}^1\cdot\hat{i}^0&\hat{k}^1\cdot\hat{j}^0&\hat{k}^1\cdot\hat{k}^0 \end{array}\right] }_{R_0^1(\psi)} \underbrace{ \left[\begin{array}{c} p_x^0 \\ p_y^0 \\ p_z^0 \end{array}\right] }_{{\bf p}^0} }$
$\qquad\Downarrow$
$\displaystyle{ R_0^1(\psi)= \left[\begin{array}{ccc} C_\psi & S_\psi & 0 \\ -S_\psi & C_\psi & 0 \\ 0 & 0 & 1 \end{array}\right] }$

●right-handed rotation about y-axis

$\displaystyle{{\bf p}^0=p_x^0\hat{i}^0+p_y^0\hat{j}^0+p_z^0\hat{k}^0}$
$\displaystyle{{\bf p}^1=p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1}$
$\qquad\Downarrow$
$\displaystyle{p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1=p_x^0\hat{i}^0+p_y^0\hat{j}^0+p_z^0\hat{k}^0}$
$\qquad\Downarrow$
$\displaystyle{\hat{i}^1\cdot (p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1)=p_x^1=p_x^0\hat{i}^1\cdot\hat{i}^0+p_y^0\hat{i}^1\cdot\hat{j}^0+p_z^0\hat{i}^1\cdot\hat{k}^0}$
$\displaystyle{\hat{j}^1\cdot (p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1)=p_y^1=p_x^0\hat{j}^1\cdot\hat{i}^0+p_y^0\hat{j}^1\cdot\hat{j}^0+p_z^0\hat{j}^1\cdot\hat{k}^0}$
$\displaystyle{\hat{k}^1\cdot (p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1)=p_z^1=p_x^0\hat{k}^1\cdot\hat{i}^0+p_y^0\hat{k}^1\cdot\hat{j}^0+p_z^0\hat{k}^1\cdot\hat{k}^0}$
$\qquad\Downarrow$
$\displaystyle{ \underbrace{ \left[\begin{array}{c} p_x^1 \\ p_y^1 \\ p_z^1 \end{array}\right] }_{{\bf p}^1} = \underbrace{ \left[\begin{array}{ccc} \hat{i}^1\cdot\hat{i}^0&\hat{i}^1\cdot\hat{j}^0&\hat{i}^1\cdot\hat{k}^0\\ \hat{j}^1\cdot\hat{i}^0&\hat{j}^1\cdot\hat{j}^0&\hat{j}^1\cdot\hat{k}^0\\ \hat{k}^1\cdot\hat{i}^0&\hat{k}^1\cdot\hat{j}^0&\hat{k}^1\cdot\hat{k}^0 \end{array}\right] }_{R_0^1(\theta)} \underbrace{ \left[\begin{array}{c} p_x^0 \\ p_y^0 \\ p_z^0 \end{array}\right] }_{{\bf p}^0} }$
$\qquad\Downarrow$
$\displaystyle{ R_0^1(\theta)= \left[\begin{array}{ccc} C_\theta & 0 & -S_\theta \\ 0 & 1 & 0 \\ S_\theta & 0 & C_\theta \end{array}\right] }$

●right-handed rotation about x-axis

$\displaystyle{{\bf p}^0=p_x^0\hat{i}^0+p_y^0\hat{j}^0+p_z^0\hat{k}^0}$
$\displaystyle{{\bf p}^1=p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1}$
$\qquad\Downarrow$
$\displaystyle{p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1=p_x^0\hat{i}^0+p_y^0\hat{j}^0+p_z^0\hat{k}^0}$
$\qquad\Downarrow$
$\displaystyle{\hat{i}^1\cdot (p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1)=p_x^1=p_x^0\hat{i}^1\cdot\hat{i}^0+p_y^0\hat{i}^1\cdot\hat{j}^0+p_z^0\hat{i}^1\cdot\hat{k}^0}$
$\displaystyle{\hat{j}^1\cdot (p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1)=p_y^1=p_x^0\hat{j}^1\cdot\hat{i}^0+p_y^0\hat{j}^1\cdot\hat{j}^0+p_z^0\hat{j}^1\cdot\hat{k}^0}$
$\displaystyle{\hat{k}^1\cdot (p_x^1\hat{i}^1+p_y^1\hat{j}^1+p_z^1\hat{k}^1)=p_z^1=p_x^0\hat{k}^1\cdot\hat{i}^0+p_y^0\hat{k}^1\cdot\hat{j}^0+p_z^0\hat{k}^1\cdot\hat{k}^0}$
$\qquad\Downarrow$
$\displaystyle{ \underbrace{ \left[\begin{array}{c} p_x^1 \\ p_y^1 \\ p_z^1 \end{array}\right] }_{{\bf p}^1} = \underbrace{ \left[\begin{array}{ccc} \hat{i}^1\cdot\hat{i}^0&\hat{i}^1\cdot\hat{j}^0&\hat{i}^1\cdot\hat{k}^0\\ \hat{j}^1\cdot\hat{i}^0&\hat{j}^1\cdot\hat{j}^0&\hat{j}^1\cdot\hat{k}^0\\ \hat{k}^1\cdot\hat{i}^0&\hat{k}^1\cdot\hat{j}^0&\hat{k}^1\cdot\hat{k}^0 \end{array}\right] }_{R_0^1(\phi)} \underbrace{ \left[\begin{array}{c} p_x^0 \\ p_y^0 \\ p_z^0 \end{array}\right] }_{{\bf p}^0} }$
$\qquad\Downarrow$
$\displaystyle{ R_0^1(\phi)= \left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & C_\phi & S_\phi \\ 0 & -S_\phi & C_\phi \end{array}\right] }$

座標系

●ドローンの運動の記述にあたっては、次のような地上座標系 ${\cal F}^i$ と機体固定座標系 ${\cal F}^b$ を用います。

${\cal F}^i$ における機体の位置は次の変数で表します。

$p_n$ : the inertial (north) position of the aircraft along $\hat{i}^i$ in ${\cal F}^i$
$p_e$ : the inertial (east) position of the aircraft along $\hat{j}^i$ in ${\cal F}^i$
$h$ : the altitude of the aircraft measured along $-\hat{k}^i$ in ${\cal F}^i$

また、 ${\cal F}^b$ は、 ${\cal F}^i$ から、次のオイラー角

$\phi$ : the roll angle defined with respect to ${\cal F}^{v2}$
$\theta$ : the pitch angle defined with respect to ${\cal F}^{v1}$
$\psi$ : the yaw angle defined with respect to ${\cal F}^{v}$

を用いて、次のように逐次回転して得られます。

$\displaystyle{{\cal F}^i\rightarrow{\cal F}^v\rightarrow{\cal F}^{v1}\rightarrow{\cal F}^{v2}\rightarrow{\cal F}^b$

ここで、 ${\cal F}^{v}$ は原点を機体重心にもつ回転前の座標系（ ${\cal F}^{i}$ に平行）を表しています。

● ${\cal F}^b$ における任意ベクトル ${\bf p}$ （後出の位置ベクトル、速度ベクトルなど）は、 ${\cal F}^{i}$ から、次のような変換を行って得られます。

$\displaystyle{ {\bf p}^b=R_v^b(\phi,\theta,\psi){\bf p}^v }$
$\displaystyle{ R_v^b(\phi,\theta,\psi)=R_{v2}^b(\phi) R_{v1}^{v2}(\theta) R_v^{v1}(\psi) }$
$\displaystyle{= \left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & C_\phi & S_\phi \\ 0 & -S_\phi & C_\phi \end{array}\right] \left[\begin{array}{ccc} C_\theta & 0 & -S_\theta \\ 0 & 1 & 0 \\ S_\theta & 0 & C_\theta \end{array}\right] \left[\begin{array}{ccc} C_\psi & S_\psi & 0 \\ -S_\psi & C_\psi & 0 \\ 0 & 0 & 1 \end{array}\right] }$
$\displaystyle{ = \left[\begin{array}{ccc} C_\theta C_\psi & C_\theta S_\psi & -S_\theta \\ S_\phi S_\theta C_\psi - C_\phi S_\psi & S_\phi S_\theta S_\psi + C_\phi C_\psi & S_\phi C_\theta \\ C_\phi S_\theta C_\psi + C_\phi S_\psi & C_\phi S_\theta S_\psi - S_\phi C_\psi & C_\phi C_\theta \end{array}\right] }$

● $\displaystyle{{\cal F}^v\rightarrow{\cal F}^{v1}$
$\displaystyle{ {\bf p}^{v1}=R_v^{v1}(\psi){\bf p}^v }$
$\displaystyle{ R_v^{v1}(\psi)= \left[\begin{array}{ccc} C_\psi & S_\psi & 0 \\ -S_\psi & C_\psi & 0 \\ 0 & 0 & 1 \end{array}\right] }$

● $\displaystyle{{\cal F}^{v1}\rightarrow{\cal F}^{v2}$
$\displaystyle{ {\bf p}^{v2}=R_{v1}^{v2}(\theta){\bf p}^{v1} }$
$\displaystyle{ R_{v1}^{v2}(\theta)= \left[\begin{array}{ccc} C_\theta & 0 & -S_\theta \\ 0 & 1 & 0 \\ S_\theta & 0 & C_\theta \end{array}\right] }$

● $\displaystyle{{\cal F}^{v2}\rightarrow{\cal F}^b$
$\displaystyle{ {\bf p}^{b}=R_{v2}^b(\phi){\bf p}^{v2} }$
$\displaystyle{ R_{v2}^b(\phi)= \left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & C_\phi & S_\phi \\ 0 & -S_\phi & C_\phi \end{array}\right] }$

● ${\cal F}^b$ における任意ベクトル ${\bf p}$ （後出の位置ベクトル、速度ベクトルなど）の時間変化がどのように表されるかを調べます。

これは、次式で表されます。

$\displaystyle{\frac{d}{dt}_i{\bf p}=\frac{d}{dt}_b{\bf p}+\omega_{b/i}\times{\bf p} }$

実際

$\displaystyle{{\bf p}+\delta{\bf p}=(1-C_{(-\delta \phi)})({\bf p}\cdot \hat{s}) \hat{s}+C_{(-\delta \phi)}{\bf p} -S_{(-\delta \phi)}(\hat{s}\times{\bf p}) }$
$\displaystyle{\delta{\bf p}\simeq \delta \phi (\hat{s}\times{\bf p}) \Rightarrow \frac{\delta{\bf p}}{\delta t}\simeq\frac{\delta \phi}{\delta t} (\hat{s}\times{\bf p}) }$
$\qquad\Downarrow$
$\displaystyle{\frac{d}{dt}_i{\bf p}\simeq \underbrace{\dot{\phi} \hat{s}}_{\omega_{b/i}}\times{\bf p} }$

運動方程式

${\cal F}^b$ における運動は次の変数を用いて表されます。

$u$ : the body frame velocity measured along $\hat{i}^b$ in ${\cal F}^b$
$v$ : the body frame velocity measured along $\hat{j}^b$ in ${\cal F}^b$
$w$ : the body frame velocity measured along $\hat{k}^b$ in ${\cal F}^b$
$p$ : the roll rate measured along $\hat{i}^b$ in ${\cal F}^b$
$q$ : the pitch rate measured along $\hat{j}^b$ in ${\cal F}^b$
$r$ : the yaw rate measured along $\hat{k}^b$ in ${\cal F}^b$

これらが ${\cal F}^i$ からどう表されるか、すなわち次の変数の微分とどう関係するか（キネマティックス）を調べます。

$p_n$ : the inertial (north) position of the quadrotor along $\hat{i}^i$ in ${\cal F}^i$
$p_e$ : the inertial (east) position of the quadrotor along $\hat{j}^i$ in ${\cal F}^i$
$h$ : the altitude of the aircraft measured along $-\hat{k}^i$ in ${\cal F}^i$
$\phi$ : the roll angle defined with respect to ${\cal F}^{v2}$
$\theta$ : the pitch angle defined with respect to ${\cal F}^{v1}$
$\psi$ : the yaw angle defined with respect to ${\cal F}^{v}$

● ${\cal F}^i$ と ${\cal F}^b$ における速度ベクトルについては

$\displaystyle{ \left[\begin{array}{c} u \\ v \\ w \end{array}\right] = \underbrace{ \left[\begin{array}{ccc} C_\theta C_\psi & C_\theta S_\psi & -S_\theta \\ S_\phi S_\theta C_\psi - C_\phi S_\psi & S_\phi S_\theta S_\psi + C_\phi C_\psi & S_\phi C_\theta \\ C_\phi S_\theta C_\psi + C_\phi S_\psi & C_\phi S_\theta S_\psi - S_\phi C_\psi & C_\phi C_\theta \end{array}\right] }_{R_v^b} \left[\begin{array}{c} \dot{p}_n \\ \dot{p}_e \\ -\dot{h} \end{array}\right] }$

の関係が成り立ちます。したがって、

$\displaystyle{ \left[\begin{array}{c} \dot{p}_n \\ \dot{p}_e \\ -\dot{h} \end{array}\right] = \underbrace{ \left[\begin{array}{ccc} C_\theta C_\psi &S_\phi S_\theta C_\psi - C_\phi S_\psi &C_\phi S_\theta C_\psi + C_\phi S_\psi \\ C_\theta S_\psi &S_\phi S_\theta S_\psi + C_\phi C_\psi &C_\phi S_\theta S_\psi - S_\phi C_\psi \\ -S_\theta &S_\phi C_\theta &C_\phi C_\theta \end{array}\right] }_{R_b^v=(R_v^b)^T} \left[\begin{array}{c} u \\ v \\ w \end{array}\right] }$

● ${\cal F}^i$ と ${\cal F}^b$ における角速度ベクトルについては

$\displaystyle{ \left[\begin{array}{c} p \\ q \\ r \end{array}\right] = \underbrace{R_{v2}^b(\dot{\phi})}_{I} \left[\begin{array}{c} \dot{\phi} \\ 0 \\ 0 \end{array}\right] + R_{v2}^b(\phi)\underbrace{R_{v1}^{v2}(\dot{\theta})}_{I} \left[\begin{array}{c} 0 \\ \dot{\theta} \\ 0 \end{array}\right] + R_{v2}^b(\phi)R_{v1}^{v2}(\theta)\underbrace{R_{v}^{v1}(\dot{\psi})}_{I} \left[\begin{array}{c} 0 \\ 0 \\ \dot{\psi} \end{array}\right] }$
$\displaystyle{ = \left[\begin{array}{c} \dot{\phi} \\ 0 \\ 0 \end{array}\right] + \left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & C_\phi & S_\phi \\ 0 & -S_\phi & C_\phi \end{array}\right] \left[\begin{array}{c} 0 \\ \dot{\theta} \\ 0 \end{array}\right] + \left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & C_\phi & S_\phi \\ 0 & -S_\phi & C_\phi \end{array}\right] \left[\begin{array}{ccc} C_\theta & 0 & -S_\theta \\ 0 & 1 & 0 \\ S_\theta & 0 & C_\theta \end{array}\right] \left[\begin{array}{c} 0 \\ 0 \\ \dot{\psi} \end{array}\right] }$

の関係が成り立ちます。したがって、

$\displaystyle{ \left[\begin{array}{c} p \\ q \\ r \end{array}\right]= \left[\begin{array}{ccc} 1 & 0 & -S_\theta \\ 0 & C_\phi & S_\phi C_\theta \\ 0 & -S_\phi & C_\phi C_\theta \end{array}\right] \left[\begin{array}{c} \dot{\phi} \\ \dot{\theta} \\ \dot{\psi} \end{array}\right] }$

すなわち、次の関係式を得ます。

●次に並進運動（ダイナミックス）は次式で表されます。

$\displaystyle{m\frac{d}{dt}_i{\bf v}=\frac{d}{dt}_b{\bf v}+\omega_{b/i}\times{\bf v}=\frac{1}{m}{\bf f} }$
$\displaystyle{ \frac{d}{dt}_b \left[\begin{array}{c} u\\ v\\ w \end{array}\right] + \left[\begin{array}{c} p\\ q\\ r \end{array}\right] \times \left[\begin{array}{c} u\\ v\\ w \end{array}\right] =\frac{1}{m} \left[\begin{array}{c} f_x\\ f_y\\ f_z \end{array}\right] }$

$\displaystyle{ \left[\begin{array}{c} \dot{u}\\ \dot{v}\\ \dot{w} \end{array}\right] =- \underbrace{ \left[\begin{array}{ccc} 0 & -r & q\\ r & 0 &-p\\ -q & p & 0 \end{array}\right] \left[\begin{array}{c} u\\ v\\ w \end{array}\right] }_{\left[\begin{array}{c} qw-rv \\ ru-pw\\ pv-qu \end{array}\right]} +\frac{1}{m} \left[\begin{array}{c} f_x\\ f_y\\ f_z \end{array}\right] }$

●また回転運動（ダイナミックス）は次式で表されます。

$\displaystyle{\frac{d}{dt}_i{\bf h}=\frac{d}{dt}_b{\bf h}+\omega_{b/i}\times{\bf h}={\bf m} }$
$\displaystyle{{\bf h}=J\omega_{b/i} }$
$\displaystyle{J= \left[\begin{array}{ccc} \int(y^2+z^2)dm & -\int xydm & -\int xzdm \\ -\int yxdm & \int(x^2+z^2)dm & -\int yzdm \\ -\int zxdm & -\int zydm & \int(x^2+y^2)dm \end{array}\right] =\left[\begin{array}{ccc} J_x & 0 & 0 \\ 0 & J_y & 0 \\ 0 & 0 & J_z \end{array}\right] }$

$\displaystyle{ \frac{d}{dt}_bJ \left[\begin{array}{c} p\\ q\\ r \end{array}\right] + \left[\begin{array}{c} p\\ q\\ r \end{array}\right] \times J \left[\begin{array}{c} p\\ q\\ r \end{array}\right] = \left[\begin{array}{c} \tau_\phi\\ \tau_\theta\\ \tau_\psi \end{array}\right] }$
$\displaystyle{ \left[\begin{array}{c} J_x\dot{p}\\ J_y\dot{q}\\ J_z\dot{r} \end{array}\right] + \underbrace{ \left[\begin{array}{ccc} 0 & -r & q\\ r & 0 &-p\\ -q & p & 0 \end{array}\right] \left[\begin{array}{c} J_xp\\ J_yq\\ J_zr \end{array}\right] }_{\left[\begin{array}{c} (J_z-J_y)qr \\ (J_x-J_z)pr\\ (J_y-J_x)pq \end{array}\right]} = \left[\begin{array}{c} \tau_\phi\\ \tau_\theta\\ \tau_\psi \end{array}\right] }$

$\displaystyle{ \left[\begin{array}{c} \dot{p}\\ \dot{q}\\ \dot{r} \end{array}\right] =- \underbrace{ \left[\begin{array}{ccc} 0 & 0 & \frac{J_z-J_y}{J_x}q\\ \frac{J_x-J_z}{J_y} r & 0 & 0\\ 0 & \frac{J_y-J_x}{J_z}p & 0 \end{array}\right] \left[\begin{array}{c} p\\ q\\ r \end{array}\right] }_{ \left[\begin{array}{c} \frac{J_z-J_y}{J_x}qr \\ \frac{J_x-J_z}{J_y}pr \\ \frac{J_y-J_x}{J_z}pq \end{array}\right] } + \left[\begin{array}{c} \frac{1}{J_x}\tau_\phi\\ \frac{1}{J_y}\tau_\theta\\ \frac{1}{J_z}\tau_\psi \end{array}\right] }$

外力

４つのモータへの回転指令を $\delta_f,\delta_r,\delta_b,\delta_\ell$ とします。このとき、各モータによる推力とトルクは次式で表されます。

$\displaystyle{F_f=k_1\delta_f}$
$\displaystyle{F_r=k_1\delta_r}$
$\displaystyle{F_b=k_1\delta_b}$
$\displaystyle{F_\ell=k_1\delta_\ell}$

$\displaystyle{\tau_f=k_2\delta_f}$
$\displaystyle{\tau_r=k_2\delta_r}$
$\displaystyle{\tau_b=k_2\delta_b}$
$\displaystyle{\tau_\ell=k_2\delta_\ell}$

また、４つのモータの合推力と各軸回りのトルクは次式で表されます。

$\displaystyle{F=F_f+F_r+F_b+F_\ell}$
$\displaystyle{\tau_\phi=\ell(F_\ell-F_r)}$
$\displaystyle{\tau_\theta=\ell(F_f-F_b)}$
$\displaystyle{\tau_\psi=\tau_r+\tau_\ell-\tau_f-\tau_b}$

すなわち

$\displaystyle{ \left[\begin{array}{c} F\\ \tau_\phi\\ \tau_\theta\\ \tau_\psi \end{array}\right] = \underbrace{ \left[\begin{array}{cccc} k_1 & k_1 & k_1 & k_1 \\ 0 & -\ell k_1 & 0 & \ell k_1 \\ \ell k_1 & 0 & -\ell k_1 & 0 \\ -k_2 & k_2 & -k_2 & k_2 \end{array}\right] }_{{\cal M}} \left[\begin{array}{c} \delta_f\\ \delta_r\\ \delta_b\\ \delta_\ell \end{array}\right] }$

また、並進運動方程式の外力項は次式となります。

$\displaystyle{ \left[\begin{array}{c} f_x\\ f_y\\ f_z \end{array}\right] = R_v^b \left[\begin{array}{c} 0\\ 0\\ mg \end{array}\right] + \left[\begin{array}{c} 0\\ 0\\ -F \end{array}\right] =mg \left[\begin{array}{ccc} -S_\theta \\ S_\phi C_\theta \\ C_\phi C_\theta \end{array}\right] + \left[\begin{array}{c} 0\\ 0\\ -F \end{array}\right] }$

このときドローンの運動方程式に基づく状態方程式は次式となります。

$\displaystyle{ \left[\begin{array}{c} \dot{p}_n \\ \dot{p}_e \\ -\dot{h} \\ \hline \dot{\phi} \\ \dot{\theta} \\ \dot{\psi} \\ \hline \dot{u}\\ \dot{v}\\ \dot{w} \\ \hline \dot{p}\\ \dot{q}\\ \dot{r} \end{array}\right] = \left[\begin{array}{ccc|ccc} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array}\right.}$
$\displaystyle{ \left.\begin{array}{ccc|ccc} C_\theta C_\psi &S_\phi S_\theta C_\psi - C_\phi S_\psi &C_\phi S_\theta C_\psi + C_\phi S_\psi & 0 & 0 & 0\\ C_\theta S_\psi &S_\phi S_\theta S_\psi + C_\phi C_\psi &C_\phi S_\theta S_\psi - S_\phi C_\psi & 0 & 0 & 0\\ -S_\theta &S_\phi C_\theta &C_\phi C_\theta & 0 & 0 & 0\\ \hline 0 & 0 & 0 & 1 & S_\phi T_\theta & C_\phi T_\theta \\ 0 & 0 & 0 & 0 & C_\phi & -S_\phi \\ 0 & 0 & 0 & 0 & \frac{S_\phi}{C_\theta} & \frac{C_\phi}{C_\theta}\\ \hline 0 & r &- q & 0 & 0 & 0\\ r- & 0 &p & 0 & 0 & 0\\ q & -p & 0 & 0 & 0 & 0\\ \hline 0 & 0 & 0 & 0 & 0 & -\frac{J_z-J_y}{J_x}q\\ 0 & 0 & 0 & -\frac{J_x-J_z}{J_y} r & 0 & 0\\ 0 & 0 & 0 & 0 & -\frac{J_y-J_x}{J_z}p & 0 \end{array}\right] }$
$\displaystyle{ \times \left[\begin{array}{c} p_n \\ p_e \\ -h \\ \hline \phi \\ \theta \\ \psi \\ \hline u\\ v\\ w \\ \hline p\\ q\\ r \end{array}\right] + \left[\begin{array}{cccc} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \frac{1}{m} & 0 & 0 & 0 \\ \hline 0 & \frac{1}{J_x} & 0 & 0 \\ 0 & 0 & \frac{1}{J_y} & 0 \\ 0 & 0 & 0 & \frac{1}{J_z} \\ \end{array}\right] \left[\begin{array}{cccc} k_1 & k_1 & k_1 & k_1 \\ 0 & -\ell k_1 & 0 & \ell k_1 \\ \ell k_1 & 0 & -\ell k_1 & 0 \\ -k_2 & k_2 & -k_2 & k_2 \end{array}\right]^{-1} \left[\begin{array}{c} \delta_f\\ \delta_r\\ \delta_b\\ \delta_\ell \end{array}\right] + \left[\begin{array}{c} 0 \\ 0 \\ 0 \\ \hline 0 \\ 0 \\ 0 \\ \hline -gS_\theta \\ gS_\phi C_\theta \\ gC_\phi C_\theta \\ \hline 0 \\ 0 \\ 0 \end{array}\right] }$

これを見ると、 $A$ 行列に状態変数が絡んでおり、複雑なダイナミクスをもっていることが分かります。

運動方程式の簡単化

次の仮定のもとで、運動方程式の簡単化を行ないます。

仮定１：　 $qr$ 、 $pr$ 、 $pq$ は小さく無視できる
仮定２：　 $\phi$ 、 $\theta$ は小さく無視できる

まず、仮定１のもとで、次の近似が可能です。

$\displaystyle{ \left[\begin{array}{c} \dot{p}\\ \dot{q}\\ \dot{r} \end{array}\right] =- \left[\begin{array}{c} \frac{J_z-J_y}{J_x}qr \\ \frac{J_x-J_z}{J_y}pr \\ \frac{J_y-J_x}{J_z}pq \end{array}\right] + \left[\begin{array}{c} \frac{1}{J_x}\tau_\phi\\ \frac{1}{J_y}\tau_\theta\\ \frac{1}{J_z}\tau_\psi \end{array}\right] \simeq \left[\begin{array}{c} \frac{1}{J_x}\tau_\phi\\ \frac{1}{J_y}\tau_\theta\\ \frac{1}{J_z}\tau_\psi \end{array}\right] }$

このとき、仮定２をおけば次の近似が可能です。

$\displaystyle{ \left[\begin{array}{c} \dot{\phi} \\ \dot{\theta} \\ \dot{\psi} \end{array}\right] = \left[\begin{array}{ccc} 1 & S_\phi T_\theta & C_\phi T_\theta \\ 0 & C_\phi & -S_\phi \\ 0 & \frac{S_\phi}{C_\theta} & \frac{C_\phi}{C_\theta} \end{array}\right] \left[\begin{array}{c} p \\ q \\ r \end{array}\right] \simeq \left[\begin{array}{c} p \\ q \\ r \end{array}\right] }$
$\displaystyle{ \left[\begin{array}{c} \ddot{\phi} \\ \ddot{\theta} \\ \ddot{\psi} \end{array}\right] \simeq \left[\begin{array}{c} \dot{p}\\ \dot{q}\\ \dot{r} \end{array}\right] \simeq \left[\begin{array}{c} \frac{1}{J_x}\tau_\phi\\ \frac{1}{J_y}\tau_\theta\\ \frac{1}{J_z}\tau_\psi \end{array}\right] }$

これは、機体の傾きを独立して制御できることを示しています。

さらに次のような近似が可能です。

$\displaystyle{ \left[\begin{array}{c} \dot{p}_x \\ \dot{p}_y \\ \dot{p}_z \end{array}\right] = {R}_{v1}^b \left[\begin{array}{c} u \\ v \\ w \end{array}\right] }$
$\displaystyle{ \left[\begin{array}{c} \ddot{p}_x \\ \ddot{p}_y \\ \ddot{p}_z \end{array}\right] = \underbrace{\dot{R}_{v1}^b}_{\simeq0} \left[\begin{array}{c} {u}\\ {v}\\ {w} \end{array}\right] + {R}_{v1}^b \left[\begin{array}{c} \dot{u}\\ \dot{v}\\ \dot{w} \end{array}\right] \simeq {R}_{v1}^b \left[\begin{array}{c} \dot{u}\\ \dot{v}\\ \dot{w} \end{array}\right] }$
$\displaystyle{ \left[\begin{array}{c} \ddot{p}_x \\ \ddot{p}_y \\ \ddot{p}_z \end{array}\right] \simeq \left[\begin{array}{ccc} C_\theta & S_\phi S_\theta & C_\phi S_\theta \\ 0 & C_\phi & -S_\phi \\ -S_\theta &S_\phi C_\theta &C_\phi C_\theta \end{array}\right] (g \left[\begin{array}{c} -S_\theta \\ S_\phi C_\theta \\ C_\phi C_\theta \end{array}\right] +\frac{1}{m} \left[\begin{array}{c} 0\\ 0\\ -F \end{array}\right]) }$
$\displaystyle{ = \left[\begin{array}{c} 0 \\ 0 \\ g \end{array}\right] -\frac{1}{m}F \left[\begin{array}{ccc} C_\phi S_\theta \\ -S_\phi \\ C_\phi C_\theta \end{array}\right] \simeq \left[\begin{array}{c} 0 \\ 0 \\ g-\frac{1}{m}F \end{array}\right] }$

結局、ドローンの運動方程式は、次式のように簡単化されます。

$\displaystyle{m\ddot{p}_z=mg-(F_f+F_r+F_b+F_\ell)}$
$\displaystyle{J_x\ddot{\phi}=\ell(F_\ell-F_r)}$
$\displaystyle{J_y\ddot{\theta}=\ell(F_f-F_b)}$
$\displaystyle{J_z\ddot{\psi}=\tau_r+\tau_\ell-\tau_f-\tau_b}$

いま，重力補償のために $mg=4k_f\delta^*$ から決まる $\delta^*$ を各操作入力に前もって加えておくことにすると、次式が成り立ちます。

$\displaystyle{m\ddot{p}_z=4k_f\delta^*-(k_f\delta_f+k_f\delta_r+k_f\delta_b+k_f\delta_\ell)}$
$\displaystyle{=-k_f(\delta_f-\delta^*)-k_f(\delta_r-\delta^*)-k_f(\delta_b-\delta^*)-k_f(\delta_\ell-\delta^*))}$

$\displaystyle{J_x\ddot{\phi}=\ell(k_f\delta_\ell-k_f\delta_r) =\ell(k_f(\delta_\ell-\delta^*)-k_f(\delta_r-\delta^*))}$

$\displaystyle{J_y\ddot{\theta}=\ell(k_f\delta_f-k_f\delta_b)=\ell(k_f(\delta_f-\delta^*)-k_f(\delta_b-\delta^*))}$

$\displaystyle{J_z\ddot{\psi}=(k_\tau\delta_r+k_\tau\delta_\ell)-(k_\tau\delta_f+k_\tau\delta_b)}$
$\displaystyle{=(k_\tau(\delta_r-\delta^*)+k_\tau(\delta_\ell-\delta^*))-(k_\tau(\delta_f-\delta^*)+k_\tau(\delta_b-\delta^*))}$

このときドローンの線形状態方程式は次式となります。

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{p}_z \\ \dot{\phi} \\ \dot{\theta} \\ \dot{\psi} \\\hline \ddot{z} \\ \ddot{\phi} \\ \ddot{\theta} \\ \ddot{\psi} \\ \end{array}\right] }_{\dot{x}} = \underbrace{ \left[\begin{array}{cccc|cccc} 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\\hline 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array}\right] }_{A} % \left[\begin{array}{c|c} % 0_{4\times4} & I_4 \\\hline % 0_{4\times4} & 0_{4\times4} % \end{array}\right] \underbrace{ \left[\begin{array}{c} {p_z} \\ {\phi} \\ {\theta} \\ {\psi} \\\hline \dot{z} \\ \dot{\phi} \\ \dot{\theta} \\ \dot{\psi} \\ \end{array}\right] }_{x} }$
$\displaystyle{+ \underbrace{ \left[\begin{array}{cccc} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\\hline -\frac{k_f}{m} & -\frac{k_f}{m} & -\frac{k_f}{m} & -\frac{k_f}{m} \\ 0 & -\frac{\ell k_f}{J_x} & 0 & \frac{\ell k_f}{J_x} \\ \frac{\ell k_f}{J_y} & 0 & -\frac{\ell k_f}{J_y} & 0 \\ -\frac{k_\tau}{J_z} & \frac{k_\tau}{J_z} & -\frac{k_\tau}{J_z} & \frac{k_\tau}{J_z} \\ \end{array}\right] }_{B} % \left[\begin{array}{c} % 0_{4\times4} \\ % B_2 \\ % \end{array}\right] \underbrace{ \left[\begin{array}{c} \delta_f-\delta^* \\ \delta_r-\delta^* \\ \delta_b-\delta^* \\ \delta_\ell-\delta^* \\ \end{array}\right] }_{u} }$