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

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

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

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

次の例題を考える。

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


図7.1

解答 運動エネルギー T=\frac{1}{2}J\dot{\theta}^2J=\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)の右辺を1次近似して

\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}^2J=\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)の右辺を1次近似して

\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]の階数が次数4に等しいことを確認すればよい。これは,\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++]
のコマンドを与えればよい。その結果,階数はすべて4と計算され,可制御性が成り立つ。\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_rT_\theta,K_\thetaは,それぞれ閉ループ系における台車位置と角度の時間応答の速さ(どれくらいの時間である値に到達するか)を表している。たとえば,台車を5秒間に50cm移動させるとき,周期T_\ellの振子は3度以内の振れに抑えることが妥当と考えられるときは,T_r=5,K_r=0.5T_\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}

を求め,次式からFF_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. 定値外乱を抑制する

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

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

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

したがって,これまでの知識だけでは多くの現実問題を解くことはできない。実際の問題では,必ずといっていいほど,さまざまな「定値外乱」の影響を受ける状況のもとで,「設定値」の保持や,ときには変更を要求されるからである。このような問題を解決するためには,これまでの制御方式に加えて,まず「積分動作」の導入による構造的な改良が必要となる。本章では,定値外乱の影響の除去や,設定値の変更に対応できるように,積分動作を入れて安定化を行う方策について学ぶ。その安定化に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.5w=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)=1w=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)=0w=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) }

を達成したい。ここで,定数ベクトルrm個の設定値からなる。もし制御目的(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] }

の関係を満足しているはずである。したがって,どのようなwrに対しても,x_\inftyu_\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設計

<

これまで,閉ループ系(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}})は可観測対とする。

さらに,FF_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>0V>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_SFF_IHから,つぎを構成する。

\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)の妥当性について説明しなさい。

【本章のねらい】
・1次系,2次系の状態フィードバックの最適設計を,リッカチ方程式を解いて行う。
・ リッカチ方程式を解く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. }

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

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

を設定し,これを最小化する問題を解くことである。ただし,Q>0R>0QRは正定行列。テキスト「線形システム制御入門」の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 }

まず,1次系の場合を考える。

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

\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=1K=1とする。このとき,つぎの重み係数をもつ評価関数を最小にするfを決定せよ。
(1) q=1, r=1
(2) q=1, r=0.5
(3) q=1, r=0.1
また,x(0)=1のとき,閉ループ系の時間応答をMATLABでシミュレーションせよ。

つぎに,2次系の場合を考える。

例題5.2 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が2つ,つぎに第3式より\pi_2が2つ,さらに第2式より\pi_1が1つ定まり,つぎのように4組の解をもつ。すなわち

\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 つぎの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 2次系

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

を得る。まず,第1式より\pi_3

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

と求まる。つぎに,第3式より\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}) }

でなければならない。さらに,第2式より\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 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) }

について,例題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_1q_2の与え方の指針を得ることができる(テキスト「線形システム制御入門」の第5章演習問題【1】参照。また(??),(??)の導出についてはそれぞれ第5章例題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;

ここで,4行目でハミルトン行列の固有値問題を解き,安定固有値の添字を5行目の 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 3次系

\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 4次系

\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 オブザーバベースト・コントローラの最適設計

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


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

ここで,新しい入力wvがそれぞれW>0V>0の平方根行列(X>0(\ge0)に対しYY=Xを満足する行列Y>0(\ge0)X^{\frac{1}{2}}で表す。)
により重み付けられて,n次系の入力側(B'を介して)と出力側に設置されている。また新しい出力z=C'xと入力uが取り出されており,
それぞれQ>0R>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がつぎの手順で決定できる
(この手順はLQG(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 1次系\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=1r=1\sigma=1\rho=1と選んで構成せよ。
解答 CAREに,a=0b=1q=1r=1を代入して,\Pi^2=1\Pi>0より,\Pi=1
したがって,状態フィードバックゲインは,f=\frac{1}{r^2}b\Pi=1。FAREに,a=0c=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 1次系\dot{x}(t)=ax(t)+bu(t)\,(a=-1,b=1)に対するオブザーバベースト・コントローラを,
例題5.5のCAREにおいてq=1r=1,FAREにおいて\sigma=1\rho=0.2と選んで構成せよ。

例題5.6 LQG制御則を設計するための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_1f_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=-1b=1q=1r=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=-1c=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 状態オブザーバと可観測性

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

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)の第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 2次系

\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 2次系

\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 1次系

\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-bfa-hcとなる。

演習4.2 例題4.3において,a=-1,b=1,c=,f=4,x(0)=1のとき,状態フィードバックによる閉ループ系とオブザーバベースト・コントローラによる閉ループ系の応答を比較したい。h=2h=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 (\lambdaAのすべての不安定固有値)
条件D2: Cv=0,\ Av=\lambda v \Rightarrow v=0 (\lambdaAのすべての不安定固有値)

これらの条件の一つが成り立つとき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 (\lambdaAのすべての固有値)
条件O5: Cv=0,\ Av=\lambda v \Rightarrow v=0 (\lambdaAのすべての固有値)
条件O6: (A^T,C^T)は可制御対

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

例題4.4 つぎのA行列とC行列をもつ2次系\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で,システムの次数2と等しい。したがって,この2次系は可観測である。

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

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

演習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,出力数は2なので,状態オブザーバを用いた場合は4次元であるが,恒等関数オブザーバを用いれば2次元に,線形関数オブザーバを用いれば1次元に減らすことができる。)。

例題4.5 1入力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}であり,条件(??)の第2式と第4式が満足されることを示せばよい。実際,

\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の2次系に対する恒等関数オブザーバを,その固有値が\lambda=-2となるように求めよ。この恒等オブザーバに対して例題4.2と同様のシミュレーションを行え。

例題4.6 例題4.5の1入力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)の第2式と第4式が満足されることを示せばよい。実際,

\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の2次系に対して,状態フィードバックu(t)=-y(t)-2\dot{y}(t)を考える。このとき,この状態フィードバックを実施する1次元線形関数オブザーバを,その固有値が\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次系は可観測である。

(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)

したがって,この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

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

【演習4.5】 A_{21}=A_{22}=0B_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}=0B_2=1K_1=-1K_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=-2L=-\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 状態フィードバックと可制御性

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

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に戻ることを意味するからである。

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

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

\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'をもつ1次系を表している。

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

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

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

例題3.2 減衰係数\zetaと固有角周波数\omega_n^2をもつ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をもつ2次系を表している。

演習3.3 例題2.6で得た図上(ty平面)に,望ましい第1番目のオーバーシュートの頂点の座標(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次系に対する状態フィードバックの設計法を考える。まず,1入力系の場合を考え,つぎの条件を仮定する(1入力系の可制御性行列は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 2次系

\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 つぎの2次系に対する状態フィードバック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'_nAの固有値に等しくないとし,m次元ベクトルg_1,\cdots,g_nは,上の逆行列が存在する範囲で適切に指定するものとする。これから,多入力系の場合,状態フィードバックは,固有値を指定しただけでは,一意に定まらないことがわかる。

例題3.4 2入力2次系

\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の2つの状態フィードバックは,公式(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 (\lambdaAのすべての不安定固有値)
条件S2: B^Tw=0,\ A^Tw=\lambda w \Rightarrow w=0 (\lambdaAのすべての不安定固有値)

これらの条件の一つが成り立つとき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 (\lambdaAのすべての固有値)
条件C5: B^Tw=0,\ A^Tw=\lambda w \Rightarrow w=0 (\lambdaAのすべての固有値)

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

例題3.5 つぎのA行列とB行列をもつ2次系\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で,システムの次数2と等しい。したがって,この2次系は可制御である。

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

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

演習3.6 つぎのA行列とB行列をもつ3次系\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

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

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

例題3.6 例題3.5のA行列とB行列をもつ2次系\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次系は可安定である。

(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)

したがって,この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

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

演習3.8 演習3.6のA行列とB行列をもつ3次系\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行列をもつ3次系\dot{x}(t)=Ax(t)+Bu(t)の可安定性を判定せよ。

演習問題の解答

演習3.1 T=1,K=1T'=1/10,K'=1より,f=\frac{1}{1}(\frac{1}{1/0.1}-\frac{1}{1})=9g=\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で,システムの次数3より小さい。したがって,この3次系は不可制御である。

(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で,システムの次数3と等しい。
したがって,この3次系は可制御である。

演習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

したがって,この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次系は可安定である。

演習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 漸近安定性

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

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 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)
明らかに,1次系が漸近安定である条件は,a<0である。

演習2.1 つぎの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 つぎの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]}

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

\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.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のように,1次系が零固有値をもつ場合(a=0,積分器),応答は発散することはない。一方,演習2.2(3)が示すように,2次系で2つの零固有値をもつ場合(2つの積分器が直列結合),入力側の積分器の初期値が零でないならばこれを積分する出力側の積分器の応答は発散する。このため一般には零固有値は不安定固有値とみなす)。

