制御系CAD

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

CT72 2次系のLQ制御

投稿日時: 2023年12月30日投稿者: cacsd

Home Work 7.2

[P]　よく多変数制御というけど、制御対象の運動の複数の自由度を同時に操ることを意味しているのであれば、本当にできるのかなと思ってしまう。

[C]　制御対象の可制御性と可観測性を前提とし、２次形式評価関数を最小化する状態FBを求めるLQ制御は、多変数制御のための強力なツールとなっているらしいよ。

[M]　その導出は、基本的には１入力系の場合と同じアプローチでできて、最終的にはリッカチ方程式という行列の２次方程式を解くことになるんだ。　

Flipped Classroom 7.2
[1] 仮に $q_2=0$ としても、閉ループ２次系の減衰係数は $1/sqrt(2)$ となっており、十分減衰がかかっていることを意味しています。
[2] 著者は、プログラムopt.mを用いて、多くの制御問題を解いてきました。世界一重要なツールと言っても過言ではないと思っています。

CT71 １次系のLQ制御

投稿日時: 2023年12月30日投稿者: cacsd

Home Work 7.1

[P]　平衡状態から外れたとき、これを戻すためには、いろいろな「軌道（時系列）」があると思うけど、何かリーズナブルな決め方はあるのかな？

[C]　状態変数の振る舞い（２乗面積）と操作変数の振る舞い（２乗面積）のトレードオフを計ることがよく行われているよ。

[M]　トレードオフの意味だけど、高校で、次の関数の最小値を求める問題があったね。

$\displaystyle{(1)\quad y=x+\frac{1}{x} }$

右辺第１項は増加関数、第２項は減少関数だから、グラフを描いてみると、最小値が生まれているよね。

Flipped Classroom 7.1
[1]

評価関数はスカラ値を取るので、定数倍してもOKだから、 $q^2$ で割ると、

$\displaystyle{(1)\quad \frac{J}{q^2}=\int_0^\infty (x^2(t)+(\frac{r}{q})^2u^2(t))dt =\int_0^\infty (x^2(t)+\frac{1}{\rho^2}u^2(t))dt }$

これから $1/\rho$ は操作のためのコストとみなせます。すなわちコスト $1/\rho$ がより小さくなると、操作はより大きな振る舞いが許されます。 $\rho\rightarrow\infty$ のときはcheap controlと呼ばれます。
　
[2] 　

$\displaystyle{(1)\quad \begin{array}{l} (\lambda I_2- \left[\begin{array}{cc} a & -r^{-2}b^2\\ -q^2 & -a \end{array}\right]) =(\lambda-a)(\lambda+a)-q^2r^{-2}b^2=\lambda^2-a^2-q^2r^{-2}b^2=0 \end{array} }$

を解くと

$\displaystyle{(2)\quad \lambda=\pm\sqrt{a^2+\frac{q^2}{r^2}b^2} }$

$\lambda<0$ の場合の固有ベクトルは

$\displaystyle{(3)\quad %\begin{array}{l} \left[\begin{array}{cc} a & -r^{-2}b^2\\ -q^2 & -a \end{array}\right] \left[\begin{array}{cc} v_1\\ v_2 \end{array}\right] =-\sqrt{a^2+\frac{q^2}{r^2}b^2} \left[\begin{array}{cc} v_1\\ v_2 \end{array}\right] %\end{array} }$

すなわち

