制御系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. 安定化を適切に行う

投稿日時: 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)の妥当性について説明しなさい。

前章までに，状態フィードバックと状態オブザーバを組み合わせた安定化コントローラ(オブザーバベース\,コントローラ)を，どのような条件のもとで構成できるかについて検討し，それが可制御性と可観測性(可安定性と可検出性)であることを述べた。そこでつぎに，状態フィードバックと状態オブザーバの構成を合理的に行うことを考えたい。

まず，状態フィードバックについては，代表的な方法として最小2乗規範に基づく方法がある。これは，状態フィードバックによる閉ループ系における時間応答の2乗面積が最小となるように，状態フィードバックゲインを定めるものである。ここで，状態変数の時間的振る舞いの2乗面積ばかりでなく，どのような入力を用いて閉ループ系を構成したかを評価するために，入力変数の時間的振る舞いの2乗面積についても同時に考えることが肝要である。こうして得られる制御方式を，状態フィードバックによる「LQ制御」と呼ぶ。この状態フィードバックゲインは，「リッカチ方程式」と呼ばれる行列に関する2次方程式を解いて求めることができる。

つぎに，状態オブザーバについてであるが，オブザーバベース\,コントローラによる閉ループ系における時間応答の2乗面積の和を最小化する立場から考える。この値は，理想的な状態フィードバックだけの場合と比べて，状態オブザーバを用いるので増加(劣化)する。そこで，この増分が小さくなるようにオブザーバゲインを定めることが行われる。

本章で述べる制御方式は「最適制御」と呼ばれるものの一つで，歴史的には「LQG理論」として確率システムの制御の立場から研究されてきた。しかし，ここでは確定的な議論にとどめて，主要な結果を説明した。

4. センサが足りないが大丈夫か

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

前章で，状態フィードバックによって安定化可能である条件(可安定性)や閉ループ系の固有値を任意に設定できる条件(可制御性)を学んだ。ただし，すでに指摘しておいたように，状態フィードバックという制御手段は，すべての状態変数にセンサを付けること( $y=x$ )が前提となる。しかしながら，これはいつも満足されているとは限らないので，状態変数を推定することも考えておかなければならない。現代制御は，その手段として「状態オブザーバ」を提供している。これは一つの動的システムで，これを状態空間表現で表すとき，アナログ回路や計算機で実現される。その入力は制御対象の入力と出力であり，また出力は制御対象の状態の推定値(時間の経過とともに状態の真値に近ずく)である。したがって，状態オブザーバの出力に状態フィードバックのゲイン行列をかけて，状態フィードバックを近似的に実現する方法，すなわち「オブザーバベース\,コントローラ」がよく採用される。
\par
それでは，状態オブザーバは，制御対象の状態空間表現が与えられるとき，いつでも構成可能なのであろうか。答えはノーである。状態オブザーバが構成可能な条件は「可検出性」として知られている。可検出性は，前章で学んだ可安定性と密接な関係(双対関係)にある。これについて学ぶのが本章の一つの目標である。可安定性を強めた概念として可制御性があったように，可検出性を強めた概念として「可観測性」があり，これについても本章で学ぶ。最後に，2～4章の重要な結論として，制御対象の安定化可能条件が可安定性と可検出性であり，その手段がオブザーバベース\,コントローラであることを述べる。

3. いつ安定化できるのか

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

前章で，状態方程式から安定性を調べる方法を学んだ。すなわち，対象が安定であるとは $A$ 行列の固有値が複素左半平面内にあるときであった。それでは，仮に対象が不安定なとき，どのような制御方法で，どのような条件のもとで安定化を図ることができるのだろうか。ここで，不安定な対象とは，自然に存在する物としては考えにくいが，人工物としては，倒立振子や磁気吸引浮上，方向不安定な各種ビークル，核融合炉のプラズマ閉じ込めなどに見られる。また，制御手段としては，「状態フィードバック」というものを考える。これは，すべての状態変数の線形結合を各入力にフィードバックするもので，各状態変数の変化を即座に入力に反映させることができるという意味で，最も強力な手段である。ただし，すべての状態変数にセンサを付けること( $y=x$ )が前提となるので，しばしば理想的な手段である。

さて，状態フィードバックによって安定化可能である条件は「可安定性」と呼ばれるが，これを直接明らかにするのは容易ではない。そこで，状態フィードバックによって $A$ 行列のすべての固有値を複素平面のどこにでも移動できる条件に等価な「可制御性」に着目する。これにより，線形性が保たれる限りにおいては，どのように速い時間応答でも実現できることになる。これに対し，可安定性は， $A$ 行列の不安定な固有値のみを複素左半平面内に移動できればよく，安定な固有値は移動できなくともよい。このような可安定性と可制御性の微妙な違いを，最終的には数学的理論展開を追って理解してほしい。

特異値分解(singular-value decomposition)

——————————————————————————————————————–
定理　 $A$ をサイズ $m\times n$ の行列( ${\rm rank}\,A=k$ )とする。このとき，サイズ $m\times m$ の直交行列 $U$ とサイズ $n\times n$ の直交行列 $V$ が存在して

(1) $\begin{equation*} A=U\Sigma V^T \end{equation*}$

が成り立つ。ここで，サイズ $m\times n$ の行列 $\Sigma$ は，次を満たす（ $X(m\times n)$ は， $X$ がサイズ $n\times n$ の行列であることを示す( $\times$ は行数と列数の区切)。）。

(2) $\begin{equation*} \Sigma= \left[\begin{array}{cc} \Sigma_1(k\times k) & 0\,(k\times n-k)\\ 0\,(m-k\times k) & 0\,(m-k\times n-k) \end{array}\right] \end{equation*}$

(3) $\begin{equation*} \Sigma_1={\rm diag}\{\sigma_1,\cdots,\sigma_k\} \quad (\sigma_1\ge\cdots\ge\sigma_k>0) \end{equation*}$

を満足するものである。
証明　 $A^TA\ge 0$ だから，仮定 ${\rm rank}\,A=k$ より， $A^TA$ は $k$ 個の正固有値と $n-k$ 個の零固有値をもち，互いに直交する固有ベクトルをもつ。そこで， $A^TA$ の固有値の正の平方根を，大きい順に， $\sigma_1\ge\cdots\ge\sigma_k>\sigma_{k+1}=\cdots=\sigma_{n}=0$ のように表し，対応する固有ベクトル $v_1,\cdots,v_n$ を $v_i^Tv_i=1\ (i=j),\ 0\ (i\ne j)$ を満足するようにとることができる。いま， $\Sigma_1$ を式3のように，また

(4) $\begin{equation*} V= [\begin{array}{cc} V_1 & V_2 \end{array}] = [\begin{array}{ccc|ccc} v_1 & \cdots & v_k & v_{k+1}& \cdots & v_n \end{array}] \end{equation*}$

とおくと， $V$ は直交行列となり，つぎが成り立つ。

(5) $\begin{equation*} A^TAV_1=V_1\Sigma_1^2,\quad \end{equation*}$

(6) $\begin{equation*} A^TAV_2=0 \end{equation*}$

式6の第2式の左から， $V_2^T$ をかけて

(7) $\begin{equation*} V_2^TA^TAV_2=0\Rightarrow AV_2=0 \end{equation*}$

また， $U_1=AV_1\Sigma_1^{-1}$ とおくと，式6の第1式から $U_1^TU_1=I_{k}$ を得る。そこで， $U_2$ を $U= [\begin{array}{cc} U_1 & U_2 \end{array}]$ が直交行列となるように選ぶと

(8) $\begin{equation*} U^TAV= \left[\begin{array}{cc} U_1^TAV_1 & U_1^TAV_2\\ U_2^TAV_1 & U_2^TAV_2 \end{array}\right] = \left[\begin{array}{cc} \Sigma_1 & 0\\ U_2^TU_1\Sigma_1 & 0 \end{array}\right] =\Sigma \end{equation*}$

が成り立つ。これより，式1を得る。
——————————————————————————————————————–

——————————————————————————————————————–
例題　行列 $A= \left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 1 \end{array}\right]$ の特異値分解は

(9) $\begin{equation*} A= \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right] }_{U} \underbrace{ \left[\begin{array}{ccc} \sqrt{2} & 0 & 0 \\ 0 & 1 & 0 \end{array}\right] }_{\Sigma} \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0 \\ \frac{1}{\sqrt{2}} & 0 & -\frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} \end{array}\right]^T }_{V^T} \end{equation*}$

のように与えられることを確かめなさい。
解答　 $A^TA$ のサイズは $3\times 3$ であるが， $AA^T$ のサイズは $2\times 2$ であるので， $AA^T$ を計算すると $\left[\begin{array}{cc} 1 & 0 \\ 0 & 2 \end{array}\right]$ となる（サイズ $m\times n$ の行列 $A$ の特異値を手計算で求めるには， $A^TA$ と $AA^T$ のが同じ非零固有値をもつことから，サイズの小さいほうの固有値計算を行えばよい。）。これから， $AA^T$ の固有値は $1,2$ で，その正の平方根 $1,\sqrt{2}$ が特異値であるが，定理??の $\Sigma_1$ の対角成分の特異値は大きい順に並べる約束であるので

(10) $\begin{equation*} \underbrace{ \left[\begin{array}{cc} 1 & 0 \\ 0 & 2 \end{array}\right] }_{AA^T} \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right] }_{U} = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right] }_{U} \underbrace{ \left[\begin{array}{cc} 2 & 0 \\ 0 & 1 \end{array}\right] }_{\Sigma_1^2} \end{equation*}$

のように， $U$ と $\Sigma_1$ が決まる。つぎに， $V$ については，視察によって

(11) $\begin{equation*} \underbrace{ \left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 1 \\ 0 & 1 & 1 \end{array}\right] }_{A^TA} \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0 \\ x & 0 & -y \\ y & 0 & x \end{array}\right] }_{V} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0 \\ x & 0 & -y \\ y & 0 & x \end{array}\right] }_{V} \underbrace{ \left[\begin{array}{ccc} 2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{array}\right] }_{{\rm diag}\{\Sigma_1^2,0\}} \end{equation*}$

としてよい。ここで $\left[\begin{array}{ccc}x \\y \end{array}\right]$ と $\left[\begin{array}{ccc}-y \\x \end{array}\right]$ が直交しており， $x^2+y^2=1$ の制約がある。式(11)から $x=y$ が出て， $x=y=\frac{1}{\sqrt{2}}$ と定まる。
——————————————————————————————————————–

それでは，階段化アルゴリズムについて述べる。

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

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

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

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

(12) $\begin{equation*} B^{(j)}=U^{(j)}\Sigma^{(j)}V^{(j)}\,^T \end{equation*}$

(13) $\begin{equation*} \Sigma^{(j)}= \left[\begin{array}{cc} \Sigma_1^{(j)} & 0 \\ 0 & 0 \end{array}\right] \end{equation*}$

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

ステップ3 座標変換

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

(14) $\begin{equation*} W^{(j)}= \left[\begin{array}{cc} I_s & 0 \\ 0 & U^{(j)} \end{array}\right] \end{equation*}$

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

(15) $\begin{equation*} A^{(j+1)}=W^{(j)}\,^TA^{(j)}W^{(j)},\ T^{(j+1)}=T^{(j)}W^{(j)} \end{equation*}$

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

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

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

——————————————————————————————————————–

——————————————————————————————————————–
例題
階段化アルゴリズムを用いて，つぎの可安定性を調べなさい。

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

解答階段化アルゴリズムの各ステップは，つぎのように計算される。

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

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

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

(19) $\begin{equation*} A^{(2)}= \left[\begin{array}{cc|c} 0 & 0 & 1 \\ 0 & 0 & 0 \\ \hline 0 & 0 & -1 \end{array}\right] \end{equation*}$

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

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

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

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

まず，可安定性と可制御性は座標変換により不変であることを示す。
——————————————————————————————————————–
定理 $(A,B)$ が可安定対[可制御対]ならば， $(TAT^{-1},TB)$ も可安定対[可制御対]である。
証明 $n$ 次系 $\dot{x}(t)=Ax(t)+Bu(t)$ に対する $u(t)=-Fx(t)$ による安定化可能性は，座標変換後の $\dot{x}'(t)=TAT^{-1}x'(t)+TBu(t)$ に対する $u(t)=-FT^{-1}x'(t)$ による安定化可能性と等価である。実際

(21) $\begin{equation*} {\rm det}(\lambda-TAT^{-1}+TB\cdot FT^{-1})= {\rm det}(\lambda-A+BF) \end{equation*}$

が成り立ち，閉ループ系の固有値は，座標変換によって不変だからである。

可制御性については，シルベスターの公式

${\rm rank}\,X+{\rm rank}\,Y-n\le{\rm rank}\,XY\le{\rm min}\{{\rm rank}\,X,{\rm rank}\,Y\}$ ( $n$ は $X$ の列数)

を用いて

(22) $\begin{equation*} \begin{array}{ll} {\rm rank}\, [\begin{array}{cccc} TB & TAT^{-1}TB & \cdots & (TAT^{-1})^{n-1}TB \end{array}] \\ ={\rm rank}\, [\begin{array}{cccc} B & AB &\cdots & A^{n-1}B \end{array}] \end{array} \end{equation*}$

が成り立つことから，定理の主張を得る。
——————————————————————————————————————–

この定理に基づいて，可安定性や可制御性の判別をしやすいように座標変換行列 $T$ を選ぶことを考える。すなわち， $n$ 次系 $\dot{x}(t)=Ax(t)+Bu(t)$ に対して，つぎを満足するように座標変換できるかどうかを調べる。

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

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

