制御系CAD

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

を行うと，閉ループ系

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

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

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

となる。第2項は

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

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

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

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

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

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

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

図6.2 定値外乱の影響

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

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

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

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

の時間応答は

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

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

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

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

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

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

6.1.2 積分動作を加える

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

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

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

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

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

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

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

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

を合わせて

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

を構成する。式(16)に

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

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

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

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

図6.3 積分動作を入れる

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

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

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

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

と計算できるから

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

と表される。したがって

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

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

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

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

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

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

図6.4 定値外乱の抑制

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

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

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

6.1.3 設定値を変更する

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

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

に対して，制御目的

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

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

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

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

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

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

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

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

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

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

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

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

を代入すると

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

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

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

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

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

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

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

また，つぎが成り立つ。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

を用いて

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

の出力を用いて

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

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

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

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

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

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

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

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

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

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

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

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

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

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

を行えば

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

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

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

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

すなわち

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

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

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

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

6.3 LQI設計

5 LQG制御

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

偏差系

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

を，状態フィードバック

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

ただし

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

に対して，評価関数

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

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

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

の解 $\Pi$ を求めると

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

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

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

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

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

例題5.2　２次系

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

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

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

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

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

を要素ごとに整理して

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

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

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

ここで，(1)だけが， $\pi_1,\pi_2>0,\ \pi_1\pi_2-\pi_3^2>0$ を満たす
（ $\left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]>0\Leftrightarrow \pi_1>0,\ \pi_1\pi_2-\pi_3^2>0$ 。このとき， $\pi_2>0$ 。）。したがって

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

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

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

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

例題5.3　２次系

$\displaystyle{ \left\{\begin{array}{l} %\underbrace{ \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right] %}_{\dot{x}(t)} = %\underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 0 & -2\zeta\omega_n \end{array}\right] %}_{A} %\underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] + %\underbrace{ \left[\begin{array}{c} 0 \\ \omega_n^2 \end{array}\right] %}_{B} \\ y(t)= %\underbrace{ \left[\begin{array}{cc} 1 & 0 \end{array}\right] %}_{C} %\underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] \end{array}\right. }$

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

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

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

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

となることを示せ。
解答　リッカチ方程式

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

を要素ごとに整理して

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

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

$\displaystyle{ \pi_3=\pm r\omega_n^{-2}q }$

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

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

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

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

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

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

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

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

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

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

演習5.3　２次系

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

によって得る。

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

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

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

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

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

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

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

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

eig(A-B*F)

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

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

演習5.5　３次系

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

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

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

演習5.6　４次系

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

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

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

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

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

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

ここで，新しい入力 $w$ と $v$ がそれぞれ $W>0$ と $V>0$ の平方根行列（ $X>0(\ge0)$ に対し $YY=X$ を満足する行列 $Y>0(\ge0)$ を $X^{\frac{1}{2}}$ で表す。）
により重み付けられて， $n$ 次系の入力側（ $B'$ を介して）と出力側に設置されている。また新しい出力 $z=C'x$ と入力 $u$ が取り出されており，
それぞれ $Q>0$ と $R>0$ の平方根行列により重み付けられている。

この $n$ 次系

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

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

$\displaystyle{(5.13)\quad \left\{\begin{array}{l} \dot{\hat{x}}(t)=(A-HC-BF)\hat{x}(t)+Hy(t) \\ u(t)=-F\hat{x}(t) \end{array}\right. }$

による閉ループ系は，次式で表される。

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

このインパルス応答行列

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

に関する評価規範として

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

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

Step 1　行列方程式

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

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

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

Step 2　行列方程式

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

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

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

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

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

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

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

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

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

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

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

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

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

演習問題の解答

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

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

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

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

$\begin{array}{ll} \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ 0 & -1 \end{array}\right] + \left[\begin{array}{cc} 0 & 0 \\ 1 & -1 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]- \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \\ \times \left[\begin{array}{c} 0 \\ 1 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] + \left[\begin{array}{c} 1 \\ 0 \end{array}\right] \left[\begin{array}{cc} 1 & 0 \end{array}\right] = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \end{array}$

を要素ごとに整理して

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

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

$F= \left[\begin{array}{cc} 0 & 1 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]= \left[\begin{array}{cc} \pi_3 & \pi_2 \end{array}\right]= \left[\begin{array}{cc} 1 & -1+\sqrt{3} \end{array}\right]$
%[f,p]=opt([0 1;0 -1],[0;1],[1 0],1,1)