$\displaystyle{(4)\quad \left\{\begin{array}{l} av_1-r^{-2}b^2v_2=-\sqrt{a^2+\frac{q^2}{r^2}b^2}v_1\\ -q^2v_1-av_2=-\sqrt{a^2+\frac{q^2}{r^2}b^2}v_2 \end{array}\right. }$

$v_1=1$ の場合(と規格化する場合)は

$\displaystyle{(4)\quad %\left\{\begin{array}{l} a-r^{-2}b^2v_2=-\sqrt{a^2+\frac{q^2}{r^2}b^2} \Rightarrow v_2=\frac{a+\sqrt{a^2+\frac{q^2}{r^2}b^2}}{r^{-2}b^2} %\end{array}\right. }$

したがって

$\displaystyle{(5)\quad %\left\{\begin{array}{l} f=r^{-2}b^2v_2v_1^{-1}=r^{-2}b^2\frac{a+\sqrt{a^2+\frac{q^2}{r^2}b^2}}{r^{-2}b^2} =a+\sqrt{a^2+\frac{q^2}{r^2}b^2} %\end{array}\right. }$

CT64 補遺6

投稿日時: 2023年12月30日投稿者: cacsd

n次系の周波数応答

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

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

このとき、正弦波入力

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

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

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

を、零状態応答の式

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

に代入して

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

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

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

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

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

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

CT63 ２次系のステップ応答

投稿日時: 2023年12月30日投稿者: cacsd

Home Work 6.3

[P]　物理の言葉ではないと思うけど、無定位系とよばれる制御対象があるよ。モータやビークルのように、復原項がない場合は、位置が定まらないという意味らしい。無定位系では、位置についてステップ応答の実験ができないよね。

[C]　単位フィードバックを行えば、ステップ応答の実験ができるのだけど、フィードバックの知識がないと気が付かないかな。

[M]　位置と速度を状態変数にもつ無定位系

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

に対して、単位フィードバック

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

によって復原項を入れると、次のような２次系になるからだね。

$\displaystyle{(3)\quad \begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t)\\ \dot{x}_2(t) \end{array}\right] = \left[\begin{array}{cc} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] + \left[\begin{array}{cc} 0 \\ \omega_n^2 \end{array}\right] u(t)\\ y(t)= \left[\begin{array}{cc} 1 & 0 \end{array}\right] \left[\begin{array}{c} x_1(t)\\ x_2(t) \end{array}\right] \end{array} }$

Flipped Classroom 6.3
[1] 　

[2]

$1^\circ$ $\zeta>1$ のとき、インパルス応答を積分して、ステップ応答は次のように計算されます。

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

$2^\circ$ $\zeta=1$ のとき、インパルス応答を積分して、ステップ応答は次のように計算されます。

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

$3^\circ$ $\zeta<1$ のとき、インパルス応答を積分して、公式

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

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

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

[3] 　

ヒントより、 $T_p=\frac{\pi}{\lambda_I}$ だから

$\displaystyle{(8)\quad \begin{array}{l} S(T_p)=1-\frac{\omega_n}{\lambda_I}e^{\lambda_RT_p}\sin(\lambda_IT_p+\phi)\\ =1-\frac{\omega_n}{\lambda_I}e^{\lambda_R\frac{\pi}{\lambda_I}}\sin(\lambda_I\frac{\pi}{\lambda_I}+\phi)\\ =1-\frac{\omega_n}{\lambda_I}e^{\frac{\lambda_R}{\lambda_I}\pi}\sin(\pi+\phi) =1+\frac{\omega_n}{\lambda_I}e^{\frac{\lambda_R}{\lambda_I}\pi}\sin(\phi)\\ =1+\frac{\omega_n}{\lambda_I}e^{\frac{\lambda_R}{\lambda_I}\pi}\frac{\lambda_I}{\omega_n} =1+e^{\frac{-\zeta\omega_n}{\omega_n\sqrt{1-\zeta^2}}\pi}=1+e^{\frac{-\zeta\pi}{\sqrt{1-\zeta^2}}} \end{array} }$

[4]　

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

[5] 　

CT62 ２次系のインパルス応答

投稿日時: 2023年12月30日投稿者: cacsd

Home Work 6.2

[P]　一般にインパルス応答を実験で求めようとすると大変そうだね。構造試験ではインパルス・ハンマーという道具が使用されるらしいけど、やはりインパルス応答が調べられているんだね。

[M]　インパルス応答の定義式(6.24)をよく見てみると、

$\displaystyle{(6.24)\quad G(t)=C\exp(At)B=C\exp(At)\left[b_1 \cdots b_m\right] }$

となっていて、これは $x(0)=b_1,\cdots,b_m$ の場合の零入力応答を並べたものと言えるね。そうすると、数学的には、インパルス応答は初期値を $B$ の各列ベクトルに設定したときの応答ということになるね。

[C]　この考え方は、非線形シミュレータで適当な初期値を与えたいときに役に立ちそうだね。

Flipped Classroom 6.2
[1] 　

$\displaystyle{(1)\quad \begin{array}{l} (\lambda I_2- \left[\begin{array}{cc} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right]) =\lambda(\lambda+2\zeta\omega_n)+\omega_n^2=\lambda^2+2\zeta\omega_n\lambda+\omega_n^2=0 \end{array} }$

を解くと

$\displaystyle{(2)\quad \lambda=\frac{1}{2}(-2\zeta\omega_n\pm \sqrt{(2\zeta\omega_n)^2-4\omega_n^2}) =-\zeta\omega_n\pm \omega_n\sqrt{\zeta^2-1} }$

これは次のように場合分けされます。

$1^\circ$ $\zeta>1$ のとき、 $\lambda_1=-\zeta\omega_n+ \omega_n\sqrt{\zeta^2-1}, \lambda_2=-\zeta\omega_n- \omega_n\sqrt{\zeta^2-1}$

$2^\circ$ $\zeta=1$ のとき、 $\lambda_1=\lambda_2=-\omega_n$

$3^\circ$ $\zeta<1$ のとき、 $\lambda_1=-\zeta\omega_n+ j\omega_n\sqrt{1-\zeta^2}, \lambda_2=-\zeta\omega_n- j\omega_n\sqrt{1-\zeta^2}$

[2]

$1^\circ$ $\zeta>1$ のとき

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

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

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

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

$2^\circ$ $\zeta=1$ のとき、 $\lambda=-\omega_n$ とおいて

$\displaystyle{(6.28.2)\quad \underbrace{ \left[\begin{array}{cc} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] }_{A} = \underbrace{ \left[\begin{array}{cc} 1 & 1\\ \lambda & \lambda+1 \end{array}\right] }_{V} \underbrace{ \left[\begin{array}{cc} \lambda & 1\\ 0 & \lambda \end{array}\right] }_{\Lambda} \underbrace{ \left[\begin{array}{cc} 1 & 1\\ \lambda & \lambda+1 \end{array}\right]^{-1} }_{V^{-1}}}$

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

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

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

$3^\circ$ $\zeta<1$ のとき、 $\lambda_R=-\zeta\omega_n, \lambda_I=\omega_n\sqrt{1-\zeta^2}$ とおいて

$\displaystyle{(6.28.3)\quad \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] }_{A} = \underbrace{ \left[\begin{array}{cc} 1 & 0 \\ \lambda_R & \lambda_I \end{array}\right] }_{V} \underbrace{ \left[\begin{array}{cc} \lambda_R & \lambda_I \\ -\lambda_I & \lambda_R \end{array}\right] }_{\Lambda} \underbrace{ \frac{1}{\lambda_I} \left[\begin{array}{cc} \lambda_I & 0 \\ -\lambda_R & 1 \end{array}\right] }_{V^{-1}}}$

