Human interface Section, P&I Lab, Titech
実時間物理シミュレーション技術
東京工業大学精密工学研究所
長谷川晶一
Human interface Section, P&I Lab, Titech
内容
物理シミュレーションの意義
なぜ,今物理シミュレーションなのか?
物理シミュレーションでできること
シミュレータの特徴と選び方
接触力計算法,高速化法
制御
シミュレータ内の物体の動かし方
デモ
我々が開発している物理シミュレータ付きVR開発環境
Springheadと力覚インタフェースSPIDARのデモ
なぜ物理シミュレーションが必要
か?
Human interface Section, P&I Lab, Titech
入力に対する反応の多様さ
→ゲーム世界の多様さ・楽しさ
入力の進化
アナログパッド,画像,力覚 → 選択肢の増加
たくさんの反応を用意しなければならない
 従来の延長
 シミュレーション
場合分けの爆発
たくさんの手間と時間
自動的に多様でリアルな反応
従来の新技術と同じこと
 多量の2次元画像
 ストーリの分岐
→
→
3DCG
束構造に
Human interface Section, P&I Lab, Titech
物理シミュレーションが役立つ例
動きの生成
例:歩行
シミュレーションなし: 足の動きと体の動きを別々に計算
シミュレーションあり: 足を動かすと自然と体が動く
当たり判定
なし: 当たり判定専用計算
あり: ポリゴンモデルの接触位置が毎ステップ求まる
ダメージ計算
なし: 当たったもの,当たり方毎に値を用意
あり: 加わった力の大きさからダメージを計算
効果音
なし: イベントごとに音を用意
あり: 加わった力,材質から自動的に音を発生
Human interface Section, P&I Lab, Titech
物理シミュレーションの役割
ゲーム世界の物体が
現実世界と同様に自律的に動く
計測と制御
従来のゲームでも物理は使われていた?
従来は,特別に,物理法則を作りこんでいた.
物理シミュレータでは,物体はすべて
物理法則にしたがう.
Human interface Section, P&I Lab, Titech
物理シミュレータを使うには
物理シミュレータの中身を知る
物理シミュレータの選定
自分のゲームでは何が大事か?
Havok, MathEngine? (商用), Tokamak(無料)
Open Dynamics Engine, Springhead (オープンソース)
パフォーマンスを引き出す
スピードと精度のトレードオフ
精度同士のトレードオフ
ボトルネックの特定
Human interface Section, P&I Lab, Titech
物理シミュレーションとは
物理法則(現実世界)は微分方程式で記述できる
たとえば
質点の運動
剛体の運動
流体の運動
f  mx
 
mv  f , Iω
シミュレーションは微分方程式の数値解の計算
f  mx
運動方程式:
f
m
f t
差分方程式にすると: x (t  t )  x (t ) 
m
x(t  t )  x(t )  x (t )t
順番に求めて行く:
x(t )  x (t )t  x(0)
x(2t )  ...
x(3t )  ......
x
…
x(0) x(t) x(2t)
Human interface Section, P&I Lab, Titech
剛体運動のシミュレーション
剛体の運動の話だけにします.
剛体
硬いもの,変形しないもの
剛体だと考えてよいものが多い
積み木,ボール,ロボット,人体・・・
剛体でないもの
スポンジ,粘土,水・・・
Human interface Section, P&I Lab, Titech
剛体の運動
v: 速度
ω:角速度
m: 質量
I: 慣性テンソル
f: 外力
(すべて絶対座標系)
運動方程式
mv  f
 τ
Iω
mv(t  t )  mv(t )  ft
I(t  t ) ω(t  t )  I(t ) ω(t )  τt
f  0, τ  0 ならば,速度一定・角運動量一定
Human interface Section, P&I Lab, Titech
剛体に働く力
重力→ f=mg… 定数
バネ→ f=kx… 位置に比例
拘束力
mg
kx
力の大きさは不明
剛体同士の位置・速度関係が決まっている
 蝶番:2物体の相対位置が一定
 抗力:2物体が互いに侵入しない
 静止摩擦力:物体が滑らない