「式(23)のように座標変換できれば， $A_k$ が安定行列のとき可安定性が成り立ち，そうでないときは可安定性は成り立たない。また，式(23)のように座標変換できなければ，可制御性が成り立っている。」

つぎの結果は，階段化アルゴリズムを用いて得られる。
——————————————————————————————————————–
定理
サイズ $n\times n$ の行列 $A$ とサイズ $n\times m$ の行列 $B$ に対して，つぎのように変換する直交行列 $T$ が存在する。

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

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

(25) $\begin{equation*} {\rm rank}\,B_j=m_j\quad (j=1,\cdots,k-1) \end{equation*}$

(26) $\begin{equation*} {\rm rank}\,B_k=m_k\ または\ 0 \end{equation*}$

(27) $\begin{equation*} m\ge m_1 \ge \cdots \ge m_{k-1} \ge m_k \end{equation*}$

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

条件C0:　式(24)において， $B_k\ne 0$

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

条件S0:　式(24)において， $B_k=0$ のとき， $A_k$ は安定行列
——————————————————————————————————————–

$B_k=0$ のとき，式(23)のような形式が得られていることは，式\(ref{eq3.4.3})の $T$ を上の $T^{-1}=T^T$ と読み替えれば明らかであろう。また，可安定性と可制御性と相違は，前者が $B_k\ne0$ ばかりでなく $B_k=0$ を許すところにあることに注意する。可制御性の条件C0と可安定性の条件S0の妥当性は，あとで示す。

それでは，まず可制御性について，新しい条件を追加した定理を述べる。
——————————————————————————————————————–
定理　可制御性が成り立つための等価な条件

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

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

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

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

条件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$ のすべての固有値)

ここで，条件C4を満足する $A$ の固有値 $\lambda$ は可制御固有値，満足しない固有値は不可制御固有値と呼ばれる。
——————————————————————————————————————–

前節までに，定義 $\Leftrightarrow$ 条件C1 $\Leftrightarrow$ 条件C2 $\Rightarrow$ 条件C3を示したので，以下では，条件C3 $\Rightarrow$ 条件C0 $\Rightarrow$ 条件C4 $\Leftrightarrow$ 条件C5 $\Rightarrow$ 条件C2を示す。

——————————————————————————————————————–
証明＜条件C3 $\Rightarrow$ 条件C0＞
$B_k=0$ とすると，式(24)を用いて

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

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

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

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

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

(30) $\begin{equation*} \left[\begin{array}{cc} T^TB & T^TAT-\lambda I_n \end{array}\right] = T^T \left[\begin{array}{cc} B & A-\lambda I_n \end{array}\right] \left[\begin{array}{cc} I_n & 0 \\ 0 & T \end{array}\right] \end{equation*}$

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

(31) $\begin{equation*} {\rm rank}\, [\begin{array}{cc} T^TB & T^TAT-\lambda I_n \end{array}] = {\rm rank}\, [\begin{array}{cc} B & A-\lambda I_n \end{array}] =n \end{equation*}$

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

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

(32) $\begin{equation*} w^T\, [\begin{array}{cc} B & A-\lambda I_n \end{array}] =0 \Rightarrow w =0 \end{equation*}$

すなわち

(33) $\begin{equation*} B^Tw=0,\ A^Tw=\lambda w \Rightarrow w=0 \end{equation*}$

と等価である。

証明　＜条件C5 $\Rightarrow$ 条件C2＞
条件C5より

(34) $\begin{equation*} w^TA^{i-1}B=0\quad (i=0,1,2,\cdots) \end{equation*}$

を得て

(35) $\begin{equation*} w^T\, [\begin{array}{ccccc} B & AB & \cdots & A^{n-1}B \end{array}] =0 \Rightarrow w =0 \end{equation*}$

が成り立つ。これは条件C2を意味する。
——————————————————————————————————————–

それでは，可安定性が成り立つための条件について考えよう。可制御であれば可安定であるので，問題は，不可制御で可安定の場合をどのように特徴付けるかである。不可制御であれば条件C0が不成立なので， $B_k=0$ となり， $A_k$ の固有値だけは状態フィードバックによってどうすることもできない。このとき可安定であるためには， $A_k$ が安定行列であることを前提とするしかない。これは条件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$ \quad (〃)

証明　定義 $\Leftrightarrow$ 条件S0だから，条件S1 $\Leftrightarrow$ 条件S2との等価性を示せばよい。そのために， $B_k=0$ のとき，式(29)の行列のランクが落ちるのは

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

(37) $\begin{equation*} {\rm rank}\, \left. \left[\begin{array}{c} B_k^T \\ A_k^T-\lambda I_{m_k} \end{array}\right] \right|_{B_k=0} = {\rm rank}\, \left[\begin{array}{c} 0 \\ A_k^T-\lambda I_{m_k} \end{array}\right] \le m_k \end{equation*}$

の部分で， $\lambda$ が $A_k$ の固有値に一致した場合である。したがって，条件S0を保証するためには， $\lambda$ に $A_k$ のすべての不安定固有値を入れて，式(37)を調べ，列フルランクとなればよい。すなわち，条件C4 $\Leftrightarrow$ 条件C5において， $\lambda$ は $A_k$ のすべての不安定固有値という条件をつければよい。
——————————————————————————————————————–

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

例えば，可制御性は状態フィードバックにより不変であることを示そう。

——————————————————————————————————————–
定理　 $(A,B)$ が可制御対ならば $(A-BF,B)$ も可制御対である。
証明条件C5より，任意の $\lambda$ に対して

(38) $\begin{equation*} B^Tw=0,\ A^Tw=\lambda w \Rightarrow\ w=0 \end{equation*}$

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

(39) $\begin{equation*} B^Tw=0,\ (A-BF)^Tw=\lambda w \Rightarrow\ w=0 \end{equation*}$

を示せばよい。実際， $B^Tw=0$ より

(40) $\begin{equation*} (A-BF)^Tw=A^Tw-F^TB^Tw=A^Tw=\lambda w \end{equation*}$

を得るが，このような $w$ は，仮定(38)より， $w=0$ でなければならない。
——————————————————————————————————————–

2. 状態方程式から安定性を調べる

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

いま対象が物理的な釣り合いの状態にあるとき，何らかの原因でこれが乱されたとしよう。このとき，対象に働きかけて，元の物理的な釣り合いの状態に速やかに戻すことは，最も基本的な制御目的として知られている。もし，対象自身に復元力が備わっていて，元の状態に戻ろうとする性質があれば，対象は「安定」(正確には「漸近安定」)と呼ばれる。このことを状態方程式$\dot{x}=Ax+Bu$から，どのように調べるかが，本章のテーマである。

そのために，まず状態方程式を解くことから始める。その解は，与えられた入力のもとでの状態の時間的振る舞い，すなわち「時間応答」を表している。これは，「零入力応答」と「零状態応答」の和として表される。ここで，零入力応答とは，入力を零すなわち$u=0$として，初期状態$x(0)\ne0$から出発して得られる時間応答であり，零状態応答とは，初期状態$x(0)=0$のもとで適当な入力が与えられるときの時間応答である。対象が物理的な釣り合いの状態にあって，これが何らかの原因で乱されたときの時間的振る舞いは，零入力応答に対応する。そこで，1次系$\dot{x}=ax$に対しては，$a<0$ならば安定である，また，$n$次系$\dot{x}=Ax$に対しては，「正方行列$A$の固有値の実部がすべて負であるならば安定」と判定すればよいことを学ぶ。

一方，零状態応答の中で重要なのは「インパルス応答」と呼ばれるもので，これが与えられると，どのような入力が与えられても，たたみこみ積分を行って出力が求まることを理解してほしい。

1. 物理法則から状態方程式を導く

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

私たちが制御したい対象が物である場合，その振る舞いは物理法則により記述される。一般に，動的な振る舞いを表す物理法則は微分方程式で表される。微分方程式というと尻込みする人もいるかもしれないが，最初は，ニュートンの運動第2法則 $m\ddot{x}=F$ に代表されるような，簡単な線形の微分方程式から始めればよい。ただし，これ1本で済むというわけにはいかない。通常は複数本の微分方程式が出てきて，やっかいなことにそれらがお互いに関連し合っている。どんな問題を解く場合も対象の性質をよく分析することが第一歩であるから，連立された微分方程式を扱える統一的表現法があると望ましい。そこで，対象を記述するいくつかの微分方程式を，1階の連立微分方程式としてまとめて，行列表現することがよく行われる。これを制御対象の「状態方程式」と呼ぶ。また，その解変数を「状態変数」という。

さて，制御を行うには対象の動的な振る舞い(状態変数の時間変化)を観測することが必要であり，そのためにセンサが用いられる。ここで，現実には，すべての状態変数にセンサを配置できるわけではなく，通常は状態変数の一部がセンサにより観測される。センサにより観測される変数と状態変数の間の関係式を，制御対象の「出力方程式」と呼ぶ。

本章では，動的な振る舞いが連立線形微分方程式で記述されるいくつかの制御対象に対して，状態方程式と出力方程式をペアにした「状態空間表現」を求める。これを図示した「ブロック線図」は，計算機シミュレーションを行う際にも有用である。

1.1 状態空間表現の導出例
1.1.1 ペースメーカ

高齢化社会の到来に伴い，より優れた福祉・医療機器の開発が工学分野の大きなテーマの一つとなっている。図1.1に示すのは，心臓のペースメーカの簡単な原理図である。これは，まず左側の閉回路でコンデンサへの充電を行い，つぎにスイッチを切り替えてできる右側の閉回路で放電を行うという動作を周期的に繰り返すことにより，心臓のペースメーカの役割を果たそうとするものである。ここでは，状態方程式を導く最初の例として，このようなＲＣ回路における充電と放電について考える。

そのために，キルヒホッフの電圧則より，左側閉回路と右側閉回路の回路方程式を考えると，それぞれ

(1) $\begin{equation*} R_1\underbrace{C\dot{v}(t)}_{i(t)}+v(t)=u(t) \end{equation*}$

(2) $\begin{equation*} R_2\underbrace{C\dot{v}(t)}_{i(t)}+v(t)=0 \end{equation*}$

となる。ここで， $C$ はコンデンサの容量， $R_1$ は充電回路の抵抗， $R_2$ は心臓の等価抵抗， $v(t)$ は時刻 $t$ におけるコンデンサ両端の電圧， $u(t)$ は時刻 $t$ におけるバッテリの電圧である。

図1.1　心臓のペースメーカ

式(1)は，すでに， $v(t)$ に関する1階の線形微分方程式であるので，両辺を $R_1C$ で割って，つぎの状態方程式を得る。この解変数 $v(t)$ を状態変数と呼ぶ。

(3) $\begin{equation*} \dot{v}(t)=-\frac{1}{R_1C}v(t)+\frac{1}{R_1C}u(t) \end{equation*}$

状態方程式(3)を図1.2のように図示し，これを状態方程式に基づくブロック線図と呼ぶ。この描き方のポイントは，式(3)の右辺を表すのに加え合わせ記号○を用いることと，また $\dot v$ を積分して $v$ を得て右辺と左辺を関連付けていることである。なお，加え合わせにおけるプラス符号は省略することが多い。

図1.2 ペースメーカの充電回路のブロック線図

このブロック線図から，外部より与えられる入力変数 $u(t)$ が，状態変数 $v(t)$ の微分値に影響を与え， $v(t)$ が外部に取り出されることが見てとれる。状態変数は1個であるので，式(3)で表される動的システムを1次システム(first-order system)または1次系 $^*$ と呼ぶ。

同様に，式(2)から得られる状態方程式は

(4) $\begin{equation*} \dot{v}(t)=-\frac{1}{R_2C}v(t) \end{equation*}$

であり，これによるブロック線図は図1.3のように示される。

図1.3 ペースメーカの放電回路のブロック線図

微分方程式(4)の解が

(5) $\begin{equation*} \displaystyle{v(t)=e^{-\frac{1}{R_2C}t}v(0)} \end{equation*}$

と与えられることはよいであろう(式(4)に代入して確かめよ)。状態方程式(4)は入力変数をもたないが，状態変数の初期値によって，状態変数の時間的振る舞いが現れる。この意味で，1次系(4)は自励系(autonomous system) 自由系(unforced system)と呼ばれる。つぎのシミュレーション例 $^{**}$ をみてみよう。

シミュレーション1.1 式(5)で表されるコンデンサ電圧 $v〔{\rm V}〕$ の時間的振る舞いを， $R_2C=1$ ， $v(0)=1\,{\rm V}$ の場合について図1.4に示す。

図1.4 コンデンサ放電時の電圧変化

問1.1 図1.4において，時刻 $t=1,3,5$ における $v$ の値を

(6) $\begin{equation*} e^{x}=1+x+\frac{1}{2}x^2+\cdots+\frac{1}{k!}x^k+\cdots \end{equation*}$

によって近似計算しなさい。

*系はsystemの訳語。ここでは「××システム」を簡潔に「××系」と書く。
**本書では，時間応答のコンピュータによるシミュレーション(simulation)の欄を設けた。最終的には時間応答の数学的理解が大切であるが，まずは，なぜそのような時間的振る舞いが現れるのかを物理的イメージをもって考えながら，典型的な時間応答に親しみをもってほしい。なお，本書の数値計算については演習問題の【4】を参照のこと。

1.1.2 教室のドア

教室で物の動きを実感できるものに，図1.5に示すようなばねとダンパ $^*$ からなる緩衝装置を付けたドアがある。これは，開いたドアをできるだけ速やかに静かに閉めるためのものである。

図1.5 緩衝装置をつけたドア