$\displaystyle{(3.29)\quad \exp(\Lambda t)= e^{\lambda_R t} \left[\begin{array}{cc} \cos\lambda_I t & \sin\lambda_I t\\ \sin\lambda_I t & \cos\lambda_I t \end{array}\right]}$

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

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

CT61 １次系の応答

投稿日時: 2023年12月30日投稿者: cacsd

Home Work 6.1

[P]　１次系の実用例としては、心臓のペースメーカがあるね。これはRC電気回路で、コンデンサへの充放電を繰り返しているよ。

[M]　コンデンサの放電はいわゆる零入力応答（初期値応答）だけど、一般の外部入力があるときは非斉次微分方程式となって、数学的には面倒だね。でもテキストの説明は比較的わかりやすいかな。

[C]　１次系と言えば時定数、その説明には２通り（初期値応答またはステップ応答を用いて）あるかな。船舶分野では、１次系を「野本モデル」とよび、世界的に通用するらしいよ。その時定数は船が真っすぐ進むかどうか（進路安定性）の判断に使われ、なんと時定数が負の場合もあるというから驚きだね。

Flipped Classroom 6.1
[1] 　(6.9)を用いた場合

$\displaystyle{(1)\quad \begin{array}{l} y(t)=ce^{at}x(0)+\int_0^t ce^{a\tau}b\cdot u(t-\tau)d\tau\\ =1\cdot e^{-\frac{1}{T}t}\cdot 0+\int_0^t 1\cdot e^{-\frac{1}{T}(\tau)}\frac{K}{T}\cdot 1d\tau\\ =\frac{K}{T}\int_0^t e^{-\frac{1}{T}\tau}d\tau =\frac{K}{T}\frac{1}{-\frac{1}{T}}\left[e^{-\frac{1}{T}\tau}\right]_0^t\\ =-K(e^{-\frac{1}{T}t}-1) =K(1-e^{-\frac{1}{T}t}) \end{array} }$

