制御系CAD

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

LQG制御

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

LQG制御…Homework

[1]　次のようなオブザーバベース・コントローラによる閉ループ系を考えます。ここで，新しい入力 $w$ と $v$ がそれぞれ $W>0$ と $V>0$ の平方根行列（ $X>0(\ge0)$ に対し $YY=X$ を満足する行列 $Y>0(\ge0)$ を $X^{\frac{1}{2}}$ で表す）により重み付けられて， $n$ 次系の入力側（ $B'$ を介して）と出力側に設置されています。また新しい出力 $z=C'x$ と入力 $u$ が取り出されており，それぞれ $Q>0$ と $R>0$ の平方根行列により重み付けられています。

図１　LQG制御系設計の枠組み

可制御かつ可観測な制御対象

$\displaystyle{(1)\quad \left\{\begin{array}{ll} \dot{x}(t)=Ax(t)+Bu(t)+B'W^{\frac{1}{2}}w(t)&(x(t)\in{\bf R}^n,u(t),w(t)\in{\bf R}^m)\\ y(t)=Cx(t)+V^{\frac{1}{2}}v(t)&(y(t),v(t)\in{\bf R}^p)\\ z(t)=C'x(t)&(z(t)\in{\bf R}^p) \end{array}\right. }$

を安定化するオブザーバベース・コントローラ

$\displaystyle{(2)\quad \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. }$

を２次形式評価関数

$\displaystyle{(3)\quad J=\int_0^\infty(z^T(t)Qz(t)+u^T(t)Ru(t))\,dt\quad (Q>0,R>0) }$

を最小にするように、 $F$ と $H$ を決定する問題を考えます。これらは

$\displaystyle{(4)\quad \boxed{{\bf CARE}:\Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C'^TQC'=0} }$

$\displaystyle{(5)\quad \boxed{{\bf FARE}:\Gamma A^T+A\Gamma-\Gamma C^TV^{-1}C\Gamma+B'WB'^T=0} }$

を満足する $\Pi>0$ と $\Gamma>0$ を用いて、次のように与えられます。

$\displaystyle{(6)\quad \boxed{F=R^{-1}B^T\Pi} }$

$\displaystyle{(7)\quad \boxed{H=(V^{-1}C\Gamma)^T} }$

このような制御方式をLQG制御と呼びます。

[2]　上の安定な閉ループ系は次式で表されます。

$\displaystyle{(8)\quad \boxed{\begin{array}{l} \left[\begin{array}{cc} \dot{x}(t) \\ \dot{\hat{x}}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A & -BF \\ HC&A-HC-BF \end{array}\right] }_{A_{CL1}} \left[\begin{array}{cc} {x}(t) \\ \hat{x}(t) \end{array}\right] + \underbrace{ \left[\begin{array}{cc} B'W^{\frac{1}{2}} & 0 \\ 0 & HV^{\frac{1}{2}} \end{array}\right] }_{B_{CL1}} \left[\begin{array}{cc} {w}(t) \\ {v}(t) \end{array}\right]\\ \left[\begin{array}{cc} {y'}(t) \\ {u'}(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} Q^{\frac{1}{2}}C' & 0 \\ 0 & -R^{\frac{1}{2}}F \end{array}\right] }_{C_{CL1}} \left[\begin{array}{cc} {x}(t) \\ \hat{x}(t) \end{array}\right] \end{array}} }$

これに座標変換

$\displaystyle{(9)\quad \left[\begin{array}{cc} {x}(t) \\ {e}(t) \end{array}\right] = \left[\begin{array}{cc} I_n & 0 \\ -I_n & I_n \end{array}\right] \left[\begin{array}{cc} {x}(t) \\ \hat{x}(t) \end{array}\right],\quad \left[\begin{array}{cc} {x}(t) \\ \hat{x}(t) \end{array}\right] = \left[\begin{array}{cc} I_n & 0 \\ I_n & I_n \end{array}\right] \left[\begin{array}{cc} {x}(t) \\ {e}(t) \end{array}\right] }$

を行うと

