２次系のステップ応答

２次系のステップ応答…Homework

[0]　次の $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.}$

いま行列 $B$ の第 $j$ 列を $B_j$ 、行列 $C$ の第 $i$ 列を $C_i$ で表すとき、第 $j$ 番目の入力 $u_j$ から第 $i$ 番目の出力 $y_i$ までの時間応答は

$\displaystyle{(2)\quad y_i(t)=\underbrace{C_i\exp(At)x(0)}_{zero-input\ response} +\underbrace{\int_0^tG_{ij}(t-\tau)u_j(\tau)d\tau}_{zero-state\ response}}$

のように表されます。ここで、 $G_{ij}(t)$ は(1)のインパルス応答行列の $(i,j)$ 要素で、次のように与えられます。

$\displaystyle{(3)\quad G_{ij}(t)=C_i\exp(At)B_j}$

このときステップ入力

$\displaystyle{(4)\quad u_j(t)= \left\{\begin{array}{ll} 0 & (t\le0)\\ 1 & (t>0) \end{array}\right.}$

に対する零状態応答はステップ応答と呼ばれ、インパルス応答を積分して

$\displaystyle{(5)\quad y_i(t)=\int_0^tG_{ij}(\tau)d\tau}$

のように計算されます。

[1] さて、次のようにパラメタライズされた漸近安定な１入力１出力２次系の状態空間表現を考えます。

$\displaystyle{(6a)\quad \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{(6b)\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$ とした２次系(7)に対するインパル応答を(6)から求めます。そのために、行列 $A$ の固有値 $\lambda(A)=\{\lambda_1,\lambda_2\}$ を計算すると、 $\zeta>1$ 、 $\zeta<1$ 、 $\zeta=1$ に応じて、それぞれ次に示すようになります。

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

図1 $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{(8a)\quad G(t)=\frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}(e^{\lambda_2t}-e^{\lambda_1t})}$
$\displaystyle{(8b)\quad \boxed{S(t)=1+\frac{1}{\lambda_2-\lambda_1}(\lambda_1e^{\lambda_2t}-\lambda_2e^{\lambda_1t})}}$

実際、インパルス応答を積分して、ステップ応答は次のように計算されます。

$(9)\quad \begin{array}{l} \displaystyle{S(t)=\int_0^tG(\tau)d\tau=\int_0^t\frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}(e^{\lambda_2\tau}-e^{\lambda_1\tau})d\tau}\\ \displaystyle{= \frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}\left[\frac{e^{\lambda_2\tau}}{\lambda_2}-\frac{e^{\lambda_1\tau}}{\lambda_1}\right]_0^t}\\ \displaystyle{=\frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}(\frac{e^{\lambda_2t}}{\lambda_2}-\frac{e^{\lambda_1t}}{\lambda_1})-\frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}(\frac{1}{\lambda_2}-\frac{1}{\lambda_1})}\\ \displaystyle{=1+\frac{1}{\lambda_2-\lambda_1}(\lambda_1e^{\lambda_2t}-\lambda_2e^{\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{(10a)\quad G(t)=\frac{\omega_n^2}{\lambda_I}e^{\lambda_Rt}\sin\lambda_I t}$
$\displaystyle{(10b)\quad \boxed{S(t)=1-\frac{\omega_n}{\lambda_I}e^{\lambda_Rt}\sin(\lambda_I t+\tan^{-1}\frac{\lambda_I}{-\lambda_R})}}$

実際、インパルス応答を積分して、公式

$\displaystyle{(11)\quad \int e^{ax}\sin bx dx=\frac{e^{ax}}{a^2+b^2}(a\sin bx-b \cos bx)}$

を用いて、ステップ応答は次のように計算されます。

$(12)\quad \begin{array}{l} \displaystyle{S(t)=\int_0^tG(\tau)d\tau=\int_0^t\frac{\omega_n^2}{\lambda_I}e^{\lambda_R\tau}\sin\lambda_I \tau d\tau}\\ \displaystyle{=\frac{\omega_n^2}{\lambda_I}\left[ \frac{e^{\lambda_R\tau}}{\lambda_R^2+\lambda_I^2}(\lambda_R\sin\lambda_I\tau-\lambda_I\cos\lambda_I\tau)\right]_0^t}\\ \displaystyle{=\frac{\omega_n^2}{\lambda_I}\frac{1}{\omega_n^2} (e^{\lambda_Rt}(\lambda_R\sin\lambda_It-\lambda_I\cos\lambda_It)+\lambda_I)}\\ \displaystyle{=1-\frac{\omega_n}{\lambda_I}e^{\lambda_Rt} (\sin\lambda_It\times\frac{-\lambda_R}{\omega_n}+\cos\lambda_It\times\frac{\lambda_I}{\omega_n})}\\ \displaystyle{=1-\frac{\omega_n}{\lambda_I}e^{\lambda_Rt} (\sin\lambda_It\cos\phi+\cos\lambda_It\sin\phi)}\\ \displaystyle{=1-\frac{\omega_n}{\lambda_I}e^{\lambda_Rt}\sin(\lambda_It+\phi) \quad(\phi=\tan^{-1}\frac{\lambda_I}{-\lambda_R})} \end{array}$

[5] ２次系のステップ応答(3)

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

$\displaystyle{(13a)\quad G(t)=\lambda^2te^{\lambda t}}$
$\displaystyle{(13b)\quad \boxed{S(t)=1+(\lambda t-1)e^{\lambda t}}}$

実際、インパルス応答を積分して、ステップ応答は次のように計算されます。

$(14)\quad \begin{array}{l} \displaystyle{S(t)=\int_0^tG(\tau)d\tau=\int_0^t\lambda^2\tau e^{\lambda\tau}d\tau}\\ \displaystyle{=\lambda^2\left[\tau \frac{1}{\lambda}e^{\lambda\tau}\right]_0^t -\lambda^2\int_0^t \frac{1}{\lambda}e^{\lambda\tau}d\tau}\\ \displaystyle{=\lambda te^{\lambda t}-\lambda\left[\frac{1}{\lambda}e^{\lambda\tau}\right]_0^t}\\ \displaystyle{=1+(\lambda t-1)e^{\lambda t}} \end{array}$

２次振動系の同定…Homework

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

図2 ２次系のステップ応答の比較

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

$\displaystyle{(15)\quad S(t)=1-\frac{\omega_n}{\lambda_I}e^{\lambda_Rt}\sin(\lambda_I t+\phi)\quad(\phi=\tan^{-1}\frac{\lambda_I}{-\lambda_R})}$

の第１オーバーシュートについて、次が成り立ちます。

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

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

$(17)\quad \begin{array}{l} \displaystyle{\dot{S}(t)=-\frac{\omega_n}{\lambda_I}(\lambda_Re^{\lambda_Rt}\sin(\lambda_It+\phi)+e^{\lambda_Rt}\cos(\lambda_It+\phi)\times\lambda_I)=0}\\ \displaystyle{\Rightarrow \sin(\lambda_It+\phi)\times\frac{-\lambda_R}{\omega_n}-\cos(\lambda_It+\phi)\times\frac{\lambda_I}{\omega_n}=0}\\ \displaystyle{\Rightarrow \sin(\lambda_It+\phi)\cos\phi-\cos(\lambda_It+\phi)\sin\phi=0}\\ \displaystyle{\Rightarrow \sin(\lambda_It+\phi-\phi)=0}\\ \displaystyle{\Rightarrow \sin\lambda_It=0}\\ \displaystyle{\Rightarrow \lambda_IT_p=\pi} \end{array}$

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

$\displaystyle{(18)\quad T_p=\frac{\pi}{\lambda_I}=\frac{\pi}{\omega_n\sqrt{1-\zeta^2}}}$

この時刻 $T_p$ における $S(t)$ の値は

$\displaystyle{(19)\quad S(T_p)=1-\frac{\omega_n}{\lambda_I}e^{\lambda_R\frac{\pi}{\lambda_I}}\sin(\lambda_I\frac{\pi}{\lambda_I}+\phi) =1-\underbrace{\frac{\omega_n}{\lambda_I}}_{\frac{1}{\sin\phi}}e^{\frac{\lambda_R}{\lambda_I}\pi}\underbrace{\sin(\pi+\phi)}_{-\sin\phi}}$

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

$\displaystyle{(20)\quad S(T_p)=1+e^{\frac{\lambda_R}{\lambda_I}\pi}=1+\underbrace{\exp(-\frac{\zeta\pi}{\sqrt{1-\zeta^2}})}_{p_0} }$

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

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

演習A23…Flipped Classroom

図2のステップ応答#1を描き、第１オーバーシュートの頂点の座標 $(T_p,1+p_0)$ をマウスをクリックして取得し、減衰係数と固有角周波数 $(\zeta,\omega_n)$ を求めるプログラムを作成せよ。

MATLAB

%a23.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;
 step(sys1,sys2,sys3,t) 
 grid
 title('Step Resp of 2nd-order System')
 xlabel('time')
 ylabel('y(t)')
 legend('zeta=0.01','zeta=0.707','zeta=1')
%-----
 [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
 hold on
 plot(Tp,1+p0,'o')
 A1a=[0 1;-wn^2 -2*zeta]; 
 sys1a=ss(A1a,B,C,D);
 step(sys1a,t) 
%-----
%eof

SCILAB

制御系CAD

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

２次系のステップ応答