ロボットシミュレーション用のルンゲ・クッタ法プログラムの検討（２）２次振動系のステップ応答　その２

はじめに

こんにちは、こうへいです。

前回、２次遅れ系のステップ応答の計算式を紹介しました。今回は２次遅れ系の代表例であるバネマスダンパ系のステップ応答を計算して、ルンゲ・クッタ法の結果と比較したいと思います。

kouhei-no-homepage.hatenablog.com

バネマスダンパ系の伝達関数と２次遅れ系の伝達関数

質量 ${M}$ 、ダンピング係数 ${\zeta}$ 、バネ定数 ${D}$ とすると、バネマスダンパ系の伝達関数は次のような式になります。

${ \displaystyle G(s)=\frac{1}{M s^2 + D s + K_{s}} }$

２次遅れ系の伝達関数と比較するために、２次の項の係数 ${M}$ で分母と分子を割ると

${ \displaystyle G(s)=\frac{\frac{1}{M}}{s^2 + \frac{D}{M} s + \frac{K_{s}}{M}} }$

となります。これと以下の２次遅れ系の伝達関数と比較すると

${ \displaystyle G(s)=\frac{K \omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2} }$

ゲイン ${K}$ 、減衰係数 ${\zeta}$ 、自然角振動数 ${\omega_n}$ が以下のように求まります。

${ \displaystyle K=\frac{1}{K_s}\\ \zeta=\frac{D}{2\sqrt{K_s M}}\\ \omega_n=\sqrt{\frac{K_s}{M}} }$

以上の結果と、以下の前回の投稿で示した２次遅れ系のステップ応答を表す式

${\zeta > 1}$ ：異なる二つの実数解になる場合のステップ応答

${ \displaystyle y(t)=K \left[ 1- \frac{1}{2 \sqrt{\zeta^2 - 1}} e^{- \zeta \omega_n t} \left\{ \left( \zeta + \sqrt{\zeta^2-1} \right) e^{\omega_n \sqrt{\zeta^2-1} t } \\ - \left( \zeta - \sqrt{\zeta^2-1}\right) e^{-\omega_n \sqrt{\zeta^2-1} t } \right\} \right] }$

${\zeta ＜ 1}$ ：共役複素解になる場合のステップ応答

${ \displaystyle y(t)=K \left\{ 1- \frac{1}{\sqrt{1-\zeta^2}} e^{- \zeta \omega_n t} \sin \left( \omega_n \sqrt{1-\zeta^2} t + \phi \right) \right\}\\ \phi=\tan^{-1} \frac{\sqrt{1-\zeta^2}}{\zeta} }$

${\zeta = 1}$ ：重解になる場合のステップ応答

${ \displaystyle y(t)=K \left\{ 1 - e^{-\omega_n t}( \omega_n t +1) \right\} }$

以上から、ゲインと減衰係数と自然角振動数と計算時間を引数として２次遅れ系のステップ応答を求める関数を作成しました。

サンプルではバネマスダンパ系の質量、ダンピング係数、バネ定数からゲインと減衰係数と自然角振動数を算出してステップ応答関数に与えています。

これでルンゲ・クッタ法のバネマスダンパ系の単位ステップ応答の計算結果と比較する準備ができました。

比較計算プログラム

github.com

計算結果

以上のプログラムを用いて、計算をしました。バネマスダンパ系のパラメータですが

D=3(Ns/m)、Ks=10(N/m)、M=1(kg)

です。

この場合

${K=0.1}$ (m)、 ${\zeta=0.475}$ 、 ${\omega_n=3.162}$ (rad/s)となります。自然角振動数の単位をHzに変換すると大体0.5Hzです。

まず、刻み幅が0.1秒の場合です。

f:id:kouhei_ito:20200104232929p:plain — 刻み幅0.1秒の場合

青い線がルンゲ・クッタ法による数値計算の結果、オレンジ色の線が解析解です。刻み幅が0.1秒では荒すぎて、誤差が大きようです。

次に、刻み幅を0.01秒にしました。

f:id:kouhei_ito:20200104234228p:plain — 刻み幅0.01秒の場合

0.01秒でかなり誤差も小さくなりましたが、動きが急なところでは誤差が目立ちます。

次に、刻み幅を0.001秒にしました。

f:id:kouhei_ito:20200104234644p:plain — 刻み幅0.001秒の場合

ここまで刻み幅を小さくすると、このスケールではほぼ誤差が見えません。

今回のサンプルでは刻み幅を0.001秒にすれば満足のいく結果となりましたが、この値はケースにより異なりますのでシステムの応答の速さ等が判る場合は参考にして決めるといいと考えます。

おわりに

ルンゲ・クッタ法と解析解を比較することができました。今回作成したルンゲ・クッタ関数は正常に動作しているようです。

今回のコードではメインの中でループを作って関数を呼び出しています。場当たり的な再利用を考えないコードならば結果さえ出れば良いので問題だとは思いませんが、今後再利用を視野に入れると、これが果たして使いやすいのか疑問があります。今回のサンプルのように単純に一つの物体のシミュレーションならば目立ちませんが複数の物体の計算をさせる場合等を考慮してコードを見直して行こうと思います。

加えて２次遅れ系のステップ応答を計算する関数もできたので自前Toolboxが揃って行くかもしれません。

次回は１７１７モータをシミュレーションします。