制御系CAD

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

6. 定値外乱を抑制する

投稿日時: 2022年8月18日投稿者: cacsd

これまで，安定性の定義，安定化の可能性とその方策を学んできた。その結果，私たちは，閉ループ系の零入力応答を適切に安定化する，すなわち，ある初期状態から出発して速やかに零状態(平衡状態)にもっていくことはできる。この零入力応答は，適当な強さのインパルス入力に対する零状態応答に等しいことを思い出してほしい。

例えば，一定速度で回転しているモータを考え，これに瞬間的な負荷がかかって，その速度が乱されたとする。このときは，状態フィードバックによってもとの速度に復帰させることができる。しかし，モータで何か対象物を回転させるときは，モータ軸への持続的な負荷が生じ，極端な場合はモータは回転することさえできないであろう。同様に，飛行機，船舶，車などのビークルに強い向い風が吹けば定速を維持することは困難となるであろう。

また，モータやビークルでは，速度の設定値を変えることが制御の大切な役割なのに，そのための手段はまったく学んでいない。

したがって，これまでの知識だけでは多くの現実問題を解くことはできない。実際の問題では，必ずといっていいほど，さまざまな「定値外乱」の影響を受ける状況のもとで，「設定値」の保持や，ときには変更を要求されるからである。このような問題を解決するためには，これまでの制御方式に加えて，まず「積分動作」の導入による構造的な改良が必要となる。本章では，定値外乱の影響の除去や，設定値の変更に対応できるように，積分動作を入れて安定化を行う方策について学ぶ。その安定化にLQ制御を適用するとき，その制御方式を簡潔に表すために「LQI制御」と呼ぶ。

6.1 1次系における定値外乱の影響を抑制する
6.1.1 定値外乱の影響

1章で述べた直流モータの状態方程式(1.11)すなわち

$\displaystyle{(1)\quad \dot{\omega}(t)=-\frac{K_TK_E}{JR}\omega(t)+\frac{K_T}{JR}v(t)-\frac{1}{J}\tau_L(t) }$

を思い出そう。これは，二つの入力変数 $v(t),\,\tau_L(t)$ をもち，印加電圧 $v(t)$ は操作できるが，負荷トルク $\tau_L(t)$ は操作できない外乱入力であった。ここで，負荷トルクは時間によらずつねに一定とし，その値は測定できないとする。

以下では，この状況を一般化した次式で表される1次系を考える。

$\displaystyle{(2)\quad \dot{x}(t)=ax(t)+bu(t)+w }$

ここで， $u(t)$ は操作変数， $w$ は定値外乱(constant disturbance)と呼ばれる。そのブロック線図を図6.1に示す。

図6.1 外乱入力をもつ1次系

いま，1次系(2)に対して，状態フィードバック

$\displaystyle{(3)\quad u(t)=-fx(t) }$

を行うと，閉ループ系

$\displaystyle{(4)\quad \dot{x}(t)=(a-bf)x(t)+w }$

を得る。ここで， $a-bf<0$ とする。この時間応答は

$\displaystyle{(5)\quad x(t)=e^{(a-bf)t}x(0)+\int_0^t e^{(a-bf)(t-\tau)}w\,d\tau }$

となる。第2項は

$\displaystyle{(6)\quad w\left[\frac{1}{-(a-bf)}e^{(a-bf)(t-\tau)}\right]_0^t =\frac{w}{-(a-bf)}(1-e^{(a-bf)t}) }$

と計算できるので， $a-bf<0$ を考慮すると次式を得る。

$\displaystyle{(7)\quad x(t)\ \rightarrow\ \frac{w}{-(a-bf)} \quad (t\rightarrow\infty) }$

すなわち，定値外乱が加わると

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

とはならない。このことをシミュレーション例で確認しよう。

シミュレーション 6.1
1次系 $\dot{x}(t)=x(t)+u(t)+w$ の状態フィードバック $u(t)=-2x(t)$ による閉ループ系において， $x(0)=1.5$ ， $w=1$ が加わるときの時間応答を，定値外乱がない場合(破線)とともに図6.2に示す。

図6.2 定値外乱の影響

このような定値外乱の影響を取り除くために，式(3)の代わりに

$\displaystyle{(9)\quad u(t)=-fx(t)+v(t) }$

のような方策を用いる。このときの閉ループ系

$\displaystyle{(10)\quad \dot{x}(t)=(a-bf)x(t)+bv(t)+w }$

の時間応答は

$\displaystyle{(11)\quad x(t)=e^{(a-bf)t}x(0)+\int_0^t e^{(a-bf)(t-\tau)} \underbrace{(bv(t)+w)}_{to be 0} \,d\tau }$

となる。これより，外乱を打ち消すように

$\displaystyle{(12)\quad v(t)=-\frac{1}{b}w }$

と選べば，式(8)を達成できる。しかしながら， $w$ の値は未知としているので，式(12)は実際には使用できない。そこで，どのようにして式(12)の右辺を近似的に作り出すかが問題となる。

図6.2の下に次の問を設定していましたが，紙幅の関係で省略しました。

図6.2で， $x(\infty)=-\displaystyle{\frac{w}{a-bf}}$ であることを確かめなさい。

6.1.2 積分動作を加える

この方法の基礎となるアイデアは

$\displaystyle{(13)\quad x_I(t)=\int_0^tx(\tau)\,d\tau \rightarrow constant \quad (t\rightarrow\infty) }$

を達成できれば $x(t)\rightarrow 0 \ (t\rightarrow\infty)$ が成り立つ，すなわち，積分値が一定値に収束すれば，その被積分項は零に収束するというものである。このように何らかの変数を積分する仕組みのことを，積分動作(integral action)と呼んでいる。

そこで，式(13)を満足する $x_I(t)$ をフィードバックして，定値外乱を打ち消すことが考えられる。すなわち，式(12)の代わりに， $v(t)=-f_Ix_I(t)$ を用いて( $f_I$ はスケーリングファクタ)，結果的に

$\displaystyle{(14)\quad v(t) \rightarrow -\frac{1}{b}w \quad (t\rightarrow\infty) \quad \Leftrightarrow \quad x_I(t) \rightarrow \frac{w}{f_Ib} \quad (t\rightarrow\infty) }$

を達成することができないか検討しよう。

まず，式(2)と，式(13)の $x_I(t)$ の定義から得られる積分器(integrator)

$\displaystyle{(15)\quad \dot{x}_I(t)=x(t) \quad (x_I(0)=0) }$

を合わせて

$\displaystyle{(16)\quad \underbrace{ \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_I(t) \end{array}\right] }_{\dot{x}_E(t)} = %\underbrace{ \left[\begin{array}{cc} a & 0 \\ 1 & 0 \end{array}\right] %}_{A_E} \underbrace{ \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] }_{x_E(t)} + %\underbrace{ \left[\begin{array}{c} b \\ 0 \end{array}\right] %}_{B_E} u(t) + \underbrace{ \left[\begin{array}{c} w \\ 0 \end{array}\right] }_{w_E} }$

を構成する。式(16)に

$\displaystyle{(17)\quad u(t)=-fx(t) %\underbrace{ -f_Ix_I(t) %}_{+v(t)} =- \underbrace{ \left[\begin{array}{cc} f & f_I \end{array}\right] }_{F_E} \underbrace{ \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] }_{x_E(t)} }$

を代入すると，閉ループ系は

$\displaystyle{(18)\quad \underbrace{ \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_I(t) \end{array}\right] }_{\dot{x}_E(t)} = \underbrace{ \left[\begin{array}{cc} a-bf & -bf_I \\ 1 & 0 \end{array}\right] }_{A_{EF}} \underbrace{ \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] }_{x_E(t)} + \underbrace{ \left[\begin{array}{c} w \\ 0 \end{array}\right] }_{w_E} }$

となる。そのブロック線図を図6.3に示す。

図6.3 積分動作を入れる

さて， $A_{EF}$ は安定行列であるとする。式(18)の解は

$\displaystyle{(19)\quad x_E(t)=e^{A_{EF}t}{x_E}(0) +\int_0^t e^{A_{EF}(t-\tau)}{w_E}\,d\tau }$

となる。ここで， $A_{EF}$ が正則であることに注意して，第2項は

$\displaystyle{(20)\quad \left[-e^{A_{EF}(t-\tau)}A_{EF}^{-1}\right]_0^t{w_E} =-(I-e^{A_{EF}t})A_{EF}^{-1}{w_E} }$

と計算できるから

$\displaystyle{(21)\quad x_E(t)=e^{A_{EF}t}{x_E}(0) -(I-e^{A_{EF}t})A_{EF}^{-1}{w_E} }$

と表される。したがって

$\displaystyle{(22)\quad x_E(t)\ \rightarrow\ -A_{EF}^{-1}{w_E}\quad (t\rightarrow\infty) }$

を得る。ここで $A_{EF}^{-1}= \displaystyle{\frac{1}{bf_I}} \left[\begin{array}{cc} 0 & bf_I \\ -1 & a-bf \end{array}\right]$ だから

$\displaystyle{(23)\quad \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] \ \rightarrow\ \left[\begin{array}{c} 0 \\ \displaystyle{\frac{w}{bf_I}} \end{array}\right] \quad (t\rightarrow\infty) }$

となり，外乱の影響が抑制されていることと，式(14)が成り立つことが確認できる。また，式(17)からつぎが成り立つ。

$\displaystyle{(24)\quad u(t) \ \rightarrow\ -\frac{w}{b} \quad (t\rightarrow\infty) }$

シミュレーション6.2
1次系 $\dot{x}(t)=x(t)+u(t)+w$ の状態フィードバック $u(t)=-3x(t)-x_I(t)$ による閉ループ系において， $x(0)=1$ ， $w=1$ が加わるときの時間応答を図6.4に示す(破線は $x_I(t)$ )。

図6.4 定値外乱の抑制

図6.4の下に次の問を設定していましたが，紙幅の関係で省略しました。

図6.4で， $x_I(\infty)=\displaystyle{\frac{w}{bf_I}}$ ， $u(\infty)=-\displaystyle{\frac{w}{b}}$ であることを確かめなさい。

また，(19),(20),(21)におい $e^{A_{EF}t}=\exp(A_{EF}t)$ であることに注意してください。

6.1.3 設定値を変更する

前節では，定値外乱の加わる1次系

$\displaystyle{(25)\quad \dot{x}(t)=ax(t)+bu(t)+w\quad (b\ne0) }$

に対して，制御目的

$\displaystyle{(26)\quad x(t)\ \rightarrow\ 0 \quad (t\rightarrow\infty) }$

を達成する方法を学んだ。ここでは，任意の定数 $r$ に対して，制御目的

$\displaystyle{(27)\quad x(t)\ \rightarrow\ r \quad (t\rightarrow\infty) }$

を達成する問題を考える。この $r$ を設定値(set value)と呼ぶ（参照(reference)値，目標(target)値，指令(command)値と呼ぶこともある）。

図6.5 外乱を抑制したうえで設定値を変更する

この問題の解決のためには，外乱を抑制するための図6.4を，図6.5に示すように変更すればよい。実際，式(27)を $x(t)-r \ \rightarrow\ 0 \quad (t\rightarrow\infty)$ と書き直してみれば

$\displaystyle{(28)\quad \dot{x}_I(t)=x(t)-r \qquad (x_I(0)=0) }$

のような積分動作を導入すればよいことに気づくであろう。

式(25)と式(26)を合わせて

$\displaystyle{(29)\quad %\underbrace{ \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_I(t) \end{array}\right] %}_{\dot{x}_E(t)} = %\underbrace{ \left[\begin{array}{cc} a & 0 \\ 1 & 0 \end{array}\right] %}_{A_E} %\underbrace{ \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] %}_{x_E(t)} + %\underbrace{ \left[\begin{array}{c} b \\ 0 \end{array}\right] %}_{B_E} u(t) + %\underbrace{ \left[\begin{array}{c} w \\ -r \end{array}\right] %}_{w_E} }$

を構成する。いま，式(29)に

$\displaystyle{(30)\quad u(t)=-fx(t)-f_Ix_I(t) =- %\underbrace{ \left[\begin{array}{cc} f & f_I \end{array}\right] %}_{F_E} %\underbrace{ \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] %}_{x_E(t)} }$

を代入すると

$\displaystyle{(31)\quad %\underbrace{ \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_I(t) \end{array}\right] %}_{\dot{x}_E(t)} = \underbrace{ \left[\begin{array}{cc} a-bf & -bf_I \\ 1 & 0 \end{array}\right] }_{A_{EF}} %\underbrace{ \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] %}_{x_E(t)} + %\underbrace{ \left[\begin{array}{c} w \\ -r \end{array}\right] %}_{w_E} }$

となる。 $A_{EF}$ が安定行列のとき，6.1.2項と同様にして

$\displaystyle{(32)\quad \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] \ \rightarrow\ \underbrace{ \displaystyle{\frac{1}{bf_I}} \left[\begin{array}{cc} 0 & bf_I \\ -1 & a-bf \end{array}\right] }_{A_{EF}^{-1}} \left[\begin{array}{c} -w \\ r \end{array}\right] \quad (t\rightarrow\infty) }$

を得る。この第1行は， $x(t)\,\rightarrow\,r\ (t\rightarrow\infty)$ となっていて，制御目的(27)が成り立つ。また，第2行から

$\displaystyle{(33)\quad -f_Ix_I(t)\ \rightarrow\ -\frac{1}{b}w -\frac{a-bf}{b}r \quad (t\rightarrow\infty) }$

を得る。ここで，設定値 $r$ は既知だから，第2項はあらかじめフィードフォワードすることができる。この場合の制御方式はつぎのように表される。

$\displaystyle{(34)\quad u(t)=-fx(t)+f_I\underbrace{\int_0^t(r-x(\tau))\,d\tau}_{-x_I(t)} +\underbrace{\left(-\frac{a-bf}{b}\right)}_{f_r}r }$

また，つぎが成り立つ。

$\displaystyle{(35)\quad u(t) \ \rightarrow\ u_\infty=-\frac{1}{b}(w+ar) \quad (t\rightarrow\infty) }$

それでは，フィードフォワード項による速応性改善の効果をシミュレーションによって確かめよう。

シミュレーション6.3
1次系 $\dot{x}(t)=x(t)+u(t)+w$ の状態フィードバック $u(t)=-3x(t)-x_I(t)$ による閉ループ系において， $x(0)=0$ ， $w=0$ のとき，設定値 $r=1$ へ追従する様子を図6.6，図6.7に示す(破線は $x_I(t)$ )。

図6.6 設定値への追従(フィードフォワードなし)

図6.7 設定値への追従(フィードフォワードあり)

問6.1
図6.6と図6.7で， $x_I(\infty)=\displaystyle{\frac{w}{bf_I}+\frac{a-bf}{bf_I}r}$ ， $u(\infty)=\displaystyle{-\frac{1}{b}(w+ar)}$ であることを確かめなさい。

6.2 積分型コントローラ
6.2.1 積分動作を加えた状態フィードバック

次式で表される可制御かつ可観測で， $m$ 入力 $p$ 出力をもつ $n$ 次系を考える。

$\displaystyle{(36)\quad \dot{x}(t)=Ax(t)+Bu(t)+w }$

$\displaystyle{(37)\quad y(t)=C_Mx(t) }$

ここで，状態方程式には，操作入力 $u(t)$ のほかに，定値外乱 $w$ が加わっていること， $C$ 行列を $C_M$ と書いたことに注意する。いま，出力変数の一部やそれらの組合せからなる新しい $m$ 個の被制御変数(controlled variables)を

$\displaystyle{(38)\quad z(t)=C_Sy(t)=\underbrace{C_SC_M}_{C}x(t) }$

のように選び( $(A,C)$ は可観測対とする)，定値外乱に無関係に，制御目的

$\displaystyle{(39)\quad z(t)\ \rightarrow\ r \quad (t\rightarrow\infty) }$

を達成したい。ここで，定数ベクトル $r$ は $m$ 個の設定値からなる。もし制御目的(39)が物理的に可能とすると，ある状態 $x_\infty$ と入力 $u_\infty$ が確定し

$\displaystyle{(40)\quad \left[\begin{array}{c} -w \\ r \end{array}\right] = \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] \left[\begin{array}{c} x_\infty \\ u_\infty \end{array}\right] }$

の関係を満足しているはずである。したがって，どのような $w$ と $r$ に対しても， $x_\infty$ と $u_\infty$ が定まるように，被制御変数(38)を

$\displaystyle{(41)\quad {\rm rank}\, \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S}=n+m }$

が成り立つように選ぶものとする。

図6.8 積分動作を加えた状態フィードバックによる閉ループ系

さて，制御目的(39)を達成するために，図6.8に示すような，つぎの積分動作を加えた状態フィードバックを考える。

$\displaystyle{(42)\quad u(t)=-Fx(t) -F_I\underbrace{\int_0^t(z(\tau)-r)\,d\tau}_{x_I(t)} }$

ここで，第2項は積分動作を表している。このように $x_I(t)$ を定義すると

$\displaystyle{(43)\quad \dot{x}_I(t)=z(t)-r \qquad (x_I(0)=0) }$

を得る。式(36)と式(43)を合わせて

$\displaystyle{(44)\quad %\underbrace{ \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_I(t) \end{array}\right] %}_{\dot{x}_E(t)} = %\underbrace{ \left[\begin{array}{cc} A & 0 \\ C & 0 \end{array}\right] %}_{A_E} %\underbrace{ \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] %}_{x_E(t)} + %\underbrace{ \left[\begin{array}{c} B \\ 0 \end{array}\right] %}_{B_E} u(t) + %\underbrace{ \left[\begin{array}{c} w \\ -r \end{array}\right] %}_{w_E} }$

を得る。式(44)に，式(42)すなわち

$\displaystyle{(45)\quad u(t)=- %\underbrace{ \left[\begin{array}{cc} F & F_I \end{array}\right] %}_{F_E} %\underbrace{ \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] %}_{x_E(t)} }$

を代入すると，閉ループ系は，つぎのように表される。

$\displaystyle{(46)\quad %\underbrace{ \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_I(t) \end{array}\right] %}_{\dot{x}_E(t)} = \underbrace{ \left[\begin{array}{cc} A-BF & -BF_I \\ C & 0 \end{array}\right] }_{A_{EF}} %\underbrace{ \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] %}_{x_E(t)} + %\underbrace{ \left[\begin{array}{c} w \\ -r \end{array}\right] %}_{w_E} }$

いま， $A_{EF}$ は安定行列であるとする。このとき， $A_{EF}$ は正則であり，つぎのように書ける。

$\displaystyle{(47)\quad A_{EF} = \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S} - \left[\begin{array}{cc} B \\ 0 \end{array}\right] \left[\begin{array}{cc} F & F_I+I_m \end{array}\right] }$

よって， $A_{EF}$ の逆行列は，公式

$(P-XQ^{-1}Y)^{-1}=P^{-1}+P^{-1}X(Q-YP^{-1}X)^{-1}YP^{-1}$

を用いて

$\displaystyle{(48)\quad \begin{array}{ll} A_{EF}^{-1} =& S^{-1}+ S^{-1} \left[\begin{array}{cc} B \\ 0 \end{array}\right] \left(I_m- \left[\begin{array}{cc} F & F_I+I_m \end{array}\right] S^{-1} \left[\begin{array}{cc} B \\ 0 \end{array}\right] \right)^{-1} \\[5mm] &\times \left[\begin{array}{cc} F & F_I+I_m \end{array}\right] S^{-1} \end{array} }$

ここで， $\displaystyle{S^{-1} \left[\begin{array}{cc} B \\ 0 \end{array}\right] = \left[\begin{array}{cc} 0 \\ I_m \end{array}\right] }$ に注意して，整理すると

$\displaystyle{(49)\quad A_{EF}^{-1} = \left[\begin{array}{cc} I_n & 0 \\ -F_I^{-1}F & -F_I^{-1} \end{array}\right] S^{-1} }$

のように計算される。ところで，式(46)から

$\displaystyle{(50)\quad \left[\begin{array}{c} x(t) \\ x_I(t) \end{array}\right] \ \rightarrow\ A_{EF}^{-1} \left[\begin{array}{c} -w \\ r \end{array}\right] \quad (t\rightarrow\infty) }$

を得る。この第1ブロック行は，式(40)より

$\displaystyle{(51)\quad x(t)\ \rightarrow\ \left[\begin{array}{cc} I_n & 0 \end{array}\right] S^{-1} \left[\begin{array}{c} -w \\ r \end{array}\right] =x_\infty \quad (t\rightarrow\infty) }$

となって，式(40)の第2ブロック行から制御目的(39)が成り立つ。また，式(50)の第2ブロック行から

$\displaystyle{(52)\quad -F_Ix_{I}(t)\ \rightarrow\ \left[\begin{array}{cc} F & I_m \end{array}\right] S^{-1} \left[\begin{array}{c} -w \\ r \end{array}\right] = Fx_\infty+u_\infty \quad (t\rightarrow\infty) }$

を得る。ここで，設定値 $r$ は既知だから， $r$ に関係した項 $\left[\begin{array}{cc} F & I_m \end{array}\right] S^{-1} \left[\begin{array}{c} 0 \\ r \end{array}\right]$ をフィードフォワードして，速応性を改善できる。すなわち，制御目的(39)を達成する制御方式は

$\displaystyle{(53)\quad u(t)=-Fx(t) +F_I\underbrace{\int_0^t(r-z(\tau))\,d\tau}_{-x_I(t)} +F_rr }$

のように表され，ここで， $F_r$ はつぎのように決定できる。

$\displaystyle{(54)\quad F_r= \left[\begin{array}{cc} F & I_m \end{array}\right] S^{-1} \left[\begin{array}{c} 0 \\ I_m \end{array}\right] }$

6.2.2 積分動作を加えたオブザーバベース\,コントローラ

式(42)の代わりに，実際には，状態オブザーバ

$\displaystyle{(55)\quad \dot{\hat{x}}(t)=(A-HC_M)\hat{x}(t)+Hy(t)+Bu(t) }$

の出力を用いて

$\displaystyle{(56)\quad u(t)=-F\hat{x}(t)-F_Ix_I(t) }$

を実施することになる。ここでの積分動作を加えたオブザーバベース\,コントローラは，式(56)を式(55)に代入した

$\displaystyle{(57)\quad \dot{\hat{x}}(t)=(A-HC_M-BF)\hat{x}(t)-BF_Ix_I(t)+Hy(t) }$

と，式(43)を合わせて，つぎのように表される。

$\displaystyle{(58)\quad \begin{array}{rl} \left[\begin{array}{c} \dot{\hat{x}}(t) \\ \dot{x}_I(t) \end{array}\right] =& \underbrace{ \left[\begin{array}{cc} A-HC_M-BF & -BF_I \\ 0(m\times n) & 0(m\times m) \end{array}\right] }_{A_K} \left[\begin{array}{c} \hat{x}(t) \\ x_I(t) \end{array}\right]\\[10mm] &+ \underbrace{ \left[\begin{array}{cc} H & 0(n\times m)\\ C_S & -I_m \end{array}\right] }_{B_K} \left[\begin{array}{c} y(t) \\ r \end{array}\right] \end{array} }$

$\displaystyle{(59)\quad u(t)= \underbrace{- \left[\begin{array}{cc} F & F_I \end{array}\right] }_{C_K} \left[\begin{array}{c} \hat{x}(t) \\ x_I(t) \end{array}\right] }$

これによる閉ループ系は，式(56)を式(36)に代入した

$\displaystyle{(60)\quad \dot{x}(t)=Ax(t)-BF\hat{x}(t)-BF_Ix_I(t)+w(t) }$

と，式(43)，(57)を合わせて

$\displaystyle{(61)\quad \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_I(t) \\ \dot{\hat{x}}(t) \end{array}\right] = \left[\begin{array}{ccc} A & -BF_I & -BF \\ C & 0 & 0 \\ HC_M & -BF_I & A-HC_M-BF \end{array}\right] \left[\begin{array}{c} x(t) \\ x_I(t) \\ \hat{x}(t) \end{array}\right] + \left[\begin{array}{c} w \\ -r \\ 0 \end{array}\right] }$

のように表される。そのブロック線図を図6.9に示す。

図6.9積分動作を加えたオブザーバベース\,コントローラによる閉ループ系

いま，閉ループ系(61)に，座標変換

$\displaystyle{(62)\quad \left[\begin{array}{c} x(t) \\ x_I(t) \\ e(t) \end{array}\right] = \left[\begin{array}{ccc} I_n & 0 & 0 \\ 0 & I_m & 0 \\ -I_n & 0 & I_n \end{array}\right] \left[\begin{array}{c} x(t) \\ x_I(t) \\ \hat{x}(t) \end{array}\right] }$

を行えば

