3. RUI を用いた情報世界へのインタラクションシステム
3.6. 力覚提示装置としての RUI システム
3.6.2. エアホッケー物理モデル
3.6.2.3. ペナルティ法
ペナルティ法に関して詳しく説明する．
パックが壁に向かって運動を行い，壁に衝突する場合を考える．先に説明した通り，シ ミュレーションにおいてパックが壁に衝突した場合，パックは壁にめり込んでおり，パッ クの端と壁表面の位置関係は下図のようになる．
図 3.13 ペナルティ法
侵入量 l
パック端 壁表面
反力 F
3.6.2.4. 剛体の反発運動
並進・回転運動を行っている剛体であるパックが，壁やマレットに衝突した場合には反 発運動を行う．まずパックが壁に衝突した時の反発運動に関して説明する．
パックは壁への衝突時に，壁の垂直方向にはペナルティ法による反力を受け，平行方向 には摩擦力を受けて運動を行う．今回は摩擦力としては，剛体の速度に比例した力，つま り粘性抵抗を与えることで摩擦を表現する．
図 3.14 剛体の反発運動
図 3.14のようにパックが２次元平面上で並進・回転運動を行って，壁に対して衝突・反 発を行う場合を考える．パックが壁に衝突（侵入）した時刻をtとし，このときのパック の位置を(Xp(t),Yp(t))，速度を(Vx(t),Vy(t))，回転変位をθ(t)，角速度をω(t)，パックの質 量をMp，半径をr，壁表面の位置をYwとおく．また，ペナルティ法におけるバネ係数を
k，ダンパ係数をc，壁との摩擦としての粘性係数をuとおく．
このときの反力Fy(t)と摩擦力Fx(t)は次のように与えられる．
) ( }
) ) ( {(
)
( t k Yp t r Yw c Vy t
Fy = − − − − ×
…(3.7))}
( )
( { )
( t u Vx t r t
Fx = − × − ω
…(3.8)この式(3.7)，(3.8)と前出の式(3.1)～(3.6)を用いると，tからdt秒後の時刻t+dtでのパッ クの位置，速度，回転の変位，角速度は次のように与えられる．
dt t Vx t Xp dt t
Xp ( + ) = ( ) + ( )
…(3.9)Yp ( t + dt ) = Yp ( t ) + Vy ( t ) dt
…(3.10)Y 軸
衝 突 前 の パック
壁 パック
摩擦力
衝 突 後 の パック
X 軸
反力
Mp dt t t Fx Vx dt t
Vx ( )
) ( )
( + = +
…(3.11)dt Mp
t t Fy Vy dt t
Vy ( )
) ( )
( + = +
…(3.12)ω(t)dt
θ(t) dt)
θ(t + = +
…(3.13)2 ) 2 ( 1 ) ) (
( )
( dt I mr
I t t rFx dt
t + = ω − =
ω
…(3.14)この式(3.9)～(3.14)を用いて，パックが壁に衝突した場合の反発運動の計算を行う．
次に，パックとマレットが衝突した時の反発運動に関して以下に説明する．
パック（円）がマレット（円）に衝突すると，パックはマレットとの中心を結んだ中心 線方向に反力を受け，また，中心線と垂直の方向に摩擦力を受けて運動を行う．反力の計 算にはペナルティ法を用いて計算を行い，摩擦力は，パックと壁の反発運動時と同様に粘 性抵抗による力とする．
図 3.15 パックとマレットの反発
図 3.15のように，パックがマレットに衝突した場合を考える．パックがマレットに衝突 した時刻をtとし，このときのパックの中心位置を(Xp(t),Yp(t))，速度を(Vpx(t),Ypy(t))， 角速度をω(t)，マレットの中心位置を(Xm(t),Ym(t))，速度を(Vmx(t),Ymy(t))とおく．また，
摩擦力
X 軸 Y 軸
侵入量
θ(t)
θ(t)
マレット パック
パ ッ ク の 運動方向
マレットの 運動方向 反力
φ (t)
パックとマレット間で生じる反力Fcと摩擦力Ftは次のように求めていく．
まず，角度θ(t)，パックとマレットの中心位置の関係を求めると次のようになる．
))2 ( ) ( 2 ( )) ( ) ( (
) ( ) )) (
( cos(
t Ym t Yp t
Xm t Xp
t Xm t t Xp
− +
−
= −
θ …(3.15)
))2 ( ) ( 2 ( )) ( ) ( (
) ( ) )) (
( sin(
t Ym t Yp t
Xm t Xp
t Ym t t Yp
− +
−
= −
θ …(3.16)
次に，パックの中心線方向，つまり反力方向の速度成分Vpc(t)，中心線と垂直方向，つ まり摩擦力方向の速度成分Vpt(t)と，マレットの反力方向の速度成分Vmc(t)，摩擦力方向 の速度成分Vmt(t)は次のようになる．なお摩擦力方向は，パックの時計回りの方向（図で あればX軸負，Y軸正方向）を正とする．
)) ( 2 / cos(
)) ( cos(
)
(t Vpx t Vpy t
Vpc = × θ + × π −θ …(3.17) ))
( 2 / sin(
)) ( sin(
)
(t Vpx t Vpy t
Vpt =− × θ + × π −θ …(3.18) ))
( 2 / cos(
)) ( cos(
)
(t Vmx t Vmy t
Vmc = × θ + × π −θ …(3.19) ))
( 2 / sin(
)) ( sin(
)
(t Vmx t Vmy t
Vmt =− × θ + × π −θ …(3.20)
式(3.15)～(3.20)を用いると，反力Fcと摩擦力Ftは次のように求められる．
)) ( ) ( ( )}
( )}
( ) ( { )}
( ) ( { { )
(t k Xp t Xm t 2 Ypt Ymt 2 Rp Rm c Vpc t Vmct
Fc =− × − + − − + − × − …(3.21)
)) ( )
( )
( ( )
( t u Vpt t Rp t Vmt t
Ft = × + × ω −
…(3.22)これより，パックがマレットから受けるX軸方向の力Fpx(t)とY軸方向の力Fpy(t)は次 のようになる．
)) ( 2 / cos(
) ( )) ( cos(
)
( t Fc t Ft t t
Fpx = × θ + × π + θ
…(3.23))) ( 2 / sin(
) ( )) ( sin(
)
( t Fc t Ft t t
Fpy = × θ + × π + θ
…(3.24)式(3.15)，(3.16)，(3.22)，(3.23)，(3.24)を用いて，時刻tからdt秒後の時刻t+dtでのパッ クの位置，速度，回転変位，角速度は以下のように求められる．
dt t Vpx t
Xp dt t
Xp( + )= ( )+ ( ) …(3.25) dt
t Vpy t Yp dt t
Yp( + )= ( )+ ( ) …(3.26) Mp dt
t t Fpx Vpx dt
t
Vpx ( )
) ( )
( + = + …(3.27) Mp dt
t t Fpy Vpy dt
t
Vpy ( )
) ( )
( + = + …(3.28) dt
t t dt
t ) ( ) ( )
(
θ ω
θ
+ = + …(3.29) 2 ) ( 1 ) ) (( )
( dt Ip Mpr2
Ip t Ft t Rp
dt
t+ =
ω
+ × =ω
…(3.30)パックがマレットと衝突している間に関しては，式(3.25)～(3.30)を用いてパックの運動 計算を行う．
今回のエアホッケーの物理モデルでは，RUIを通して力覚提示を行うため，パックがマ レットに衝突した際に，マレットがパックから受ける力を求める必要がある．マレットが パックから力を受ける方向は，マレットがCGアバタの両腕の先端にあり，マレットの中 心が円を描いて運動していることを考慮し，その円の接線方向とする．
図 3.16 マレットが受ける力の方向
マレットの中心が描く円の中心を
( a , b )
，半径をr
とおくと，マレットの中心位置が)) ( ), (
( Xm t Ym t
であるときの，マレットの中心を通る円の接線の方程式は次のように与え られる．b b t Ym a r b x t Ym
a t
y Xm +
+ −
− −
− −
= ( ) ( )
) (
)
(
2・・・(3.31)
これより，この接線がX軸となす角度α(t)については次のような関係式が得られる．
)2 ) ( 2 ( ) ) ( (
) )) (
( cos(
b t Ym a t Xm
b t t Ym
− +
−
= −
α …(3.32),
)2 ) ( 2 ( ) ) ( (
) )) (
( sin(
b t Ym a t Xm
a t t Xm
− +
−
− −
=
α …(3.33)
式(3.23)，(3.24)，(3.32)，(3.33)より，マレットがパックから受ける力Fm(t)は次式のよう になる．なお，マレットの中心が描く円の右回転方向を正とする．
反力
Y軸
X 軸 θ
マレット θ
パック
摩擦力
α マレットが 受ける力
マレットの
中心軌道
3.6.2.5. エアホッケーアプリケーションの仕様
エアホッケーアプリケーションにおいて各定数は次のように設定した．
・運動間隔 0.001 [s]
・エアホッケーの台の仕様 横幅 1200 [mm] 縦幅 2000 [mm]
・パックの仕様に関する定数 質量 40 [g] 半径 60 [mm]
・マレットの仕様に関する定数 質量 120 [g] 半径 100 [mm]
・パックが壁に衝突した際の運動に関する定数
使用するペナルティ法における定数 バネ係数 1000 [N/m] ダンパ係数 10[ N/m⋅s−1] 摩擦における粘性係数 0.0001[ N/m⋅s−1]
・パックがマレットに衝突した際の運動に関する定数
使用するペナルティ法における定数 バネ係数 1000 [N/m] ダンパ係数 10[ N/m⋅s−1] 摩擦における粘性係数 0.0001[ N/m⋅s−1]
操作者の体感として，上記にあるようなマレットやパック，エアホッケーの台を使用し て，エアホッケーを行っているという感覚を得ることができる．
図 3.17 エアホッケーアプリケーション
3.6.3. 運動処理周期の短縮化
物体同士が衝突した際の反力計算にペナルティ法を使用したが，ペナルティ法を使用す る上で，プログラム内で運動計算を行う周期を短くする必要がある．何故なら，運動を行 う時間間隔が長いと物体の侵入量が大きくなり過ぎ，物体に加える反力も大きくなり過ぎ てしまうために反発の表現が不安定となるからである．よって，運動処理周期を短縮する 必要がある．
運動処理周期の短縮化の目安としては，安定的に物体同士の接触を表現できる最低限の 目安としての値である1[kHz]（運動の時間間隔が1[ms]）とする[37]．また，運動処理周期 のサンプリングレートが1[kHz]よりも高くなるほど，より硬い物体の力覚を提示すること が可能となる．これを受けた力覚提示装置に関する既存の研究として，より高い更新レー トを実現させることを目的として赤羽ら[38]は，PCと力覚提示装置の中間層にハプティッ クコントローラを配置することで 10[kHz]という高速制御を実現し，より硬い物体の力覚 提示に成功した．これを考慮し，1[kHz]より一桁上回る10[kHz]への処理周期の短縮化を目 標とし，最低限でも1[kHz]を目指す．
短縮化の手法としては2種類提案する．なお，一般的な環境下での実装を目指し，リア
ルタイムOSは用いずWindows系OS上で使える手法を考案，実装し，実験を行った．
実験には，OSは全てWindows2000で，メモリ，ビデオカードやCPU数の異なる複数の PCを使用した．PC名はCPUが二つの場合Dual1，CPUが一つの場合はSingle1というよ うにおいた．
表 3.2 短縮化手法１で使用したPC
PC名 CPU メモリ ビデオカード
Single1 Celeron 1.4GHz 382.48 MB 無し
36.00 MB
Single2 Celeron 1.4GHz 383.48 MB GeForce4 Ti4800
124.38 MB
Dual1 P-Ⅲ 800MHz×2 255.48 MB GeForce2 MX 100/200
59.12 MB
Dual2 Xeon 3.06GHz×2 2046.96 MB ELSA GLADIAC FX 736Ultra
123.96 MB
・短縮化手法1
短縮化手法1として，計算と描画の処理を分担・並列に行えるように，マルチスレッド
り，FullScreenサイズかWindowサイズかで描画方式が異なるため，２パターンの描画形式 で計測を行った．また，プログラムに関しては，DualCPUのPCでは各CPUで描画と計算 を分担して実行した．
図 3.18 短縮化手法1における各処理の流れ
各 PC における計算処理に掛かった時間の計測結果において，平均処理時間，平均サン プリングレートを表 3.3に示す．使用PC名において(F)はFullScreenサイズ，(W)はWindow サイズで描画を行ったことを示す．また計測したデータにおいて，1 ループにかかった時 間が設定した時間の倍の 200[μs]以上となった場合をエラーとみなし，そのエラーの総数 を総データ数で割った結果をエラー率とした．また，力覚提示に最低限必要なサンプリン グレートを 1[kHz]として，ループに 1[ms]以上掛かった場合は別に算出した．その結果を 図 3.19に示す．取得した総データ数は31997である．
表 3.3 短縮化手法1の結果（計算処理時間）
PC名 平均処理時間[ms] 標準偏差 平均サンプリングレート[kHz]
Dual2(F) 0.1008 0.8634
× 10
−2 9.921Dual2(W) 0.1009 0.8647
× 10
−2 9.914Dual1(F) 0.1017 0.2817
× 10
−1 9.830Dual1(W) 0.1017 0.1343
× 10
−1 9.830Single2(F) 3.637 14.85 0.2749
Single2(W) 3.824 14.88 0.2615
Single1(F) 0.1122 0.1305 8.910
Single1(W) 3.631 15.42 0.2754
図 3.19 短縮化手法1の結果（エラー率）
結果より，この短縮化手法1では，SingleCPUのPCではなく，計算と描画処理を二つの CPUで分担したDualCPUのPCを使用した場合に，ほぼ10[kHz]のサンプリングレートで 計算処理が行えることが確認できた．SingleCPUにおいても，平均的には9[kHz]のサンプ リングレートで処理が行える場合もあるが，処理時間にばらつきが見られる．このばらつ きによって，物体の運動が発散してしまう恐れが残る．
・短縮化手法2
二つ目の手法は手法1よりも，より一般的な環境下での実装を考えて，SingleCPUのPC での運動処理周期の短縮化を可能とするものを目指した．手法としては，プログラム内に おいて描画を更新する際に，前回の描画更新からの経過時間を計算し，その経過時間に相 当する回数だけまとめて運動処理を行うものである．この手法2により，物理モデルの運 動に関しては1[kHz]の短縮化を可能とした．
実際に，短縮化手法 1 で運動処理周期のレートが落ちた「SingleCPU の PC における
Windowサイズでのプログラム」に関して，表 3.4に示すPCを用いて実際の運動計算処理
に掛かった時間を計測した．
表 3.4 短縮化手法2で使用したPC
PC名 OS CPU メモリ ビデオカード
DesktopPC1 Windous2000 Celeron 1.4GHz 383.48MB GeForce FX 5600
125.50MB
DesktopPC2 Windous2000 Celeron 1.4GHz 383.48MB GeForce4 Ti 4800
124.38MB DesktopPC3 Windous2000 Pentium 4 2.4GHz 511.48MB GeForce4 Ti 4200
124.38MB
NotePC WindowsXp Pentium 4-M 2.5GHz 511.34MB GeForceFX Go5600
61.50MB
表 3.5 短縮化手法2の結果
PC名 利用データ数 平均処理時間[ms] 平均サンプリングレート[kHz]
DesktopPC1 40000 1.000 0.9999
DesktopPC2 40000 1.000 1.000
DesktopPC3 40000 1.000 1.000
NotePC 40000 1.000 0.996
計測の結果を表 3.5に示す．性能の異なるSingleCPUのPCにおいて，運動処理に掛か る平均時間はほぼ 1[ms]となり，平均してほぼ 1[kHz]のサンプリングレートで運動処理が 行われることを確認できた．
力覚提示は描画の更新レートに依存しており，本システムでは一般的なレートである 60[Hz]とした．ここで，力覚提示は力覚演算の連続性とその提示タイミングの２つの要素 に分けられる．力覚演算の連続性は先の赤羽らの論文に報告されているように対象物体の 硬さの正確な提示に必要であり，1[kHz]以上にて制御を行うことが必要である．しかしな がら，その提示タイミングは視覚や聴覚等他の感覚統合に関わる問題であり，視覚提示の 更新レートによる60分の1秒程度のずれがあっても違和感は少なくなっている．
また，短縮化手法2とペナルティ法を使用した物理モデル内で，物体同士の衝突・反発 運動に関しての精度を測定する実験を行った．
実験内容は，一次元運動のみを行う二つの円盤を用意し，静止している一つの円盤Aに