$\displaystyle{(10)\quad \boxed{\begin{array}{l} \left[\begin{array}{cc} \dot{x}(t) \\ \dot{e}(t) \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A-BF & -BF \\ 0 & A-HC \end{array}\right] }_{A_{CL2}} \left[\begin{array}{cc} {x}(t) \\ {e}(t) \end{array}\right] +\underbrace{ \left[\begin{array}{cc} B'W^{\frac{1}{2}} & 0 \\ -B'W^{\frac{1}{2}} & HV^{\frac{1}{2}} \end{array}\right] }_{B_{CL2}} \left[\begin{array}{cc} {w}(t) \\ {v}(t) \end{array}\right]\\ \left[\begin{array}{cc} {y'}(t) \\ {u'}(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} Q^{\frac{1}{2}}C' & 0 \\ -R^{\frac{1}{2}}F & -R^{\frac{1}{2}}F \end{array}\right] }_{C_{CL2}} \left[\begin{array}{cc} {x}(t) \\ {e}(t) \end{array}\right] \end{array}} }$

となります。このとき閉ループ系のインパルス応答は次式で与えられます。

$\displaystyle{(11)\quad G_{CL}(t) =C_{CL1} \exp(A_{CL1}t) B_{CL1} =C_{CL2} \exp(A_{CL2}t) B_{CL2} }$

● $i=1,2$ に対して、評価関数 $J$ の総和は次式で与えられます。

$(12)\quad \begin{array}{l} \displaystyle{\int_0^\infty{\rm tr}(G_{CL}^T(t)G_{CL}(t))\,dt\quad (Q>0,R>0)}\\ \displaystyle{={\rm tr}(\int_0^\infty B_{CLi}^T\exp(A_{CLi}^Tt)C_{CLi}^TC_{CLi}\exp(A_{CLi}t)B_{CLi}\,dt)}\\ \displaystyle{={\rm tr}(B_{CLi}^T\int_0^\infty \exp(A_{CLi}^Tt)\underbrace{C_{CLi}^TC_{CLi}}_{Q_{CLi}}\exp(A_{CLi}t)\,dt\,B_{CLi})}\\ \displaystyle{={\rm tr}(B_{CLi}^T\underbrace{\int_0^\infty \exp(A_{CLi}^Tt)Q_{CLi}\exp(A_{CLi}t)\,dt}_{\Pi}\,B_{CLi})}\\ \displaystyle{={\rm tr}(\Pi B_{CLi}B_{CLi}^T)} \end{array} }$

ここで $\Pi>0$ は次式を満足します。

$\displaystyle{(13)\quad \Pi A_{CLi}+A_{CLi}^T\Pi+Q_{CLi}=0 }$

ラグランジュの未定定数法を適用するために、次の評価関数を考えます。

$\displaystyle{(14)\quad J'={\rm tr}(\Pi B_{CLi}B_{CLi}^T)+{\rm tr}((\Pi A_{CLi}+A_{CLi}^T\Pi+Q_{CLi})\Gamma) }$

これを最小化する場合の必要条件は、次式となります。

$\displaystyle{(15)\quad \frac{\partial J'}{\partial\Gamma}=0,\ \frac{\partial J'}{\partial\Pi}=0,\ \frac{\partial J'}{\partial F}=0,\ \frac{\partial J'}{\partial H}=0 }$

以下では、次の分割を考えます。

$\displaystyle{(16)\quad \Pi=\left[\begin{array}{cc} \Pi_{1}&\Pi_{3} \\ \Pi_{3}^T&\Pi_{2} \end{array}\right],\ \Gamma=\left[\begin{array}{cc} \Gamma_{1}&\Gamma_{3} \\ \Gamma_{3}^T&\Gamma_{2} \end{array}\right] }$

[3]　 $i=1$ の場合について、必要条件を一つ一つ調べていきます。

● $\frac{\partial J'}{\partial\Gamma}=0$

$\displaystyle{(17)\quad \frac{\partial J}{\partial\Gamma}=\Pi A_{CL1}+A_{CL1}^T\Pi+Q_{CL1}=0 }$

において

$\displaystyle{(18)\quad \begin{array}{l} \Pi A_{CL1}= \left[\begin{array}{cc} \Pi_{1}&\Pi_{3} \\ \Pi_{3}^T&\Pi_{2} \end{array}\right] \left[\begin{array}{cc} A & -BF \\ HC&A-HC-BF \end{array}\right]\\ = \left[\begin{array}{cc} \Pi_{1}A+\Pi_{3}HC & -\Pi_{1}BF+\Pi_{3}(A-HC-BF)\\ \Pi_{3}^TA+\Pi_{2}HC & -\Pi_{3}^TBF+\Pi_{2}(A-HC-BF) \end{array}\right] \end{array} }$

$\displaystyle{(19)\quad \begin{array}{l} A_{CL1}^T\Pi =(\Pi A_{CL1})^T\\ =\left[\begin{array}{cc} A^T\Pi_{1}+C^TH^T\Pi_{3}^T&\\ -F^TB^T\Pi_{1}+(A-HC-BF)^T\Pi_{3}^T & \end{array}\right.\\ \left.\begin{array}{cc} A^T\Pi_{3}+C^TH^T\Pi_{2}\\ -F^TB^T\Pi_{3}+(A-HC-BF)^T\Pi_{2} \end{array}\right] \end{array} }$

$\displaystyle{(20)\quad Q_{CL1}= \left[\begin{array}{cc} C'^TQC' & 0 \\ 0 & F^TRF \end{array}\right] }$

を代入して、次を得ます。

$\displaystyle{(21)\quad \Pi_{1}A+\Pi_{3}HC+A^T\Pi_{1}+C^TH^T\Pi_{3}^T+C'^TQC'=0 }$

$\displaystyle{(22)\quad -\Pi_{1}BF+\Pi_{3}(A-HC-BF)+A^T\Pi_{3}+C^TH^T\Pi_{2}=0 }$

$\displaystyle{(23)\quad \begin{array}{l} -\Pi_{3}^TBF+\Pi_{2}(A-HC-BF) \\ -F^TB^T\Pi_{3}+(A-HC-BF)^T\Pi_{2}+F^TRF=0 \end{array} }$

(22)+(23)より

$\displaystyle{(24)\quad \begin{array}{l} (\Pi_{3}+\Pi_{2})(A-HC)+(A-BF)^T(\Pi_{3}+\Pi_{2})\\ +\underbrace{F^TRF-(\Pi_{1}+\Pi_{3}+\Pi_{3}^T+\Pi_{2})BF}_{0\ (assumed)}=0\\ \Rightarrow \Pi_{3}+\Pi_{2}=0\Rightarrow \Pi_{3}=-\Pi_{2}\Rightarrow \Pi_{3}^T=\Pi_{3} \end{array} }$

(23)より

$\displaystyle{(25)\quad \begin{array}{l} \Pi_{2}(A-HC)+(A-HC)^T\Pi_{2}+F^TRF\\ \underbrace{-(\Pi_{3}^T+\Pi_{2})BF-F^TB^T(\Pi_{3}+\Pi_{2})}_{0}=0\\ \Rightarrow \Pi_{2}>0 \end{array} }$

(21)+(22)より

$\displaystyle{(26)\quad \begin{array}{l} (\Pi_{1}+\Pi_{3})(A-BF)+A^T(\Pi_{1}+\Pi_{3})\\ +\underbrace{C^TH^T(\Pi_{3}^T+\Pi_{2})}_{0}+C'^TQC'+F^TRF-F^TB^T(\Pi_{1}+\Pi_{3})=0 \end{array} }$

$\displaystyle{(27)\quad \begin{array}{l} (\Pi_{1}+\Pi_{3})(A-BF)+(A-BF)^T(\Pi_{1}+\Pi_{3})+C'^TQC'+F^TRF=0\\ \Rightarrow \Pi_{1}+\Pi_{3}=\Pi_{1}-\Pi_{2}>0 \end{array} }$

● $\frac{\partial J'}{\partial\Pi}=0$

$\displaystyle{(28)\quad \frac{\partial J}{\partial\Pi}=B_{CL1}B_{CL1}^T+\Gamma A_{CL1}^T+A_{CL1}\Gamma=0 }$

において

$\displaystyle{(29)\quad \begin{array}{l} A_{CL1}\Gamma= \left[\begin{array}{cc} A & -BF \\ HC&A-HC-BF \end{array}\right] \left[\begin{array}{cc} \Gamma_{1}&\Gamma_{3} \\ \Gamma_{3}^T&\Gamma_{2} \end{array}\right]\\ =\left[\begin{array}{cc} A\Gamma_{1} -BF\Gamma_{3}^T &A\Gamma_{3} -BF\Gamma_{2}\\ HC\Gamma_{1}+(A-HC-BF)\Gamma_{3}^T &HC\Gamma_{3}+(A-HC-BF)\Gamma_{2} \end{array}\right] \end{array} }$

$\displaystyle{(30)\quad \begin{array}{l} \Gamma A_{CL1}^T =(A_{CL1}\Gamma)^T\\ =\left[\begin{array}{cc} \Gamma_{1}A^T -\Gamma_{3}F^TB^T & \Gamma_{1}C^TH^T+\Gamma_{3}(A-HC-BF)^T\\ \Gamma_{3}^TA^T -\Gamma_{2}F^TB^T & \Gamma_{3}^TC^TH^T+\Gamma_{2}(A-HC-BF)^T \end{array}\right] \end{array} }$

$\displaystyle{(31)\quad B_{CL1}B_{CL1}^T =\left[\begin{array}{cc} B'WB'^T & 0 \\ 0 & HVH^T \end{array}\right] }$

を代入して、次を得ます。

$\displaystyle{(32)\quad A\Gamma_{1} -BF\Gamma_{3}^T+\Gamma_{1}A^T -\Gamma_{3}F^TB^T+B'WB'^T=0 }$

$\displaystyle{(33)\quad A\Gamma_{3} -BF\Gamma_{2}+\Gamma_{1}C^TH^T+\Gamma_{3}(A-HC-BF)^T=0 }$

$\displaystyle{(34)\quad \begin{array}{l} HC\Gamma_{3}+(A-HC-BF)\Gamma_{2}\\ +\Gamma_{3}^TC^TH^T+\Gamma_{2}(A-HC-BF)^T+HVH^T=0 \end{array} }$

(33)-(34)より

$\displaystyle{(35)\quad \begin{array}{l} (A-HC)(\Gamma_{3}-\Gamma_{2})+(\Gamma_{3}-\Gamma_{2})(A-BF)^T\\ +\underbrace{(\Gamma_{1}-\Gamma_{3}-\Gamma_{3}^T+\Gamma_{2})C^TH^T-HVH^T}_{0\ (assumed)}=0\\ \Rightarrow \Gamma_{3}-\Gamma_{2}=0\Rightarrow \Gamma_{3}=\Gamma_{2}\Rightarrow \Gamma_{3}^T=\Gamma_{3} \end{array} }$

(34)より

$\displaystyle{(36)\quad \begin{array}{l} (A-BF)\Gamma_{2}+\Gamma_{2}(A-BF)^T+H^TVH\\ +\underbrace{HC(\Gamma_{3}-\Gamma_{2})+(\Gamma^T_{3}-\Gamma_{2})C^TH^T}_{0}=0\Rightarrow \Gamma_{2}>0 \end{array} }$

(32)-(33)より

$\displaystyle{(37)\quad \begin{array}{l} (\Gamma_{1}-\Gamma_{3})(A-HC)^T+A(\Gamma_{1}-\Gamma_{3})\\ +\underbrace{BF(\Gamma_{3}^T-\Gamma_{2})}_{0}+B'WB'^T+HVH^T-HC(\Gamma_{1}-\Gamma_{3})=0 \end{array} }$

$\displaystyle{(38)\quad \begin{array}{l} (\Gamma_{1}-\Gamma_{3})(A-HC)^T+(A-HC)(\Gamma_{1}-\Gamma_{3})+B'WB'^T+HVH^T=0\\ \Rightarrow \Gamma_{1}-\Gamma_{3}=\Gamma_{1}-\Gamma_{2}>0 \end{array} }$

●準備１

$\displaystyle{(39)\quad \Pi B_{CL1}B_{CL1}^T= \left[\begin{array}{cc} \Pi_{1}&\Pi_{3} \\ \Pi_{3}^T&\Pi_{2} \end{array}\right] \left[\begin{array}{cc} B'WB'^T & 0 \\ 0 & HVH^T \end{array}\right] }$

(1) $\begin{eqnarray*} \end{eqnarray*}$

$\displaystyle{(40)\quad {\rm tr}(\Pi B_{CL1}B_{CL1}^T)={\rm tr}(\Pi_{1}B'WB'^T+\Pi_{2}HVH^T) }$

$\displaystyle{(41)\quad \frac{\partial {\rm tr}(\Pi B_{CL1}B_{CL1}^T)}{\partial H}=2\Pi_{2}HV }$

●準備２

$\displaystyle{(42)\quad \begin{array}{l} \Pi A_{CL1}\Gamma= \left[\begin{array}{cc} \Pi_{1}&\Pi_{3} \\ \Pi_{3}^T&\Pi_{2} \end{array}\right] \left[\begin{array}{cc} A & -BF \\ HC&A-HC-BF \end{array}\right] \Gamma\\ =\left[\begin{array}{cc} \Pi_{1}A+\Pi_{3}HC & -\Pi_{1}BF+\Pi_{3}(A-HC-BF)\\ \Pi_{3}^TA+\Pi_{2}HC & -\Pi_{3}^TBF+\Pi_{2}(A-HC-BF) \end{array}\right] \left[\begin{array}{cc} \Gamma_{1}&\Gamma_{3} \\ \Gamma_{3}^T&\Gamma_{2} \end{array}\right] \end{array} }$

$\displaystyle{(43)\quad \begin{array}{l} {\rm tr}(\Pi A_{CL1}\Gamma)\\ ={\rm tr}(\Pi_{1}A\Gamma_{1}+\Pi_{3}HC\Gamma_{1} -\Pi_{1}BF\Gamma_{3}^T+\Pi_{3}(A-HC-BF)\Gamma_{3}^T)\\ +{\rm tr}(\Pi_{3}^TA\Gamma_{3}+\Pi_{2}HC\Gamma_{3} -\Pi_{3}^TBF\Gamma_{2}+\Pi_{2}(A-HC-BF)\Gamma_{2}) \end{array} }$

$\displaystyle{(44)\quad \frac{\partial {\rm tr}(\Pi A_{CL1}\Gamma)}{\partial F} =-B^T\Pi_1\Gamma_{3}-B^T\Pi_3^T\Gamma_3-B^T\Pi_{3}\Gamma_{2}-B^T\Pi_{2}\Gamma_{2} }$

$\displaystyle{(45)\quad \frac{\partial {\rm tr}(\Pi A_{CL1}\Gamma)}{\partial H}=\Pi_{3}^T\Gamma_{1}C^T-\Pi_{3}^T\Gamma_{3}C^T+\Pi_{2}\Gamma_{3}^TC^T-\Pi_{2}\Gamma_{2}C^T }$

●準備３

$\displaystyle{(46)\quad \begin{array}{l} A_{CL1}^T\Pi\Gamma=(\Pi A_{CL1})^T\Gamma=\\ \left[\begin{array}{cc} A^T\Pi_{1}+C^TH^T\Pi_{3}^T&\\ -F^TB^T\Pi_{1}+(A-HC-BF)^T\Pi_{3}^T & \end{array}\right.\\ \left.\begin{array}{cc} A^T\Pi_{3}+C^TH^T\Pi_{2}\\ -F^TB^T\Pi_{3}+(A-HC-BF)^T\Pi_{2} \end{array}\right] \left[\begin{array}{cc} \Gamma_{1}&\Gamma_{3} \\ \Gamma_{3}^T&\Gamma_{2} \end{array}\right] \end{array} }$

$\displaystyle{(47)\quad \begin{array}{l} {\rm tr}(A_{CL1}^T\Pi\Gamma)\\ ={\rm tr}(A^T\Pi_{1}\Gamma_{1}+C^TH^T\Pi_{3}^T\Gamma_{1}+A^T\Pi_{3}\Gamma_{3}^T+C^TH^T\Pi_{2}\Gamma_{3}^T)\\ +{\rm tr}( -F^TB^T\Pi_{1}\Gamma_{3}+(A-HC-BF)^T\Pi_{3}^T\Gamma_{3} \\ -F^TB^T\Pi_{3}\Gamma_{2}+(A-HC-BF)^T\Pi_{2}\Gamma_{2}) \end{array} }$

$\displaystyle{(48)\quad \frac{\partial {\rm tr}(\Pi A_{CL1}\Gamma)}{\partial F}=-B^T\Pi_1\Gamma_{3}-B^T\Pi_3^T\Gamma_3-B^T\Pi_{3}\Gamma_{2}-B^T\Pi_{2}\Gamma_{2} }$

$\displaystyle{(49)\quad \frac{\partial {\rm tr}(\Pi A_{CL1}\Gamma)}{\partial H}=\Pi_{3}^T\Gamma_{1}C^T-\Pi_{3}^T\Gamma_{3}C^T+\Pi_{2}\Gamma_{3}^TC^T-\Pi_{2}\Gamma_{2}C^T }$

●準備４

$\displaystyle{(50)\quad Q_{CL1}^T\Gamma= \left[\begin{array}{cc} C^TQC & 0 \\ 0 & F^TRF \end{array}\right] \left[\begin{array}{cc} \Gamma_{1}&\Gamma_{3} \\ \Gamma_{3}^T&\Gamma_{2} \end{array}\right] }$

$\displaystyle{(51)\quad {\rm tr}(Q_{CL1}\Gamma)={\rm tr}(C^TQC\Gamma_{1}+F^TRF\Gamma_{2}) }$

$\displaystyle{(52)\quad \frac{\partial {\rm tr}(Q_{CL1}\Gamma)}{\partial F}=2RF\Gamma_{2} }$

● $\frac{\partial J'}{\partial F}=0$

$\displaystyle{(53)\quad \begin{array}{l} \frac{\partial J}{\partial F}=2RF\Gamma_{2}-2B^T(\Pi_1+\underbrace{\Pi_3^T}_{-\Pi_2})\underbrace{\Gamma_3}_{\Gamma_2}-2B^T\underbrace{(\Pi_{3}+\Pi_{2})}_{0}\Gamma_{2}=0\\ \qquad\Downarrow\Pi=\Pi_1-\Pi_2>0\\ 2(RF-B^T\Pi)\Gamma_{2}=0\\ \qquad\Downarrow\Gamma_2>0\\ \boxed{F=R^{-1}B^T\Pi}\\ \qquad\Downarrow\Pi(A-BF)+(A-BF)^T\Pi+C'^TQC'+F^TRF=0\\ \boxed{\Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C'^TQC'=0} \end{array} }$

● $\frac{\partial J'}{\partial H}=0$

$\displaystyle{(54)\quad \begin{array}{l} \frac{\partial J}{\partial H}=2\Pi_{2}HV+2\underbrace{\Pi_{3}^T}_{-\Pi_{2}}(\Gamma_{1}-\underbrace{\Gamma_{3}}_{\Gamma_{2}})C^T+2\Pi_{2}\underbrace{(\Gamma_{3}^T-\Gamma_{2})}_{0}C^T=0\\ \qquad\Downarrow\Gamma=\Gamma_1-\Gamma_2>0\\ 2\Pi_{2}(HV-\Gamma C^T)=0\\ \qquad\Downarrow\Pi_2>0\\ \boxed{H=(V^{-1}C\Gamma)^T}\\ \qquad\Downarrow(A-HC)\Gamma+\Gamma(A-HC)^T+B'WB'^T+HVH^T=0\\ \boxed{\Gamma A^T+A\Gamma-\Gamma C^TV^{-1}C\Gamma+B'WB'^T=0} \end{array} }$

補遺　上述の議論では、次についての検討が必要です。

検討事項１　(13)の妥当性
検討事項２　(24),(35)における仮定の妥当性
検討事項３　 $i=2$ の場合の導出
検討事項４　十分性の証明

LQ制御の応用

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

LQ制御の応用…Homework

[1]　倒立振子は制御なしには平衡状態を保てないので、制御技術の有効性を示すのに、よく研究されています。次図は様々な倒立振子を示しています。

図1 様々な倒立振子

これらの倒立振子について、非線形シミュレータを開発し、平衡状態まわりの挙動を表す線形状態方程式を得ていました。これらに基づいて、以下では、LQ制御系を設計します。

[2] CIP

平衡状態 $\theta^*=0$ 周りの挙動を表す線形状態方程式は次式で与えられます。

$\displaystyle{(2.1)\quad \frac{d}{dt}\left[\begin{array}{c} r(t)\\ \theta(t)\\ \dot{r}(t)\\ \dot{\theta}(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cccc} 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ 0 & -\frac{3gm}{4M+m} & 0 & 0\\ 0 & \frac{3(M+m)g}{(4M+m)\ell} & 0 & 0\\ \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} r(t)\\ \theta(t)\\ \dot{r}(t)\\ \dot{\theta}(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0\\ 0\\ \frac{4}{4M+m}\\ \frac{3}{(4M+m)\ell} \end{array}\right] }_{B} \underbrace{F(t)}_{u(t)} }$

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

$\displaystyle{(2.2)\quad u(t)= \underbrace{- \left[\begin{array}{cccc} f_1 & f_2 & f_3 & f_4 \end{array}\right] }_{-F} \underbrace{ \left[\begin{array}{c} r(t)\\ \theta(t)\\ \dot{r}(t)\\ \dot{\theta}(t) \end{array}\right] }_{x(t)} }$

を、評価関数

$\displaystyle{(2.3)\quad J=\int_0^\infty (\underbrace{\frac{1}{M_r^2}}_{q_1^2}r^2(t)+\underbrace{\frac{1}{M_\theta^2}}_{q_2^2}\theta^2(t)+\rho^2\underbrace{\frac{1}{M_F^2}}_{R}u^2(t))\,dt }$

$\displaystyle{(2.4)\quad |r(t)|<M_r,\ |\theta(t)|<M_\theta,\ |F(t)|<M_F }$

を最小にするように設計します。ここで、(2.4)は各変数が平衡状態周りで取り得る範囲を表しています。これらの逆数を重み係数に用いる方法をInverse Weighting法と呼ぶことがあります。 $\rho$ は状態変数と入力変数の間のバランスをとる役割を果たします。また、上の評価関数はあくまで一例であり、他にも候補があり、その選択は設計者に任せられています。

MATLAB

%cCIP_lq.m
%-----
 clear all, close all
 global mc m ell g th0
 mc=1; m=0.1; ell=0.2; g=9.8;
 th0=3/180*pi %input('th(0) = <0,180> ')/180*pi;
%-----
 A=[zeros(2,2) eye(2);zeros(2,4)];
 A(3,1)=0; 
 A(3,2)=-3*m*g/(m+4*mc); 
 A(4,1)=0; 
 A(4,2)=3*(m+mc)*g/((m+4*mc)*ell);  
 B=zeros(4,1);
 B(3)=4/(m+4*mc);
 B(4)=-3/((m+4*mc)*ell); 
%-----
 C=eye(2,4);
 Mr=0.5; Mth=3/180*pi; MF=1; rho=1;
 Q=diag([1/Mr^2,1/Mth^2]); R=rho^2*1/MF^2;
 [F,pcl]=opt(A,B,C,Q,R)
 C=diag([100,180/pi])*C;
 sys=ss(A-B*F,B,[C;-F],[]);
 t=0:0.05:5;
 x0=[0;th0;0;0];
 initial(sys,x0,t), grid
%-----
 xs=[0;0;0;0];
 us=0;
 sim("nCIP_lq")
%-----
%end

[3] CIP2

平衡状態 $\theta^*=0$ と平衡入力 $F^*=(M+m)g\sin\alpha$ 周りの挙動を表す線形状態方程式は次式で与えられます。

$(3.1)\quad \begin{array}{l} \displaystyle{\frac{d}{dt}\left[\begin{array}{c} r(t)\\ \theta(t)\\ \dot{r}(t)\\ \dot{\theta}(t) \end{array}\right] = \underbrace{\left[\begin{array}{cccc} 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ 0 & -\frac{6\cos\alpha mg}{8M+(5-3\cos2\alpha)m} & 0 & 0\\ 0 & \frac{6(M+m)g}{(8M+(5-3\cos2\alpha)m)\ell} & 0 & 0\\ \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} r(t)\\ \theta(t)\\ \dot{r}(t)\\ \dot{\theta}(t) \end{array}\right] }_{x(t)}}\\ \displaystyle{+ \underbrace{\left[\begin{array}{c} 0\\ 0\\ \frac{8}{8M+(5-3\cos2\alpha)m}\\ \frac{6\cos\alpha}{(8M+(5-3\cos2\alpha)m)\ell} \end{array}\right] }_{B} \underbrace{(F(t)-F^*)}_{u(t)}} \end{array}$

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

$\displaystyle{(3.2)\quad u(t)= \underbrace{- \left[\begin{array}{cccc} f_1 & f_2 & f_3 & f_4 \end{array}\right] }_{-F} \underbrace{ \left[\begin{array}{c} r(t)\\ \theta(t)\\ \dot{r}(t)\\ \dot{\theta}(t) \end{array}\right] }_{x(t)} }$

を、評価関数

$\displaystyle{(3.3)\quad J=\int_0^\infty (\underbrace{\frac{1}{M_r^2}}_{q_1^2}r^2(t)+\underbrace{\frac{1}{M_\theta^2}}_{q_2^2}\theta^2(t)+\rho^2\underbrace{\frac{1}{M_F^2}}_{R}u^2(t))\,dt }$

$\displaystyle{(3.4)\quad |r(t)|<M_r,\ |\theta(t)|<M_\theta,\ |F(t)-F^*|<M_F }$

を最小にするように設計します。

MATLAB

%cCIP2_lq.m
%-----
 clear all, close all
 global mc m ell g alpha th0
 mc=1; m=0.1; ell=0.2; g=9.8;
 alpha=3/180*pi;
 th0=3/180*pi %input('ths   = <0,180> ')/180*pi;
%-----
 A=[zeros(2,2) eye(2);zeros(2,4)];
 A(3,1)=0; 
 A(3,2)=-6*cos(alpha)*m*g/(8*mc+(5-3*cos(2*alpha))*m); 
 A(4,1)=0; 
 A(4,2)=6*(m+mc)*g/((8*mc+(5-3*cos(2*alpha))*m)*ell);  
 B=zeros(4,1);
 B(3)=8/(8*mc+(5-3*cos(2*alpha))*m);
 B(4)=-6*cos(alpha)/((8*mc+(5-3*cos(2*alpha))*m)*ell); 
%-----
 C=eye(2,4);
 Mr=0.5; Mth=3/180*pi; MF=1; rho=1;
 Q=diag([1/Mr^2,1/Mth^2]); R=rho^2*1/MF^2;
 [F,pcl]=opt(A,B,C,Q,R)
 C=diag([100,180/pi])*C;
 sys=ss(A-B*F,B,[C;-F],[]);
 t=0:0.05:5;
 x0=[0;th0;0;0];
 initial(sys,x0,t), grid
%-----
 xs=[0;0;0;0];
 us=(mc+m)*g*sin(alpha);
 sim("nCIP2_lq")
%-----
%end

●いまどれくらいレールが傾けられたか（ $\alpha$ の値）は未知とします。この場合 $F^*=(M+m)g\sin\alpha$ の設定ができないので、これを零として制御系のシミュレーションを実行してみます。そうすると安定化された倒立位置がずれてしまいます（平衡状態の移動）。これを防ぐにはLQI制御を適用する必要があります。

[4] AIP

平衡状態 $\theta_1^*=0,\ \theta_2^*=0$ と平衡入力 $\tau^*=-(m_1+2m_2)\ell_1g\sin\theta_1^*=0$ 周りの挙動を表す線形状態方程式は次式で与えられます。

$(4.1)\quad \begin{array}{l} \displaystyle{\frac{d}{dt}\left[\begin{array}{c} \theta_1(t)-\theta_1^*\\ \theta_2(t)\\ \dot{\theta}_1(t)\\ \dot{\theta}_2(t) \end{array}\right] =}\\ \displaystyle{\underbrace{\left[\begin{array}{cccc} 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ \frac{3(m_1+2m_2)g}{(4m_1+3m_2)\ell_1} & -\frac{9m_2g}{2(4m_1+3m_2)\ell_1} & 0 & 0\\ -\frac{9(m_1+2m_2)g}{2(4m_1+3m_2)\ell_1} & \frac{9m_2g}{(4m_1+3m_2)\ell_2} & 0 & 0\\ \end{array}\right] }_{A} \left[\begin{array}{c} \theta_1(t)-\theta_1^*\\ \theta_2(t)\\ \dot{\theta}_1(t)\\ \dot{\theta}_2(t) \end{array}\right]\\ \displaystyle{+ \underbrace{\left[\begin{array}{c} 0\\ 0\\ \frac{6}{2(4m_1+3m_2)\ell_1^2}\\ -\frac{9}{2(4m_1+3m_2)\ell_1\ell_2} \end{array}\right] }_{B} \underbrace{(\tau(t)-\tau^*)}_{u(t)}} \end{array}$

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

$\displaystyle{(4.2)\quad u(t)= \underbrace{- \left[\begin{array}{cccc} f_1 & f_2 & f_3 & f_4 \end{array}\right] }_{-F} \underbrace{ \left[\begin{array}{c} \theta_1(t)-\theta_1^*\\ \theta_2(t)\\ \dot{\theta}_1(t)\\ \dot{\theta}_2(t) \end{array}\right] }_{x(t)} }$

を、評価関数

$\displaystyle{(4.3)\quad J=\int_0^\infty (\underbrace{\frac{1}{M_{\theta_1}^2}}_{q_1^2}(\theta_1(t)-\theta_1^*)^2(t)+\underbrace{\frac{1}{M_{\theta_2}^2}}_{q_2^2}\theta^2(t)+\rho^2\underbrace{\frac{1}{M_\tau^2}}_{R}u^2(t))\,dt }$

$\displaystyle{(4.4)\quad |\theta_1(t)-\theta_1^*|<M_{\theta_1},\ |\theta_2(t)|<M_{\theta_2},\ |\tau(t)-\tau^*|<M_\tau }$

を最小にするように設計します。

MATLAB

%cAIP_lq.m
%-----
 clear all, close all
 global m1 m2 ell1 ell2  g th10 th20
 m1=0.1; m2=0.2; ell1=0.1; ell2=0.2; g=9.8;
 th10=0/180*pi %input('th1(0) = <0,180> ')/180*pi;
 th20=-3/180*pi %input('th2(0) = <0,180> ')/180*pi;  
%----- th10=0 における線形モデル
 A=[zeros(2,2) eye(2);zeros(2,4)];
 A(3,1)=3*(2*m2+m1)*g/((3*m2+4*m1)*ell1); 
 A(3,2)=-9*m2*g/(2*(3*m2+4*m1)*ell1); 
 A(4,1)=-9*(2*m2+m1)*g/(2*(3*m2+4*m1)*ell2); 
 A(4,2)=3*(3*m2+m1)*g/((3*m2+4*m1)*ell2);  
 B=zeros(4,1);
 B(3)=3/((3*m2+4*m1)*ell1^2);
 B(4)=-9/(2*(3*m2+4*m1)*ell1*ell2);
%-----
 C=eye(2,4);
 Mth1=3/180*pi; Mth2=3/180*pi; MF=1; rho=1;
 Q=diag([1/Mth1^2,1/Mth2^2]); R=rho^2*1/MF^2;
 [F,pcl]=opt(A,B,C,Q,R)
 C=diag([180/pi,180/pi])*C;
 sys=ss(A-B*F,B,[C;-F],[]);
 t=0:0.01:1;
 x0=[th10;th20;0;0];
 initial(sys,x0,t), grid
%-----
 xs=[0;0;0;0];
 us=-(m1+2*m2)*ell1*g*sin(th10)
 sim("nAIP_lq")
%-----
%end

●平衡状態 $\theta_1^*\ne 0,\ \theta_2^*=0$ に移動させるためには、LQI制御を適用する必要があります。

[5] PIP

平衡状態 $\theta_1^*=0,\ \theta_2^*=0$ 周りの挙動を表す線形状態方程式は次式で与えられます。

$(5.1)\quad \begin{array}{l} \displaystyle{\frac{d}{dt}\left[\begin{array}{c} r(t)\\ \theta_1(t)\\ \theta_2(t)\\ \dot{r}(t)\\ \dot{\theta}_1(t)\\ \dot{\theta}_2(t) \end{array}\right] =}\\ \displaystyle{\underbrace{\left[\begin{array}{cccccc} 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 1\\ 0 & -\frac{3m_1g}{4M+m_1+m_2} & -\frac{3m_2g}{4M+m_1+m_2} & 0 & 0 & 0\\ 0 & \frac{3(4m+4m_1+m_2)g}{4(4M+m_1+m_2)\ell_1} & \frac{9m_2g}{4(4M+m_1+m_2)\ell_1} & 0 & 0& 0\\ 0 & \frac{9gm_1}{4(4M+m_1+m_2)\ell_2} & \frac{3(4m+m_1+4m_2)g}{4(4M+m_1+m_2)\ell_2} & 0 & 0& 0\\ \end{array}\right]}_{A} \left[\begin{array}{c} r(t)\\ \theta_1(t)\\ \theta_2(t)\\ \dot{r}(t)\\ \dot{\theta}_1(t)\\ \dot{\theta}_2(t) \end{array}\right]}\\ \displaystyle{+ \underbrace{ \left[\begin{array}{c} 0\\ 0\\ 0\\ \frac{4}{4M+m_1+m_2}\\ -\frac{3}{(4M+m_1+m_2)\ell_1}\\ -\frac{3}{(4M+m_1+m_2)\ell_2} \end{array}\right] }_{B} \underbrace{F(t)}_{u(t)} \end{array}$

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

$\displaystyle{(5.2)\quad u(t)= \underbrace{- \left[\begin{array}{cccccc} f_1 & f_2 & f_3 & f_4 & f_5 & f_6 \end{array}\right] }_{-F} \underbrace{ \left[\begin{array}{c} r(t)\\ \theta_1(t)\\ \theta_2(t)\\ \dot{r}(t)\\ \dot{\theta}_1(t)\\ \dot{\theta}_2(t) \end{array}\right] }_{x(t)} }$

を、評価関数

$\displaystyle{(5.3)\quad J=\int_0^\infty (\underbrace{\frac{1}{M_r^2}}_{q_1^2}r^2(t)+\underbrace{\frac{1}{M_{\theta_1}^2}}_{q_2^2}\theta_1^2(t)+\underbrace{\frac{1}{M_{\theta_2}^2}}_{q_2^2}\theta^2(t)+\rho^2\underbrace{\frac{1}{M_F^2}}_{R}u^2(t))\,dt }$

$\displaystyle{(5.4)\quad |r(t)|<M_r,\ |\theta_1(t)|<M_{\theta_1},\ |\theta_2(t)|<M_{\theta_2},\ |F(t)|<M_F }$

を最小にするように設計します。

MATLAB

%cPIP_lq.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; 
 th10=0/180*pi %input('th1(0) = <0,180> ')/180*pi;
 th20=-1/180*pi %input('th2(0) = <0,180> ')/180*pi;  
%-----
 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); 
%------
 C=eye(3,6);
 Mr=0.5; Mth1=3/180*pi; Mth2=3/180*pi; MF=1; rho=1;
 Q=diag([1/Mr^2,1/Mth1^2,1/Mth2^2]); R=rho^2*1/MF^2;
 [F,pcl]=opt(A,B,C,Q,R);
 C=diag([100,180/pi,180/pi])*C;
 sys=ss(A-B*F,B,[C;-F],[]);
 t=0:0.05:5;
 x0=[0;th10;th20;0;0;0];
 initial(sys,x0,t), grid
%-----
 xs=[0;0;0;0;0;0];
 us=0;
 sim("nPIP_lq")   
%-----
%eof

[6] DIP

平衡状態 $\theta_1^*=0,\ \theta_2^*=0$ 周りの挙動を表す線形状態方程式は次式で与えられます。

$(6.1)\quad \begin{array}{l} \displaystyle{ \frac{d}{dt}\left[\begin{array}{c} r(t)\\ \theta_1(t)\\ \theta_2(t)\\ \dot{r}(t)\\ \dot{\theta}_1(t)\\ \dot{\theta}_2(t) \end{array}\right] =}\\ \displaystyle{ \left[\begin{array}{cccccc} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & -\frac{3m_1g}{Mm_1+m_1^2+(3M+m1)m_2} & -\frac{3m_2g}{Mm_1+m_1^2+(3M+m1)m_2} \\ 0 & \frac{3(4Mm_1+4m_1^2+3m_2^2+3(18M+13m_1)m_2)g}{(4Mm_1+m_1^2+(3M+m1)m_2)\ell_1} & \frac{9(2M+m_1)m_2g}{(4Mm_1+m_1^2+(3M+m1)m_2)\ell_1} \\ 0 & \frac{9(2Mm_1+m_1^2+3(2M+m_1)m_2)g}{(4Mm_1+m_1^2+(3M+m_1)m_2)\ell_2} & \frac{3(4Mm_1+m_1^2+12(3M+m_1)m_2)g}{(4Mm_1+m_1^2+(3M+m_1)m_2)\ell_2} \\ \end{array}\right.}\\ \displaystyle{\left.\begin{array}{cccccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\\ 0 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0\\ \end{array}\right] \left[\begin{array}{c} r(t)\\ \theta_1(t)\\ \theta_2(t)\\ \dot{r}(t)\\ \dot{\theta}_1(t)\\ \dot{\theta}_2(t) \end{array}\right]}\\ \displaystyle{+ \underbrace{\left[\begin{array}{c} 0\\ 0\\ 0\\ \frac{4m_1+3m_2}{4Mm_1+m_1^2+(3M+m1)m_2}\\ -\frac{3(2m_1+m2)}{2Mm_1+m_1^2+(3M+m1)m_2}\\ \frac{3m_1}{2(4Mm_1+m_1^2+(3M+m_1)m_2)\ell_2} \end{array}\right] }_{B} \underbrace{F(t)}_{u(t)}} \end{array}$

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

$\displaystyle{(6.2)\quad u(t)= \underbrace{- \left[\begin{array}{cccccc} f_1 & f_2 & f_3 & f_4 & f_5 & f_6 \end{array}\right] }_{-F} \underbrace{ \left[\begin{array}{c} r(t)\\ \theta_1(t)\\ \theta_2(t)\\ \dot{r}(t)\\ \dot{\theta}_1(t)\\ \dot{\theta}_2(t) \end{array}\right] }_{x(t)} }$

を、評価関数

$\displaystyle{(6.3)\quad J=\int_0^\infty (\underbrace{\frac{1}{M_r^2}}_{q_1^2}r^2(t)+\underbrace{\frac{1}{M_{\theta_1}^2}}_{q_2^2}\theta_1^2(t)+\underbrace{\frac{1}{M_{\theta_2}^2}}_{q_2^2}\theta^2(t)+\rho^2\underbrace{\frac{1}{M_F^2}}_{R}u^2(t))\,dt }$

$\displaystyle{(6.4)\quad |r(t)|<M_r,\ |\theta_1(t)|<M_{\theta_1},\ |\theta_2(t)|<M_{\theta_2},\ |F(t)|<M_F }$

を最小にするように設計します。

MATLAB

%cDIP_lq.m
%-----
 clear all, close all
 global mc m1 m2 ell1 ell2 g th10 th20
 mc=1; m=0.1; m1=m; m2=m; ell=0.15; ell1=ell; ell2=ell; g=9.8;
 th10=1/180*pi %input('th1(0) = <0,180> ')/180*pi;
 th20=0 %input('th2(0) = <0,180> ')/180*pi;   
 %-----
 A=[zeros(3,3) eye(3);zeros(3,6)];
 A(4,1)=0; 
 A(4,2)=-18*m*g/(2*m+7*mc); 
 A(4,3)=3*m*g/(4*m+14*mc);  
 A(5,1)=0; 
 A(5,2)=(15*m+12*mc)*g/((2*m+7*mc)*ell);  
 A(5,3)=-(9*m+18*mc)*g/((8*m+28*mc)*ell);   
 A(6,1)=0; 
 A(6,2)=-(9*m+18*mc)*g/((2*m+7*mc)*ell);  
 A(6,3)=(15*m+48*mc)*g/((8*m+28*mc)*ell);   
 B=zeros(6,1);
 B(4)=7/(2*m+7*mc);
 B(5)=-9/((4*m+14*mc)*ell);
 B(6)=3/((4*m+14*mc)*ell);
%------
 C=eye(3,6);
 Mr=0.5; Mth1=3/180*pi; Mth2=3/180*pi; MF=1; rho=1;
 Q=diag([1/Mr^2,1/Mth1^2,1/Mth2^2]); R=rho^2*1/MF^2;
 [F,pcl]=opt(A,B,C,Q,R);
 C=diag([100,180/pi,180/pi])*C;
 sys=ss(A-B*F,B,[C;-F],[]);
 t=0:0.05:5;
 x0=[0;th10;th20;0;0;0];
 initial(sys,x0,t), grid
%-----
 xs=[0;0;0;0;0;0];
 us=0;
 sim("nDIP_lq")   
%-----
%eof

LQ制御

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

LQ制御…Homework

[1]　可制御な制御対象

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

を安定化する状態フィードバック

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

の決定法を考えます。一つの方法は，閉ループ系

$\displaystyle{(3)\quad %\left\{\begin{array}{l} \dot{x}(t)=(A-BF)x(t) \\ %z(t)=Cx(t) %\end{array}\right. }$

の時間応答に関する評価規範として，２次形式評価関数

$(4)\quad \boxed{\begin{array}{l} \displaystyle{J=\int_0^\infty (z^T(t)Qz(t)+u^T(t)Ru(t))\,dt\quad (Q>0,R>0)}\\ \displaystyle{z(t)=Cx(t)\quad (z(t)\in{\bf R}^q)} \end{array}}$

を設定し，これを最小化する問題を解くことです。ただし、 $(A,C)$ は可観測対とします。これによる状態フィードバックのゲイン行列 $F$ は，リッカチ方程式

$\displaystyle{(5)\quad {\boxed{\Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0}} }$

の解 $\Pi>0$ を用いて，次式で与えられます。

$\displaystyle{(6)\quad {\boxed{F=R^{-1}B^T\Pi}} }$

この証明はNoteに示しています。

[2]　以下では、代表的な２次系に対して、評価関数を設定して、ゲイン行列 $F$ を求めてみます。その際、リッカチ方程式は４つの解候補を持ちますが、 $\Pi>0$ の条件、すなわち

$\displaystyle{(7)\quad \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]>0\Leftrightarrow \pi_1>0,\ \pi_1\pi_2-\pi_3^2>0 }$

を用いて（ $\pi_2>0$ ）、解を１つに絞ることに注意してください。

●いま２次系（２重積分器）

$\displaystyle{(8)\quad \boxed{\left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] + \left[\begin{array}{c} 0 \\ 1 \end{array}\right] u(t)} }$

を安定化する状態フィードバック $u(t)=-f_1x_1(t)-f_2x_2(t)$ を，評価関数

$\displaystyle{(9)\quad J=\int_0^\infty (x_1^2(t)+x_2^2(t)+u^2(t))\,dt }$

を最小にするように求めると

$\displaystyle{(10)\quad \boxed{ f_1=1,\quad f_2=\sqrt{3}\right) }}$

となります。実際、リッカチ方程式

$\displaystyle{(11)\quad \begin{array}{l} \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right] + \left[\begin{array}{cc} 0 & 0 \\ 1 & 0 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]\\ - \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{c} 0 \\ 1 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] + \left[\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right]\\ = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \end{array} }$

を要素ごとに整理して

$\displaystyle{(12)\quad \left\{ \begin{array}{l} -\pi_3^2+1=0\\ \pi_1-\pi_2\pi_3=0\\ 2\pi_3-\pi_2^2+1=0 \end{array} \right. }$

を得る。これは，まず第1式より $\pi_3$ が２つ，つぎに第3式より $\pi_2$ が２つ，さらに第2式より $\pi_1$ が１つ定まり，つぎのように４組の解をもつ。すなわち

$\displaystyle{(13)\quad \left\{\begin{array}{llll} \pi_3= 1 & \Rightarrow \left\{\begin{array}{l} \pi_2= \sqrt{3}\\ \pi_2=-\sqrt{3} \end{array}\right. & \left.\begin{array}{l} \Rightarrow\pi_1= \sqrt{3}\\ \Rightarrow\pi_1=-\sqrt{3} \end{array}\right. & \left.\begin{array}{ll} (*) &○\\ (**) &× \end{array}\right. \\ \pi_3=-1 & \Rightarrow \left\{\begin{array}{llll} \pi_2= j \\ \pi_2=-j \end{array}\right. & \hspace{-3mm} \left.\begin{array}{l} \Rightarrow\pi_1=-j\\ \Rightarrow\pi_1=j \end{array}\right. & \left.\begin{array}{ll} (***) &×\\ (****) &× \end{array}\right. \end{array}\right. }$

ここで，(*)だけが， $\pi_1,\pi_2>0,\ \pi_1\pi_2-\pi_3^2>0$ を満たします。したがって

$\displaystyle{(14)\quad \left[\begin{array}{cc} f_1 & f_2 \end{array}\right]= \left[\begin{array}{cc} 0 & 1 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]= \left[\begin{array}{cc} \pi_3 & \pi_2 \end{array}\right]= \left[\begin{array}{cc} 1 & \sqrt{3} \end{array}\right] }$

例題1　２次系（無定位系）

$\displaystyle{(15)\quad \boxed{\left\{\begin{array}{l} %\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 \\ 0 & -2\zeta\omega_n \end{array}\right] %}_{A} %\underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] + %\underbrace{ \left[\begin{array}{c} 0 \\ \omega_n^2 \end{array}\right] %}_{B} \\ 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] \end{array}\right.} }$

を安定化する状態フィードバック $u(t)=-f_1x_1(t)-f_2x_2(t)$ を，評価関数

$\displaystyle{(16)\quad J=\int_0^\infty (q^2y^2(t)+r^2u^2(t))\,dt }$

を最小にするように求めると

$\displaystyle{(17)\quad \boxed{ f_1=\frac{q}{r},\quad f_2=\frac{2}{\omega_n}\left(-\zeta+\sqrt{\zeta^2+\frac{1}{2}\frac{q}{r}}\right) }}$

また

$\displaystyle{(18)\quad J=\int_0^\infty (q_1^2x_1^2(t)+q_2^2x_2^2(t)+r^2u^2(t))\,dt }$

を最小にするように求めると

$\displaystyle{(19)\quad \boxed{\left\{\begin{array}{l} f_1=\frac{q_1}{r} \\ f_2=\frac{2}{\omega_n}\left(-\zeta+\sqrt{\zeta^2+\frac{1}{2}\frac{q_1}{r}+ \left(\frac{q_2}{r}\right)^2\left(\frac{\omega_n}{2}\right)^2 }\right) \end{array}\right.} }$

例題2　２次系

$\displaystyle{(20)\quad \boxed{\left\{\begin{array}{l} %\underbrace{ \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right] %}_{\dot x} = %\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} + %\underbrace{ \left[\begin{array}{c} 0 \\ \omega_n^2 \end{array}\right] %}_{B} u(t) \\ 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] \end{array}\right.} }$

を安定化する状態フィードバック $u(t)=-f_1x_1(t)-f_2x_2(t)$ を，評価関数

$\displaystyle{(16)\quad J=\int_0^\infty (q^2y^2(t)+r^2u^2(t))\,dt }$

を最小にするように求めると

$\displaystyle{(21)\quad \boxed{\left\{\begin{array}{l} f_1=-1+\sqrt{1+\left(\frac{q}{r}\right)^2}\\ f_2=-\frac{2}{\omega_n}\left(-\zeta+\sqrt{\zeta^2-\frac{1}{2}+\frac{1}{2}\sqrt{1+\left(\frac{q}{r}\right)^2}}\right) \end{array}\right.} }$

また

$\displaystyle{(18)\quad J=\int_0^\infty (q_1^2x_1^2(t)+q_2^2x_2^2(t)+r^2u^2(t))\,dt }$

を最小にするように求めると

$\displaystyle{(22)\quad \boxed{\left\{\begin{array}{l} f_1=-1+\sqrt{1+\left(\frac{q_1}{r}\right)^2}\\ f_2=-\frac{2}{\omega_n}\left(-\zeta+\sqrt{\zeta^2-\frac{1}{2}+\frac{1}{2}\sqrt{1+\left(\frac{q_1}{r}\right)^2}+\left(\frac{q_2}{r}\right)^2\left(\frac{\omega_n}{2}\right)^2}\right) \end{array}\right.} }$

[3]　計算機でリッカチ方程式解くには、ハミルトン行列と呼ばれる

$\displaystyle{(23)\quad \boxed{M=\left[\begin{array}{cc} A & -BR^{-1}B^T \\ C^TQC & -A^T \end{array}\right]} }$

を考えます。このハミルトン行列の固有値は、実軸に対称ばかりでなく、虚軸にも対称となるという性質を持っています。これらのうち安定な固有値と対応する固有ベクトルを、次のように求めます。

$\displaystyle{(24)\quad \begin{array}{l} \underbrace{ \left[\begin{array}{cc} A & -BR^{-1}B^T \\ -C^TQC & -A^T \end{array}\right]}_{M(2n\times 2n)} \underbrace{ \left[\begin{array}{c} V_1 \\ V_2 \end{array}\right]}_{V^-(2n\times n)}\\ = \underbrace{ \left[\begin{array}{c} V_1 \\ V_2 \end{array}\right]}_{V^-(2n\times n)} \underbrace{ {\rm diag}\{\lambda_1,\cdots,\lambda_n\} }_{\Lambda^-(n\times n)} \end{array} }$

これから

$\displaystyle{(25)\quad \boxed{\Pi=V_2V_1^{-1}} }$

のように求められます。

●ハミルトン行列を経由してリッカチ方程式解くためのコードは次のようになります。

MATLAB

%opt.m
%---------------------------------------
 function [F,p]=opt(A,B,C,Q,R)
%---------------------------------------
 n=size(A,1);
 w2=R\B';
 [v,p]=eig([A -B*w2;-C'*Q*C -A']);
 p=diag(p);
 [dummy,index]=sort(real(p));
 p=p(index(1:n));
 x=v(1:n,index(1:n));
 y=v(n+1:n+n,index(1:n));
 X=real(y/x);
 F=w2*X;
%---------------------------------------
%eof

SCILAB

function [F,p]=opt(A,B,C,Q,R)
 w2=R\B';
 [v,p]=spec([A -B*w2;-C'*Q*C -A']); p=diag(p);
 [dummy,index]=gsort(real(p));
 n=size(A,1); j=index(n+1:n+n);
 p=p(j); V1=v(1:n,j); V2=v(n+1:n+n,j);
 X=real(V2/V1); F=w2*X;
endfunction

MATLABでは、関数lqrがリッカチ方程式解くために準備されています。

演習A52…Flipped Classroom

$\displaystyle{ (26)\ %\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 \\ 0 & -1 \end{array}\right] %}_{A} %\underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] %}_{x(t)} + %\underbrace{ \left[\begin{array}{c} 0 \\ 1 \end{array}\right] %}_{B} u(t) }$