$\displaystyle{(63)\quad \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_I(t) \\\hline \dot{e}(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc|c} A-BF & -BF_I & -BF \\ C & 0 & 0 \\\hline 0 & 0 & A-HC_M \end{array}\right] }_{ A_{EF}'= \left[\begin{array}{c|c} A_{EF} & - \left[\begin{array}{cc} B \\ 0 \end{array}\right] F \\[5mm] \hline 0 & \widehat{A} \end{array}\right] } \left[\begin{array}{c} x(t) \\ x_I(t) \\\hline e(t) \end{array}\right] + \left[\begin{array}{c} w \\ -r \\\hline -w \end{array}\right] }$

となり，閉ループ系の固有値は，積分動作を加えた状態フィードバックだけの閉ループ系の固有値と状態オブザーバの固有値からなる。

ここで， $A_{EF}$ は安定行列であるとする。このとき，式(63)より

$\displaystyle{(64)\quad \left[\begin{array}{c} x(t) \\ x_I(t) \\\hline e(t) \end{array}\right] \ \rightarrow\ \underbrace{ \left[\begin{array}{c|c} A_{EF}^{-1} & -A_{EF}^{-1} \left[\begin{array}{cc} B \\ 0 \end{array}\right] F\widehat{A}^{-1} \\[5mm]\hline 0 & \widehat{A}^{-1} \end{array}\right] }_{A_{EF}'\,^{-1}} \left[\begin{array}{c} -w \\ r \\\hline w \end{array}\right] \quad (t\rightarrow\infty) }$

すなわち

$\displaystyle{(65)\quad x(t)\ \rightarrow\ x_\infty \quad (t\rightarrow\infty) }$

$\displaystyle{(66)\quad -F_Ix_{I}(t)\ \rightarrow\ F(x_\infty+e_\infty)+u_\infty \quad (t\rightarrow\infty) }$

$\displaystyle{(67)\quad e(t)\ \rightarrow\ e_\infty=\widehat{A}^{-1}w \quad (t\rightarrow\infty) }$

を得る。したがって，定値外乱が存在するときは状態オブザーバに関して定常偏差が残るにもかかわらず，制御目的(39)が成り立つことがわかる。また，フィードフォーワードゲイン $F_r$ は，式(66)の導出過程から，式(54)と同様に決めてよい。

6.3 LQI設計

5 LQG制御

投稿日時: 2022年8月18日投稿者: cacsd

これまで，閉ループ系(18)，(31)，(46)，(63)において， $A_{EF}$ は安定行列であるとしていた。ここでは，これを満足させるための具体的手段として，前章で学んだLQ制御を使うことを考えたい。LQ制御の議論における閉ループ系は自励系(入力をもたない系)を前提にしていたが，本章における閉ループ系は入力をもつことに注意が必要である。この前提を満足させるために，定常状態との差をとって得られる偏差系(error system)が用いられる。

制御目的(39)が達成されたとき成り立つ式(40)より

$\displaystyle{(68)\quad \left[\begin{array}{c} 0 \\ 0 \end{array}\right] = \left[\begin{array}{cc} A & 0 \\ C & 0 \end{array}\right] \left[\begin{array}{c} x_\infty \\ x_{I\infty} \end{array}\right] + \left[\begin{array}{c} B \\ 0 \end{array}\right] u_\infty + \left[\begin{array}{c} w \\ -r \end{array}\right] }$

を得る( $x_{I\infty}$ は定数ベクトル)。まず，(44)から式(68)を引いて，つぎの偏差系を得る。

偏差系E1：
$\displaystyle{(69)\quad \frac{d}{dt} %\underbrace{ \left[\begin{array}{c} x(t)-x_\infty \\ x_I(t)-x_{I\infty} \end{array}\right] %}_{x_E(t)-x_{E\infty}} = \underbrace{ \left[\begin{array}{cc} A & 0 \\ C & 0 \end{array}\right] }_{A_{E1}} %\underbrace{ \left[\begin{array}{c} x(t)-x_\infty \\ x_I(t)-x_{I\infty} \end{array}\right] %}_{x_E(t)-x_{E\infty}} + \underbrace{ \left[\begin{array}{c} B \\ 0 \end{array}\right] }_{B_{E1}} (u(t)-u_\infty) }$

この両辺を微分すれば，状態変数の中の定数ベクトルを除くことができて

偏差系E2：
$\displaystyle{(70)\quad \frac{d}{dt} %\underbrace{ \left[\begin{array}{c} {\dot x}(t) \\ z(t)-r \end{array}\right] %}_{\dot{x}_E(t)} = \underbrace{ \left[\begin{array}{cc} A & 0 \\ C & 0 \end{array}\right] }_{A_{E2}} %\underbrace{ \left[\begin{array}{c} {\dot x}(t) \\ z(t)-r \end{array}\right] %}_{\dot{x}_E(t)} + \underbrace{ \left[\begin{array}{c} B \\ 0 \end{array}\right] }_{B_{E2}} {\dot u}(t) }$

を得る。さらに，式(36)と(38)をまとめた

$\displaystyle{(71)\quad \left[\begin{array}{c} {\dot x}(t)-w \\ z(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S} \left[\begin{array}{c} x(t) \\ u(t) \end{array}\right] }$

から式(40)を引いて，つぎの関係式が成り立つ。

$\displaystyle{(72)\quad \left[\begin{array}{c} {\dot x}(t) \\ z(t)-r \end{array}\right] = \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S} \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] }$

これに基づいて，偏差系E2に座標変換を行えば

偏差系E3：
$\displaystyle{(73)\quad \frac{d}{dt} %\underbrace{ \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] %}_{{\tilde x}_E(t)-{\tilde x}_{E\infty}} = \underbrace{ \left[\begin{array}{cc} A & B \\ 0 & 0 \end{array}\right] }_{A_{E3}} %\underbrace{ \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] %}_{{\tilde x}_E(t)-{\tilde x}_{E\infty}} + \underbrace{ \left[\begin{array}{c} 0 \\ I_m \end{array}\right] }_{B_{E3}} {\dot u}(t) }$

を得る。ここで，つぎの関係式を用いた。

$\displaystyle{(74)\quad \underbrace{ \left[\begin{array}{cc} A & 0 \\ C & 0 \end{array}\right] }_{A_{E2}} \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S} = \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{cc} A & B \\ 0 & 0 \end{array}\right] }_{A_{E3}} }$

$\displaystyle{(75)\quad \underbrace{ \left[\begin{array}{c} B \\ 0 \end{array}\right] }_{B_{E2}} = \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S} \underbrace{ \left[\begin{array}{c} 0 \\ I_m \end{array}\right] }_{B_{E3}} }$

$\displaystyle{(76)\quad \underbrace{ \left[\begin{array}{cc} 0 & I_m \end{array}\right] }_{C_{E2}} \underbrace{ \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] }_{S} = \underbrace{ \left[\begin{array}{cc} C & 0 \end{array}\right] }_{C_{E3}} }$

問6.2
$(A_{E3},B_{E3})$ は可制御対であることを示しなさい。

問6.3
$(A_{E3},C_{E3})$ は可観測対であることを示しなさい。

例えば，偏差系E1に対して，状態フィードバック

$\displaystyle{(77)\quad u(t)-u_\infty=- \left[\begin{array}{cc} F & F_I \end{array}\right] \left[\begin{array}{c} x(t)-x_\infty \\ x_I(t)-x_{I\infty} \end{array}\right] }$

を用いて，閉ループ系を安定化すれば， $A_{EF}= \left[\begin{array}{cc} A-BF & -BF_I \\ C & 0 \end{array}\right]$ を安定行列とすることができる。これは，偏差系E2に対しても同様である。

以下に，偏差系E3をLQ制御により安定化して，積分動作を加えたオブザーバベース\,コントローラを構成する手順を示す。

アルゴリズム6.1　＜LQI設計＞

入力パラメータ1
$A(n\times n),\,B(n\times m),\,C_M(p\times n)$
ただし， $m\le p$ ， $(A,B)$ は可制御対， $(A,C_M)$ は可観測対

出力パラメータ
$A_K(n+m\times n+m),\,B_K(n+m\times p+m),\,C_K(m\times n+m)$

ステップ1　被制御変数の決定

$\left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right]$ が正則となるように( $C=C_SC_M$ )，セレクタ行列 $C_S(m\times p)$ を決める（一般に，多入力多出力系の場合，どの操作変数でどの被制御変数を制御するのかについて，物理的に実現可能な1対1対応を考えることが重要である。その際，被制御変数はフィードバックされるので観測量の中から選ばれなけばならない）。

ステップ2　偏差系の安定化

偏差系

$\displaystyle{(78)\quad \frac{d}{dt} %\underbrace{ \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] %}_{{\tilde x}_E(t)-{\tilde x}_{E\infty}} = \underbrace{ \left[\begin{array}{cc} A & B \\ 0 & 0 \end{array}\right] }_{A_{E}} %\underbrace{ \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] %}_{{\tilde x}_E(t)-{\tilde x}_{E\infty}} + \underbrace{ \left[\begin{array}{c} 0 \\ I_m \end{array}\right] }_{B_{E}} {\dot u}(t) }$

を，状態フィードバック

$\displaystyle{(79)\quad {\dot u}(t)=- \left[\begin{array}{cc} K & K_I \end{array}\right] \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] }$

によるLQ制御で安定化する。その際，評価関数としては

$\displaystyle{(80)\quad \int_0^\infty \left( \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right]^T Q_E \left[\begin{array}{c} x(t)-x_\infty \\ u(t)-u_\infty \end{array}\right] +{\dot u}^T(t)R_E{\dot u}(t)\right)\,dt }$

を用いる。ただし， $(A_E,Q_E^{\frac{1}{2}})$ は可観測対とする。

さらに， $F$ と $F_I$ を，次式から計算する。

$\displaystyle{(81)\quad \left[\begin{array}{cc} F & F_I \end{array}\right] = \left[\begin{array}{cc} K & K_I \end{array}\right] \left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right]^{-1} }$

ステップ3　オブザーバゲイン $H$ の決定

例えば，アルゴリズム5.2のステップ2を用いる。
行列 $B'$ を， $(A,B')$ が可制御対となるように選び，重み行列 $W>0$ と $V>0$ を指定し

$\displaystyle{(82)\quad \Gamma A^T+A\Gamma-\Gamma C_M^TV^{-1}C_M\Gamma+B'WB'^T=0 }$

を解いて， $\Gamma>0$ を求め，オブザーバゲイン $H$ をつぎのように定める。

$\displaystyle{(83)\quad H=V^{-1}C_M\Gamma }$

ステップ4　LQIコントローラの構成

$C_S$ ， $F$ ， $F_I$ ， $H$ から，つぎを構成する。

$\displaystyle{(84)\quad \dot{x}_K(t)=A_Kx_K(t)+B_K \left[\begin{array}{c} y(t) \\ r \end{array}\right] }$

$\displaystyle{(85)\quad u(t)=C_Kx_K(t) }$

ただし

$\displaystyle{(86)\quad A_K= \left[\begin{array}{cc} A-HC_M-BF & -BF_I \\ 0\,(m\times n) & 0\,(m\times m) \end{array}\right] }$

$\displaystyle{(87)\quad B_K= \left[\begin{array}{cc} H & 0\,(n\times m)\\ C_S & -I_m \end{array}\right] }$

$\displaystyle{(88)\quad C_K=- \left[\begin{array}{cc} F & F_I \end{array}\right] }$

この積分動作を加えたオブザーバベース\,コントローラを簡潔にLQIコントローラ}(LQI controller)，これによる制御方式をLQI制御(LQ control with integral action)と呼ぶ。

問6.4
式(81)の妥当性について説明しなさい。

【本章のねらい】
・１次系，２次系の状態フィードバックの最適設計を，リッカチ方程式を解いて行う。
・リッカチ方程式を解くMファイルを理解し，これを用いて高次系の状態フィードバックの最適設計を行なう。
・オブザーバベースト・コントローラの最適設計を，リッカチ方程式CAREを解いて状態フィードバックを，またリッカチ方程式FAREを解いてオブザーバゲインを求めて行う。

5.1 状態フィードバックの最適設計

可制御な $m$ 入力 $p$ 出力 $n$ 次線形系（ $n$ 次系）

$\displaystyle{(5.1)\quad \left\{\begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t) \\ y(t)=Cx(t) \end{array}\right. }$

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

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

の決定法を改めて考える。一つの方法は，閉ループ系

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

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

$\displaystyle{(5.4)\quad \int_0^\infty (x^T(t)C^TQCx(t)+u^T(t)Ru(t))\,dt }$

を設定し，これを最小化する問題を解くことである。ただし， $Q>0$ ， $R>0$ （ $Q$ ， $R$ は正定行列。テキスト「線形システム制御入門」の2.2.3節を参照。），
また， $(A,C)$ は可観測対とする。これによる状態フィードバックのゲイン行列は，リッカチ方程式

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

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

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

まず，１次系の場合を考える。

例題5.1 時定数 $T$ と定常ゲイン $K$ をもつ１次系

$\displaystyle{ \dot{x}(t)=-\frac{1}{T}x(t)+\frac{K}{T}u(t) }$

に対して，評価関数

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

を最小にするように，状態フィードバック $u=-fx$ を決定せよ。
解答　リッカチ方程式

$\displaystyle{ \frac{1}{r^2}(\frac{K}{T})^2\Pi^2-2(-\frac{1}{T})\Pi-q^2=0 }$

の解 $\Pi$ を求めると

$\displaystyle{ \Pi=\frac{(-\frac{1}{T})\pm\sqrt{(-\frac{1}{T})^2-\frac{1}{r^2}(\frac{K}{T})^2(-q^2)}}{\frac{1}{r^2}(\frac{K}{T})^2} =\frac{r^2}{K^2}T\left(-1\pm\sqrt{1+\left(\frac{q}{r}\right)^2K^2}\right) }$

$\Pi>0$ より，求めるゲイン $f$ は

$\displaystyle{ f=\frac{1}{r^2}\frac{K}{T}\Pi=\frac{1}{K}\left(-1+\sqrt{1+\left(\frac{q}{r}\right)^2K^2}\right) }$

演習5.1　例題5.1において， $T=1$ ， $K=1$ とする。このとき，つぎの重み係数をもつ評価関数を最小にする $f$ を決定せよ。
(1) $q=1, r=1$ ，
(2) $q=1, r=0.5$ ，
(3) $q=1, r=0.1$
また， $x(0)=1$ のとき，閉ループ系の時間応答をMATLABでシミュレーションせよ。

つぎに，２次系の場合を考える。

例題5.2　２次系

$\displaystyle{ \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{ \int_0^\infty (x_1^2(t)+x_2^2(t)+u^2(t))\,dt }$

を最小にするように求めよ。
解答　リッカチ方程式

$\displaystyle{ \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{ \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{ \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} (1) &○\\ (2) &× \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} (3) &×\\ (4) &× \end{array}\right. \end{array}\right. }$

ここで，(1)だけが， $\pi_1,\pi_2>0,\ \pi_1\pi_2-\pi_3^2>0$ を満たす
（ $\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{ \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] }$

演習5.2　つぎの２次系について，例題5.2と同じ問題設定で解け。

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

例題5.3　２次系

$\displaystyle{ \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{ \int_0^\infty (q^2y^2(t)+r^2u^2(t))\,dt }$

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

$\displaystyle{ 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{ \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] \nonumber \end{array} }$

を要素ごとに整理して

$\displaystyle{ \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{ \pi_3=\pm r\omega_n^{-2}q }$

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

$\displaystyle{ \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{ \pi_2=r^2\omega_n^{-3}(-2\zeta+\sqrt{4\zeta^2\pm 2r^{-1}q}) }$

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

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

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

$\displaystyle{ \pi_1=r\omega_n^{-1}q\sqrt{2r^{-1}q+4\zeta^2} }$
$\displaystyle{ \pi_2=r^2\omega_n^{-3}(-2\zeta+\sqrt{2r^{-1}q+4\zeta^2}) }$
$\displaystyle{ \pi_3=r\omega_n^{-2}q }$

でなければならない。このとき $\pi_1\pi_2-\pi_3^2>0$ も満足される。したがって

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

演習5.3　２次系

$\displaystyle{ %\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) }$

について，例題5.3と同じ問題設定で解くと

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

のように与えられることを示せ。

ちなみに，例題5.3において，評価関数

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

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

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

となる。演習5.3においては

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

となる。これらから重み係数 $q_1$ と $q_2$ の与え方の指針を得ることができる（テキスト「線形システム制御入門」の第５章演習問題【1】参照。また(??)，(??)の導出についてはそれぞれ第５章例題5.1と演習問題【2】参照。）。

さて，リッカチ方程式を計算機を用いて解くときは，ハミルトン行列 $M$ の安定固有値と対応する固有ベクトルを

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

のように求めて，状態フィードバックゲインを

$\displaystyle{(5.11)\quad F=R^{-1}B^T\underbrace{V_2V_1^{-1}}_{\Pi} }$

によって得る。

例題5.4　リッカチ方程式 $\Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0$ を解き，状態フィードバックゲイン $F=R^{-1}B^T\Pi$ と閉ループ系の固有値 $\lambda(A-BF)$ を求めるMファイル

[F,p]=opt(A,B,C,Q,R)

を準備し，例題5.2を解け。
解答 Mファイル opt.m の作成例をつぎに示す。

%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;

ここで，４行目でハミルトン行列の固有値問題を解き，安定固有値の添字を５行目の index に得て，リッカチ方程式の解 X を得ている。

例題5.2を解くには，MATLABにつぎのコマンドを与えればよい。

%lqr.m
A=[0 1;0 0]; B=[0; 1]; C=[1 0]; Q=1; R=1;
[F,p]=opt(A,B,C,Q,R);

また， p が閉ループ系の固有値 $\lambda(A-BF)$ に等しいことは，

eig(A-B*F)

としてみれば確認できる。

演習5.4　演習5.2を，例題5.4にならって，MATLABで解け。

演習5.5　３次系

$\displaystyle{ \left\{\begin{array}{l} %\underbrace{ \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \\ \dot{x}_3(t) \end{array}\right] %}_{\dot{x}(t)} = \left[\begin{array}{cccc} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \end{array}\right] %\underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \\ x_3(t) \end{array}\right] %}_{x(t)} + \left[\begin{array}{cc} 0 & 0 \\ 1 & -1 \\ 0 & 1 \\ \end{array}\right] \left[\begin{array}{c} u_1(t) \\ u_2(t) \end{array}\right] \\ \left[\begin{array}{c} y_1(t) \\ y_2(t) \end{array}\right] = \left[\begin{array}{cccc} 1 & 0 & 0 \\ 0 & 0 & 1 \\ \end{array}\right] %\underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \\ x_3(t) \end{array}\right] %}_{x(t)} \end{array}\right. }$

を安定化する状態フィードバック $u(t)=-f_1x_1(t)-f_2x_2(t)-f_3x_3(t)$ を，つぎの評価関数を最小にするように求めよ。

$\displaystyle{ J=\int_0^\infty(y_1^2(t)+y_2^2(t)+u_1^2(t)+u_2^2(t))\,dt }$

演習5.6　４次系

$\displaystyle{ \left\{\begin{array}{l} %\underbrace{ \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \\ \dot{x}_3(t) \\ \dot{x}_4(t) \end{array}\right] %}_{\dot{x}(t)} = \left[\begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \end{array}\right] %\underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \\ x_3(t) \\ x_4(t) \end{array}\right] %}_{x(t)} + \left[\begin{array}{c} 0 \\ 0 \\ 1 \\ 0 \end{array}\right] u(t) \\ y(t) = \left[\begin{array}{cccc} 0 & 1 & 0 & 0 \end{array}\right] %\underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \\ x_3(t) \\ x_4(t) \end{array}\right] %}_{x(t)} \end{array}\right. }$

を安定化する状態フィードバック $u(t)=-f_1x_1(t)-f_2x_2(t)-f_3x_3(t)-f_4x_4(t)$ を，つぎの評価関数を最小にするように求めよ。

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

5.2 オブザーバベースト・コントローラの最適設計

オブザーバベースト・コントローラによる閉ループ系の時間応答に関する評価するために，つぎのブロック線図を考える。

オブザーバベースト・コントローラによる閉ループ系の評価

ここで，新しい入力 $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$ の平方根行列により重み付けられている。

この $n$ 次系

$\displaystyle{(5.12)\quad \left\{\begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t)+B'W^{\frac{1}{2}}w(t) \\ y(t)=Cx(t)+V^{\frac{1}{2}}v(t)\\ z(t)=C'x(t) \end{array}\right. }$

に対して，オブザーバベースト・コントローラ

$\displaystyle{(5.13)\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{(5.14)\quad \left\{\begin{array}{rcl} \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_{CL}} \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_{CL}} \left[\begin{array}{cc} {w}(t) \\ {v}(t) \end{array}\right] \\ \left[\begin{array}{cc} {z'}(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_{CL}} \left[\begin{array}{cc} {x}(t) \\ \hat{x}(t) \end{array}\right] \end{array}\right. }$

このインパルス応答行列

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

に関する評価規範として

$\displaystyle{(5.16)\quad J=\int_0^\infty{\rm tr}(G_{CL}^T(t)G_{CL}(t))\,dt }$

を設定し，これを最小化する問題を解くことができる（ ${\rm tr}$ は行列のトレース）。その詳細は割愛するが，
状態フィードバックのゲイン行列 $F$ と状態オブザーバのゲイン行列 $H$ がつぎの手順で決定できる
（この手順はＬＱＧ(Linear Quadratic Gaussian)制御問題の解法と同等である。）。

Step 1　行列方程式

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

の解 $\Pi>0$ を求めて，つぎの状態フィードバックのゲイン行列 $F$ を定める。

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

Step 2　行列方程式

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

の解 $\Gamma>0$ を求めて，つぎの状態オブザーバのゲイン行列 $H$ を定める。

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

例題5.5　１次系 $\dot{x}(t)=ax(t)+bu(t)\,(a=0,b=1)$ に対するオブザーバベースト・コントローラを，

$\displaystyle{ {\bf CARE}:\ -r^{-2}b^2\Pi^2+2a\Pi+q^2=0 }$
$\displaystyle{ {\bf FARE}:\ -\rho^{-2}c^2\Gamma^2+2a\Gamma+\sigma^2=0 }$

において $q=1$ ， $r=1$ ， $\sigma=1$ ， $\rho=1$ と選んで構成せよ。
解答　CAREに， $a=0$ ， $b=1$ ， $q=1$ ， $r=1$ を代入して， $\Pi^2=1$ 。 $\Pi>0$ より， $\Pi=1$ 。
したがって，状態フィードバックゲインは， $f=\frac{1}{r^2}b\Pi=1$ 。FAREに， $a=0$ ， $c=1$ ， $\sigma=1$ ， $\rho=1$ を代入して，
$\Gamma^2=1$ 。 $\Gamma>0$ より， $\Gamma=1$ 。したがって，オブザーバゲインは， $h=\frac{1}{\rho^2}c\Gamma=1$ 。
以上から，オブザーバベースト・コントローラが，つぎのように得られる。

$\displaystyle{ \left\{\begin{array}{l} \dot{\hat{x}}(t)=(a-hc-bf)\hat{x}(t)+hy(t)=-2\hat{x}(t)+y(t) \\ u(t)=-f\hat{x}(t) =-\hat{x}(t) \end{array}\right. }$

演習5.7　１次系 $\dot{x}(t)=ax(t)+bu(t)\,(a=-1,b=1)$ に対するオブザーバベースト・コントローラを，
例題5.5のCAREにおいて $q=1$ ， $r=1$ ，FAREにおいて $\sigma=1$ ， $\rho=0.2$ と選んで構成せよ。

例題5.6　ＬＱＧ制御則を設計するためのMファイルを作成せよ。

[AK,BK,CK,pK,pcare,pfare] = optobc(A,B,C,CC,Q,R,BB,W,V)

解答　例題5.4で作成したMファイル opt.m を用いる。

%optobc.m
function [AK,BK,CK,pK,pcare,pfare] = optobc(A,B,C,CC,Q,R,BB,W,V)
[F,pcare]=opt(A,B,CC,Q,R);
[H,pfare]=opt(A’,C’,BB’,W,V); H=H’;
AK=A-H*C-B*F; BK=H; CK=F; pK=eig(AK);

演習5.8　演習5.7を，例題5.6で作成したMファイル optobc.m を用いて解け。

演習問題の解答

【演習5.1】
(1) $f=-1+\sqrt{2}$ ，
(2) $f=-1+\sqrt{3}$ ，
(3) $f=-1+\sqrt{11}$ 。

MATLABによるシミュレーションはつぎのように行えばよい。

%lqr1.m
T=1; K=1; a=-1/T; b=K/T; c=1; t=0:0.1:5; x0=1;
s0=ss(a,b,c,0); y0=initial(s0,x0,t);
f1=-1+sqrt(2); s1=ss(a-b*f1,b,c,0); y1=initial(s1,x0,t); u1=-f1*y1;
f2=-1+sqrt(3); s2=ss(a-b*f2,b,c,0); y2=initial(s2,x0,t); u2=-f2*y2;
f3=-1+sqrt(11); s3=ss(a-b*f3,b,c,0); y3=initial(s3,x0,t); u3=-f3*y3;
figure(1),subplot(121),plot(t,y0,t,y1,t,y2,t,y3),grid,title(‘y’)
figure(1),subplot(122),plot(t,u1,t,u2,t,u3),grid,title(‘u’)

【演習5.2】
(1)リッカチ方程式

$\begin{array}{ll} \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ 0 & -1 \end{array}\right] + \left[\begin{array}{cc} 0 & 0 \\ 1 & -1 \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] \\ \times \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}{c} 1 \\ 0 \end{array}\right] \left[\begin{array}{cc} 1 & 0 \end{array}\right] = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \end{array}$

を要素ごとに整理して

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

を得る。 $\pi_1>0,\ \pi_1\pi_2-\pi_3^2>0$ を満たす解は， $\pi_1= \sqrt{3},\,\pi_2=-1+\sqrt{3},\,\pi_3=1$ 。
したがって

$F= \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 & -1+\sqrt{3} \end{array}\right]$
%[f,p]=opt([0 1;0 -1],[0;1],[1 0],1,1)

(2)リッカチ方程式

$\begin{array}{ll} \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ -1 & 0 \end{array}\right] + \left[\begin{array}{cc} 0 & -1 \\ 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] \\ \times \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}{c} 1 \\ 0 \end{array}\right] \left[\begin{array}{cc} 1 & 0 \end{array}\right] = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \end{array}$

を要素ごとに整理して

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

を得る。 $\pi_1>0,\ \pi_1\pi_2-\pi_3^2>0$ を満たす解は， $\pi_1= \sqrt{-4+2\sqrt{2}},\,\pi_2=\sqrt{-2+\sqrt{2}},\,\pi_3=-1+\sqrt{2}$ 。したがって

$F= \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{2} & \sqrt{-2+\sqrt{2}} \end{array}\right]$
%[f,p]=opt([0 1;-1 0],[0;1],[1 0],1,1)

【演習5.3】
リッカチ方程式

$\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] \\ \times \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}$