ちなみに、(6.8)を用いた場合

$\displaystyle{(2)\quad \begin{array}{l} y(t)=ce^{at}x(0)+\int_0^t ce^{a(t-\tau)}b\cdot u(\tau)d\tau\\ =1\cdot e^{-\frac{1}{T}t}\cdot 0+\int_0^t 1\cdot e^{-\frac{1}{T}(t-\tau)}\frac{K}{T}\cdot 1d\tau\\ =\frac{K}{T}e^{-\frac{1}{T}t} \int_0^t e^{\frac{1}{T}\tau}d\tau =\frac{K}{T}e^{-\frac{1}{T}t} \frac{1}{\frac{1}{T}}\left[e^{\frac{1}{T}\tau}\right]_0^t\\ =Ke^{-\frac{1}{T}t}(e^{\frac{1}{T}t}-1) =K(1-e^{-\frac{1}{T}t}) \end{array} }$

[2] 　 $t=0$ における接線の式は、傾きが $\dot{y}(0)$ 、切片が $y(0)$ だから

$\displaystyle{(3)\quad \dot{y}(0)t+y(0)=K\frac{1}{T}t+0=\frac{K}{T}t }$

これが $K$ の値をとる時刻は $t=T$ となります。また、 $y(\infty)=K$ だから

$\displaystyle{(4)\quad \frac{K(1-e^{-\frac{1}{T}T})}{K}=1-\frac{1}{e}=1-1/2.718=0.632 }$

[3]　図6.2の〇で示す点をクリックすればよいと言えます。

CT54 補遺5

投稿日時: 2023年12月30日投稿者: cacsd

まず穴埋め問題の解答を示します。太字がキーワードなので、確認してください。

[1] 制御対象の状態空間表現が与えられているとする。これは「状態方程式」 ${\dot{x}(t)=Ax(t)+Bu(t)}$ と観測方程式 ${y(t)=Cx(t)}$ からなる。ここで、 ${x(t),u(t),y(t)}$ はそれぞれ次数 ${n,m,p}$ をもつベクトル、 ${A,B,C}$ はそれぞれサイズ ${n\times n,n\times m,p\times n}$ をもつ行列である。また、初期状態「 ${x(0)=0}$ 」は、平衡状態を表すように選ばれているものとする。

