はじめに
あるシステムが存在するとき、入力に対して出力は一般的に遅れる。これはある意味で常識かつ経験則的な考え方とも言える。というのは、仮に出力が入力よりも進んでいたら、現在の出力を検知することができれば未来の入力を検知することができてしまい、未来予知となってしまうからである。したがって、一般的なシステムは『遅れ系』であることが多い。今回はRLC回路のような振動解を持つ『二次遅れ系』のシステムに対して単位ステップを入力した場合、どのような応答をするのかをPythonを用いて考察する。
二次遅れ系
二次遅れ系は、RLC直列回路のような振動する系のモデルであり以下の伝達関数で定義される。
G(s) =\frac{\omega_n^2}{s^2+2\zeta \omega_n +\omega^2_n}
ただし、$\zeta$は減衰係数で$\omega_n$は固有角速度である。
ここで、$0.3<\zeta<1.5$という領域で$\zeta$を動かしたとき、単位ステップ応答はどのように変化するのかを確認する。
部分分数展開とラプラス逆変換
単位ステップ応答Y(s)は以下のように表すことができる。
Y(s)=G(s)U(s)=\frac{\omega_n^2}{s^2+2\zeta \omega_n +\omega^2_n} \frac{1}{s}=\frac{1}{s}-\frac{s^2+2\zeta \omega_n}{s^2+2\zeta \omega_n +\omega^2_n}
となる。
$\frac{s^2+2\zeta \omega_n}{s^2+2\zeta \omega_n +\omega^2_n}$に対して、更に部分分数展開をすることが可能かを考える。（実数のみではなく複素数まで範囲を拡張する）
ここで、分母＝0と考えて分母を因数分解できるかを考えた場合、
\lambda_\pm = -\zeta \omega_n \pm \sqrt{\zeta^2-1}\omega_n
を解として$(s-\lambda_+)(s-\lambda_-)$と分母は因数分解できる。
\frac{s^2+2\zeta \omega_n}{s^2+2\zeta \omega_n +\omega^2_n}=\frac{a}{s-\lambda_+}+\frac{b}{s-\lambda_-}
とおく。（a,bは実数もしくは複素数の定数である）
(分母)=0の二次方程式が重解を持たない場合と持つ場合で場合分けをする。（重解を保つ場合は極限を取ることで重解を持たない場合より推定する方法もあるが、、）
重解を持たない場合
a=\frac{\zeta+ \sqrt{\zeta^2-1}}{2\sqrt{\zeta^2-1}},b=\frac{-\zeta+ \sqrt{\zeta^2-1}}{2\sqrt{\zeta^2-1}}
と与えられる。したがって、出力の時間応答は、
y(t)=1-(a e^{\lambda_+ t}+ b e^{\lambda_- t})
と表すことができる。
重解を持つ場合
Y(s)は、以下のように$\frac{1}{s},\frac{1}{s+\omega},\frac{\omega_n}{(s+\omega_n)^2}$の和で表すことができる。
Y(s)=\frac{\omega_n^2}{(s+\omega_n)^2}\frac{1}{s}=\frac{1}{s}-\frac{1}{s+\omega}+\frac{\omega_n}{(s+\omega_n)^2}
したがって、
y(t)=1-e^{-\omega_n t}(1+\omega_n t)
と表すことができる。
プログラム
上記の考察をもとにして、減衰係数$\zeta$を変化させた場合において$y(t)$の時間応答を調べるプログラムを書いた。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
# 要素分割数
n=1000
# 時刻はt=0からt=t0まで
# 減衰係数はjita_minからjita_maxまで変化するものとする。
t0=0.6
jita_min=0.3
jita_max=1.5
t=np.zeros((n,n))
y=np.zeros((n,n))
jita=np.zeros((n,n))
#print(X)
# 固有角速度
omega_n=60
for i in range(n):
for k in range(n):
t[i][k]= (i/n)*t0
jita[i][k]=jita_min+(jita_max-jita_min)*k/n
# 解が異なる2つの実数解のとき
if jita[i][k] > 1:
a= (jita[i][k]+((jita[i][k])**2-1)**0.5)/(2*((jita[i][k])**2-1)**0.5)
b= (-jita[i][k]+((jita[i][k])**2-1)**0.5)/(2*((jita[i][k])**2-1)**0.5)
y[i][k]= 1-(a* np.exp(-(jita[i][k]-(jita[i][k]**2-1)**0.5)*omega_n*t[i][k]) +b* np.exp(-(jita[i][k]+1j*(jita[i][k]**2-1)**0.5)*omega_n*t[i][k]))
# 解が重解のとき
elif jita[i][k]==1:
y[i][k]=1-(np.exp(-omega_n*t[i][k]))*(1+omega_n*t[i][k])
# 解が異なる2つの虚数解のとき
else:
a= (jita[i][k]+1j*((-jita[i][k])**2+1)**0.5)/(1j*2*((-jita[i][k])**2+1)**0.5)
b= (-jita[i][k]+1j*((-jita[i][k])**2+1)**0.5)/(1j*2*((-jita[i][k])**2+1)**0.5)
y[i][k]= 1-(a* np.exp(-(jita[i][k]-1j*(-jita[i][k]**2+1)**0.5)*omega_n*t[i][k]) +b* np.exp(-(jita[i][k]+1j*(-jita[i][k]**2+1)**0.5)*omega_n*t[i][k]))
plt.title("減衰係数と単位ステップ応答")
plt.xlabel('時間(s)')
plt.ylabel('応答値')
plt.contourf(t,y,jita)
plt.colorbar(label="減衰係数")
plt.savefig("減衰係数と単位ステップ応答.png")
plt.show()
このプログラムを実行すると以下のような画像が出力される。
このように、減衰係数を大きくするほど収束までに振動する強度が小さくなる。
まとめ
今回は、二次遅れ系の単位ステップ応答について部分分数展開を用いたラプラス逆変換を行うことで特性をグラフ化した。具体的には、特性方程式が重解を持たない場合と持つ場合で分けることによって解を表現した。ただし、Pythonでは根号の中身が負になる計算を上手くcomplex型を通して扱うことができなかったので、3パターンに場合分けすることで解いた。
参考文献