ｎ次系の時間応答

ｎ次系の時間応答…Homework

[1] 次の $m$ 入力 $p$ 出力 $n$ 次系の状態空間表現を考えます。

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

これから、 $n$ 次系の時間応答の表現式

$\displaystyle{(2)\quad y(t)=C\exp(At)x(0)+\int_0^tC\exp(A(t-\tau))Bu(\tau)d\tau}$

を得ます。ここで、インパルス応答行列

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

を定義しますと（ $G(t)\in{\rm\bf R}^{p\times m}$ ）、上式は

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

すなわち、 $n$ 次系の時間応答は、零入力応答（第１項）と零状態応答（第２項）の和となります。

いま、(1)をラプラス変換して、 $x(0)=0$ とすると（零状態応答だけを考えます）

$\displaystyle{(5a)\quad \hat{x}(s)-\underbrace{x(0)}_0=A\hat{x}(s)+B\hat{u}(s) \Rightarrow \hat{x}(s)=(sI_n-A)^{-1}B\hat{u}(s) }$

$\displaystyle{(5b)\quad \hat{y}(s)=C\hat{x}(s) \Rightarrow \hat{y}(s)=\underbrace{C(sI_n-A)^{-1}B}_{\hat{G}(s)}\hat{u}(s) }$

ここで現れた

$\displaystyle{(6)\quad \boxed{\hat{G}(s)=C(sI_n-A)^{-1}B}}$

を伝達関数行列と呼びます。

● $n$ 次系の時間応答についての相互関係は次図のように表されます。多変数系の場合、 $i$ 番目の入力から $j$ 番目の出力までのステップ応答と周波数応答を考えます。

図1 多変数系時間等応答の相互関係

２次系の時間応答…Homework

[2] さて、次のようにパラメタライズされた漸近安定な１入力１出力２次系の状態空間表現を考えます。 $\zeta>0$ 、 $\omega_n$ 、 $K$ はそれぞれ減衰係数、固有角周波数、定常ケインと呼ばれます。

$\displaystyle{(7a)\quad \boxed{\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\\ -\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(t)} + \underbrace{ \left[\begin{array}{cc} 0 \\ K\omega_n^2 \end{array}\right] }_{B} u(t)} }$

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

以下では、簡単のため $K=1$ とした２次系に対するインパルス応答(3)を求めます。そのために、行列 $A$ の固有値 $\lambda(A)=\{\lambda_1,\lambda_2\}$ を計算すると、 $\zeta>1$ 、 $\zeta<1$ 、 $\zeta=1$ に応じて、それぞれ次に示すようになります。

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

図2 $A=\left[\begin{array}{cc} 0 & 1\\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right]$ の固有値分布

[3] ２次系のインパルス応答(1)

$\zeta>1$ のとき、 $\lambda_1,\lambda_2=-\zeta\omega_n\pm\omega_n\sqrt{\zeta^2-1}$ とおくと、インパルス応答は次式で与えられます。

$\displaystyle{(9)\quad \boxed{G(t)=\frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}(e^{\lambda_2t}-e^{\lambda_1t})}}$

実際

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

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

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

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

[4] ２次系のインパルス応答(2)

$\zeta<1$ のとき、 $\lambda_1,\lambda_2=\underbrace{-\zeta\omega_n}_{\lambda_R}\pm j\underbrace{\omega_n\sqrt{1-\zeta^2}}_{\lambda_I}$ とおくと、インパルス応答とステップ応答はそれぞれそれぞれ次式で与えられます。

$\displaystyle{(12)\quad \boxed{G(t)=\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\sin\lambda_I t}}$

実際

$\displaystyle{(13a)\quad \underbrace{ \left[\begin{array}{cc} 0 & 1 \\ -\omega_n^2 & -2\zeta\omega_n \end{array}\right] }_{A} = \underbrace{ \left[\begin{array}{cc} 1 & 0 \\ \lambda_R & \lambda_I \end{array}\right] }_{V} \underbrace{ \left[\begin{array}{cc} \lambda_R & \lambda_I \\ -\lambda_I & \lambda_R \end{array}\right] }_{\Lambda} \underbrace{ \frac{1}{\lambda_I} \left[\begin{array}{cc} \lambda_I & 0 \\ -\lambda_R & 1 \end{array}\right] }_{V^{-1}}}$
$\displaystyle{(13b)\quad \exp(\Lambda t)= e^{\lambda_R t} \left[\begin{array}{cc} \cos\lambda_I t & \sin\lambda_I t\\ \sin\lambda_I t & \cos\lambda_I t \end{array}\right]}$

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

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