$\displaystyle{ (27)\ %\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 \\ -1 & 0 \end{array}\right] %}_{A} %\underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] %}_{x(t)} + %\underbrace{ \left[\begin{array}{c} 0 \\ 1 \end{array}\right] %}_{B} u(t) }$

Note A52-1 行列による微分　

いま、任意の行列 $X$ の $(i,j)$ 要素を $x_{ij}=[X]_{ij}$ で表すとき、スカラ関数 $f$ を行列変数 $X$ の各要素で微分して得られる行列を $\displaystyle{[\frac{\partial f}{\partial X}]_{ij}=\frac{\partial f}{\partial x_{ij}}}$ で定義します。このとき、行列のトレースについて、次が成り立ちます。

$\displaystyle{(1)\ {\rm tr}AB={\rm tr}BA}$
$\displaystyle{(2)\ \frac{\partial}{\partial X}{\rm tr}AXB=A^TB^T}$
$\displaystyle{(3)\ \frac{\partial}{\partial X}{\rm tr}AX^TB=BA}$
$\displaystyle{(4)\ \frac{\partial}{\partial X}{\rm tr}AXBX^T=A^TXB^T+AXB}$

実際、

$\displaystyle{(1)\ \sum_{i}[AB]_{ii}=\sum_{i}\sum_{j}a_{ij}b_{ji} =\sum_{j}\sum_{i}b_{ji}a_{ij}=\sum_{j}[BA]_{jj}}$