例題2.3 つぎの行列Aをもつ2次系\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つの固有値が実数で負であるので漸近安定である。

(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と求められる。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つの固有値の実部が負(正)であるので漸近安定である(ない)

演習2.3 例題2.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参照)。

以下では,簡単のために,1入力1出力の場合を考える。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' }

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

まず,次式で表される1次系を考える。

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

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

解答 インパルス応答は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次系のステップ応答を求めよ。

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

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

\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 つぎの1次系のステップ応答を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)を使って読み取れ。

つぎに,次式で表される2次系を考える。

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

したがって,図上で第1番目のオーバーシュートの頂点の座標(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 つぎの2次系のステップ応答を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の図上で,第1番目のオーバーシュートの頂点の座標をginput(1)を使って読み込み,固有角周波数と減衰係数を同定せよ。

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

例題2.7 つぎの1次系の時間応答を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 つぎの1次系のボード線図を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の2次系のボード線図をMATLABで図示せよ。

最後に,自由系の時間応答(初期値応答)の計算法を示す。

例題2.9 つぎの1次自由系の時間応答を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)の2次系のインパルス応答を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つとも実数で負だから漸近安定。
(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。2つとも実数で負だから漸近安定。

演習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の固有値は2つとも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.したがって,横線1-\frac{1}{e}を図示すればよい。

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

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

演習2.6

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

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

演習2.7

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

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

演習2.8

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

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

演習2.9

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

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

1 状態空間表現

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

ここで,A行列,B行列,C行列,D行列を参照するには,それぞれ

sys.a,sys.b,sys.c,sys.d

を用いる。たとえばA行列の(2,2)要素を-0.1に変更するには,つぎのコマンドを与えればよい。

sys.a(2,2)=-0.1

また,B行列,D行列が存在しない自由系 sys0 を定義するには,つぎのコマンドを与えればよい。

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

もし状態方程式だけを指定する場合は,C=[]; D=[]; または C=eye(size(A)); D=0; とする。

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

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

ここで,A行列,B行列,C行列,D行列を参照するには,それぞれ

sys.a,sys.b,sys.c,sys.d

を用いる。たとえばA行列の(2,2)要素を-0.1に変更するには,つぎのコマンドを与えればよい。

A(2,2)=-0.1; sys.a=A

また,B行列,D行列が存在しない自由系 sys0 を定義するには,つぎのコマンドを与えればよい。

//state_space_cont.m
A=[0 1;0 -0.01]; B=[]; C=[1 0]; D=[];
sys0=syslin(A,B,C,D)

もし状態方程式だけを指定する場合は,C=[]; D=[]; または C=eye(size(A)); D=zeros(size(A,1),size(B,2)); とする。

さて,状態空間表現のブロック線図を作成することにより,どのように入力変数から状態変数を経て出力変数につながるかを表すことができる。これは,Simulinkなどを用いた時間応答シミュレーションに役に立つ。

例題1.3 つぎの1次系の状態空間表現をブロック線図で表せ。

\displaystyle{ \left\{\begin{array}{l} \dot{x}(t)=-\frac{1}{T}x(t)+\frac{K}{T}u(t) \\ y(t)=cx(t) \end{array}\right. }

解答 積分器の出力をxとすると,その入力は\dot{x}であることに注意し,状態方程式の左辺\dot{x}と右辺-\frac{1}{T}x+\frac{K}{T}uを等しくなるよう描いて,つぎのブロック線図を得る(加え合せ点○で符号を表示しない場合は+とみなす。また,●は引き出し点を表す。)。

演習1.5 つぎの1次系の状態空間表現をブロック線図で表せ。

\displaystyle{ \left\{\begin{array}{l} \dot{x}(t)=u(t) \\ y(t)=x(t)+u(t) \end{array}\right. }

例題1.4 つぎの2次系の状態空間表現をブロック線図で表せ。

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

解答 \dot{x}_1(t)=x_2(t),\ \dot{x}_2(t)=-\omega_n^2x_1(t)-2\zeta\omega_nx_2(t)+\omega_n^2u(t),\ y(t)=x_1(t)より,つぎのブロック線図を得る。

演習1.6 つぎの2次系の状態空間表現をブロック線図で表せ。

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

演習1.7 演習1.2 で得られた状態空間表現について,ブロック線図を描け。

1.2 状態空間表現の座標変換と直列結合

n次系の状態空間表現

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

に対して,座標変換

\displaystyle{(1.3)\quad x'(t)=Tx(t)\quad({\rm det}T\ne0) }

を行うと,つぎの状態空間表現を得る。

\displaystyle{(1.4)\quad \left\{\begin{array}{l} \dot{x}'(t)=A'x'(t)+B'u(t) \\ y(t)=C'x'(t) \end{array}\right. }

ただし,変換後の係数行列は次式で与えられる。

\displaystyle{(1.5)\quad \left\{\begin{array}{l} A'=TAT^{-1} \\ B'=TB \\ C'=CT^{-1} \\ \end{array}\right. }

例題1.5 2次系の状態空間表現

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

に対して,つぎの座標変換を行え。

\displaystyle{ \left[\begin{array}{c} x_1'(t) \\ x_2'(t) \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ 1 & -a_1 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] }

解答 状態空間表現に

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

を代入して

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

この状態方程式の左から

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

をかけて,つぎの状態空間表現を得る。

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

演習1.8 2次系の状態空間表現

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

に対して,つぎの座標変換を行え。

\displaystyle{ \left[\begin{array}{c} x_1'(t) \\ x_2'(t) \end{array}\right] = \left[\begin{array}{cc} 1 & 0 \\ -\zeta\omega_n & \omega_n\sqrt{1-\zeta^2} \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] }

演習1.9 2次系の状態空間表現

\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right] = \left[\begin{array}{cc} a_{11} & a_{12} \\ a_{21} & a_{22} \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] + \left[\begin{array}{c} b_{1} \\ b_{2} \end{array}\right] u(t) \nonumber \\ y(t) = \underbrace{ \left[\begin{array}{cc} c_1 & c_2 \end{array}\right] }_C \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] \quad(c_1\ne0) \end{array}\right. }

に対して,どのような座標変換

\displaystyle{ \left[\begin{array}{c} x_1'(t) \\ x_2'(t) \end{array}\right] = \underbrace{ \left[\begin{array}{cc} t_{11} & t_{12} \\ t_{21} & t_{22} \end{array}\right] }_T \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] }

を行なえば,出力方程式をつぎの形にできるか。

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

さて,つぎの2つの状態空間表現を考える。

\displaystyle{(1.6)\quad S_1: \left\{\begin{array}{l} \dot{x}_1(t)=A_1x_1(t)+B_1u_1(t) \\ y_1(t)=C_1x_1(t)+D_1u_1(t) \end{array}\right. }

\displaystyle{(1.7)\quad S_2: \left\{\begin{array}{l} \dot{x}_2(t)=A_2x_2(t)+B_2u_2(t) \\ y_2(t)=C_2x_2(t)+D_2u_2(t) \end{array}\right. }

S_1の入力変数の数とS_2の出力変数の数が等しいとき,S_2の出力をS_1の入力に直接結合することを考える。その状態空間表現は次式で与えられる。

\displaystyle{(1.8)\quad \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \end{array}\right] = \left[\begin{array}{cc} A_1 & B_1C_2 \\ 0 & A_2 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] + \left[\begin{array}{c} B_1D_2 \\ B_2 \end{array}\right] u_2(t) \\ y_1(t)= \left[\begin{array}{cc} C_1 & D_1C_2 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] +D_1D_2u_2(t) \end{array}\right. }

例題1.6 1次系の状態空間表現

\displaystyle{ \left\{\begin{array}{l} \dot{x}(t)=-\frac{1}{T}x(t)+\frac{K}{T}u(t) \\ y(t)=x(t) \end{array}\right. }

の入力に,むだ時間を1次系として近似したときの状態空間表現

\displaystyle{ \left\{\begin{array}{l} \dot{x}_L(t)=-\frac{2}{L}x_L(t)+\frac{2}{L}u_L(t) \\ y_L(t)=2x_L(t)-u_L(t) \end{array}\right. }

の出力を結合して得られるシステムの状態空間表現を求めよ。
解答 まず,状態方程式は

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

u(t)=y_L(t)=2x_L(t)-u_L(t)を代入し,

\displaystyle{ \dot{x}_L(t)=-\frac{2}{L}x_L(t)+\frac{2}{L}u_L(t) }

と合わせて

\displaystyle{ \left[\begin{array}{c} \dot{x}(t) \\ \dot{x}_L(t) \end{array}\right] = \left[\begin{array}{ccc} -\frac{1}{T} & 2\frac{K}{T} \\ 0 & -\frac{2}{L} \end{array}\right] \left[\begin{array}{c} x(t) \\ x_L(t) \end{array}\right] + \left[\begin{array}{c} -\frac{K}{T} \\ \frac{2}{L} \end{array}\right] u_L(t) }

のように得られる。また,出力方程式は

\displaystyle{ y(t)= \left[\begin{array}{ccc} 1 & 0 \end{array}\right] \left[\begin{array}{c} x(t) \\ x_L(t) \end{array}\right] }

のように得られる。

演習1.10 2次系

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

の入力に,むだ時間を2次系として近似したときの状態空間表現

\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_{L1}(t) \\ \dot{x}_{L2}(t) \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ -\frac{12}{L^2} & -\frac{6}{L} \end{array}\right] \left[\begin{array}{c} x_{L1}(t) \\ x_{L2}(t) \end{array}\right] + \left[\begin{array}{c} 0 \\ 1 \end{array}\right] u_L(t) \nonumber \\ y_L(t) = \left[\begin{array}{cc} 0 & -\frac{12}{L} \end{array}\right] \left[\begin{array}{c} x_{L1}(t) \\ x_{L2}(t) \end{array}\right] +u_L(t) \end{array}\right. }

の出力を結合して得られるシステムの状態空間表現を求めよ。

MATLABを用いて座標変換を行うために,たとえば演習1.8で,\omega_n=1,\zeta=-0.01の場合,つぎのコマンドを与えればよい。

%coordinate_transformation.m
om=1; zeta=0.01;
A=[0 1;-om^2 -2*zeta*om]; B=[0;om^2]; C=[1 0]; D=0;
sys1=ss(A,B,C,D);
T=[1 0;-zeta*om om*sqrt(1-zeta^2)];
sys2=ss2ss(sys1,inv(T))

SCILABを用いて座標変換を行うために,たとえば演習1.8で,\omega_n=1,\zeta=-0.01の場合,つぎのコマンドを与えればよい。

//coordinate_transformation.m
om=1; zeta=0.01;
A=[0 1;-om^2 -2*zeta*om]; B=[0;om^2]; C=[1 0]; D=0;
sys1=syslin(‘c’,A,B,C,D);
T=[1 0;-zeta*om om*sqrt(1-zeta^2)];
sys2=ss2ss(sys1,T)

つぎに,MATLABを用いて直列結合を行うためには,たとえば演習1.10でL=1の場合,つぎのコマンドを与えればよい。

%serial_connection.m
A1=[0 1;0 0]; B1=[0;1]; C1=[1 0]; D1=0;
sys1=ss(A1,B1,C1,D1);
L=1;
A2=[0 1;-12/L^2 -6/L]; B2=[0;1]; C2=[0 -12/L]; D2=1;
sys2=ss(A2,B2,C2,D2);
sys=sys1*sys2

ここで, sys1 と sys2 の順番に注意する。

最後に,状態空間表現から伝達関数行列表現(8章の(8.2)式を参照。)を求めるには,コマンド tf を用いる。たとえば,上のむだ時間を2次系として近似したときの状態空間表現 sys2 の伝達関数を求めるには,コマンド tf(sys2) を与えればよい。

つぎに,SCILABを用いて直列結合を行うためには,たとえば演習1.10でL=1の場合,つぎのコマンドを与えればよい。

//serial_connection.m
A1=[0 1;0 0]; B1=[0;1]; C1=[1 0]; D1=0;
sys1=syslin(‘c’,A1,B1,C1,D1);
L=1;
A2=[0 1;-12/L^2 -6/L]; B2=[0;1]; C2=[0 -12/L]; D2=1;
sys2=syslin(‘c’,A2,B2,C2,D2);
sys=sys1*sys2

ここで, sys1 と sys2 の順番に注意する。

最後に,状態空間表現から伝達関数行列表現(8章の(8.2)式を参照。)を求めるには,コマンド ss2tf を用いる。たとえば,上のむだ時間を2次系として近似したときの状態空間表現 sys2 の伝達関数を求めるには,コマンド ss2tf(sys2) を与えればよい。

演習問題の解答

演習1.1 \omega=\dot{\theta}とおくと,\dot{\theta}=\omega, \dot{\omega}=-\frac{3g}{4\ell}\theta。したがって,つぎの状態方程式を得る。

\displaystyle{ \left[\begin{array}{c} \dot{\theta} \\ \dot{\omega} \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ -\frac{3g}{4\ell} & 0 \end{array}\right] \left[\begin{array}{c} \theta \\ \omega \end{array}\right] }

演習1.2 v_1=\dot{x}_1, v_2=\dot{x}_2とおくと,\dot{x}_1=v_1, \dot{v}_1=-(x_1-x_2)+u\dot{x}_2=v_2, \dot{v}_2=-(x_2-x_1)。したがって,つぎの状態方程式を得る。

\displaystyle{ \left[\begin{array}{c} \dot{x}_1 \\ \dot{v}_1 \\ \dot{x}_2 \\ \dot{v}_2 \end{array}\right] = \left[\begin{array}{cccc} 0 & 1 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & -1 & 0 \\ \end{array}\right] \left[\begin{array}{c} x_1 \\ v_1 \\ x_2 \\ v_2 \end{array}\right] + \left[\begin{array}{c} 0 \\ 1 \\ 0 \\ 0 \end{array}\right] u }

または

\displaystyle{ \left[\begin{array}{c} \dot{x}_1 \\ \dot{x}_2 \\ \dot{v}_1 \\ \dot{v}_2 \end{array}\right] = \left[\begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \\ \end{array}\right] \left[\begin{array}{c} x_1 \\ x_2 \\ v_1 \\ v_2 \end{array}\right] + \left[\begin{array}{c} 0 \\ 0 \\ 1 \\ 0 \end{array}\right] u }

演習1.3 一定のトルク\tau^*のもとで一定の角速度\omega^*を得ている。これらは運動方程式\dot{\omega}=-{\omega}+\tauを満足するので,0=-\omega^*+\tau^*。これを運動方程式から辺々引き算して,つぎの1次系としての状態空間表現を得る。

\displaystyle{ \left\{\begin{array}{l} \underbrace{\frac{d}{dt}(\omega-\omega^*)}_{\dot x} =\underbrace{-1}_{a}\underbrace{(\omega-\omega^*)}_{x} +\underbrace{1}_{b} \underbrace{(\tau-\tau^*)}_{u} \\ \underbrace{\omega-\omega^*}_{y} =\underbrace{1}_{c}\underbrace{\omega-\omega^*}_{x} \end{array}\right. }

演習1.4 x_1={y}, x_2=\dot{y}とおくと,\dot{x}_1=x_2, \dot{x}_2=-a_1x_2-a_2x_1+u。したがって,つぎの状態空間表現を得る。

\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1 \\ \dot{x}_2 \\ \end{array}\right] = \left[\begin{array}{cccc} 0 & 1 \\ -a_2 & -a_1 \\ \end{array}\right] \left[\begin{array}{c} x_1 \\ x_2 \\ \end{array}\right] + \left[\begin{array}{c} 0 \\ 1 \\ \end{array}\right] u \\ y = \left[\begin{array}{cccc} 1 & 0 \end{array}\right] \left[\begin{array}{c} x_1 \\ x_2 \\ \end{array}\right] \end{array}\right. }

または

\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_2 \\ \dot{x}_1 \\ \end{array}\right] = \left[\begin{array}{cccc} -a_1 & -a_2 \\ 1 & 0 \\ \end{array}\right] \left[\begin{array}{c} x_2 \\ x_1 \\ \end{array}\right] + \left[\begin{array}{c} 1 \\ 0 \\ \end{array}\right] u \\ y = \left[\begin{array}{cccc} 0 & 1 \end{array}\right] \left[\begin{array}{c} x_2 \\ x_1 \\ \end{array}\right] \end{array}\right. }

演習1.5 \dot{x}=u,\ y=x+uより,つぎのブロック線図を得る。

演習1.6 \dot{x}_1=x_2,\ \dot{x}_2=-2\zeta\omega_nx_2+\omega_n^2u,\ y=c_2x_2より,つぎのブロック線図を得る。

演習1.7 \dot{v}_1=-(x_1-x_2)+u,\ \dot{v}_2=-(x_2-x_1)より,つぎのブロック線図を得る。

演習1.8 状態空間表現に

\displaystyle{ \left[\begin{array}{c} x_1 \\ x_2 \end{array}\right] = \left[\begin{array}{cc} 1 & 0 \\ -\zeta\omega_n & \omega_n\sqrt{1-\zeta^2} \end{array}\right] \left[\begin{array}{c} x_1' \\ x_2' \end{array}\right] }

を代入して

\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{cc} 1 & 0 \\ -\zeta\omega_n & \omega_n\sqrt{1-\zeta^2} \end{array}\right] \left[\begin{array}{c} \dot{x}_1' \\ \dot{x}_2' \end{array}\right] \nonumber \\ = \left[\begin{array}{cc} 0 & 1 \\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] \left[\begin{array}{cc} 1 & 0 \\ -\zeta\omega_n & \omega_n\sqrt{1-\zeta^2} \end{array}\right] \left[\begin{array}{c} x_1' \\ x_2' \end{array}\right] + \left[\begin{array}{c} 0 \\ 1 \end{array}\right] u \nonumber \\ y= \left[\begin{array}{cc} 1 & 0 \end{array}\right] \left[\begin{array}{cc} 1 & 0 \\ -\zeta\omega_n & \omega_n\sqrt{1-\zeta^2} \end{array}\right] \left[\begin{array}{c} x_1' \\ x_2' \end{array}\right] \end{array}\right. }