を要素ごとに整理して

$\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_1>0,\ \pi_1\pi_2-\pi_3^2>0$ を満たすものを選ぶと

$\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.$

を得る。これらを用いて， $f_1$ と $f_2$ の表現式が，例題??と同様にして得られる。

【演習5.4】
%lqr2.m
A=[0 1;0 -1]; B=[0;1]; C=[1 0]; Q=1; R=1;
[F,p]=opt(A,B,C,Q,R);
sys=ss(A-B*F,B,eye(2,2),0); x0=[1;0];
t=0:0.1:10; x=initial(sys,x0,t)’; y=C*x; u=-F*x;
figure(2),subplot(121),plot(t,y),grid,title(‘y’)
figure(2),subplot(122),plot(t,u),grid,title(‘u’)
(2)の場合，A=[0 1;-1 0] と書き換える。

【演習5.5】
%lqr3.m
A=[0 1 0;0 0 0;0 0 -1]; B=[0 0;1 -1;0 1]; C=[1 0 0;0 0 1];
[F,p]=opt(A,B,C,1,1)

【演習5.6】
%lqr4.m
A=[0 0 1 0;0 0 0 1;-1 1 0 0;1 -1 0 0]; B=[0;0;1;0]; C=[0 1 0 0];
[F,p]=opt(A,B,C,1,1)

【演習5.7】
CARE: $-\frac{1}{r^2}b^2\Pi^2+2a\Pi+q^2=0$ に， $a=-1$ ， $b=1$ ， $q=1$ ， $r=1$ を代入して， $\Pi^2+2\Pi-1=0$ 。 $\Pi>0$ より $\Pi=-1+\sqrt{2}$ 。
したがって，状態フィードバックゲインは， $f=\frac{1}{r^2}b\Pi=-1+\sqrt{2}$ 。

FARE: $-\frac{1}{\rho^2}c^2\Gamma^2+2a\Gamma+\sigma^2=0$ に， $a=-1$ ， $c=1$ ， $\sigma=1$ ， $\rho=0.2$ を代入して， $25\Gamma^2+2\Gamma-1=0$ 。 $\Gamma>0$ より， $\Gamma=\frac{1}{25}(-1+\sqrt{26})$ 。
したがって，オブザーバゲインは， $h=\frac{1}{\rho^2}c\Gamma=-1+\sqrt{26}$ 。

以上から，オブザーバベースト・コントローラが，つぎのように得られる。

$\left\{\begin{array}{l} \dot{\hat{x}}(t)=(a-hc-bf)\hat{x}(t)+hy(t)=(1-\sqrt{2}-\sqrt{26})\hat{x}(t)+(-1+\sqrt{26})y \\ u(t)=-f\hat{x}(t)=-(-1+\sqrt{2})\hat{x}(t) \end{array}\right.$

【演習5.8】
%lqr5.m
[A

4 状態オブザーバと可観測性

投稿日時: 2022年8月18日投稿者: cacsd

【本章のねらい】
・状態オブザーバを構成する。
・可観測性と可検出性を判定する。

4.1 状態オブザーバ

いま制御対象は平衡状態にあるとし，何らかの要因でこれが乱されたとき，速やかに元の平衡状態に戻す手段として， $m$ 入力 $p$ 出力 $n$ 次元線形系（ $n$ 次系）

$\displaystyle{(4.1)\quad \left\{\begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t) \\ y(t)=Cx(t) \end{array}\right. }$

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

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

を考えた。しかしながら，現実には状態変数をすべて計測できる場合は少ない。したがって，状態フィードバックは実際には実施できるとは限らない。そこで，状態オブザーバとよばれる $n$ 次系

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

が考案されている。ここで，サイズ $n\times p$ の行列 $H$ が設計パラメータである。このブロック線図をつぎに示す。

実際，(4.3)から(4.1)の第１式を辺々引き算すると

$\displaystyle{(4.4)\quad \underbrace{\dot{\hat{x}}(t)-\dot{x}(t)}_{\dot{e}(t)}=(A-HC)\underbrace{(\hat{x}(t)-x(t))}_{e(t)} }$

ここで，行列 $A-HC$ が安定行列であれば

$\displaystyle{(4.5)\quad for\ any\ e(0)\ne0,\ e(t)\rightarrow 0\quad(t\rightarrow\infty) }$

となって， $n$ 次系(4.3)の状態 $\hat{x}(t)$ が $n$ 次系(4.1)の状態 ${x}(t)$ に漸近していく。したがって，行列 $A-HC$ が安定行列となるように状態オブザーバのゲイン行列 $H$ をどう求めるか問題となる。

一つのアプローチは，つぎの仮想的な $n$ 次系

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

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

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

を求めることである。実際，閉ループ系は

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

となって，行列 $(A-HC)^T$ を安定行列，よって行列 $A-HC$ を安定行列とすることができる。

したがって，前章の状態フィードバックの設計法をそのまま援用できるが，実際には次章のLQG制御問題として解く場合が多い。

例題4.1　２次系

$\displaystyle{ \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 & 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}{cc} 0 \\ 1 \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] }_{x(t)} \end{array}\right. }$

に対する状態オブザーバ $\dot{\hat{x}}(t)=(A-HC)\hat{x}(t)+Hy(t)+Bu(t)$ を，行列 $A-HC$ の固有値が $-1,-2$ となるように構成せよ。
解答　行列 $A$ の特性多項式は

$\displaystyle{ {\rm det}(\lambda I_2-A)= {\rm det}\left[\begin{array}{cc} \lambda & -1 \\ 0 & \lambda \end{array}\right] =\lambda^2+\underbrace{0}_{a_1}\lambda+\underbrace{0}_{a_2} }$

行列 $A-HC$ の特性多項式は

$\displaystyle{ {\rm det}(\lambda I_2-A+HC)=(\lambda-(-1))(\lambda-(-2))= \lambda^2+\underbrace{3}_{a'_1}\lambda+\underbrace{2}_{a'_2} }$

これらから，オブザーバゲイン $H$ は，つぎのように計算される。

$\displaystyle{ H^T= \left[\begin{array}{cc} 2-0 & 3-0 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right]^{-1} \left[\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right]^{-1} =\left[\begin{array}{cc} 3 & 2 \end{array}\right] }$

したがって，求める状態オブザーバ $\dot{\hat{x}}=(A-HC)\hat{x}+Hx+Bu$ は

$\displaystyle{ \left[\begin{array}{c} \dot{\hat{x}}_1 \\ \dot{\hat{x}}_2 \end{array}\right] = ( \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right]- \left[\begin{array}{c} 3 \\ 2 \end{array}\right] \left[\begin{array}{cc} 1 & 0 \end{array}\right] ) \left[\begin{array}{c} \hat{x}_1 \\ \hat{x}_2 \end{array}\right] + \left[\begin{array}{c} 3 \\ 2 \end{array}\right] y + \left[\begin{array}{c} 0 \\ 1 \end{array}\right] u }$

$\displaystyle{ = \left[\begin{array}{cc} -3 & 1 \\ -2 & 0 \end{array}\right] \left[\begin{array}{c} \hat{x}_1 \\ \hat{x}_2 \end{array}\right] + \left[\begin{array}{c} 3 \\ 2 \end{array}\right] y + \left[\begin{array}{c} 0 \\ 1 \end{array}\right] u }$

演習4.1　２次系

$\displaystyle{ \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 & 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}{cc} 0 \\ 1 \end{array}\right] }_{B} u(t)\\ y(t)= \underbrace{ \left[\begin{array}{cc} 1 & 1 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \end{array}\right] }_{x(t)} \end{array}\right. }$

に対する状態オブザーバ $\dot{\hat{x}}(t)=(A-HC)\hat{x}(t)+Hy(t)+Bu(t)$ を，行列 $A-HC$ の固有値が $-2,-2$ となるように構成せよ。

例題4.2　例題4.1において， $x_1(0)=1,x_2(0)=0.5$ のときの零入力応答 $x(t)$ を，状態オブザーバの出力 $\hat{x}(t)=x(t)+e(t)$ が追従する様子をMATLABでシミュレーションせよ。
解答 $x(t)=\exp(At)x(0)$ と $\hat{x}(t)=x(t)+\exp((A-HC)t)(\hat{x}(0)-x(0))$ を重ねて描くためのMファイルは，たとえばつぎのように書ける。

%obs_err.m
A=[0 1;0 0]; C=[1 0]; H=[3;2];
sys1=ss(A,zeros(2,1),eye(2,2),0); x0=[1;0.5];
sys2=ss(A-H*C,zeros(2,1),eye(2,2),0); xh0=[0;0];
t=0:0.1:5; x=initial(sys1,x0,t); e=initial(sys2,xh0-x0,t); xh=x+e;
figure(1),subplot(121),plot(t,x(:,1),t,xh(:,1)),axis([0 5 0 4]),grid
figure(1),subplot(122),plot(t,x(:,2),t,xh(:,2)),axis([0 5 0 2]),grid

演習4.2　演習4.1において， $x_1(0)=1,x_2(0)=0.5$ のときの零入力応答 $x(t)$ を，状態オブザーバの出力 $\hat{x}(t)=x(t)+e(t)$ が追従する様子をMATLABでシミュレーションせよ。

さて， $n$ 次系(4.1)に対する状態フィードバック(4.2)を，状態オブザーバ(4.3)の出力(＝状態)を用いて

$\displaystyle{(4.9)\quad u(t)=-F\hat{x}(t) }$

のように実施するとき，オブザーバベースト・コントローラは $n$ 次系

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

となる。このとき，閉ループ系はつぎのように表される( $e(t)={\hat x}(t)-x(t)$ )。

$\displaystyle{(4.11)\quad \left[\begin{array}{c} \dot{x}(t) \\ \dot{e}(t) \end{array}\right]= \left[\begin{array}{ccc} A-BF & -BF \\ 0 & A-HC \end{array}\right] \left[\begin{array}{c} x(t) \\ e(t) \end{array}\right] }$

このように閉ループ系の $A$ 行列の固有値の集合は，状態フィードバックによる $A-BF$ の固有値の集合と状態オブザーバによる $A-HC$ の固有値の集合の和となる（テキスト「線形システム制御入門」の 4.4節参照）。

例題4.3　１次系

$\displaystyle{ \left\{\begin{array}{l} \dot{x}(t)=ax(t)+bu(t) \\ y(t)=cx(t) \end{array}\right. }$

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

$\displaystyle{ u(t)=-fx(t) }$

と，状態オブザーバ

$\displaystyle{ \dot{\hat{x}}(t)=(a-hc)\hat{x}(t)+hy(t)+bu(t)\quad（a-hc<0） }$

を考える。このときオブザーバベースト・コントローラ

$\displaystyle{ \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$ 行列の固有値を求めよ。
解答オブザーバベースト・コントローラによる閉ループ系は

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

すなわち

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

で表される。この行列 $A_F$ の固有値を求めると

$\displaystyle{ \begin{array}{lll} &&{\rm det}(\lambda I_2-A_F) ={\rm det} \left[\begin{array}{cc} \lambda-a & bf \\ -hc & \lambda-a+hc+bf \end{array}\right]\\ &&=(\lambda-a)(\lambda-a+hc+bf)+bfhc\\ &&=\lambda^2-(a-bf+a-hc)\lambda+(a-bf)(a-hc)\\ &&=(\lambda-(a-bf))(\lambda-(a-hc))=0 \end{array} }$

より， $a-bf$ と $a-hc$ となる。

演習4.2　例題4.3において， $a=-1,b=1,c=,f=4,x(0)=1$ のとき，状態フィードバックによる閉ループ系とオブザーバベースト・コントローラによる閉ループ系の応答を比較したい。 $h=2$ と $h=9$ の場合について，それらの応答をMATLABでシミュレーションせよ。

4.2 可観測性と可検出性

どのような $n$ 次系に対しても状態オブザーバが求まるわけではない。状態オブザーバが構成可能な条件を可検出性という。また，可検出性の十分条件である可観測性の条件も知られている。これらの定義と等価な条件をまとめてお（テキスト「線形システム制御入門」の定理4.1，定理4.2参照）。

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

定義D0:　状態オブザーバを構成可能
条件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$ のすべての不安定固有値)

これらの条件の一つが成り立つとき $n$ 次系は可検出， $(A,B)$ は可検出対という。

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

定義O0:　任意有限時間の入力と出力から，初期状態を一意に決定可能
条件O1: $\displaystyle{\int_0^t \exp(A^T\tau)C^TC\exp(A\tau)\,d\tau>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)$ は可制御対

これらの条件の一つが成り立つとき $n$ 次系は可観測， $(A,B)$ は可観測対という。

例題4.4　つぎの $A$ 行列と $C$ 行列をもつ２次系 $\dot{x}(t)=Ax(t)+Bu(t),y(t)=Cx(t)$ の可観測性を，可観測性行列の階数を求めて判定せよ。

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

$\displaystyle{ (2)\quad A= \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right] ,\ C= \left[\begin{array}{cc} 0 & 1 \end{array}\right] }$

$\displaystyle{ (3)\quad A= \left[\begin{array}{cc} 0 & 1 \\ 0 & -1 \end{array}\right] ,\ C= \left[\begin{array}{cc} 1 & 1 \end{array}\right] }$
解答
(1) 可観測性行列は， $\left[\begin{array}{c} C \\ CA \end{array}\right]= \left[\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right]$ 。この階数は２で，システムの次数２と等しい。したがって，この２次系は可観測である。

(2) 可観測性行列は， $\left[\begin{array}{c} C \\ CA \end{array}\right]= \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right]$ 。この階数は１で，システムの次数２と等しくない。したがって，この２次系は可観測でない。

(3) 可観測性行列は， $\left[\begin{array}{c} C \\ CA \end{array}\right]= \left[\begin{array}{cc} 1 & 1 \\ 0 & 0 \end{array}\right]$ 。この階数は１で，システムの次数２と等しくない。したがって，この２次系は可観測でない。

演習4.4　例題4.4における可観測性の判定を， $A$ 行列の固有値に基づいて行なえ。

4.3　状態オブザーバの低次元化

これまで，状態オブザーバの出力 $\hat{x}$ に状態フィードバックのゲイン行列 $F$ をかけてオブザーバベースト・コントローラを構成した。しかしながら，最終的に必要なのは，状態 $x$ の推定値ではなく，その線形関数 $u=Kx=-Fx$ の推定値であることから，線形関数オブザーバ

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

が考案されている。ここで，サイズ $q\times n$ の適当な行列 $U$ を用いて

