１次系の時間応答

投稿日時: 投稿者:

１次系の時間応答…Homework

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

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

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

\displaystyle{(2)\quad y(t)=ce^{at}x(0)+\int_0^tce^{a(t-\tau)}bu(\tau)d\tau}

を得ます（Note A02-1参照）。ここで、インパルス応答と呼ばれる

\displaystyle{(3)\quad \boxed{g(t)=ce^{at}b}}

を定義しますと、上式は

\displaystyle{(4)\quad \boxed{y(t)=\underbrace{ce^{at}x(0)}_{zero-input\ response} +\underbrace{\int_0^tg(t-\tau)u(\tau)d\tau}_{zero-state\ response}} }

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

●実際、a=-1,b=1,c=1、x(0)=1、u(t)=\sin(10t)についてy(t)を求めてみると次のグラフが得られます。


図1 どれが零状態応答か？

[2] (1)への入力uとして次図の時間関数を考えます。


図2 ステップ入力とインパルス入力(T_s\rightarrow0)

左図においてT_s\rightarrow0としたものがステップ入力、

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

です。これを(2)に代入すると、

\displaystyle{(6)\quad y(t)=\int_0^tg(t-\tau)\cdot 1 d\tau=\int_{t}^0g(t')(-dt')=\int_0^tg(t')dt'=\int_0^tg(\tau)d\tau}

より、ステップ応答が次のように表されます。

\displaystyle{(7)\quad \boxed{y(t)=\int_0^tg(\tau)d\tau}}

すなわち、インパルス応答を積分してステップ応答が得られます。逆に、ステップ応答を微分してインパルス応答が得られることも分かります。

図2の右図は左図の微分\dot{u}を表していますが、T_s\rightarrow0としたものがインパルス入力（デルタ関数）です。実は(3)で定義したインパルス応答はこのインパルス入力に対する応答です。このことを理解するために、いま１次系の状態方程式を微分して

\displaystyle{(8)\quad \ddot{x}(t)=a\dot{x}(t)+b\dot{u}(t)}

を得ます。これは１次系に入力\dot{u}を与えた場合の応答は\dot{x}(t)であることを示しています。図2の左図でT_s\rightarrow0としたステップ入力に対する応答がステップ応答です。図2の右図でT_s\rightarrow0としたインパルス入力はステップ入力の微分ですから、その応答はステップ応答の微分として得られ、インパルス応答となります。すなわち、ステップ応答を微分するとインパルス応答が得られることが分かります。

[3] 特別な３種類の入力（ステップ入力、インパルス入力、正弦波入力）に対する零状態応答（ステップ応答、インパルス応答、正弦波応答）を考え、これらの相互関係を調べます。


図3 １次系の時間応答の相互関係

一般に、「システムを知る」ということは「どんな入力を与えられても出力を推測できる」ことを意味します。したがって、平衡状態（x(0)=0）にある場合、インパルス応答さえ分かれば、畳み込み積分を行って出力（零状態応答）が計算できます。この意味で、インパルス応答は重要なのです。また、インパルス応答自体の計算は、x(0)=bとしたときの零入力応答を求めて行います。ステップ応答が分かれば、これを微分してインパルス応答が分かります。また、インパルス応答をラプラス変換したものが、正弦波応答と密接に関係しています。このことはあとで詳しく述べます。

[4] 次のようにパラメータライズされた漸近安定な１次系を考えます。T、Kはそれぞれ時定数、定常ゲインと呼ばれます。

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

このインパルス応答は次式で与えられます。

\displaystyle{(10)\quad g(t)=\frac{K}{T}e^{-\frac{1}{T}t}}

これを積分して

\displaystyle{(11)\quad y(t)=\frac{K}{T}\int_0^te^{-\frac{1}{T}\tau}d\tau=\frac{K}{T}\left[\frac{e^{-\frac{1}{T}\tau}}{-\frac{1}{T}}\right]_0^t =-K(e^{-\frac{1}{T}t}-1)=K(1-e^{-\frac{1}{T}t}) }

すなわちステップ応答は次式で表されます。

\displaystyle{(12)\quad \boxed{y(t)=K(1-e^{-\frac{1}{T}t})}}

ここで、t=Tのときは

\displaystyle{(13)\quad y(T)=K(1-\frac{1}{e})=0.632K}

となります。また、t=0での接線は

\displaystyle{(14)\quad y(t)=\dot{x}(0)t=\frac{K}{T}t}

ですから、x(T)=Kとなります。この時定数Tは、最終値Kの63.2%まで到達する時間、または初期速度で到達した場合の時間（初期時刻における接線がKと交わる時刻）として特徴付けられます。

●実際、T=0.5,1,2、K=1について、ステップ応答x(t)を求めてみると次のグラフが得られます。最終値K=1に到達する時間は大体5Tであることが分かります。


図4 １次系のステップ応答から時定数を求める補助線とは？

[5] 時定数は平衡状態に復帰する時間の目安にもなります。漸近安定な１次自由系を次のようにパラメータライズします。

\displaystyle{(15)\quad \dot{x}(t)=-\frac{1}{T}x(t)}

この解は

\displaystyle{(16)\quad x(t)=e^{-\frac{1}{T}t}x(0)}

ですから、t=Tのときは

\displaystyle{(17)\quad x(T)=\frac{1}{e}x(0)=0.367x(0)}

となります。また、t=0での接線は

\displaystyle{(18)\quad x(t)=\dot{x}(0)t+x(0)=-\frac{1}{T}x(0)t+x(0)}

ですから、x(T)=0となります。この時定数Tは、初期値の36.7%まで復帰する時間、または初期速度で復帰した場合の時間（初期時刻における接線が0と交わる時刻）として特徴付けられます。

●実際、時刻t=0において平衡状態x=0が乱されx(0)=1になったとし、T=0.5,1,2についてx(t)を求めてみると次のグラフが得られます。平衡状態に復帰する時間は大体5Tであることを表します。

<
図5 １次自由系において平衡状態に復帰する時間（時定数）を定める補助線とは？

演習A02-1…Flipped Classroom

1^\circ 図1の時間応答を描くプログラムを作成せよ。

MATLAB 

%a01_1.m
clear all, close all
a=-1; b=1; c=1;
sys=ss(a,b,c,[]);
t=0:0.01:10;
u=sin(10*t);
lsim(sys,u,t,1), hold on
lsim(sys,u,t,0)
lsim(sys,u*0,t,1)
grid
title('Time Response of 1st-order System')
xlabel('time')
ylabel('x(t)')
legend('total resp','zero-state resp','zero-input resp')
%eof
SCILAB 
coming soon

2^\circ 図4のステップ応答#1を描き、適当な補助線を引いて、交点の座標から時定数を求めるプログラムを作成せよ。

MATLAB 
%a02_2.m
clear all, close all
T1=0.5; T2=1; T3=2; K=1;
a1=-1/T1; a2=-1/T2; a3=-1/T3;
b1=K/T1; b2=K/T2; b3=K/T3; 
c=1;
sys1=ss(a1,b1,c,[]);
sys2=ss(a2,b2,c,[]);
sys3=ss(a3,b3,c,[]);
t=0:0.01:5;
step(sys1,sys2,sys3,t)
grid
hold on
plot([0,5],(1-1/exp(1))*[1 1])
T=ginput(3)
title('Step Responses of 1st-order System')
xlabel('time')
ylabel('x(t)')
legend('T1=0.5','T2=1','T3=2','1-1/e')
%eof
SCILAB 
coming soon

[6] １次系(1)をラプラス変換（Note A02-3参照）して、x(0)=0とすると（零状態応答を考えているので）

\displaystyle{(19)\quad s\hat{x}(s)-\underbrace{x(0)}_0=a\hat{x}(s)+b\hat{u}(s) \Rightarrow \hat{x}(s)=(s-a)^{-1}b\hat{u}(s) }

\displaystyle{(20)\quad \hat{y}(s)=c\hat{x}(s) \Rightarrow \hat{y}(s)=\underbrace{c(s-a)^{-1}b}_{\hat{G}(s)}\hat{u}(s) }

ここで現れた

\displaystyle{(21)\quad \boxed{\hat{g}(s)=c(s-a)^{-1}b}}

を伝達関数と呼びます。伝達関数G(s)にs=j\omegaを代入したG(j\omega)を周波数伝達関数と呼び、正弦波入力に対する時間応答、すなわち周波数応答を特徴づけることができます。

●漸近安定な１次系(9)の伝達関数は次式で与えられます。

\displaystyle{(22)\quad \hat{g}(s)=\frac{K}{Ts+1} }

演習A02-2...Flipped Classroom

1^\circ １次系のステップ応答を次式を逆ラプラス変換して求めなさい。

\displaystyle{(23)\quad \hat{y}(s)=\underbrace{\frac{K}{Ts+1}}_{\hat{g}(s)}\frac{1}{s}=K(\frac{1}{s}-\frac{T}{Ts+1}) }

2^\circ 次のような無駄時間と呼ばれる特殊な入出力関係を考えます。

\displaystyle{(24)\quad y(t)=u(t-L)\quad(L>0) }

これをラプラス変換して無駄時間の伝達関数は次式をとなることを示しなさい。

\displaystyle{(25)\quad \boxed{\hat{g}(s)=e^{-sL}} }

ヒント：

(26)\quad \begin{array}{l} \displaystyle{\hat{y}(s)=\int_0^\infty u(t-L)e^{-st}\,dt =\int_{-L}^\infty u(\tau)e^{-s(\tau+L)}\,d\tau}\\ \displaystyle{=\underbrace{e^{-sL}}_{\hat{g}(s)}\underbrace{\int_0^\infty u(\tau)e^{-s\tau}\,d\tau}_{\hat{u}(s)}} \end{array}

[7] 漸近安定な１次系(9)を考えます。正弦波入力に対する零状態応答すなわち周波数応答は、Note A02-2を参照して

\displaystyle{(27)\quad y(t)=\frac{K}{\sqrt{1+(\omega T)^2}}(\sin(\omega t-\tan^{-1}\omega T)+e^{-\frac{1}{T}t}\sin(\tan^{-1}\omega T))}

ここで、t\rightarrow\inftyとすると

\displaystyle{(28)\quad \boxed{y(t)\simeq\underbrace{\frac{K}{\sqrt{1+(\omega T)^2}}}_{|\hat{g}(j\omega)|} \sin(\omega t\underbrace{-\tan^{-1}\omega T}_{\angle\hat{g}(j\omega)})} }
ただし
\displaystyle{(29)\quad \hat{g}(s)=\frac{K}{Ts+1} }

これは、入力が正弦波のときは、時間が十分立てば、出力も正弦波となることを示しています。その振幅と位相はそれぞれ\hat{g}(j\omega)の絶対値と偏角となっています。

|\hat{g}(j\omega)|をゲイン、\angle\hat{g}(j\omega)を位相と呼びます。このゲイン線図と位相線図をペアにして片対数グラフにプロットしたものをボード線図と呼びます。ゲインはdb値(20{\log_{10}|\hat{g}(j\omega)|)、位相はdeg値(\frac{180}{\pi}\angle\hat{g}(j\omega))です。

●実際、T=0.5,1,2について、ボード線図を描いてみると次のグラフが得られます


図6 １次系のボード線図から時定数を求める補助線とは？

ちなみに

\displaystyle{(30)\quad |\hat{G}(j\omega_b)|=\frac{|\hat{G}(j0)|}{\sqrt{2}}}

を満足する\omega_bは

\displaystyle{(31)\quad \omega_b=\frac{1}{T}}

となり、帯域幅と呼びます。実際、上図で-3dBの水平線との交点が1/Tであることが確かめられます。

演習A02-3...Flipped Classroom

1^\circ 図1のゲイン曲線#1を描き、3dB下がる点の座をマウスをクリックして取得し、時定数を求めるプログラムを作成せよ。

MATLAB 
%a02_3.m
clear all, close all
T1=0.5; T2=1; T3=2; K=1;
a1=-1/T1; a2=-1/T2; a3=-1/T3;
b1=K/T1; b2=K/T2; b3=K/T3; 
c=1;
sys1=ss(a1,b1,c,[]);
sys2=ss(a2,b2,c,[]);
sys3=ss(a3,b3,c,[]);
sys0=ss([],[],[],1/sqrt(2));
w=logspace(-1,1);
bode(sys1,sys2,sys3,sys0,w)
grid
W=ginput(3)
title('Bode Diagrams of 1st-order System')
xlabel('Freq')
ylabel('Gain')
%ylabel('Phase')
legend('T1=0.5','T2=1','T3=2','-3dB')
%eof
SCILAB 
coming soon

Note A02-1　状態方程式の解

●(１次自由系を表す)微分方程式

\displaystyle{(1)\quad \dot{x}(t)=ax(t)}

の解は

\displaystyle{(2)\quad x(t)=e^{at}x(0)}

と表されます。これが解であることは元の微分方程式に代入すればすぐに確かめられ、また次のようにして導くことができます。

元の微分方程式を

\displaystyle{(3)\quad \dot{x}(t)-ax(t)=0}

と書いて、左から積分因数と呼ばれるe^{-at}をかけると

\displaystyle{(4)\quad e^{-at}\dot{x}(t)-ae^{-at}x(t)=0}

すなわち

\displaystyle{(5)\quad \frac{d}{dt}(e^{-at}x(t))=0}

したがって、cを定数として

\displaystyle{(6)\quad e^{-at}x(t)=c\ \Rightarrow\ x(t)=e^{at}c}

ここで、t=0とおくとcは初期値x(0)に等しいので

\displaystyle{(7)\quad \boxed{x(t)=e^{at}x(0)}}

と表されます。

●(１次系の状態方程式を表す)微分方程式

\displaystyle{(8)\quad \dot{x}(t)=ax(t)+bu(t)}

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

\displaystyle{(9)\quad \frac{d}{dt}(e^{-at}x(t))=e^{-at}bu(t)}

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

\displaystyle{(10)\quad \left[e^{-at}x(t)\right]_0^{t'}=\int_0^{t'}e^{-at}bu(t)dt}

すなわち

\displaystyle{(11)\quad e^{-at'}x(t')-x(0)=\int_0^{t'}e^{-at}bu(t)dt}

したがって

\displaystyle{(12)\quad x(t')=e^{at'}x(0)+e^{at'}\int_0^{t'}e^{-at}bu(t)dt}

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

\displaystyle{(13)\quad \boxed{x(t)=e^{at}x(0)+\int_0^te^{a(t-\tau)}bu(\tau)d\tau}}

が得られます。

Note A02-2　正弦波応答

(１次系の状態方程式を表す)微分方程式

\displaystyle{(1)\quad \dot{x}(t)=ax(t)+bu(t) }

に対して、正弦波入力

\displaystyle{(2)\quad u(t)=\sin\omega t}

に対する解を求めます。これは公式

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

を用いて、次のように得られます。

(4)\quad \begin{array}{l} \displaystyle{x(t)=\int_0^te^{a(t-\tau)}b \sin\omega\tau\,d\tau=e^{at}b\int_0^te^{-a\tau} \sin\omega\tau\,d\tau}\\ \displaystyle{=e^{at}b\left[\frac{e^{-a\tau}}{a^2+\omega^2}(-a\sin\omega\tau-\omega\cos\omega\tau)\right]_0^t}\\ \displaystyle{=\frac{e^{at}b}{a^2+\omega^2}(e^{-at}(-a\sin\omega t-\omega\cos\omega t)+\omega)}\\ \displaystyle{=\frac{b}{\sqrt{a^2+\omega^2}}(\sin\omega t\frac{-a}{\sqrt{a^2+\omega^2}}-\cos\omega t\frac{\omega}{\sqrt{a^2+\omega^2}}+e^{+at}\frac{\omega}{\sqrt{a^2+\omega^2}})}\\ \displaystyle{=\frac{b}{\sqrt{a^2+\omega^2}}(\sin\omega t\cos\theta-\cos\omega t\sin\theta+e^{+at}\sin\theta)}\\ \displaystyle{=\frac{b}{\sqrt{a^2+\omega^2}}(\sin(\omega t-\theta)+e^{+at}\sin\theta)\quad (\theta=\tan^{-1}\frac{\omega}{-a})}\\ \displaystyle{=\frac{b}{\sqrt{a^2+\omega^2}}(\sin(\omega t+\phi)-e^{+at}\sin\phi)\quad (\phi=-\theta)} \end{array}

Note A02-3　ラプラス変換

時間関数x(t)のラプラス変換は次式によって定義されます。

\displaystyle{(1)\quad \boxed{\hat{x}(s)=\int_0^\infty x(t)e^{-st}\,dt}}

Note A02-4　dB値 20{\log_{10}(\cdot)