ロボットシミュレーション用のルンゲ・クッタ法プログラムの検討（１）ルンゲ・クッタ関数

はじめに

趣味のマイクロマウス活動において、ずーっと思い描いている妄想に、マイクロマウスシミュレータがあります。

ずいぶん昔に斜め探索の探索アルゴリズムの検討のために探索シュミレータを作ったことがありますが、本当の意味でのマイクロマウスの動作シミュレーションを作ってみたいです。 github.com

最近のお気に入り言語はPythonなのでこれを使って作ってみようかと思っています。そこで微分方程式を解く場面を想定して、まずはルンゲ・クッタ法を用いて微分方程式を解く汎用の関数の開発を始めてみました。

最近はMatlab＆Simulinkも気軽に使える素晴らしい状態になったのでシミュレーションを手軽にしたい場合はそれらのプログラムを使うことをお勧めします。

理系的な戯れとしては、その辺のところを楽しむのが趣旨なので、わざわざ苦労します。

でもね、勝負に勝ちたいと思った場合は手段は選んではいけないのだと思います。ずいぶん前にYariYari2000と言う恥ずかしい名前の４WSのロボットを作ったんですけどロボットのステアリングが前後２箇所とも切れたので制御的には多入力他出力系の制御問題を解く現代制御理論的なことをやりました。この時は、Matlabに頼らざる得なかったので使いました。PID制御までならMatlab等使わずになんとかなりそうですが、無いと辛いと思います。

ロボット競技の楽しみ方もそれぞれなので、その辺りも自分でやるのも楽しいかなとも思います。

ちなみにMatlabのシミュレーションのデフォルトはode45と言うアルゴリズムになっています。これは４次と５次のルンゲ・クッタ法で計算して、その差を見ながら、刻み幅を自動的に調整して計算精度を保証しようと言うものらしいです。２０年以上前からこれが用いられていますので、マイナーな更新があるとは思いますが、大変信頼性のある方法なのだと考えられます。

ルンゲ・クッタ法　Runge-Kutta Methods

ルンゲ・クッタ法についてはいろいろなところで説明がなされているので、内容についてはそれらに任せたいと思います。このメソッドはドイツの学者カール・ルンゲ先生とマルティン・クッタ先生が共同で開発したメソッドです。お二人ともゲッティンゲン大学 で研究をされていて、航空工学を少し勉強した者としては、ゲッティンゲン風洞を思い出すと同時に、クッタ先生はあのクッタ・ジュコーフスキーの定理のクッタ先生だと気が付き更に感慨深いものがあります。

いつの日か戯れとしてルンゲ・クッタ法について考えてみたいなと思います。

ルンゲ・クッタ法によって何ができるかと言うと「常微分方程式を数値的に解ける」と言うことです。

ルンゲ＝クッタ法 - Wikipedia

ここで取り上げるのは４次ルンゲ・クッタ法です。結論の数式だけ以下に示すと

求めたい時間の関数を ${x}$ として、その導関数が以下だとします。

${\displaystyle \frac{dx}{dt}=f(t,x) }$

すると、初期値を適切に設定して、以下の式を順に計算して求めていくと、最後に求めたい変数 ${x}$ の次の値がもとまります。

${\displaystyle k_1=h f(t_n, x_n)\\ k_2=h f(t_n+\frac{h}{2}, x_n+\frac{k_1}{2} )\\ k_3=h f(t_n+\frac{h}{2}, x_n+\frac{k_2}{2} )\\ k_4=h f(t_n+h, x_n+k_3 )\\ }$

${\displaystyle x_{n+1}=x_n+\frac{k_1 + 2 k_2 + 2 k_3 + k_4}{6} }$

以上で一回分の計算がおわり、求められた値を使って、同じことを繰り返すのが数値計算的に微分方程式を解く流れになります。

汎用的なルンゲ・クッタ関数を作るには

引数として関数を渡す

ルンゲ・クッタ法の部分を関数化しようとした際に、様々な運動方程式や回路方程式を表す導関数に対応しなければなりません。そのためには関数を引数として渡すと言うことを実現しなければなりませんがPythonでは、普通に関数名を引数として渡すだけです。

引数が可変でも対応できないとならない

導関数が使用する変数は様々な導関数によって種類や数が違うと考えられます。そこで、Pythonでの可変引数の実現方法ですが関数の引数リストの最後にアスタリスク＊をつけた変数を置きます。Cのポインタと似ています。

サンプルプログラム

以上のことを踏まえてバネマスダンパ系のステップ応答を計算するプログラムを作成してみました。一応コードにコメントをつけて使い方等は判るように書いたつもりです。計算が正しいかは次回の投稿で２次振動系のステップ応答について取り上げたいと思いますので、そちらで検証します。

github.com

計算結果

以下に計算結果のグラフを示します。雰囲気的にはうまく計算できているような気がします。

おわりに

趣味のロボット作りに生かせるロボットシミュレーションプロジェクトについて書いてみました。

どう書けば使いやすいコードになっていくのか使って行かないと判らないと思いますので、これから少しづつ書いていけたらと思います。一応のターゲットはマイクロマウスです。マイクロマウスのコードはほぼC言語またはC++で書かれているのでPythonでシミュレーションすることが本当にいいのかは疑問を持っていたところです。一番良いのはC言語で書かれたマウスのコードとPythonで書かれたシミュレータがリンクして動かせればと思います。やはりこれは、マウスハードの仮想環境とか作らないとダメなのかと考えると完成は遥か遠くに感じます。とりあえずはマウスのコードがそのまま試せると言うのは次の目標にして、アルゴリズムの検討や制御パラメータの調整に生かせるものにすると言うのが遠い将来の目標です。

次回は計算結果が正しいのかを検討するために、２次振動系のステップ応答の解析解を紹介しつつ検討したいと思います。さらに、その次はファルハーバーの1717モータのシミュレーションをしてみます。

次回以降の予定