このドアの運動は回転運動であるが，話しをわかりやすくするため，図1.6に示すような等価な直線運動として調べてみよう。その出発点は，ニュートンの運動第2法則

(7) $\begin{equation*} M\ddot{r}(t)=F(t) \end{equation*}$

である。ここで， $M$ はドアの質量， $r(t)$ は時刻 $t$ におけるドアの変位， $F(t)$ は時刻 $t$ においてドアに働く力であり

(8) $\begin{equation*} F(t)=-D\dot{r}(t)-Kr(t)+u(t) \end{equation*}$

のように表すことができる。ここで，ダンパが第1項の力を，ばねが第2項の力を与える。 $u$ は人がドアに与える力である。式(7)と式(8)より

(9) $\begin{equation*} M\ddot{r}(t)=-D\dot{r}(t)-Kr(t)+u(t) \end{equation*}$

を得る。

図1.6 ドアの簡単なモデル

これは2階の線形微分方程式であるが， $v(t)=\dot{r}(t)$ を定義すると

(10) $\begin{equation*} \dot{r}(t)=v(t) \end{equation*}$

(11) $\begin{equation*} \dot{v}(t)=-\frac{D}{M}v(t)-\frac{K}{M}r(t)+\frac{1}{M}u(t) \end{equation*}$

のような1階の連立線形微分方程式で表される。これらを行列表示すると

(12) $\begin{equation*} \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 u(t) \end{equation*}$

のような状態方程式を得る $^{**}$ 。ここで，状態変数は $r(t)$ と $v(t)$ ，入力変数は $u(t)$ である。また，図1.7のようなブロック線図が得られる。

図1.7 ドアのブロック線図

さて，2個の状態変数のうち，ドアの変位 $r$ の $c_r$ 倍の電圧 $y_r$ ，すなわち

(13) $\begin{equation*} y_r(t)=c_r r(t) \end{equation*}$

を得るセンサはあるが，ドアの速度を計測するセンサはないものとする。このとき， $y_r$ を出力変数と呼ぶ。これは，つぎの出力方程式により表される。

(14) $\begin{equation*} \underbrace{y_r(t)}_{y(t)}= \underbrace{ \left[\begin{array}{cc} c_r & 0 \end{array}\right] }_C \underbrace{ \left[\begin{array}{c} r(t) \\ v(t) \end{array}\right] }_{x(t)} \end{equation*}$

以上から，ドアに対して，状態方程式(12)と出力方程式(14)からなる2次系(second-order system)としての状態空間表現を得た。

シミュレーション 式(12)において, $M=1\,{\rm kg}$ , $K=1\,{\rm N/m}$ ， $r(0)=0.5\,{\rm m}$ ， $v(0)=0\,{\rm m/s}$ ， $u(t)=0\,{\rm N}$ のとき， $D=5,2,0.1\,{\rm N\cdot s/m}$ の三つの場合について，ドア開度 $r〔{\rm m}〕$ の時間的振る舞いを図1.8に示す。

図1.8 ドア開度の時間的振る舞い

問1.2 図1.8の三つの時間応答に対応して，ドアはそれぞれどのように閉まるか説明しなさい。

*ばねとダンパの特性値を調整するためのねじを回すことにより行われる。
**本書では， $\underbrace{\Delta}_{\circ}$ のように書いて，△を○で定義・表記する(△は○に等しいとする)。

1.1.3 直流モータ

代表的なアクチュエータとしてモータがある。例えば図1.9に示すのは，ロボットアームを駆動する直流モータである。

図1.9　直流モータ

このモデルは図1.10のように表される。

図1.10　直流モータのモデル

このとき，つぎが成り立つ。

(15) $\begin{equation*} J\dot{\omega}(t)=K_Ti(t)-\tau_L(t) \end{equation*}$

(16) $\begin{equation*} L\dot{i}(t)=-Ri(t)-K_E\omega(t)+v(t) \end{equation*}$

ここで，式(15)は機械系としての運動方程式であるが，電流による発生トルクの項 $K_T\,i$ を含む。 $K_T$ はトルク定数と呼ばれる。また，式(16)は電気系としての回路方程式であるが，角速度 $\omega=\dot{\theta}$ による逆起電力の項 $K_E\,\omega$ を含む。 $K_E$ は逆起電力定数と呼ばれる。このように，モータは機械系と電気系の混合系という特徴をもつ。式(15)と式(16)に

(17) $\begin{equation*} \dot{\theta}(t)=\omega(t) \end{equation*}$

を加えたものを行列表示すると

(18) $\begin{eqnarray*} %\begin{equation} %\renewcommand{\arraystretch}{1.3} %\begin{array}{rl} \underbrace{ \left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & J & 0 \\ 0 & 0 & L \end{array}\right] }_{E} \left[\begin{array}{c} \dot{\theta}(t) \\ \dot{\omega}(t) \\ \dot{i}(t) \end{array}\right] & = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & K_T \\ 0 & -K_E & -R \end{array}\right] }_{A} \left[\begin{array}{c} \theta(t) \\ \omega(t) \\ i(t) \end{array}\right] \nonumber\\ & + \underbrace{ \left[\begin{array}{cc} 0 & 0 \\ 0 & -1 \\ 1 & 0 \end{array}\right] }_{B} \left[\begin{array}{cc} v(t) \\ \tau_L(t) \end{array}\right] %\end{array} %\end{equation} \end{eqnarray*}$

となる $^*$ 。この左から， $E^{-1}$ をかけて

(19) $\begin{equation*} \underbrace{ \left[\begin{array}{c} \dot{\theta}(t) \\ \dot{\omega}(t) \\ \dot{i}(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & \frac{K_T}{J} \\ 0 & -\frac{K_E}{L} & -\frac{R}{L} \end{array}\right] }_{E^{-1}A} \underbrace{ \left[\begin{array}{c} \theta(t) \\ \omega(t) \\ i(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{cc} 0 & 0 \\ 0 & -\frac{1}{J} \\ \frac{1}{L} & 0 \end{array}\right] }_{E^{-1}B} \underbrace{ \left[\begin{array}{cc} v(t) \\ \tau_L(t) \end{array}\right] }_{u(t)} \end{equation*}$

のような状態方程式を得る。状態方程式(19)は二つの入力変数 $v,\,\tau_L$ をもち， $v$ は操作できるが， $\tau_L$ は操作できない外乱であることに注意してほしい。

問1.3 式(19)を用いて，直流モータのブロック線図を描きなさい。

さて，この直流モータに対しては，角度 $\theta$ の $c_\theta$ 倍の電圧 $y_\theta$ と，角加速度 $\omega$ の $c_\omega$ 倍の電圧 $y_\omega$ が測れるものとすると，出力方程式は

(20) $\begin{equation*} \underbrace{ \left[\begin{array}{c} y_\theta(t) \\ y_\omega(t) \end{array}\right] }_{y(t)}= \underbrace{ \left[\begin{array}{ccc} c_\theta & 0 & 0 \\ 0 & c_\omega & 0 \end{array}\right] }_C \underbrace{ \left[\begin{array}{c} \theta(t) \\ \omega(t) \\ i(t) \end{array}\right] }_{x(t)} \end{equation*}$

となる。このように，直流モータに対して，状態方程式(19)と出力方程式(20)からなる3次系(third-order system)としての状態空間表現を得た。つぎに一つのシミュレーション例を示す。

シミュレーション1.3 式(19)において， $J=1.5\times10^{-5}\,{\rm kg\cdot m^2}$ ， $K_T=0.03\,{\rm N\cdot m/A}$ ， $L=0.001\,{\rm H}$ ， $R=1\,{\rm \Omega}$ ， $K_E=0.05\,{\rm V\cdot s/rad}$ のとき， $\tau_L(t)=0\,{\rm N}$ ， $v(t)=1\,{\rm V}$ ， $x(0)=0$ の場合の $\omega〔{\rm rad/s}〕$ と $i〔{\rm A}〕$ の時間的振る舞いを図1.11に示す

図1.11　直流モータの時間応答

ところで，私たちは物理的な感覚として，機械的な動きと電気的な動きでは速さが格段に違うことを知っている。直流モータは機械系と電気系の混合系であることを述べたが，制御目的は位置制御や速度制御のように機械系に関わるのが普通であるので，状態変数としては $\theta$ と $\omega$ だけでよさそうである。式(16)をみると，直流モータの電気的時定数( $i$ の時定数 $^{**}$ )は

(21) $\begin{equation*} T_e=\frac{L}{R} \end{equation*}$

で与えられ，上の例では $T_e=0.001\,{\rm s}$ である。ところが，図1.11からわかるように， $\omega$ の時定数は約 $0.01\,{\rm s}$ である。したがって，電流は角速度に比べて10倍速く落ち着くので，式(16)の左辺を零とおいてみよう。すなわち

(22) $\begin{equation*} 0=-Ri(t)-K_E\omega(t)+v(t) \end{equation*}$

これから $i$ を求めて，式(15)に代入してみると

(23) $\begin{equation*} \dot{\omega}(t)=-\frac{K_TK_E}{JR}\omega(t) +\frac{K_T}{JR}v(t) -\frac{1}{J}\tau_L(t) \end{equation*}$

を得る。ここで， $\omega$ の時定数

(24) $\begin{equation*} T_m=\frac{RJ}{K_TK_E} \end{equation*}$

は直流モータの機械的時定数と呼ばれている。上の例で計算してみると $T_m=0.01\,{\rm s}$ である。したがって，もし，直流モータの電気的時定数が機械的時定数に比べて十分小さい場合(経験則は $T_e<0.1\,T_m$ )は，式(17)と式(23)を合わせて，つぎの状態方程式をもつ2次系としてよい。

(25) $\begin{equation*} \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{1}{T_m} \end{array}\right] }_A \underbrace{ \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right] }_{x(t)}+ \underbrace{ \left[\begin{array}{cc} 0 & 0 \\ \frac{1}{T_mK_E} & -\frac{1}{J} \end{array}\right] }_B \underbrace{ \left[\begin{array}{c} v(t)\\ \tau_L(t) \end{array}\right] }_{u(t)} \end{equation*}$

式(19)と比較すると，状態空間表現の次数を1だけ減らしたことになる。
これは，モデルの低次元化の一例である。

低次元化の過程を図1.12～図1.14に示しておく。

図1.12　式(1.19)に基づく低次元化前のブロック線図

図1.13　式(1.22)を用いた低次元化中のブロック線図

図1.14　式(1.22)を用いた低次元化中のブロック線図

*式(18)は，式(19)のように物理パラメータどうしの演算を含まず，それらの変動の影響を考察するのに便利な形式であり，ディスクリプタ形式の状態方程式と呼ばれる。
**ここでは，2.1.3項で学ぶ時定数の知識を前提にしている。

1.2 状態空間表現へのモデリング

前節で述べた内容は，モデリング(modeling)という枠組でとらえられる。これは，何らかの対象に対して，目的に応じた科学的分析を通してそのモデルを得ようとするものである。私たちが扱う対象は，入力と出力をもつ動的システム $^*$ (dynamical system)であり，そのモデリングを行うのは動的な振る舞いの制御のためであるので，制御対象(controlled object)と呼ばれる。ここでの科学的分析とは，おもに物理法則から線形定係数常微分方程式 $^{**}$ で表される運動方程式や回路方程式の導出と，これに含まれるパラメータ値の決定をいう。そして，モデルの一つとして，状態空間表現(state-space description)という数理モデル $^{***}$ (mathematical model)を考え，その構築例を示してきた。

本書では，状態空間表現を次式のように書く。

(26) $\begin{equation*} \dot{x}(t)=Ax(t)+Bu(t) \end{equation*}$

(27) $\begin{equation*} y(t)=Cx(t)+Du(t) \end{equation*}$

ここで， $x(t)$ は $n$ 個の状態変数(state variables)をもつベクトル， $u(t)$ は $m$ 個の入力変数(input variables)をもつベクトル， $y(t)$ は $p$ 個の出力変数(output variables)をもつベクトルである。また， $A$ ， $B$ ， $C$ ， $D$ は，それぞれサイズ $n\times n$ ， $n\times m$ ， $p\times n$ ， $p\times m$ の実行列であり，それぞれ $A$ 行列，\mbox{\boldmath $B$ }行列， $C$ 行列， $D$ 行列と呼ばれる。状態方程式(state equation)(26)と出力方程式}(output equation)(27)からなる状態空間表現を， $m$ 入力 $p$ 出力をもつ $n$ 次系( $n$ th-order system)と参照する。 $Du(t)$ の項は直達項(direct part)を表し，存在しないこともある。特に，1次系のときは， $\dot{x}(t)=ax(t)+bu(t)$ などと書く。

*動的システムは，微分方程式・差分方程式のどちらで記述されるかによって連続時間系・離散時間系，重ね合わせの原理が成り立つか否かによって線形系・非線形系，常微分方程式か偏微分方程式かによって集中定数系・分布定数系，係数パラメータの時間依存性によって時変系・時不変系，入出力が確率過程であるか否かによって決定系・確率系などに分類される。
**非線形系の場合の取り扱いは7章で述べる。1～6章までは線形時不変系のみを扱う。
***他の数理モデルとして伝達関数表現がある。状態空間表現と伝達関数表現の間の相互関係については8章で述べる。
****他のアプローチとして，入力と出力の時系列データからモデリングを行うシステム同定がある。

1.3 状態空間表現の座標変換

状態空間表現を見やすくする一つの手段として，座標変換(coordinate transformation)があるので，これについて説明しよう。
いま， $n$ 次系

(28) $\begin{equation*} \dot{x}(t)=Ax(t)+Bu(t) \end{equation*}$

