機械学習による弾道シミュレーションのパラメータ推定

弾道シミュレーションをPythonで実装し、シミュレーションに用いるパラメータを観測値に合わせて調整する方法を解説します。パラメータ調整には機械学習を使用しますが、ここでは微分を直接計算せず、パラメータに摂動を加えて誤差変化を観察する原始的なアプローチをとります。

題材として、ゴルフボールの弾道計算を取り扱います。この記事で解説している内容は、GolfTaught Memoの弾道・飛距離計算機能を実装するために使用されたものです。コードの全体像は、GolfTaught Memoのリポジトリ内のJupyter notebookを参照してください。

回転するボールの弾道シミュレーション

ゴルフボールの運動方程式を設定します。ボールの位置・速度を変数として、その他の決定困難な量はパラメータとして抽出します。空気（＝流体）中の物体は様々な要素の影響を受けて運動しますが、ここではそれらを定数と仮定して可能な限り統合・簡略化することにします。

ボールの運動方程式

飛翔中のゴルフボールの運動方程式を以下のように表現します。

$$\frac{dv_y}{dt}=C_L v_x \ – g \ – k_1 v_y \ – k_2 v_y |v_y|\tag{1}$$

$$\frac{dv_x}{dt}=-C_L v_y \ – k_1 v_x \ – k_2 v_x |v_x|\tag{2}$$

ここで、 \(g=9.8m/s^2\) です。

空気の抵抗力

ボールに働く空気抵抗には、粘性抵抗と慣性抵抗が存在します。それぞれ速度の1乗・2乗に比例するため、係数 \(k_1, k_2\) を用いて表現しています。

マグヌス効果による揚力

回転する球と空気の相互関係から、マグヌス効果による揚力が生まれます。揚力はボールの速度に比例するとし、その係数を \(C_L\) とします。ここでは便宜的に、（正確な表現ではありませんが）この係数をマグヌス係数と呼ぶことにします。

\(C_L\) にはボールの回転数 \(n\) が影響します。ボールの回転数は飛翔中に減衰していくことが想定されるため、これを指数関数を用いて表現します。

$$n(t)=C\exp(-k_t\cdot t)$$

ここで \(C\) は定数であり、ボールの回転数の初期値に相当します。 \(C_L\) には他にも空気の流速や密度が影響しますが、これらはすべて定数であると仮定し、 \(C\) も含めた1つの定数にまとめて

$$C_L(t)=k_m n(t)=k_m\exp(-k_t\cdot t)\tag{3}$$

と表現します。

弾道シミュレーションの実装（Python）

以上の運動方程式にもとづいて、ゴルフボールの弾道シミュレーションを実装します。

def simulate(bs, la, params=[0.1, 0.001, 0.1, 0.1], dt=0.01):
    k1, k2, km, kt = params
    v0 = mph2ms(bs)
    vy = v0 * np.sin(deg2rad(la))
    vx = v0 * np.cos(deg2rad(la))
    xy = np.zeros_like(vy) + 1e-7
    xx = np.zeros_like(vx)
    traj = np.zeros((len(v0), 10000, 2))
    for t in range(10000):
        update = xy > 0
        if update.sum() == 0:
            traj = traj[:, :t, :]
            break
        vy_u, vx_u = vy[update].copy(), vx[update].copy()
        Magnus_coeff = km * np.exp(-kt * t * dt)
        vy[update] += (Magnus_coeff * vx_u - 9.8 - k1 * vy_u - k2 * vy_u * np.abs(vy_u)) * dt
        vx[update] -= (Magnus_coeff * vy_u + k1 * vx_u + k2 * vx_u * np.abs(vx_u)) * dt
        xy[update] += vy_u * dt
        xx[update] += vx_u * dt
        traj[:, t, :] = np.array([xx, xy]).T
    return traj

simulate関数は弾道をシミュレートし、ボールの軌跡を返す関数です。

弾道シミュレーションのパラメータ推定

前章の結果から、このシミュレーションは \(k_1, k_2, k_m, k_t\) の4つのパラメータを持つことがわかりました。ここでは、原始的な機械学習法を用いてこれらのパラメータを推定します。

ゴルフボールの弾道データの取得

ゴルフ競技では最大14本のクラブを使用することができます。これらのクラブは、ボールの打ち出し角を決めるロフト角の大きさが異なるため、それらを同じ力で振っても、ボールの飛距離が異なります。具体的には、ロフト角が大きいほど球が高く上がり、その分飛距離が落ちます。逆に、ドライバーと呼ばれるロフトの小さいクラブでは、低い弾道でボールが遠くまで飛びます。

①どのクラブで、②どのくらいの力（＝クラブのヘッドスピード）でボールを打つと、(A)どのくらい高く、(B)どのくらい遠くまでボールが飛ぶのかというデータは、たとえばTranckman Universityから取得することができます（データの再配布にあたるため、Githubリポジトリでは、取得したデータgolf_flight.xlsxの共有はしていません）。

弾道データを線形回帰で補間する

Trackman Universityからは、ドライバー（ロフト角10度）・6番アイアン（ロフト角30度）・ピッチングウェッジ（ロフト角45度）で打ったボールについての

1. 打ち出し初速（Ball Speed）
2. 打ち出し角（Launch Angle）
3. 打ち出し点から着弾点までの飛距離（Carry）
4. ボールの最高到達点（Height）
5. 着弾後のボールの転がりも含めた総飛距離（Total）

を取得することができます。これらとクラブのヘッドスピード（Club Speed）の関係をプロットすると、（変曲点を含む）線形モデル

$$\mathrm{Parameter_i}=a_i\times\mathrm{Club Speed}+b_i$$

に当てはめられることがわかります。

なお、Runとは、ボールが着弾した後に転がった距離で、Total – Carryで計算できます。

線形モデルの係数 \(a_i, b_i\) はクラブのロフト角によって変化しますが、これもさらに線形回帰を用いて、以下のように表すことができます。

$$a_i=a_\mathrm{slope}\times\mathrm{Loft Angle}+a_\mathrm{intercept}$$

$$b_i=b_\mathrm{slope}\times\mathrm{Loft Angle}+b_\mathrm{intercept}$$

同じ人が異なるクラブを振った時のヘッドスピードは、その人の6番アイアンのクラブのヘッドスピード \(\pm\alpha\) で計算できることを仮定すると、取得した3種類のクラブ以外でボールを打った時の情報を補間することができます。

飛距離と高度からシミュレーションパラメータを推定する

以上より、各クラブを色々なヘッドスピードで振った時のボールの情報が得られます。このとき、それぞれの初速（Ball Speed）と角度（Launch Angle）でボールを打ち出すと、それぞれの飛距離（Carry）と高度（Height）の弾道になるはずです。シミュレーションがこれらのデータを再現できるように、4種類のパラメータ \(k_1, k_2, k_m, k_t\) を調整します。

誤差関数の設定

実際の飛距離・高度と、シミュレーションによる飛距離・高度の差は、それぞれの平均2乗誤差の重みづけ和で評価します。

def loss(carry, height, traj):
    return np.mean((traj[:, -1, 0] - carry)**2 + ((traj[:, :, 1].max(axis=1) - height) * 10)**2)

ここで、飛距離が最大300m程度であるのに対し、高度は30m程度なので、誤差のバラツキは飛距離のほうが \(10^2\) 倍程度大きくなります。これは誤差評価において飛距離が重視され、高度のズレが無視されやすいことを意味します。そこで、ここでは2乗計算の前に、高度の