$\displaystyle{(2)\ [\frac{\partial}{\partial X}{\rm tr}AXB]_{ij} =\frac{\partial}{\partial x_{ij}}\sum_{k}[AXB]_{kk} =\frac{\partial}{\partial x_{ij}}\sum_{k}\sum_{i,j}a_{ki}x_{ij}b_{jk}}$
$\displaystyle{=\sum_{k}b_{jk}a_{ki}=[BA]_{ji}=[A^TB^T]_{ij}}$

$\displaystyle{(3)\ [\frac{\partial}{\partial X}{\rm tr}AX^TB]_{ij} =\frac{\partial}{\partial x_{ij}}\sum_{k}[AX^TB]_{kk} =\frac{\partial}{\partial x_{ij}}\sum_{k}\sum_{i,j}a_{ki}x_{ji}b_{jk}}$
$\displaystyle{=\sum_{k}b_{ik}a_{kj}=[BA]_{ij}}$

$\displaystyle{(4)\ \frac{\partial}{\partial X}{\rm tr}AXBX^T =\frac{\partial}{\partial X}{\rm tr}(AXB)X^T +\frac{\partial}{\partial X}{\rm tr}AX(BX^T)}$
$\displaystyle{=AXB+A^TXB^T}$

Note A52-2 LQ制御問題の解法　

