皆さんこんにちは。AmedTech代表の天野です。
前回はMicrosoft Excelを使って伝達関数をシミュレーションする方法をご紹介させていただきました。今回はその続きとして二次振動性要素のシミュレーションをご紹介いたします。
2階常微分の差分化
前回微分を差分化することによって表計算でも微分を計算できることをお示しいたしました。これは2階常微分にも拡張が可能です。
微分を差分に変換する方法の復習ですが、前回ご紹介したように現在と未来の値を使う方法（前方差分法）、現在と過去の値を使う方法（後方差分法）、過去と未来の値を使い方法の三種類がありました。ここでは過去の値を使う差分法（後方差分法）を利用します。
$$\frac{dx}{dt}=\frac{x(t)-x(t-1)}{\Delta t}\tag{1}$$
これを2階微分に拡張すると、次に様になります。
$$\frac{d^2x}{dt^2}=\frac{d\left.(\frac{dx}{dt}\right.)}{dt}=\frac{\frac{x_{n+1}-x_{n}}{\Delta t}-\frac{x_{n}-x_{n-1}}{\Delta t}}{\Delta t}=\frac{x_{n+1}+x_{n-1}-2x_{n}}{\Delta t^2}\tag{2}$$
ただしここで気を付けなくてはならないのは未来の値x(n+1)が入ってしまっていることです。これは、エクセルで計算する際には2つ前の値を利用することで代用します。すなわち
$$\frac{d^2x}{dt^2}=\frac{d\left.(\frac{dx}{dt}\right.)}{dt}=\frac{\frac{x_{n}-x_{n-1}}{\Delta t}-\frac{x_{n-1}-x_{n-2}}{\Delta t}}{\Delta t}=\frac{x_{n}+x_{n-2}-2x_{n-1}}{\Delta t^2}\tag{2.1}$$
とすることで（正確ではありませんが）十分実用になります。実際にやってみましょう。
前回は1次遅れ要素をエクセルでシミュレーションしましたが、今回は2次振動性要素をシミュレーションしてみます。2次振動性要素の伝達関数は
$$G(s)=\frac{\omega_{n}^2}{s^2+2\zeta\omega_{n}s+\omega_{n}^2}\tag{3}$$
と表されますが、定常ゲインを１とする次の表現を用いことが多いです。
$$G(s)=\frac{1}{1+2\zeta(\frac{s}{\omega_{n}})+(\frac{s}{\omega_{n}})^2}\tag{4}$$
\(\omega_{n}\)を固定角周波数\(\zeta\)をダンピングレシオといいます。
では上式(3)をエクセルでシミュレーションしてみましょう。まずは左辺に出力x(t)右辺に入力y(t)にまとめます。（覚えてますか？入力はy(t)出力はx(t)でしたね。）
$$G(s)=\frac{X(s)}{Y(s)}=\frac{\omega_{n}^2}{s^2+2\zeta\omega_{n}s+\omega_{n}^2}$$
$$s^2X(s)+2\zeta\omega_{n}sX(s)+\omega_{n}^2X(s)=\omega_{n}^2Y(s)$$
ここに上式(2.1)を第1項に、式(1)を第2項に適用します。この結果が次の式ですが、周波数領域X(s)から時間領域x(t)に変わっていることに注意してください。
$$\frac{x(t)-2x(t-1)+x(t-2)}{dt^2}+2\zeta\omega_{n}\frac{x(t)-x(t-1)}{dt}+\omega_{n}^2x(t)=\omega_{n}^2y(t)$$
次にこれを左辺にx(t)、右辺にy(t)、x(t-1)、x(t-2)をまとめます。
$$\left.(\frac{1}{dt^2}+\frac{2\zeta\omega_{n}}{dt}+\omega_{n}^2\right.)x(t)=\omega_{n}^2y(t)+\left.(\frac{2}{dt^2}+\frac{2\zeta\omega_{n}}{dt}\right.)x(t-1)-\frac{1}{dt^2}x(t-2)$$
ここで左辺は
$$\frac{1+2\zeta\omega_{n}dt+\omega_{n}^2dt^2}{dt^2}x(t)$$
となるので、最終的には
$$x(t)=\frac{dt^2}{1+2\zeta\omega_{n}dt+\omega_{n}^2dt^2}\biggl(\omega_{n}^2y(t)+\Bigl(\frac{2+2\zeta\omega_{n}dt}{dt^2}\Bigl)x(t-1)-\frac{1}{dt^2}x(t-2)\biggl)\tag{5}$$
となります。y(t)には入力値を、x(t-1)、x(t-2)にはそれぞれ1つ前と2つ前の出力値を入れます。エクセルではこのようになります。
G列の各数字にはF列にある名前を付けてください。（メニュー→数式→名前の定義）
D列は式（5）を入れます。
D6=dt^2/(1+2*Zeta*Wn*dt+Wn^2*dt^2)*(Wn^2*B6-1/dt^2*D4+(2/dt^2+2*Zeta*Wn/dt)*D5)
C列には教科書に載っている以下の2次振動性要素のステップ応答の式を入れます。ステップ応答は\(\zeta\)の大きさによって式が変わるので、ここではif文を使って場合分けをしています。
$$0<\zeta<1のとき\hspace{10pt} h(t)=1-\frac{\mathrm{e}^{-\zeta\omega_{n}t}}{\sqrt{1-\zeta^{2}}}\sin(\sqrt{1-\zeta^{2}}\omega_{n}t+\tan^{-1}\frac{\sqrt{1-\zeta^2}}{\zeta})$$
$$\zeta=1のとき\hspace{10pt}h(t)=1-\mathrm{e}^{-\omega_{n}t}(1+\omega_{n}t)$$
$$\zeta>1のとき\hspace{10pt}h(t)=1-\frac{\zeta+\sqrt{\zeta^2-1}}{2\sqrt{\zeta^2-1}}\mathrm{e}^{-(\zeta-\sqrt{\zeta^2-1})\omega_{n}t}+\frac{\zeta-\sqrt{\zeta^2-1}}{2\sqrt{\zeta^2-1}}\mathrm{e}^{-(\zeta+\sqrt{\zeta^2-1})\omega_{n}t}$$
C6=IF(Zeta<=0,”Error”,IF(Zeta<1,1-(EXP(-Zeta*Wn*A6)/SQRT(1-Zeta^2))*SIN(SQRT(1-Zeta^2)*Wn*A6+ATAN(SQRT(1-Zeta^2)/Zeta)),IF(Zeta=1,1-EXP(-Wn*A6)*(1+Wn*A6),1-((Zeta+SQRT(Zeta^2-1))/(2*SQRT(Zeta^2-1))*EXP(-(Zeta-SQRT(Zeta^2-1))*Wn*A6))+((Zeta-SQRT(Zeta^2-1))/(2*SQRT(Zeta^2-1))*EXP(-(Zeta+SQRT(Zeta^2-1))*Wn*A6)))))
E列には差分の二乗を計算させています。この結果を次に示します。（ただしこの例ではdt=0.2、\(\omega_{n}\)=0.2でシミュレーションしています。）青が入力でステップ入力、赤が厳密解、緑が差分法による計算結果、そして紫が厳密解と差分解の差の二乗でスケールは右側です。まずは過剰補償(Over damping)でζ=1.2の時。
次に適正補償(Critical damping)でζ=1.0の時。
これら2つのグラフでは厳密解と差分法による解の違いは最大で0.025％程度ですが…
ダンピングファクターが0.8,0.5,0.2と小さくなるにつれて誤差が広がり、0.2では0.1％を超えます。とはいえ、この程度であれば今後出てくるPID制御のシミュレーションには十分使えます。特にMathematicaやMATLABなどの科学計算ソフトを使わずに、多くのPCにインストールされている表計算ソフトであるMicrosoft Excelを使って微分方程式の解を比較的簡便に得られるのは、結構使い勝手が良いのではないかと思います。