$\displaystyle{(4.13)\quad \left\{\begin{array}{l} \hat{A}:\ stable\ matrix\\ UA=\hat{A}U+\hat{B}C\\ UB=\hat{J}\\ K=\hat{C}U+\hat{D}C \end{array}\right. }$

を満足させることができれば

$\displaystyle{(4.14)\quad \underbrace{\dot{\hat{x}}(t)-U\dot{x}(t)}_{\dot{e}(t)}=\hat{A}\underbrace{(\hat{x}(t)-Ux(t))}_{e(t)} }$

が成り立ち，これから

$\displaystyle{(4.15)\quad \hat{x}(t)\rightarrow U\hat{x}(t)\quad(t\rightarrow\infty) }$

したがって

$\displaystyle{(4.16)\quad z(t)\rightarrow u(t)=K\hat{x}(t)\quad(t\rightarrow\infty) }$

を得る（テキスト「線形システム制御入門」の 4.3節参照）。特に， $K=I_n$ の場合は恒等関数オブザーバと呼ばれる。線形関数オブザーバを使用する利点は，状態オブザーバ(4.3)の次元数は制御対象の次元数 $n$ と同じであったが，これを減らすことができる点である（たとえば，倒立振子の次元数は４，出力数は２なので，状態オブザーバを用いた場合は４次元であるが，恒等関数オブザーバを用いれば２次元に，線形関数オブザーバを用いれば１次元に減らすことができる。）。

例題4.5　１入力 $p$ 出力 $2p$ 次系

$\displaystyle{ \left\{\begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot{y}(t) \\ \ddot{y}(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{cc} 0 & I_p \\ A_{21} & A_{22} \end{array}\right] }_A \underbrace{ \left[\begin{array}{c} y(t) \\ \dot{y}(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{cc} 0 \\ B_2 \end{array}\right] }_B u(t) \\ y(t)= \underbrace{ \left[\begin{array}{cc} I_p & 0 \end{array}\right] }_C \underbrace{ \left[\begin{array}{c} y(t) \\ \dot{y}(t) \end{array}\right] }_{x(t)} \end{array}\right. }$

に対して， $p$ 次元の恒等関数オブザーバの一つは次式で与えられることを示せ。

$\displaystyle{ \left\{\begin{array}{l} \dot{\hat{x}}(t)= \underbrace{(A_{22}-L)}_{\hat{A}}\hat{x}(t) +\underbrace{(A_{21}+(A_{22}-L)L)}_{\hat{B}}y(t) +\underbrace{B_{2}}_{\hat{J}}u(t)\\ z(t)= \underbrace{ \left[\begin{array}{c} 0 \\ I_p \end{array}\right] }_{\hat{C}} \hat{x}(t) + \underbrace{ \left[\begin{array}{c} I_p \\ L \end{array}\right] }_{\hat{D}} y(t) \end{array}\right. }$

ここで， $L$ は $\hat{A}$ を安定行列とするサイズ $p\times p$ の適当な行列である。
解答　恒等関数オブザーバの場合は $K=I_{2p}$ であり，条件(??)の第２式と第４式が満足されることを示せばよい。実際，

$\displaystyle{ U= \left[\begin{array}{cc} -L & I_p \end{array}\right] }$

と選べば，次式が成り立つ。

$\displaystyle{ \underbrace{ \left[\begin{array}{cc} A_{21} & A_{22}-L \\ I_p & 0\\ 0 & I_p \end{array}\right] }_{ \left[\begin{array}{cc} UA \\ I_{2p} \end{array}\right]} = \underbrace{ \left[\begin{array}{cc} A_{22}-L & A_{21}+(A_{22}-L)L\\ 0 & I_p\\ I_p & L \end{array}\right] }_{ \left[\begin{array}{cc} \hat{A} & \hat{B} \\ \hat{C} & \hat{D} \end{array}\right]} \underbrace{ \left[\begin{array}{cc} -L & I_p \\ I_p & 0 \end{array}\right] }_{ \left[\begin{array}{cc} U \\ C \end{array}\right]} }$

演習4.5 例題4.1の２次系に対する恒等関数オブザーバを，その固有値が $\lambda=-2$ となるように求めよ。この恒等オブザーバに対して例題4.2と同様のシミュレーションを行え。

例題4.6　例題4.5の１入力 $p$ 出力 $2p$ 次系に対して，線形関数

$\displaystyle{ u(t)= \underbrace{ \left[\begin{array}{cc} K_1 & K_2 \end{array}\right] }_K \underbrace{ \left[\begin{array}{c} y(t) \\ \dot{y}(t) \end{array}\right] }_{x(t)} }$

を推定する線形関数オブザーバの一つは，つぎに $(p+1)$ 入力 $1$ 出力{\bf 1次系}として与えられることを示せ。

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

ただし， $L=A_{22}-\lambda I_p$ と定める。
解答　条件(4.13)の第２式と第４式が満足されることを示せばよい。実際，

$\displaystyle{ U=K_2 \left[\begin{array}{cc} -L & I_p \end{array}\right] }$

と選べば，次式が成り立つ。

$\displaystyle{ \underbrace{ \left[\begin{array}{cc} K_2A_{21} & K_2(A_{22}-L) \\ K_1 & K_2 \end{array}\right] }_{ \left[\begin{array}{cc} UA \\ K \end{array}\right]} = \underbrace{ \left[\begin{array}{cc} \lambda & K_2(A_{21}+\lambda L) \\ 1 & K_1+K_2L \end{array}\right] }_{ \left[\begin{array}{cc} \hat{A} & \hat{B} \\ \hat{C} & \hat{D} \end{array}\right]} \underbrace{ \left[\begin{array}{cc} -K_2L & K_2 \\ I_p & 0 \end{array}\right] }_{ \left[\begin{array}{cc} U \\ C \end{array}\right]} }$

演習4.6　例題4.1の２次系に対して，状態フィードバック $u(t)=-y(t)-2\dot{y}(t)$ を考える。このとき，この状態フィードバックを実施する１次元線形関数オブザーバを，その固有値が $\lambda=-2$ となるように求めよ。また，その出力が $u(t)$ を追跡するシミュレーションを行え。

演習問題の解答

【演習4.1】　行列 $A$ の特性多項式は

${\rm det}(\lambda I_2-A)= {\rm det}\left[\begin{array}{cc} \lambda & -1 \\ 0 & \lambda \end{array}\right] =\lambda^2+\underbrace{0}_{a_1}\lambda+\underbrace{0}_{a_2}$

行列 $A-HC$ の特性多項式は

${\rm det}(\lambda I_2-A+HC)=(\lambda-(-2))^2= \lambda^2+\underbrace{4}_{a'_1}\lambda+\underbrace{4}_{a'_2}$ \\
これらから，オブザーバゲイン $H$ は，つぎのように計算される。\\
\hspace{5mm}
$H^T= \left[\begin{array}{cc} 4-0 & 4-0 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right]^{-1} \left[\begin{array}{cc} 1 & 0 \\ 1 & 1 \end{array}\right]^{-1} =\left[\begin{array}{cc} 0 & 4 \end{array}\right]$

したがって，求める状態オブザーバ $\dot{\hat{x}}=(A-HC)\hat{x}+Hx+Bu$ は

$\left[\begin{array}{c} \dot{\hat{x}}_1 \\ \dot{\hat{x}}_2 \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ -4 & -4 \end{array}\right] \left[\begin{array}{c} \hat{x}_1 \\ \hat{x}_2 \end{array}\right] + \left[\begin{array}{c} 0 \\ 4 \end{array}\right] y + \left[\begin{array}{c} 0 \\ 1 \end{array}\right] u$

【演習4.2】

%obs_err2.m
A=[0 1;0 0]; C=[1 1]; H=[0;4];
sys1=ss(A,zeros(2,1),eye(2,2),0); x0=[1;0.5];
sys2=ss(A-H*C,zeros(2,1),eye(2,2),0); xh0=[0;0];
t=0:0.1:5; x=initial(sys1,x0,t); e=initial(sys2,xh0-x0,t); xh=x+e;
figure(2),subplot(121),plot(t,x(:,1),t,xh(:,1)),axis([0 5 0 4]),grid
figure(2),subplot(122),plot(t,x(:,2),t,xh(:,2)),axis([0 5 0 2]),grid

【演習4.3】

%observer_based_controller.m
a=-1; b=1; c=1; f=4; sys=ss(a-b*f,b,c,0);
h1=9; A1=[a-b*f -b*f;0 a-h1*c]; sys1=ss(A1,zeros(2,1),[1 0],0);
h2=2; A2=[a-b*f -b*f;0 a-h2*c]; sys2=ss(A2,zeros(2,1),[1 0],0);
t=0:0.1:3; x0=1; xh0=0; X0=[x0; xh0-x0];
y=initial(sys,x0,t); y1=initial(sys1,X0,t); y2=initial(sys2,X0,t);
figure(3),plot(t,y,t,y1,t,y2),grid
legend(‘h=0′,’h=9′,’h=2’)

【演習4.4】
(1) $A$ 行列の固有値は $\lambda_1=\lambda_2=0$ 。

${\rm rank}\, \left[\begin{array}{c} C \\ A-\lambda_i I_2 \end{array}\right]= {\rm rank} \left[\begin{array}{ccc} 1 & 0 \\ 0 & 1 \\ 0 & 0 \end{array}\right]= 2\ \ (i=1,2)$

したがって，この２次系は可観測である。

(2) $A$ 行列の固有値は $\lambda_1=\lambda_2=0$ 。

${\rm rank}\, \left[\begin{array}{c} C \\ A-\lambda_i I_2 \end{array}\right]= {\rm rank} \left[\begin{array}{ccc} 0 & 1 \\ 0 & 1 \\ 0 & 0 \end{array}\right]= 1\ \ (i=1,2)$

したがって，この２次系は可観測ではない。

(3) $A$ 行列の固有値は $\lambda_1=0, \lambda_2=-1$ 。

${\rm rank}\, \left[\begin{array}{c} C \\ A-\lambda_1 I_2 \end{array}\right]= {\rm rank} \left[\begin{array}{ccc} 1 & 1 \\ 0 & 1 \\ 0 & -1 \end{array}\right]= 2$

${\rm rank}\, \left[\begin{array}{c} C \\ A-\lambda_2 I_2 \end{array}\right]= {\rm rank} \left[\begin{array}{ccc} 1 & 1 \\ -1 & 1 \\ 0 & 0 \end{array}\right]= 1$

したがって，この２次系は可観測ではない。

【演習4.5】　 $A_{21}=A_{22}=0$ ， $B_2=1$ だから，例題4.1の恒等オブザーバは，次式となる。

$\left\{\begin{array}{l} \dot{\hat{x}}(t)= \underbrace{-L}_{A_{22}-L} \hat{x}(t) \underbrace{-L^2}_{A_{21}+(A_{22}-L)L} y(t)+ \underbrace{1}_{B_{2}} u(t)\\ z(t)= \left[\begin{array}{c} 0 \\ 1 \end{array}\right] \hat{x}(t) + \left[\begin{array}{c} 1 \\ L \end{array}\right] y(t) \end{array}\right.$

ただし， $L=2$ 。また，恒等オブザーバの出力 $\hat{x}$ が状態 $x(t)$ を追跡していく様子は，つぎのMファイルを用いてシミュレーションできる。

%obs_err3.m
A=[0 1;0 0]; B=[0;1]; C=[1 0];
A21=0; A22=0; B2=1; L=2; U=[-L 1];
AH=A22-L; BH=A21-(A22-L)*L; CH=[0;1]; DH=[1;L]; JH=B2;
sys1=ss(A,zeros(2,1),eye(2,2),0); x0=[1;0.5];
sys2=ss(AH,0,1,0); xh0=0;
t=0:0.1:5;
x=initial(sys1,x0,t)’; y=C*x;
e=initial(sys2,xh0-U*x0,t)’; xh=U*x+e; z=CH*xh+DH*y;
figure(4),subplot(121),plot(t,x(1,:),t,z(1,:)),axis([0 5 0 4]),grid
figure(4),subplot(122),plot(t,x(2,:),t,z(2,:)),axis([0 5 0 2]),grid

【演習4.6】　 $A_{21}=A_{22}=0$ ， $B_2=1$ ， $K_1=-1$ ， $K_2=-2$ だから，例題4.1の関数オブザーバは，次式となる。

$\left\{\begin{array}{l} \dot{\hat{x}}(t)= \lambda \hat{x}(t) +\underbrace{2\lambda^2}_{K_2(A_{21}+\lambda L)} y(t) \underbrace{-2}_{K_2B_2} u(t)\\ z(t)= \hat{x}(t) \underbrace{-(1-2\lambda)}_{K_1+K_2L} u(t) \end{array}\right.$

ただし， $\lambda=-2$ ， $L=-\lambda$ 。また，関数オブザーバの出力 $z(t)$ が状態フィードバック $u(t)$ を追跡していく様子は，つぎのMファイルを用いてシミュレーションできる。

%obs_err4.m
A=[0 1;0 0]; B=[0;1]; C=[1 0]; F=[1 2];
A21=0; A22=0; B2=1; K1=-1; K2=-2;
lambda=-1; L=A22-lambda; U=K2*[-L 1];
AH=lambda; BH=K2*(A21+lambda*L); JH=K2*B2; CH=1; DH=K1+K2*L;
AA=[A-B*F B*CH;0 0 AH]; BB=[B; 0]; CC=[-F 0];
sys1=ss(A-B*F,zeros(2,1),-F,0); x0=[1;0.5];
sys2=ss(AA,BB,CC,0); xh0=0; X0=[x0;xh0-U*x0];
t=0:0.1:5; u1=initial(sys1,x0,t); u2=initial(sys2,X0,t);
figure(5),plot(t,u1,t,u2),axis([0 5 -2 2]),grid
legend(‘sf’,’obc’)

3 状態フィードバックと可制御性

投稿日時: 2022年8月18日投稿者: cacsd

【本章のねらい】
・状態フィードバックを設計する。
・可制御性と可安定性を判定する。

3.1 状態フィードバック

いま制御対象は平衡状態にあるとし，何らかの要因でこれが乱されたとき，適当な手段を用いて，速やかに元の平衡状態に戻したい。そのような手段の一つとして， $m$ 入力 $p$ 出力 $n$ 次元線形系（ $n$ 次系）

$\displaystyle{(3.1)\quad \left\{\begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t) \\ y(t)=Cx(t) \end{array}\right. }$

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

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

を考える。このとき，(3.2)式を(3.1)式に代入して，閉ループ系

$\displaystyle{(3.3)\quad \left\{\begin{array}{l} \dot{x}(t)=\underbrace{(A-BF)}_{A_F}x(t) \\ y(t)=Cx(t) \end{array}\right. }$

を得る。このブロック線図を次に示す。

上の制御目的が達成されるためには，閉ループ系の $A$ 行列，すなわち，行列 $A_F=A-BF$ が安定行列となるように，状態フィードバックのゲイン行列 $F$ を決めればよい。実際，

(3.4)「任意の $x(0)\ne0$ に対して， $x(t)=\exp(A-BF)x(0)\rightarrow 0\ (t\rightarrow\infty)$ 」

が成り立ち，これは平衡状態 $x=0$ に戻ることを意味するからである。

まず，１次系の状態フィードバックの例を考える。

例題3.1　時定数 $T$ と定常ゲイン $K$ をもつ１次系

$\displaystyle{ \dot{x}(t)=-\frac{1}{T}x(t)+\frac{K}{T}u(t) }$

に対して，新しい入力 $v(t)$ をもつ状態フィードバック

$\displaystyle{ u(t)=-\underbrace{\frac{T}{K}(\frac{1}{T'}-\frac{1}{T})}_fx(t)+\underbrace{\frac{T}{K}\frac{K'}{T'}}_gv(t) }$

を行うと閉ループ系の時定数は $T'$ ，定常ゲインは $K'$ となることを示せ。
解答　フィードバック式を状態方程式の $u(t)$ に代入すると，閉ループ系は

$\displaystyle{ \dot{x}(t)=-\frac{1}{T}x(t)+\frac{K}{T}(-\frac{T}{K}(\frac{1}{T'}-\frac{1}{T})x(t)+\frac{T}{K}\frac{K'}{T'}v(t))=-\frac{1}{T'}x(t)+\frac{K'}{T'}v(t) }$

となる。これは時定数 $T'$ と定常ゲイン $K'$ をもつ１次系を表している。

演習3.1　１次系 $\dot{x}(t)=-x(t)+u(t)$ に対するフィードバック $u(t)=-fx(t)+gv(t)$ を，閉ループ系の時定数が $1/10$ ，定常ゲインが $1$ となるように
定めよ。

演習3.2　例題2.5で得た図上( $t$ – $x$ 平面)に，望ましい閉ループ系の時定数 $T'$ ，定常ゲイン $K'$ を表す点 $(T',K')$ を指定し，ginputを使って読み込み，これを達成するフィードバック $u(t)=-fx(t)+gv(t)$ を定めよ。

つぎに，２次系の状態フィードバックの例を考える。

例題3.2　減衰係数 $\zeta$ と固有角周波数 $\omega_n^2$ をもつ２次系

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

に対して，新しい入力 $v(t)$ をもつ状態フィードバック

$\displaystyle{ u(t)=- \underbrace{ \left[\begin{array}{cc} \frac{1}{\omega_n^2}(\omega_n'\,^2-\omega_n^2) & \frac{2}{\omega_n^2}(\zeta'\omega_n'-\zeta\omega_n) \end{array}\right] }_{F} \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] +\underbrace{\frac{\omega_n'\,^2}{\omega_n^2}}_Gv(t) }$

を行うと，閉ループ系の減衰係数は $\zeta'$ と固有角周波数は $\omega_n'\,^2$ となることを示せ。
解答　 $u(t)=-Fx(t)+Gv(t)$ を状態方程式 $\dot{x}(t)=Ax(t)+Bu(t)$ に代入すると，閉ループ系は $\dot{x}(t)=(A-BF)x(t)+BGv(t)$ となる。ここで

$\displaystyle{ \begin{array}{lll} &&A-BF= \left[\begin{array}{cc} 0 & 1 \\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] - \left[\begin{array}{c} 0 \\ \omega_n^2 \end{array}\right]\\ && \times \left[\begin{array}{cc} \frac{1}{\omega_n^2}(\omega_n'\,^2-\omega_n^2) & \frac{2}{\omega_n^2}(\zeta'\omega_n'-\zeta\omega_n) \end{array}\right] =\left[\begin{array}{cc} 0 & 1 \\ -\omega_n'\,^2 & -2\zeta'\omega_n' \end{array}\right] \end{array} }$

$\displaystyle{ BG= \left[\begin{array}{c} 0 \\ \omega_n^2 \end{array}\right] \frac{\omega_n'\,^2}{\omega_n^2} = \left[\begin{array}{c} 0 \\ \omega_n'\,^2 \end{array}\right] }$

これは減衰係数 $\zeta'$ と固有角周波数 $\omega_n'\,^2$ をもつ２次系を表している。

演習3.3　例題2.6で得た図上( $t$ – $y$ 平面)に，望ましい第１番目のオーバーシュートの頂点の座標 $(T_p',y(T_p'))$ を指定し，ginputを使って読み込み，対応する減衰係数 $\zeta'$ と固有角周波数 $\omega_n'\,^2$ を(2.22)式を使って求め，これを達成するフィードバック $u(t)=-f_1x_1(t)-f_2x_2(t)+gv(t)$ を定めよ。

さて， $n$ 次系に対する状態フィードバックの設計法を考える。まず，１入力系の場合を考え，つぎの条件を仮定する（１入力系の可制御性行列は $n$ 次の正方行列となり，本条件はその正則性を意味する。）。

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

また，行列 $A$ の固有値を $\lambda_1,\cdots,\lambda_n$ ，行列 $A_F=A-BF$ の固有値を $\lambda'_1,\cdots, \lambda'_n$ とするとき，，それぞれの特性多項式を次式で表す。

$\displaystyle{(3.6)\quad {\rm det}(\lambda I_n-A)&=&(\lambda-\lambda_1)\cdots(\lambda-\lambda_n) \nonumber\\ &=&\lambda^n+a_1\lambda^{n-1}+\cdots+a_n }$

$\displaystyle{(3.7)\quad {\rm det}(\lambda I_n-A_F)&=&(\lambda-\lambda'_1)\cdots(\lambda-\lambda'_n)\nonumber\\ &=&\lambda'\,^n+a'_1\lambda'\,^{n-1}+\cdots+a'_n }$

このとき，閉ループ系の $A$ 行列 $A_F=A-BF$ の固有値を，指定された安定固有値 $\lambda'_1,\cdots, \lambda'_n$ に設定する状態フィードバックのゲイン行列 $F$ は

$\displaystyle{(3.8)\quad \begin{array}{lll} F&=&\left[\begin{array}{cccc} a_n'-a_n & \cdots & a_2'-a_2 & a_1'-a_1 \end{array}\right] \nonumber\\ &\times& \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_2 & a_1 & \cdots & 0 & 0 \\ a_1 & 1 & \cdots & 0 & 0 \\ 1 & 0 & \cdots & 0 & 0 \end{array}\right]^{-1} \nonumber\\ &\times& \left[\begin{array}{cccc} B & AB &\cdots & A^{n-1}B \end{array}\right]^{-1} \end{array} }$

または

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

で与えられる（テキスト「線形システム制御入門」の3.3節参照）。(3.8})式と(3.9})式を比較すると，前者は $A$ の特性多項式の係数の計算を，後者は $A$ のべき乗計算を必要とすることに注意する。

例題3.3　２次系

$\displaystyle{ \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 & 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) }$

に対する状態フィードバック $u=-Fx$ を，行列 $A-BF$ の固有値が，つぎのものとなるように求めよ。
(1) $\lambda'_1=-1,\ \lambda'_2=-2$
(2) $\lambda'_1=-1+j,\ \lambda'_2=-1-j$
解答　 $A$ 行列の特性多項式は

$\displaystyle{ {\rm det}(\lambda I_2-A)= {\rm det}\left[\begin{array}{cc} \lambda & -1 \\ 0 & \lambda \end{array}\right] =\lambda^2+\underbrace{0}_{a_1}\lambda+\underbrace{0}_{a_2} }$

(1) 行列 $A-BF$ の特性多項式は

$\displaystyle{ {\rm det}(\lambda I_2-A+BF)=(\lambda-(-1))(\lambda-(-2))= \lambda^2+\underbrace{3}_{a'_1}\lambda+\underbrace{2}_{a'_2} }$

したがって，ゲイン行列 $F$ は，つぎのように計算される。

$\displaystyle{ F= \left[\begin{array}{cc} 2-0 & 3-0 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right]^{-1} \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right]^{-1} =\left[\begin{array}{cc} 2 & 3 \end{array}\right] }$

(2) 行列 $A-BF$ の特性多項式は

$\displaystyle{ {\rm det}(\lambda I_2-A+BF)=(\lambda-(-1+j))(\lambda-(-1-j))= \lambda^2+\underbrace{2}_{a'_1}\lambda+\underbrace{2}_{a'_2} }$

したがって，ゲイン行列 $F$ は，つぎのように計算される。

$\displaystyle{ F= \left[\begin{array}{cc} 2-0 & 2-0 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right]^{-1} \left[\begin{array}{cc} 0 & 1 \\ 1 & -1 \end{array}\right]^{-1} =\left[\begin{array}{cc} 2 & 2 \end{array}\right] }$

演習3.4　つぎの２次系に対する状態フィードバック $u=-Fx$ を，閉ループ系の $A$ 行列の固有値が $\lambda'_1=\lambda'_2=-1$ となるように求めよ。
(1) $\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)$

(2) $\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)$

最後に多入力をもつ $n$ 次系に対して状態フィードバックのゲイン行列を求めることを考える。 $A-BF$ の固有値 $\lambda'_1,\cdots,\lambda'_n$ に対応する固有ベクトルを $v_1,\cdots,v_n$ とするとき，次式が成り立つ。

$\displaystyle{(3.10)\quad (A-BF)v_i=\lambda'_iv_i\quad (i=1,\cdots,n) }$

ここで

$\displaystyle{(3.11)\quad g_i=Fv_i\quad (i=1,\cdots,n) }$

とおくと

$\displaystyle{(3.12)\quad (A-\lambda'I_n)v_i=Bg_i\quad (i=1,\cdots,n) }$

これから， $F$ を次式で求めることが考えられる。

$\displaystyle{(3.13)\quad \begin{array}{lll} F&=& \left[\begin{array}{cccc} g_1 & \cdots & g_n \end{array}\right] \nonumber\\ &\times& \left[\begin{array}{cccc} (A-\lambda_1'I_n)^{-1}Bg_1 &\cdots & (A-\lambda_n'I_n)^{-1}Bg_n \end{array}\right]^{-1} \end{array} }$

ただし， $\lambda'_1,\cdots,\lambda'_n$ は $A$ の固有値に等しくないとし， $m$ 次元ベクトル $g_1,\cdots,g_n$ は，上の逆行列が存在する範囲で適切に指定するものとする。これから，多入力系の場合，状態フィードバックは，固有値を指定しただけでは，一意に定まらないことがわかる。

例題3.4　２入力２次系

$\displaystyle{ \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 & 0 \\ 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}{cc} 1 & 1 \\ 1 & -1 \end{array}\right] }_{B} u(t) }$

に対するつぎの状態フィードバックによる閉ループ系おける行列 $A-BF$ の固有値を求めよ。

(1) $u(t)=- \underbrace{ \left[\begin{array}{cc} 1 & 1 \\ 1 & -1 \end{array}\right] }_{F} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \end{array}\right] }_{{x}(t)}$

(2) $u(t)=- \underbrace{ \left[\begin{array}{cc} 1 & 0 \\ 1 & -2 \end{array}\right] }_{F} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \end{array}\right] }_{{x}(t)}$
解答
(1) 閉ループ系における行列 $A-BF$ は

$A-BF=\left[\begin{array}{cc} 0 & 0 \\ 0 & -1 \end{array}\right] - \left[\begin{array}{cc} 1 & 1 \\ 1 & -1 \end{array}\right] \left[\begin{array}{cc} 1 & 1 \\ 1 & -1 \end{array}\right] = \left[\begin{array}{cc} -2 & 0 \\ 0 & -3 \end{array}\right]$

この特性多項式は

${\rm det}( \lambda I_2- \left[\begin{array}{cc} -2 & 0 \\ 0 & -3 \end{array}\right] ) ={\rm det} \left[\begin{array}{cc} \lambda+2 & 0 \\ 0 & \lambda+3 \end{array}\right]=(\lambda+2)(\lambda+3)$

したがって，行列 $A-BF$ の固有値は， $-2,-3$ 。

(2) 閉ループ系の $A$ 行列は

$A-BF=\left[\begin{array}{cc} 0 & 0 \\ 0 & -1 \end{array}\right] - \left[\begin{array}{cc} 1 & 1 \\ 1 & -1 \end{array}\right] \left[\begin{array}{cc} 1 & 0 \\ 1 & -2 \end{array}\right] = \left[\begin{array}{cc} -2 & 2 \\ 0 & -3 \end{array}\right]$

この特性多項式は

${\rm det}( \lambda I_2- \left[\begin{array}{cc} -2 & 2 \\ 0 & -3 \end{array}\right] ) ={\rm det} \left[\begin{array}{cc} \lambda+2 & -2 \\ 0 & \lambda+3 \end{array}\right]=(\lambda+2)(\lambda+3)$

したがって，行列 $A-BF$ の固有値は， $-2,-3$ 。

演習3.5　例題3.4の２つの状態フィードバックは，公式(3.13})において

(1) $\left[\begin{array}{cc} g_1 & g_2 \end{array}\right]= \left[\begin{array}{cc} 1 & 1 \\ 1 & -1 \\ \end{array}\right]$

(2) $\left[\begin{array}{cc} g_1 & g_2 \end{array}\right]= \left[\begin{array}{cc} 1 & 1 \\ 1 & 2 \\ \end{array}\right]$

と指定して得られることを，MATLABを用いて確かめよ。

3.2 可制御性と可安定性

どのような $n$ 次系に対しても，閉ループ系を安定化をする状態フィードバックが求まるわけではない。その条件を可安定性という。また，(3.6})は，可安定性の十分条件である可制御性の条件として知られている。これらの定義と等価な条件をまとめておく（テキスト「線形システム制御入門」の定理3.5，定理3.6参照）。

【可安定性の定義とその等価な条件】
定義S0:　状態フィードバックにより安定化可能
条件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$ のすべての不安定固有値)

これらの条件の一つが成り立つとき $n$ 次系は可制御， $(A,B)$ は可制御対という。

【可制御性の定義とその等価な条件】
定義C0:　任意初期状態を，任意有限時間内に，任意状態に移動可能
条件C1:　 $\displaystyle{\int_0^t \exp(A\tau)BB^T\exp(A^T\tau)\,d\tau>0 \quad (\forall t>0)}$
条件C2:　 ${\rm rank}\, \underbrace{ \left[\begin{array}{cccc} B & AB & \cdots & A^{n-1}B \end{array}\right] }_{controllability\ matrix} =n$
条件C3:　 $F$ を選んで， $A-BF$ の固有値を任意に設定可能
条件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$ のすべての固有値)

これらの条件の一つが成り立つとき $n$ 次系は可制御， $(A,B)$ は可制御対という。

例題3.5　つぎの $A$ 行列と $B$ 行列をもつ２次系 $\dot{x}(t)=Ax(t)+Bu(t)$ の可制御性を，可制御性行列の階数を求めて判定せよ。

(1) $A= \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right] ,\ B= \left[\begin{array}{c} 0 \\ 1 \end{array}\right]$

(2) $A= \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right] ,\ B= \left[\begin{array}{c} 1 \\ 0 \end{array}\right]$

(3) $A= \left[\begin{array}{cc} 0 & 1 \\ 0 & -1 \end{array}\right] ,\ B= \left[\begin{array}{c} 0 \\ 1 \end{array}\right]$

解答
(1) 可制御性行列は， $\left[\begin{array}{cc} B & AB \end{array}\right]= \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right]$ である。この階数は２で，システムの次数２と等しい。したがって，この２次系は可制御である。

(2) 可制御性行列は， $\left[\begin{array}{cc} B & AB \end{array}\right]= \left[\begin{array}{cc} 1 & 0 \\ 0 & 0 \end{array}\right]$ である。この階数は１で，システムの次数２と等しくない。したがって，この２次系は可制御でない。

(3) 可制御性行列は， $\left[\begin{array}{cc} B & AB \end{array}\right]= \left[\begin{array}{cc} 0 & 1 \\ 1 & -1 \end{array}\right]$ である。この階数は２で，システムの次数２と等しい。したがって，この２次系は可制御である。

演習3.6　つぎの $A$ 行列と $B$ 行列をもつ３次系 $\dot{x}(t)=Ax(t)+Bu(t)$ の可制御性を，可制御性行列の階数を求めて判定せよ。

(1) $A= \left[\begin{array}{ccc} 0 & 1 & 0 \\ 0 & -1 & 1 \\ 0 & 0 & -1 \end{array}\right] ,\ B= \left[\begin{array}{c} 0 \\ 1 \\ 0 \end{array}\right]$

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

MATLABを用いて可制御性を判定するには，たとえば例題3.5(3)の $A$ 行列と $B$ 行列に対しては，つぎのコマンドを与えればよい（tolは零判定基準でデータの誤差を考慮して決める。省略すればデフォルト値が用いられる。）。

%controllability_check.m
A=[0 1;0 -1]; B=[0;1]; n=size(A,1); r=eig(A), tol=0.01;
for i=1:n, rank([B A-r(i)*eye(n)],tol)==n, end

ここで，４行目の結果がすべて真であれば，可制御である。

演習3.7　上のコマンドを用いて演習3.6の $A$ 行列と $B$ 行列をもつ３次系 $\dot{x}(t)=Ax(t)+Bu(t)$ の可制御性を判定せよ。

例題3.6　例題3.5の $A$ 行列と $B$ 行列をもつ２次系 $\dot{x}(t)=Ax(t)+Bu(t)$ の可安定性を判定せよ。
解答
(1) $A$ 行列の固有値は $\lambda_1=\lambda_2=0$ 。ともに不安定固有値。

${\rm rank} \left[\begin{array}{cc} B & A-\lambda_i I_2 \end{array}\right]= {\rm rank} \left[\begin{array}{ccc} 0 & 0 & 1 \\ 1 & 0 & 0 \end{array}\right]= 2\ \ (i=1,2)$

したがって，この２次系は可安定である。

(2) $A$ 行列の固有値は $\lambda_1=\lambda_2=0$ 。ともに不安定固有値。

${\rm rank} \left[\begin{array}{cc} B & A-\lambda_i I_2 \end{array}\right]= {\rm rank} \left[\begin{array}{ccc} 1 & 0 & 1 \\ 0 & 0 & 0 \end{array}\right]= 1\ \ (i=1,2)$

したがって，この２次系は可安定ではない。

(3) $A$ 行列の固有値は $\lambda_1=0,\ \lambda_2=-1$ 。 $\lambda_1$ のみ不安定固有値。

${\rm rank} \left[\begin{array}{cc} B & A-\lambda_1 I_2 \end{array}\right]= {\rm rank} \left[\begin{array}{ccc} 0 & 0 & 1 \\ 1 & 0 & -1 \end{array}\right]= 2$

したがって，この２次系は可安定である。

演習3.8　演習3.6の $A$ 行列と $B$ 行列をもつ３次系 $\dot{x}(t)=Ax(t)+Bu(t)$ の可安定性を判定せよ。

MATLABを用いて可安定性を判定するには，たとえば例題3.5(3)の $A$ 行列と $B$ 行列に対しては，つぎのコマンドを与えればよい。

%stabilizability_check.m
A=[0 1;0 -1]; B=[0;1]; n=size(A,1); r=eig(A), tol=0.01;
for i=1:n
if real(r(i))>=0, rank([B A-r(i)*eye(n)],tol)==n, end
end

演習3.9　上のコマンドを用いて演習3.6の $A$ 行列と $B$ 行列をもつ３次系 $\dot{x}(t)=Ax(t)+Bu(t)$ の可安定性を判定せよ。

演習問題の解答

演習3.1　 $T=1,K=1$ ， $T'=1/10,K'=1$ より， $f=\frac{1}{1}(\frac{1}{1/0.1}-\frac{1}{1})=9$ ， $g=\frac{1}{1}\frac{1}{0.1}=10$ 。

演習3.2　たとえば，つぎのMファイルを実行すればよい。

%sf1.m
T1=1; K1=1; a1=-1/T1; b1=K1/T1; sys1=ss(a1,b1,1,0);
t=0:0.1:5; step(sys1,t); x=ginput(1); T2=x(1); K2=x(2);
f=T1/K1*(1/T2-1/T1); g=T1/K1*K2/T2;
a2=a1-b1*f; b2=b1*g; sys2=ss(a2,b2,1,0);
hold on, step(sys2,t)

演習3.3　たとえば，つぎのMファイルを実行すればよい。

%sf2.m
w1=1; z1=0.01; A1=[0 1;-w1^2 -2*z1*w1],; B1=[0;w1^2]; C=[1 0];
sys1=ss(A1,B1,C,0);
t=0:0.05:10; step(sys1,t); x=ginput(1); Tp=x(1); p0=x(2)-1;
w2=sqrt(log(p0)^2+pi^2)/Tp; z2=sqrt(log(p0)^2/(log(p0)^2+pi^2));
Kp=(w2^2-w1^2)/w1^2, Kd=(2/w1^2)*(z2*w2-z1*w1)
F=[Kp Kd]; G=w2^2/w1^2;
A2=A1-B1*F; B2=B1*G; sys2=ss(A2,B2,C,0);
hold on, step(sys2,t)

演習3.4　行列 $A-BF$ の特性多項式は

${\rm det}(\lambda I_2-A+BF)=(\lambda-(-1))^2= \lambda^2+\underbrace{2}_{a'_1}\lambda+\underbrace{1}_{a'_2}$

(1) $A$ 行列の特性多項式は
${\rm det}(\lambda I_2-A)= {\rm det}\left[\begin{array}{cc} \lambda & -1 \\ 0 & \lambda+1 \end{array}\right] =\lambda(\lambda+1) =\lambda^2+\underbrace{1}_{a_1}\lambda+\underbrace{0}_{a_2}$

したがって，ゲイン行列 $F$ は，つぎのように計算される。

$F= \left[\begin{array}{cc} 1-0 & 2-1 \end{array}\right] \left[\begin{array}{cc} 1 & 1 \\ 1 & 0 \end{array}\right]^{-1} \left[\begin{array}{cc} 0 & 1 \\ 1 & -1 \end{array}\right]^{-1} =\left[\begin{array}{cc} 1 & 1 \end{array}\right]$

(2) $A$ 行列の特性多項式は

${\rm det}(\lambda I_2-A)= {\rm det}\left[\begin{array}{cc} \lambda & -1 \\ 1 & \lambda \end{array}\right] =\lambda^2+1 =\lambda^2+\underbrace{0}_{a_1}\lambda+\underbrace{1}_{a_2}$

したがって，ゲイン行列 $F$ は，つぎのように計算される。

$F= \left[\begin{array}{cc} 1-1 & 2-0 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right]^{-1} \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right]^{-1} =\left[\begin{array}{cc} 0 & 2 \end{array}\right]$