(29) $\begin{equation*} y(t)=Cx(t) \end{equation*}$

に対して，つぎの座標変換を行いたい。

(30) $\begin{equation*} x'(t)=Tx(t) \end{equation*}$

ただし， $T$ は正則とする。式(30)を式(28)に代入すると

(31) $\begin{equation*} \dot{x}(t)=T^{-1}\dot{x'}(t) \end{equation*}$

に注意して

(32) $\begin{equation*} T^{-1}\dot{x'}(t)=AT^{-1}x'(t)+Bu(t) \end{equation*}$

%すなわち

(33) $\begin{equation*} \dot{x'}(t)=\underbrace{TAT^{-1}}_{A'}x'(t)+\underbrace{TB}_{B'}u(t) \end{equation*}$

となる。また，式(30)を式(29)に代入すると

(34) $\begin{equation*} y(t)=\underbrace{CT^{-1}}_{C'}x'(t) \end{equation*}$

となる。この結果を，参照しやすいようにつぎにまとめておく。

定理1.1 $n$ 次系 $\dot{x}(t)=Ax(t)+Bu(t),\ y(t)=Cx(t)$ に対して，座標変換 $x'(t)=Tx(t)$ を行うと，新しい $n$ 次系は次式で表される。

(35) $\begin{equation*} \dot{x'}(t)=A'x'(t)+B'u(t) \end{equation*}$

(36) $\begin{equation*} y(t)=C'x'(t) \end{equation*}$

ただし

(37) $\begin{equation*} A'=TAT^{-1},\ B'=TB,\ C'=CT^{-1} \end{equation*}$

例題1.1 直流モータの状態方程式(25)において， $\tau_L$ を零とおくと

(38) $\begin{equation*} \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{1}{T_m} \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}{T_mK_E} \end{array}\right] }_B \underbrace{ v(t) }_{u(t)} % \end{equation*}$

である。これに対して，座標変換

(39) $\begin{equation*} \underbrace{ \left[\begin{array}{c} x'_1(t) \\ x'_2(t) \end{array}\right] }_{x'(t)} = \underbrace{ \left[\begin{array}{cc} 1 & T_m \\ 0 & -T_m \end{array}\right] }_T \underbrace{ \left[\begin{array}{c} \theta(t) \\ \omega(t) \end{array}\right] }_{x(t)} % \end{equation*}$

を行うと，新しい状態方程式は

(40) $\begin{equation*} \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 & -\frac{1}{T_m} \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} \frac{1}{K_E} \\ -\frac{1}{K_E} \end{array}\right] }_{B'} \underbrace{ v(t) }_{u(t)} \end{equation*}$

となることを示しなさい。

解答座標変換後の $A$ 行列と $B$ 行列は，定理1.1を用いて

(41) $\begin{equation*} \underbrace{ \left[\begin{array}{cc} 1 & T_m \\ 0 & -T_m \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ 0 & -\frac{1}{T_m} \end{array}\right] \left[\begin{array}{cc} 1 & 1 \\ 0 & -\frac{1}{T_m} \end{array}\right] }_{TAT^{-1}} = \underbrace{ \left[\begin{array}{cc} 0 & 0 \\ 0 & -\frac{1}{T_m} \end{array}\right] }_{A'} \end{equation*}$

(42) $\begin{equation*} \underbrace{ \left[\begin{array}{cc} 1 & T_m \\ 0 & -T_m \end{array}\right] \left[\begin{array}{cc} 0 \\ \frac{1}{T_mK_E} \end{array}\right] }_{TB} = \underbrace{ \left[\begin{array}{c} \frac{1}{K_E} \\ -\frac{1}{K_E} \end{array}\right] }_{B'} \end{equation*}$

のように得られる。

ここで，2次系の状態方程式が，二つの1次系の状態方程式

(43) $\begin{equation*} \dot{x}'_1(t)=\frac{1}{K_E}v(t),\quad \dot{x}'_2(t)=-\frac{1}{T_m}x'_2(t)-\frac{1}{K_E}v(t) \end{equation*}$

に分離されており，入力から状態変数への影響の考察をしやすくなっていることに注意してほしい。

1.4 状態空間表現の直列結合

制御対象の状態空間表現を求める際に，図1.15に示すように，二つの部分システムの状態空間表現を求めておいて，これらを直列結合(serial connection)する場合がある。このときの結合システムの状態空間表現を求めることを考える。

図1.15 直列結合( $u_1=y_2$ )

まず，その結果を定理の形で示そう。

定理1.2 二つの状態空間表現

(44) $\begin{equation*} \dot x_1(t)=A_1x_1(t)+B_1u_1(t) \\ \end{equation*}$

(45) $\begin{equation*} y_1(t)=C_1x_1(t)+D_1u_1(t) \end{equation*}$

および

(46) $\begin{equation*} \dot x_2(t)=A_2x_2(t)+B_2u_2(t) \end{equation*}$

(47) $\begin{equation*} y_2(t)=C_2x_2(t)+D_2u_2(t) % \end{equation*}$

に対して， $u_1=y_2$ のように直列結合した場合の状態空間表現は

(48) $\begin{equation*} \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}{cc} B_1D_2 \\ B_2 \end{array}\right] u_2(t) \end{equation*}$

(49) $\begin{equation*} 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{equation*}$

のように得られる。

証明 $\dot x_1(t)=A_1x_1(t)+B_1u_1(t)$ と $y_1(t)=C_1x_1(t)+D_1u_1(t)$ に， $u_1(t)=y_2(t)=C_2x_2(t)+D_2u_2(t)$ を代入して

(50) $\begin{equation*} \dot x_1(t)=A_1x_1(t)+B_1(C_2x_2(t)+D_2u_2(t)) \end{equation*}$

(51) $\begin{equation*} y_1(t)=C_1x_1(t)+D_1(C_2x_2(t)+D_2u_2(t)) \end{equation*}$

となる。第1式と $\dot x_2(t)=A_2x_2(t)+B_2u_2(t)$ をまとめたものと，第2式から，定理の結果を得る。

例題1.2 2次系の制御対象

(52) $\begin{equation*} \dot x(t)= \left[\begin{array}{cc} 0 & 1 \\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] x(t)+ \left[\begin{array}{cc} 0 \\ K\omega_n^2 \end{array}\right] u(t) \end{equation*}$

(53) $\begin{equation*} y(t)=\left[\begin{array}{cc} 1 & 0 \end{array}\right]x(t) \end{equation*}$

に対して( $x(t)$ は2次元ベクトル)，1次系のアクチュエータ

(54) $\begin{equation*} \displaystyle{\dot{x}_a(t) =-\frac{1}{T_a}x_a(t) +\frac{K_a}{T_a}u_a(t)} \end{equation*}$

(55) $\begin{equation*} y_a(t)=x_a(t) \end{equation*}$

を， $u=y_a$ のように直列結合した場合の状態空間表現を求めなさい。

解答定理1.2を用いて，直列結合の状態空間表現として

(56) $\begin{equation*} \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_a(t) \end{array}\right]= \left[\begin{array}{cc|c} 0 & 1 & 0\\ -\omega_n^2 & -2\zeta\omega_n & K\omega_n^2 \\\hline 0 & 0 & -\frac{1}{T_a} \end{array}\right] \left[\begin{array}{c} x(t) \\ x_a(t) \end{array}\right]+ \left[\begin{array}{cc} 0 \\ 0 \\\hline \frac{K_a}{T_a} \end{array}\right] u_a(t) \end{equation*}$

(57) $\begin{equation*} y(t)= \left[\begin{array}{cc|c} 1 & 0 & 0 \end{array}\right] \left[\begin{array}{c} x(t) \\ x_a(t) \end{array}\right] \end{equation*}$

が得られる $^*$ 。

問1.4 例題1.2の直列結合の状態空間表現を，状態ベクトルが $[x_a(t)\ x(t)]^T$ となるように求めなさい。

*ここで， $A$ 行列の縦線と横線， $B$ 行列の横線は，状態ベクトルの要素 $x$ ， $x_a$ のサイズに適合するように引かれている。

演習問題

【1】いろいろな計測装置の基礎となる電気回路の一つにブリッジ回路がある。
例えば，図1.16に示すブリッジ回路 $^*$ を考えてみよう。この回路方程式は

(58) $\begin{equation*} L\dot{i}(t)=-\frac{R_1+R_3}{R_1+R_2}R_2i(t)-\frac{R_1+R_3}{R_1+R_2}v(t) +u(t)-R_3C\dot{v}(t) \end{equation*}$

(59) $\begin{equation*} C\dot{v}(t)=-\left(\frac{R_2}{R_1+R_2}-\frac{R_4}{R_3+R_4}\right)i(t) -\left(\frac{1}{R_1+R_2}+\frac{1}{R_3+R_4}\right)v(t) \end{equation*}$

で与えられる。いま，ブリッジ条件

(60) $\begin{equation*} R_1R_4=R_2R_3 \end{equation*}$

が成り立つとして，つぎの状態方程式を導出しなさい。

(61) $\begin{equation*} \underbrace{ \left[\begin{array}{c} \dot{i}(t) \\ \dot{v}(t) \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} -\frac{1}{L}(\frac{R_1R_2}{R_1+R_2}+\frac{R_3R_4}{R_3+R_4})\ \ 0 \\ 0\ \ -\frac{1}{C}(\frac{1}{R_1+R_2}+\frac{1}{R_3+R_4}) \end{array}\right] }_A \underbrace{ \left[\begin{array}{c} i(t) \\ v(t) \end{array}\right] }_{x(t)}+ \underbrace{ \left[\begin{array}{c} \frac{1}{L} \\ 0 \end{array}\right] }_B u(t) \end{equation*}$

この状態方程式に基づいて，平衡ブリッジ回路のブロック線図を描きなさい。

図1.16　ブリッジ回路

【2】さまざまな柔軟構造物の制振問題は，重要な制御のテーマである。
その特徴は，図1.17に示す連結台車 $^{**}$ にもみられる。この運動方程式は

(62) $\begin{equation*} m_1\ddot{x}_1(t)=-k(x_1(t)-x_2(t))+f_1(t)+w_1(t) \end{equation*}$

(63) $\begin{equation*} m_2\ddot{x}_2(t)=-k(x_2(t)-x_1(t))+w_2(t) \end{equation*}$

で与えられる。ここで， $m_1$ と $m_2$ はそれぞれ台車1と台車2の質量， $k$ はばね定数である。このとき，つぎの状態方程式を導出しなさい。

(64) $\begin{equation*} \begin{array}{rcl} \underbrace{ \left[\begin{array}{c} \dot{x_1}(t) \\ \dot{x_2}(t) \\ \dot{v_1}(t) \\ \dot{v_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) \\ v_1(t) \\ v_2(t) \end{array}\right] }_{x(t)} \\ & & \hspace*{-1zw} + \underbrace{ \left[\begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & 0 \\ \frac{1}{m_1} & 0 & \frac{1}{m_1} \\ 0 & \frac{1}{m_2} & 0 \end{array}\right] }_B \underbrace{ \left[\begin{array}{c} w_1(t)\\ w_2(t)\\ f_1(t) \end{array}\right] }_{u(t)} \end{array} \end{equation*}$

この状態方程式に基づいて，連結台車のブロック線図を描きなさい。

図1.17 連結台車

【3】式23で表される直流モータにおいて，一定入力 $v^*$ ，一定負荷 $\tau_L^*$ のもとで，一定角速度 $\omega^*$ の平衡状態が達成されているものとする。この平衡状態を基準とする直流モータの時間的振る舞いを表す状態方程式を示しなさい。

【4】本書におけるすべての数値計算は，対話型の行列計算環境である学生版MATLAB $^{***}$ を用いて行っている。また，すべての時間応答のグラフは，(非線形)微分方程式による対話型シミュレーション環境である学生版SIMULINK $^{****}$ を用いて得ている。時間応答のシミュレーションのためには，状態方程式のブロック線図を描くことが必要となる。例えば，心臓のペースメーカのブロック線図(図1.3)を得たとすると，SIMULINKでは，これを図1.18のようにほぼそのままの構成で，対話型操作により表現する。ブロックIntegratorの初期値とブロックGainの値を設定し，微分方程式のソルバーの種類，サンプリング周期，シミュレーション時間などを設定すれば，ブロックScopeに図1.1の時間応答を直ちにみることができる。時系列データの処理やグラフ化はMATLABで行える。

MATLABとSIMULINKが手元にあれば，シミュレーション1.3と同一条件下で，直流モータの低次元化後の状態方程式25による角速度の応答を，低次元化前の状態方程式19によるものと比較しなさい。

図1.18 SIMULINKによる微分方程式のブロック表現

*高橋・有本：回路網とシステム理論，コロナ社 (1974)のpp.65 $\sim$ 66から引用。
**B.Wie, D.２.Bernstein : Benchmark Problems for Robust Control Design, ACC Proc. pp.2047 $\sim$ 2048 (1992) から引用。
***The Student Edition of MATLAB-Version\,5 User’s Guide, Prentice Hall (1997)
****The Student Edition of SIMULINK-Version\,2 User’s Guide, Prentice Hall (1998)

まえがき

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

本書は，動く物の制御に関心をもつ工学系の学生の皆さんに，状態方程式に基づく制御について説明するために書かれている。言うまでもないことだが，何らかの動く物を制御しようとすれば，それは物理法則に逆らってはできない。したがって，制御の出発点は，対象の動きを支配する運動方程式や回路方程式であり，これは一般には非線形の微分方程式で表される。どのような制御を行うかについてはさまざまなものが考えられるが，最も基本的なものは，ある物理的な釣り合いの状態(平衡状態)を保持することであり，この目的は安定化と呼ばれている。そのためには，すべての状態変数の変化をセンサによって即座にとらえて対象を操作する仕組み(状態フィードバック)と，またセンサの不足に備えて状態変数を推定する工夫(状態オブザーバ)が基本となる。