可制御かつ可観測な $n$ 次系

$\displaystyle{(1)\quad \dot{x}(t)=Ax(t)+Bu(t),\ y(t)=Cx(t) }$

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

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

による閉ループ系

$\displaystyle{(3)\quad \begin{array}{l} \dot{x}(t)=(A-BF)x(t),\ z(t)=Cx(t)\\ A_F=A-BF:\ stable\ matrix \end{array} }$

に対して、評価関数

$\displaystyle{(4)\quad J=\int_0^\infty(z^T(t)Qz(t)+u^T(t)Ru(t))\,dt }$

を最小化するように $F$ を決める問題を考えます。

閉ループ系における状態の振る舞いは次式で与えられます。

$\displaystyle{(5)\quad x(t)=\exp({A_Ft})x(0) }$

ここで、１次系の場合は初期状態は $x(0)\ne0$ であればよかったのですが、一般の場合はインパルス応答となるように $B$ の列ベクトル $B^{(i)}$ を考えます。各インパルス応答

$\displaystyle{(6)\quad x^{(i)}(t)=\exp({A_Ft})B^{(i)} }$

に対する評価関数 $J$ の総和は

$(7)\quad \begin{array}{l} \displaystyle{\sum\int_0^\infty x^{(i)}\,^T(t)(C^TQC+F^TRF)x^{(i)}(t)\,dt}\\ \displaystyle{=\sum B^{(i)}\,^T\underbrace{\int_0^\infty \exp(A_F^Tt)(C^TQC+F^TRF)\exp(A_Ft)\,dt}_{\Pi}B^{(i)}}\\ \displaystyle{=\sum B^{(i)}\,^T\Pi B^{(i)}={\rm tr}\ \Pi BB^T} \end{array} }$

と書けます。いま $\Gamma$ をラグランジュの未定定数として

$\displaystyle{(8)\quad J'={\rm tr}\ \Pi x(0)x^T(0)+{\rm tr}\ ((\Pi A_F+A_F^T\Pi+C^TQC+F^TRF)\Gamma) }$

を最小化する問題を考えます。ここで、制約条件は、リャプノフ方程式と呼ばれる

$\displaystyle{(9)\quad \Pi A_F+A_F^T\Pi+C^TQC+F^TRF=0 }$

ですが、 $\Pi>0$ は $A_F$ が安定行列を意味することに注意します。

そこで、必要条件として次を得ます。

$\displaystyle{(11)\quad \frac{\partial J'}{\partial\Pi}&=&BB^T+\Gamma A_F^T+A_F\Gamma=0\Rightarrow \Gamma>0 }$
$\displaystyle{(12)\quad \frac{\partial J'}{\partial F}&=&2(RF-B^T\Pi)\Gamma=0\Rightarrow F=R^{-1}B^T\Pi }$
$\displaystyle{(13)\quad \frac{\partial J'}{\partial\Gamma}&=&\Pi A_F+A_F^T\Pi+C^TQC+F^TRF=0 }$

ここで、第２式から得られる $F=R^{-1}B^T\Pi$ を第３式に代入して

$\displaystyle{(14)\quad \begin{array}{l} \Pi (A-BR^{-1}B^T\Pi)+(A-BR^{-1}B^T\Pi)^T\Pi+C^TQ\nonumber\\ +(R^{-1}B^T\Pi)^TRR^{-1}B^T\Pi=0 \end{array} }$

すなわち、リッカチ方程式と呼ばれる $\Pi$ の行列方程式

$\displaystyle{(15)\quad \Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0 }$

を得ます。これから $\Pi>0$ を求めて、 $F$ は

$\displaystyle{(16)\quad F=R^{-1}B^T\Pi }$

のように得られます。このような制御方式をLQ制御と呼びます。

一方、十分性の議論は次のように行われます。まず、被積分項は次のように表すことができます。

$\displaystyle{(17)\quad \begin{array}{lll} &&y^TQy+u^TRu=(u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x)\nonumber\\ &&-\underbrace{\frac{d}{dt}x^T\Pi x}_{2x^T\Pi\dot{x}} \end{array} }$

実際、右辺に $\dot{x}=Ax+Bu$ を代入し、リッカチ方程式を用いると

$\displaystyle{(18)\quad \begin{array}{lll} &&{\rm RHS}=u^TRu+2x^T\Pi Bu+x^T\Pi BR^{-1}B^T\Pi x-2x^T\Pi(Ax+Bu)\nonumber\\ &&=u^TRu+x^T(\Pi A+A^T\Pi+C^TQC)x-2x^T\Pi Ax={\rm LHS} \end{array} }$

したがって、上記の両辺を積分して

$\displaystyle{(19)\quad J=\int_0^\infty (u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x)\,dt-\left[x^T\Pi x\right]_0^\infty }$

を得ます。ここで、 $x(t)=\exp(A_Ft)x(0)\rightarrow 0\ (t\rightarrow\infty)$ を前提とするので

$\displaystyle{(20)\quad J=\int_0^\infty (u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x)\,dt+x^T(0)\Pi x(0) }$

を得ます。これから $u=-R^{-1}B^T\Pi x$ が評価関数を最小化することが分かります。

Note A53-3a 例題１(17)の導出　
　
リッカチ方程式

$\displaystyle{(1)\quad \begin{array}{l} \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ 0 & -2\zeta\omega_n \end{array}\right] + \left[\begin{array}{cc} 0 & 0 \\ 1 & -2\zeta\omega_n \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]\\ - \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{c} 0 \\ \omega_n^2 \end{array}\right] r^{-2} \left[\begin{array}{cc} 0 & \omega_n^2 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]\\ + \left[\begin{array}{c} 1 \\ 0 \end{array}\right] q^{2} \left[\begin{array}{cc} 1 & 0 \end{array}\right] % \left[\begin{array}{cc} % q_1^2 & 0 \\ % 0 & q_2^2 % \end{array}\right] = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \end{array}} }$

を要素ごとに整理して

$\displaystyle{(2)\quad \left\{ \begin{array}{l} -r^{-2}\omega_n^4\pi_3^2+q^2=0\\ \pi_1-2\zeta\omega_n\pi_3-r^{-2}\omega_n^4\pi_2\pi_3=0\\ 2\pi_3-4\zeta\omega_n\pi_2-r^{-2}\omega_n^4\pi_2^2=0 \end{array}\right.} }$

を得る。まず，第１式より $\pi_3$ が

$\displaystyle{(3)\quad \pi_3=\pm r\omega_n^{-2}q }$

と求まる。つぎに，第３式より $\pi_2$ は

$\displaystyle{(4)\quad \pi_2= %\frac{-b\pm\sqrt{b^2-ac}}{a} r^2\omega_n^{-4}(-2\zeta\omega_n\pm\sqrt{(2\zeta\omega_n)^2\pm 2r^{-2}\omega_n^4r\omega_n^{-2}q}) }$

となるが， $\pi_2>0$ より

$\displaystyle{(5)\quad \begin{array}{l} \pi_2=r^2\omega_n^{-3}(-2\zeta+\sqrt{4\zeta^2+ 2r^{-1}q})\\ \pi_3=r\omega_n^{-2}q \end{array} }$

でなければならない。さらに，第２式より $\pi_1$ は

$\displaystyle{(6)\quad \begin{array}{l} \pi_1=(2\zeta\omega_n+r^{-2}\omega_n^4r^2\omega_n^{-3}(-2\zeta+\sqrt{4\zeta^2+ 2r^{-1}q}))r\omega_n^{-2}q\\ =r\omega_n^{-1}q\sqrt{4\zeta^2+ 2r^{-1}q \end{array} }$

となる（ $\pi_1>0$ ）。すなわち(1)の解として

$\displaystyle{(7)\quad \left\{\begin{array}{l} \pi_1=r\omega_n^{-1}q\sqrt{2r^{-1}q+4\zeta^2}}\\ \pi_2=r^2\omega_n^{-3}(-2\zeta+\sqrt{2r^{-1}q+4\zeta^2})}\\ \pi_3=r\omega_n^{-2}q} \end{array}\right.} }$

を得ます（このとき $\pi_1\pi_2-\pi_3^2>0$ も満足されます）。したがって

$\displaystyle{(8)\quad \begin{array}{l} \left[\begin{array}{cc} f_1 & f_2 \end{array}\right]=r^{-2} \left[\begin{array}{cc} 0 & \omega_n^2 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]=r^{-2}\omega_n^2 \left[\begin{array}{cc} \pi_3 & \pi_2 \end{array}\right]\\ = \left[\begin{array}{cc} \frac{q}{r} & \frac{2}{\omega_n}\left(-\zeta+\sqrt{\zeta^2+\frac{1}{2}\frac{q}{r}}\right) \end{array}\right] \end{array}} }$

Note A53-3b 例題１(19)の導出

リッカチ方程式

$\displaystyle{(1)\quad \begin{array}{l} \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ 0 & -2\zeta\omega_n \end{array}\right] + \left[\begin{array}{cc} 0 & 0 \\ 1 & -2\zeta\omega_n \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]\\ - \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{c} 0 \\ \omega_n^2 \end{array}\right] r^{-2} \left[\begin{array}{cc} 0 & \omega_n^2 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]\\ + \left[\begin{array}{cc} q_1^2 & 0 \\ 0 & q_2^2 \end{array}\right] = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \end{array}} }$

