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

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

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

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

3 なぜ物理シミュレーションが必要か？ 入力に対する反応の多様さ →ゲーム世界の多様さ・楽しさ 入力の進化 アナログパッド，画像，力覚
入力に対する反応の多様さ →ゲーム世界の多様さ・楽しさ 入力の進化 アナログパッド，画像，力覚 ユーザの選択肢が増加 たくさんの反応を用意しなければならない 従来の延長 たくさんの手間と時間 フレーム問題 シミュレーション 自動的に多様でリアルな反応

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

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

6 物理シミュレータの選定 いくつもの物理シミュレータ いくつかの方式と特徴 物理シミュレータの中身を知る必要がある．
Havok， MathEngine, Tokamak （商用） Open Dynamics Engine, Springhead （オープンソース） いくつかの方式と特徴 スピードと精度のトレードオフ 何の精度を重視するか 物理シミュレータの中身を知る必要がある．

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

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

9 剛体の運動 運動方程式 ならば，速度一定・角運動量一定 m: 質量 I: 慣性テンソル f: 外力 r: 外力の作用点の位置 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 解析法１（関節） f A pA B pB 拘束： 運動方程式： からfを求められる

13 解析法２（抗力） fn 拘束： B 運動方程式: A 抗力は，反発だけ： pB pA
利点： 物体同士が侵入しない． 欠点： 遅い，跳ね返り係数を考慮していない．

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

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

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

17 ペナルティ法（面接触） 面で接触する場合 正確な接触力の計算 力は接触部分全体から発生 抗力・摩擦力の分布モデルを考える
静止摩擦←→動摩擦は一度に切り替わるとする 発生する力を接触部全体について積分 意外と簡単な式になる 解析法に比べ，とくに摩擦力が正確になる

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

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

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

21 ペナルティ法は高速 60 ペナルティ法（Springhead） 50 解析法（Open Dynamics Engine） 40 30
計算時間[ms] 20 10 5 10 15 ブロック数

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

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

24 解析法の高速化 Featherstoneの方法 とても速いが， 巨大行列を作らない 葉から根に向かって１つずつ拘束力を求める
構造を変えるのに多少時間がかかる 輪があると計算できない →抗力計算には向かない

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

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

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

28 解析法の高速化 ループを含めた高速化方法 山根克 98 「構造可変なリンク系の動力学計算法とヒューマンフィギュアの運動計算」 Open HRP (Humanoid Robotics Platform)に 使われているようです （なぜかOpen HRPは非常に重い）

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

30 Springheadの選択 接触力計算はペナルティ法 ロボットなど関節を持つ物はFeatherstone法
高速性，安定性 接触領域に分布モデルを考える→摩擦の精度 多少剛体同士が侵入するが，あきらめる． ロボットなど関節を持つ物はFeatherstone法 構造が変化しないので非常に高速． ループのある機構は少ない． ペナルティ法とは簡単に組み合わせられる． シミュレーション法・高速化法の特性を考慮して， シミュレータを選んでください．

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

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

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

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

35 PD制御の性質と調節 k b バネ係数kとダンパ係数bでPD制御の性質が決まる バネ・ダンパと考えられるので， m バネ係数 k
大きいほど早く目標に近づこうとする 小さいほど柔らかい ＝力が加わったときに動きやすい ダンパ係数b 負だと振動がどんどん大きくなる 0だと振動が止まらない 大きいほど振動しにくい 大きいほど撃力が伝わる kが大きいほど早く動く． b>0ならば，収束する（振動が収まる）． 　 連続系ではこうなるが，シミュレータは離散系 そううまくはいかない． m 目標位置 k b

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 k 2 b=k 実際に良く使うのはこの辺

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

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

42

43

44

45

46

47

48

49

50

51

52

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

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

55 凸形状 凸形状の便利な性質 凸形状 非凸形状 凸形状 距離が極小となる点が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) 凸形状 非凸形状

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

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

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