ところで，対象の時間的振る舞いは，ある平衡状態の近くでは，重ね合わせの原理がほぼ成り立ち，線形の微分方程式で十分近似できることが多い。これをラプラス変換して得られる伝達関数を前提にして，周波数領域で制御方式を考えるのが「古典制御」である。これに対し，状態空間法と呼ばれる本書のアプローチは，運動方程式や回路方程式から状態方程式を導いて，時間領域で制御方式を考えるもので，「現代制御」として知られている。近年は，古典制御と現代制御は相補的な関係を一段と強め，より高度に発展させられている。本シリーズでは，古典制御については，「フィードバック制御入門」で説明されている。また，その編集委員会によって，現代制御の線形理論について説明する本書のタイトルは，「線形システム制御入門」と与えられ，解析より設計の話にやや重点をおくことが意図されている。

以上のことをふまえて，本書では，現代制御の考え方を説明するためのストーリをつぎのように組み立ててみた。

1章　物理法則から状態方程式を導く
2章　状態方程式から安定性を調べる
3章　いつ安定化できるのか(可制御性と状態フィードバックの話)
4章　センサが足りないが大丈夫か(可観測性と状態オブザーバの話)
5章　安定化を適切に行う(最適制御の話)
6章　定値外乱を抑制する
7章　非線形の運動方程式から始める
8章　伝達関数から状態方程式を導く

ここで，1章から6章までは，2.3節の一部を除いて，対象が線形微分方程式で表される線形システムだけを扱い，7章で非線形システムへの対応について，8章で伝達関数と状態方程式の間の関係について述べた。したがって，微分積分，線形代数，力学，電気回路に関する基礎知識があれば，本書の大部分は古典制御よりも先に読めるであろう。つぎに，記述にあたっては，1次系について導入を行ったうえで，その結果を一般の高次系に拡張するというスタイルをとった。特に2次系についての説明は丁寧に行った。さらに，いくつかの節には$^*$印を付けて，最初はスキップしてもよいことを示した。

本書を，半期(約12時間)の現代制御のテキストとして使う場合は，1～5章と8章をじっくり勉強する基礎重視のコースと，1～7章で1次系と2次系に関する議論をおもにフォローし，高次系についての結果を受け入れていく応用重視のコースの2通りが考えられる。後者においては，実際の制御系設計で必要となる数値計算の指針を演習問題の形で示したので参考にしていただきたい。

おわりに，編集委員の先生方には，著者の拙い原稿に何度も目を通していただき，数多くの貴重なコメントをいただいた。このご支援に対し，心より感謝の意を表したい。それにもかかわらず，著者の力不足のため，現代制御のテキストとして至らないところがあることを恐れている。読者の率直なご意見を，{\tt kajiwara@nams.kyushu-u.ac.jp}まで，お寄せいただければ幸いである。

2000年3月
梶原宏之

8 最小実現問題

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

【本章のねらい】
・伝達関数行列に対応する状態空間表現（実現）を求める。
・実現の可制御可観測な部分系を取り出して，与えられた伝達関数行列に対応する最小実現を求める。

8.1 実現問題

状態空間表現

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

から伝達関数行列表現

$\displaystyle{(8.2)\quad \left\{\begin{array}{l} Y(s)=G(s)U(s)\\ G(s)=C(sI_n-A)^{-1}B+D \end{array}\right. }$

を求める計算を，つぎのドイルの記法を用いて表す。

$\displaystyle{(8.3)\quad \left[\begin{array}{c|c} A & B \\\hline C & D \end{array}\right] =C(sI_n-A)^{-1}B+D }$

ここでは，逆に伝達関数行列表現(8.2)から状態空間表現(8.1)を求める問題、~~すなわち~~実現問題を考える。すなわち

$\displaystyle{(8.4)\quad G(s)= \left[\begin{array}{c|c} A & B \\\hline C & D \end{array}\right] }$

を満足する行列 $A$ ， $B$ ， $C$ ， $D$ を求めたい。この状態空間表現は伝達関数行列表現の実現と呼ばれるが，一意には決まらず，座標変換と状態空間の次元数などに関して自由度がある。次式は座標変換によって伝達関数行列は不変であることを表している。

$\displaystyle{(8.5)\quad \left[\begin{array}{c|c} A & B \\\hline C & D \end{array}\right]= \left[\begin{array}{c|c} TAT^{-1} & TB \\\hline CT^{-1} & D \end{array}\right] }$

また，１入力１出力系の場合は， $G(s)=G^T(s)$ だから，次式が成り立つ。

$\displaystyle{(8.6)\quad \left[\begin{array}{c|c} A & B \\\hline C & D \end{array}\right]= \left[\begin{array}{c|c} A^T & C^T \\\hline B^T & D \end{array}\right] }$

以下では，まず１入力１出力系の場合を考える。１次系について，次式が成り立つ。

$\displaystyle{(8.7)\quad \frac{b}{s+a}+d = \left[\begin{array}{c|c} -a & b \\\hline 1 & d \end{array}\right]= \left[\begin{array}{c|c} -a & 1 \\\hline b & d \end{array}\right] }$

２次系について，次式が成り立つ。

$\displaystyle{ \frac{b_1s+b_2}{s^2+a_1s+a_2}+d }$

$\displaystyle{(8.8)\quad = \left[\begin{array}{cc|c} -a_1 & -a_2 & 1 \\ 1 & 0 & 0 \\\hline b_1 & b_2 & d \end{array}\right] = \left[\begin{array}{cc|c} 0 & 1 & 0 \\ -a_2 & -a_1 & 1 \\\hline b_2 & b_1 & d \end{array}\right] }$

$\displaystyle{(8.9)\quad = \left[\begin{array}{cc|c} -a_1 & 1 & b_1 \\ -a_2 & 0 & b_2 \\\hline 1 & 0 & d \end{array}\right] = \left[\begin{array}{cc|c} 0 & -a_2 & b_2 \\ 1 & -a_1 & b_1 \\\hline 0 & 1 & d \end{array}\right] }$

ここで，(8.8)式の２つの実現は，

$\displaystyle{ T= \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right] }$

について，(8.5)式の関係にある。また，(8.9)式は(8.8)式を(8.6)式の意味で転置したものである。

一般に，１入力１出力 $n$ 次線形系の伝達関数がプロパー $^*$ な有理関数で与えられるとき，これに対応する状態空間表現の例を次に示す $^{**}$ 。

$\displaystyle{(8.10)\quad G(s)=\frac{b_1s^{n-1}+\cdots+b_n}{s^n+a_1s^{n-1}+\cdots+a_n}+d }$
$\displaystyle{ = \left[\begin{array}{cccc|c} -a_1 & \cdots & -a_{n-1} & -a_n & 1 \\ 1 & \cdots & 0 & 0 & 0 \\ \vdots & \ddots & \vdots & \vdots & 0 \\ 0 & \cdots & 1 & 0 & 0 \\ \hline b_1 & \cdots & b_{n-1} & b_n & d \\ \end{array}\right] = \left[\begin{array}{cccc|c} 0 & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & 0 \\ 0 & 0 & \cdots & 1 & 0 \\ -a_n & -a_{n-1} & \cdots & -a_1 & 1 \\ \hline b_n & b_{n-1} & \cdots & b_1 & d \\ \end{array}\right] }$
$\displaystyle{ = \left[\begin{array}{cccc|c} -a_1 & 1 & \cdots & 0 & b_1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ -a_{n-1} & 0 & \cdots & 1 & b_{n-1} \\ -a_n & 0 & \cdots & 0 & b_n \\\hline 1 & 0 & \cdots & 0 & 0 \end{array}\right]\ = \left[\begin{array}{cccc|c} 0 & \cdots & 0 & -a_n & b_n \\ 1 & \cdots & 0 & -a_{n-1} & b_{n-1} \\ \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & \cdots & 1 & -a_1 & b_1 \\\hline 0 & \cdots & 0 & 1 & 0 \end{array}\right] }$

*分子多項式の次数が分母多項式の次数を超えない場合をいう。
**『線形システム制御入門』の8.1節のファディーブの公式を用いて証明できる。

例題8.1 むだ時間要素 $u(t)=v(t-L)$ の伝達関数は $e^{-Ls}$ で与えられ，しばしば次式で近似される。

$\displaystyle{ e^{-Ls}\simeq\frac{1-0.5Ls}{1+0.5Ls} \nonumber }$

このとき，対応する状態空間表現の一つを求めよ。

解答
$\displaystyle{ \frac{1-0.5Ls}{1+0.5Ls}=\frac{2}{1+0.5Ls}-1= \left[\begin{array}{c|c} -\frac{2}{L} & \frac{2}{L} \\\hline 2 & -1 \end{array}\right] }$

演習8.1 むだ時間伝達関数のつぎの近似式に対応する状態空間表現の一つを求めよ。

$\displaystyle{ e^{-Ls}\simeq\frac{1-\frac{1}{2}Ls+\frac{1}{12}L^2s^2}{1+\frac{1}{2}Ls+\frac{1}{12}L^2s^2} \nonumber }$

以上の１入力１出力系の実現を利用して，多入力出力系の場合の実現の一つを求める方法を考える。

$m$ 入力 $p$ 出力系の伝達関数行列

$\displaystyle{(8.11)\quad G(s)= \left[\begin{array}{cccc} G_{11}(s) & \cdots & G_{1m}(s) \\ \vdots & \ddots & \vdots \\ G_{p1}(s) & \cdots & G_{pm}(s) \end{array}\right] }$

の各要素の実現を

$\displaystyle{(8.12)\quad G_{ij}(s)= \left[\begin{array}{c|c} A_{ij} & B_{ij} \\\hline C_{ij} & D_{ij} \end{array}\right] \ (i=1,\cdots,p;\,j=1,\cdots,m) }$

とするとき， $G(s)$ の実現は

$\displaystyle{(8.13)\quad \left[\begin{array}{cccc} G_{11}(s) & \cdots & G_{1m}(s) \\ \vdots & \ddots & \vdots \\ G_{p1}(s) & \cdots & G_{pm}(s) \end{array}\right] = \left[\begin{array}{ccc|ccc} A_1 & \cdots & 0 & B_{1} \\ \vdots & \ddots & \vdots & \vdots \\ 0 & \cdots & A_p & B_{p} \\\hline C_{1} & \cdots & C_{p} & D \end{array}\right] }$

で与えられる。ここで

$\displaystyle{(8.14)\quad A_i= \left[\begin{array}{ccc} A_{i1} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & A_{im} \end{array}\right] \ (i=1,\cdots,p) }$

は，サイズ $(n_{i1}+\cdots+n_{im})\times (n_{i1}+\cdots+n_{im})$ の行列

$\displaystyle{(8.15)\quad B_i= \left[\begin{array}{ccc} B_{i1} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & B_{im} \end{array}\right] \ (i=1,\cdots,p) }$

は，サイズ $(n_{i1}+\cdots+n_{im})\times m$ の行列

$\displaystyle{(8.16)\quad C_i= \left[\begin{array}{ccc} C_{i1} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & C_{im} \end{array}\right] \ (i=1,\cdots,p) }$

は，サイズ $p\times (n_{i1}+\cdots+n_{im})$ の行列

$\displaystyle{(8.17)\quad D= \left[\begin{array}{ccc} D_{11} & \cdots & D_{1m} \\ \vdots & \ddots & \vdots \\ D_{p1} & \cdots & D_{pm} \end{array}\right] \ (i=1,\cdots,p) }$

は，サイズ $p\times m$ の行列である。

ちなみに， $m$ 入力 $p$ 出力伝達関数行列

$\displaystyle{(8.18)\quad G(s) =\frac{1}{s^n+a_1s^{n-1}+\cdots+a_n} (G_1s^{n-1}+\cdots+G_n)+D }$

の可制御な実現の一つは

$\displaystyle{(8.19)\quad G(s)&=& \left[\begin{array}{cccc|c} -a_1I_m & \cdots & -a_{n-1}I_m & -a_nI_m & I_m \\ I_m & \cdots & 0 & 0 & 0 \\ \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & \cdots & I_m & 0 & 0 \\\hline G_1 & \cdots & G_{n-1} & G_n & D \end{array}\right] }$
$\displaystyle{(8.20)\quad &=& \left[\begin{array}{cccc|c} 0 & I_m & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & I_m & 0 \\ -a_nI_m & -a_{n-1}I_m & \cdots & -a_1I_m & I_m \\\hline G_n & G_{n-1} & \cdots & G_1 & D \end{array}\right] }$

また，可観測な実現の一つは

$\displaystyle{(8.21)\quad G(s)&=& \left[\begin{array}{cccc|c} -a_1I_p & I_p & \cdots & 0 & G_1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ -a_{n-1}I_p & 0 & \cdots & I_p & G_{n-1} \\ -a_nI_p & 0 & \cdots & 0 & G_n \\\hline I_p & 0 & \cdots & 0 & D \end{array}\right]\ }$
$\displaystyle{(8.22)\quad &=& \left[\begin{array}{cccc|c} 0 & \cdots & 0 & -a_nI_p & G_n \\ I_p & \cdots & 0 & -a_{n-1}I_p & G_{n-1} \\ \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & \cdots & I_p & -a_1I_p & G_1 \\\hline 0 & \cdots & 0 & I_p & D \end{array}\right]\ }$

で与えられる。

例題8.2 つぎの１入力２出力伝達関数行列の実現を一つ求めよ。