この状態方程式の左から

\displaystyle{ \left[\begin{array}{cc} 1 & 0 \\ -\zeta\omega_n & \omega_n\sqrt{1-\zeta^2} \end{array}\right]^{-1} = \frac{1}{\omega_n\sqrt{1-\zeta^2}} \left[\begin{array}{cc} \omega_n\sqrt{1-\zeta^2} & 0 \\ \zeta\omega_n & 1 \end{array}\right] }

をかけて,つぎの状態空間表現を得る。

\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1' \\ \dot{x}_2' \end{array}\right] = \left[\begin{array}{cc} -\zeta\omega_n & \omega_n\sqrt{1-\zeta^2} \\ -\omega_n\sqrt{1-\zeta^2} & -\zeta\omega_n \end{array}\right] \left[\begin{array}{c} x_1' \\ x_2' \end{array}\right] + \left[\begin{array}{c} 0 \\ \frac{1}{\omega_n\sqrt{1-\zeta^2}} \end{array}\right] u \nonumber \\ y= \left[\begin{array}{cc} 1 & 0 \end{array}\right] \left[\begin{array}{c} x_1' \\ x_2' \end{array}\right] \end{array}\right. }

演習1.9

\displaystyle{ \underbrace{ \left[\begin{array}{cc} 1 & 0 \end{array}\right] }_{C'} \underbrace{ \left[\begin{array}{cc} t_{11} & t_{12} \\ t_{21} & t_{22} \end{array}\right] }_T = \underbrace{ \left[\begin{array}{cc} c_1 & c_2 \end{array}\right] }_C \quad(Tは正則) }

を満足させればよいので,たとえばt_{11}=c_1(\ne0), t_{12}=c_2, t_{21}=0, t_{22}=1と選べばよい。

演習1.10

\displaystyle{ \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \\\hdashline \dot{x}_{L1}(t) \\ \dot{x}_{L2}(t) \end{array}\right] = \left[\begin{array}{cc:cc} 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & -\frac{12}{L} \\\hdashline 0 & 0 & 0 & 1 \\ 0 & 0 & -\frac{12}{L^2} & -\frac{6}{L} \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \\ \hdashline x_{L1}(t) \\ x_{L2}(t) \end{array}\right] + \left[\begin{array}{c} 0 \\ 1 \\ \hdashline 0 \\ 1 \end{array}\right] u_L(t) \\ y(t)= \left[\begin{array}{cc:cc} 1 & 0 & 0 & 0 \end{array}\right] \left[\begin{array}{c} x_1(t) \\ x_2(t) \\ \hdashline x_{L1}(t) \\ x_{L2}(t) \end{array}\right] \end{array}\right. }

UV Control

●Drone becomes to be very popular tool to take a video from the sky. It is a surprise that we can control its motion with 6 degrees of freedom so smartly. What makes it possible?

In general, a drone has 4 propellers as shown. This means that 4 degrees of freedom can be controlled. They should be the heave z and the roll \phi, the pitch \theta, the yaw \psi. Assume that rotational speeds of 4 propellers can be set to the specified values, \omega_{1c}, \omega_{2c}, \omega_{3c}, \omega_{4c}.

Then we can produce 1 force F_z and 3 torques \tau_\phi, \tau_\theta, \tau_\psi to control the heave z and the roll \phi, the pitch \theta, the yaw \psi respectively by combining the 4 propellers’ thrusts as follows:

\left[\begin{array}{c} F_z \\ \tau_\phi \\ \tau_\theta \\ \tau_\psi \end{array}\right]= \left[\begin{array}{cccc} k_1 & k_1 & k_1 & k_1 \\ 0 & -\ell k_1 & 0 & \ell k_1 \\ \ell k_1 & 0 & -\ell k_1 & 0 \\ -k_2 & k_2 & -k_2 & k_2 \end{array}\right] \left[\begin{array}{c} \omega_{1c} \\ \omega_{2c} \\ \omega_{3c} \\ \omega_{4c}  \end{array}\right].

This means that the MIMO dynamics is decoupled into 4 SISO dynamics. We will employ a PID controller corresponding to each SISO dynamics. Therefore it is the static decoupling to make the control system design easier.

●Many underwater vehicles have been developed for various purposes. For example, we can find them as follows.

Note that 4 to 8 thrustors are equipped to realize the static decoupling. The number depends on motion necessary to achieve the mission of UVs. Most control systems will be designed in the same way as drone.

●Prof. Koterayama and Prof. Nakamura developed the following UV called as DELTA.

The mission of DELTA is for wide-area survey. It works in a towing mode for wade scanning, and in a self-propulsive mode for local investigation. The detail is discussed in the paper:

Masahiko NAKAMURA, Hiroyuki KAJIWARA, Wataru KOTERAYAMA:
Development of an ROV operated both as towed and self-propulsive vehicle
Ocean Engineering, vol.28, no.1, pp.1-43, 2000.

DELTA has only two thrustors. So only the heave z and the yaw \psi are controlled. However it is necessary to control the roll \phi and the pitch \theta to keep the equilibrium state. Such a control system can’t be realized by PID control. The following control system is designed by LQG approach.

For any UV, the collision avoidance is an important function. For example, a deep diving will be necessary to avoid an obstacle as follows.

In order to change the depth in DELTA, there is no way except decreasing the thrust. This means the large speed variation, by which hydrodynamic force will be varied largely. So we need a robust controller or a scheduling controller. Such a control scheme is depicted in the following.

Ship Control

●For the directional control of a ship, consider the following NOMOTO model.

\displaystyle{ \dot{r}(t)=-\frac{1}{T(t)}r(t)+\frac{K(t)}{T(t)}\delta(t-L)+w(t)} }