拘束力の計算が難しい
fn
ft
Human interface Section, P&I Lab, Titech
拘束力の求め方
例:球関節の拘束.
2物体が 点pAと点pBで繋がっている
拘束の式: pA = pB
A
pA
この式を満たすように,関節に働く力を求める
解析法:運動方程式と拘束の式を連立させて解く
 David Baraff 89-93など
 Havok,Tokamak, Open Dynamics Engine
ペナルティ法:拘束違反に応じた力を加える
 昔からいろいろな用途に使われてきた
 Springhead
 接触判定時に拘束違反の量(侵入量)を調べる必要がある
pB
B
Human interface Section, P&I Lab, Titech
解析法1(関節)
B
f
A
rB
pA pB
 B  p
 A  0
拘束: p
(rA,  θ A  (p A  rA )  θ A  (p A  rA ))
 (r  θ  (p  r )  θ  (p  r ))  0
B,
rA
IA
mB
M
B
B
B
B
B
Br  b  0
運動方程式:
 mA






B
1
 rA  

  



 θ A    (p A  rA )  
 r   
f
1
 B  







I B  θ B   (p B  rB )  
r 
C f
Af  b  0
計算量はo(n3)
Human interface Section, P&I Lab, Titech
解析法2(抗力)
 B  p
 A )  n  0
拘束: (p
B
pB
fn
運動方程式: Mr  Cfn
af  b  0
pA
A
(Br  b)  n  0
抗力は,反発だけ:
これを満たす f を見つける→
利点:
欠点:
0 -a/b
f 0
線形計画法,2次計画法・・・
物体同士が侵入しない.
遅い,跳ね返り係数を考慮していない.
f
Human interface Section, P&I Lab, Titech
解析法3(摩擦力)
クーロンの摩擦モデル
f S  0 f N (静止摩擦)
f d   f N (動摩擦)
拘束:
B
fN
抗力は,反発力:
fS
A
運動方程式:
 B  p
 A )  n  0
