物理シミュレーションと数値積分

最近、グランツーリスモ7のプレイヤーから、物理シミュレーションが破綻する現象が報告されています。特にサンババスで極端な破綻が起きるとのことで、空中に60万km/h以上で打ち出される事例の報告までされています。報告を見る限り、リアにバラストを積み、車高を下げることで発生するようです。

一般に物理エンジンは数値積分を内部的に使っています。物理エンジンでシミュレーションする系に固いレートのバネが入っていると計算が発散しやすく、 GT7で報告されているような車が宙に突如放り出されるような挙動に陥りやすくなります。 GT7でも同じことが起きているとまでは断言できませんが、簡単なモデルを使って「再現」をしてみましょう。

サスペンションモデル

車両のサスペンションを含む最も単純な物理運動表現として、車体をタイヤとサスペンションに単純化したモデルを導入します。

このモデルでは、タイヤのバネ定数 $ k_1 \text{$N/m$} $ や、ハブを起点に評価したホイールレート $ k_2 \text{$N/m$} $ を代表的なパラメータとして取り扱います。

システム方程式は

$$ \frac{d}{dt} \begin{pmatrix} x_1 \\ x_2 \\ \dot{x}_1 \\ \dot{x}_2 \end{pmatrix} = \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -\frac{k_1 + k_2}{m_t} & \frac{k_2}{m_t} & -\frac{c_2}{m_t} & \frac{c_2}{m_t} \\ \frac{k_2}{m_b} & -\frac{k_2}{m_b} & \frac{c_2}{m_b} & -\frac{c_2}{m_b} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ \dot{x}_1 \\ \dot{x}_2 \end{pmatrix} $$

GT7で発散が報告されているパラメータを元に下記の式からシミュレーションで使う仮の値を概算で設定します。

車体重量 $ m_b $ : 1,000 kg
タイヤ総重量 $ m_t $ : 100 kg
サスの固有振動数 $ f_2 $ : 1.5 Hz
サスの減衰比 $ \zeta_2 $ : 0.25

このほか、タイヤのバネ定数 $ k_1 $ をおよそ $ 2.0 \cdot 10^6 \text{N/m} $ とします。

システム方程式で使うパラメータはバネマスダンパー系の特性値の公式を使って概算値として求めます。

$$ \begin{aligned} k_2 &= m_b (2 \pi f_2)^2 \\ c_2 &= 2 m_b \zeta_2 (2 \pi f_2) \end{aligned} $$

物理シミュレーション

物理シミュレーションでは、現実世界の物理現象を数値的に再現するために、微分方程式を解く必要があります。これに数値積分法が使用されます。原理が簡単な前進オイラー法と、物理シミュレーションエンジンでポピュラーな半陰オイラー法を試します。

物理モデルが行列 $ A \in \mathrm{GL}_{2n}(\mathbb{R}) $ を使って、微分方程式

$$ \dot{x} = \frac{d}{dt} \begin{pmatrix} X \\ V \end{pmatrix} = A x $$

で与えられたとき、数値積分で逐次的にステップ $ n $ の $ x_n $ を求めます。

本件のモデルに対しては

$$ A = \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -2 \cdot 10^4 & 8.9 \cdot 10^2 & -4.7 \cdot 10 & 4.7 \cdot 10 \\ 8.9 \cdot 10 & -8.9 \cdot 10 & 4.7 & -4.7 \\ \end{pmatrix} $$

です。

前進オイラー法

前進オイラー法は、次の状態を現在の状態からタイムステップを使って計算するシンプルな数値積分法です。数式で表すと次のようになります。

$$ x_{n+1} = x_n + \Delta t A x_n $$

先程のサスペンションモデルをこの数値積分で解いてみましょう。

通常、ゲーム向けのシミュレーションではシミュレーションのタイムステップと画面のリフレッシュレートを一致させますから、GT7のリフレッシュレート120 Hzをひとつの水準にします。また、対比のため、1kHzでのシミュレーションもしてみましょう。

初期値として車両が加速トルクをうけて、準静的に1 cmスクワットし、これが開放された瞬間を取ります ($ x_0^\intercal= [0,-0.01,0,0] $ )。

