octaveの微分代数方程式のソルバー　DASPK solverに関して

ご無沙汰しておりました。
我がチームの設計の方も一段落つき、設計としてもっとも忙しい時期は超えたなという感がある季節になりました。
およそ2週間後に迫る荷重試験が、これまた設計にとってもっとも過酷な試験ではありますが、それさえクリア出来れば二年目の設計としての試験飛行に突入するのを待つのみです。

さて、今回は最近取り組んでいる人力飛行機のフライトシミュレータの紹介と、そのために用いたoctave関数「daspk」の紹介をしておこう。

これが鳥人間コンテストにおける発進と旋回をシミュレートしたものである。操舵はoctaveの割り込み入力関数kbhitを用いてリアルタイム(笑)で入力する。
左上に出ているグラフの上が機速対時間　下が出力対時間である。緑の方がパイロット出力だ。右下のグラフは状態のグラフ。上が縦で、下が横・方向だ。詳しい中身は少し機体運動学を勉強すれば読めるので割愛する。

…いかに人力飛行機の操縦が難しいか分かる。UIが整っていないおかげで、三次元位置のグラフと、文字で表示される機体姿勢から操舵を入力しているのだが、旋回を終えたあたりで激しく高度が上下している。そして高度的には着水している。（これでもまだ上手くなった)　
風の影響をまだ入れていないため、横・方向について旋回までは全く直線を維持している。風無しでこれか…というのが正直な感想。

人力飛行機の操縦の難しさは有次元安定微係数から考えれば理論立てて述べることは出来るのだが、この話は長くなるので割愛。

これがUI。ざっと調べたけれどもoctave単体ではジョイスティック（可変抵抗)の値をリアルタイムで読んで操舵入力に反映するのは難しそうだ。出来たとしても、octaveの貧弱なグラフィックスではたぶん処理落ちしてしまう。

そこで、他の汎用言語に移りたいと思うわけで、現在その勉強を始めている。

しかし、それを大きく阻む壁がある。これはまさにoctaveの強い魅力という意味でもあるのだが。

「微分代数方程式のソルバー」である。

連立常微分方程式は4次のRunge-Kutta法など多くのソルバーもあるし、理論もそこまで難しくない。そもそも解くだけならそこら中にソルバーのアルゴリズムは公開されている。

しかし、こうした常微分方程式のソルバーが解けるのは、
$\frac{dx}{dy}=...$
という形のみで構成された連立常微分方程式のみである。

私はそこまで数学に詳しくないので、これだけで微分代数方程式も解けると言われてしまえばこの記事のお役はここで終わりである。
しかし、言い訳がましくなってしまうが学部二年で人力飛行機の設計をしなければならない場合が多い人力飛行機界において、「手軽に解けること」は非常に重要だと考える。

そこで登場するのが「daspk」なのである。

そもそも、微分代数方程式とは何よ。という方もいると思うので、軽く説明。わたしもちゃんと分かっていないので何とも言えないが、上に述べたような形に直せない(直すのがとっても面倒くさい)微分方程式だとかってに思っている。たとえば、
$(\frac{d}{dt}-X_{u})u-X_\alpha\alpha+(W_0\frac{d}{dt}+gcos\theta_0)\theta=X_{\delta t}\delta_{t}$
$-Z_u+(U_0\frac{d}{dt}-Z_{\alpha})\alpha-(U_0+Z_q)\frac{d}{dt}-gsin\theta_0)\theta=Z_{\delta e}\delta_{e}+Z_{\delta t}\delta_{t}$
$-M_u+(M_{\dot{\alpha}}\frac{d}{dt}+M_\alpha)\alpha+({\frac{d}{dt}}^2-M_q\frac{d}{dt})\theta=M_{\delta e}\delta_{e}+M_{\delta t}\delta_{t}$
$\frac{d}{dt}\theta=q$
(縦の微少擾乱方程式,加藤寛一郎,大屋昭男,柄沢研二,航空機力学入門より)
といった具合の方程式である。

daspkはその汎用性が高いのに、octaveの常微分方程式のソルバー「lsode」より存在感が薄いように思われる。
daspkの使い方は以下の通り

まず解く関数を定義する。
function res=yuukivel_sample(x,dx,t, ....(使う定数たち))
res(1)=dx-(x+tan(x)+spline(〜)...
res(2)=〜
〜
endfunction

といった感じだ。微分代数方程式は0=(右辺)として記述する。daspkは定義された関数においてres=0となるように値を推測して計算を進めていくらしい。もちろん数値計算であるから、非線形でも構わない。splineを関数内に入れることだって出来る。（もちろん計算は遅くなる)

呼び出しは
[x dx]=daspk(@(x,dx,t)yuukivel_sample(x,dt,t, ...(使う定数たち)),x0,dx0,t)
のたった一行。@以下は無名関数といって、「yuukivel_sample関数の変数は@内ですよ」と述べている。関数名の後は関数のパラメータであり、そのあとはxの初期値(x0)、xの微分値の初期値(dx0)、値を求めたい時刻tの値となる。

ね。簡単でしょ。慣れてしまえばlsodeの登場する幕はなくなる気がする。ただし、どちらが精度が良いのかについては正直、分からない。

私は、daspkを用いて微少擾乱の近似を弱めたものと、オイラー角と位置を求める微分方程式を関数内に突っ込んで解かせている。式の行数は15行にもなる。（まだ完成していないので無駄な行もあったりする）その結果が冒頭に上げたシミュレーション結果である。

機体運動解析は、機体の挙動をみるという純粋に解析的な用途だけでなく、パイロットの操縦訓練、設計のパラメータ決めの決定打になりうる。私個人としては、運動解析が手軽になることによって尾翼容積や動ファクタ比の値で尾翼周りの設計をする時代が終焉を迎えてくれればなぁ、と思ったりもする。もちろん、安定微係数の推算値が現実と大きく異なることもありうる。桁の剛性なども機体挙動にもかかわってくる。そのあたりをいかに意識して運動解析を使うかもまた、設計者の腕の見せ所なのかなぁとつらつらと考える。