(2)リッカチ方程式

$\begin{array}{ll} \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ -1 & 0 \end{array}\right] + \left[\begin{array}{cc} 0 & -1 \\ 1 & 0 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]- \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \\ \times \left[\begin{array}{c} 0 \\ 1 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] + \left[\begin{array}{c} 1 \\ 0 \end{array}\right] \left[\begin{array}{cc} 1 & 0 \end{array}\right] = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \end{array}$

を要素ごとに整理して

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

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

$F= \left[\begin{array}{cc} 0 & 1 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]= \left[\begin{array}{cc} \pi_3 & \pi_2 \end{array}\right]= \left[\begin{array}{cc} -1+\sqrt{2} & \sqrt{-2+\sqrt{2}} \end{array}\right]$
%[f,p]=opt([0 1;-1 0],[0;1],[1 0],1,1)

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

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

を要素ごとに整理して

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

を得る。 $\pi_1>0,\ \pi_1\pi_2-\pi_3^2>0$ を満たすものを選ぶと

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

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

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

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

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

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

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

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

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

【演習5.8】
%lqr5.m
[A

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

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

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

4.1 状態オブザーバ

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

例題4.1　２次系

$\displaystyle{ \left\{\begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \\ \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \end{array}\right] }_{x(t)}+ \underbrace{ \left[\begin{array}{cc} 0 \\ 1 \end{array}\right] }_{B} u(t)\\ y(t)= \underbrace{ \left[\begin{array}{cc} 1 & 0 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \end{array}\right] }_{x(t)} \end{array}\right. }$

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

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

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

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

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

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

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

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

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

演習4.1　２次系

$\displaystyle{ \left\{\begin{array}{l} \underbrace{ \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \\ \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \end{array}\right] }_{x(t)}+ \underbrace{ \left[\begin{array}{cc} 0 \\ 1 \end{array}\right] }_{B} u(t)\\ y(t)= \underbrace{ \left[\begin{array}{cc} 1 & 1 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \end{array}\right] }_{x(t)} \end{array}\right. }$

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

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

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

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

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

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

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

$\displaystyle{(4.10)\quad \left\{\begin{array}{l} \dot{\hat{x}}(t)=(A-HC-BF)\hat{x}(t)+Hy(t) \\ u(t)=-F\hat{x}(t) \end{array}\right. }$

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

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

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

例題4.3　１次系

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

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

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

と，状態オブザーバ

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

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

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

による閉ループ系の $A$ 行列の固有値を求めよ。
解答オブザーバベースト・コントローラによる閉ループ系は

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

すなわち

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

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

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

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

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

4.2 可観測性と可検出性

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

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

定義D0:　状態オブザーバを構成可能
条件D1:　 ${\rm rank}\, \left[\begin{array}{c} C \\ A-\lambda I_n \end{array}\right]=n$ 　( $\lambda$ は $A$ のすべての不安定固有値)
条件D2:　 $Cv=0,\ Av=\lambda v \Rightarrow v=0$ 　( $\lambda$ は $A$ のすべての不安定固有値)

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

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

定義O0:　任意有限時間の入力と出力から，初期状態を一意に決定可能
条件O1: $\displaystyle{\int_0^t \exp(A^T\tau)C^TC\exp(A\tau)\,d\tau>0 \quad (\forall t>0)}$
条件O2:　 ${\rm rank}\, \underbrace{ \left[\begin{array}{c} C \\ CA \\ \vdots\\ CA^{n-1} \end{array}\right] }_{observability\ matrix} =n$
条件O3:　 $H$ を選ぶことにより， $A-HC$ の固有値を任意に設定可能
条件O4:　 ${\rm rank} \left[\begin{array}{c} C \\ A-\lambda I_n \end{array}\right]=n$ 　( $\lambda$ は $A$ のすべての固有値)
条件O5:　 $Cv=0,\ Av=\lambda v \Rightarrow v=0$ 　( $\lambda$ は $A$ のすべての固有値)
条件O6:　 $(A^T,C^T)$ は可制御対

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

が成り立ち，これから

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

したがって

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

演習問題の解答

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

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

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

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

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

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

【演習4.2】

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

【演習4.3】

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 状態フィードバック

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$\displaystyle{(3.8)\quad \begin{array}{lll} F&=&\left[\begin{array}{cccc} a_n'-a_n & \cdots & a_2'-a_2 & a_1'-a_1 \end{array}\right] \nonumber\\ &\times& \left[\begin{array}{ccccc} a_{n-1} & a_{n-2} & \cdots & a_1 & 1 \\ a_{n-2} & a_{n-3} & \cdots & 1 & 0 \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ a_2 & a_1 & \cdots & 0 & 0 \\ a_1 & 1 & \cdots & 0 & 0 \\ 1 & 0 & \cdots & 0 & 0 \end{array}\right]^{-1} \nonumber\\ &\times& \left[\begin{array}{cccc} B & AB &\cdots & A^{n-1}B \end{array}\right]^{-1} \end{array} }$

または

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

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

例題3.3　２次系

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \\ \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \end{array}\right] }_{{x}(t)}+ \underbrace{ \left[\begin{array}{c} 0 \\ 1 \end{array}\right] }_{B} u(t) }$

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

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

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

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

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

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

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

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

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

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

