実時間物理シミュレーション技術 東京工業大学精密工学研究所 長谷川晶一.

Presentation on theme: "実時間物理シミュレーション技術 東京工業大学精密工学研究所 長谷川晶一."— Presentation transcript:

1 実時間物理シミュレーション技術 東京工業大学精密工学研究所 長谷川晶一

2 内容 物理シミュレーションの意義 シミュレータの特徴と選び方 制御 デモ なぜ，今物理シミュレーションなのか？
物理シミュレーションでできること シミュレータの特徴と選び方 接触力計算法，高速化法 制御 シミュレータ内の物体の動かし方 デモ 我々が開発している物理シミュレータ付きVR開発環境Springheadと力覚インタフェースSPIDARのデモ

3 なぜ物理シミュレーションが必要か？ 入力に対する反応の多様さ →ゲーム世界の多様さ・楽しさ 入力の進化
入力に対する反応の多様さ →ゲーム世界の多様さ・楽しさ 入力の進化 アナログパッド，画像，力覚 → 選択肢の増加 たくさんの反応を用意しなければならない 従来の延長 場合分けの爆発 たくさんの手間と時間 シミュレーション 自動的に多様でリアルな反応 従来の新技術と同じこと 多量の2次元画像 → ３DCG ストーリの分岐 → 束構造に 例を挙げても良いかも アナログが活きてるゲーム フレーム問題解消例 オトギリソウ・かまいたちの夜 → 木構造を束構造にして解消

4 物理シミュレーションが役立つ例 動きの生成 例：歩行 当たり判定 ダメージ計算 効果音
動きの生成 例：歩行 シミュレーションなし： 足の動きと体の動きを別々に計算 シミュレーションあり： 足を動かすと自然と体が動く 当たり判定 なし： 当たり判定専用計算 あり： ポリゴンモデルの接触位置が毎ステップ求まる ダメージ計算 なし： 当たったもの，当たり方毎に値を用意 あり： 加わった力の大きさからダメージを計算 効果音 なし： イベントごとに音を用意 あり： 加わった力，材質から自動的に音を発生 DEMO カブトムシの例 物理なしとあり 効果音の自動生成の例

5 物理シミュレーションの役割 ゲーム世界の物体が 現実世界と同様に自律的に動く 計測と制御 従来のゲームでも物理は使われていた？
従来は，特別に，物理法則を作りこんでいた． 物理シミュレータでは，物体はすべて 物理法則にしたがう．

6 物理シミュレータを使うには 物理シミュレータの中身を知る 物理シミュレータの選定 パフォーマンスを引き出す 自分のゲームでは何が大事か？
Havok， MathEngine? （商用）, Tokamak（無料） Open Dynamics Engine, Springhead （オープンソース） パフォーマンスを引き出す スピードと精度のトレードオフ 精度同士のトレードオフ ボトルネックの特定 TOKAMAK・Open Dynamics Engine・Springheadをデモする

7 物理シミュレーションとは … 物理法則(現実世界)は微分方程式で記述できる シミュレーションは微分方程式の数値解の計算 m たとえば
質点の運動 剛体の運動 流体の運動 運動方程式： m f 差分方程式にすると： x 順番に求めて行く： … x(0) x(Dt) x(2Dt)

8 剛体運動のシミュレーション 剛体の運動の話だけにします． 剛体 硬いもの，変形しないもの 剛体だと考えてよいものが多い 剛体でないもの
積み木，ボール，ロボット，人体・・・ 剛体でないもの スポンジ，粘土，水・・・

9 剛体の運動 運動方程式 ならば，速度一定・角運動量一定 m: 質量 I: 慣性テンソル f: 外力 v: 速度 ω:角速度
(すべて絶対座標系) 運動方程式 ならば，速度一定・角運動量一定

10 剛体に働く力 kx mg fn ft 重力→ f=mg… 定数 バネ→ f=kx… 位置に比例 拘束力 拘束力の計算が難しい
力の大きさは不明 剛体同士の位置・速度関係が決まっている 蝶番：2物体の相対位置が一定 抗力：2物体が互いに侵入しない 静止摩擦力：物体が滑らない 拘束力の計算が難しい kx mg fn ft

11 拘束力の求め方 例：球関節の拘束． ２物体が 点pAと点pBで繋がっている B 拘束の式： pA = pB
この式を満たすように，関節に働く力を求める 解析法：運動方程式と拘束の式を連立させて解く David Baraff 89-93など Havok，Tokamak, Open Dynamics Engine ペナルティ法：拘束違反に応じた力を加える 昔からいろいろな用途に使われてきた Springhead 接触判定時に拘束違反の量（侵入量）を調べる必要がある B pB pA A

