システムの制御設計で2次系伝達関数のモデルを扱う際にどのような制御特性で設計するか考えることがあるとがあると思います。
そういった際に減衰係数や固有角振動数を変えることでどのような振る舞いをするのか把握しておくと設計を進めやすくなります。
本記事では、2次系伝達関数モデルの減衰係数を変更した際の挙動の違いを、ステップ応答波形とボード線図で比較します。比較のためのシミュレーションは無料で扱える Python-Control のライブラリを用いて計算し、matplotlib で波形を出力します。
Python-Control 基本機能
Python-Control の基本は下記の記事を参照ください。
2次系伝達関数モデル
s領域における2次系の伝達関数は次式になります。
ζ は減衰係数、ωn は固有角振動数です。
$$ G(s) = \frac{\omega_{n}^2}{s + 2 \zeta \omega_{n} s + \omega_{n}^2 } $$
ステップ応答波形描画
まず演算に必要なライブラリをインポートします。
次に2次系モデルを扱いやすくするため tf_2nd_order( ) で関数化します。
なお関数の第一引数は ζ 、第二引数は ωn とします。
ステップ応答を出力する際は numpy の関数 linspace( ) を使って時間パラメータを作成しています。下記のコードでは 0sec から 2sec まで 1000 分割の分解能になります。
そして伝達関数モデルと時間を関数 step( ) の引数に与えてステップ応答の出力を得ます。
# 数値演算ライブラリ
import numpy as np
# 制御工学ライブラリ
from control import matlab
# グラフ描画
from matplotlib import pyplot as plt
# 2次系伝達関数モデル作成の関数
def tf_2nd_order(zeta, wn):
num = [wn*wn]
den = [1, 2*zeta*wn, wn*wn]
G = matlab.tf(num, den)
return G
# 伝達関数モデル作成
G2nd = tf_2nd_order(0.6, 10.0)
# 伝達関数モデル出力
print(G2nd)
# ステップ応答
t = np.linspace(0, 2.0, 1000)
(output, t) = matlab.step(G2nd, t)
# 出力パラメータ
plt.plot(t, output)
plt.axhline(1.0, ls=":", color="blue")
# グラフ表示設定
plt.rcParams['font.family'] = 'Times New Roman' # 全体のフォント
plt.title('2nd order step', fontsize=10) # グラフタイトル
plt.xlabel('time [sec]', fontsize=10) # x軸ラベル
plt.ylabel('Output', fontsize=10) # y軸ラベル
plt.xlim([0, 2.0]) # x軸範囲
plt.ylim([0, 1.2]) # y軸範囲
plt.tick_params(labelsize = 10) # 軸ラベルの目盛りサイズ
plt.tight_layout() # ラベルがきれいに収まるよう表示
plt.grid() # グリッド表示
plt.show() # グラフの表示結果
100
----------------
s^2 + 12 s + 100
応答特性パラメータ比較
パラメータに対する挙動変化
減衰係数を変化させてステップ応答波形を比較します。
減衰係数 ζ を変化させると波形の振動が変化します。
1＜ ζ は非振動的波形、0 ≦ ζ ≦ 1 は振動的波形になります。
つまり ζ が小さくなるほど振動性が強くなります。
固有角振動数 ωn が大きくなると応答は早く、小さくなると応答は遅くなります。
減衰係数を変えた際の振動性の比較をしたいので ωn は固定値とします。
Pythonソースコード
行数が増えるので先ほどとは書き方を変えています。
固定値のパラメータはソースコードの初めに大文字で定義しています。
減衰係数を変化させるために for 文を使うことにします。
数式を ζ＝i×0.2 とし、i は 0 から 10 まで変化させ ζ は 0 から 2.0 まで 0.2 刻みで変化させます。
for文を回した後にグラフ表示設定を行い最後にグラフを表示させます。
ついでにボード線図も出力します。
#------------ ライブラリインポート ------------#
# 数値演算ライブラリ
import numpy as np
# 制御工学ライブラリ
from control import matlab
# グラフ描画
from matplotlib import pyplot as plt
#------------ 固定パラメータ ------------#
# ステップ応答時間設定
TIME_START = 0.0
TIME_END = 2.5
TIME_STEP = 1000
#------------ 関数定義 ------------#
# 2次系伝達関数モデル作成の関数
def tf_2nd_order(zeta, wn):
num = [wn*wn]
den = [1, 2*zeta*wn, wn*wn]
G = matlab.tf(num, den)
return G
#------------ Main ------------#
# 伝達関数モデル作成 & ステップ応答導出
for i in range(0, 11):
param_zeta = round(i*0.2,1)
param_wn = 10.0
G2nd = tf_2nd_order(param_zeta, param_wn)
# ステップ応答
t = np.linspace(TIME_START, TIME_END, TIME_STEP)
(output, t) = matlab.step(G2nd, t)
# 波形出力パラメータ
# (ls:線種, lw:線太さ, alpha:透過度, color:色指定, marker:マーカー形状指定, ms:マーカーサイズ, label:凡例名)
label_name = "ζ = " + str(param_zeta)
plt.plot(t, output, alpha=0.7, label=label_name)
# 収束値を点線表示
plt.axhline(1.0, ls=":", color="blue")
# グラフ表示設定
plt.rcParams['font.family'] = 'Times New Roman' # 全体のフォント
plt.legend(loc="upper right",fontsize=8) # 凡例の表示（2：位置は第二象限）
plt.title('2nd order step', fontsize=10) # グラフタイトル
plt.xlabel('time [sec]', fontsize=10) # x軸ラベル
plt.ylabel('Output', fontsize=10) # y軸ラベル
plt.xlim([0, TIME_END]) # x軸範囲
plt.ylim([0, 2.0]) # y軸範囲
plt.tick_params(labelsize = 10) # 軸ラベルの目盛りサイズ
plt.tight_layout() # ラベルがきれいに収まるよう表示
plt.grid() # グリッド表示
plt.show() # グラフの表示
# 伝達関数モデル作成 & ボード線図
for i in range(0, 11):
param_zeta = round(i*0.2,1)
param_wn = 10.0
G2nd = tf_2nd_order(param_zeta, param_wn)
# ボード線図作成
label_name = "ζ = " + str(param_zeta)
matlab.bode(G2nd, label=label_name)
# グラフ表示設定
plt.legend(loc="upper right",fontsize=8) # 凡例の表示（2：位置は第二象限）
plt.show() # グラフの表示
2次系伝達関数モデルステップ応答シミュレーション
1＜ ζ は非振動的波形、0 ≦ ζ ≦ 1 は振動的波形になっていることが分かります。
ζ が小さくなるにつれて振動性が強くなっています。
このことから2次系のシステム制御設計においてオーバーシュートさせたくない場合は 1 ＜ ζ で設計するとよいです。
また、ボード線図から -40dB で減衰して位相は -180° になることが分かります。ζ が大きいほど、位相は緩やかに低下していきます。
おわりに
本記事では2次系伝達関数のステップ応答をシミュレーションしました。
1＜ ζ は非振動的波形、0 ≦ ζ ≦ 1 は振動的な波形になることが分かりました。
2次系の制御設計において、ステップ入力信号に対して出力を振動させたくなければ 1＜ ζ で設計すれば良いということが分かります。