を要素ごとに整理して、上と同様にして導出されます。

Note A52-3c 例題２(21)の導出　

リッカチ方程式

$\displaystyle{(1)\quad \begin{array}{ll} \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] + \left[\begin{array}{cc} 0 & -\omega_n^2 \\ 1 & -2\zeta\omega_n \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]\\ - \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{c} 0 \\ \omega_n^2 \end{array}\right] r^{-2} \left[\begin{array}{cc} 0 & \omega_n^2 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]\\ + \left[\begin{array}{c} 1 \\ 0 \end{array}\right] q^2 \left[\begin{array}{cc} 1 & 0 \end{array}\right] = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \end{array}}$

を要素ごとに整理して

$\displaystyle{(2)\quad \left\{\begin{array}{l} -2\omega_n^2\pi_3-r^{-2}\omega_n^4\pi_3^2+q^2=0 \\ \pi_1-2\zeta\omega_n\pi_3-\omega_n^2\pi_2-r^{-2}\omega_n^4\pi_2\pi_3=0 \\ 2\pi_3-4\zeta\omega_n\pi_2-r^{-2}\omega_n^4\pi_2^2=0 \end{array}\right.$

を得る。まず，第１式より $\pi_3$ が

$\displaystyle{(3)\quad \pi_3=r^2\omega_n^{-2}(-1\pm\sqrt{1+r^{-2}q^2}) }$

と求まる。つぎに，第３式より $\pi_2$ は

$\displaystyle{(4)\quad \begin{array}{l} \pi_2= %\frac{-b\pm\sqrt{b^2-ac}}{a} a=r^{-2}\omega_n^4, b=2\zeta\omega_n, c=-2\pi_3 r^2\omega_n^{-4}(-2\zeta\omega_n\pm\sqrt{(2\zeta\omega_n)^2\pm 2r^{-2}\omega_n^4r^2\omega_n^{-2}(-1\pm\sqrt{1+r^{-2}q^2})})\\ =r^2\omega_n^{-3}(-2\zeta\pm\sqrt{4\zeta^2\pm 2(-1\pm\sqrt{1+r^{-2}q^2})}) \end{array}} }$

となるが， $\pi_2>0$ より

$\displaystyle{(5)\quad \begin{array}{l} \pi_2=r^2\omega_n^{-3}(-2\zeta+\sqrt{4\zeta^2+ 2(-1+\sqrt{1+r^{-2}q^2})})}\\ \pi_3=r^2\omega_n^{-2}(-1+\sqrt{1+r^{-2}q^2}) \end{array} }$

でなければならない。さらに，第２式より $\pi_1$ は

$\displaystyle{(6)\quad \begin{array}{l} \pi_1=(2\zeta\omega_n+r^{-2}\omega_n^4\pi_2)\pi_3+\omega_n^2\pi_2 \\ =(2\zeta\omega_n+r^{-2}\omega_n^4r^2\omega_n^{-3}(-2\zeta+\sqrt{4\zeta^2+ 2(-1+\sqrt{1+r^{-2}q^2})})})\\ \times r^2\omega_n^{-2}(-1+\sqrt{1+r^{-2}q^2})\\ +\omega_n^2r^2\omega_n^{-3}(-2\zeta+\sqrt{4\zeta^2+ 2(-1+\sqrt{1+r^{-2}q^2})})}\\ =r^2\omega_n^{-1}(2\zeta+(-2\zeta+\sqrt{4\zeta^2+ 2(-1+\sqrt{1+r^{-2}q^2})})})(-1+\sqrt{1+r^{-2}q^2})\\ +r^2\omega_n^{-1}(-2\zeta+\sqrt{4\zeta^2+ 2(-1+\sqrt{1+r^{-2}q^2})})}\\ =r^2\omega_n^{-1}\sqrt{4\zeta^2+ 2(-1+\sqrt{1+r^{-2}q^2})}(-1+\sqrt{1+r^{-2}q^2})\\ +r^2\omega_n^{-1}(-2\zeta+\sqrt{4\zeta^2+2(-1+\sqrt{1+r^{-2}q^2})})}\\ =r^2\omega_n^{-1}(-2\zeta+\sqrt{1+r^{-2}q^2}\sqrt{4\zeta^2-2+2\sqrt{1+r^{-2}q^2}}) \end{array} }$

となる（ $\pi_1>0$ ）。すなわち(1)の解として

$\displaystyle{(7)\quad \left\{\begin{array}{l} \pi_1=r^2\omega_n^{-1}(-2\zeta+\sqrt{1+r^{-2}q^2} \sqrt{4\zeta^2-2+2\sqrt{1+r^{-2}q^2}}) \\ \pi_2=r^2\omega_n^{-3}(-2\zeta+\sqrt{4\zeta^2-2+2\sqrt{1+r^{-2}q^2}}) \\ \pi_3=r^2\omega_n^{-2}(-1+\sqrt{1+r^{-2}q^2}) \end{array}\right.}$

を得ます（このとき $\pi_1\pi_2-\pi_3^2>0$ も満足されます）。したがって

$\displaystyle{(8)\quad \begin{array}{l} \left[\begin{array}{cc} f_1 & f_2 \end{array}\right]=r^{-2} \left[\begin{array}{cc} 0 & \omega_n^2 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]=r^{-2}\omega_n^2 \left[\begin{array}{cc} \pi_3 & \pi_2 \end{array}\right]\\ = \left[\begin{array}{cc} -1+\sqrt{1+\left(\frac{q}{r}\right)^2} & -\frac{2}{\omega_n}\left(-\zeta+\sqrt{\zeta^2-\frac{1}{2}+\frac{1}{2}\sqrt{1+\left(\frac{q}{r}\right)^2}}\right) \end{array}\right] \end{array}} }$

Note A52-3d 例題２(22)の導出　

リッカチ方程式

$\displaystyle{(1)\quad \begin{array}{l} \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] + \left[\begin{array}{cc} 0 & -\omega_n^2 \\ 1 & -2\zeta\omega_n \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]\\ - \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{c} 0 \\ \omega_n^2 \end{array}\right] r^{-2} \left[\begin{array}{cc} 0 & \omega_n^2 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]\\ + \left[\begin{array}{cc} q_1^2 & 0 \\ 0 & q_2^2 \end{array}\right] = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \end{array}} }$

を要素ごとに整理して、上と同様にして導出されます。

補遺　上述の議論では、次についての検討が必要です。

検討事項１　(9)の妥当性

リャプノフ方程式

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

リャプノフ方程式…Homework

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

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

これに出力方程式

$\displaystyle{(2)\quad y(t)=Cx(t) }$

を考慮する場合も含めると、漸近安定性の判定法は次のようにまとめられます。

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

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

条件A0:　 $\exp(At)\rightarrow 0\quad(t\rightarrow\infty)$

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

条件A2:　 $\exists \Pi>0: \Pi A+A^T\Pi+I_n=0$

条件A3:　 $\exists \Pi>0: \boxed{\Pi A+A^T\Pi+C^TC=0}$ 　ただし、 $(A,C)$ は可観測対
——————————————————————————————————————–

ここで、条件A2、条件A3における行列方程式はリャプノフ方程式と呼ばれます。
　＜定義DA $\Leftrightarrow$ 条件A0 $\Leftrightarrow$ 条件A1＞
はすでに示しています。以下では、
　＜定義DA $\Rightarrow$ 条件A2 $\Rightarrow$ 条件A1＞および
　＜定義DA $\Rightarrow$ 条件A3 $\Rightarrow$ 条件A1＞
を示して全部の等価性を主張します。

●＜定義DA $\Rightarrow$ 条件A2 $\Rightarrow$ 条件A1＞　定義DAより、(1)の解 $x(t)=\exp(At)x(0)$ の二乗面積は有界となり、

$\displaystyle{(3)\quad \int_0^\infty \underbrace{x^T(t)x(t)}_{V(x(t))}dt=x^T(0)\underbrace{\int_0^\infty\exp(A^Tt)\exp(At)dt}_{\Pi}x(0)>0 }$

から、 $\Pi>0$ が定まります。これは次のようにリャプノフ方程式を満足します。

$(4)\quad \begin{array}{l} \Pi A+A^T\Pi\\ \displaystyle{=\int_0^\infty\exp(A^Tt)\exp(At)dt A+A^T \int_0^\infty\exp(A^Tt)\exp(At)dt}\\ \displaystyle{=\int_0^\infty\frac{d}{dt}(\exp(A^Tt)\exp(At))dt=\left[\exp(A^Tt)\exp(At)\right]_0^\infty\\ =\exp(A^T\infty)\underbrace{\exp(A\infty)}_{0}-\exp(A^T0)\underbrace{\exp(A0)}_{I_n}=-I_n}\\ \end{array} }$

次に、条件A2が成り立つとします。 $A$ が ${\rm Re}(\lambda)\ge0$ を満たす固有値 $\lambda$ をもつとし、これに対応する固有ベクトルを $v\ne0$ とします。このとき

$\displaystyle{(5)\quad \begin{array}{l} 0=v^H(\Pi A+A^T\Pi+I_n)v\\ =v^H\Pi Av+v^HA^T\Pi v+v^Hv\\ =\lambda v^H\Pi v+\lambda^*v^H\Pi v+v^Hv\\ =2\underbrace{{\rm Re}(\lambda)}_{\ge0} \underbrace{v^H\Pi v}_{>0}+\underbrace{v^Hv}_{>0}>0 \end{array} }$

となって矛盾。したがって条件A1を得ます。

●安定行列 $A=\left[\begin{array}{cc} 0 & 1 \\ -2 & -3 \end{array}\right]$ に対して条件A2を確かめます。リャプノフ方程式は

$\displaystyle{(6)\quad \underbrace{ \left[\begin{array}{cc} p_1 & p_3 \\ p_3 & p_2 \end{array}\right] }_{P} \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -2 & -3 \end{array}\right] }_{A} + \underbrace{ \left[\begin{array}{cc} 0 & -2 \\ 1 & -3 \end{array}\right] }_{A^T} \underbrace{ \left[\begin{array}{cc} p_1 & p_3 \\ p_3 & p_2 \end{array}\right] }_{P} =- \underbrace{ \left[\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right] }_{I_2} }$

となり、次の連立1次方程式に書き直すことができます。

$\displaystyle{(7)\quad \left[\begin{array}{ccc} 0 & 0 & -4 \\ 1 & -2 & -3 \\ 0 & -6 & 2 \end{array}\right] \left[\begin{array}{c} p_1 \\ p_2 \\ p_3 \end{array}\right] = \left[\begin{array}{c} -1 \\ 0 \\ -1 \end{array}\right] }$

これより、次を得ます。

$\displaystyle{(8)\quad P= \left[\begin{array}{cc} 1.25 & 0.25 \\ 0.25 & 0.25 \end{array}\right]>0 }$

●＜定義DA $\Rightarrow$ 条件A3 $\Rightarrow$ 条件A1＞　定義DAより、(1),(2)の解 $y(t)=C\exp(At)x(0)$ の二乗面積は有界となり、

$\displaystyle{(9)\quad \int_0^\infty \underbrace{y^T(t)y(t)}_{V(x(t))}dt=x^T(0)\underbrace{\int_0^\infty\exp(A^Tt)C^TC\exp(At)dt}_{\Pi}x(0)>0 }$

から、 $\Pi>0$ が定まります。ここで、 $(A,C)$ は可観測対の条件が、可観測性グラミアン $W_o(t)>0$ を保証し、 $\Pi=W_o(\infty)>0$ が成り立ちます。これは次のようにリャプノフ方程式を満足します。

$(10)\quad \begin{array}{l} \Pi A+A^T\Pi\\ \displaystyle{=\int_0^\infty\exp(A^Tt)C^TC\exp(At)dt A+A^T \int_0^\infty\exp(A^Tt)C^TC\exp(At)dt}\\ \displaystyle{=\int_0^\infty\frac{d}{dt}(\exp(A^Tt)C^TC\exp(At))dt=\left[\exp(A^Tt)C^TC\exp(At)\right]_0^\infty\\ =\exp(A^T\infty)C^TC\underbrace{\exp(A\infty)}_{0}-\exp(A^T0)C^TC\underbrace{\exp(A0)}_{I_n}=-C^TC}\\ \end{array} }$

次に、条件A3が成り立つとします。 $A$ が ${\rm Re}(\lambda)\ge0$ を満たす固有値 $\lambda$ をもつとし、これに対応する固有ベクトルを $v\ne0$ とします。このとき

$\displaystyle{(11)\quad \begin{array}{l} 0=v^H(\Pi A+A^T\Pi+C^TC)v\\ =v^H\Pi Av+v^HA^T\Pi v+v^HC^TCv\\ =\lambda v^H\Pi v+\lambda^*v^H\Pi v+v^HC^TCv\\ =2\underbrace{{\rm Re}(\lambda)}_{\ge0} \underbrace{v^H\Pi v}_{>0}+\underbrace{v^HC^TCv}_{\ge0}\ge0 \end{array} }$

となって矛盾。したがって条件A1を得ます。

Note A51 ラウス・フルビッツの安定判別法

次式で表される $n$ 次自由系を考えます。

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

この漸近安定性は

$\displaystyle{(102)\quad \forall x(0)\ne0: x(t)=\exp(At)x(0)\,\rightarrow\,0\quad(t\rightarrow\infty) }$

で定義され， $A$ 行列のすべての固有値の実部が負であることが必要十分条件であることはよく知られています。 $A$ 行列の中に不確かなパラメータを含む場合，漸近安定性を保証する範囲を求めるために，Routh-Hurwitzの安定判別法が有用です。これは， $A$ 行列の特性多項式