120Hzの水準では車体の振幅が次第に拡大していき、ほんの3秒の間にはね飛ぶようになります(グラフ上段)。1kHzの水準ではダンパーの効果もあり、振幅がゼロに収束します。厳密解(実車体の動き)はこちらに近いはずです。

グラフ下段に示した通り、120Hzのシミュレーション条件ではシステムへ仕事の供給がないのに系のエネルギーが指数的に増えていきます。

半陰オイラー法

半陰オイラー法は、前進オイラー法の安定性を改善するため、速度の更新結果に基づいて位置を更新します。

$$ \begin{aligned} V_{n+1} &= V_n + \Delta t \dot{V_n} \\ X_{n+1} &= X_n + \Delta t V_{n+1} \\ &= X_n + \Delta t V_n + (\Delta t)^2 \dot{V_n} \end{aligned} $$

これを折り込み、前進オイラー法と同様に1ステップで計算するには…

$$ x_{n+1} = x_n + \Delta t (I + \Delta t \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} ) A x_n $$

この式で下のように $ \tilde{A} $ を取れば実装上は前進オイラー法と同じです。

$$ \tilde{A} = (I + \Delta t \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} ) A $$

これを使って修正した行列で同様にサスペンションモデルをシミュレーションすると120Hzでも発散しないことがわかります。

本シミュレーションの条件では発生しませんでしたが、より厳しいパラメータを与えると、この半陰オイラー法を採用しても発散します。GT7も、一般の物理シミュレーション同様に、この数値積分法を使っていると考えられますので、剛体同士の接触などが起きていたのだろうと思います。

数値積分の安定域

シンプルな常微分方程式

$$ \begin{aligned} \dot{y} &= \lambda y \\ y(0) &= 1 \\ \lambda & \in \{ \text{Re}(x)<0 | x \in \mathbb{C}\} \end{aligned} $$

を考えます。この常微分方程式の解は、$ y = exp(\lambda t) $ です。この方程式を前進オイラー法で解くことを考えると、タイムステップ $ \Delta t $ に対して、数値解は

$$ \begin{aligned} y_{n+1} &= y_n + \Delta t \lambda y_n \\ &= (1+\Delta t \lambda) y_n \\ &= (1+\Delta t \lambda)^n y_0. \end{aligned} $$

このとき、$ \lambda $ の定義から $ t \rightarrow \infty $ で $ y = 0 $ ですから、数値解の収束条件は $ | 1 + \Delta t \lambda | < 1 $ とわかります。

同様に

$$ B = \begin{pmatrix} \lambda_{1} & & \\ & \ddots & \\ & & \lambda_{2n} \end{pmatrix} $$

を考え、$ \dot{y} = By $ について考察すると、すべての $ \lambda_k $ が前記の収束条件を満たしている必要性があります。

また、仮に$ A $ が $ B = P^{-1}AP $ と対角化可能であれば、$ B $ の対角に並ぶ要素は $ A $ の固有値です。数値積分の収束について議論する場合は、Aの固有値に注目すれば良さそうです。

本ケースの考察

サスペンションモデルの $ A $ の固有値を求めると $ -23.7±141.8i, -2.2±9.0i $であることから、120Hz ($ \Delta t = \frac{1}{120} $) の条件で、(絶対値が大きい方の固有値で) $ | 1+ \frac{-23.7+141.8i}{120}| > 1 $ となり、収束条件を満たしません。

前進オイラー法では、収束条件を満たすには $ \Delta t < \frac{1}{433} $ (433 Hz以上)である必要があります。実際、1kHzの条件では収束していました。

次に半陰オイラー法について $ \tilde{A} $ の固有値を求めると $ -110.8±91.7i, -2.5±9.0i $です。収束条件は $ \tilde{A} $ について、前進オイラー法と同じ議論は成り立つことから、収束条件を確認すると、 $ | 1+ \frac{-110.8+91.7i}{120}| \approx 0.77 < 1 $ で、120Hzの数値積分でも収束条件を満たすことがわかります。