演習3.5　たとえば，つぎのMファイルを実行すればよい。

%sf_minputs.m
A=[0 0;0 -1]; B=[1 1;1 -1]; r1=-2; r2=-3;
disp(‘(1)’), X1=[1 1;1 -1];
V1=[((A-r1*eye(2))\B)*X1(:,1) ((A-r2*eye(2))\B)*X1(:,2)];
F1=X1/V1, AF1=A-B*F1, ev1=eig(AF1)
disp(‘(2)’), X2=[1 1;1 2];
V2=[((A-r1*eye(2))\B)*X2(:,1) ((A-r2*eye(2))\B)*X2(:,2)];
F2=X2/V2, AF2=A-B*F2, ev2=eig(AF2)

演習3.6

(1) 可制御性行列は，
$\left[\begin{array}{ccc} B & AB & A^2B \end{array}\right]= \left[\begin{array}{ccc} 0 & 1 & -1 \\ 1 & -1 & 1 \\ 0 & 0 & 0 \end{array}\right]$ 。この階数は２で，システムの次数３より小さい。したがって，この３次系は不可制御である。

(2) 可制御性行列は，
$\left[\begin{array}{ccc} B & AB & A^2B \end{array}\right]= \left[\begin{array}{cccccc} 0 & 0 & 1 &-1 &-1 & 1\\ 1 &-1 &-1 & 1 & 0 & 0\\ 0 & 1 & 0 & 2 & 0 & 4 \end{array}\right]$ 。この階数は３で，システムの次数３と等しい。
したがって，この３次系は可制御である。

演習3.7　Mファイル{\tt controllability\_check.m}のデータ{\tt A,B}の定義を次のように書き換える。

(1) A=[0 1 0;0 -1 1;0 0 -1]; B=[0;1;0];
(2) A=[0 1 0;-1 -1 0;0 0 2]; B=[0 0;1 -1;0 1];

演習3.8
(1) $A$ 行列の固有値は $\lambda_1=0,\lambda_2=\lambda_3=-1$ 。 $\lambda_1$ のみ不安定固有値。

${\rm rank} \left[\begin{array}{cc} B & A-\lambda_1 I_3 \end{array}\right]= {\rm rank} \left[\begin{array}{cccc} 0 & 0 & 1 & 0 \\ 1 & 0 & -1 & 1 \\ 0 & 0 & 0 & -1 \end{array}\right]= 3$

したがって，この３次系は可安定である。

(2) $A$ 行列の固有値は $\lambda_1,\lambda_2=(-1\pm j\sqrt{3})/2,\ \lambda_3=2$ 。 $\lambda_3$ のみ不安定固有値。\\

${\rm rank} \left[\begin{array}{cc} B & A-\lambda_3 I_3 \end{array}\right]= {\rm rank} \left[\begin{array}{ccccc} 0 & 0 & -2 & 1 & 0 \\ 1 & -1 & 0 & -2 & 1 \\ 0 & 1 & 0 & 0 & -2 \end{array}\right]= 3$

したがって，この３次系は可安定である。

演習3.9　Mファイル{\tt stabilizability\_check.m}のデータ{\tt A,B}の定義を，次のように書き換える。

(1) A=[0 1 0;0 -1 1;0 0 -1]; B=[0;1;0];
(2) A=[0 1 0;-1 -1 0;0 0 2]; B=[0 0;1 -1;0 1];

2 漸近安定性

投稿日時: 2022年8月18日投稿者: cacsd

【本章のねらい】
・状態空間表現を用いて漸近安定性の判定を行う。
・状態空間表現を用いて時間応答の計算を行う。

2.1 漸近安定性

いま制御対象は平衡状態（物理的な釣り合いの状態）にあるとし，何らかの要因でこれが乱されたとき，その振舞いはつぎの $m$ 入力 $p$ 出力 $n$ 次元線形系（ $n$ 次系）の状態空間表現によって表されるものとする。

$\displaystyle{(2.1)\quad \left\{\begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t) \\ y(t)=Cx(t) \end{array}\right. }$

このとき，もし元の平衡状態に戻るならば，その平衡状態は漸近安定，または $n$ 次系(??)は漸近安定という。平衡状態は零状態 $x=0$ で表し，これを保持する入力は零入力 $u=0$ となるように状態空間表現を得ておくと， $n$ 次系(??)の漸近安定性は， $u=0$ の場合の状態方程式の解が $0$ に収束するかどうかで決まる。 $n$ 次系(??)の零入力応答，すなわち $n$ 次自由系

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

の時間応答 $x(t)$ について

$\displaystyle{(2.3)\quad \forall x(0)\ne0:\ x(t)\rightarrow 0 \quad (t\rightarrow\infty) }$

となれば，漸近安定である。もし

$\displaystyle{(2.4)\quad \exists x(0)\ne0:\ x(t)\rightarrow \hspace{-3.5mm}/\hspace{2.5mm} 0 \quad (t\rightarrow\infty) }$

ならば，漸近安定ではない，すなわち不安定である。

例題2.1　１次自由系 $\dot{x}(t)=ax(t),\ x(0)\ne0$ の時間応答を調べ，漸近安定となる条件を求めよ。
解答　 $\dot{x}(t)=ax(t)$ の解は， $x(t)=e^{at}x(0)$ と表わされる。
(1) $a<0$ ならば， $t\rightarrow\infty$ のとき， $x(t)\rightarrow0$ 。
(2) $a>0$ ならば， $t\rightarrow\infty$ のとき， $x(t)\rightarrow\infty$ 。
(3) $a=0$ ならば， $x(t)=x(0)$ 。
明らかに，１次系が漸近安定である条件は， $a<0$ である。

演習2.1　つぎの１次系が漸近安定となる定数 $f$ の範囲を求めよ。
(1) $\dot{x}(t)=(1-f)x(t)$ 　　(2) $\dot{x}(t)=(-1-2f)x(t)$

● $n$ 次自由系 $\dot{x}(t)=Ax(t)$ の時間応答は

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

で与えられる（テキスト「線形システム制御入門」の定理2.2参照）。ここで， $\exp(At)$ は

$\displaystyle{(2.6)\quad \exp(At)=I_n+At+\frac{1}{2}(At)^2+\cdots+\frac{1}{k!}(At)^k+\cdots }$

で定義される。たとえば

$\displaystyle{(2.7)\quad \exp \left(\left[\begin{array}{cc} \lambda_1 & 0\\ 0 & \lambda_2 \end{array}\right]t\right) = \left[\begin{array}{cc} e^{\lambda_1t} & 0 \\ 0 & e^{\lambda_2t} \end{array}\right] }$
$\displaystyle{(2.8)\quad \exp \left(\left[\begin{array}{cc} \lambda_R & \lambda_I \\ \lambda_I & \lambda_R \end{array}\right]t\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] }$
$\displaystyle{(2.9)\quad \exp \left(\left[\begin{array}{cc} \lambda & 1 \\ 0 & \lambda \end{array}\right]t\right) =e^{\lambda t} \left[\begin{array}{cc} 1 & t \\ 0 & 1 \end{array}\right] }$

のように計算される（テキスト「線形システム制御入門」の定理2.4参照）。これらを用いて，次の例題を考える。

例題2.2　つぎの２次自由系の時間応答を求め，漸近安定性を判定せよ。

$\displaystyle{(1) \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right]= \left[\begin{array}{cc} -1 & 0 \\ 0 & -2 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right]}$

$\displaystyle{(2) \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right]= \left[\begin{array}{cc} 1 & 1 \\ -1 & 1 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right]}$

$\displaystyle{(3) \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right]= \left[\begin{array}{cc} -1 & 1 \\ 0 & -1 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right]}$

解答
(1) 式(??)を用いて，時間応答は次式となる。

$\displaystyle{\left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] = \exp \left(\left[\begin{array}{cc} -1 & 0 \\ 0 & -2 \end{array}\right]t\right) \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] = \left[\begin{array}{cc} e^{-t} & 0 \\ 0 & e^{-2t} \end{array}\right] \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] = \left[\begin{array}{c} e^{-t}x_1(0)\\ e^{-2t}x_2(0) \end{array}\right]}$

これより，つぎが成り立つ。

$\displaystyle{\forall \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] \ne \left[\begin{array}{c} 0 \\ 0 \end{array}\right]:\ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] \rightarrow \left[\begin{array}{c} 0 \\ 0 \end{array}\right] \quad (t\rightarrow\infty)}$

したがって，漸近安定である。

(2) 式(??)を用いて，時間応答は次式となる。

$\displaystyle{\left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] = \exp \left(\left[\begin{array}{cc} 1 & 1 \\ -1 & 1 \end{array}\right]t\right) \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] = e^{t} \left[\begin{array}{cc} \cos t & \sin t \\ -\sin t & \cos t \end{array}\right] \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] = \left[\begin{array}{c} e^{t}(x_1(0)\cos t +x_2(0)\sin t )\\ e^{t}(-x_1(0)\sin t +x_2(0)\cos t ) \end{array}\right]}$

これより，つぎが成り立つ。

たとえば $\left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] = \left[\begin{array}{c} 0 \\ 1 \end{array}\right]$ に対して， $\left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] \rightarrow \hspace{-3.5mm}/\hspace{2.5mm} \left[\begin{array}{c} 0 \\ 0 \end{array}\right] \quad (t\rightarrow\infty)$

したがって，漸近安定でない。

(3) 式(??)を用いて，時間応答は次式となる。

$\displaystyle{\left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] = \exp \left(\left[\begin{array}{cc} -1 & 1 \\ 0 & -1 \end{array}\right]t\right) \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] = e^{-t} \left[\begin{array}{cc} 1 & t \\ 0 & 1 \end{array}\right] \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] = \left[\begin{array}{c} e^{-t}(x_1(0)+t x_2(0))\\ e^{-t}x_2(0) \end{array}\right]}$

これより，つぎが成り立つ。

したがって，漸近安定である。

演習2.2　つぎの２次自由系の時間応答を求め，漸近安定性を判定せよ。

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

$\displaystyle{(2) \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right]= \left[\begin{array}{cc} -1 & 1 \\ -1 & -1 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right]$

$(3) \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]}$

一般に， $n$ 次系(??)が漸近安定であるための必要十分条件は行列 $A$ のすべての固有値の実部が負であることである（テキスト「線形システム制御入門」の定理2.6参照）。すべての固有値の実部が負である行列を安定行列と呼ぶ。また，実部が負の固有値を安定固有値，実部が非負の固有値を不安定固有値と呼ぶ（例題2.1のように，１次系が零固有値をもつ場合（ $a=0$ ，積分器），応答は発散することはない。一方，演習2.2(3)が示すように，２次系で２つの零固有値をもつ場合（２つの積分器が直列結合），入力側の積分器の初期値が零でないならばこれを積分する出力側の積分器の応答は発散する。このため一般には零固有値は不安定固有値とみなす）。

例題2.3　つぎの行列 $A$ をもつ２次系 $\dot{x}(t)=Ax(t)$ の漸近安定性を，行列 $A$ の固有値を求めて判定せよ。

$\displaystyle{(1)\ A=\left[\begin{array}{cc} 0 & 1 \\ -1 & -2 \end{array}\right]}$
$\displaystyle{(2)\ A=\left[\begin{array}{cc} 0 & 1 \\ 0 & -1 \end{array}\right]}$
$\displaystyle{(3)\ A=\left[\begin{array}{cc} 0 & 1 \\ -1 & 1 \end{array}\right]}$

解答

(1)　行列 $A$ の固有値は，特性方程式

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

の解として， $\lambda_1=-1,\lambda_2=-1$ と求められる。２つの固有値が実数で負であるので漸近安定である。

(2)　行列 $A$ の固有値は，特性方程式

$\displaystyle{\det(\lambda I_2-A)= \det\left[\begin{array}{cc} \lambda & -1 \\ 0 & \lambda+1 \end{array}\right] =\lambda(\lambda+1)=0}$

の解として， $\lambda_1=0,\lambda_2=-1$ と求められる。１つの固有値が零であるので漸近安定でない。

(3)　行列 $A$ の固有値は，特性方程式

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

の解として， $\lambda_1=\frac{-1+ j\sqrt{3}}{2}, \lambda_2=\frac{-1- j\sqrt{3}}{2}$ （ $\lambda_1=\frac{+1+ j\sqrt{3}}{2}, \lambda_2=\frac{+1- j\sqrt{3}}{2}$ ）と求められる。２つの固有値の実部が負（正）であるので漸近安定である（ない）。

演習2.3　例題2.2と演習2.2の２次系の漸近安定性を，行列 $A$ の固有値を求めて判定せよ。

MATLABを用いて漸近安定性を判定するには，たとえば例題2.3(3)に対しては，つぎのコマンドを与えればよい。

%stability_check.m
A=[0 1;1 -1]; n=size(A,1); %行列データの定義と次数の取得
poles=eig(A) %行列の固有値の計算
sum(real(poles)<0)==n %論理式による漸近安定性の判定

SCILABを用いて漸近安定性を判定するには，たとえば例題2.3(3)に対しては，つぎのコマンドを与えればよい。

//stability_check.m
A=[0 1;1 -1]; n=size(A,1); //行列データの定義と次数の取得
poles=spec(A) //行列の固有値の計算
sum(real(poles)<0)==n //論理式による漸近安定性の判定

2.2 時間応答

$n$ 次系(??)の 時間応答は，次式で表わされる。

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

ただし，

$\displaystyle{(2.11)\quad G(t)=C\exp(At)B }$

すなわち，時間応答は零入力応答と零状態応答の和となる（テキスト「線形システム制御入門」の定理2.3参照）。

以下では，簡単のために，１入力１出力の場合を考える。 $G(t)$ は インパルス応答と呼ばれる。また，入力を $u(t)=1$ と与える場合の零状態応答

$\displaystyle{(2.12)\quad y(t)=\int_0^tG(t-\tau)\,d\tau=-\int_t^0G(\tau')\,d\tau'=\int_0^tG(\tau')\,d\tau' }$

は，ステップ応答と呼ばれる。したがって，ステップ応答はインパルス応答を積分して，逆にインパルス応答はステップ応答を微分して得られる。

まず，次式で表される１次系を考える。

$\displaystyle{(2.13)\quad \dot{x}(t)=-\frac{1}{T}x(t)+\frac{K}{T}u(t) }$

例題2.4　１次系(??)のステップ応答を求めよ。

解答インパルス応答は $g(t)=\frac{K}{T}e^{-\frac{1}{T}t}$ だから，ステップ応答は

$\displaystyle{x(t)=\int_0^t\frac{K}{T}e^{-\frac{1}{T}(t-\tau)}\cdot1\,d\tau=\frac{K}{T}e^{-\frac{1}{T}t}\left[Te^{\frac{1}{T}\tau}\right]_0^t =Ke^{-\frac{1}{T}t}(e^{-\frac{1}{T}t}-1)}$
$\displaystyle{=K(1-e^{-\frac{1}{T}t})}$

で表わされる。または，インパルス応答を直接積分して

$\displaystyle{x(t)=\int_0^t\frac{K}{T}e^{-\frac{1}{T}\tau}\,d\tau=\frac{K}{T}\left[-Te^{-\frac{1}{T}\tau}\right]_0^t=K(1-e^{-\frac{1}{T}t})}$

演習2.4　つぎの１次系のステップ応答を求めよ。

(1) $\dot{x}(t)=-x(t)+u(t)$ 　　(2) $\dot{x}(t)=x(t)+u(t)$

１次系(??)のステップ応答

$\displaystyle{(2.14)\quad x(t)=K(1-e^{-\frac{1}{T}t}) }$

において， $t\rightarrow\infty$ のとき $x(t)\rightarrow K$ となるので， $K$ は定常ゲインとよばれる。また， $T$ は時定数とよばれ，次式により特徴づけられる。

$\displaystyle{(2.15)\quad x(T)=K(1-\frac{1}{e})=0.632K }$

$\displaystyle{(2.16)\quad \dot{x}(0)=\frac{K}{T} }$

すなわち，時定数は，応答が定常値 $K$ の63.2%に到達する時刻，または応答の初期時刻における接線が定常値 $K$ に到達する時刻として求められる。

例題2.5　つぎの１次系のステップ応答をMATLABで計算し，図示せよ。

$\dot{x}(t)=-x(t)+u(t)$

解答　MATLABに，つぎのコマンドを与えればよい。

%step_resp1.m
sys=ss(-1,1,1,0); %状態空間表現のデータ構造体の定義
t=0:0.1:5; %ステップ応答を計算する時刻の定義
step(sys,t), grid %ステップ応答の計算と図示

SCILABに，つぎのコマンドを与えればよい。

//step_resp1.m
sys=syslin(‘c’,-1,1,1,0); //状態空間表現のデータ構造体の定義
t=0:0.1:5; //ステップ応答を計算する時刻の定義
plot(t,csim(‘step’,t,sys)), ,mtlb_grid //ステップ応答の計算と図示

演習2.5　例題2.5の図上に，時定数を求めるための補助線を引き，交点の座標をginput(1)を使って読み取れ。

つぎに，次式で表される２次系を考える。

$\displaystyle{(2.17)\quad \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] }_{x} \end{array}\right. }$

ここで， $A$ 行列の固有値はつぎのように計算される。

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

これらに対応して， $A$ 行列の固有値分解は，次式で与えられる。

$\displaystyle{(2.19)\quad 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}} \quad (\zeta>1) }$

$\displaystyle{(2.20)\quad 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}} \quad (\zeta<1) }$

$\displaystyle{(2.21)\quad 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}} \quad (\zeta<1)(\zeta=1) }$

これらに対応して，インパルス応答

$\displaystyle{(2.22)\quad G(t)=C\exp(At)B=CV\exp(\Lambda t)V^{-1}B }$

は，次式で与えられる。

$\displaystyle{(2.23)\quad \left\{\begin{array}{ll} G(t)=\frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}(e^{\lambda_2t}-e^{\lambda_1t}) & (\zeta>1)\\ G(t)=\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\sin\lambda_I t & (\zeta<1) \\ G(t)=\lambda^2te^{\lambda t} & (\zeta=1) \end{array}\right. }$

これらに対応して，ステップ応答は，次式で与えられる。

$\displaystyle{(2.24)\quad \left\{\begin{array}{ll} y(t)=1+\frac{1}{\lambda_2-\lambda_1}(\lambda_1e^{\lambda_2t}-\lambda_2e^{\lambda_1t}) & (\zeta>1)\\ y(t)=1-\frac{\omega_n}{\lambda_I}e^{\lambda_Rt}\sin(\lambda_I t-\tan^{-1}\frac{\lambda_I}{\lambda_R}) & (\zeta<1) \\ y(t)=1+(\lambda t-1)e^{\lambda t} & (\zeta=1) \end{array}\right. }$

特に， $\zeta<1$ の場合のステップ応答は，つぎのように与えられる。

$\displaystyle{(2.25)\quad y(t)=1-{1\over\sqrt{1-\zeta^2}}\exp(-\zeta\omega_nt)\sin(\omega_n\sqrt{1-\zeta^2}\,t+\phi) }$

ただし

$\displaystyle{(2.26)\quad \phi=\tan^{-1}{\sqrt{1-\zeta^2}\over\zeta} }$

このとき，ステップ応答の行き過ぎ時間 $T_p$ と行き過ぎ量 $p_0=y(T_p)-1$ は，次式で表される。

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

したがって，図上で第１番目のオーバーシュートの頂点の座標 $(T_p,y(T_p))$ を求めれば，固有角周波数 $\omega_n$ と減衰係数 $\zeta$ を，つぎのように計算できる。

$\displaystyle{(2.28)\quad (\omega_n,\zeta)=\left(\frac{\sqrt{(\ln p_0)^2+\pi^2}}{T_p},\frac{|\ln p_0|}{\sqrt{(\ln p_0)^2+\pi^2}}\right) }$

例題2.6　つぎの２次系のステップ応答をMATLABで計算し，図示せよ。

$\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right]= \left[\begin{array}{cc} 0 & 1 \\ -1 & -0.02 \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)\\ 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}\right. }$

解答 MATLABに，つぎのコマンドを与えればよい。

%step_resp2.m
A=[0 1;-1 -0.02]; B=[0;1]; C=[1 0];
sys=ss(A,B,C,0); %0によって適合するサイズの零行列をD行列に設定
t=0:0.1:10; step(sys,t), grid

SCILABに，つぎのコマンドを与えればよい。

//step_resp2.m
A=[0 1;-1 -0.02]; B=[0;1]; C=[1 0];
sys=syslin(‘c’,A,B,C);
t=0:0.1:10; plot2d(t,csim(‘step’,t,sys)), mtlb_grid

演習2.6　例題2.6の図上で，第１番目のオーバーシュートの頂点の座標をginput(1)を使って読み込み，固有角周波数と減衰係数を同定せよ。

さて，もう一つ重要な $n$ 次系(??)の時間応答として，正弦波入力を与える場合の零状態応答（周波数応答）がある。

例題2.7　つぎの１次系の時間応答をMATLABで計算し，図示せよ。

$\displaystyle{ \dot{x}(t)=-x(t)+u(t),\ x(0)=1 }$

ただし， $u(t)=\sin10t$ とする。

解答 MATLABに，つぎのコマンドを与えればよい。

%sin_resp.m
sys=ss(-1,1,1,0); x0=1;
t=0:0.05:10; u=sin(10*t); %入力を指定
y=lsim(sys,u,t,x0); %時間応答の計算
plot(t,y), grid %時間応答の図示

SCILABに，つぎのコマンドを与えればよい。。

//sin_resp.m
sys=syslin(‘c’,-1,1,1,0); x0=1;
t=0:0.05:10; u=sin(10*t); //入力を指定
y=csim(u,t,sys,x0); //時間応答の計算
plot2d(t,y), mtlb_grid //時間応答の図示

演習2.7　例題2.7の図上に，零入力応答と零状態応答を重ねて表示せよ。

この例が示すように， $n$ 次系(??)に対して正弦波入力 $u(t)=\sin\omega t$ を与える場合の零状態応答は十分時間が経つと振幅と位相が変化した次式となる。

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

ここで，周波数伝達関数 $G(j\omega)$ は

$\displaystyle{(2.30)\quad G(j\omega)=C(j\omega I_n-A)^{-1}B }$

で与えられる。古典制御で用いられるボード線図は，角周波数毎にゲイン $20\log_{10}|G(j\omega)|$ と位相 $\angle G(j\omega)$ を対にしたものである。

例題2.8　つぎの１次系のボード線図をMATLABで図示せよ。

$\dot{x}(t)=-x(t)+u(t)$

解答　MATLABに，つぎのコマンドを与えればよい。

%bode_diag.m
sys=ss(-1,1,1,0);
w={1e-1,1e1}; %角周波数範囲の指定
bode(sys,w), grid %ボード線図の計算と図示

SCILABに，つぎのコマンドを与えればよい。

//bode_diag.m
sys=syslin(‘c’,-1,1,1,0);
fmin=1e-1/(2*%pi); fmax=1e1/(2*%pi); //周波数範囲の指定
bode(sys,fmin,fmax) //ボード線図の計算と図示

演習2.8　例題2.6の２次系のボード線図をMATLABで図示せよ。

最後に，自由系の時間応答（初期値応答）の計算法を示す。

例題2.9　つぎの１次自由系の時間応答をMATLABで計算し，図示せよ。

$\dot{x}(t)=-x(t),\ x(0)=1$

解答　MATLABに，つぎのコマンドを与えればよい。

%free_resp.m
sys=ss(-1,0,1,0); x0=1;
t=0:0.1:5;
initial(sys,x0,t), grid %初期値応答の計算と図示

SCILABに，つぎのコマンドを与えればよい。

//free_resp.m
sys=syslin(‘c’,-1,0,1,0); x0=1;
t=0:0.1:5; u=zeros(1,length(t)); //零入力を指定
y=csim(u,t,sys,x0); //初期値応答の計算
plot2d(t,y), mtlb_grid //初期値応答の図示

演習2.9　例題(2.6)の２次系のインパルス応答をMATLABで図示せよ。

演習問題の解答

演習2.1

(1) 漸近安定であるための条件は， $1-f<0$ 。したがって， $f>1$ 。
(2) 漸近安定であるための条件は， $-1-2f<0$ 。したがって， $f>-\frac{1}{2}$ 。

演習2.2
(1) (??)式を用いて，時間応答は次式となる。

$\displaystyle{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] = \exp \left(\left[\begin{array}{cc} -1 & 0 \\ 0 & 0 \end{array}\right]t\right) \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] \nonumber\\ = \left[\begin{array}{cc} e^{-t} & 0 \\ 0 & 1 \end{array}\right] \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] = \left[\begin{array}{c} e^{-t}x_1(0)\\ x_2(0) \end{array}\right] }$

これより，つぎが成り立つ。

$\displaystyle{ \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] = \left[\begin{array}{c} 0 \\ 1 \end{array}\right] \Rightarrow \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] \rightarrow \left[\begin{array}{c} 0 \\ 0 \end{array}\right] \quad (t\rightarrow\infty) }$

したがって，漸近安定でない。

(2) (??)式を用いて，時間応答は次式となる。

$\displaystyle{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] = \exp \left(\left[\begin{array}{cc} -1 & 1 \\ -1 & -1 \end{array}\right]t\right) \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] \nonumber\\ = e^{-t} \left[\begin{array}{cc} \cos t & \sin t \\ -\sin t & \cos t \end{array}\right] \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] = \left[\begin{array}{c} e^{-t}(x_1(0)\cos t +x_2(0)\sin t )\\ e^{-t}(-x_1(0)\sin t +x_2(0)\cos t ) \end{array}\right] }$

これより，つぎが成り立つ。

$\displaystyle{ \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] \ne \left[\begin{array}{c} 0 \\ 0 \end{array}\right] \Rightarrow \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] \rightarrow \left[\begin{array}{c} 0 \\ 0 \end{array}\right] \quad (t\rightarrow\infty) }$

したがって，漸近安定である。

(3) (??)式を用いて，時間応答は次式となる。