$\displaystyle{(103)\quad {\rm det}(\lambda I_n-A)=\underbrace{a_0}_{1}\lambda^n+a_1\lambda^{n-1}+\cdots+a_{n-1}\lambda+a_n }$

の係数 $a_1,a_2,\cdots,a_n$ から得られるRouth表

$\displaystyle{(104)\quad \boxed{\begin{array}{c|rrrrrr} s^{n} & a_0 & a_2 & a_4 & a_6 & \cdots \\ s^{n-1} & a_1 & a_3 & a_5 & a_7 & \cdots \\ s^{n-2} & \frac{a_1a_2-a_0a_3}{a_1}=b_1 & \frac{a_1a_4-a_0a_5}{a_1}=b_2 & \frac{a_1a_6-a_0a_7}{a_1}=b_3 & \cdots & \cdots \\ s^{n-3} & \frac{b_1a_3-a_1b_2}{b_1}=c_1 & \frac{b_1a_5-a_1b_3}{b_1}=c_2 & \cdots & \cdots & \vdots \\ s^{n-4} & \frac{c_1b_2-b_1c_2}{c_1}=d_1 & \frac{c_1b_3-b_1c_3}{c_1}=d_2 & \cdots & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ s^{0} & & & & & \end{array}} }$

の第１列の要素 $R_1=a_1,R_2=b_1,R_3=c_1,R_4=d_1,\cdots$ がすべて正( $a_0=1$ )であれば漸近安定，または係数 $a_1,a_2,a_3,\cdots$ がすべて正( $a_0=1$ )かつHurwitz行列式

$\displaystyle{(105)\quad \boxed{\begin{array}{l} D_1=a_1 \\ D_2={\rm det}< \left[\begin{array}{cc} a_1 & a_3 \\ a_0 & a_2 \end{array}\right] \\ D_3={\rm det} \left[\begin{array}{ccc} a_1 & a_3 & a_5 \\ a_0 & a_2 & a_4 \\ 0 & a_1 & a_3 \end{array}\right] \\ \vdots \\ D_{n-1}={\rm det} \left[\begin{array}{ccccc} a_1 & a_3 & a_5 & \cdots & a_{2n-3} \\ a_0 & a_2 & a_4 & \cdots & a_{2n-4} \\ 0 & a_1 & a_3 & \cdots & a_{2n-5} \\ 0 & a_0 & a_2 & \cdots & a_{2n-6} \\ 0 & 0 & a_1 & \cdots & a_{2n-7} \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & a_{n-1} \\ \end{array}\right] \end{array}} }$

がすべて正であれば漸近安定と判定します。両者の間には，次が成り立ちます。

$\displaystyle{(106)\quad R_1=\Delta_1,\ R_2=\frac{\Delta_2}{\Delta_1},\cdots,R_n=\frac{\Delta_n}{\Delta_{n-1}} }$

●たとえば、 $n=2,3$ の場合は、次の条件となります。

$\lambda^2+a_1\lambda+a_2$ の場合は、 $a_1>0,a_2>0$

$\lambda^3+a_1\lambda^2+a_2\lambda+a_3$ の場合は、 $\displaystyle{a_1>0,a_2>\frac{a_3}{a_1},a_3>0}$

証明

●Routh-Hurwitzの安定判別法の時間領域における証明については，P.C.Parksによって，自由系が漸近安定であるための必要十分条件

$\displaystyle{(107)\quad \exists P>0: PA+A^TP=-H^TH\le0 }$

に基づくものが示されています（ただし， $(A,H)$ は可観測対）。彼は

$\displaystyle{(108)\quad A=\left[\begin{array}{cccccc} 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1\\ -a_n & -a_{n-1} & -a_{n-2} & \cdots & -a_1 \end{array}\right] }$

を，ある正則行列 $T$ で相似変換した行列

$\displaystyle{(109)\quad B=\left[\begin{array}{ccccccc} 0 & 1 & \cdots & 0 & 0 & 0 \\ -b_n & 0 & \cdots & 0 & 0 & 0 \\ 0 & -b_{n-1} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots& \vdots\\ 0 & 0 & \cdots & -b_3 & 0 & 1 \\ 0 & 0 & \cdots & 0 & -b_2 & -b_1 \end{array}\right]=TAT^{-1} }$

に対して

$\displaystyle{(110)\quad \exists \Pi>0: \Pi B+B^T\Pi=-H^TH\le0 }$

を満足する $\Pi$ と $H$ を，次のように求めています。

$\displaystyle{(111)\quad \Pi=\left[\begin{array}{ccccccc} b_1b_2\cdots b_n & 0 & \cdots & 0 \\ 0 & b_1b_2\cdots b_{n-1} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & b_1 \end{array}\right] }$

$\displaystyle{(112)\quad H=\left[\begin{array}{cccc} 0 & \cdots & 0 & \sqrt{2}b_1 \end{array}\right] }$

実際， $\Pi B$ を計算してみると

$\displaystyle{(113)\quad \left[\begin{array}{ccccccc} b_1b_2\cdots b_n & 0 & \cdots & 0 \\ 0 & b_1b_2\cdots b_{n-1} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & b_1 \end{array}\right] \left[\begin{array}{ccccccc} 0 & 1 & \cdots & 0 & 0 & 0 \\ -b_n & 0 & \cdots & 0 & 0 & 0 \\ 0 & -b_{n-1} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots& \vdots\\ 0 & 0 & \cdots & -b_3 & 0 & 1 \\ 0 & 0 & \cdots & 0 & -b_2 & -b_1 \end{array}\right] }$
$\displaystyle{ = \left[\begin{array}{ccccccc} 0 & b_1b_2\cdots b_n & \cdots & 0 & 0 & 0 \\ -b_1b_2\cdots b_n & 0 & \cdots & 0 & 0 & 0 \\ 0 & -b_1b_2\cdots b_{n-1} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots& \vdots\\ 0 & 0 & \cdots & -b_1b_2b_3 & 0 & b_1b_2 \\ 0 & 0 & \cdots & 0 & -b_1b_2 & -b_1^2 \end{array}\right] }$

となります。したがって

$\displaystyle{(114)\quad &&\Pi B+B^T\Pi= \left[\begin{array}{ccccccc} 0 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & -2b_1^2 \end{array}\right] =- \underbrace{ \left[\begin{array}{cccc} 0 \\ \vdots \\ 0 \\ \sqrt{2}b_1 \end{array}\right] }_{H^T} \underbrace{ \left[\begin{array}{cccc} 0 & \cdots & 0 & \sqrt{2}b_1 \end{array}\right] }_{H} \nonumber }$

が得られます。明らかに， $b_n\ne0$ のとき

$\displaystyle{(115)\quad {\rm rank} \left[\begin{array}{ccccccc} B-\lambda I_n \\ H \end{array}\right] = {\rm rank} \left[\begin{array}{ccccccc} -\lambda & 1 & \cdots & 0 & 0 & 0 \\ -b_n & -\lambda & \cdots & 0 & 0 & 0 \\ 0 & -b_{n-1} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots& \vdots\\ 0 & 0 & \cdots & -b_3 & -\lambda& 1 \\ 0 & 0 & \cdots & 0 & -b_2 & -b_1-\lambda \\ 0 & 0 & \cdots & 0 & 0 & \sqrt{2}b_1 \end{array}\right]=n }$

となり， $(B,H)$ は可観測対です。また，関係式

$\displaystyle{(116)\quad &&b_1=\Delta_1,\ b_2=\frac{\Delta_2}{\Delta_1},\ b_3=\frac{\Delta_3}{\Delta_1\Delta_2},\ b_4=\frac{\Delta_1\Delta_4}{\Delta_2\Delta_3},\ \cdots,\ b_n=\frac{\Delta_{n-3}\Delta_n}{\Delta_{n-2}\Delta_{n-1}} }$

が成り立つことが示されており， $\Pi$ の正定性からRouth-Hurwitzの安定判別法の妥当性が出ます。

●以下では、２つの行列

$\displaystyle{(117)\quad A=\left[\begin{array}{cccccc} 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1\\ -a_n & -a_{n-1} & -a_{n-2} & \cdots & -a_1 \end{array}\right] }$

$\displaystyle{(118)\quad B=\left[\begin{array}{ccccccc} 0 & 1 & \cdots & 0 & 0 & 0 \\ -b_n & 0 & \cdots & 0 & 0 & 0 \\ 0 & -b_{n-1} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots& \vdots\\ 0 & 0 & \cdots & -b_3 & 0 & 1 \\ 0 & 0 & \cdots & 0 & -b_2 & -b_1 \end{array}\right] }$

に対して

$\displaystyle{(119)\quad B=TAT^{-1} }$

を満足する座標変換行列 $T$ を導出します。ここで

$\displaystyle{(120)\quad {\rm det}(\lambda I_n-B)={\rm det}(\lambda I_n-A) =\lambda^n+a_1\lambda^{n-1}+\cdots+a_{n-1\lambda+}a_n }$

が成り立ちます。

さて、次の関係式はよく知られています。

$\displaystyle{(121)\quad A=V_A\Lambda V_A^{-1} }$

ただし

$\displaystyle{(122)\quad \Lambda=\left[\begin{array}{ccccccc} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \lambda_n \end{array}\right] }$

$\displaystyle{(123)\quad V_A=\left[\begin{array}{ccccccc} 1 & 1 & \cdots & 1 \\ \lambda_1 & \lambda_2 & \cdots & \lambda_n \\ \vdots & \vdots & \ddots & \vdots \\ \lambda_1^{n-1} & \lambda_2^{n-1} & \cdots & \lambda_n^{n-1} \end{array}\right] }$

ここで，簡単のために， $A$ の固有値 $\{\lambda_1,\lambda_2,\cdots,\lambda_n\}$ は互いに相異なるものとします。このとき

$\displaystyle{(124)\quad B=V_B\Lambda V_B^{-1} }$

を満足する $V_B$ を求めれば，座標変換行列 $T$ は

$\displaystyle{(125)\quad T=V_BV_A^{-1} }$

のように表されます。

(119)に関する補足

● $n=2$ のとき

$\displaystyle{ A=\left[\begin{array}{cccccc} 0 & 1 \\ -a_2 & -a_1 \end{array}\right] }$
$\displaystyle{ B=\left[\begin{array}{cccccc} 0 & 1 \\ -b_2 & -b_1 \end{array}\right] }$

ですから，自明な場合

$\displaystyle{ b_1=a_1=\Delta_1,\ b_2=a_2=\frac{\Delta_2}{\Delta_1} }$

で， $V_B=V_A$ ，すなわち $T=I_2$ となります。

● $n=3$ のとき

$\displaystyle{ A=\left[\begin{array}{cccccc} 0 & 1 & 0 \\ 0 & 0 & 1 \\ -a_3 & -a_2 & -a_1 \end{array}\right] }$
$\displaystyle{ B=\left[\begin{array}{cccccc} 0 & 1 & 0 \\ -b_3 & 0 & 1 \\ 0 & -b_2 & -b_1 \end{array}\right] }$

となります。いま

$\displaystyle{ \left[\begin{array}{cccccc} 0 & 1 & 0 \\ -b_3 & 0 & 1 \\ 0 & -b_2 & -b_1 \end{array}\right] \left[\begin{array}{c} v_1 \\ v_2 \\ v_3 \end{array}\right]= \lambda \left[\begin{array}{c} v_1 \\ v_2 \\ v_3 \end{array}\right] }$

すなわち