$\displaystyle{ G(s)= \left[\begin{array}{c} \frac{1}{s+1}\\ \frac{s+1}{s^2+3s+2} \end{array}\right] }$

解答

$\displaystyle{ \frac{1}{s+1} = \left[\begin{array}{c|c} -1 & 1 \\\hline 1 & 0 \end{array}\right] }$

$\displaystyle{ \frac{s+1}{s^2+3s+2} = \left[\begin{array}{cc|c} -3 & -2 & 1 \\ 1 & 0 & 0 \\\hline 1 & 1 & 0 \end{array}\right] }$

だから

$\displaystyle{ G(s)= \left[\begin{array}{ccc|c} -1 & 0 & 0 & 1 \\ 0 & -3 & -2 & 1\\ 0 & 1 & 0 & 0 \\\hline 1 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 \end{array}\right] }$

となる。

演習8.2 つぎの２入力１出力伝達関数行列の実現を一つ求めよ。

$\displaystyle{ G(s)= \left[\begin{array}{cc} \frac{s+2}{s^2+3s+2} & \frac{1}{s+2} \end{array}\right] }$

8.2 最小実現

前節で求めた伝達関数行列の実現は，最小次元数をもつとは限らない。ここでは，最小次元数をもつ実現，すなわち最小実現を求める方法を考える。

一般に，不可制御かつ不可観測な状態空間表現(8.1)は適当な座標変換により，つぎの正準構造をもつように変換できることが知られている。

$\displaystyle{(8.23)\quad \left[\begin{array}{c|c} TAT^{-1} & TB \\\hline CT^{-1} & D \end{array}\right] = \left[\begin{array}{cccc|c} A_1 & 0 & X_{13} & 0 & B_1 \\ X_{21} & A_2 & X_{23} & X_{24} & B_2 \\ 0 & 0 & A_3 & 0 & 0 \\ 0 & 0 & X_{43} & A_4 & 0 \\\hline C_1 & 0 & C_3 & 0 & 0 \end{array}\right] }$

ここで，正方行列 $A_1$ ， $A_2$ ， $A_3$ ， $A_4$ の次数は一意に定まり， $(A_1,B_1)$ は可制御対， $(A_1,C_1)$ は可観測対である。この正準構造のブロック線図を図8.1に示す。

図8.1　正準構造のブロック線図

さて，(8.23)を表す伝達関数行列を $G(s)$ とおくと，つぎが成り立つ。

$\displaystyle{(8.24)\quad \left[\begin{array}{c|c} A_1 & B_1 \\\hline C_1 & 0 \end{array}\right]=G(s) }$

これが最小実現であることが知られている。すなわち，最小実現は，可制御かつ可観測な部分系であり，入力から出力までの伝達特性を表している。

最小実現の計算は，まず，正準構造の可制御な部分系

$\displaystyle{(8.25)\quad \left[\begin{array}{cc|c} A_1 & 0 & B_1 \\ X_{21} & A_2 & B_2 \\\hline C_1 & 0 & 0 \end{array}\right]=G(s) }$

を求め，この可観測な部分系として求めるか，または，正準構造の可観測な部分系

$\displaystyle{(8.26)\quad \left[\begin{array}{cc|c} A_1 & X_{13} & B_1 \\ 0 & A_3 & 0 \\\hline C_1 & C_3 & 0 \end{array}\right]=G(s) }$

を求め，この可制御な部分系として求めればよい。したがって，可制御な部分系または可観測な部分系の計算が基礎となる。

一般に，つぎのように変換する直交行列 $T$ が存在する。

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

ここで， $B_1,\cdots,B_{k-1}$ は横長の形状となり，行フルランク $m_1 \ge \cdots \ge m_{k-1}$ をもつ。また， $B_k$ は横長で，行フルランク $m_k$ をもつか零行列のどちらかである。 $(A,B)$ は，前者の場合は可制御対，後者の場合は不可制御対である。

そのような直交行列 $T$ と $m_1,\cdots,m_k$ を求める階段化アルゴリズムを実行するMファイルの作成例をつぎに示す。

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

ここで，tol は，7行で有効階数 m(j) を決めるために適切に選ばれる零判定基準である。

同様に，つぎのように変換する直交行列 $T$ が存在する。

$\displaystyle{(8.28)\quad \left[\begin{array}{c|c} T^TAT & T^TB \\\hline CT & D \end{array}\right] = \left[\begin{array}{ccccc|c} A_1 & C_2 & \cdots & 0 & 0 & X \\ X & A_2 & \ddots & \vdots & \vdots & \vdots \\ \vdots & \vdots & \ddots & C_{k-1} & 0 & X \\ X & X & \cdots & A_{k-1} & C_k & X \\ X & X & \cdots & X & A_k　 & X \\\hline C_1 & 0 & \cdots & 0 & 0 & D \end{array}\right] }$

ここで， $C_1,\cdots,C_{k-1}$ は縦長の形状となり，列フルランク $p_1 \ge \cdots \ge p_{k-1}$ をもつ。また， $C_k$ は縦長で，列フルランク $p_k$ をもつか零行列のどちらかである。 $(A,C)$ は，前者の場合，可観測対，後者の場合，不可観測対である。

そのような直交行列を求める階段化アルゴリズムを実行するMファイルの作成例をつぎに示す。

%staircase2.m
function [T,p]=staircase2(A,C,tol)
[T,p]=staircase(A’,C’,tol)

例題8.3 例題8.2で得た実現の可制御かつ可観測な部分系を求め，最小実現を得よ。

解答 Mファイル

%min_realization.m
A=[-1 0 0;0 -3 -2; 0 1 0]; B=[1;1;0]; C=[1 0 0;0 1 1]; D=[0;0];
[T,m]=staircase(A,B,0.01);
AA=T’*A*T; BB=T’*B; CC=C*T; n=sum(m);
sys=ss(AA(1:n,1:n),BB(1:n:),CC(:1:n),D); tf(sys)

を実行して，

$\displaystyle{ T= \left[\begin{array}{ccc} -0.7071 & 0.5774 & -0.4082 \\ -0.7071 & -0.5774 & 0.4082 \\ 0 & 0.5774 & 0.8165 \end{array}\right], \ m_1=1,\ m_2=1,\ m_3=0 }$

を得る。このとき

$\displaystyle{ \left[\begin{array}{c|c} T^TAT & T^TB \\\hline CT & D \end{array}\right] = \left[\begin{array}{ccc|c} -2.0000 & -0.0000 & 1.7321 & -1.4142 \\ -1.2247 & -1.0000 & 2.1213 & 0 \\ 0 & 0 & -1.0000 & 0 \\\hline -0.7071 & 0.5774 & -0.4082 & 0 \\ -0.7071 & 0 & 1.2247 & 0 \\ \end{array}\right] }$

となる。この可制御部分を取り出すと

$\displaystyle{ G(s)= \left[\begin{array}{cc|c} -2.0000 & -0.0000 & -1.4142\\ -1.2247 & -1.0000 & 0\\\hline -0.7071 & 0.5774 & 0\\ -0.7071 & 0 & 0\\ \end{array}\right] }$

となり，しかもこれは可観測であるので最小実現である。

演習8.3 演習8.2で得た実現の可制御かつ可観測な部分系を求め，最小実現を得よ。

演習問題の解答

【演習8.1】
$\frac{1-\frac{1}{2}Ls+\frac{1}{12}L^2s^2}{1+\frac{1}{2}Ls+\frac{1}{12}L^2s^2} =\frac{-\frac{12}{L}s}{s^2+\frac{6}{L}Ls+\frac{12}{L^2}}+1$

から，実現の一つは

$\left[\begin{array}{cc|c} 0 & 1 & 0 \\ -\frac{12}{L^2} & -\frac{6}{L} & 1 \\\hline 0 & -\frac{12}{L} & 1 \end{array}\right]$

である。

【演習8.2】
$\frac{s+2}{s^2+3s+2} = \left[\begin{array}{cc|c} -3 & -2 &1 \\ 1 & 0 & 0 \\\hline 1 & 2 & 0 \end{array}\right]$ ， $\frac{1}{s+2} = \left[\begin{array}{c|c} -2 & 1 \\\hline 1 & 0 \end{array}\right]$

だから

$G(s)= \left[\begin{array}{ccc|cc} -3 & -2 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & -2 & 0 & 1\\\hline 1 & 2 & 1& 0 & 0 \end{array}\right]$

を得る。

【演習8.3】Mファイル

%min_realization2.m
A=[-3 -2 0;1 0 0; 0 0 -2]; B=[1 0;0 0;0 1]; C=[1 2 1]; D=[0 0];
[T,m]=staircase(A,B,0.01);
AA=T’*A*T; BB=T’*B; CC=C*T; n=sum(m);
sys=ss(AA(1:n,1:n),BB(1:n:),CC(:1:n),D); tf(sys)

を実行して，

$T= \left[\begin{array}{ccc} -0.4082 & -0.1826 & -0.8944 \\ -0.8165 & -0.3651 & 0.4472 \\ -0.4082 & 0.9129 & 0.0000 \\ \end{array}\right]$ ， $m_1=1,\, m_2=1,\, m_3=0$

を得る。このとき

$\left[\begin{array}{c|c} T^TAT & T^TB \\\hline CT & D \end{array}\right] = \left[\begin{array}{ccc|cc} -1.1667 & 0.3727 & -0.0000 & -0.4082 & -0.4082 \\ 0.3727 & -1.8333 & -0.0000 & -0.1826 & 0.9129 \\ -2.7386 & -1.2247 & -2.0000 & -0.8944 & 0.0000 \\\hline -2.4495 & 0 & -0.0000 & 0 & 0 \\ \end{array}\right]$

となる。この可観測部分を取り出すと

$G(s)= \left[\begin{array}{cc|cc} -1.1667 & 0.3727 & -0.4082 & -0.4082 \\ 0.3727 & -1.8333 & -0.1826 & 0.9129 \\\hline -2.4495 & 0 & 0 & 0 \\ \end{array}\right]$

となり，しかもこれは可制御であるので最小実現である。

7 非線形システムの線形化

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

【本章のねらい】
・振子を例に，MAXIMA $^*$ を用いて，非線形運動方程式を計算し，非線形状態方程式を求める。また，平衡状態のまわりで線形化して，線形状態方程式を求め，安定性と可制御性の判定を行う。
・振子の線形状態方程式を用いて，積分動作を加えた状態フィードバックのLQI設計を行い，これを非線形状態方程式に適用して，安定化シミュレーションを行う。

* http://ja.wikipedia.org/wiki/Maximaを参照。

7.1 非線形システムのモデリングの例

次の例題を考える。

例題7.1 図7.1の振子を考える。MAXIMAを用いて，非線形運動方程式を計算し，これから非線形状態方程式を求めよ。さらに，平衡状態のまわりで線形化して，線形状態方程式を求めよ。

図7.1

解答運動エネルギー $T=\frac{1}{2}J\dot{\theta}^2$ （ $J=\frac{1}{3}m\ell^2$ ），ポテンシャルエネルギー $V=mg\cos\theta$ から，ラグランジュ関数 $L=T-V$ を得て，これをラグランジュ方程式

$\displaystyle{\frac{d}{dt}\left(\frac{\partial L}{\partial \dot\theta}\right)-\frac{\partial L}{\partial\theta}=0}$

に代入する。この計算をMAXIMAを用いて行うためには，コマンド

/*pen*/
dth:’diff(th(t),t);
ddth:’diff(th(t),t,2);
x:ell*sin(th(t));
y:ell*cos(th(t));
J:(1/3)*m*ell^2;
T:(1/2)*m*(diff(x,t)^2+diff(y,t)^2)+(1/2)*J*dth^2;
V:m*g*ell*cos(th(t));
L:T-V;
LE:diff(diff(L,dth),t)-diff(L,th(t))=0,trigreduce;
sol:solve(LE,ddth);
f:matrix([dth],[rhs(sol[1])]);

を与えればよい $^*$ 。すなわち，非線形運動方程式は次式により得られる。

$\displaystyle{\ddot{\theta}=\frac{3g}{4\ell}\sin\theta}$

これより，非線形状態方程式

$\underbrace{ \left[\begin{array}{c} \dot{\theta}\\ \ddot{\theta} \end{array}\right] }_{\dot \xi} = \underbrace{ \left[\begin{array}{c} \dot{\theta}\\ \frac{3g}{4\ell}\sin\theta \end{array}\right] }_{f(\xi)}$

が求められる。平衡状態は

$f(\xi^*)=0 \Rightarrow \xi^* =\left[\begin{array}{c} \theta^* \\ 0 \end{array}\right]$ ( $\theta^*=0,\pi$ )

のように求められる。この平衡状態 $\xi^*$ のまわりで， $f(\xi)$ の右辺を１次近似して

$\underbrace{ \left[\begin{array}{c} f_1(\xi)\\ f_2(\xi) \end{array}\right] }_{f(\xi)} \simeq \underbrace{\left. \left[\begin{array}{cc} 0 & 1 \\ \frac{\partial f_2}{\partial\theta} & \frac{\partial f_2}{\partial\dot{\theta}} \end{array}\right]\right|_{{\theta=\theta^*}\atop{\dot{\theta}=0\ }} }_{A} \underbrace{ \left[\begin{array}{c} \theta-\theta^*\\ \dot{\theta} \end{array}\right] }_{x=\xi-\xi^*}$

を得る。この計算をMAXIMAを用いて行うためには，コマンド

A:matrix([diff(f[1,1],th(t)),diff(f[1,1],dth)],
[diff(f[2,1],th(t)),diff(f[2,1],dth)]);