演習3.4　つぎの２次系に対する状態フィードバック $u=-Fx$ を，閉ループ系の $A$ 行列の固有値が $\lambda'_1=\lambda'_2=-1$ となるように求めよ。
(1) $\underbrace{ \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \\ \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 0 & -1 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \end{array}\right] }_{{x}(t)}+ \underbrace{ \left[\begin{array}{c} 0 \\ 1 \end{array}\right] }_{B} u(t)$

(2) $\underbrace{ \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \\ \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -1 & 0 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \end{array}\right] }_{{x}(t)}+ \underbrace{ \left[\begin{array}{c} 0 \\ 1 \end{array}\right] }_{B} u(t)$

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

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

ここで

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

とおくと

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

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

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

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

例題3.4　２入力２次系

$\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \\ \end{array}\right] }_{\dot{x}(t)}= \underbrace{ \left[\begin{array}{cc} 0 & 0 \\ 0 & -1 \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \end{array}\right] }_{{x}(t)}+ \underbrace{ \left[\begin{array}{cc} 1 & 1 \\ 1 & -1 \end{array}\right] }_{B} u(t) }$

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

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

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

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

この特性多項式は

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

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

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

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

この特性多項式は

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

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

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

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

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

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

3.2 可制御性と可安定性

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

【可安定性の定義とその等価な条件】
定義S0:　状態フィードバックにより安定化可能
条件S1:　 ${\rm rank}\, \left[\begin{array}{cc} B & A-\lambda I_n \end{array}\right] =n$ ( $\lambda$ は $A$ のすべての不安定固有値)
条件S2:　 $B^Tw=0,\ A^Tw=\lambda w \Rightarrow w=0$ ( $\lambda$ は $A$ のすべての不安定固有値)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

演習問題の解答

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

演習3.6

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

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

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

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

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

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

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

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

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

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

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

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

2 漸近安定性

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

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

2.1 漸近安定性

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

で定義される。たとえば

$\displaystyle{(2.7)\quad \exp \left(\left[\begin{array}{cc} \lambda_1 & 0\\ 0 & \lambda_2 \end{array}\right]t\right) = \left[\begin{array}{cc} e^{\lambda_1t} & 0 \\ 0 & e^{\lambda_2t} \end{array}\right] }$
$\displaystyle{(2.8)\quad \exp \left(\left[\begin{array}{cc} \lambda_R & \lambda_I \\ \lambda_I & \lambda_R \end{array}\right]t\right) =e^{\lambda_R t} \left[\begin{array}{cc} \cos \lambda_I t & \sin \lambda_I t \\ -\sin \lambda_I t & \cos \lambda_I t \end{array}\right] }$
$\displaystyle{(2.9)\quad \exp \left(\left[\begin{array}{cc} \lambda & 1 \\ 0 & \lambda \end{array}\right]t\right) =e^{\lambda t} \left[\begin{array}{cc} 1 & t \\ 0 & 1 \end{array}\right] }$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

解答

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 時間応答

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

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

ただし，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

は，次式で与えられる。

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

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

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

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

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

ただし

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

演習問題の解答

演習2.1

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

演習2.4

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

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

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

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

演習2.6

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

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

演習2.7

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

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

演習2.8

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

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

演習2.9

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

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