[5] ２次系のインパルス応答(3)

$\zeta=1$ のとき、 $\lambda_1=\lambda_2=\lambda=-\zeta\omega_n=-\omega_n$ とおくと、インパルス応答とステップ応答はそれぞれそれぞれ次式で与えられます。

$\displaystyle{(15)\quad \boxed{G(t)=\lambda^2te^{\lambda t}}}$

実際

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

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

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

２次振動系の同定…Homework

[7] ２次系(7)において、 $\zeta=0.01,1/\sqrt{2},1$ 、 $\omega_n=1$ 、 $K=1$ の場合、インパルス応答のシミュレーション例を示します。

図3 ２次系のインパルス応答の比較

$\zeta<1$ の場合の(7)で表される２次系を２次振動系と呼びます。そのインパルス応答

$\displaystyle{(18)\quad G(t)=\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\sin\lambda_I t}$

の第k+1番目 $(k=0,1,\cdots)$ の頂点の座標について、次が成り立ちます。

$\displaystyle{(19)\quad \boxed{(t_k,G(t_k))=\left(\frac{k\pi+\phi}{\omega_n\sqrt{1-\zeta^2}},\omega_n\exp(-\frac{\zeta(k\pi+\phi)}{\sqrt{1-\zeta^2}})\right)} }$

実際、 $\dot{G}(t)=0$ を計算すると

$(20)\quad \begin{array}{l} \displaystyle{\dot{G}(t)=\frac{\omega_n^2}{\lambda_I}(\lambda_Re^{\lambda_Rt}\sin\lambda_It+e^{\lambda_Rt}\cos\lambda_It\times\lambda_I)=0}\\ \displaystyle{\Rightarrow \sin\lambda_It\times\frac{-\lambda_R}{\omega_n}-\cos\lambda_It\times\frac{\lambda_I}{\omega_n}=0}\\ \displaystyle{\Rightarrow \sin\lambda_It\cos\phi-\cos\lambda_It\sin\phi=0,\quad\phi=\tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta}}\\ \displaystyle{\Rightarrow \sin(\lambda_It-\phi)=0}\\ \displaystyle{\Rightarrow \lambda_It_k-\phi=k\pi\quad(k=0,1,\cdots)} \end{array}$

となって、次式を得ます。

$\displaystyle{(21)\quad t_k=\frac{k\pi+\phi}{\lambda_I}=\frac{k\pi+\phi}{\omega_n\sqrt{1-\zeta^2}}$

この時刻 $t_k$ における $G(t)$ の値は

$\displaystyle{(22)\quad \begin{array}{l} G(t_k)=\frac{\omega_n^2}{\lambda_I}e^{\lambda_R\frac{k\pi+\phi}{\lambda_I}}\sin(\lambda_I\frac{\phi}{\lambda_I})\\ =\omega_n\underbrace{\frac{\omega_n}{\lambda_I}}_{\frac{1}{\sin\phi}}e^{\frac{\lambda_R}{\lambda_I}(k\pi+\phi)}\underbrace{\sin(k\pi+\phi)}_{(-1)^k\sin\phi}\\=(-1)^k\omega_ne^{\frac{\lambda_R}{\lambda_I}(k\pi+\phi)} \end{array} }$

となって、次式を得ます。

$\displaystyle{(23)\quad G(t_k)=(-1)^k\omega_ne^{\frac{\lambda_R}{\lambda_I}(k\pi+\phi)} =(-1)^k\omega_n\exp(-\frac{\zeta}{\sqrt{1-\zeta^2}}(k\pi+\phi)) }}$

●さて、(19)からそれぞれ次の関係式を得ます。

$\displaystyle{(24)\quad t_{k+2}-t_k=\frac{(k+2)\pi+\phi}{\omega_n\sqrt{1-\zeta^2}}-\frac{k\pi+\phi}{\omega_n\sqrt{1-\zeta^2}} =\frac{2\pi}{\omega_n\sqrt{1-\zeta^2}}=T}$

$\displaystyle{(25)\quad \frac{G(t_{k+2})}{G(t_k)} =\exp(-\frac{\zeta((k+2)\pi+\phi)}{\sqrt{1-\zeta^2}} +\frac{\zeta(k\pi+\phi)}{\sqrt{1-\zeta^2}}) =\exp(-\frac{2\pi\zeta}{\sqrt{1-\zeta^2}})=\Lambda}$