$\displaystyle{ \left\{\begin{array}{l} v_2=\lambda v_1 \\ v_3=\lambda v_2 +b_3v_1\\ \lambda v_3+b_1v_3+b_2v_2=0 \end{array}\right. }$

となる $v_1=1,v_2,v_3$ を求めると

$\displaystyle{ \left[\begin{array}{c} v_1 \\ v_2 \\ v_3 \end{array}\right]= \left[\begin{array}{l} 1 \\ \lambda \\ \lambda^2+b_3 \end{array}\right] = \underbrace{ \left[\begin{array}{cccccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ b_3 & 0 & 1 \end{array}\right]}_{\Gamma} \left[\begin{array}{c} 1 \\ \lambda \\ \lambda^2 \end{array}\right] }$

ただし

$\displaystyle{ \lambda^3+b_1\lambda^2+(b_2+b_3)\lambda+b_1b_3=0 }$

これは

$\displaystyle{ \lambda^3+a_1\lambda^2+a_2\lambda+a_3=0 }$

と一致しなければならないので

$\displaystyle{ b_1=a_1,\ b_2=a_2-\frac{a_3}{a_1},\ b_3=\frac{a_3}{a_1} }$

となります。また，Hurwitz行列式を用いて

$\displaystyle{ b_1=\Delta_1,\ b_2=\frac{\Delta_2}{\Delta_1},\ b_3=\frac{\Delta_3}{\Delta_1\Delta_2} }$

と表されるます。実際

以上から， $V_B=\Gamma V_A$ ，すなわち $T=\Gamma$ となります。

● $n=4$ のとき

$\displaystyle{ A=\left[\begin{array}{cccccc} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -a_4 & -a_3 & -a_2 & -a_1 \end{array}\right] }$

$\displaystyle{ B=\left[\begin{array}{cccccc} 0 & 1 & 0 & 0 \\ -b_4 & 0 & 1 & 0 \\ 0 & -b_3 & 0 & 1 \\ 0 & 0 & -b_2 & -b_1 \end{array}\right] }$

となります。いま

$\displaystyle{ \left[\begin{array}{cccccc} 0 & 1 & 0 & 0 \\ -b_4 & 0 & 1 & 0 \\ 0 & -b_3 & 0 & 1 \\ 0 & 0 & -b_2 & -b_1 \end{array}\right] \left[\begin{array}{c} v_1 \\ v_2 \\ v_3 \\ v_4 \end{array}\right]= \lambda \left[\begin{array}{c} v_1 \\ v_2 \\ v_3 \\ v_4 \end{array}\right] }$

すなわち

$\displaystyle{ \left\{\begin{array}{l} v_2=\lambda v_1 \\ v_3=\lambda v_2 +b_4v_1\\ v_4=\lambda v_3 +b_3v_2\\ \lambda v_4+b_1v_4+b_2v_3=0 \end{array}\right. }$

となる $v_1=1,v_2,v_3,v_4$ を求めると

$\displaystyle{ \left[\begin{array}{c} v_1 \\ v_2 \\ v_3 \\ v_4 \end{array}\right]= \left[\begin{array}{l} 1 \\ \lambda \\ \lambda^2+b_4 \\ \lambda^3+(b_3+b_4)\lambda \end{array}\right] = \underbrace{ \left[\begin{array}{cccccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ b_4 & 0 & 1 & 0 \\ 0 & b_3+b_4 & 0 & 1 \end{array}\right]}_{\Gamma} \left[\begin{array}{c} 1 \\ \lambda \\ \lambda^2 \\ \lambda^3 \end{array}\right] }$

ただし

$\displaystyle{ \lambda^4+b_1\lambda^3+(b_2+b_3+b_4)\lambda^2+b_1(b_3+b_4)\lambda+b_2b_4=0 }$

これは

$\displaystyle{ \lambda^4+a_1\lambda^3+a_2\lambda^2+a_3\lambda+a_4=0 }$

と一致しなければならないので

$\displaystyle{ b_1=a_1,\ b_2=a_2-\frac{a_3}{a_1},\ b_3=\frac{a_3}{a_1}-\frac{a_4}{a_2-\frac{a_3}{a_1}},\ b_4=\frac{a_4}{a_2-\frac{a_3}{a_1}} }$

となります。また，Hurwitz行列式を用いて

$\displaystyle{ b_1=\Delta_1,\ b_2=\frac{\Delta_2}{\Delta_1},\ b_3=\frac{\Delta_3}{\Delta_1\Delta_2},\ b_4=\frac{\Delta_1\Delta_4}{\Delta_2\Delta_3} }$

と表されます。実際

以上から， $V_B=\Gamma V_A$ ，すなわち $T=\Gamma$ となります。

● $n=5$ のとき
$\displaystyle{ A=\left[\begin{array}{cccccc} 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ -a_5 & -a_4 & -a_3 & -a_2 & -a_1 \end{array}\right] }$

$\displaystyle{ B=\left[\begin{array}{cccccc} 0 & 1 & 0 & 0 & 0 \\ -b_5 & 0 & 1 & 0 & 0 \\ 0 & -b_4 & 0 & 1 & 0 \\ 0 & 0 & -b_3 & 0 & 1 \\ 0 & 0 & 0 & -b_2 & -b_1 \end{array}\right] }$

となります。いま

$\displaystyle{ \left[\begin{array}{cccccc} 0 & 1 & 0 & 0 & 0 \\ -b_5 & 0 & 1 & 0 & 0 \\ 0 & -b_4 & 0 & 1 & 0 \\ 0 & 0 & -b_3 & 0 & 1 \\ 0 & 0 & 0 & -b_2 & -b_1 \end{array}\right] \left[\begin{array}{c} v_1 \\ v_2 \\ v_3 \\ v_4 \\ v_5 \end{array}\right]= \lambda \left[\begin{array}{c} v_1 \\ v_2 \\ v_3 \\ v_4 \\ v_5 \end{array}\right] }$

すなわち

$\displaystyle{ \left\{\begin{array}{l} v_2=\lambda v_1 \\ v_3=\lambda v_2 +b_5v_1\\ v_4=\lambda v_3 +b_4v_2\\ v_5=\lambda v_4 +b_3v_3\\ \lambda v_5+b_1v_5+b_2v_4=0 \end{array}\right. }$

となる $v_1=1,v_2,v_3,v_4,v_5$ を求めると

$\displaystyle{ \left[\begin{array}{c} v_1 \\ v_2 \\ v_3 \\ v_4 \\ v_5 \end{array}\right]= \left[\begin{array}{l} 1 \\ \lambda \\ \lambda^2+b_5 \\ \lambda^3+(b_4+b_5)\lambda \\ \lambda^4+(b_3+b_4+b_5)\lambda^2+b_3b_5 \end{array}\right] }$
$\displaystyle{ = \underbrace{ \left[\begin{array}{cccccc} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ b_5 & 0 & 1 & 0 & 0 \\ 0 & b_4+b_5 & 0 & 1 & 0 \\ b_3b_5 & 0 & b_3+b_4+b_5 & 0 & 1 \end{array}\right]}_{\Gamma} \left[\begin{array}{c} 1 \\ \lambda \\ \lambda^2 \\ \lambda^3 \\ \lambda^4 \end{array}\right] \nonumber }$

ただし

$\displaystyle{ \lambda^5+b_1\lambda^4+(b_2+b_3+b_4+b_5)\lambda^3+b_1(b_3+b_4+b_5)\lambda^2+(b_2(b_4+b_5)+b_3b_5)\lambda+b_1b_3b_5=0\nonumber }$

これは

$\displaystyle{ \lambda^5+a_1\lambda^4+a_2\lambda^3+a_3\lambda^2+a_4\lambda+a_5=0 }$

と一致しなければならないので

$\displaystyle{ &&b_1=a_1,\ b_2=a_2-\frac{a_3}{a_1},\ b_3=\frac{a_3}{a_1}-\frac{a_4-\frac{a_5}{a_1}}{b_2},\ b_4=\frac{a_4-\frac{a_5}{a_1}}{b_2}-\frac{a_5}{a_1b_3},\ b_5=\frac{a_5}{a_1b_3} \nonumber }$

となります。また，Hurwitz行列式を用いて

$\displaystyle{ b_1=\Delta_1,\ b_2=\frac{\Delta_2}{\Delta_1},\ b_3=\frac{\Delta_3}{\Delta_1\Delta_2},\ b_4=\frac{\Delta_1\Delta_4}{\Delta_2\Delta_3},\ b_5=\frac{\Delta_2\Delta_5}{\Delta_3\Delta_4} }$

と表される。実際

$\displaystyle{ \begin{array}{l} b_1=\Delta_1=a_1 \\ b_2=\frac{\Delta_2}{\Delta_1}=\frac{a_1a_2-a_3}{a_1}=a_2-\frac{a_3}{a_1} \\ b_3=\frac{\Delta_3}{\Delta_1\Delta_2}=\frac{\Delta_2a_3-(a_1a_4-a_5)a_1}{a_1\Delta_2}=\frac{a_3}{a_1}-\frac{a_4-\frac{a_5}{a_1}}{b_2} \\ b_4=\frac{\Delta_1\Delta_4}{\Delta_2\Delta_3}=\frac{1}{b_2}(a_4-\frac{a_5}{a_1}-\frac{a_5(a_2-\frac{a_3}{a_1})}{R_3})=\frac{a_4-\frac{a_5}{a_1}}{b_2}-\frac{a_5}{a_1b_3} \\ b_5=\frac{\Delta_2\Delta_5}{\Delta_3\Delta_4}=\frac{\Delta_2\Delta_4a_5}{\Delta_3\Delta_4}=\frac{a_5}{a_1b_3} \end{array} }$

以上から， $V_B=\Gamma V_A$ ，すなわち $T=\Gamma$ となる。

●ここまで調べれば、 $n=i$ の場合の(119)から $n=i+1$ の場合の(119)を示すことができる思います（あとでやってみます）。

可検出性と可観測性

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

可検出性と可観測性…Homework

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

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

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

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

で与えられました。ここで、オブザーバゲイン $H$ を $A-HC$ が安定行列になるように選ぶ必要がありました。でもこれはいつでも可能であるわけではありません。すなわち $H$ を $A-HC$ が安定行列になるように選べる場合に状態オブザーバが構成可能となり、このとき可検出性(detectability)が成り立つと言います。

$A-HC$ を安定行列とするためには、双対システム

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

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

$\displaystyle{(4)\quad v(t)=-H^Tw(t)}$

による閉ループ系

$\displaystyle{(5)\quad \dot{w}(t)=(A^T-C^TH^T)w(t)=(A-HC)^Tw(t)}$

を安定化して求めます。ということは、可検出性の条件は双対システムの可安定性と等価になります。すなわち可安定性の条件において

$\displaystyle{(6)\quad A\leftarrow A^T,\ B\leftarrow C^T}$

のように置き換えて、次が得られます。

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

定義DD:　状態オブザーバを構成可能

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

条件D2:　 $Cv=0,\ Av=\lambda v \Rightarrow v=0$ 　( $\lambda$ は $A$ のすべての不安定固有値)

可検出性の判定は行列 $A$ と $C$ を用いて行われるので、可検出性が成り立つとき、対 $(A,C)$ は可検出対という言い方をします。

可安定性よりも強い条件として可制御性がありましたが、これに対応する概念は可観測性(observability)と言われ、これと等価な条件は次のようにまとめられます。

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

定義DO:　任意有限時間 $[0,t]$ 上の入出力データ $u(\tau),y(\tau)\tau\in[0,t]$ から，初期状態 $x(0)$ を一意に決定できる

条件O1:　 $\displaystyle{\underbrace{\int_0^t \exp(A^T\tau)C^TC\exp(A\tau)\,d\tau}_{W_o(t)}>0 \quad (\forall t>0)}$

条件O2:　 ${\rm rank}\, \underbrace{ \left[\begin{array}{c} C \\ CA \\ \vdots\\ CA^{n-1} \end{array}\right] }_{observability\ matrix} =n$

条件O3:　 $H$ を選んで， $A-HC$ の固有値を任意に設定可能

条件O4:　 ${\rm rank} \left[\begin{array}{c} C \\ A-\lambda I_n \end{array}\right]=n$ 　( $\lambda$ は $A$ のすべての固有値)

条件O5:　 $Cv=0,\ Av=\lambda v \Rightarrow v=0$ 　( $\lambda$ は $A$ のすべての固有値)

条件O6:　 $(A^T,C^T)$ は可制御対

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

上で述べたように、可観測性の議論は、可制御性の議論において、 $A$ を $A^T$ に、 $B$ を $C^T$ に置き換えて行えばよいので、次が成り立ちます。

条件O6 $\Leftrightarrow$ 条件O1 $\Leftrightarrow$ 条件O2 $\Leftrightarrow$ 条件O3 $\Leftrightarrow$ 条件O4 $\Leftrightarrow$ 条件O5

●＜定義DO $\Rightarrow$ 条件O1＞

適当な入力のもとでの出力は次式で表されます。

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

ここで、零状態応答は計算できますので、これを出力から引いた $z(t)$ を次のようにおきます。

$\displaystyle{(9)\quad \underbrace{y(t)-\int_0^tC\exp(A\tau)Bu(t-\tau)\,d\tau}_{z(t)}=C\exp(At)x(0)}$

したがって、入出力データ $u(\tau),z(\tau)\tau\in[0,t]$ から，初期状態 $x(0)$ を一意に決定できるかが問題となります。

このとき、 $W_o(t)>0$ を否定して、矛盾を導きます。

$(10)\quad \begin{array}{l} \displaystyle{\overline{(\forall v\ne0:v^TW_o(t)v>0)}}\\ \displaystyle{\Rightarrow \exists v\ne0:w^TW_o(t)v=0}\\ \displaystyle{\Rightarrow \exists v\ne0,\forall \tau\le t: C\exp(A\tau)v=0}\\ \displaystyle{\Rightarrow x(0)=v\Rightarrow z(t)=0} \end{array}$

これは、入出力データ $u(\tau),z(\tau)\tau\in[0,t]$ から，初期状態 $x(0)$ を一意に決定できることにはならないことを示しています。

●＜条件O1 $\Rightarrow$ 定義DO＞

上と同様に(9)を考えますと、入出力データ $u(\tau),z(\tau)\tau\in[0,t]$ から，初期状態 $x(0)$ は次式で計算できます。

$\displaystyle{(11)\quad x(0)=W_o^{-1}(t)\int_0^t\exp(A^T\tau)C^Tz(\tau)\,d\tau}$

実際、左辺は、次のように $x(0)$ となります。

$\displaystyle{(12)\quad (\int_0^t\exp(A^T\tau)C^TC\exp(A\tau)\,d\tau)^{-1}\int_0^t\exp(A^T\tau)C^TC\exp(At)x(0)\,d\tau=x(0) }$

以上で、可観測性 $\Leftrightarrow$ O1 $\Leftrightarrow$ O2が示されました。

Note A43 可安定性と可制御性

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

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

定義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$ は可制御固有値、満足しない固有値を不可制御固有値と呼びます。

最小次元状態オブザーバ

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

を得ます。ここで

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

を仮定し、 $\hat{J}$ を

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

とおけば

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

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

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

を仮定すると

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

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

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

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

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

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

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

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

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

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

● モータの状態方程式

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

と出力方程式

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

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

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

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

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

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

図1 状態推定の例

演習 A42…Flipped Classroom

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

MATLAB

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

SCILAB

coming soon

Note A42　線形関数オブザーバ

制御対象の状態空間表現

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

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

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

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

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

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

から

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

を得ます。また

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

状態オブザーバ

投稿日時: 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$ のすべての固有値)

[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)$ も可制御対であることは次のように証明されます。