を与えればよい。すなわち，線形状態方程式は，次式のように求められる。

$\underbrace{ \frac{d}{dt} \left[\begin{array}{c} \theta-\theta^*\\ \dot{\theta} \end{array}\right] }_{\dot x} = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ \frac{3g}{4\ell}\cos\theta^* & 0 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} \theta-\theta^*\\ \dot{\theta} \end{array}\right] }_{x}$

この線形状態方程式は，平衡状態近傍では，

$\sin\theta\simeq\cos\theta^*(\theta-\theta^*)$

と近似できるので，これを非線形運動方程式に代入し

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

を得て，これから直接求めることもできる。

演習7.1 例題7.1において， $\ell=0.25{\rm [m]},g=9.8{\rm [m/s^2]}$ のとき，つぎの初期値について，Simulinkを用いて，非線形状態方程式と安定な線形状態方程式による時間応答 $\theta(t)\,(0\le t\le10)$ を比較せよ。
(1) $\theta(0)=(3/180)\pi,\dot{\theta}(0)=0$
(2) $\theta(0)=(177/180)\pi,\dot{\theta}(0)=0$

例題7.2 図7.2の台車によって駆動される振子を考える。MAXIMAを用いて，非線形運動方程式を計算し，これから非線形状態方程式を求めよ。さらに，平衡状態のまわりで線形化して，線形状態方程式を求めよ。

図7.2

解答台車の運動エネルギーは $T_{cart}=\frac{1}{2}M\dot{r}^2$ ，振子の運動エネルギーは
$T_{pen}=\frac{1}{2}m(\dot{x}^2+\dot{y}^2)+\frac{1}{2}J\dot{\theta}^2$ （ $J=\frac{1}{3}m\ell^2$ ），
であり，振子のポテンシャルエネルギーは $U_{pen}=mg\cos\theta$ 。これらから，
ラグランジュ関数 $L=T_{cart}+T_{pen}-U_{pen}$ を得て，これをラグランジュ方程式

$\displaystyle{\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{r}}\right)-\frac{\partial L}{\partial r}=F}$

$\displaystyle{\frac{d}{dt}\left(\frac{\partial L}{\partial \dot\theta}\right)-\frac{\partial L}{\partial\theta}=0}$

に代入する。この計算をMAXIMAを用いて行うためには，コマンド

/*pend*/
dr:’diff(r(t),t); ddr:’diff(r(t),t,2);
dth:’diff(th(t),t); ddth:’diff(th(t),t,2);
T0:(1/2)*M*diff(r(t),t)^2;
x1:r(t)+ell*sin(th(t)); y1:ell*cos(th(t));
J1:(1/3)*m*ell^2;
T1:(1/2)*m*(diff(x1,t)^2+diff(y1,t)^2)+(1/2)*J1*dth^2;
V1:m*g*y1;
L:T0+T1-V1;
LE1:diff(diff(L,dr),t)-diff(L,r(t))=F,trigreduce,ratsimp;
LE2:diff(diff(L,dth),t)-diff(L,th(t))=0,trigreduce,ratsimp;
sol:solve([LE1,LE2],[ddr,ddth]);

を与えればよい。この結果から，非線形運動方程式

$\underbrace{ \left[\begin{array}{cc} M+m & m\ell\cos\theta \\ m\ell\cos\theta & \frac{4}{3}m\ell^2 \end{array}\right] }_{M(\theta)} \underbrace{ \left[\begin{array}{c} \ddot{r}\\ \ddot{\theta} \end{array}\right] }_{\ddot{\xi}_1} + \underbrace{ \left[\begin{array}{cc} 0 & -m\ell\dot{\theta}\sin\theta \\ 0 & 0 \end{array}\right] }_{C(\theta)} \underbrace{ \left[\begin{array}{c} \dot{r}\\ \dot{\theta} \end{array}\right] }_{\dot{\xi}_1}$ \\
\hspace{5mm}
$+ \underbrace{ \left[\begin{array}{c} 0 \\ -m\ell g\sin\theta \end{array}\right] }_{G(\theta)} = \underbrace{ \left[\begin{array}{c} F\\ 0 \end{array}\right] }_{\zeta}$

が得られる。これより，非線形状態方程式

$\underbrace{ \left[\begin{array}{c} \dot{\xi}_1\\ \ddot{\xi}_1 \end{array}\right] }_{\dot \xi} = \underbrace{ \left[\begin{array}{c} \dot{\xi}_1\\ -M(\theta)^{-1}(\zeta-C(\theta)\dot{\xi}_1-G(\theta)) \end{array}\right] }_{f(\xi,\zeta)}$

が求められる。平衡状態は， $f(\xi^*,\zeta^*)=0$ より， $\dot{\xi}_1=0$ ，したがって， $G(\theta)=\zeta$ を満足することから

$\xi^* =\left[\begin{array}{c} 0 \\ \theta^* \\ 0 \\ 0 \end{array}\right]$ ( $\theta^*=0,\pi$ )，
$\zeta^* =\left[\begin{array}{c} F^* \\ 0 \end{array}\right]$ ( $F^*=0$ )

のように求められる。この平衡状態近傍で， $f(\xi,\zeta)$ の右辺を１次近似して

$\underbrace{ \left[\begin{array}{c} \dot{r}\\ \dot{\theta}\\ f_3(r,\theta,\dot{r},\dot{\theta})\\ f_4(r,\theta,\dot{r},\dot{\theta})\\ \end{array}\right] }_{f(\xi,\zeta)} \simeq \underbrace{\left. \left[\begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \frac{\partial f_3}{\partial r} & \frac{\partial f_3}{\partial\theta} & \frac{\partial f_3}{\partial\dot{r}} & \frac{\partial f_3}{\partial\dot{\theta}} \\ \frac{\partial f_4}{\partial r} & \frac{\partial f_4}{\partial\theta} & \frac{\partial f_4}{\partial\dot{r}} & \frac{\partial f_4}{\partial\dot{\theta}} \end{array}\right] \right|_{{{{{r=0}\atop{\theta=\theta^*}}\atop{\dot{r}=0}}\atop{\dot{\theta}=0}}\atop{F^*=0}} }_{A} \underbrace{ \left[\begin{array}{c} r \\ \theta-\theta^* \\ \dot{r} \\ \dot{\theta} \end{array}\right] }_{x=\xi-\xi^*}$
$+ \underbrace{\left. \left[\begin{array}{c} 0 \\ 0 \\ \frac{\partial f_3}{\partial F} \\ \frac{\partial f_4}{\partial F} \end{array}\right] \right|_{{{{{r=0}\atop{\theta=\theta^*}}\atop{\dot{r}=0}}\atop{\dot{\theta}=0}}\atop{F^*=0}} }_{B} \underbrace{ F }_{u}$

を得る。この計算をMAXIMAを用いて行うためには，コマンド

f:matrix([dr],[dth],[rhs(sol[1][1])],[rhs(sol[1][2])]);
a31:diff(f[3,1],r(t)); a32:diff(f[3,1],th(t));
a33:diff(f[3,1],dr); a34:diff(f[3,1],dth);
a41:diff(f[4,1],r(t)); a42:diff(f[4,1],th(t));
a43:diff(f[4,1],dr); a44:diff(f[4,1],dth);
b3:diff(f[3,1],F); b4:diff(f[4,1],F);
A:matrix([0,0,1,0],[0,0,0,1],[a31,a32,a33,a34],
[a41,a42,a43,a44]); B:matrix([0],[0],[b3],[b4]);
A1:A,th(t)=0,F=0,trigreduce,ratsimp;
B1:B,th(t)=0,F=0,trigreduce,ratsimp;
A2:A,th(t)=%pi,F=0,trigreduce,ratsimp;
B2:B,th(t)=%pi,F=0,trigreduce,ratsimp;

を与えればよい。この結果から， $\theta^*=0$ のときの線形状態方程式は

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

となる。また， $\theta^*=\pi$ のときの線形状態方程式は，{\tt th(t)=\%pi}
% \hspace{5mm}{\tt A2:A,th(t)=\%pi,F=0,trigreduce,ratsimp; }\\
% \hspace{5mm}{\tt B2:B,th(t)=\%pi,F=0,trigreduce,ratsimp;}\\
として

$\underbrace{ \frac{d}{dt} \left[\begin{array}{c} r \\ \theta-\pi \\ \dot{r} \\ \dot{\theta} \end{array}\right] }_{\dot x} = \underbrace{ \left[\begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & -\frac{3mg}{4M+m} & 0 & 0 \\ 0 & -\frac{3(M+m)g}{(4M+m)\ell} & 0 & 0 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} r \\ \theta-\pi \\ \dot{r} \\ \dot{\theta} \end{array}\right] }_{x} + \underbrace{ \left[\begin{array}{cccc} 0 \\ 0 \\ \frac{4}{4M+m} \\ \frac{3}{(4M+m)\ell} \\ \end{array}\right] }_{B} \underbrace{ F }_{u}$ \\
のように求められる。これらは，平衡状態近傍での近似式
$\sin\theta\simeq\cos\theta^*(\theta-\theta^*),\dot{\theta}^2=0$
を非線形運動方程式に代入して得られる

$\left[\begin{array}{cc} M+m & m\ell\cos\theta^* \\ m\ell\cos\theta^* & \frac{4}{3}m\ell^2 \end{array}\right] \left[\begin{array}{c} \ddot{r}\\ \ddot{\theta} \end{array}\right] + \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \left[\begin{array}{c} \dot{r}\\ \dot{\theta} \end{array}\right]$
$+ \left[\begin{array}{c} 0 \\ -m\ell g\cos\theta^*(\theta-\theta^*) \end{array}\right] = \left[\begin{array}{c} F\\ 0 \end{array}\right]$ \\
から直接求めたものと一致する。

演習7.2 例題7.2の振子を軸支した台車を，図7.3のように傾斜角 $\alpha$ の台上に置いた。MAXIMAを用いて，非線形運動方程式を計算し，これから非線形状態方程式を求めよ。さらに，平衡状態のまわりで線形化して，線形状態方程式を求めよ。

図7.3

例題7.3 例題7.2の台車によって駆動される振子について，そこで求めた状態方程式に基づいて，MAXIMAを用いて，漸近安定性と可制御性について調べよ。

解答まず，安定性を調べるために， $A$ 行列の固有値を計算する。これは，MAXIMAに，のコマンドを与えればよい。

lambda:eigenvalues(A1);

その結果， $\theta^*=0$ の場合， $\pm \sqrt{\frac{3(M+m)g}{(4M+m)\ell}}$ となり，また， $\theta^*=\pi$ の場合， $\pm j\sqrt{\frac{3(M+m)g}{(4M+m)\ell}}$ となる。どちらとも，漸近安定とはいえない。

つぎに，可制御性を， $A$ 行列の固有値に基づいて調べるためには，各固有値 $\lambda_i\,(i=1,\cdots,4)$ について，行列 $\left[\begin{array}{cc} B & \lambda I_n -A \end{array}\right]$ の階数が次数４に等しいことを確認すればよい。これは， $\theta^*=0$ の場合，MAXIMAに

rank(addcol(A1-lambda[1][1]*ident(4),B1));
rank(addcol(A1-lambda[1][2]*ident(4),B1));
rank(addcol(A1-lambda[2][1]*ident(4),B1));
rank(addcol(A1-lambda[2][2]*ident(4),B1));

% For[i = 1, i $\le$ 4,
% W=MapThread[Append,{A-lambda[[i]] IdentityMatrix[4],B}];
% Print[MatrixRank[W]]; i++]
のコマンドを与えればよい。その結果，階数はすべて４と計算され，可制御性が成り立つ。 $\theta^*=\pi$ の場合も同様である。

演習7.3 例題7.3を， $M=1{\rm [kg]},m=0.1{\rm [kg]},\ell=0.25{\rm [m]},g=9.8{\rm [m/s^2]}$ のとき，MATLABにより検討せよ。

7.2 非線形系に対する線形制御の適用

つぎの例題を考える。

例題7.4 例題7.2の台車によって駆動される倒立振子について，台車の位置 $r$ と振子の角度 $\theta$ が計測できるものとする。このとき，振子を倒立状態に保ちながら，台車を指定位置 $r^*=0.5{\rm [m]}$ まで移動させるための，積分動作を加えた状態フィードバックを求めよ。ただし， $M=1{\rm [kg]},m=0.1{\rm [kg]},\ell=0.25{\rm [m]},g=9.8{\rm [m/s^2]}$ とする。

解答まず，出力方程式として

$\underbrace{ \left[\begin{array}{c} r \\ \theta \end{array}\right] }_{y_M} = \underbrace{ \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{array}\right] }_{C_M} \underbrace{ \left[\begin{array}{c} r \\ \theta \\ \dot{r} \\ \dot{\theta} \end{array}\right] }_{x}$

を得る。また，台車の位置 $r$ は

$\underbrace{ r }_{z} = \underbrace{ \left[\begin{array}{cccc} 1 & 0 \end{array}\right] }_{C_S} \underbrace{ \left[\begin{array}{c} r \\ \theta \end{array}\right] }_{y_M} = \underbrace{ \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \end{array}\right] }_{C=C_SC_M} \underbrace{ \left[\begin{array}{c} r \\ \theta \\ \dot{r} \\ \dot{\theta} \end{array}\right] }_{x}$

のように表される。このとき，正方行列

$S=\left[\begin{array}{cc} A & B \\ C & 0 \end{array}\right] = \left[\begin{array}{cccc:c} 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0\\ 0 & -\frac{3mg}{4M+m} & 0 & 0 & \frac{3}{4M+m}\\ 0 & \frac{3(M+m)g}{(4M+m)\ell} & 0 & 0 & -\frac{3}{(4M+m)\ell} \\ \hdashline 1 & 0 & 0 & 0 & 0 \end{array}\right]$