(p
Mr  Cf
f n  0
摩擦力の条件: f  n / | f |  / 1   2
場合分けを無くすため 0   とすることが多い
Human interface Section, P&I Lab, Titech
ペナルティ法(抗力)
 拘束を解かない.
Af  b  0
 拘束を侵した分だけ罰(ペナルティ)として力を加える.
.
侵入量 d,相対速度 d
f  Kd  Bd
バネ ダンパ
 繰り返すうちに拘束が満たされる...はず.
力が直接決まるので,計算量はo(n)
接触部にバネとダンパを入れたと考えられる
利点:
欠点:
高速,跳ね返り係数を考慮できる.
物体同士が多少侵入してしまう.
Human interface Section, P&I Lab, Titech
ペナルティ法(摩擦力)
静止摩擦:ずれに比例した力を加える
接触の履歴を利用
クーロンの摩擦モデルをそのまま利用できる
f d   f N (動摩擦)
f S  0 f N (静止摩擦)
静止摩擦
動摩擦 静止摩擦
Human interface Section, P&I Lab, Titech
ペナルティ法(面接触)
面で接触する場合
力は接触部分全体から発生
正確な接触力の計算
分布バネダンパモデルを考える
発生する力を接触部全体について積分
意外と簡単な式になる
解析法に比べ,とくに摩擦力が正確
接触部分の形を求める必要
Human interface Section, P&I Lab, Titech
計算量と高速化
ここまでで,物理の話は終わりです.
関節を持つ剛体
剛体同士の接触
のシミュレーションができるようになりました.
ここからは,リアルタイム動作に必要な,
計算量と高速化について話します.
Human interface Section, P&I Lab, Titech
解析法の計算量と接触数
接触している
多物体が
場合多くの拘束が働く.
関節で繋がれる
f1n
f2n
Af  b  0
f 0
行列Aの次元=抗力 fi の数
方程式を解くので,行列の次数nに対して,計算量はo(n3)
立方体の面接触1つ:4次元
10個つむと:40次元
Human interface Section, P&I Lab, Titech
解析法の高速化
 Aを早く解く工夫
 ODEの場合
Af  b  0
f 0
 巨大な行列Aを作らず,
2物体単位で計算する
 繰り返し計算する
 繰り返し回数が少ないと精度が落ちる (ペナルティ法に近い)
巨大行列Aを解いた場合
2体単位で,繰り返し
解く場合
Open Dynamics Engineのドキュメントより
Human interface Section, P&I Lab, Titech
構造を利用した,解析法の高速化
Articulated Body
多数の剛体を関節でつないだもの
動物や人間の体,ロボットなど
構造に特徴
剛体が輪になっていない→木構造
輪になっている例 全部繋がっている例
Human interface Section, P&I Lab, Titech
解析法で求めると...
3
2
B1
A
D
Aの運動方程式: M ArA  CA1 (f1 )
1の拘束の式: B1ArA  B1BrB  b1  0
Bの・・・・・・・・・・・・・・・・ ここで, rやf は6次元ベクトル
r  vx
v y
vz
 x
 y
 z T
並べてみると
MA





B
 1A



C A1
MB
CB1
MC
CC2
MD
B1B
B 2B
B 2C
B 3C
CB2
B 3D
 rA   0 
  

 rB   0 
CC3  rC   0 
  

C D3  rD    0 
 f    b 
 1   1 
 f 2    b 2 
  

f

b
 3   3 
0だらけ,すかすか,sparse行列
普通に解くより早い方法があるのでは?
Human interface Section, P&I Lab, Titech
Featherstoneの方法
巨大行列を普通に解かない
根から葉に向かって1つずつ拘束力を求める
とても速い O(n)
Human interface Section, P&I Lab, Titech
Featherstoneの方法
1. 全体を一つの剛体だと考え,加
速度を求める.
 剛体Aの加速度が求まる
3
2
D
B
1
A
2. AとBから先の2つの剛体と考え
1. 関節1に働く力を求める.
2. Bの加速度を求める
3
2
B
1
3. BとCから先 〃
A
D
Human interface Section, P&I Lab, Titech
Featherstoneの方法
先ほどのsparse行列で考えると,
3
2
B1
A
D
1. 葉から根に向かってMを合成
MA





B
 1A



C A1
MB
CB1
MC
CC2
MD
B1B
B 2B
B 2C
B 3C
CB2
B 3D
2.加速度rA を求める.
 rA   0 
  



r
0
 B  

CC3  rC   0 
  

C D3  rD    0 
 f    b  f 3を消去
 1   1 
 f 2    b 2 
  

f

b
 3   3  rDを消去
Human interface Section, P&I Lab, Titech
Featherstoneの方法
3.根から葉に向かって,fi , r を求める
MA





B
 1A



C A1
MB
CB1
MC
CC2
MD
B1B
B 2B
B 2C
B 3C
CB2
B 3D
 rA   0 
  

 rB   0 
CC3  rC   0 
  

C D3  rD    0 
 f    b 
 1   1 
 f 2    b 2 
  

f

b
 3   3 
rAが求まったので,
f1を求める
Human interface Section, P&I Lab, Titech
ペナルティ法は高速
バネダンパモデルから力が求まる
f   kx  bx
60
計算時間[ms]
50
ペナルティ法(Springhead)
解析法(Open Dynamics Engine)
40
30
20
10
0
0
5
10
ブロック数
15
Human interface Section, P&I Lab, Titech
解析法とペナルティ法
解析法
1ステップで拘束と評価関数を満たす力を計算
ステップが大きく取れる.1ステップの計算は多い.
摩擦や跳ね返り係数などは難しい
→ 動きの精度低 or 評価関数化難,計算量増大
繰り返し計算による高速化:1ステップに何度も計算
ペナルティ法
バネダンパモデル→1ステップの計算は速い.
侵入量∝Δt なので,ステップ数が多くなる.
摩擦・跳ね返り係数なども簡単にモデル化できる.
シミュレーション法・高速化法の特性を考慮して,
シミュレータを選んでください.
Human interface Section, P&I Lab, Titech
Springhead
我々が開発している物理シミュレータ
開発の動機
力覚インタフェースに使いたい
超高速更新(>300Hz),安定性重視
シミュレータでロボット対戦(ロボコン)をやりたい
2台以上
フィールドに障害物
押し合い → 摩擦が重要
Human interface Section, P&I Lab, Titech
Springheadの選択
接触力計算はペナルティ法
高速性,安定性
摩擦の精度
多少剛体同士が侵入するが,気にしない
ロボットなど関節を持つ物はFeatherstone法
構造が変化しないので非常に高速.
ループのある機構は少ない.
ペナルティ法とは簡単に組み合わせられる.
Human interface Section, P&I Lab, Titech
ゲームの作成
シミュレータの限界
限界を見極めてモデリングする必要がある
安定性の限界
極端に重いものと軽いもの
• 衝突時に速度が大きくなりすぎる
質量と慣性モーメントの比率がおかしいもの
• 回転速度が極端に大きくなることがある
極端に大きさが異なるもの
• 衝突判定の精度が問題となる
計算速度の限界
物体数,ポリゴン数,時間刻み
Human interface Section, P&I Lab, Titech
ゲームの作成
物理以外に必要なもの
ゲームのルール
• 物体の位置や速度で判定
ダメージの計算
• 接触力や関節に働く力の計測
物理シミュレーションを活かして作ると効果的
人物,動物,車などの動き
自然の動きではなく,意図を持った動き
→ コントロール=制御が必要
Human interface Section, P&I Lab, Titech
物体の制御
シミュレータ内の物体は,速度・慣性を持つ.
自動的に運動する
強制的に位置を指定
不自然な動き,非常に大きな力が発生
シミュレーションの意味がない
やわらかい動き
物理を無視せず,実世界と同じく力を加えて動かす.
=制御をする
Human interface Section, P&I Lab, Titech
PD制御
 目標位置に質点を持っていくには?
誤差に応じた力を加えてやる
例えば f=k(xg – x)では?
mx  k ( xg  x) 
m
x0
f
xg
初期位置 目標位置
x  ( x0  xg ) cos( k / m t )  xg
運動方程式: f  mx
振動し続ける
誤差に比例:バネ
Proportional制御
 f  k ( xg  x)  b( x g  x ) では?
x  ( x0  xg ) e (b / 2)t cos( k / m t )  xg
減衰し,振動が止まる
微分に比例:ダンパ
Differential制御
Human interface Section, P&I Lab, Titech
PD制御の性質と調節
 バネ係数kとダンパ係数bでPD制御の性質が決まる
 バネ・ダンパと考えられるので,
k
 バネ係数 k
 大きいほど早い.小さいほど柔らかい.
 ダンパ係数b
 負だと振動がどんどん大きくなる
 0だと振動が止まらない
 大きいほど振動しにくい
目標位置
m
b
x  ( x0  xg ) e(b / 2)t cos( k / m t )  xg
 b>0ならば,収束する(振動が収まる).
 kが大きいほど早く動く. 振動の周期は 2 k / m
 シミュレーションは時間が離散 少なくとも周期<2Δtは無理
Δt
Human interface Section, P&I Lab, Titech
シミュレータ上でのPD制御
シミュレータ上でのPD制御
時間が離散
m
位置 x
速度 v
0
目標位置
Human interface Section, P&I Lab, Titech
シミュレータ上でのPD制御
安定性の確認
行列A
An →0 ならば,このPD制御は安定になる
安定:そのうち目標位置でとまる
Human interface Section, P&I Lab, Titech
シミュレータ上でのPD制御
安定性の確認(つづき)
|Aの固有値| < 1 ならば, An →0 なので,
Aの固有値を求めてみる
Aの固有値
Human interface Section, P&I Lab, Titech
シミュレータ上でのPD制御
~~
安定な k,b の範囲, 早く静止する k,b を求める
~
b
~ ~
b=k
2
実際に良く使うのはこの辺
~
k
Human interface Section, P&I Lab, Titech
無次元化について

この数値は何?
この数値は質量m,ステップΔtに依らない?
Human interface Section, P&I Lab, Titech
PD制御の実験
~~
いろいろな k, b での挙動をご覧ください
~
~~
b
b=k
2
1.5
~
0.5
2
k
Human interface Section, P&I Lab, Titech
ペナルティ法とPD制御
ペナルティ法もバネ・ダンパモデル
.
侵入量 r,相対速度 r
f  Kr  Br  c
バネ・ダンパ係数は?
Springheadでは,
バネ ダンパ
では多物体が接触したときに安定しない.
を使用
Human interface Section, P&I Lab, Titech
まとめ
物理シミュレータの役割
ゲーム世界の多様性,多様な反応の実現
シミュレータの中身と選定
解析法とペナルティ法
 私はペナルティ法が好きですが,目的次第です.
高速化法
 行列Aをいかに早く解くか(解かずに済ますか)
やわらかく動かすための制御
PD制御
バネ・ダンパ係数の決め方,無次元化
Human interface Section, P&I Lab, Titech
謝辞
このような機会を与えて下さった,
星野准一先生,新清士様,IGDAの皆様
一緒に開発したSpringhead開発チーム
田崎勇一,岡田直樹,市川宙,白井暁彦,
藤井伸旭,田上信一郎
ご清聴ありがとうございました.
Human interface Section, P&I Lab, Titech
参照
デモとすべてのソースはSpringheadのWebで公開
http://springhead.info/
この資料や参考文献リストもおく予定
ご意見,ご質問,お問い合わせを
[email protected]
[email protected]
WebのWiki・掲示板・バグトラッカー
でお待ちしております
Human interface Section, P&I Lab, Titech
Human interface Section, P&I Lab, Titech
Human interface Section, P&I Lab, Titech
Human interface Section, P&I Lab, Titech
衝突
実世界の物体は
互いに侵入しない.
跳ね返る.
再現するためには
衝突検知
剛体間に働く力の決定
Human interface Section, P&I Lab, Titech
階層化による高速化
Model Hierarchy:
各ノードが子ノードを含む Bounding Volume を持つ.
 Bouding Volume は,球,直方体などで表される.
階層構造の一番したのノード(葉)が多面体モデルを持
つ.
Bouding Volumeが球で子が2つの場合の例:
Human interface Section, P&I Lab, Titech
階層化による高速化
上の階層から
順に判定
枝刈りできる
→高速
Human interface Section, P&I Lab, Titech
衝突検知
衝突の検出
衝突しているかどうかすべての剛体について調
べる.
剛体の形状は,多面体で表現されている.
衝突検知を簡単にするため,凸形状に分割しておく
Human interface Section, P&I Lab, Titech
凸形状
凸形状の便利な性質
凸形状 距離が極小となる点が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)
Human interface Section, P&I Lab, Titech
衝突検知(GJK)
凸形状A上の点から,
凸形状B上の点へのベ
クトルを
原点を始点に並べると
ベクトルの終点の集合
も凸形状になる
Human interface Section, P&I Lab, Titech
衝突検知(GJK2)
1
V0 : 凸形状内の任意の点
Wi :OViとOWiの内積が最小の点
Vi :三角形Wi-2 Wi-1 Wi内の点で原点に
一番近い点
2
3
4
Human interface Section, P&I Lab, Titech
接触解析
凸形状の便利な性質2
凸形状の交差部分も凸形状
交差部分の形状が簡単に求まる
Half space
representation
D. E. Muller and F.P.Preparata:
“Finding the intersection of two convex” (1978)
Human interface Section, P&I Lab, Titech
接触解析
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.
Dual Transformation
O
n
O
1
n
Human interface Section, P&I Lab, Titech
接触解析
Finding the intersection of two convex(2)
Half space
representation
Dual transform
Vertex of intersection
Dual transform
Human interface Section, P&I Lab, Titech
接触量の積分
The intersection is convex polyhedron
Intersection = upper bound – lower bound
Intersection
=
ー
Contact normal
We can integrate penalty by each face.
d+da
d
d+db
d+db
b
d+da
a
o
Human interface Section, P&I Lab, Titech
力とトルクの計算
 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
The computations for the results are simple
Human interface Section, P&I Lab, Titech
Human interface Section, P&I Lab, Titech
Featherstoneの方法
 リンクiとより葉側のリンク全体
の運動方程式が,
fi  Miri  Zi
のように,書ける.
Mi: リンクiとより葉側のリンクの
質量・慣性による項
Zi: 外力,重力,コリオリ力による項
Mi,Ziは加速度 r によらない
1. 葉のMiは既知なので,
葉から根に向かってMiを
計算する.
fi
2. 根ではfiは0なので,運動方程
式から加速度r を求める.
3. 根から葉に向かって,fi , r
を求める
リンクi
葉のリンク
ri
葉のリンク
リンクi とより葉側のリンク
根のリンク
ダウンロード

発表資料 - 東京工業大学 長谷川晶一研究室