次の記事(成功しました！)
線形回帰分析を用いたロケットの頂点時間の導出②[C言語,Arduino]
1. 概要
この記事は，時間と気圧高度のデータから水ロケットの軌道を推定し，頂点時間を計算する手法を解説したものです．結論から言うと，今回作成したプログラムでは水ロケットのような短時間で頂点高度に到達するものでは上手く行きませんでした．
しかしハイブリッドロケットなど長時間飛行するようなロケットなどではうまくいくことが期待されるので，とりあえずメモとして残します．
もうしばらくやる予定なのでまた今度もう少し良いモデルを考えて実験します．
GitHubのリポジトリ作っときましたので参考までに．コードはここに投げていきます．
3.式の導出
まずロケットの挙動がどうなっているか，ロケットのシミュレーションを行います．今回はOpenRocketを用いて行います．ただ，エンジンにペットボトルは選択できないので適当に作ったロケットで見てみます．
自由落下運動ならば，高度 \(y\) をプロットしてみると放物線運動となり， \[ y = a t^2 + b t + c \] という式となることが予想されます．実際に，上のグラフを見てもエンジン燃焼終了後は概ね放物線運動をしていることが読み取れるかと思います．
放物線運動(二次関数)のパラメータは3つなので，3点のデータがあれば連立方程式を解いて \(a,b,c\) は求まり，頂点に到達する時刻 \(t=-\frac{b}{2a}\) を計算することができます．
しかし，センサーから読み取られたデータはノイズがある程度存在するため，正しい結果が得られない可能性があります．そこで今回は線形回帰分析の知識を使って， \(a,b,c\) を推定してみます．
センサーから得られた高度と時間のデータ \({(t_1,y_1),...,(t_n,y_n)}\) があります．ここでこれらのデータが放物線 \[ y = a t^2 + b t + c \] にノイズ \(e_i\) が乗ったものとして考えていきます．
ノイズ \(e_i\) は元の関数との差なので \[e_i=a t_i^2 + b t_i +c - y_i\] と表す事ができます．ここでノイズの二乗和を評価関数 \(J\) として定義します．(普通の最小二乗法をします) \[ J = \frac{1}{2} \sum_{i=1}^{n} (a t_i^2 + b t_i +c - y_i)^2 \]
この評価関数を最小とするようなパラメータ \(a,b,c\) が，データに一番フィットするような放物線のパラメータと成るはずです．評価関数 \(J\) を\(a,b,c\) について偏微分したものが0となる連立方程式を考えます．
\[ \begin{eqnarray*} \frac{\partial J}{\partial a} &=& \sum_{i=1}^n(a t_i^2 + bt_i + c - y_i )t_i^2 = 0 \\ \frac{\partial J}{\partial b} &=& \sum_{i=1}^n(a t_i^2 + bt_i + c - y_i )t_i = 0 \\ \frac{\partial J}{\partial c} &=& \sum_{i=1}^n(a t_i^2 + bt_i + c - y_i )= 0 \\ \end{eqnarray*} \]
行列形式にすれば， \[ ( \sum_{i=1}^n \left[ \begin{array}{rrr} t_i^4 & t_i^3 & t_i^2 \\ t_i^3 & t_i^2 & t_i \\ t_i^2 & t_i & 1 \end{array} \right] ) \left[ \begin{array}{rrr} a \\ b \\ c \end{array} \right] = ( \sum_{i=1}^n\left[ \begin{array}{rrr} t_i^2 y_i \\ t_i y_i \\ y_i \end{array} \right] ) \] となりこの連立方程式をとけば評価関数 \(J\) について最適なパラメータ \(a,b,c\) がわかります． (※データ数が少なかったり運が悪いと正則行列でない場合があるのでそうならないようにお祈りしましょう(?))
また，地球の地表付近で観測された放物線運動ならば \(a\) は 重力加速度 \(g\) を用いればおよそ \({-\frac{1}{2}g}\) 付近であるはずです． 計算した結果のパラメータが正しいかどうかを判定するのに使えるかと思います
3.シミュレーション
どの程度データ数があれば十分かや，実装したプログラムが正しいかどうかを確認する意味も含めてシミュレーションを行いました．シミュレーションでは高度のデータにノイズとして±0.5mの乱数を加えています．また，あとで示すArduinoのプログラムではループ時間が約0.03秒であったので，それに従ってシミュレーションを行いました．
プログラムはこんな感じ．REPマクロとか修正してなかったけど大丈夫でしょう．余分にライブラリincludeしてるけど気にしないでください．
結果はこんな感じです．Markdownはいいぞ……．
上のグラフを見てない人も見た人も，なんかヤヴァい値が見えますね．
データの個数が100個以下ではうまくいくときもあれば，上手く行かないときもあるというような非常に怖い結果が出ています．500個のデータが有ればどれも上手く言っているようですが，少なくとも水ロケットでは使い物にならなそうです．
なぜなのかはデータをプロットしてやれば明確で，実際に三番目のテストケースのデータの個数が100の場合のやつをプロットしてみるとこうなります．
時間が短すぎてもはや直線で，放物線には見えないですね，これでは値が意味不明になるのもわかります．こればっかりはどうしようもないので，次はもう少し制約の強いモデルを考えて見ますので次回の記事を読んでもらうことを期待します……．
4.Arduinoのコード
一応ペットボトルロケットでやってみましたが上手く行かなかったです．だいぶ長くなるのでGitHubのリンクを貼っときます．
センサーはmpu9250とbmp180です．加速度のノルムが一定数以上になった場合，サンプリング周期0.03sで高度データを取得，30個データがたまった時点で計算を行います．
もう一度いいますがこれは上手く行かなかったやつなので，また別のプログラムを作る予定ですが，数学モデルが変わるのでこの負の遺産を残しておきます．
5.結論とか
上手く行かなかったのは多分制約が弱すぎるからなので，とりあえずa=-4.9の制約を入れてどうなるかやって，それでだめなら \[ \frac{dy}{dt} = -gt + b \] を基底関数としてやったりしてみます．
追記：
次の記事(成功しました！)
線形回帰分析を用いたロケットの頂点時間の導出②[C言語,Arduino]
追記：
次の記事(成功しました！)
線形回帰分析を用いたロケットの頂点時間の導出②[C言語,Arduino]
0 件のコメント:
コメントを投稿