は正則である。したがって，次式が定まる。

$\left[\begin{array}{c} x_\infty \\ u_\infty \end{array}\right] = \left[\begin{array}{c} r_\infty \\ \theta_\infty \\ \dot{r}_\infty \\ \dot{\theta}_\infty \\\hdashline u_\infty \end{array}\right] =S^{-1} \left[\begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\\hdashline r^* \end{array}\right] = \left[\begin{array}{c} r^* \\ 0 \\ 0 \\ 0 \\\hdashline 0 \end{array}\right]$

つぎに，偏差系

$\frac{d}{dt} \underbrace{ \left[\begin{array}{c} x-x_\infty \\ u-u_\infty \end{array}\right] }_{x_{E3}} = \underbrace{ \left[\begin{array}{cccc} A & B \\ 0 & 0 \end{array}\right] }_{A_{E3}} \underbrace{ \left[\begin{array}{c} x-x_\infty \\ u-u_\infty \end{array}\right] }_{x_{E3}} + \underbrace{ \left[\begin{array}{cccc} 0 \\ 1 \\ \end{array}\right] }_{B_{E3}} \dot{u}$

に対する評価関数を，たとえば

$\displaystyle{J=\int_0^\infty\{\frac{1}{K_r^2}(r-r^*)^2+\frac{1}{K_\theta^2}\theta^2+\frac{T_r^2}{K_r^2}\dot{r}^2+\frac{T_\theta^2}{K_\theta^2}\theta^2}$ $\displaystyle{+\frac{1}{K_u^2}u^2+\frac{T_u^2}{K_u^2}\dot{u}^2\}\,dt}$

のように設定する。ここで， $T_r,K_r$ と $T_\theta,K_\theta$ は，それぞれ閉ループ系における台車位置と角度の時間応答の速さ（どれくらいの時間である値に到達するか）を表している。たとえば，台車を5秒間に50cm移動させるとき，周期 $T_\ell$ の振子は3度以内の振れに抑えることが妥当と考えられるときは， $T_r=5,K_r=0.5$ と $T_\theta=(1/4)T_\ell,K_\theta=(5/180)\pi$ と与える。つぎに， $T_u,K_u$ は駆動系の時間応答の速さを考慮して選ぶ。

制御目的を達成する積分動作を加えた状態フィードバック

$u=-Fx-F_Ix_I$ 　ただし， $x_I=\int_0^t(r(\tau)-r^*)\,d\tau$

は，上の評価関数を最小にするように，状態フィードバック

$\dot{u}=- \left[\begin{array}{cc} K & K_I \end{array}\right]x_{E3}$

を求め，次式から $F$ と $F_I$ を定めればよい。

$\left[\begin{array}{cc} F & F_I \end{array}\right] = \left[\begin{array}{cc} K & K_I \end{array}\right]S^{-1}$

以上の計算を，MATLABで行うためには

%pend.m
global M m ell g th0
M=1; m=0.1; ell=0.25; g=9.8; th0=3/180*pi;
E=[M+m m*ell;m*ell (4/3)*m*ell^2];
A21=-E\[0 0;0 -m*ell*g]; A=[zeros(2) eye(2);A21 zeros(2)]
B2=E\[1;0]; B=[zeros(2,1);B2]
CM=[1 0 0 0;0 1 0 0]; CS=[1 0]; C=CS*CM;
AE=[A B;zeros(1,5)]; BE=[zeros(4,1); 1]; CE=eye(5);
Tcart=0.1; Kr=0.5;
Tpend=(1/4)*(2*pi*sqrt(4/3*ell/g)); Kth=5/180*pi;
Tamp=0.1; Kamp=5;
q1=1/Kr; q2=1/Kth; q3=Tcart/Kr; q4=Tpend/Kth; q5=1/Kamp;
QE=diag([q1 q2 q3 q4 q5].^2); RE=(Tamp/Kamp)^2;
[FE,pAE]=opt(AE,BE,CE,QE,RE)
FE=FE/[A B;C 0]; F=FE(:,1:4); FI=FE(:,5);
ACL=[A-B*F -B*FI;C 0];
BCL=[zeros(4,1);-0.5];
CCL=[CM zeros(2,1);-F -FI];
DCL=[zeros(2,1);0];
sys=ss(ACL,BCL,CCL,DCL);
step(sys)

を与えればよい。

もちろん，これらは重み係数の初期値であり，実際には調整を伴うことはいうまでもない。

演習7.4 演習7.3のように，振子を軸支した台車を傾斜角 $\alpha=(5/180)\pi{\rm [rad]}$ の台上に置く場合，例題7.4と同様の設計を行え。

例題7.5 例題7.4で設計した積分動作を加えた状態フィードバックを，非線形ダイナミクスを表す非線形状態方程式と結合して，閉ループ系の時間応答をシミュレーションせよ。

解答まず，非線形状態方程式を，つぎのS-functionで記述する。

%spend.m
function [sys,x0]=spend(t,state,input,flag)
global M m ell g th0
if abs(flag)==1
u =input(1);
r =state(1); th =state(2); dr =state(3); dth=state(4);
Mp=[M+m m*ell*cos(th); m*ell*cos(th) (4/3)*m*ell^2];
Cp=[0 -m*ell*dth*sin(th);0 0];
Gp=[0; -m*ell*g*sin(th)];
sys=[dr;dth;Mp\(-Cp*[dr;dth]-Gp+[u;0])];
elseif flag==3,
sys=state(1:2);
elseif flag==0
sys=[4;0;2;1;0;0];
x0=[0;th0;0;0];
else
sys=[];
end

つぎに，Simulink上に

図7.4

のように，このS-functionブロックと例題??で設計した積分動作を加えた状態フィードバックを接続して，シミュレーションを行う。

演習7.5 演習7.4で設計した積分動作を加えた状態フィードバックを，非線形ダイナミクスを表す非線形状態方程式と結合して，閉ループ系の時間応答をシミュレーションせよ。

演習問題の解答

【演習7.1】Simulinkを用いて，次図のブロック線図を作成する。

上半分が非線形運動方程式 $\ddot{\theta}=\frac{3g}{4\ell}\sin\theta$ を，下半分が平衡状態 $\theta^*=\pi$ まわりの線形運動方程式 $\ddot{\theta}=-\frac{3g}{4\ell}(\theta-\theta^*)$ を表している。角度を得る積分器の初期値は，非線形運動方程式の場合は $\theta(0)$ を，線形運動方程式の場合は $\theta(0)-\pi$ を与えることに注意する。すなわち，(1)では，非線形運動方程式の場合は $\theta(0)=(3/180)\pi$ を，線形運動方程式の場合は $\theta(0)-\pi=(-177/180)\pi$ を与える。また，(2)では，非線形運動方程式の場合は $\theta(0)=(177/180)\pi$ を，線形運動方程式の場合は $\theta(0)-\pi=(-3/180)\pi$ を与える。

【演習7.2】非線形運動方程式を求める計算をMAXIMAを用いて行うためには

/*pend2*/
dr:’diff(r(t),t); ddr:’diff(r(t),t,2);
dth:’diff(th(t),t); ddth:’diff(th(t),t,2);
x0:r(t)*cos(alpha); y0:r(t)*sin(alpha);
T0:(1/2)*M*(diff(x0,t)^2+diff(y0,t)^2);
V0:M*g*y0;
x1:x0+ell*sin(th(t)); y1:y0+ell*cos(th(t));
J1:(1/3)*m*ell^2;
T1:(1/2)*m*(diff(x1,t)^2+diff(y1,t)^2)+(1/2)*J1*dth^2;
V1:m*g*y1;
L:T0+T1-V0-V1;
LE1:diff(diff(L,dr),t)-diff(L,r(t))=F,trigreduce,ratsimp;
LE2:diff(diff(L,dth),t)-diff(L,th(t))=0,trigreduce,ratsimp;
sol:solve([LE1,LE2],[ddr,ddth]);

のコマンドを与えればよい。この結果から，非線形運動方程式
$\underbrace{ \left[\begin{array}{cc} M+m & m\ell\cos(\theta+\alpha) \\ m\ell\cos(\theta+\alpha) & \frac{4}{3}m\ell^2 \end{array}\right] }_{M(\theta)} \underbrace{ \left[\begin{array}{c} \ddot{r}\\ \ddot{\theta} \end{array}\right] }_{\ddot{\xi}_1} + \underbrace{ \left[\begin{array}{cc} 0 & -m\ell\dot{\theta}\sin(\theta+\alpha) \\ 0 & 0 \end{array}\right] }_{C(\theta)} \underbrace{ \left[\begin{array}{c} \dot{r}\\ \dot{\theta} \end{array}\right] }_{\dot{\xi}_1}$
$+ \underbrace{ \left[\begin{array}{c} 0 \\ -m\ell g\sin\theta \end{array}\right] }_{G(\theta)} = \underbrace{ \left[\begin{array}{c} F-(M+m)g\sin\alpha\\ 0 \end{array}\right] }_{\zeta}$
が得られる。これより，平衡状態は
$\xi^* =\left[\begin{array}{c} 0 \\ \theta^* \\ 0 \\ 0 \end{array}\right]$ ( $\theta^*=0,\pi$ )，
$\zeta^* =\left[\begin{array}{c} F^* \\ 0 \end{array}\right]$ ( $F^*=(M+m)g\sin\alpha$ )
のように求まる。 $\theta^*=0$ のときの線形状態方程式は
$\underbrace{ \frac{d}{dt} \left[\begin{array}{c} r \\ \theta \\ \dot{r} \\ \dot{\theta} \end{array}\right] }_{\dot x} = \underbrace{ \left[\begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & -\frac{6mg\cos\alpha}{8M+(5-3\cos2\alpha)m} & 0 & 0 \\ 0 & \frac{6(M+m)g}{(8M+(5-3\cos2\alpha)m)\ell} & 0 & 0 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} r \\ \theta \\ \dot{r} \\ \dot{\theta} \end{array}\right] }_{x}$ \\
$+ \underbrace{ \left[\begin{array}{cccc} 0 \\ 0 \\ \frac{8}{8M+(5-3\cos2\alpha)m} \\ -\frac{6\cos\alpha}{(8M+(5-3\cos2\alpha)m)\ell} \\ \end{array}\right] }_{B} \underbrace{ (F-(M+m)g\sin\alpha) }_{u}$

【演習7.3】

%pend.m
M=1; m=0.1; ell=0.25; g=9.8;
E=[M+m m*ell;m*ell (4/3)*m*ell^2];
A21=-E\[0 0;0 -m*ell*g]; A=[zeros(2) eye(2);A21 zeros]
B2=E\[1;0]; B=[zeros(2,1);B2]
r=eig(A)
tol=1e-6;
for i=1:4, k=rank([B A-r(i)*eye(4)],tol), end

【演習7.4】次を実行し，あとは例題7.4と同様に行えばよい。

%pend2.m
M=1; m=0.1; ell=0.25; g=9.8; alpha=5/180*pi;
a32=-(3*m*g*cos(alpha))/(7*M+(11-3*cos(2*alpha))/2*m);
a42=(3*(M+m)*g)/(7*M+(11-3*cos(2*alpha))/2*m)/ell;
b3=7/(7*M+(11-3*cos(2*alpha))/2*m);
b4=-3/(7*M+(11-3*cos(2*alpha))/2*m)/ell;
A21=[0 a32;0 a42]; A=[zeros(2) eye(2);A21 zeros(2)]
B2=[b3;b4]; B=[zeros(2,1);B2]

【演習7.5】まず，非線形ダイナミクスを，つぎのS-functionで記述する。

%spend2.m
function [sys,x0]=spend2(t,state,input,flag)
M=1; m=0.1; ell=0.25; g=9.8; th0=0; alpha=5/180*pi;
if abs(flag)==1
u=input(1);
r=state(1); th=state(2); dr=state(3); dth=state(4);
Mp=[M+m m*ell*cos(th+alpha); m*ell*cos(th+alpha) (7/3)*m*ell^2];
Cp=[0 -m*ell*dth*sin(th+alpha);0 0];
Gp=[0; -m*ell*g*sin(th)];
sys=[dr;dth;Mp $\backslash$ (-Cp*[dr;dth]-Gp+[u-(M+m)*g*sin(alpha);0])];
elseif flag==3
sys=state(1:2);
elseif flag==0
sys=[4;0;2;1;0;0]; x0=[0;th0;0;0];
else
sys=[];
end

Simulink上に，このS-functionブロックと，演習7.4で設計した積分動作を加えた状態フィードバックを記述して，シミュレーションを行う。詳細は割愛する。

6. 定値外乱を抑制する

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

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

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) }$

以下では，この状況を一般化した次式で表される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.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 }$

図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) }$

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

を合わせて

を構成する。式(16)に

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

となる。そのブロック線図を図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]$ だから

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

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

図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 外乱を抑制したうえで設定値を変更する

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

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

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

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

を代入すると

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

を得る。この第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) }$

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

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

の関係を満足しているはずである。したがって，どのような $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)を合わせて

を得る。式(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)} }$

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

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

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

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

を用いて

ここで， $\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{(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{(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)を合わせて

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

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

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

を行えば

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

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

すなわち

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

6.3 LQI設計

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

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

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

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

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

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

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

問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] }$

以下に，偏差系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　被制御変数の決定

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

偏差系

を，状態フィードバック

$\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制御で安定化する。その際，評価関数としては

を用いる。ただし， $(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] }$

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