PR

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

AI便利帳
Sponsored

弾道シミュレーションを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乗計算の前に、高度の誤差に10倍の重みづけを行っています。

このloss関数を誤差関数として、機械学習によって、シミュレーションのパラメータを推定します。

（原始的な）機械学習の原理

4種類のパラメータ \(k_1, k_2, k_m, k_t\) は、弾道に対し複雑な影響を与えます。上図は6番アイアンを平均的なヘッドスピードで振る設定で、4種類のパラメータをそれぞれ10％増減させた場合の弾道の変化をプロットしたものです。

それぞれのパラメータを適切な方向に変化させると、誤差を減らすことができます。通常の機械学習では、パラメータを更新する方向を判定するために微分を用いますが、最終結果がシミュレーションに依存する場合は、微分を計算することができません。そこで、パラメータを少しだけ変えた場合の誤差の変動によって微分を近似します。

$$\frac{dL}{dk_i}\simeq\frac{L(k_i+\Delta k_i)-L(k_i)}{\Delta k_i}$$

ここで、 \(L(k_i)\) はパラメータ \(k_i\) でシミュレーションをした場合の誤差（それ以外のパラメータは固定）をあらわします。この考え方にもとづき、それぞれのパラメータを以下のように更新します。

$$k_i\leftarrow k_i-\alpha\times k_i\times\mathrm{max}\left(-1, \mathrm{min}\left(1, \frac{L(k_i+\Delta k_i)-L(k_i)}{c}\right)\right)$$

ここで、分母に \(\Delta k_i\) ではなく定数 \(c\) を用いるのは、微小数 \(\Delta k_i\) で割ることによる急な桁数増加を防ぐためです。また、 \(\alpha\) は学習率に相当します。ここではmax, minで更新量に範囲を設けているため、「最大でパラメータの \(\alpha\) 倍を増減させる」ことを意味します。

np.random.seed(42)
lr = 0.001
k1, k2, km, kt = 0.1, 0.001, 0.3, 0.3
ls = 0
ls_save = []
for i in range(3000):
    if i == 2000:
        lr = 0.0001
    bs, la, ca, he, _ = flight_stats(np.random.randint(9, 60, size=64), np.random.randint(60, 100, size=64))
    L = loss(ca, he, simulate(bs, la, [k1, k2, km, kt]))
    
    L_k1 = loss(ca, he, simulate(bs, la, [k1 * 1.001, k2, km, kt]))
    L_k2 = loss(ca, he, simulate(bs, la, [k1, k2 * 1.001, km, kt]))
    L_km = loss(ca, he, simulate(bs, la, [k1, k2, km * 1.001, kt]))
    L_kt = loss(ca, he, simulate(bs, la, [k1, k2, km, kt * 1.001]))
    
    k1 -= max(-1, min(1, (L_k1 - L) / 10)) * k1 * lr
    k2 -= max(-1, min(1, (L_k2 - L) / 10)) * k2 * lr
    km -= max(-1, min(1, (L_km - L) / 10)) * km * lr
    kt -= max(-1, min(1, (L_kt - L) / 10)) * kt * lr
    
    ls += L
    if (i + 1) % 100 == 0:
        ls_mean = ls / 100
        print("iter:{:>5} | loss={:>9.3f} | k1={:.3e}, k2={:.3e}, km={:.3e}, kt={:.3e}".format(i+1, ls_mean, k1, k2, km, kt))
        ls_save.append(ls_mean)
        ls = 0

以上による原始的な機械学習は3000イテレーション程度で収束し、

  • \(k_1=0.1376\)
  • \(k_2=0.0011\)
  • \(k_m=0.2284\)
  • \(k_t=0.3643\)

程度になることがわかりました。この結果を用いてシミュレーションを行った結果が下図です。

まとめ

この記事では、飛翔体の運動方程式を設定し、Pythonでシミュレーションを行いました。今回はクラブの種類によるボールの初期スピン量の違いを考慮しませんでしたが、これを式に含めることで精度が向上する可能性があります。

また、現実の観測値に合わせてシミュレーションのパラメータを決定する方法についても解説しました。誤差の計算に複雑な過程を必要とする場合、パラメータによる微分を直接計算することはできませんが、機械学習の原始形態ともいえる方法を利用することが可能です。

もっと知りたいこと、感想を教えてください！

Sponsored