where r(t)=\dot{\psi}(t) and

\displaystyle{ T(t)=\frac{L}{U(t)}T',\ K(t)=\frac{U(t)}{L}K' \quad (U_1\le U(t)\le U_2) }

We assume the time lag L to take account of the delay to get the rudder effect. We use the same notation for the ship length L assuming no confusion. The time constant T and gain constant K are varied according to the variation of forward speed U. From the following equation, it is clear that the rudder effect is proportional to the square of the forward speed.

\displaystyle{ \dot{r}(t)=-\underbrace{\left(\frac{U(t)}{L}\right)\frac{1}{T'}}_{a(t)=\frac{1}{T(t)}}r(t) +\underbrace{\left(\frac{U(t)}{L}\right)^2\frac{K'}{T'}}_{b(t)=\frac{K(t)}{T(t)}}\delta(t-L)+w(t) }

Given a nominal forward speed U^*, the following equation holds.

\displaystyle{ \dot{r}(t)=-\underbrace{\left(\frac{U(t)}{U^*}\right)\frac{1}{T^*}}_{a(t)=\frac{1}{T(t)}}r(t) +\underbrace{\left(\frac{U(t)}{U^*}\right)^2\frac{K^*}{T^*}}_{b(t)=\frac{K(t)}{T(t)}}\delta(t-L)+w(t) }

where

\displaystyle{ T^*=\frac{L}{U^*}T',\ K^*=\frac{U^*}{L}K' \quad (U_1\le U^*\le U_2) }

Zig-Zag Test: Neglecting the time lag L and the forward speed variation, a simulation example of the Zig-Zag test is given as follows.

From this response, we used to obtain both time constant T and gain constant K. However, the response is not always the step response by which a linear system is completely characterized. As a ship has dynamics of an astatic (no-fixed-position) system, the Zig-Zag test is a compromise way for its identification.

We should notice that if we apply the unity feedback as shown in the above we can obtain the step response and identify (\zeta,\omega_n) from (T_p, 1+p_0) by the following formula.

\displaystyle{ (T_p, 1+p_0)\ \Rightarrow\  (\zeta,\omega_n)= \left( \sqrt{\frac{(\log p_0)^2}{(\log p_0)^2+\pi^2}}, \frac{\sqrt{(\log p_0)^2+\pi^2}}{T_p} \right) }

PID control: Under the nominal forward speed U^*, neglecting the time lag L, consider the following NOMOTO model.

\displaystyle{ \dot{r}(t)=-\frac{1}{T^*}r(t)+\frac{K^*}{T^*}\delta(t)+w(t) }

Then PID control is presented as follows.

\displaystyle{ \delta(t)=K_P(\psi_c-\psi(t))+K_D\frac{d}{dt}(\psi_c-\psi(t))+K_I\int_0^t(\psi_c-\psi(\tau))\,d\tau }

Although we used to tune the P gain K_P, the D gain K_D, the I gain K_I on the site, here we want try to determine them based on the NOMOTO model. The closed-loop system by PID control is given by

Here assuming \dot{\psi}_c=0, we have

\displaystyle{ \delta(t)=-K_P\psi(t)-K_Dr(t)+K_P\psi_c+K_I\int_0^t(\psi_c-\psi(\tau))\,d\tau }

Therefore, we have another configuration of PID control system.

Letting \zeta=\frac{1}{2\sqrt{T^*K^*}},\omega_n=\sqrt{\frac{K^*}{T^*}}, the following state equation is obtained.

\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ 0 & -2\zeta\omega_n \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ \omega_n^2 \end{array}\right] }_{B} \underbrace{\delta(t)}_{u(t)} + \left[\begin{array}{c} 0 \\ w \end{array}\right] }

