このシステムの時間応答ってどんな感じなんだろか、
このフィルタの周波数特性ってどのようになっているのだろうか、と気になることは多いと思います。
今回は身近なPythonを用いてシステムのステップ応答や周波数特性を求める方法を書いていきます。
・pythonと制御工学の勉強をしている方
・本を読んでもいまいちコードと理論が結びつかない方
モデル
今回は以下のような単純なモデルを例にプログラムを書いていきます。
このモデルの微分方程式は図１に示した通りですが、ここから伝達関数を求めていきます。
まずラプラス変換すると、
$$ms^{2}X(s)+csX(s)+kX(s)=F(s)$$
となります。X(s)でくくると
$$X(s)\{ms^{2}+cs+k\}=F(s)$$
ここで出力/入力の形にします。入力が力で出力が変位ですので、
$$\frac{Output}{Input}=\frac{X(s)}{F(s)}=\frac{1}{ms^{2}+cs+k}$$
以上が今回使うモデルになります。
P制御
P制御は入力と出力の偏差に比例してゲインをかけます。
ブロック線図は図2のようになります。
ここから、$K_p$を変えてどのような応答が出るのかを調べます。
プログラム
今回は以下のパラメータを想定します。
$m=0.01[Kg],c=0.015[N/m/s],k=1[N/m]$
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from control.matlab import *
plt.rcParams["font.size"] = 20
"""
linestyleを付与する関数
yieldを見つけるとジェネレータと判断され、linestyle[lineID]が返される。
]ineIDは最初0だが、「lineID = (lineID+1)%len(linestyle)」の書き方で
1,2,3,0,1,2,,,とループする。
"""
def linestyle_generator():
linestyle = ['-','--','-.',':']
lineID=0
while True:
yield linestyle[lineID]
lineID = (lineID+1)%len(linestyle)
"""
x軸のラベルとy軸のラベル、そして凡例をつけるための関数
"""
def plot_set(fig_ax,*args):
fig_ax.set_xlabel(args[0])
fig_ax.set_ylabel(args[1])
fig_ax.grid(ls=':')
if len(args)==3:
fig_ax.legend(loc=args[2])
k=1
c = 1.5e-2
m = 1.0e-2
P = tf([0,1],[m,c,k])
ref = 30 #deg
Kp = (0.5,1,2)
LS = linestyle_generator()
fig,ax = plt.subplots(figsize=(10,10))
for i in range(3):
K = tf([0,Kp[i]],[0,1])
# K = tf([Kp[i]],[1])でもOK
Gyr = feedback(P*K,1))
y,t = step(Gyr,np.arange(0,10,0.01))
pltargs={'ls':next(LS),'label':'$K_p$'+str(Kp[i])}
ax.plot(t,y*ref,**pltargs)
ax.axhline(ref,color = 'k',linewidth = 0.5)
plot_set(ax,'t[s]','y[m]','best')
これによって時間応答がわかります。
伝達関数が変わったときは
y,t = step(Gyr,np.arange(0,10,0.01))ここのGyrを別の文字にして、その文字にtf([],[])を入れれば任意のステップ応答がわかります。
P制御だけですと定常偏差が生まれてしまいます。
その原理は以下の記事で説明しています。
bode線図
bode線図に含まれる情報量はとてもリッチだそうです。。
伝達関数やシステムのモデルさえわかってしまえば、いつでもbode線図を出力できる！
そう思えたら非常に安心感があります。
プログラム
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from control.matlab import *
plt.rcParams["font.size"] = 20
def linestyle_generator():
linestyle = ['-','--','-.',':']
lineID=0
while True:
yield linestyle[lineID]
lineID = (lineID+1)%len(linestyle)
def plot_set(fig_ax,*args):
fig_ax.set_xlabel(args[0])
fig_ax.set_ylabel(args[1])
fig_ax.grid(ls=':')
if len(args)==3:
fig_ax.legend(loc=args[2])
def bodeplot_set(fig_ax,*args):
# ゲイン線図のグリッドとy軸ラベルの設定
fig_ax[0].grid(which='both',ls=':')
fig_ax[0].set_ylabel('gain[dB]')
# 位相線図のグリッドとx軸y軸ラベルの設定
fig_ax[1].grid(which='both',ls=':')
fig_ax[1].set_xlabel('$\omega$[rad/s]')
fig_ax[1].set_ylabel('Phase[deg]')
if(len(args))>0:
fig_ax[1].legend(loc=args[0]) #引数が一つ以上：ゲイン線図に表示
if len(args)>1:
fig_ax[0].legend(loc=args[1])
LS = linestyle_generator()
fig,ax = plt.subplots(2,1,figsize=(10,10))
k=1
c = 1.5e-2
m = 1.0e-2
P = tf([0,1],[m,c,k])
ref = 30 #deg
Kp = (0.5,1,2)
for i in range(len(Kp)):
K = tf([0,Kp[i]],[0,1])
Gyr = feedback(P*K,1)
gain, phase ,w = bode(Gyr,logspace(-1,3),plot = False)
pltargs = {'ls':next(LS),'label':'$K_P$='+str(Kp[i])}
ax[0].semilogx(w,20*np.log10(gain),**pltargs)
ax[1].semilogx(w,phase*180/np.pi,**pltargs)
bodeplot_set(ax,'lower left')
参考文献にも記しますが「Pythonによる制御工学入門」という本を参考にしています。
このプログラムですと、例えばハイパスフィルタなどの時定数をいくつにしたらいいかを検討する際,
簡単に周波数特性を確認できます。例えば以下のようにWhのタプルと、HPFの伝達関数、そしてループするrange(len(Wh))のところを変えればOKです。
#関数の定義は省略
LS = linestyle_generator()
fig,ax = plt.subplots(2,1,figsize=(10,10))
Wh = (10,50,100,150,200)
for i in range(len(Wh)):
HPF = tf([1,0],[1,Wh[i]])
gain, phase ,w = bode(HPF,logspace(-1,3),plot = False)
pltargs = {'ls':next(LS),'label':'$Wh$='+str(Wh[i])}
ax[0].semilogx(w,20*np.log10(gain),**pltargs)
ax[1].semilogx(w,phase*180/np.pi,**pltargs)
bodeplot_set(ax,'lower left')
#fig.savefig('HPF.PDF') 保存するとき用
参考文献
南さんはYouTubeもアップされており、大学院の講義でわからないところがあった際、
動画を拝見しお世話になりました。
コメント