$\displaystyle{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] = \exp \left(\left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right]t\right) \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] \nonumber\\ = \left[\begin{array}{cc} 1 & t \\ 0 & 1 \end{array}\right] \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] = \left[\begin{array}{c} x_1(0)+t x_2(0)\\ x_2(0) \end{array}\right] }$

これより，つぎが成り立つ。

$\displaystyle{ \left[\begin{array}{c} x_1(0) \\ x_2(0) \end{array}\right] = \left[\begin{array}{c} 0 \\ 1 \end{array}\right] \Rightarrow \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] \rightarrow \hspace{-3.5mm}/\hspace{2.5mm} \left[\begin{array}{c} 0 \\ 0 \end{array}\right] \quad (t\rightarrow\infty) }$

したがって，漸近安定でない。

演習2.3　例題1.2について:
(1) $\det(\lambda I_2-A)=(\lambda+1)(\lambda+2)=0$ より，行列 $A$ の固有値は $-1,-2$ 。２つとも実数で負だから漸近安定。
(2) $\det(\lambda I_2-A)=(\lambda-1)^2+1=\lambda^2-2\lambda+2=0$ より，行列 $A$ の固有値は $1\pm j$ 。実部が正だから漸近安定ではない。
(3) $\det(\lambda I_2-A)=(\lambda+1)^2=0$ より，行列 $A$ の固有値は $-1,-1$ 。２つとも実数で負だから漸近安定。

演習1.2について:
(1) $\det(\lambda I_2-A)=(\lambda+1)\lambda=0$ より，行列 $A$ の固有値は $-1,0$ 。零固有値をもつので漸近安定はでない。
(2) $\det(\lambda I_2-A)=(\lambda+1)^2+1=\lambda^2+2\lambda+2=0$ より，行列 $A$ の固有値は $-1\pm j$ 。実部が負だから漸近安定。
(3) $\det(\lambda I_2-A)=\lambda^2=0$ より，行列 $A$ の固有値は２つとも $0$ 。零固有値をもつので漸近安定でない。

演習2.4

(1) $x(t)=\int_0^te^{-(t-\tau)}\,d\tau=e^{-t}\left[e^{\tau}\right]_0^t=e^{-t}(e^{t}-1)=1-e^{-t}$
(2) $x(t)=\int_0^te^{+(t-\tau)}\,d\tau=e^{t}\left[-e^{-\tau}\right]_0^t=e^{t}(-e^{-t}+1)=-1+e^{t}$

演習2.5　定常ゲインは１．したがって，横線 $1-\frac{1}{e}$ を図示すればよい。

%step_resp1_cont.m
hold on %重ね書きの準備
plot([0 5],(1-exp(-1))*[1 1]) %レベル1-1/e=0.632の表示
w=ginput(1)

//step_resp1_cont.m
plot2d([0 5],(1-exp(-1))*[1 1]) %レベル1-1/e=0.632の表示
w=locate(1)

演習2.6

%step_resp2_cont.m
[Tp,p0]=ginput(1); p0=p0-1;
zeta=sqrt(log(p0)^2/(log(p0)^2+pi^2));
wn=sqrt(log(p0)^2+pi^2)/Tp;

//step_resp2_cont.m
w=locate(1); Tp=w(1); p0=w(2)-1;
zeta=sqrt(log(p0)^2/(log(p0)^2+%pi^2));
wn=sqrt(log(p0)^2+%pi^2)/Tp;

演習2.7

%sin_resp_cont.m
x0=1; u=zeros(1,length(t)); y1=lsim(sys,u,t,x0); %零入力応答の計算
x0=0; u=sin(10*t); y2=lsim(sys,u,t,x0); %零状態応答の計算
hold on, plot(t,y1,’g’,t,y2,’r’)

//sin_resp_cont.m
x0=1; u=zeros(1,length(t)); y1=csim(u,t,sys,x0); %零入力応答の計算
x0=0; u=sin(10*t); y2=csim(u,t,sys,x0); %零状態応答の計算
hold on, plot2d(t,y1,’g’,t,y2,’r’)

演習2.8

%bode_diag2.m
A=[0 1;-1 -0.02]; B=[0;1]; C=[1 0]; sys=ss(A,B,C,0);
w={1e-1,1e1}; bode(sys,w), grid

//bode_diag2.m
A=[0 1;-1 -0.02]; B=[0;1]; C=[1 0]; sys=syslin(‘c’,A,B,C,0);
fmin=1e-1/(2*%pi); fmax=1e1/(2*%pi);
bode(sys,fmin,fmax)

演習2.9

%impulse_resp.m
A=[0 1;-1 -0.02]; B=[0;1]; C=[1 0]; sys=ss(A,B,C,0);
t=0:0.1:5; initial(sys,B,t), grid

//impulse_resp.m
A=[0 1;-1 -0.02]; B=[0;1]; C=[1 0]; sys=syslin(‘c’,A,B,C,0);
t=0:0.1:5; u=zeros(1,length(t));
y=csim(u,t,sys,B);

1 状態空間表現

投稿日時: 2022年8月18日投稿者: cacsd

【本章のねらい】
・状態空間表現を導き，その構造をブロック線図で表す。
・状態空間表現に対する操作として座標変換と直列結合を行う。

1.1 状態空間表現の導出とブロック線図

古典制御では，主に１入力１出力をもつ線形ダイナミカルシステム（線形系）の入出力特性に注目し，これを伝達関数で表した。一方，現代制御では，多入出力をもつ線形系において，「まず入力が状態に影響を及ぼし，つぎに状態の一部が出力として現れる」と考え，前者を状態方程式で，後者を出力方程式で表す。すなわち，状態方程式と出力方程式をペアにした状態空間表現

$\displaystyle{(1.1)\quad \left\{\begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t) \\ y(t)=Cx(t)+Du(t) \end{array}\right. }$

を求めることが出発点となる。ここで，実ベクトル $x(t),\ u(t),\ y(t)$ は時刻 $t$ におけるそれぞれ $n$ 個の状態変数， $m$ 個の入力変数， $p$ 個の出力変数を要素にもち，状態，入力，出力を表す。 $n$ を次数または次元数と呼ぶ。また， $A,\ B,\ C,\ D$ は，それぞれサイズ $n\times n,\ n\times m,\ p\times n,\ p\times m$ の実行列である（実は，状態空間表現の係数行列には定着した呼び名がない。以下では，行列 $A,B,C,D$ を，それぞれ $A$ 行列， $B$ 行列， $C$ 行列， $D$ 行列と呼ぶことがある。）。

以下では，(1.1)式の状態空間表現で表される $m$ 個の入力変数， $p$ 個の出力変数， $n$ 個の状態変数をもつ線形系を，簡単に $m$ 入力 $p$ 出力n次系と呼ぶ。これには，つぎのような特別な場合が含まれる。

ここで，(1)は直達項 $Du(t)$ がない場合である。(2)と(4)は入力を考えない自由系(unforced system)の場合で，n次自由系と呼ぶ。(3)と(4)は状態方程式だけを考える場合で， $y=x$ とみなす。

ここでは，制御対象のダイナミクスは連立された線形常微分方程式で表されるとし，これから連立１階線形微分方程式の行列表現を導いて状態方程式を得る。
一方，状態変数の一部がどのように取り出されるかを出力方程式で表す。

例題1.1　つぎの運動方程式を考える。

$\displaystyle{ M\ddot{r}(t)+D\dot{r}(t)+Kr(t)=f(t) }$

ここで， $r(t)$ は変位， $f(t)$ は外力， $M$ は質量， $D$ は減衰係数， $K$ はバネ定数である。 $f(t)$ を入力変数とし， $r(t)$ と $v(t)=\dot{r}(t)$ を状態変数とする２次系の状態方程式を導出せよ。
解答　運動方程式から，つぎの連立１階微分方程式を得る。

$\displaystyle{ \left\{\begin{array}{l} \dot{r}(t)=v(t) \\ \dot{v}(t)=-\frac{D}{M}v(t)-\frac{K}{M}r(t)+\frac{1}{M}f(t) \end{array}\right. }$

これを行列表現して， $r(t)$ と $v(t)$ を状態変数， $f(t)$ を入力変数とする，つぎの状態方程式を得る。

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{r}(t) \\ \dot{v}(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -\frac{K}{M} & -\frac{D}{M} \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} r(t) \\ v(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ \frac{1}{M} \end{array}\right] }_{B} \underbrace{f(t)}_{u(t)} }$

演習1.1　つぎの運動方程式から２次自由系の状態方程式を導出せよ。

$\displaystyle{ \ddot{\theta}(t)=-\frac{3g}{4\ell}\theta(t) }$

演習1.2　つぎの連成した運動方程式から， $u(t)$ を入力変数とする４次系の状態方程式を導出せよ。

$\displaystyle{ \left\{\begin{array}{l} \ddot{x}_1(t)=-(x_1(t)-x_2(t))+u(t) \\ \ddot{x}_2(t)=-(x_2(t)-x_1(t)) \end{array}\right. }$

例題1.2　つぎの運動方程式を考える。

$\displaystyle{ J\ddot{\theta}(t)+D\dot{\theta}(t)=\tau(t) }$

ここで， $\theta(t)$ は回転角， $\tau(t)$ はトルク， $J$ は慣性モーメント， $D$ は減衰係数である。 $\tau(t)$ を入力変数とし， $\theta(t)$ と $\omega(t)=\dot{\theta}(t)$ を状態変数とする２次系の状態方程式を導出せよ。また，次の場合の出力方程式を示せ。
(1) $\theta(t)$ を計測できる場合
(2) $\omega(t)$ を計測できる場合
解答　運動方程式から，つぎの連立１階微分方程式を得る。

$\displaystyle{ \left\{\begin{array}{l} \dot{\theta}(t)=\omega(t) \\ \dot{\omega}(t)=-\frac{D}{J}\omega(t)+\frac{1}{J}\tau(t) \end{array}\right. }$

これを行列表現すると， $\theta(t)$ と $\omega(t)$ を状態変数， $\tau(t)$ を入力変数とする，つぎの状態方程式を得る。

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

(1) $\theta(t)$ を計測できる場合の出力方程式として，次式を得る。

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

(2) $\dot{\theta}(t)$ を計測できる場合の出力方程式として，次式を得る。

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

演習1.3　例題1.2の(2)を考える。角速度を $\omega(t)=\dot{\theta}(t)$ とおくと，運動方程式は $\dot{\omega}(t)+{\omega}(t)=\tau(t)$ となる。いま一定のトルク $\tau^*$ のもとで，一定の角速度 $\omega^*$ を得ているとする。このとき， $x(t)=\omega(t)-\omega^*$ を状態変数， $u(t)=\tau(t)-\tau^*$ を入力変数， $y(t)=\omega(t)-\omega^*$ を出力変数とする状態空間表現を導出せよ（一般に，零状態 $x=0$ が平衡状態を表すように状態空間表現を導出する。）。

演習1.4　つぎの２階線形常微分方程式から， $u(t)$ を入力変数， $y(t)$ を出力変数とする状態空間表現を導出せよ。

$\displaystyle{ \ddot{y}(t)+a_1\dot{y}(t)+a_2y(t)=u(t) }$

MATLABを用いて状態空間表現 sys を定義するには，たとえば例題1.2の(1)で， $J=1,D=-0.01$ の場合，つぎのコマンドを与えればよい。

%state_space.m
A=[0 1;0 -0.01]; B=[0;1]; C=[1 0]; D=0;
sys=ss(A,B,C,D)

ここで， $A$ 行列， $B$ 行列， $C$ 行列， $D$ 行列を参照するには，それぞれ

sys.a，sys.b，sys.c，sys.d

を用いる。たとえば $A$ 行列の $(2,2)$ 要素を-0.1に変更するには，つぎのコマンドを与えればよい。

sys.a(2,2)=-0.1

また， $B$ 行列， $D$ 行列が存在しない自由系 sys0 を定義するには，つぎのコマンドを与えればよい。

%state_space_cont.m
A=[0 1;0 -0.01]; B=[]; C=[1 0]; D=[];
sys0=ss(A,B,C,D)

もし状態方程式だけを指定する場合は，C=[]; D=[]; または C=eye(size(A)); D=0; とする。

SCILABを用いて状態空間表現 sys を定義するには，たとえば例題1.2の(1)で， $J=1,D=-0.01$ の場合，つぎのコマンドを与えればよい。

//state_space.m
A=[0 1;0 -0.01]; B=[0;1]; C=[1 0]; D=0;
sys=syslin(A,B,C,D)

ここで， $A$ 行列， $B$ 行列， $C$ 行列， $D$ 行列を参照するには，それぞれ

sys.a，sys.b，sys.c，sys.d

を用いる。たとえば $A$ 行列の $(2,2)$ 要素を-0.1に変更するには，つぎのコマンドを与えればよい。

A(2,2)=-0.1; sys.a=A

また， $B$ 行列， $D$ 行列が存在しない自由系 sys0 を定義するには，つぎのコマンドを与えればよい。

//state_space_cont.m
A=[0 1;0 -0.01]; B=[]; C=[1 0]; D=[];
sys0=syslin(A,B,C,D)

もし状態方程式だけを指定する場合は，C=[]; D=[]; または C=eye(size(A)); D=zeros(size(A,1),size(B,2)); とする。

さて，状態空間表現のブロック線図を作成することにより，どのように入力変数から状態変数を経て出力変数につながるかを表すことができる。これは，Simulinkなどを用いた時間応答シミュレーションに役に立つ。

例題1.3　つぎの１次系の状態空間表現をブロック線図で表せ。

$\displaystyle{ \left\{\begin{array}{l} \dot{x}(t)=-\frac{1}{T}x(t)+\frac{K}{T}u(t) \\ y(t)=cx(t) \end{array}\right. }$

解答　積分器の出力を $x$ とすると，その入力は $\dot{x}$ であることに注意し，状態方程式の左辺 $\dot{x}$ と右辺 $-\frac{1}{T}x+\frac{K}{T}u$ を等しくなるよう描いて，つぎのブロック線図を得る（加え合せ点○で符号を表示しない場合は＋とみなす。また，●は引き出し点を表す。）。

演習1.5　つぎの１次系の状態空間表現をブロック線図で表せ。

$\displaystyle{ \left\{\begin{array}{l} \dot{x}(t)=u(t) \\ y(t)=x(t)+u(t) \end{array}\right. }$

例題1.4　つぎの２次系の状態空間表現をブロック線図で表せ。

$\displaystyle{ \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} c_1 & 0 \end{array}\right] %}_{C} %\underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] %}_{x} \end{array}\right. }$

解答　 $\dot{x}_1(t)=x_2(t),\ \dot{x}_2(t)=-\omega_n^2x_1(t)-2\zeta\omega_nx_2(t)+\omega_n^2u(t),\ y(t)=x_1(t)$ より，つぎのブロック線図を得る。

演習1.6　つぎの２次系の状態空間表現をブロック線図で表せ。

$\displaystyle{ \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 \\ 0 & -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} 0 & c_2 \end{array}\right] %}_{C} %\underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] %}_{x} \end{array}\right. }$

演習1.7　演習1.2 で得られた状態空間表現について，ブロック線図を描け。

1.2 状態空間表現の座標変換と直列結合

$n$ 次系の状態空間表現

$\displaystyle{(1.2)\quad \left\{\begin{array}{l} \dot{x}(t)=Ax(t)+Bu(t) \\ y(t)=Cx(t) \end{array}\right. }$

に対して，座標変換

$\displaystyle{(1.3)\quad x'(t)=Tx(t)\quad({\rm det}T\ne0) }$

を行うと，つぎの状態空間表現を得る。

$\displaystyle{(1.4)\quad \left\{\begin{array}{l} \dot{x}'(t)=A'x'(t)+B'u(t) \\ y(t)=C'x'(t) \end{array}\right. }$

ただし，変換後の係数行列は次式で与えられる。

$\displaystyle{(1.5)\quad \left\{\begin{array}{l} A'=TAT^{-1} \\ B'=TB \\ C'=CT^{-1} \\ \end{array}\right. }$

例題1.5　２次系の状態空間表現

$\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right] = \left[\begin{array}{cc} 0 & -a_2 \\ 1 & -a_1 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] + \left[\begin{array}{c} 1 \\ 0 \end{array}\right] u(t) \nonumber \\ y(t)= \left[\begin{array}{cc} 0 & 1 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] \end{array}\right. }$

に対して，つぎの座標変換を行え。

$\displaystyle{ \left[\begin{array}{c} x_1'(t) \\ x_2'(t) \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ 1 & -a_1 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] }$

解答　状態空間表現に

$\displaystyle{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ 1 & -a_1 \end{array}\right]^{-1} \left[\begin{array}{c} x_1'(t) \\ x_2'(t) \end{array}\right] = \left[\begin{array}{cc} a_1 & 1 \\ 1 & 0 \end{array}\right] \left[\begin{array}{c} x_1'(t) \\ x_2'(t) \end{array}\right] }$

を代入して

$\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{cc} a_1 & 1 \\ 1 & 0 \end{array}\right] \left[\begin{array}{c} \dot{x}_1'(t) \\ \dot{x}_2'(t) \end{array}\right] = \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}{c} x_1'(t) \\ x_2'(t) \end{array}\right] + \left[\begin{array}{c} 1 \\ 0 \end{array}\right] u \nonumber \\ y(t)= \left[\begin{array}{cc} 0 & 1 \end{array}\right] \left[\begin{array}{cc} a_1 & 1 \\ 1 & 0 \end{array}\right] \left[\begin{array}{c} x_1'(t) \\ x_2'(t) \end{array}\right] \end{array}\right. }$

この状態方程式の左から

$\displaystyle{ \left[\begin{array}{cc} a_1 & 1 \\ 1 & 0 \end{array}\right]^{-1} = \left[\begin{array}{cc} 0 & 1 \\ 1 & -a_1 \end{array}\right] }$

をかけて，つぎの状態空間表現を得る。

$\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1'(t) \\ \dot{x}_2'(t) \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ -a_2 & -a_1 \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) \nonumber \\ 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}\right. }$

演習1.8　２次系の状態空間表現

$\displaystyle{ \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] %}_{x} \end{array}\right. \quad(\zeta<1) }$

に対して，つぎの座標変換を行え。

$\displaystyle{ \left[\begin{array}{c} x_1'(t) \\ x_2'(t) \end{array}\right] = \left[\begin{array}{cc} 1 & 0 \\ -\zeta\omega_n & \omega_n\sqrt{1-\zeta^2} \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] }$

演習1.9　２次系の状態空間表現

$\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right] = \left[\begin{array}{cc} a_{11} & a_{12} \\ a_{21} & a_{22} \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] + \left[\begin{array}{c} b_{1} \\ b_{2} \end{array}\right] u(t) \nonumber \\ y(t) = \underbrace{ \left[\begin{array}{cc} c_1 & c_2 \end{array}\right] }_C \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] \quad(c_1\ne0) \end{array}\right. }$

に対して，どのような座標変換

$\displaystyle{ \left[\begin{array}{c} x_1'(t) \\ x_2'(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} t_{11} & t_{12} \\ t_{21} & t_{22} \end{array}\right] }_T \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] }$

を行なえば，出力方程式をつぎの形にできるか。

$\displaystyle{ y(t) = \underbrace{ \left[\begin{array}{cc} 1 & 0 \end{array}\right] }_{C'} \left[\begin{array}{c} x'_1(t) \\ x'_2(t) \end{array}\right] }$

さて，つぎの２つの状態空間表現を考える。

$\displaystyle{(1.6)\quad S_1: \left\{\begin{array}{l} \dot{x}_1(t)=A_1x_1(t)+B_1u_1(t) \\ y_1(t)=C_1x_1(t)+D_1u_1(t) \end{array}\right. }$

$\displaystyle{(1.7)\quad S_2: \left\{\begin{array}{l} \dot{x}_2(t)=A_2x_2(t)+B_2u_2(t) \\ y_2(t)=C_2x_2(t)+D_2u_2(t) \end{array}\right. }$

$S_1$ の入力変数の数と $S_2$ の出力変数の数が等しいとき， $S_2$ の出力を $S_1$ の入力に直接結合することを考える。その状態空間表現は次式で与えられる。

$\displaystyle{(1.8)\quad \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right] = \left[\begin{array}{cc} A_1 & B_1C_2 \\ 0 & A_2 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] + \left[\begin{array}{c} B_1D_2 \\ B_2 \end{array}\right] u_2(t) \\ y_1(t)= \left[\begin{array}{cc} C_1 & D_1C_2 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] +D_1D_2u_2(t) \end{array}\right. }$

例題1.6　１次系の状態空間表現

$\displaystyle{ \left\{\begin{array}{l} \dot{x}(t)=-\frac{1}{T}x(t)+\frac{K}{T}u(t) \\ y(t)=x(t) \end{array}\right. }$

の入力に，むだ時間を１次系として近似したときの状態空間表現

$\displaystyle{ \left\{\begin{array}{l} \dot{x}_L(t)=-\frac{2}{L}x_L(t)+\frac{2}{L}u_L(t) \\ y_L(t)=2x_L(t)-u_L(t) \end{array}\right. }$

の出力を結合して得られるシステムの状態空間表現を求めよ。
解答　まず，状態方程式は

$\displaystyle{ \dot{x}(t)=-\frac{1}{T}x(t)+\frac{K}{T}u(t) }$

に $u(t)=y_L(t)=2x_L(t)-u_L(t)$ を代入し，

$\displaystyle{ \dot{x}_L(t)=-\frac{2}{L}x_L(t)+\frac{2}{L}u_L(t) }$

と合わせて

$\displaystyle{ \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_L(t) \end{array}\right] = \left[\begin{array}{ccc} -\frac{1}{T} & 2\frac{K}{T} \\ 0 & -\frac{2}{L} \end{array}\right] \left[\begin{array}{c} x(t) \\ x_L(t) \end{array}\right] + \left[\begin{array}{c} -\frac{K}{T} \\ \frac{2}{L} \end{array}\right] u_L(t) }$

のように得られる。また，出力方程式は

$\displaystyle{ y(t)= \left[\begin{array}{ccc} 1 & 0 \end{array}\right] \left[\begin{array}{c} x(t) \\ x_L(t) \end{array}\right] }$

のように得られる。

演習1.10 ２次系

$\displaystyle{ \left\{\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 & 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) \nonumber \\ 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}\right. }$

の入力に，むだ時間を２次系として近似したときの状態空間表現

$\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_{L1}(t) \\ \dot{x}_{L2}(t) \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ -\frac{12}{L^2} & -\frac{6}{L} \end{array}\right] \left[\begin{array}{c} x_{L1}(t) \\ x_{L2}(t) \end{array}\right] + \left[\begin{array}{c} 0 \\ 1 \end{array}\right] u_L(t) \nonumber \\ y_L(t) = \left[\begin{array}{cc} 0 & -\frac{12}{L} \end{array}\right] \left[\begin{array}{c} x_{L1}(t) \\ x_{L2}(t) \end{array}\right] +u_L(t) \end{array}\right. }$

の出力を結合して得られるシステムの状態空間表現を求めよ。

MATLABを用いて座標変換を行うために，たとえば演習1.8で， $\omega_n=1,\zeta=-0.01$ の場合，つぎのコマンドを与えればよい。

%coordinate_transformation.m
om=1; zeta=0.01;
A=[0 1;-om^2 -2*zeta*om]; B=[0;om^2]; C=[1 0]; D=0;
sys1=ss(A,B,C,D);
T=[1 0;-zeta*om om*sqrt(1-zeta^2)];
sys2=ss2ss(sys1,inv(T))

SCILABを用いて座標変換を行うために，たとえば演習1.8で， $\omega_n=1,\zeta=-0.01$ の場合，つぎのコマンドを与えればよい。

//coordinate_transformation.m
om=1; zeta=0.01;
A=[0 1;-om^2 -2*zeta*om]; B=[0;om^2]; C=[1 0]; D=0;
sys1=syslin(‘c’,A,B,C,D);
T=[1 0;-zeta*om om*sqrt(1-zeta^2)];
sys2=ss2ss(sys1,T)

つぎに，MATLABを用いて直列結合を行うためには，たとえば演習1.10で $L=1$ の場合，つぎのコマンドを与えればよい。

%serial_connection.m
A1=[0 1;0 0]; B1=[0;1]; C1=[1 0]; D1=0;
sys1=ss(A1,B1,C1,D1);
L=1;
A2=[0 1;-12/L^2 -6/L]; B2=[0;1]; C2=[0 -12/L]; D2=1;
sys2=ss(A2,B2,C2,D2);
sys=sys1*sys2

ここで， sys1 と sys2 の順番に注意する。

最後に，状態空間表現から伝達関数行列表現（8章の(8.2)式を参照。）を求めるには，コマンド tf を用いる。たとえば，上のむだ時間を２次系として近似したときの状態空間表現 sys2 の伝達関数を求めるには，コマンド tf(sys2) を与えればよい。

つぎに，SCILABを用いて直列結合を行うためには，たとえば演習1.10で $L=1$ の場合，つぎのコマンドを与えればよい。

//serial_connection.m
A1=[0 1;0 0]; B1=[0;1]; C1=[1 0]; D1=0;
sys1=syslin(‘c’,A1,B1,C1,D1);
L=1;
A2=[0 1;-12/L^2 -6/L]; B2=[0;1]; C2=[0 -12/L]; D2=1;
sys2=syslin(‘c’,A2,B2,C2,D2);
sys=sys1*sys2

ここで， sys1 と sys2 の順番に注意する。

最後に，状態空間表現から伝達関数行列表現（8章の(8.2)式を参照。）を求めるには，コマンド ss2tf を用いる。たとえば，上のむだ時間を２次系として近似したときの状態空間表現 sys2 の伝達関数を求めるには，コマンド ss2tf(sys2) を与えればよい。

演習問題の解答

演習1.1　 $\omega=\dot{\theta}$ とおくと， $\dot{\theta}=\omega, \dot{\omega}=-\frac{3g}{4\ell}\theta$ 。したがって，つぎの状態方程式を得る。

$\displaystyle{ \left[\begin{array}{c} \dot{\theta} \\ \dot{\omega} \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ -\frac{3g}{4\ell} & 0 \end{array}\right] \left[\begin{array}{c} \theta \\ \omega \end{array}\right] }$

演習1.2　 $v_1=\dot{x}_1, v_2=\dot{x}_2$ とおくと， $\dot{x}_1=v_1, \dot{v}_1=-(x_1-x_2)+u$ ， $\dot{x}_2=v_2, \dot{v}_2=-(x_2-x_1)$ 。したがって，つぎの状態方程式を得る。

$\displaystyle{ \left[\begin{array}{c} \dot{x}_1 \\ \dot{v}_1 \\ \dot{x}_2 \\ \dot{v}_2 \end{array}\right] = \left[\begin{array}{cccc} 0 & 1 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & -1 & 0 \\ \end{array}\right] \left[\begin{array}{c} x_1 \\ v_1 \\ x_2 \\ v_2 \end{array}\right] + \left[\begin{array}{c} 0 \\ 1 \\ 0 \\ 0 \end{array}\right] u }$

または

$\displaystyle{ \left[\begin{array}{c} \dot{x}_1 \\ \dot{x}_2 \\ \dot{v}_1 \\ \dot{v}_2 \end{array}\right] = \left[\begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \\ \end{array}\right] \left[\begin{array}{c} x_1 \\ x_2 \\ v_1 \\ v_2 \end{array}\right] + \left[\begin{array}{c} 0 \\ 0 \\ 1 \\ 0 \end{array}\right] u }$

演習1.3　一定のトルク $\tau^*$ のもとで一定の角速度 $\omega^*$ を得ている。これらは運動方程式 $\dot{\omega}=-{\omega}+\tau$ を満足するので， $0=-\omega^*+\tau^*$ 。これを運動方程式から辺々引き算して，つぎの１次系としての状態空間表現を得る。

$\displaystyle{ \left\{\begin{array}{l} \underbrace{\frac{d}{dt}(\omega-\omega^*)}_{\dot x} =\underbrace{-1}_{a}\underbrace{(\omega-\omega^*)}_{x} +\underbrace{1}_{b} \underbrace{(\tau-\tau^*)}_{u} \\ \underbrace{\omega-\omega^*}_{y} =\underbrace{1}_{c}\underbrace{\omega-\omega^*}_{x} \end{array}\right. }$

演習1.4　 $x_1={y}, x_2=\dot{y}$ とおくと， $\dot{x}_1=x_2, \dot{x}_2=-a_1x_2-a_2x_1+u$ 。したがって，つぎの状態空間表現を得る。

$\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1 \\ \dot{x}_2 \\ \end{array}\right] = \left[\begin{array}{cccc} 0 & 1 \\ -a_2 & -a_1 \\ \end{array}\right] \left[\begin{array}{c} x_1 \\ x_2 \\ \end{array}\right] + \left[\begin{array}{c} 0 \\ 1 \\ \end{array}\right] u \\ y = \left[\begin{array}{cccc} 1 & 0 \end{array}\right] \left[\begin{array}{c} x_1 \\ x_2 \\ \end{array}\right] \end{array}\right. }$