Then PID control is a state feed back with an integral action as follows.

\displaystyle{ \underbrace{\delta(t)}_{u(t)} =- \underbrace{ \left[\begin{array}{cc} K_P & K_D \end{array}\right] }_{F} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \end{array}\right] }_{x(t)} +\underbrace{K_P}_{G} \underbrace{\psi_c}_{v} +K_I \underbrace{\int_0^t(\psi_c-\psi(\tau))\,d\tau}_{x_I(t)} }

Therefore, after specifying the desired first overshoot (T_p',1+p_0'), we calculate the corresponding to (\zeta',\omega_n') and then determine PID gains by

\displaystyle{ K_P=\frac{\omega_n'^2}{\omega_n^2},\  K_D=\frac{2}{\omega_n^2}(\zeta'\omega_n'-\zeta\omega_n),\  K_I=\omega_I\frac{\omega_n'^2}{\omega_n^2} }

Here we try \omega_I=\frac{w_n'}{10} at first, and then tune it by observing the disturbance rejection.

The CLPS by PID control is presented by

\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \\ \dot{x}_I(t) \end{array}\right] }_{\dot{x}_{CL}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ -\omega_n'^2 & -2\zeta'\omega_n' & \omega_I\omega_n'^2 \\ -1 & 0 & 0 \end{array}\right] }_{A_{CL}} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ x_I(t) \end{array}\right] }_{x_{CL}(t)} + \underbrace{ \left[\begin{array}{cc} 0 & 0\\ \omega_n'^2 & 1\\ 1 & 0 \end{array}\right] }_{B_{CL}} \left[\begin{array}{c} \psi_c \\ w \end{array}\right] }

Here, assuming \dot{w}=0, we differentiate it and obtain

\displaystyle{ \underbrace{ \left[\begin{array}{c} \ddot{\psi}(t) \\ \ddot{r}(t) \\ \ddot{x}_I(t) \end{array}\right] }_{\ddot{x}_{CL}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ -\omega_n'^2 & -2\zeta'\omega_n' & \omega_I\omega_n'^2 \\ -1 & 0 & 0 \end{array}\right] }_{A_{CL}} \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \\ \dot{x}_I(t) \end{array}\right] }_{\dot{x}_{CL}(t)} }

From the CLPS stability, we have

\displaystyle{ \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ x_I(t) \end{array}\right] }_{x_{CL}} \rightarrow \underbrace{ \left[\begin{array}{ccc} 0 & -1 & 0\\ \omega_n'^2 & 2\zeta'\omega_n' & -\omega_I\omega_n'^2 \\ 1 & 0 & 0 \end{array}\right]^{-1} }_{-A_{CL}^{-1}} \underbrace{ \left[\begin{array}{cc} 0 \\ \omega_n'^2\psi_c +w\\ \psi_c \end{array}\right] }_{B_{CL}} = \left[\begin{array}{cc} \psi_c\\ 0\\ \frac{w}{\omega_I\omega_n'^2} \end{array}\right] }

This means that the integral action is estimating the unknown disturbance.

The following simulation shows the CLPS responses by PD control. It is observed that the disturbance rejection is not sufficient because of the lack of integral action.