12 解析法１（関節） B 拘束： f rB A pB pA rA 運動方程式： 計算量はo(n3)

13 解析法２（抗力） fn 拘束： B 運動方程式: A 抗力は，反発だけ： pB pA -a/b f
抗力は，反発だけ： これを満たす f を見つける→ 線形計画法,2次計画法・・・ 利点： 物体同士が侵入しない． 欠点： 遅い，跳ね返り係数を考慮していない．

14 解析法３（摩擦力） クーロンの摩擦モデル 拘束： 運動方程式: B 抗力は，反発力: 摩擦力の条件： A
fN 抗力は，反発力: 摩擦力の条件： fS A 場合分けを無くすため とすることが多い

15 ペナルティ法（抗力） . 侵入量 d，相対速度 d 拘束を解かない． 拘束を侵した分だけ罰(ペナルティ)として力を加える．
繰り返すうちに拘束が満たされる．．．はず． 力が直接決まるので，計算量はo(n) 接触部にバネとダンパを入れたと考えられる 侵入量 d，相対速度 d . バネ ダンパ 利点： 高速，跳ね返り係数を考慮できる． 欠点： 物体同士が多少侵入してしまう．

16 ペナルティ法（摩擦力） 静止摩擦：ずれに比例した力を加える 接触の履歴を利用 クーロンの摩擦モデルをそのまま利用できる 静止摩擦 動摩擦

17 ペナルティ法（面接触） 面で接触する場合 正確な接触力の計算 力は接触部分全体から発生 分布バネダンパモデルを考える
発生する力を接触部全体について積分 意外と簡単な式になる 解析法に比べ，とくに摩擦力が正確 接触部分の形を求める必要

18 計算量と高速化 ここまでで，物理の話は終わりです． ここからは，リアルタイム動作に必要な， 計算量と高速化について話します． 関節を持つ剛体
剛体同士の接触 のシミュレーションができるようになりました． ここからは，リアルタイム動作に必要な， 計算量と高速化について話します．

19 方程式を解くので，行列の次数nに対して，計算量はo(n3)
解析法の計算量と接触数 接触している 関節で繋がれる 多物体が 場合多くの拘束が働く． f1n f2n 行列Aの次元＝抗力 fi の数 方程式を解くので，行列の次数nに対して，計算量はo(n3) 立方体の面接触１つ：４次元 １０個つむと：４０次元

