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

はじめに

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

今回は、前回に続き開発中のルンゲ・クッタ関数を使ったバネマスダンパ系のステップ応答が正しいのかどうか確認してみたいと思います。（と、思いましたが長くなってしまったので、その１とその２に分けました。）

（すみません、間違っている事を書いている可能性も多分にあるので、もし万が一ここに書かれている事を試してご迷惑をおかけしたとしても、自己責任でお願いします。）

kouhei-no-homepage.hatenablog.com

制御工学で必ず出てくる２次遅れ系

制御工学では２次振動系もしくは２次遅れ系と言う呼称で表されるシステムが登場します。例としてバネマスダンパ系やLCR回路等です。これは数値計算で解かなくても解析的に微分方程式が解けるので、数値計算が正しいかどうかをベンチマークするのにちょうど良いと思います。

便利な２次遅れ系

また、２次振動系は振動数と減衰係数の二つのパラメータ（これにゲインを加えて３つのパラメータ）でその特性を表すことができます。複雑な現象を単純化したい場合に、その次元を落としたモデルとして使うこともありますので、１次遅れ系と併せて馴染んでおくと良いものです。

手計算できる２次遅れ系

加えて、２次遅れ系については手計算で解ける限界になります。２次遅れ系の「２次」と言う言葉は２次方程式の「２次」と同じです。２次方程式は有名な解の公式や平方完成で頑張れば絶対に解けます。２次遅れ系の解析には２次方程式を解くと言う動作が入ってくるので、公式等を使えば必ず解けると言うのが大事なことなのです。

　２次遅れ系

さて、２次遅れ系について考えていきます。ブログにこの話を書こうと思ったときに、どのくらい詳しく書くか迷いました。書こうと思えばいっぱい書けるけど、ルンゲ・クッタ関数のベンチマーク的なものの解説ならば、そんなにだらだらとは書けないので、簡略化バージョンにしたいと思います。

ちなみに２次振動系の解析をするのに必要な知識として

２次方程式の解の求め方
複素数
オイラーの公式
ラプラス変換と逆変換
部分分数展開
因数定理、指数法則、etc......

僕の授業ではこれらのことを教えていくのに半年かけてます。制御の話って殆どできない・・・・でも、これらの事を疎かにすると、自分では勉強できなくなると思っています。更に「でも」ですが、これらのことって教室で聞くときっとつまらないんですよね、学校の授業ってこの学問勉強すると「おもしろそうだな！」と思わせることも大事だと思ってるので、自分で勉強できる力を付けようと思って頑張ってるんですけど、その結果モチベーションを削いでいるような気もしてとてもジレンマです。

２次遅れ系の伝達関数

ごめんなさい、いきなり２次遅れ系の伝達関数を書きます。

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

これが２次遅れ系の伝達関数と呼ばれるものです。教科書によっては $K$ が無い場合もあります。この式の中の各パラメータは

${K}$ ：ゲイン

${\zeta}$ ：減衰係数

${\omega_n}$ ：自然角周波数（振動数）（rad/s）

となります。

伝達関数について

「伝達関数とはなんぞや？」と言う説明をちゃんとするのは申し訳ないですけど省力して、なんとなくの話で済ませたいと思います。

僕たちは「時間が経過すると物事がどうなっていくか？」と言う世界を通常は意識的にも無意識にも理解しています。その時間中心の世界が直感的に理解している世界なのですが、その世界で計算しようとすると結構難しい積分計算をガリガリしていかないとならなくなります。時間の世界でガリガリ計算するのは難しい積分計算をしないとならないのでしんどいけど、「一手間」かけると計算が四則演算レベルで済んでしまう世界があります。その世界を「周波数の領域」と呼んでいます。

周波数の世界に入るための「一手間」をラプラス変換といいます。ラプラス変換という魔法により計算を簡単にしていけます。でもその世界での計算もはじめての人には簡単には思えないかもしれませんね。でもこの方法はやはり画期的なことで、制御工学を構築した人たちはやはり大変な苦労をしてここに行き着いたんだろうなと思います。

伝達関数というのはそのラプラス変換されて周波数の世界での入力と出力の比（出力／入力）を表す式です。

したがって、伝達関数に入力の関数をかけると出力を計算できます。ですので、システムの性質を表す関数で、この関数を解析するとシステムの性質を解き明かすことができ伝達関数に基づくことにより、制御システムを設計することができます。

単位ステップ応答

システムに単位ステップ入力を与えた時の出力を単位ステップ応答と呼びます。伝達関数でシステムが記述されているときは伝達関数に ${1/s}$ をかけると求められます。出力を ${Y(s)}$ とすると、単位ステップ応答を表す式は以下のようになります。

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

逆ラプラス変換

上の単位ステップ入力を表す式を逆ラプラス変換すると式が時間の世界の式に戻って時間が経ったら出力がどうなるかが判るようになります。

逆ラプラス変換するには分母の２次式を解の公式と因数定理を使って因数分解します。そのあと、部分分数展開というテクニックで１次の有利関数（分数形式で表された式）の和の形に上の式を変形します。すると、１次の有理関数を逆ラプラス変換するにはラプラス変換表と言うカンニングペーパーがあってそれを見ると一瞬にできるので、逆ラプラス変換が完了します。

２次式の因数分解

２次式の因数分解をするときに２次式＝０の２次方程式の解が ${\alpha}$ 、 ${\beta}$ だとすると、２次式は因数定理によると以下の式のように因数分解できます。（２次の項の係数を ${a}$ とします。）

${ \displaystyle a(s-\alpha)(s-\beta)=0 }$

よって、２次方程式の解がわかれば左辺の２次式は必ず因数分解できます。

ここで以下の伝達関数の分母の２次式＝０の式

${ \displaystyle s^2 + 2 \zeta \omega_n s + \omega_n^2= 0 }$

を解いた解が、異なる二つの実数解、共役複素解、重解になるかで因数分解の結果が変わります。これは高校で学習する判別式を立てると、その条件が判るのですが、判別式を立てて調べると以下のことがわかります。

${\zeta > 1}$ ：異なる二つの実数解
${\zeta ＜ 1}$ ：共役複素解
${\zeta = 1}$ ：重解

逆ラプラス変換の過程では、以上の因数分解の結果の場合わけにより計算過程が変わってくるのでそれぞれ計算を行って、それぞれの場合わけごとに異なった答えを得ます。

今回の記事ではガリガリ逆ラプラス変換する過程は省力します。授業でやると学生はゲンナリしています。とってもしんどそうです。

２次振動系のステップ応答をあらわす時間関数

それでは、以下に各場合におけるステップ応答の式を掲載します。

${\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] }$

振動しないでKに収束しそうですね。減衰係数と角振動数の大きさが収束の速さに影響するようです。

${\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} }$