または

$\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_2 \\ \dot{x}_1 \\ \end{array}\right] = \left[\begin{array}{cccc} -a_1 & -a_2 \\ 1 & 0 \\ \end{array}\right] \left[\begin{array}{c} x_2 \\ x_1 \\ \end{array}\right] + \left[\begin{array}{c} 1 \\ 0 \\ \end{array}\right] u \\ y = \left[\begin{array}{cccc} 0 & 1 \end{array}\right] \left[\begin{array}{c} x_2 \\ x_1 \\ \end{array}\right] \end{array}\right. }$

演習1.5　 $\dot{x}=u,\ y=x+u$ より，つぎのブロック線図を得る。

演習1.6　 $\dot{x}_1=x_2,\ \dot{x}_2=-2\zeta\omega_nx_2+\omega_n^2u,\ y=c_2x_2$ より，つぎのブロック線図を得る。

演習1.7　 $\dot{v}_1=-(x_1-x_2)+u,\ \dot{v}_2=-(x_2-x_1)$ より，つぎのブロック線図を得る。

演習1.8　状態空間表現に

$\displaystyle{ \left[\begin{array}{c} x_1 \\ x_2 \end{array}\right] = \left[\begin{array}{cc} 1 & 0 \\ -\zeta\omega_n & \omega_n\sqrt{1-\zeta^2} \end{array}\right] \left[\begin{array}{c} x_1' \\ x_2' \end{array}\right] }$

を代入して

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

この状態方程式の左から

$\displaystyle{ \left[\begin{array}{cc} 1 & 0 \\ -\zeta\omega_n & \omega_n\sqrt{1-\zeta^2} \end{array}\right]^{-1} = \frac{1}{\omega_n\sqrt{1-\zeta^2}} \left[\begin{array}{cc} \omega_n\sqrt{1-\zeta^2} & 0 \\ \zeta\omega_n & 1 \end{array}\right] }$

をかけて，つぎの状態空間表現を得る。

$\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1' \\ \dot{x}_2' \end{array}\right] = \left[\begin{array}{cc} -\zeta\omega_n & \omega_n\sqrt{1-\zeta^2} \\ -\omega_n\sqrt{1-\zeta^2} & -\zeta\omega_n \end{array}\right] \left[\begin{array}{c} x_1' \\ x_2' \end{array}\right] + \left[\begin{array}{c} 0 \\ \frac{1}{\omega_n\sqrt{1-\zeta^2}} \end{array}\right] u \nonumber \\ y= \left[\begin{array}{cc} 1 & 0 \end{array}\right] \left[\begin{array}{c} x_1' \\ x_2' \end{array}\right] \end{array}\right. }$

演習1.9

$\displaystyle{ \underbrace{ \left[\begin{array}{cc} 1 & 0 \end{array}\right] }_{C'} \underbrace{ \left[\begin{array}{cc} t_{11} & t_{12} \\ t_{21} & t_{22} \end{array}\right] }_T = \underbrace{ \left[\begin{array}{cc} c_1 & c_2 \end{array}\right] }_C \quad(Tは正則) }$

を満足させればよいので，たとえば $t_{11}=c_1(\ne0), t_{12}=c_2, t_{21}=0, t_{22}=1$ と選べばよい。

演習1.10

$\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \\\hdashline \dot{x}_{L1}(t) \\ \dot{x}_{L2}(t) \end{array}\right] = \left[\begin{array}{cc:cc} 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & -\frac{12}{L} \\\hdashline 0 & 0 & 0 & 1 \\ 0 & 0 & -\frac{12}{L^2} & -\frac{6}{L} \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \\ \hdashline x_{L1}(t) \\ x_{L2}(t) \end{array}\right] + \left[\begin{array}{c} 0 \\ 1 \\ \hdashline 0 \\ 1 \end{array}\right] u_L(t) \\ y(t)= \left[\begin{array}{cc:cc} 1 & 0 & 0 & 0 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \\ \hdashline x_{L1}(t) \\ x_{L2}(t) \end{array}\right] \end{array}\right. }$

UV Control

投稿日時: 2022年8月17日投稿者: cacsd

●Drone becomes to be very popular tool to take a video from the sky. It is a surprise that we can control its motion with 6 degrees of freedom so smartly. What makes it possible?

In general, a drone has 4 propellers as shown. This means that 4 degrees of freedom can be controlled. They should be the heave $z$ and the roll $\phi$ , the pitch $\theta$ , the yaw $\psi$ . Assume that rotational speeds of 4 propellers can be set to the specified values, $\omega_{1c}$ , $\omega_{2c}$ , $\omega_{3c}$ , $\omega_{4c}$ .

Then we can produce 1 force $F_z$ and 3 torques $\tau_\phi$ , $\tau_\theta$ , $\tau_\psi$ to control the heave $z$ and the roll $\phi$ , the pitch $\theta$ , the yaw $\psi$ respectively by combining the 4 propellers’ thrusts as follows:

$\left[\begin{array}{c} F_z \\ \tau_\phi \\ \tau_\theta \\ \tau_\psi \end{array}\right]= \left[\begin{array}{cccc} k_1 & k_1 & k_1 & k_1 \\ 0 & -\ell k_1 & 0 & \ell k_1 \\ \ell k_1 & 0 & -\ell k_1 & 0 \\ -k_2 & k_2 & -k_2 & k_2 \end{array}\right] \left[\begin{array}{c} \omega_{1c} \\ \omega_{2c} \\ \omega_{3c} \\ \omega_{4c} \end{array}\right]$ .

This means that the MIMO dynamics is decoupled into 4 SISO dynamics. We will employ a PID controller corresponding to each SISO dynamics. Therefore it is the static decoupling to make the control system design easier.

●Many underwater vehicles have been developed for various purposes. For example, we can find them as follows.

Note that 4 to 8 thrustors are equipped to realize the static decoupling. The number depends on motion necessary to achieve the mission of UVs. Most control systems will be designed in the same way as drone.

●Prof. Koterayama and Prof. Nakamura developed the following UV called as DELTA.

The mission of DELTA is for wide-area survey. It works in a towing mode for wade scanning, and in a self-propulsive mode for local investigation. The detail is discussed in the paper:

Masahiko NAKAMURA, Hiroyuki KAJIWARA, Wataru KOTERAYAMA:
Development of an ROV operated both as towed and self-propulsive vehicle
Ocean Engineering, vol.28, no.1, pp.1-43, 2000.

DELTA has only two thrustors. So only the heave $z$ and the yaw $\psi$ are controlled. However it is necessary to control the roll $\phi$ and the pitch $\theta$ to keep the equilibrium state. Such a control system can’t be realized by PID control. The following control system is designed by LQG approach.

For any UV, the collision avoidance is an important function. For example, a deep diving will be necessary to avoid an obstacle as follows.

In order to change the depth in DELTA, there is no way except decreasing the thrust. This means the large speed variation, by which hydrodynamic force will be varied largely. So we need a robust controller or a scheduling controller. Such a control scheme is depicted in the following.

Ship Control

投稿日時: 2022年8月17日投稿者: cacsd

●For the directional control of a ship, consider the following NOMOTO model.

$\displaystyle{ \dot{r}(t)=-\frac{1}{T(t)}r(t)+\frac{K(t)}{T(t)}\delta(t-L)+w(t)} }$

where $r(t)=\dot{\psi}(t)$ and

$\displaystyle{ T(t)=\frac{L}{U(t)}T',\ K(t)=\frac{U(t)}{L}K' \quad (U_1\le U(t)\le U_2) }$

We assume the time lag $L$ to take account of the delay to get the rudder effect. We use the same notation for the ship length $L$ assuming no confusion. The time constant $T$ and gain constant $K$ are varied according to the variation of forward speed $U$ . From the following equation, it is clear that the rudder effect is proportional to the square of the forward speed.

$\displaystyle{ \dot{r}(t)=-\underbrace{\left(\frac{U(t)}{L}\right)\frac{1}{T'}}_{a(t)=\frac{1}{T(t)}}r(t) +\underbrace{\left(\frac{U(t)}{L}\right)^2\frac{K'}{T'}}_{b(t)=\frac{K(t)}{T(t)}}\delta(t-L)+w(t) }$

Given a nominal forward speed $U^*$ , the following equation holds.

$\displaystyle{ \dot{r}(t)=-\underbrace{\left(\frac{U(t)}{U^*}\right)\frac{1}{T^*}}_{a(t)=\frac{1}{T(t)}}r(t) +\underbrace{\left(\frac{U(t)}{U^*}\right)^2\frac{K^*}{T^*}}_{b(t)=\frac{K(t)}{T(t)}}\delta(t-L)+w(t) }$

where

$\displaystyle{ T^*=\frac{L}{U^*}T',\ K^*=\frac{U^*}{L}K' \quad (U_1\le U^*\le U_2) }$

●Zig-Zag Test:　Neglecting the time lag $L$ and the forward speed variation, a simulation example of the Zig-Zag test is given as follows.

From this response, we used to obtain both time constant $T$ and gain constant $K$ . However, the response is not always the step response by which a linear system is completely characterized. As a ship has dynamics of an astatic (no-fixed-position) system, the Zig-Zag test is a compromise way for its identification.

We should notice that if we apply the unity feedback as shown in the above we can obtain the step response and identify $(\zeta,\omega_n)$ from $(T_p, 1+p_0)$ by the following formula.

$\displaystyle{ (T_p, 1+p_0)\ \Rightarrow\ (\zeta,\omega_n)= \left( \sqrt{\frac{(\log p_0)^2}{(\log p_0)^2+\pi^2}}, \frac{\sqrt{(\log p_0)^2+\pi^2}}{T_p} \right) }$

●PID control: Under the nominal forward speed $U^*$ , neglecting the time lag $L$ , consider the following NOMOTO model.

$\displaystyle{ \dot{r}(t)=-\frac{1}{T^*}r(t)+\frac{K^*}{T^*}\delta(t)+w(t) }$

Then PID control is presented as follows.

$\displaystyle{ \delta(t)=K_P(\psi_c-\psi(t))+K_D\frac{d}{dt}(\psi_c-\psi(t))+K_I\int_0^t(\psi_c-\psi(\tau))\,d\tau }$

Although we used to tune the P gain $K_P$ , the D gain $K_D$ , the I gain $K_I$ on the site, here we want try to determine them based on the NOMOTO model. The closed-loop system by PID control is given by

Here assuming $\dot{\psi}_c=0$ , we have

$\displaystyle{ \delta(t)=-K_P\psi(t)-K_Dr(t)+K_P\psi_c+K_I\int_0^t(\psi_c-\psi(\tau))\,d\tau }$

Therefore, we have another configuration of PID control system.

Letting $\zeta=\frac{1}{2\sqrt{T^*K^*}},\omega_n=\sqrt{\frac{K^*}{T^*}}$ , the following state equation is obtained.

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(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} \psi(t) \\ r(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ \omega_n^2 \end{array}\right] }_{B} \underbrace{\delta(t)}_{u(t)} + \left[\begin{array}{c} 0 \\ w \end{array}\right] }$

Then PID control is a state feed back with an integral action as follows.

$\displaystyle{ \underbrace{\delta(t)}_{u(t)} =- \underbrace{ \left[\begin{array}{cc} K_P & K_D \end{array}\right] }_{F} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \end{array}\right] }_{x(t)} +\underbrace{K_P}_{G} \underbrace{\psi_c}_{v} +K_I \underbrace{\int_0^t(\psi_c-\psi(\tau))\,d\tau}_{x_I(t)} }$

Therefore, after specifying the desired first overshoot $(T_p',1+p_0')$ , we calculate the corresponding to $(\zeta',\omega_n')$ and then determine PID gains by

$\displaystyle{ K_P=\frac{\omega_n'^2}{\omega_n^2},\ K_D=\frac{2}{\omega_n^2}(\zeta'\omega_n'-\zeta\omega_n),\ K_I=\omega_I\frac{\omega_n'^2}{\omega_n^2} }$

Here we try $\omega_I=\frac{w_n'}{10}$ at first, and then tune it by observing the disturbance rejection.

The CLPS by PID control is presented by

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \\ \dot{x}_I(t) \end{array}\right] }_{\dot{x}_{CL}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ -\omega_n'^2 & -2\zeta'\omega_n' & \omega_I\omega_n'^2 \\ -1 & 0 & 0 \end{array}\right] }_{A_{CL}} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ x_I(t) \end{array}\right] }_{x_{CL}(t)} + \underbrace{ \left[\begin{array}{cc} 0 & 0\\ \omega_n'^2 & 1\\ 1 & 0 \end{array}\right] }_{B_{CL}} \left[\begin{array}{c} \psi_c \\ w \end{array}\right] }$

Here, assuming $\dot{w}=0$ , we differentiate it and obtain

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \ddot{\psi}(t) \\ \ddot{r}(t) \\ \ddot{x}_I(t) \end{array}\right] }_{\ddot{x}_{CL}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ -\omega_n'^2 & -2\zeta'\omega_n' & \omega_I\omega_n'^2 \\ -1 & 0 & 0 \end{array}\right] }_{A_{CL}} \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \\ \dot{x}_I(t) \end{array}\right] }_{\dot{x}_{CL}(t)} }$

From the CLPS stability, we have

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ x_I(t) \end{array}\right] }_{x_{CL}} \rightarrow \underbrace{ \left[\begin{array}{ccc} 0 & -1 & 0\\ \omega_n'^2 & 2\zeta'\omega_n' & -\omega_I\omega_n'^2 \\ 1 & 0 & 0 \end{array}\right]^{-1} }_{-A_{CL}^{-1}} \underbrace{ \left[\begin{array}{cc} 0 \\ \omega_n'^2\psi_c +w\\ \psi_c \end{array}\right] }_{B_{CL}} = \left[\begin{array}{cc} \psi_c\\ 0\\ \frac{w}{\omega_I\omega_n'^2} \end{array}\right] }$

This means that the integral action is estimating the unknown disturbance.

The following simulation shows the CLPS responses by PD control. It is observed that the disturbance rejection is not sufficient because of the lack of integral action.

The following simulation shows the CLPS responses by PID control. Here $\omega_I=2\frac{w_n'}{10}$ . It is observed that the disturbance rejection is satisfied but the transient response is fractured a little bit being off expectation.

●LQI control: We will try to determine the PID gains in the framework of LQI control. Consider the following error system:

$\displaystyle{ \left[\begin{array}{c} \dot\psi(t)-\dot{\psi}_c \\ \dot r(t) \\ \dot{\delta}(t) \end{array}\right] = \left[\begin{array}{ccc} 0 & 1 & 0 \\ 0 & -2\zeta\omega_n &\omega_n^2 \\ 0 & 0 & 0 \end{array}\right] \left[\begin{array}{c} \psi(t)-\psi_c \\ r(t) \\ \delta(t) \end{array}\right] + \left[\begin{array}{c} 0 \\ 0 \\ 1 \end{array}\right] \dot{\delta} (t) }$

For the error system, we will minimize the following criterion function:

$\displaystyle{ J=\int_0^\infty(\frac{1}{M_\psi^2}(\psi(t)-\psi_c)^2+\frac{T_\psi^2}{M_\psi^2}r^2(t) +\frac{1}{M_\delta^2}\delta^2(t)+\frac{T_\delta^2}{M_\delta^2}\dot{\delta}^2(t))\,dt }$

Solving this problem, the PID gains are obtained as follows:

$\displaystyle{ \dot{\delta}(t) =- \left[\begin{array}{ccc} f_\psi & f_r & f_\delta \end{array}\right] \left[\begin{array}{c} \psi(t)-\psi_c \\ r(t) \\ \delta(t) \end{array}\right] }$
$\displaystyle{ =- \underbrace{ \left[\begin{array}{ccc} f_\psi & f_r & f_\delta \end{array}\right] \left[\begin{array}{ccc} 0 & 1 & 0 \\ 0 & -2\zeta\omega_n &\omega_n^2 \\ 1 & 0 & 0 \end{array}\right]^{-1} }_{ \left[\begin{array}{ccc} K_P & K_D & K_I \end{array}\right] } \left[\begin{array}{c} \dot\psi(t)-\dot{\psi}_c \\ \dot{r}(t) \\ \psi(t)-\psi_c \end{array}\right] }$

In the case of $M_{\psi}=1, T_{\psi}=0.5T, M_{\delta}=0.5, T_{\delta}=0.1T$ , we have the following simulation result, which seems to be very nice. Note that there is no feedforward term compared with PID control.

However, in the case of having the time lag $L=9$ , the CLPS is unstable as follows:

●LPV control:

Assuming the following speed variation, the responses of CLPS by unity feedback vary largely.

Consider the following NOMOTO model based on a nominal forward speed $U^*$ .

$\displaystyle{ \left\{\begin{array}{l} \dot{\psi}(t)=r(t) \\ \displaystyle{\dot{r}(t)=-\left(\frac{U}{U^*}\right)\frac{1}{T^*}r(t) +\left(\frac{U}{U^*}\right)^2\frac{K^*}{T^*}\delta(t)} \end{array}\right. }$

Adding the rudder dynamics,

$\displaystyle{ \dot{\delta}(t)=-\frac{1}{T_a}\delta(t)+\frac{K_a}{T_a}u(t) }$

we have the following state-space representation.

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \\ \dot{\delta}(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ 0 & -\left(\frac{U}{U^*}\right)\frac{1}{T^*} & \left(\frac{U}{U^*}\right)^2\frac{K^*}{T^*} \\ 0 & 0 & -\frac{1}{T_a} \end{array}\right] }_{A(U,U^2)} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ \delta(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ \frac{K_a}{T_a} \end{array}\right] }_{B} u(t) }$
$\displaystyle{ \underbrace{ \psi(t) }_{y(t)} = \underbrace{ \left[\begin{array}{ccc} 1 & 0 & 0 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ \delta(t) \end{array}\right] }_{x(t)} }$

Note that it is possible to have the following polytopic representation.

$\displaystyle{ A(U,U^2)=p_1(U,U^2)A_1+p_2(U,U^2)A_2+p_3(U,U^2)A_3 }$

where

$\displaystyle{A_1=A(U_1,U_1^2)}$
$\displaystyle{A_2=A(U_2,U_2^2)}$
$\displaystyle{A_3=A(\frac{U_1+U_2}{2},U_1U_2)}$ 　

Defining $U_3=\frac{U_1+U_2}{2}$ ,

$\displaystyle{ p_1(U,U^2)=\frac{1}{p_0}\det \left[\begin{array}{cc} U-U_3 & U_2-U_3 \\ U^2-U_1U_2 & U_2^2-U_1U_2 \\ \end{array}\right]}$
$\displaystyle{ p_2(U,U^2)=\frac{1}{p_0}\det \left[\begin{array}{cc} U_1-U_3 & U-U_3 \\ U_1^2-U_1U_2 & U^2-U_1U_2 \\ \end{array}\right]}$
$\displaystyle{ p_3(U,U^2)=\frac{1}{p_0}\det \left[\begin{array}{cc} U_1-U_2 & U_2-U \\ U_1^2-U_2^2 & U_2^2-U^2 \\ \end{array}\right]}$
$\displaystyle{ p_0=\det \left[\begin{array}{cc} U_1-U_2 & U_2-U_3 \\ U_1^2-U_2^2 & U_2^2-U_1U_2 \\ \end{array}\right]}$
$\displaystyle{p_1(U,U^2)+p_2(U,U^2)+p_3(U,U^2)=1}$

Based on this model, the framework of LPV control is proposed. LPV control can cover the variation included in the triangle region with the vertexes $(U_1,U_1^2)$ , $(U_2,U_2^2)$ and temporal vertex $(U_3,U_1U_2)$ .

Constitute the following 2-port system.

$\displaystyle{ P: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x} \\ \dot{x}_I \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A(U,U^2)& 0 \\ -C & 0 \end{array}\right] }_{\cal{A}(U,U^2)} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 1 \end{array}\right] }_{B_1} r + \underbrace{\left[\begin{array}{c} B \\ 0 \end{array}\right] }_{B_2} u\\ \underbrace{ \left[\begin{array}{c} y_{11} \\ y_{12} \end{array}\right] }_{y_1} = \underbrace{ \left[\begin{array}{cc} 0 &\omega_I\\ \omega_DCA(U,U^2) & 0 \end{array}\right] }_{C_1} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \end{array}\right] }_{D_{11}} r + \underbrace{ \left[\begin{array}{c} 0 \\ \omega_DCB \end{array}\right] }_{D_{12}} u\\ \underbrace{ \left[\begin{array}{c} y \\ x_I \end{array}\right] }_{y_2} = \underbrace{ \left[\begin{array}{cc} C & 0\\ 0 & 1 \end{array}\right] }_{C_2} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \end{array}\right] }_{D_{21}} r \end{array}\right.}$

$A(U,U^2)=p_1(U,U^2)A_1+p_2(U,U^2)A_2+p_3(U,U^2)A_3$

For this LPV model, design the following LPV controller.

$\displaystyle{ K_0: \left\{\begin{array}{l} \dot{x}_K=A_K(U,U^2)x_K+ \underbrace{ \left[\begin{array}{cc} B_K^{(1)}(U,U^2) & B_K^{(2)}(U,U^2) \end{array}\right] }_{B_K(U,U^2)} \left[\begin{array}{c} y \\ x_I \end{array}\right] \\ u=C_K(U,U^2)x_K + \underbrace{ \left[\begin{array}{cc} D_K^{(1)}(U,U^2) & D_K^{(2)}(U,U^2) \end{array}\right] }_{D_K(U,U^2)} \left[\begin{array}{c} y \\ x_I \end{array}\right] \end{array}\right.}$

$A_K(U,U^2)=p_1(U,U^2)A_{K1}+p_2(U,U^2)A_{K2}+p_3(U,U^2)A_{K3}$
$B_K(U,U^2)=p_1(U,U^2)B_{K1}+p_2(U,U^2)B_{K2}+p_3(U,U^2)B_{K3}$
$C_K(U,U^2)=p_1(U,U^2)C_{K1}+p_2(U,U^2)C_{K2}+p_3(U,U^2)C_{K3}$
$D_K(U,U^2)=p_1(U,U^2)D_{K1}+p_2(U,U^2)D_{K2}+p_3(U,U^2)D_{K3}$

Form this, the following output feedback is obtained.

$\displaystyle{ K: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_K \\ \dot{x}_I \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A_K(U,U^2) & B_K^{(2)}(U,U^2) \\ 0 & 0 \end{array}\right] }_{A_K(U,U^2)} \left[\begin{array}{c} x_K \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{cc} B_K^{(1)}(U,U^2) & 0\\ -1& 1 \end{array}\right] }_{B_K(U,U^2)} \left[\begin{array}{c} y \\ r \end{array}\right] \\ u= \underbrace{ \left[\begin{array}{cc} C_K(U,U^2) & D_K^{(2)}(U,U^2) \end{array}\right] }_{C_K(U,U^2)} \left[\begin{array}{c} x_K \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{cc} D_K^{(1)}(U,U^2) & 0 \end{array}\right] }_{D_K(U,U^2)} \left[\begin{array}{c} y \\ r \end{array}\right] \end{array}\right. }$