20 解析法の高速化 Aを早く解く工夫 ODEの場合 巨大な行列Aを作らず， 2物体単位で計算する 繰り返し計算する
繰り返し回数が少ないと精度が落ちる (ペナルティ法に近い） 巨大行列Aを解いた場合 ２体単位で，繰り返し 解く場合 Open Dynamics Engineのドキュメントより

21 構造を利用した，解析法の高速化 Articulated Body 構造に特徴 多数の剛体を関節でつないだもの 動物や人間の体，ロボットなど
剛体が輪になっていない→木構造 輪になっている例 全部繋がっている例

22 解析法で求めると．．． D C B A 3 2 1 並べてみると 0だらけ，すかすか，sparse行列 普通に解くより早い方法があるのでは？
１の拘束の式： Bの・・・・・・・・・・・・・・・・ 並べてみると 0だらけ，すかすか，sparse行列 普通に解くより早い方法があるのでは？

23 Featherstoneの方法 巨大行列を普通に解かない 根から葉に向かって１つずつ拘束力を求める とても速い O(n)

24 Featherstoneの方法 全体を一つの剛体だと考え，加速度を求める． AとBから先の2つの剛体と考え BとCから先 〃
関節1に働く力を求める． Bの加速度を求める BとCから先 〃 1 2 3 A B C D 1 2 3 A B C D

25 Featherstoneの方法 D C B A 先ほどのsparse行列で考えると， 3 2 1 1. 葉から根に向かってMを合成
2.加速度 を求める．

26 Featherstoneの方法 3.根から葉に向かって， を求める

27 ペナルティ法は高速 バネダンパモデルから力が求まる 5 10 15 20 30 40 50 60 計算時間[ms] ブロック数
5 10 15 20 30 40 50 60 計算時間[ms] ブロック数 ペナルティ法（Springhead） 解析法（Open Dynamics Engine）

28 解析法とペナルティ法 解析法 ペナルティ法 シミュレーション法・高速化法の特性を考慮して， シミュレータを選んでください．
1ステップで拘束と評価関数を満たす力を計算 ステップが大きく取れる．1ステップの計算は多い． 摩擦や跳ね返り係数などは難しい → 動きの精度低 or 評価関数化難，計算量増大 繰り返し計算による高速化：1ステップに何度も計算 ペナルティ法 バネダンパモデル→1ステップの計算は速い． 侵入量∝Δt なので，ステップ数が多くなる． 摩擦・跳ね返り係数なども簡単にモデル化できる． シミュレーション法・高速化法の特性を考慮して， シミュレータを選んでください．

29 Springhead 我々が開発している物理シミュレータ 開発の動機 力覚インタフェースに使いたい
超高速更新(>300Hz)，安定性重視 シミュレータでロボット対戦（ロボコン）をやりたい ２台以上 フィールドに障害物 押し合い → 摩擦が重要

30 Springheadの選択 接触力計算はペナルティ法 ロボットなど関節を持つ物はFeatherstone法 高速性，安定性 摩擦の精度
多少剛体同士が侵入するが，気にしない ロボットなど関節を持つ物はFeatherstone法 構造が変化しないので非常に高速． ループのある機構は少ない． ペナルティ法とは簡単に組み合わせられる．

31 ゲームの作成 シミュレータの限界 限界を見極めてモデリングする必要がある 安定性の限界 計算速度の限界 極端に重いものと軽いもの
衝突時に速度が大きくなりすぎる 質量と慣性モーメントの比率がおかしいもの 回転速度が極端に大きくなることがある 極端に大きさが異なるもの 衝突判定の精度が問題となる 計算速度の限界 物体数，ポリゴン数，時間刻み デモしながら

32 ゲームの作成 物理以外に必要なもの 物理シミュレーションを活かして作ると効果的 人物，動物，車などの動き ゲームのルール ダメージの計算
物体の位置や速度で判定 ダメージの計算 接触力や関節に働く力の計測 物理シミュレーションを活かして作ると効果的 人物，動物，車などの動き 自然の動きではなく，意図を持った動き → コントロール＝制御が必要

33 物体の制御 シミュレータ内の物体は，速度・慣性を持つ． 自動的に運動する 強制的に位置を指定 不自然な動き，非常に大きな力が発生
シミュレーションの意味がない やわらかい動き 物理を無視せず，実世界と同じく力を加えて動かす． ＝制御をする DEMO 強制位置指定と制御の比較

34 PD制御 f m xg x0 目標位置に質点を持っていくには？ 誤差に応じた力を加えてやる 例えば f=k(xg – x)では？
振動し続ける 　　減衰し，振動が止まる f m x0 xg 初期位置 目標位置 誤差に比例:バネ Proportional制御 微分に比例：ダンパ Differential制御

35 PD制御の性質と調節 k b バネ係数kとダンパ係数bでPD制御の性質が決まる バネ・ダンパと考えられるので， m
大きいほど早い．小さいほど柔らかい． ダンパ係数b 負だと振動がどんどん大きくなる 0だと振動が止まらない 大きいほど振動しにくい b>0ならば，収束する（振動が収まる）． kが大きいほど早く動く． 振動の周期は シミュレーションは時間が離散 少なくとも周期＜2Δtは無理 m 目標位置 k b Δt

36 シミュレータ上でのPD制御 シミュレータ上でのPD制御 時間が離散 m x v 位置 速度 目標位置
目標位置 この前に，サンプリング定理の話を入れたほうが良いかも？

37 シミュレータ上でのPD制御 安定性の確認 行列A An →0 ならば，このPD制御は安定になる 安定：そのうち目標位置でとまる

38 シミュレータ上でのPD制御 安定性の確認（つづき） |Aの固有値| < 1 ならば， An →0 なので， Aの固有値を求めてみる

39 シミュレータ上でのPD制御 ~ ~ 安定な k,b の範囲， 早く静止する k,b を求める ~ ~ ~ b b=k 2 ~ k
実際に良く使うのはこの辺 ~ k

40 無次元化について この数値は何？ この数値は質量m,ステップΔtに依らない？

41 PD制御の実験 ～ ～ いろいろな k, b での挙動をご覧ください ～ ～ ～ b b=k 2 1.5 ～ k 2 0.5

42 ペナルティ法とPD制御 . ペナルティ法もバネ・ダンパモデル バネ・ダンパ係数は？ 侵入量 r，相対速度 r
では多物体が接触したときに安定しない． Springheadでは, を使用 侵入量 r，相対速度 r . バネ ダンパ

43 まとめ 物理シミュレータの役割 シミュレータの中身と選定 やわらかく動かすための制御 ゲーム世界の多様性，多様な反応の実現
解析法とペナルティ法 私はペナルティ法が好きですが，目的次第です． 高速化法 行列Aをいかに早く解くか（解かずに済ますか） やわらかく動かすための制御 PD制御 バネ・ダンパ係数の決め方，無次元化

44 謝辞 このような機会を与えて下さった， 星野准一先生，新清士様，IGDAの皆様 一緒に開発したSpringhead開発チーム
田崎勇一，岡田直樹，市川宙，白井暁彦， 藤井伸旭，田上信一郎 ご清聴ありがとうございました．

45 参照 デモとすべてのソースはSpringheadのWebで公開 この資料や参考文献リストもおく予定 ご意見，ご質問，お問い合わせを WebのWiki・掲示板・バグトラッカー でお待ちしております ここまで終わってから，Springheadのデモをいくつか見てもらう． 同時にSPIDARも見てもらう

46

47

48

49 衝突 実世界の物体は 互いに侵入しない． 跳ね返る． 再現するためには 衝突検知 剛体間に働く力の決定

50 階層化による高速化 Model Hierarchy: Bouding Volumeが球で子が2つの場合の例：
各ノードが子ノードを含む Bounding Volume を持つ． Bouding Volume は，球，直方体などで表される． 階層構造の一番したのノード(葉)が多面体モデルを持つ． Bouding Volumeが球で子が2つの場合の例：

51 階層化による高速化 上の階層から順に判定 枝刈りできる →高速

52 衝突検知 衝突の検出 衝突しているかどうかすべての剛体について調べる． 剛体の形状は，多面体で表現されている．
衝突検知を簡単にするため，凸形状に分割しておく

53 凸形状 凸形状の便利な性質 凸形状 非凸形状 凸形状 距離が極小となる点が1点 最近傍点が簡単に求まる GJK algorithm
E. G. Gilbert, D. W. Johnson and S. S. Keerthi A Fast Procedure for Computing the Distance between Complex Objects in Three-Dimensional Space (1988) 凸形状 非凸形状

54 衝突検知（GJK） 凸形状A上の点から，凸形状B上の点へのベクトルを 原点を始点に並べると ベクトルの終点の集合も凸形状になる

55 衝突検知（GJK2） １ ２ ３ ４ V0 : 凸形状内の任意の点 Wi :OViとOWiの内積が最小の点
Vi :三角形Wi-2 Wi-1 Wi内の点で原点に 一番近い点 ２ ３ ４

56 接触解析 凸形状の便利な性質２ 凸形状の交差部分も凸形状 交差部分の形状が簡単に求まる
D. E. Muller and F.P.Preparata: “Finding the intersection of two convex” (1978) Half space representation

57 接触解析 Finding the intersection of two convex
Preparation: Dual Transformation Dual transformation transform a face into a vertex and a vertex into a face. Dual transformation’s dual transformation is original facet. Muller’s algorithm use dual transformation, So I will explain dual transformation at first. Dual transformation transform a vertex into a plane and a plane into a vertex. The dual transformation of a dual transformation is identical transformation. Dual Transformation O O

58 接触解析 Finding the intersection of two convex(2) Half space
representation Dual transform Muller’s algorithm represents two convex by half space. Half space representation is common area of half spaces, which represented by planes. So, the summation of the two half space representations results in the intersection of the two convex polyhedrons. But, the representation has redundant planes. So we have to find the minimum set of the planes which represents the intersection. To find the minimum set, the algorithm uses dual transformation. Dual transformation transform the planes into vertices. And the problem of the minimum set of the planes results in a convex hull problem, which find the minimum convex from the vertices. There is fast algorithm to solve convex hull problem. We use quick hull algorithm for this. After finding the convex hull, Muller’s algorithm uses dual transform again. Then we get the vertices and faces of the intersection. Vertex of intersection Dual transform

59 接触量の積分 ＝ ー o The intersection is convex polyhedron
Intersection = upper bound – lower bound We can integrate penalty by each face. Intersection Contact normal ＝ ー Now, we got the intersection of the two polyhedrons. So the next, we integrate the penalty over the intersection. The volume of the intersection can be represented by upper bound minus lower bound. And the boundary is composition of triangles. So we can integrate penalty one by one for each triangles. d+db d+db d d+da b o d+da a

60 The computations for the results are simple
力とトルクの計算 Each point of intersection generates penalty and it’s moment. Penalty and its moment from a triangle can be represented by: Penalty Moment of penalty We thought that each point of the intersection generate penalty force and its moment. We integrate the penalty and its moment for triangles. We can symbolically integrate them. And the result is simple. So we can integrate the penalty very fast. The computations for the results are simple

61

62 Featherstoneの方法 リンクiとより葉側のリンク全体の運動方程式が， のように，書ける．
Mi: リンクiとより葉側のリンクの 質量・慣性による項 Zi: 外力，重力，コリオリ力による項 Mi,Ziは加速度 によらない 葉のMiは既知なので， 葉から根に向かってMiを 計算する． 根ではfiは0なので，運動方程式から加速度 を求める． 根から葉に向かって， を求める 葉のリンク リンクi 葉のリンク 根側のリンク リンクi とより葉側のリンク 根のリンク