The following simulation shows the CLPS responses by PID control. Here \omega_I=2\frac{w_n'}{10}. It is observed that the disturbance rejection is satisfied but the transient response is fractured a little bit being off expectation.

LQI control: We will try to determine the PID gains in the framework of LQI control. Consider the following error system:

\displaystyle{ \left[\begin{array}{c} \dot\psi(t)-\dot{\psi}_c \\ \dot r(t) \\ \dot{\delta}(t) \end{array}\right] = \left[\begin{array}{ccc} 0 & 1 & 0 \\ 0 & -2\zeta\omega_n &\omega_n^2 \\ 0 & 0 & 0 \end{array}\right] \left[\begin{array}{c} \psi(t)-\psi_c \\ r(t) \\ \delta(t) \end{array}\right] + \left[\begin{array}{c} 0 \\ 0 \\ 1 \end{array}\right] \dot{\delta} (t) }

For the error system, we will minimize the following criterion function:

\displaystyle{ J=\int_0^\infty(\frac{1}{M_\psi^2}(\psi(t)-\psi_c)^2+\frac{T_\psi^2}{M_\psi^2}r^2(t) +\frac{1}{M_\delta^2}\delta^2(t)+\frac{T_\delta^2}{M_\delta^2}\dot{\delta}^2(t))\,dt }

Solving this problem, the PID gains are obtained as follows:

\displaystyle{ \dot{\delta}(t) =- \left[\begin{array}{ccc} f_\psi & f_r & f_\delta \end{array}\right] \left[\begin{array}{c} \psi(t)-\psi_c \\ r(t) \\ \delta(t) \end{array}\right] }
\displaystyle{ =- \underbrace{ \left[\begin{array}{ccc} f_\psi & f_r & f_\delta \end{array}\right] \left[\begin{array}{ccc} 0 & 1 & 0 \\ 0 & -2\zeta\omega_n &\omega_n^2 \\ 1 & 0 & 0 \end{array}\right]^{-1} }_{ \left[\begin{array}{ccc} K_P & K_D & K_I \end{array}\right] } \left[\begin{array}{c} \dot\psi(t)-\dot{\psi}_c \\ \dot{r}(t) \\ \psi(t)-\psi_c \end{array}\right] }

In the case of M_{\psi}=1, T_{\psi}=0.5T, M_{\delta}=0.5, T_{\delta}=0.1T, we have the following simulation result, which seems to be very nice. Note that there is no feedforward term compared with PID control.

However, in the case of having the time lag L=9, the CLPS is unstable as follows:

LPV control:

Assuming the following speed variation, the responses of CLPS by unity feedback vary largely.

Consider the following NOMOTO model based on a nominal forward speed U^*.

\displaystyle{ \left\{\begin{array}{l} \dot{\psi}(t)=r(t) \\ \displaystyle{\dot{r}(t)=-\left(\frac{U}{U^*}\right)\frac{1}{T^*}r(t) +\left(\frac{U}{U^*}\right)^2\frac{K^*}{T^*}\delta(t)} \end{array}\right. }

Adding the rudder dynamics,

\displaystyle{ \dot{\delta}(t)=-\frac{1}{T_a}\delta(t)+\frac{K_a}{T_a}u(t) }

we have the following state-space representation.

\displaystyle{ \underbrace{ \left[\begin{array}{c} \dot{\psi}(t) \\ \dot{r}(t) \\ \dot{\delta}(t) \end{array}\right] }_{\dot{x}(t)} = \underbrace{ \left[\begin{array}{ccc} 0 & 1 & 0\\ 0 & -\left(\frac{U}{U^*}\right)\frac{1}{T^*} & \left(\frac{U}{U^*}\right)^2\frac{K^*}{T^*} \\ 0 & 0 & -\frac{1}{T_a} \end{array}\right] }_{A(U,U^2)} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ \delta(t) \end{array}\right] }_{x(t)} + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \\ \frac{K_a}{T_a} \end{array}\right] }_{B} u(t) }
\displaystyle{ \underbrace{ \psi(t) }_{y(t)} = \underbrace{ \left[\begin{array}{ccc} 1 & 0 & 0 \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} \psi(t) \\ r(t) \\ \delta(t) \end{array}\right] }_{x(t)} }

Note that it is possible to have the following polytopic representation.

\displaystyle{ A(U,U^2)=p_1(U,U^2)A_1+p_2(U,U^2)A_2+p_3(U,U^2)A_3 }

where

\displaystyle{A_1=A(U_1,U_1^2)}
\displaystyle{A_2=A(U_2,U_2^2)}
\displaystyle{A_3=A(\frac{U_1+U_2}{2},U_1U_2)} 

Defining U_3=\frac{U_1+U_2}{2},

\displaystyle{ p_1(U,U^2)=\frac{1}{p_0}\det \left[\begin{array}{cc} U-U_3 & U_2-U_3 \\ U^2-U_1U_2 & U_2^2-U_1U_2 \\ \end{array}\right]}
\displaystyle{ p_2(U,U^2)=\frac{1}{p_0}\det \left[\begin{array}{cc} U_1-U_3 & U-U_3 \\ U_1^2-U_1U_2 & U^2-U_1U_2 \\ \end{array}\right]}
\displaystyle{ p_3(U,U^2)=\frac{1}{p_0}\det \left[\begin{array}{cc} U_1-U_2 & U_2-U \\ U_1^2-U_2^2 & U_2^2-U^2 \\ \end{array}\right]}
\displaystyle{ p_0=\det \left[\begin{array}{cc} U_1-U_2 & U_2-U_3 \\ U_1^2-U_2^2 & U_2^2-U_1U_2 \\ \end{array}\right]}
\displaystyle{p_1(U,U^2)+p_2(U,U^2)+p_3(U,U^2)=1}

Based on this model, the framework of LPV control is proposed. LPV control can cover the variation included in the triangle region with the vertexes (U_1,U_1^2), (U_2,U_2^2) and temporal vertex (U_3,U_1U_2).

Constitute the following 2-port system.

\displaystyle{ P:  \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x} \\ \dot{x}_I \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A(U,U^2)& 0 \\ -C & 0 \end{array}\right] }_{\cal{A}(U,U^2)} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 1 \end{array}\right] }_{B_1} r + \underbrace{\left[\begin{array}{c} B \\ 0 \end{array}\right] }_{B_2} u\\ \underbrace{ \left[\begin{array}{c} y_{11} \\ y_{12} \end{array}\right] }_{y_1} = \underbrace{ \left[\begin{array}{cc} 0 &\omega_I\\ \omega_DCA(U,U^2) & 0 \end{array}\right] }_{C_1} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \end{array}\right] }_{D_{11}} r + \underbrace{ \left[\begin{array}{c} 0 \\ \omega_DCB \end{array}\right] }_{D_{12}} u\\ \underbrace{ \left[\begin{array}{c} y \\ x_I \end{array}\right] }_{y_2} = \underbrace{ \left[\begin{array}{cc} C & 0\\ 0 & 1  \end{array}\right] }_{C_2} \left[\begin{array}{c} x \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{c} 0 \\ 0 \end{array}\right] }_{D_{21}} r \end{array}\right.}

A(U,U^2)=p_1(U,U^2)A_1+p_2(U,U^2)A_2+p_3(U,U^2)A_3

For this LPV model, design the following LPV controller.