We can obtain $H_\infty$ control by fixing as $U=U^*$ . The following simulation shows the CLPS responses by $H_\infty$ control. It is observed that the responses are slightly varied according to the three kinds of speed variation.

On the other hand, the following simulation shows the CLPS responses by the LPV control. It is observed that the responses are almost same in spite of the three kinds of speed variation.

LQG Control

投稿日時: 2022年8月17日投稿者: cacsd

1. Optimal Control for MIMO system

$\bullet$ Given an nth-order system (satisfying controllability and observability)

(1) $\begin{eqnarray*} &&\dot{x}(t)=Ax(t)+Bu(t),\ y(t)=Cx(t)\ \\ &&(x(t)\in{\rm\bf R}^n,\ u(t)\in{\rm\bf R}^m,\ y(t)\in{\rm\bf R}^p) \end{eqnarray*}$

and the stabilizing state feedback

(2) $\begin{eqnarray*} u(t)=-Fx(t), \end{eqnarray*}$

the closed-loop system is represented by

(3) $\begin{eqnarray*} &&\dot{x}(t)=\underbrace{(A-BF)}_{A_F}x(t),\ y(t)=Cx(t)\\ &&A_F=A-BF:\ stable\ matrix. \end{eqnarray*}$

The behaviors of the state and the input are given by

(4) $\begin{eqnarray*} x(t)=\exp(A_Ft)x(0) \end{eqnarray*}$

and

(5) $\begin{eqnarray*} u(t)=-F\exp(A_Ft)x(0). \end{eqnarray*}$

respectively. Then consider a problem to determine $F$ to minimize the criterion function

(6) $\begin{eqnarray*} J=\int_0^\infty(y^T(t)Qy(t)+u^T(t)Ru(t))\,dt. \end{eqnarray*}$

$\bullet$ The criterion function can be written as

(7) $\begin{eqnarray*} J&=&\int_0^\infty x^T(t)(C^TQC+F^TRF)x(t)\,dt\\ &=&x^T(0)\underbrace{\int_0^\infty \exp(A_F^Tt)(C^TQC+F^TRF)\exp(A_Ft)\,dt}_{\Pi}x(0)\\ &=&{\rm tr}\ \Pi x(0)x^T(0). \end{eqnarray*}$

Note that we have the following constraint on $\Pi$ :

(8) $\begin{eqnarray*} \Pi A_F+A_F^T\Pi+C^TQC+F^TRF=0. \end{eqnarray*}$

It is known that $\Pi>0$ holds because of the closed-loop stability. Taking account of all zero-input responses, instead of (7), we consider to minimize

(9) $\begin{eqnarray*} J={\rm tr}\ \Pi. \end{eqnarray*}$

Therefore, we will minimize

(10) $\begin{eqnarray*} J'={\rm tr}\ \Pi+{\rm tr}\ ((\Pi A_F+A_F^T\Pi+C^TQC+F^TRF)\Gamma) \end{eqnarray*}$

using undetermined multiplier $\Gamma$ and the stability constraint (8). As the necessary conditions, we have

(11) $\begin{eqnarray*} \frac{\partial J'}{\partial\Pi}=I_n+\Gamma A_F^T+A_F\Gamma=0\Rightarrow \Gamma>0, \end{eqnarray*}$

(12) $\begin{eqnarray*} \frac{\partial J'}{\partial F}=2(RF-B^T\Pi)\Gamma=0\Rightarrow F=R^{-1}B^T\Pi, \end{eqnarray*}$

(13) $\begin{eqnarray*} \frac{\partial J'}{\partial\Gamma}=\Pi A_F+A_F^T\Pi+C^TQC+F^TRF=0. \end{eqnarray*}$

Substituting $F=R^{-1}B^T\Pi$ into (13),

(14) $\begin{eqnarray*} \Pi (A-BR^{-1}B^T\Pi)+(A-BR^{-1}B^T\Pi)^T\Pi\\ +C^TQC+(R^{-1}B^T\Pi)^TRR^{-1}B^T\Pi=0 \end{eqnarray*}$

That is, we have the Algebraic Riccati Equation (ARE) on $\Pi$ :

(15) $\begin{eqnarray*} {\Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0}. \end{eqnarray*}$

By solving this, we can obtain $\Pi>0$ and then $F$ by

(16) $\begin{eqnarray*} {F=R^{-1}B^T\Pi}. \end{eqnarray*}$

The control methodology is called as LQ Control.

$\bullet$ The sufficiency is discussed as follows. The following expression can be derived:

(17) $\begin{eqnarray*} y^TQy+u^TRu=(u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x) -\underbrace{\frac{d}{dt}x^T\Pi x}_{2x^T\Pi\dot{x}}. \end{eqnarray*}$

In fact, substituting $\dot{x}=Ax+Bu$ into the right hand side and using the Riccati equation,

(18) $\begin{eqnarray*} &&{\rm RHS}=u^TRu+2x^T\Pi Bu+x^T\Pi BR^{-1}B^T\Pi x-2x^T\Pi(Ax+Bu)\\ &&=u^TRu+x^T(\Pi A+A^T\Pi+C^TQC)x-2x^T\Pi Ax={\rm LHS}. \end{eqnarray*}$

Integrating the both side,

(19) $\begin{eqnarray*} 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. \end{eqnarray*}$

As $x(t)=\exp(A_Ft)x(0)\rightarrow 0\ (t\rightarrow\infty)$ ,

(20) $\begin{eqnarray*} 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). \end{eqnarray*}$

Therefore, it is shown that $u=-R^{-1}B^T\Pi x$ minimizes the criterion function.

Exercise 1
Given $A=\left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right],\ B=\left[\begin{array}{cc} 0 \\ 1 \end{array}\right],\ C=\left[\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right], Q=\left[\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right],R=1$ , obtain the solution $\Pi= \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$ of the Riccati equation $\Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0$ and calculate $F=R^{-1}B^T\Pi$ .

Expanding the Ricatti equation

(21) $\begin{eqnarray*} && \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] \nonumber\\ &&-\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] \nonumber , \end{eqnarray*}$

we have

(22) $\begin{eqnarray*} \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. \nonumber \end{eqnarray*}$

There are the two solutions on $\pi_3$ from $-\pi_3^2+1=0$ , the two solutions on $\pi_2$ from $2\pi_3-\pi_2^2+1=0$ , the one solution $\pi_1$ from $\pi_1-\pi_2\pi_3=0$ . Therefore, we have the four kinds of solution $(\pi_1,\pi_2,\pi_3)$ as follows:

(23) $\begin{eqnarray*} \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} (1) &○\\ (2) &× \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} (3) &×\\ (4) &× \end{array}\right. \end{array}\right. . \end{eqnarray*}$

The only (1) satisfies $\pi_1,\pi_2>0,\ \pi_1\pi_2-\pi_3^2>0$ . Therefore we have

(24) $\begin{eqnarray*} F= \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] . \end{eqnarray*}$

●How to solve ARE
Given the Riccati equation $\Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0$ , consider the Hamilton matrix

(25) $\begin{eqnarray*} {M=\left[\begin{array}{cc} A & -BR^{-1}B^T \\ C^TQC & -A^T \end{array}\right]}. \end{eqnarray*}$

The eigenvalues of the Hamilton matrix with size $n\times n$ are distributed symmetrically to not only the real axis but also the imaginal axis. So there are $n$ stable eigenvalues and $n$ unstable eigenvalues. The eigenvectors corresponding to $n$ stable eigenvalues $\lambda_1,\cdots,\lambda_n$ are obtained as follows:

(26) $\begin{eqnarray*} \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{eqnarray*}$

Based on this, we can obtain $\Pi$ as

(27) $\begin{eqnarray*} \Pi=V_2V_1^{-1}. \end{eqnarray*}$

A program by SCILAB to solve the Riccati equation is given as follows.

//opt.sce
function [F,p]=opt(A,B,C,Q,R)
W=R\B’; [V,R]=spec([A -B*W;-C’*Q*C -A’]);
p=diag(R); [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=W*X;
endfunction
//eof

Solve the Riccati equation in Exercise 1 by using the above program.

A=[0 1;0 0]; B=[0;1]; C=eye(2,2);
Q=diag([1 1]); R=1;
[F,p]=opt(A,B,C,Q,R)
poles=spec(A-B*F)

Exercise 2
Consider the following spring-connected carts as a control object.

The motion is governed by

(28) $\begin{eqnarray*} \left\{\begin{array}{l} m_1\ddot{x}_1(t)=-k(x_1(t)-x_2(t))+f_1(t) \\ m_2\ddot{x}_2(t)=-k(x_2(t)-x_1(t))+f_2(t) \end{array}\right. , \end{eqnarray*}$

where $k$ is a spring constant with the range $0.5<k<2$ . The state equation and output equation are given by

(29) $\begin{eqnarray*} \underbrace{ \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \\ \ddot{x}_1(t) \\ \ddot{x}_2(t) \\ \end{array}\right] }_{\dot{x}(t)} &=& \underbrace{ \left[\begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -\frac{k}{m_1} & -\frac{k}{m_1} & 0 & 0 \\ \frac{k}{m_2} & -\frac{k}{m_2}& 0 & 0 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \dot{x}_1(t) \\ \dot{x}_2(t) \\ \end{array}\right] }_{x(t)}\\ &&+ \underbrace{ \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \\ \frac{k}{m_1} & 0 \\ 0 & \frac{k}{m_2} \end{array}\right] }_{B} \underbrace{ \left[\begin{array}{c} f_1(t) \\ f_2(t) \end{array}\right] }_{u(t)}, \end{eqnarray*}$

(30) $\begin{eqnarray*} \underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] }_{y(t)} &=&- \underbrace{ \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \dot{x}_1(t) \\ \dot{x}_2(t) \\ \end{array}\right] }_{x(t)}. \end{eqnarray*}$

The control purpose is to regulate the zero-input response under the initial condition:

(31) $\begin{eqnarray*} \left[\begin{array}{c} {x}_1(0) \\ {x}_2(0) \\ \dot{x}_1(0) \\ \dot{x}_2(0) \\ \end{array}\right] = \left[\begin{array}{c} 0 \\ 0 \\ 0 \\ \frac{k}{m_2} \end{array}\right] \end{eqnarray*}$

by using the state feedback:

(32) $\begin{eqnarray*} \underbrace{ \left[\begin{array}{c} f_1(t) \\ f_2(t) \end{array}\right] }_{u(t)} &=&- \underbrace{ \left[\begin{array}{cccc} g_{11} & g_{12} & g_{13} & g_{14} \\ g_{21} & g_{22} & g_{23} & g_{24} \end{array}\right] }_{F} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \dot{x}_1(t) \\ \dot{x}_2(t) \\ \end{array}\right] }_{x(t)}. \end{eqnarray*}$

In order to determine the state feedback gain $F$ , for $k$ fixed to an appropriate nominal value $k^*$ , we will minimize the following criterion function:

(33) $\begin{eqnarray*} J=\int_0^\infty (q_1^2x_1^2(t)+q_2^2x_2^2(t)+r_1^2f_1^2(t)+r_2^2f_2^2(t))\,dt, \end{eqnarray*}$

that is

(34) $\begin{eqnarray*} J&=&\int_0^\infty ( \underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right]^T }_{y^T(t)} \underbrace{ \left[\begin{array}{cc} q_1^2 & 0 \\ 0 & q_2^2 \ \end{array}\right] }_{Q} \underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] }_{y(t)}\\ &&+\underbrace{ \left[\begin{array}{c} f_1(t) \\ f_2(t) \end{array}\right]^T }_{u^T(t)} \underbrace{ \left[\begin{array}{cc} r_1^2 & 0 \\ 0 & r_2^2 \ \end{array}\right] }_{R} \underbrace{ \left[\begin{array}{c} f_1(t) \\ f_2(t) \end{array}\right] }_{u(t)} )\,dt. \end{eqnarray*}$

For example, the closed-loop zero-input response is simulated under the following assumptions:

1) $m_1=1$ , $m_2=1$ , $k^*=1$ ,
2) $q_1=\frac{1}{M_{x_1}}=1$ , $q_2=\frac{1}{M_{x_2}}=1$ , $r_1=\frac{\rho}{M_{f_1}}=1$ , $r_2=\frac{\rho}{M_{f_2}}=1$
( $M_{x_1}=M_{x_2}=1, M_{f_1}=M_{f_2}=1,\ \rho=1$ ).

Appendix 1

Check the following properties on matrix trace.

(35) $\begin{eqnarray*} {\rm tr}AB={\rm tr}BA, \end{eqnarray*}$

(36) $\begin{eqnarray*} \frac{\partial}{\partial X}{\rm tr}AXB=A^TB^T, \end{eqnarray*}$

(37) $\begin{eqnarray*} \frac{\partial}{\partial X}{\rm tr}AX^TB=BA, \end{eqnarray*}$

(38) $\begin{eqnarray*} \frac{\partial}{\partial X}{\rm tr}AXBX^T=A^TXB^T+AXB \end{eqnarray*}$

where $\displaystyle{{\rm tr}A=\sum_{i}a_{ii}}$ for $A=[a_{ij}]_{i,j=1,\cdots,n}$ . In fact,

(39) $\begin{eqnarray*} \sum_{i}[AB]_{ii}=\sum_{i}\sum_{j}a_{ij}b_{ji} =\sum_{j}\sum_{i}b_{ji}a_{ij}=\sum_{j}[BA]_{jj}, \end{eqnarray*}$

(40) $\begin{eqnarray*} &&[\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}\\ &&=\sum_{k}b_{jk}a_{ki}=[BA]_{ji}=[A^TB^T]_{ij}, \end{eqnarray*}$

(41) $\begin{eqnarray*} &&[\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}\\ &&=\sum_{k}b_{ik}a_{kj}=[BA]_{ij}, \end{eqnarray*}$

(42) $\begin{eqnarray*} &&\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)\\ &&=AXB+A^TXB^T. \end{eqnarray*}$

LQG Control

投稿日時: 2022年8月17日投稿者: cacsd

1. Optimal Control for a SISO system

Given a 1st-order system with an input and an output:

(1) $\begin{eqnarray*} &&\dot{x}(t)=ax(t)+bu(t),\ y(t)=x(t)\ \\ &&(x(t)\in{\rm\bf R},\ u(t)\in{\rm\bf R},\ y(t)\in{\rm\bf R}) \end{eqnarray*}$

and its stabilizing state feedback:

(2) $\begin{eqnarray*} u(t)=-fx(t), \end{eqnarray*}$

the closed-loop system is represented by

(3) $\begin{eqnarray*} \dot{x}(t)=(a-bf)x(t),\ a-bf<0. \end{eqnarray*}$

The state behavior $x(t)$ and the input behavior $u(t)$ are given by

(4) $\begin{eqnarray*} x(t)=e^{(a-bf)t}x(0) \end{eqnarray*}$

and

(5) $\begin{eqnarray*} u(t)=-fe^{(a-bf)t}x(0). \end{eqnarray*}$

respectively. Then consider a problem to determine $f$ to minimize a criterion function

(6) $\begin{eqnarray*} J=\int_0^\infty(q^2x^2(t)+r^2u^2(t))\,dt. \end{eqnarray*}$

The first term of the criterion function is calculated as

(7) $\begin{eqnarray*} J_x&=&\int_0^\infty q^2x^2(t)\,dt \nonumber\\ &=&\int_0^\infty q^2e^{2(a-bf)t}x^2(0)\,dt \nonumber\\ &=&q^2x^2(0)\left[\frac{1}{2(a-bf)}e^{2(a-bf)t}\right]_0^\infty \nonumber\\ &=&\frac{q^2x^2(0)}{2(a-bf)}\left[\underbrace{e^{2(a-bf)\infty}}_{0}-\underbrace{e^{2(a-bf)0}}_{1}\right] \nonumber\\ &=&-\frac{q^2}{2(a-bf)}x^2(0)>0\quad (a-bf<0), \end{eqnarray*}$

and the second term of the criterion function is calculated as

(8) $\begin{eqnarray*} J_u&=&\int_0^\infty r^2u^2(t)\,dt \nonumber\\ &=&\int_0^\infty r^2f^2e^{2(a-bf)t}x^2(0)\,dt \nonumber\\ &=&r^2f^2x^2(0)\left[\frac{1}{2(a-bf)}e^{2(a-bf)t}\right]_0^\infty \nonumber\\ &=&-\frac{r^2f^2}{2(a-bf)}x^2(0)>0\quad (a-bf<0). \end{eqnarray*}$

Therefore, the criterion function can be written as

(9) $\begin{eqnarray*} J=\underbrace{\frac{-q^2}{2(a-bf)}x^2(0)}_{J_x}+\underbrace{\frac{-r^2f^2}{2(a-bf)}x^2(0)}_{J_u} =\underbrace{-\frac{q^2+r^2f^2}{2(a-bf)}}_{\Pi}x^2(0). \end{eqnarray*}$

Minimizing $J$ is equivalent to minimizing

(10) $\begin{eqnarray*} \Pi=-\frac{q^2+r^2f^2}{2(a-bf)}. \end{eqnarray*}$

Differentiating by $f$ brings

(11) $\begin{eqnarray*} \frac{d\Pi}{df} &=&-\frac{2r^2f}{2(a-bf)}-(-1)\frac{q^2+r^2f^2}{2(a-bf)^2}(-b)\nonumber\\ &=&-\frac{2r^2(a-bf)f+(q^2+r^2f^2)b}{2(a-bf)^2}\nonumber\\ &=&\frac{br^2f^2-2ar^2f-q^2b}{2(a-bf)^2}. \end{eqnarray*}$

Therefore

(12) $\begin{eqnarray*} \frac{d\Pi}{df}=0 & \Rightarrow & bf^2-2af-\left(\frac{q}{r}\right)^2b=0. \end{eqnarray*}$

As $f$ must satisfy $a-bf<0$ , we have

(13) $\begin{eqnarray*} {f=\frac{1}{b}\left(a+\sqrt{a^2+\left(\frac{q}{r}\right)^2b^2}\right)}\nonumber \end{eqnarray*}$

This is uniquely determined because $\Pi$ is downward convex as follows:

(14) $\begin{eqnarray*} \frac{d^2\Pi}{df^2}&=&\frac{2br^2f-2ar^2}{2(a-bf)^2} +(-2)\frac{br^2f^2-2ar^2f-q^2b}{2(a-bf)^3}(-b)\nonumber\\ &=&\frac{-r^2}{(a-bf)} +b\frac{-2r^2(a-bf)f-br^2f^2-q^2b}{(a-bf)^3}\nonumber\\ &=&\frac{-r^2(a-bf)^2-2r^2(a-bf)bf-r^2b^2f^2-q^2b^2}{(a-bf)^3}\nonumber\\ &=&\frac{-r^2(a-bf+bf)^2-q^2b^2}{2(a-bf)^3}\nonumber\\ &=&\frac{-r^2a^2-q^2b^2}{2(a-bf)^3}>0. \end{eqnarray*}$

The closed-loop system by the $f$ is given by

(15) $\begin{eqnarray*} \dot{x}(t)=(a-bf)x(t)=-\sqrt{a^2+\left(\frac{q}{r}\right)^2b^2}x(t). \end{eqnarray*}$

In the case of $a=\frac{1}{T},b=\frac{K}{T}$ ,

(16) $\begin{eqnarray*} \dot{x}(t)&=&-\frac{1}{T}\sqrt{1+\left(\frac{q}{r}\right)^2K^2}x(t)\nonumber\\ &=&-\frac{1}{\frac{T}{\sqrt{1+\left(\frac{q}{r}\right)^2K^2}}}x(t) \end{eqnarray*}$

which means that the new time constant is shorter than the original time constant $T$ .

Exercise 1
Letting $a=-1$ , $b=1$ , $x(0)=1$ , consider the following cases:

Case#1: $\displaystyle{\frac{q}{r}=2}$
Case#2: $\displaystyle{\frac{q}{r}=4}$
Case#3: $\displaystyle{\frac{q}{r}=8}$

Then simulate the behaviors of $x(t)$ and $u(t)$ as follows.

Question:
Why do we use not only $J_x$ but also $J_u$ in the criterion function.
Answer:
Check that $J_u$ is downward convex, and takes the minimum value at $f=\frac{2b}{a}$ as follows:

(17) $\begin{eqnarray*} \frac{dJ_u}{df}&=&\frac{d}{df}\left(-\frac{r^2f^2}{2(a-bf)}x^2(0)\right)\nonumber\\ &=&\frac{-r^2x^2(0)}{2}\left( \frac{2f}{(a-bf)}+(-1)\frac{f^2}{(a-bf)^2}(-b)\right)\nonumber\\ &=&\frac{-r^2x^2(0)}{2}\frac{2(a-bf)f+f^2b}{(a-bf)^2}\nonumber\\ &=&\frac{-r^2x^2(0)}{2}\frac{(2a-bf)f}{(a-bf)^2}, \end{eqnarray*}$

(18) $\begin{eqnarray*} \frac{d^2J_u}{df^2}&=&\frac{-r^2x^2(0)}{2}\left(\frac{2a-2bf}{(a-bf)^2} +(-2)\frac{(2a-bf)f}{2(a-bf)^3}(-b)\right)\nonumber\\ &=&-r^2x^2(0)\left(\frac{1}{a-bf}+b\frac{2(a-bf)f+bf^2}{(a-bf)^3}\right)\nonumber\\ &=&-r^2x^2(0)\frac{(a-bf)^2+2(a-bf)bf+b^2f^2}{(a-bf)^3}\nonumber\\ &=&-r^2x^2(0)\frac{(a-bf+bf)^2}{2(a-bf)^3}\nonumber\\ &=&-r^2x^2(0)\frac{a^2}{2(a-bf)^3}>0. \end{eqnarray*}$

For example, letting $b=1,q=1,r=1$ , for $a=-1$ and $a=1$ , the overview of $J_x$ , $J_u$ , $J=J_x+J_u$ are drown as follows.

Here the symbol “o” shows the minimum of $J$ . Note that $J_u$ is necessary to make $J$ downward convex.

Appendix
In order to extend the above discussion to MIMO systems, we should be familiar with Lagrange’s method of undetermined multipliers. We will rewrite the above discussion by using this method as follows.

From (10), note that the constraint on $\Pi$ is given by the following Lyapnov’s equation

(19) $\begin{eqnarray*} {2(a-bf)\Pi+q^2+r^2f^2=0} \end{eqnarray*}$

Here $\Pi>0$ holds because of $a-bf<0$ . Therefore, instead of minimizing $\Pi$ , we will minimize

(20) $\begin{eqnarray*} J'=\Pi+\Gamma(2(a-bf)\Pi+q^2+r^2f^2). \end{eqnarray*}$

using undetermined multiplier $\Gamma$ and the stability constraint (19). As the necessary conditions, we have

(21) $\begin{eqnarray*} \frac{\partial J'}{\partial \Pi}=1+2(a-bf)\Gamma=0\Rightarrow\Gamma>0, \end{eqnarray*}$

(22) $\begin{eqnarray*} \frac{\partial J'}{\partial f}=(-2b\Pi+2r^2f)\Gamma=0\Rightarrow f=r^{-2}b\Pi, \end{eqnarray*}$

(23) $\begin{eqnarray*} \frac{\partial J'}{\partial \Gamma}=2(a-bf)\Pi+q^2+r^2f^2=0. \end{eqnarray*}$

Substituting $f=r^{-2}b\Pi$ into (23),

(24) $\begin{eqnarray*} 2(a-br^{-2}b\Pi)\Pi+q^2+r^2r^{-4}b^2\Pi^2=0 \end{eqnarray*}$

That is, we have the second-order equation on $\Pi$

(25) $\begin{eqnarray*} {2a\Pi-r^{-2}b^2\Pi^2+q^2=0} \end{eqnarray*}$

which is called as Ricatti equation. By solving this, $\Pi>0$ is obtained by

(26) $\begin{eqnarray*} \Pi=\frac{a+\sqrt{a^2+r^{-2}q^2b^2}}{r^{-2}b^2}, \end{eqnarray*}$

and $f$ is given by

(27) $\begin{eqnarray*} f=r^{-2}b\Pi=\frac{1}{b}\left(a+\sqrt{a^2+\left(\frac{q}{r}\right)^2b^2}\right). \end{eqnarray*}$

Lastly consider the following matrix:

(28) $\begin{eqnarray*} {M=\left[\begin{array}{cc} a & -r^{-2}b^2 \\ -q^2 & -a \end{array}\right]}. \end{eqnarray*}$

which is called as Hamilton matrix. The stable eigenvalue is given by

(29) $\begin{eqnarray*} \lambda=-\sqrt{a^2+r^{-2}b^2q^2} \end{eqnarray*}$

from

(30) $\begin{eqnarray*} {\rm det}(\lambda I_2-M)=\lambda^2-a^2-r^{-2}b^2q^2=0. \end{eqnarray*}$

The corresponding eigenvector is obtained as

(31) $\begin{eqnarray*} \left[\begin{array}{cc} v_1 \\ v_2 \end{array}\right] =\left[\begin{array}{cc} 1 \\ \frac{-a-\sqrt{a^2+r^{-2}b^2q^2}}{-r^{-2}b^2} \end{array}\right] \end{eqnarray*}$

Note that

(32) $\begin{eqnarray*} \Pi=v_2v_1^{-1}=\frac{a+\sqrt{a^2+r^{-2}q^2b^2}}{r^{-2}b^2} \end{eqnarray*}$

and

(33) $\begin{eqnarray*} f=r^{-2}b\Pi=\frac{a+\sqrt{a^2+r^{-2}b^2q^2}}{b}. \end{eqnarray*}$