[2] いま、何らかの外乱により平衡状態が乱され、 ${x(0)\ne0}$ となったとする。このとき、次第に平衡状態に戻るとき、すなわち ${x(t)\rightarrow 0\ (t\rightarrow\infty)}$ となるとき、制御対象は「漸近安定」であるという。そのための条件は、行列 ${A}$ が安定行列、すなわち行列 ${A}$ の「すべての固有値の実部が負」であることである。

[3] いま、制御対象が安定ではないとき、これを「状態フィードバック」 ${u(t)=-Fx(t)}$ により安定化したい。すなわち、閉ループ系の行列「 $A-BF$ 」を安定行列としたい。それが可能かどうか調べるためには、「可制御」性の条件「 ${{\rm rank}\ [B,AB,\cdots,A^{n-1}B]=n}$ 」が成立するかどうかをチェックすればよい。

[4] 状態フィードバックを実施するためには、全ての状態変数を、「センサ」を用いて計測する必要がある。しかし、これはいつも可能とは限らないので、状態変数の値を漸近的に推定する「状態オブザーバ」 ${\dot{\hat{x}}(t)=(A-HC)\hat{x}(t)+Hy(t)+Bu(t)}$ が用いられる。これは、 ${\hat{x}(t)\rightarrow x(t)\ (t\rightarrow\infty)}$ となるように、オフザーバゲイン ${H}$ を選択することにより達成される。それが可能かどうか調べるためには、「可観測」性の条件「 ${{\rm rank}\left[\begin{array}{c}C\\CA\\\vdots\\CA^{n-1}\end{array}\right]=n}$ 」が成立するかどうかをチェックすればよい。

[5] 状態フィードバックを状態オブザーバを用いて、 ${u(t)=-F\hat{x}(t)}$ のように実施するとき、状態空間表現 ${\dot{\hat{x}}(t)=(A-HC-BF)\hat{x}(t)+Hy(t)}$ 、 ${u(t)=-F\hat{x}(t)}$ を「オブザーバベース」コントローラと呼ぶ。これによる閉ループ系の固有値の集合は、行列「 $A-BF$ 」の固有値の集合と行列「 $A-HC$ 」の固有値の集合の和集合となることが知られている。したがって、状態フィードバックと状態オブザーバの設計は独立に行うことができる。ただ、制御動作より状態推定が速く行われなければならないので、複素平面上で、行列「 $A-HC$ 」の固有値は行列「 $A-BF$ 」の固有値のより左寄りとすることが考えられる。

【選択肢】 (a) 発展方程式、(b) 状態方程式、(c) ${x(0)\ne0}$ 、(d) ${x(0)=0}$ 、(e)安定、(f)漸近安定、(g)不安定、(h)すべての固有値の実部が負、(i)ある固有値の実部が負、(j)すべての固有値の実部が正、(k)ある固有値の実部が正、(l)状態フィードバック、(m)状態オブザーバ、(n) $A-BF$ 、(o) $A-HC$ 、(p) $A-HC-BF$ 、(q)可制御、(r)可観測、(s)可安定、(t)可検出、(u) ${{\rm rank}\ [B,AB,\cdots,A^{n-1}B]=n}$ 、(v) ${{\rm rank}\left[\begin{array}{c}C\\CA\\\vdots\\CA^{n-1}\end{array}\right]=n}$ 、(w)アクチュエータ、(x)センサ、(y)PID、(z)オブザーバベース

次に、可検出性と可観測性の証明は、とりあえず次を参照してください。

可観測性

CT53 分離定理

投稿日時: 2023年12月30日投稿者: cacsd

Home Work

[M]　数学的には、

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

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

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

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

Flipped Classroom
[1]　１次系

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

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

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

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

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

すなわち

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

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

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

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

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

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

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

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

[3] n次系

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

に対する状態オブザーバ

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

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

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

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

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

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

すなわち

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

ここで、座標変換

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

を行うと

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

したがって

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

CT52　低次元化

投稿日時: 2023年12月30日投稿者: cacsd

Home Work

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

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

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

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