\displaystyle{ K_0: \left\{\begin{array}{l} \dot{x}_K=A_K(U,U^2)x_K+ \underbrace{ \left[\begin{array}{cc} B_K^{(1)}(U,U^2) & B_K^{(2)}(U,U^2) \end{array}\right] }_{B_K(U,U^2)} \left[\begin{array}{c} y \\ x_I \end{array}\right] \\ u=C_K(U,U^2)x_K + \underbrace{ \left[\begin{array}{cc} D_K^{(1)}(U,U^2) & D_K^{(2)}(U,U^2) \end{array}\right] }_{D_K(U,U^2)} \left[\begin{array}{c} y \\ x_I \end{array}\right] \end{array}\right.}

A_K(U,U^2)=p_1(U,U^2)A_{K1}+p_2(U,U^2)A_{K2}+p_3(U,U^2)A_{K3}
B_K(U,U^2)=p_1(U,U^2)B_{K1}+p_2(U,U^2)B_{K2}+p_3(U,U^2)B_{K3}
C_K(U,U^2)=p_1(U,U^2)C_{K1}+p_2(U,U^2)C_{K2}+p_3(U,U^2)C_{K3}
D_K(U,U^2)=p_1(U,U^2)D_{K1}+p_2(U,U^2)D_{K2}+p_3(U,U^2)D_{K3}

Form this, the following output feedback is obtained.

\displaystyle{ K: \left\{\begin{array}{l} \left[\begin{array}{c} \dot{x}_K \\ \dot{x}_I \end{array}\right]= \underbrace{ \left[\begin{array}{cc} A_K(U,U^2) & B_K^{(2)}(U,U^2) \\ 0 & 0 \end{array}\right] }_{A_K(U,U^2)} \left[\begin{array}{c} x_K \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{cc} B_K^{(1)}(U,U^2) & 0\\ -1& 1 \end{array}\right] }_{B_K(U,U^2)} \left[\begin{array}{c} y \\ r \end{array}\right] \\ u= \underbrace{ \left[\begin{array}{cc} C_K(U,U^2) & D_K^{(2)}(U,U^2) \end{array}\right] }_{C_K(U,U^2)} \left[\begin{array}{c} x_K \\ x_I \end{array}\right] + \underbrace{ \left[\begin{array}{cc} D_K^{(1)}(U,U^2) & 0 \end{array}\right] }_{D_K(U,U^2)} \left[\begin{array}{c} y \\ r \end{array}\right] \end{array}\right. }

We can obtain H_\infty control by fixing as U=U^*. The following simulation shows the CLPS responses by H_\infty control. It is observed that the responses are slightly varied according to the three kinds of speed variation.


On the other hand, the following simulation shows the CLPS responses by the LPV control. It is observed that the responses are almost same in spite of the three kinds of speed variation.



LQG Control

1. Optimal Control for MIMO system

\bullet Given an nth-order system (satisfying controllability and observability)

(1)   \begin{eqnarray*} &&\dot{x}(t)=Ax(t)+Bu(t),\ y(t)=Cx(t)\ \\ &&(x(t)\in{\rm\bf R}^n,\ u(t)\in{\rm\bf R}^m,\ y(t)\in{\rm\bf R}^p) \end{eqnarray*}

and the stabilizing state feedback

(2)   \begin{eqnarray*} u(t)=-Fx(t), \end{eqnarray*}

the closed-loop system is represented by

(3)   \begin{eqnarray*} &&\dot{x}(t)=\underbrace{(A-BF)}_{A_F}x(t),\ y(t)=Cx(t)\\ &&A_F=A-BF:\ stable\ matrix. \end{eqnarray*}

The behaviors of the state and the input are given by

(4)   \begin{eqnarray*} x(t)=\exp(A_Ft)x(0) \end{eqnarray*}

and

(5)   \begin{eqnarray*} u(t)=-F\exp(A_Ft)x(0). \end{eqnarray*}

respectively. Then consider a problem to determine F to minimize the criterion function

(6)   \begin{eqnarray*} J=\int_0^\infty(y^T(t)Qy(t)+u^T(t)Ru(t))\,dt. \end{eqnarray*}

\bullet The criterion function can be written as

(7)   \begin{eqnarray*} J&=&\int_0^\infty x^T(t)(C^TQC+F^TRF)x(t)\,dt\\ &=&x^T(0)\underbrace{\int_0^\infty \exp(A_F^Tt)(C^TQC+F^TRF)\exp(A_Ft)\,dt}_{\Pi}x(0)\\ &=&{\rm tr}\ \Pi x(0)x^T(0). \end{eqnarray*}

Note that we have the following constraint on \Pi:

(8)   \begin{eqnarray*} \Pi A_F+A_F^T\Pi+C^TQC+F^TRF=0. \end{eqnarray*}

It is known that \Pi>0 holds because of the closed-loop stability. Taking account of all zero-input responses, instead of (7), we consider to minimize

(9)   \begin{eqnarray*} J={\rm tr}\ \Pi. \end{eqnarray*}

Therefore, we will minimize

(10)   \begin{eqnarray*} J'={\rm tr}\ \Pi+{\rm tr}\ ((\Pi A_F+A_F^T\Pi+C^TQC+F^TRF)\Gamma) \end{eqnarray*}

using undetermined multiplier \Gamma and the stability constraint (8). As the necessary conditions, we have

(11)   \begin{eqnarray*} \frac{\partial J'}{\partial\Pi}=I_n+\Gamma A_F^T+A_F\Gamma=0\Rightarrow \Gamma>0, \end{eqnarray*}

(12)   \begin{eqnarray*} \frac{\partial J'}{\partial F}=2(RF-B^T\Pi)\Gamma=0\Rightarrow F=R^{-1}B^T\Pi, \end{eqnarray*}

(13)   \begin{eqnarray*} \frac{\partial J'}{\partial\Gamma}=\Pi A_F+A_F^T\Pi+C^TQC+F^TRF=0. \end{eqnarray*}

Substituting F=R^{-1}B^T\Pi into (13),

(14)   \begin{eqnarray*} \Pi (A-BR^{-1}B^T\Pi)+(A-BR^{-1}B^T\Pi)^T\Pi\\ +C^TQC+(R^{-1}B^T\Pi)^TRR^{-1}B^T\Pi=0 \end{eqnarray*}

That is, we have the Algebraic Riccati Equation (ARE) on \Pi:

(15)   \begin{eqnarray*} {\Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0}. \end{eqnarray*}

By solving this, we can obtain \Pi>0 and then F by

(16)   \begin{eqnarray*} {F=R^{-1}B^T\Pi}. \end{eqnarray*}

The control methodology is called as LQ Control.

\bullet The sufficiency is discussed as follows. The following expression can be derived:

(17)   \begin{eqnarray*} y^TQy+u^TRu=(u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x) -\underbrace{\frac{d}{dt}x^T\Pi x}_{2x^T\Pi\dot{x}}. \end{eqnarray*}

In fact, substituting \dot{x}=Ax+Bu into the right hand side and using the Riccati equation,

(18)   \begin{eqnarray*} &&{\rm RHS}=u^TRu+2x^T\Pi Bu+x^T\Pi BR^{-1}B^T\Pi x-2x^T\Pi(Ax+Bu)\\ &&=u^TRu+x^T(\Pi A+A^T\Pi+C^TQC)x-2x^T\Pi Ax={\rm LHS}. \end{eqnarray*}

Integrating the both side,

(19)   \begin{eqnarray*} J=\int_0^\infty (u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x)\,dt-\left[x^T\Pi x\right]_0^\infty. \end{eqnarray*}

As x(t)=\exp(A_Ft)x(0)\rightarrow 0\ (t\rightarrow\infty),

(20)   \begin{eqnarray*} J=\int_0^\infty (u+R^{-1}B^T\Pi x)^TR(u+R^{-1}B^T\Pi x)\,dt+x^T(0)\Pi x(0). \end{eqnarray*}

Therefore, it is shown that u=-R^{-1}B^T\Pi x minimizes the criterion function.

Exercise 1
Given A=\left[\begin{array}{cc} 0 & 1 \\ 0 & 0  \end{array}\right],\  B=\left[\begin{array}{cc} 0 \\ 1  \end{array}\right],\  C=\left[\begin{array}{cc} 1 & 0 \\ 0 & 1  \end{array}\right], Q=\left[\begin{array}{cc} 1 & 0 \\ 0 & 1  \end{array}\right],R=1, obtain the solution \Pi= \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \\ \end{array}\right]>0\Leftrightarrow \pi_1>0,\ \pi_1\pi_2-\pi_3^2>0 of the Riccati equation \Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0 and calculate F=R^{-1}B^T\Pi.

Expanding the Ricatti equation

(21)   \begin{eqnarray*} && \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \end{array}\right] + \left[\begin{array}{cc} 0 & 0 \\ 1 & 0 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \nonumber\\ &&-\left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] \left[\begin{array}{c} 0 \\ 1 \end{array}\right] \left[\begin{array}{cc} 0 & 1 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right] + \left[\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right] = \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \end{array}\right] \nonumber , \end{eqnarray*}

we have

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

There are the two solutions on \pi_3 from -\pi_3^2+1=0, the two solutions on \pi_2 from 2\pi_3-\pi_2^2+1=0, the one solution \pi_1 from \pi_1-\pi_2\pi_3=0. Therefore, we have the four kinds of solution (\pi_1,\pi_2,\pi_3) as follows:

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

The only (1) satisfies \pi_1,\pi_2>0,\ \pi_1\pi_2-\pi_3^2>0. Therefore we have

(24)   \begin{eqnarray*} F= \left[\begin{array}{cc} 0 & 1 \end{array}\right] \left[\begin{array}{cc} \pi_1 & \pi_3 \\ \pi_3 & \pi_2 \end{array}\right]= \left[\begin{array}{cc} \pi_3 & \pi_2 \end{array}\right]= \left[\begin{array}{cc} 1 & \sqrt{3} \end{array}\right] . \end{eqnarray*}

How to solve ARE
Given the Riccati equation \Pi A+A^T\Pi-\Pi BR^{-1}B^T\Pi+C^TQC=0, consider the Hamilton matrix

(25)   \begin{eqnarray*} {M=\left[\begin{array}{cc} A & -BR^{-1}B^T \\ C^TQC & -A^T \end{array}\right]}. \end{eqnarray*}

The eigenvalues of the Hamilton matrix with size n\times n are distributed symmetrically to not only the real axis but also the imaginal axis. So there are n stable eigenvalues and n unstable eigenvalues. The eigenvectors corresponding to n stable eigenvalues \lambda_1,\cdots,\lambda_n are obtained as follows:

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

Based on this, we can obtain \Pi as

(27)   \begin{eqnarray*} \Pi=V_2V_1^{-1}. \end{eqnarray*}

A program by SCILAB to solve the Riccati equation is given as follows.

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

Solve the Riccati equation in Exercise 1 by using the above program.

A=[0 1;0 0]; B=[0;1]; C=eye(2,2);
Q=diag([1 1]); R=1;
[F,p]=opt(A,B,C,Q,R)
poles=spec(A-B*F)

Exercise 2
Consider the following spring-connected carts as a control object.

The motion is governed by

(28)   \begin{eqnarray*} \left\{\begin{array}{l} m_1\ddot{x}_1(t)=-k(x_1(t)-x_2(t))+f_1(t) \\ m_2\ddot{x}_2(t)=-k(x_2(t)-x_1(t))+f_2(t) \end{array}\right. , \end{eqnarray*}

where k is a spring constant with the range 0.5<k<2. The state equation and output equation are given by

(29)   \begin{eqnarray*} \underbrace{ \left[\begin{array}{c} \dot{x}_1(t) \\ \dot{x}_2(t) \\ \ddot{x}_1(t) \\ \ddot{x}_2(t) \\ \end{array}\right] }_{\dot{x}(t)} &=& \underbrace{ \left[\begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -\frac{k}{m_1} & -\frac{k}{m_1} & 0 & 0 \\ \frac{k}{m_2} & -\frac{k}{m_2}& 0 & 0   \end{array}\right] }_{A} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \dot{x}_1(t) \\ \dot{x}_2(t) \\ \end{array}\right] }_{x(t)}\\ &&+ \underbrace{ \left[\begin{array}{cc} 0 & 0 \\ 0 & 0 \\ \frac{k}{m_1} &  0 \\ 0 & \frac{k}{m_2}  \end{array}\right] }_{B} \underbrace{ \left[\begin{array}{c} f_1(t) \\ f_2(t) \end{array}\right] }_{u(t)}, \end{eqnarray*}