(24)はインパルス応答の隣合う山と山の時間間隔（または谷と谷の時間間隔） $T$ は一定であることを示しています。(25)はインパルス応答の隣合う山と山の振幅比（または谷と谷の振幅比） $\Lambda$ は一定であることを示しています。

したがって、次のように、 $(T,\Lambda)$ から $(\zeta,\omega_n)$ を得ることができます。

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

演習A22…Flipped Classroom

図2の応答#1を描き、２つの山の頂点の座標をマウスをクリックして取得し、減衰係数と固有角周波数を求めるプログラムを作成せよ。

MATLAB

%a22.m
%-----
 clear all, close all
 A1=[0 1;-1 -2*0.01]; 
 A2=[0 1;-1 -2/sqrt(2)]; 
 A3=[0 1;-1 -2*1]; 
 B=[0;1]; C=[1 0]; D=[];
 sys1=ss(A1,B,C,D);
 sys2=ss(A2,B,C,D);
 sys3=ss(A3,B,C,D);
 t=0:0.1:10;
 x0=B;
 initial(sys1,sys2,sys3,x0,t) 
 grid
 title('Impulse Resp of 2nd-order System')
 xlabel('time')
 ylabel('y(t)')
 legend('zeta=0.01','zeta=0.707','zeta=1')
%-----
 [Tp,p0]=ginput(2);
 T=Tp(2)-Tp(1);
 Lambda=p0(2)/p0(1),
 zeta=sqrt(log(Lambda)^2/(log(Lambda)^2+4*pi^2))
 wn=sqrt(log(Lambda)^2+4*pi^2)/T
 hold on
 plot(Tp(1),p0(1),'o',Tp(2),p0(2),'o')
 A1a=[0 1;-wn^2 -2*zeta]; 
 sys1a=ss(A1a,B,C,D);
 initial(sys1a,x0,t) 
%-----
%eof

SCILAB

Note A22-1　状態方程式の解

$n$ 次系の状態方程式

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

の解を求めてみます。Note A21と同様にして

$\displaystyle{(2)\quad\frac{d}{dt}(\exp(-At)x(t))=\exp(-At)Bu(t)}$

までは問題ないと思います。これを $0$ から $t'$ まで積分して

$\displaystyle{(3)\quad\left[\exp(-At)x(t)\right]_0^{t'}=\int_0^{t'}\exp(-At)Bu(t)dt}$

すなわち

$\displaystyle{(4)\quad\exp(-At')x(t')-x(0)=\int_0^{t'}\exp(-At)Bu(t)dt}$

したがって

$\displaystyle{(5)\quad x(t')=\exp(At')x(0)+\exp(At')\int_0^{t'}\exp(-At)Bu(t)dt}$

ここで、 $t$ を $\tau$ に、 $t'$ を $t$ と置き換えれば

$\displaystyle{(6)\quad\boxed{x(t)=\exp(At)x(0)+\int_0^t\exp(A(t-\tau))Bu(\tau)d\tau}}$

が得られます。

Note A22-2　重ね合わせの原理

次の $m$ 入力 $p$ 出力 $n$ 次系の状態空間表現を考えます。

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

２つの入力 $u_1$ と $u_2$ に対するそれぞれの零状態応答は

$\displaystyle{(2)\quad y_1(t)=\int_0^tC\exp(A(t-\tau))Bu_1(\tau)d\tau}$

$\displaystyle{(3)\quad y_2(t)=\int_0^tC\exp(A(t-\tau))Bu_2(\tau)d\tau}$

となります。このとき入力

$\displaystyle{(4)\quad \alpha u_1(t)+\beta u_2(t)}$

に対する零状態応答は

$\displaystyle{(5)\quad \alpha y_1(t)+\beta y_2(t)}$

となります。実際、

$(6)\quad \begin{array}{l} \displaystyle{\int_0^tC\exp(A(t-\tau))B(\alpha u_1(\tau)+\beta u_2(\tau))d\tau}\\ \displaystyle{=\alpha \int_0^tC\exp(A(t-\tau))Bu_1(\tau)d\tau +\beta \int_0^tC\exp(A(t-\tau))Bu_2(\tau))d\tau}\\ =\alpha y_1(t)+\beta y_2(t)} \end{array}$

このとき、重ね合わせの原理が成り立つ、 $n$ 次系(1)は線形系(Linear System)であると言います。重ね合わせの原理が成り立たないときは非線形系(Nonlinear System)であると言います。たとえば入力を２倍したときの時間応答が元の時間応答の２倍にならなければ非線形系であると判断できます。

制御系CAD

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

ｎ次系の時間応答