訂正図5.4を出力するProgram 52は、最新のものではありませんでした。最新のp52a.m、p52a.sceを入手してください。図5.4はOKです。

MATLAB

%MATLAB: p52a.m (2024/5/1修正)
Tm=0.01; KE=0.05;
x0=[1;10];
t=0:0.001:0.05;
[t1,x1]=ode45(@(t,x) [0 1; 0 -1/Tm]*x+[0;1/Tm/KE],t,x0);
%-----Part 1
alpha=1.2
L=-(1-alpha)/Tm;
%[t2,x2]=ode45(@(t,x) (-1/Tm-L)*x,t,-x0);
[t2,x2]=ode45(@(t,x) (-1/Tm-L)*x,t,-[-L 1]*x0);
%-----Part 2
alpha=1.8
L=-(1-alpha)/Tm;
%[t3,x3]=ode45(@(t,x) (-1/Tm-L)*x,t,-x0);
[t3,x3]=ode45(@(t,x) (-1/Tm-L)*x,t,-[-L 1]*x0);
%-----
subplot(121), plot(t,x1(:,1)), hold on
subplot(121), plot(t,x1(:,1),'-.')
subplot(121), plot(t,x1(:,1),'--')
grid, axis([0 0.05 0 2])
legend('theta','alpha=1.2','alpha=1.8');
xlabel('t'), ylabel('theta(t)')
subplot(122), plot(t,x1(:,2)), hold on
%subplot(122), plot(t,x1(:,2)+x2(:,2),'-.')
%subplot(122), plot(t,x1(:,2)+x3(:,2),'--')
subplot(122), plot(t,x1(:,2)+x2,'-.')
subplot(122), plot(t,x1(:,2)+x3,'--')
%grid, axis([0 0.05 0 30])
grid, axis([0 0.05 0 80])
legend('omega','alpha=1.2','alpha=1.8');
xlabel('t'), ylabel('omega(t)')

SCILAB

//SCILAB: p52a.sce (2024/5/1修正)
Tm=0.01; KE=0.05;
function dx=f1(t,x),dx=[0 1; 0 -1/Tm]*x+[0;1/Tm/KE], endfunction
x0=[1;10];
t0=0; t=0:0.001:0.05;
x1=ode(x0,t0,t,f1);
//-----Part 1
alpha=1.2
L=-(1-alpha)/Tm;
function dx=f2(t,x),dx=(-1/Tm-L)*x, endfunction
//x2=ode(-x0,t0,t,f2);
x2=ode(-[-L 1]*x0,t0,t,f2);
//-----Part 2
alpha=1.8
L=-(1-alpha)/Tm;
//x3=ode(-x0,t0,t,f2);
x3=ode(-[-L 1]*x0,t0,t,f2);
//-----
subplot(121),plot(t,x1(1,:))
subplot(121),plot(t,x1(1,:),'-.')
subplot(121),plot(t,x1(1,:),'--')
mtlb_grid, mtlb_axis([0 0.05 0 2])
legend(['theta';'alpha=1.2';'alpha=1.8']);
xlabel('t'), ylabel('theta(t)')
subplot(122),plot(t,x1(2,:))
//subplot(122),plot(t,x1(2,:)+x2(2,:),'-.')
//subplot(122),plot(t,x1(2,:)+x3(2,:),'--')
subplot(122),plot(t,x1(2,:)+x2,'-.')
subplot(122),plot(t,x1(2,:)+x3,'--')
//mtlb_grid, mtlb_axis([0 0.05 0 30])
mtlb_grid, mtlb_axis([0 0.05 0 80])
xlabel('t'), ylabel('omega(t)')
legend(['omega';'alpha=1.2';'alpha=1.8']);

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

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

CT51 状態オブザーバ

投稿日時: 2023年12月30日投稿者: cacsd

Home Work 5.1

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

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

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

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

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

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

だから

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

これが右辺

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

に等しいので

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

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

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