(30)   \begin{eqnarray*} \underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] }_{y(t)} &=&- \underbrace{ \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0  \end{array}\right] }_{C} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \dot{x}_1(t) \\ \dot{x}_2(t) \\ \end{array}\right] }_{x(t)}. \end{eqnarray*}

The control purpose is to regulate the zero-input response under the initial condition:

(31)   \begin{eqnarray*} \left[\begin{array}{c} {x}_1(0) \\ {x}_2(0) \\ \dot{x}_1(0) \\ \dot{x}_2(0) \\ \end{array}\right] = \left[\begin{array}{c} 0 \\ 0 \\ 0 \\ \frac{k}{m_2}  \end{array}\right] \end{eqnarray*}

by using the state feedback:

(32)   \begin{eqnarray*} \underbrace{ \left[\begin{array}{c} f_1(t) \\ f_2(t) \end{array}\right] }_{u(t)} &=&- \underbrace{ \left[\begin{array}{cccc} g_{11} & g_{12} & g_{13} & g_{14} \\ g_{21} & g_{22} & g_{23} & g_{24}  \end{array}\right] }_{F} \underbrace{ \left[\begin{array}{c} {x}_1(t) \\ {x}_2(t) \\ \dot{x}_1(t) \\ \dot{x}_2(t) \\ \end{array}\right] }_{x(t)}. \end{eqnarray*}

In order to determine the state feedback gain F, for k fixed to an appropriate nominal value k^*, we will minimize the following criterion function:

(33)   \begin{eqnarray*} J=\int_0^\infty (q_1^2x_1^2(t)+q_2^2x_2^2(t)+r_1^2f_1^2(t)+r_2^2f_2^2(t))\,dt, \end{eqnarray*}

that is

(34)   \begin{eqnarray*} J&=&\int_0^\infty ( \underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right]^T }_{y^T(t)} \underbrace{ \left[\begin{array}{cc} q_1^2 & 0 \\ 0 & q_2^2 \ \end{array}\right] }_{Q} \underbrace{ \left[\begin{array}{c} x_1(t) \\ x_2(t) \end{array}\right] }_{y(t)}\\ &&+\underbrace{ \left[\begin{array}{c} f_1(t) \\ f_2(t) \end{array}\right]^T }_{u^T(t)} \underbrace{ \left[\begin{array}{cc} r_1^2 & 0 \\ 0 & r_2^2 \ \end{array}\right] }_{R} \underbrace{ \left[\begin{array}{c} f_1(t) \\ f_2(t) \end{array}\right] }_{u(t)}  )\,dt. \end{eqnarray*}

For example, the closed-loop zero-input response is simulated under the following assumptions:

1) m_1=1, m_2=1, k^*=1,
2) q_1=\frac{1}{M_{x_1}}=1, q_2=\frac{1}{M_{x_2}}=1, r_1=\frac{\rho}{M_{f_1}}=1, r_2=\frac{\rho}{M_{f_2}}=1
(M_{x_1}=M_{x_2}=1, M_{f_1}=M_{f_2}=1,\ \rho=1).

Appendix 1

Check the following properties on matrix trace.

(35)   \begin{eqnarray*} {\rm tr}AB={\rm tr}BA, \end{eqnarray*}

(36)   \begin{eqnarray*} \frac{\partial}{\partial X}{\rm tr}AXB=A^TB^T, \end{eqnarray*}

(37)   \begin{eqnarray*} \frac{\partial}{\partial X}{\rm tr}AX^TB=BA, \end{eqnarray*}

(38)   \begin{eqnarray*} \frac{\partial}{\partial X}{\rm tr}AXBX^T=A^TXB^T+AXB \end{eqnarray*}

where \displaystyle{{\rm tr}A=\sum_{i}a_{ii}} for A=[a_{ij}]_{i,j=1,\cdots,n}. In fact,

(39)   \begin{eqnarray*} \sum_{i}[AB]_{ii}=\sum_{i}\sum_{j}a_{ij}b_{ji} =\sum_{j}\sum_{i}b_{ji}a_{ij}=\sum_{j}[BA]_{jj}, \end{eqnarray*}

(40)   \begin{eqnarray*} &&[\frac{\partial}{\partial X}{\rm tr}AXB]_{ij} =\frac{\partial}{\partial x_{ij}}\sum_{k}[AXB]_{kk} =\frac{\partial}{\partial x_{ij}}\sum_{k}\sum_{i,j}a_{ki}x_{ij}b_{jk}\\ &&=\sum_{k}b_{jk}a_{ki}=[BA]_{ji}=[A^TB^T]_{ij}, \end{eqnarray*}

(41)   \begin{eqnarray*} &&[\frac{\partial}{\partial X}{\rm tr}AX^TB]_{ij} =\frac{\partial}{\partial x_{ij}}\sum_{k}[AX^TB]_{kk} =\frac{\partial}{\partial x_{ij}}\sum_{k}\sum_{i,j}a_{ki}x_{ji}b_{jk}\\ &&=\sum_{k}b_{ik}a_{kj}=[BA]_{ij}, \end{eqnarray*}

(42)   \begin{eqnarray*} &&\frac{\partial}{\partial X}{\rm tr}AXBX^T =\frac{\partial}{\partial X}{\rm tr}(AXB)X^T +\frac{\partial}{\partial X}{\rm tr}AX(BX^T)\\ &&=AXB+A^TXB^T. \end{eqnarray*}