16_1118 超弾性体のVisually Plausibleな実時間シミュレーション：高速化・安定化とアーチファクトの除去について

日本ゴム協会・第28回ゴムの力学研究分科会における講演，2016年11月18日(金)

Transcript

1. 1 超弾性体の 超弾性体のVisually Plausible Visually Plausibleな な 実時間シミュレーション： 実時間シミュレーション： 高速化・安定化と 高速化・安定化と アーチファクトの除去について アーチファクトの除去について 九州大学大学院 九州大学大学院 工学研究院機械工学部門 工学研究院機械工学部門 准教授 准教授 菊 植 亮 菊 植 亮 き く う え りょう き く う え りょう http://rk.mech.kyushu-u.ac.jp/ http://rk.mech.kyushu-u.ac.jp/ http://www.youtube.com/kikuuwe/ http://www.youtube.com/kikuuwe/
3. 4 制御と実時間シミュレーション 制御と実時間シミュレーション  実時間シミュ レーション  ロボット制御  （通常の）シミュ レーション
4. 5 制御と実時間シミュレーション 制御と実時間シミュレーション
5. 6 お話しする内容 お話しする内容 【１】 シミュレーション研究に対する私のスタンス 【２】 超弾性体の実時間シミュレーションについて 【３】 構成則計算の高速化 【４】 時間積分の安定化とArtifactの除去 【５】 「微分包含式」にもとづく摩擦・接触モデル
6. 7 【１】 【１】 シミュレーション研究に対する シミュレーション研究に対する 私のスタンス 私のスタンス 【１】 【１】 シミュレーション研究に対する シミュレーション研究に対する 私のスタンス 私のスタンス  菊植亮：“物理シミュレーションにおける原理的課題：現実・モデル・シミュレータの関係から”，日 本ロボット学会誌，Vol.34，No.6，pp.374-377，2016．
7. 8 動機：自動車産業からのニーズ 動機：自動車産業からのニーズ  【ニーズ】 重量物の精密位置決め 作業をアシストするロボットが欲しい．  【課題】 どのように反力を与えると 位置決め効率がアップするか？  【アイデア】 ロボットで擬似的な摩擦力を 発生すればよいのでは？  現実の特定の材質の摩擦を精密に再現 するのではなく，「摩擦っぽい反力」を 制御で作る.  正確さ・忠実さより，「もっともらしさ」と 「実時間性」が重要 （2003年）
8. 9 もっともらしさ？ もっともらしさ？  現実と合うシミュレーションが正しい．  それはそうだが・・・ 「現実」って何？ 「現実」って何？ 「合う」って何？ 「合う」って何？
9. 10 例１：複数の剛体ブロック 例１：複数の剛体ブロック  収束判定の閾値  各パラメータなど  複数のブロックを バラバラ落とす シミュレーション  シミュレーション結果のブレ ⇒  表面性状  角・稜線の微視的形状  最終的にどんなふうに散らばるのが正解？ ⇒ 実験したら？ ⇒ 実験結果もブレる ⇒ 実験結果もシミュレーション結果もブレる！ 実験結果もシミュレーション結果もブレる！
10. 11 例２：陽解法 例２：陽解法 vs vs 陰解法 陰解法 陽解法 陰解法 解析解  どちらも不正確  しかし陰解法のほうが安定  バネ・質量系 陽解法⇒発散 陰解法⇒収束 不正確だが好都合． 正確さだけが規範ではない？ 不正確だが好都合． 正確さだけが規範ではない？
11. 12 よいシミュレータとは よいシミュレータとは モデル （微分方程式） モデル （微分方程式） 起こりうる 現象の集合  ◦⊂◦であればOK！  ◦はしばしば主観で決まる！  その場合は，「もっともらしさ」（Visual Plausibility）が重要 現実 現実 モデル化 実装 シミュレータ （アルゴリズム） シミュレータ （アルゴリズム） 条 件 シミュレーション 結果の集合 「ブレ」 「ブレ」
12. 13 私見 私見  シミュレータの妥当性評価において，実験結果との 比較はしばしば無意味である．  シミュレーション結果の妥当性は，その現象が必ず 発生するか（あるいは実験で発生したか）ではなく， その現象が発生しえるかで評価されるべきである．  その現象が発生しえるか否かの判定は，しばしば ユーザの主観に依存する． ＝「Visually Plausible」であればＯＫ！！！  シミュレータの不正確さは不可避であるが，その不 正確さは「現実的」側に偏ることが望ましい．
13. 14 【２】 【２】 弾性体の実時間シミュレーション 弾性体の実時間シミュレーション について について 【２】 【２】 弾性体の実時間シミュレーション 弾性体の実時間シミュレーション について について
14. 15 連続体の実時間シミュレーション 連続体の実時間シミュレーション  手術シミュレータ  訓練，術前リハーサルなど  ゲーム・教育用ソフト  一般家庭用ゲーム機  アミューズメント施設等  （子どもたちにも好評）  CG（コンピュータグラフィックス） の作成・編集  3次元形状・物体のデザイン  複雑な形状をパソコン上で 直感的に設計・編集
15. 16 【2009】 経緯 経緯  非線形FEMの新しい数式表現を発見  計算量が比較的少ない  バネ・マスモデルとの関係も明確  極端な変形でも発散しない時間積分法を考案  実時間向け線形方程式ソルバ＋非線形性補正 [Kikuuwe et al., 2009, ACM-TOG/SIGGRAPH] [菊植, 2010, VR学会][菊植ら, 2012, 計算工学講演会] [Kikuuwe, 2016, Visual Computers] 【2007】 【2012】 【2010】
16. 17 節点 位置 節点力 毎秒30回～1000回 時間積分 柔軟体シミュレーションの基本 柔軟体シミュレーションの基本 ユーザ からの入力  節点位置を更新する  陽解法  陰解法  節点に加わる力を求める  バネネットワークモデル  非線形有限要素法，など 構成則
17. 18 「構成則」について 「構成則」について  線形有限要素モデル  ◦ 剛性行列が一定  × 大変形（回転）時に 　　不自然な膨張  材料的に非線形な有限要素モデル  幾何学的に非線形性な有限要素モデル  ◦ 連続体力学と整合性  × 計算量が大  バネネットワークモデル  ◦ 計算量が比較的少  × パラメータの選定方法が不明確
18. 19  陽解法 （前進オイラー法）   計算量が少なく，プログラムも単純  不安定化しやすい．  毎秒>500サイクル  陰解法 （後退オイラー法）   計算量が大きい．連立方程式を解く必要がある．  比較的安定  毎秒<30サイクル 「時間積分」について 「時間積分」について
19. 20 関連研究は2009年時点での情報です 関連研究は2009年時点での情報です 【 【3 3】 】 構成則計算の高速化 構成則計算の高速化 【 【3 3】 】 構成則計算の高速化 構成則計算の高速化  R. Kikuuwe et al.:“An Edge-Based Computationally Efficient Formulation of Saint Venant-Kirchhoff Tetrahedral Finite Elements,” ACM Tran. Graphics, 28(1), pp.8:1-8:13, 2009.
20. 21 構成則：特に 構成則：特にFEM FEMについて について  St. Venant Kirchhoff （StVK），Mooney Rivlin， Neo-Hookean など  弾性エネルギーの定義が異なる  StVKは実時間での利用も：  O’Brien & Hodgins [1999]  Mendoza & Laugier [2003]  Debunne et al. [2001]  Picinbono et al. [2003]  （2009年時点で，CG分野では） ほとんどが陽解法での利用 [Kikuuwe et al., ACM-TOG 2009]  StVKの新しい数式表現で，計算量を削減し， 陰解法での実装を実現
21. 22 StVK (St. Venant Kirchhoff) StVK (St. Venant Kirchhoff) 材料則 材料則  最も単純な「幾何学的に非線形な超弾性構成則」  ひずみエネルギー密度：  ここで E はGreen-Lagrangeひずみ：  ちなみにCauchyひずみは：  一次元のみで考えると，「非線形バネ」に相当  ポテンシャルエネルギーは変位の四次式  力は変位の三次式 変形前 変形後 長さ 力 線形バネ 非線形バネ
22. 23 弾性体（ 弾性体（FE FEモデルとバネモデルの両方） モデルとバネモデルの両方） 節点位置 節点力 （必要） 歪みエネルギー （実装には使わない） 関数 剛性行列 （陰解法で必要） 関数 関数
23. 24 StVK StVK四面体要素 四面体要素 節点位置 節点力 （必要） 歪みエネルギー （実装には使わない） ４次多項式 剛性行列 （陰解法で必要） ２次多項式 ３次多項式
24. 25  　  　  　  　 使われていない 線形FEMの５倍の計算時間 陽解法でよく利用される e.g., Mendoza & Laugier [2003]  節点力（３次多項式）  剛性行列（２次多項式）  弾性エネルギー （４次多項式）  　 StVK StVK四面体要素：文献より 四面体要素：文献より [O'Brien & Hodgins '99] [Picinbono et al. '03] [Picinbono et al. '03] [Capell et al., SIGGRAPH'02] [Capell et al., SIGGRAPH'02]
25. 26 StVK StVK四面体要素の新しい計算法 四面体要素の新しい計算法  非線形バネモデルと「それ以外」の和として表現  従来の式：  各節点の変位から辺，面，四面体毎に節点力を計算 [Picinbono et al. '03]  新しい式：  各辺の伸縮から辺および辺のペア毎に節点力を計算 [Kikuuwe et al., 2009]
26. 27 新しい計算法の直感的な解釈 新しい計算法の直感的な解釈  四面体ごとにではなく，「辺」と「辺のペア(TSEP)」ご とに接点力（f）と剛性行列（K）を計算． 従来手法 新手法  旧手法：四面体ごとの展開された多項式の計算  新手法：メッシュ全体で因数分解された多項式の計算
27. 28 四面体の歪みエネルギー 四面体の歪みエネルギー T sym T sym sym 変形前 変形後 長さの二乗の 変化量 初期長さ  「　 」の２次式として表現  辺の長さの双二次形式  Green-Lagrange歪み  四面体の歪みエネルギー
28. 29 メッシュ全体の歪みエネルギー メッシュ全体の歪みエネルギー sym T T 対角成分 非対角成分
29. 30 対角成分：非線形バネのネットワーク 対角成分：非線形バネのネットワーク  非線形バネ  力は変位の三次式  ポテンシャルエネルギー は変位の四次式 長さ 力 線形バネ 非線形バネ T W
30. 31 それ以外：辺の対（ それ以外：辺の対（TSEP TSEP）の効果 ）の効果  TSEP = Tetrahedron-Sharing Edge Pair （四面体を共有する辺の対）  非線形FEMを非線形SNモデルと「そ れ以外」の和として表現． T W
31. 32 弾性エネルギーを偏微分すると・・・ 弾性エネルギーを偏微分すると・・・  節点力 （3次式）　  接線剛性行列 （2次式）　 辺 辺の ペア
32. 33 アルゴリズム（非線形 アルゴリズム（非線形FE FEモデル） モデル） 辺 辺 辺の ペア 辺の ペア 辺 節点位置 剛性行列 節点力  バネモデルに簡略化も可能
33. 34 バネモデルに簡略化も可能：結果 バネモデルに簡略化も可能：結果  膨張：FEM　⇔ 縮小：バネのみ  無視できる違いか否か？ → 用途による？  バネモデルのパラメータ設定に利用可能
34. 35 辺の対 辺（バネ） 競合する論文： 競合する論文： [Delingette, ISBMS'08] [Delingette, ISBMS'08]  2008年7月出版．我々の論文投稿の直後 伸縮の剛性 "tensile stiffness" 角度の剛性 "angular stiffness" 体積の剛性 "volumetric stiffness"
35. 36 FLOP FLOP数の比較 数の比較 -15%～+29% -85%～-73% -73%～-62% 非線形FEM
36. 37 計算時間の比較 計算時間の比較 (a) (b) (c) (d)  接点力と剛性係数の計算時間は3割～5割減にとどまる． メモリ配置の最適化が必要か？ -38% -48% -53% -28%
37. 38 計算時間の比較（時間積分含む） 計算時間の比較（時間積分含む） (a) (b) (c) (d)  サンプリングタイム固定，剛性計算時間短縮 → 陰解法の反復演算の 回数を増やすことが可能 → おそらく過渡応答がリアルに
38. 39 実装例（ 実装例（2008 2008年当時） 年当時） (a) 直方体 四面体数1,200 (b) 恐竜 四面体数2,458 (c) スタンフォード・バニー： 四面体数4,064 (d) スタンフォード・アルマジロ： 四面体数7,554  デスクトップパソコン(Dell社製Precision 390)  Intel Core 2 Duo E6700, 2.66 GHz  Microsoft Windows XP Professional  Microsoft Visual C++ 2005 Academic Edition  ＋力覚提示装置  SensAble社製Phantom Omni  いまはもっとなめらかに動きます！
39. 40 論文： 論文： ACM Trans. Graphics, 2009 ACM Trans. Graphics, 2009
40. 41 そこそこ引用されています そこそこ引用されています
41. 42 業界の反応 業界の反応(1) (1)  菊植の手法をGPUに実装して，約５０倍高速化 できたらしい．（ただし陽解法．剛性行列無し）
42. 43 業界の反応 業界の反応(2) (2)  単眼カメラ画像からの弾性物体の形状再現に利用  Haochine et al., 2015, IEEE Trans. Visualization and Computer Graphics
43. 45  既存の非線形FEMの新しい表現  昔からある式を眺めていたら，もっと簡単に表せるという ことに気づいたという話  モデル（構成則）自体は新しくない！！  学術的な利点：  双二次形式に基づいた新しい表現  バネ・ネットワークモデルとの明確な関係  実用的な利点：  従来手法よりも計算が速い  ただし，陰解法ではこの効果はそれほど大きくない．  陰解法で実時間で実装したこと自体が当時としては珍しい？  バネ・ネットワークモデルを用いるときに，バネ定数の決 定のための参照モデルとして利用可能 「構成則の高速化」のまとめ 「構成則の高速化」のまとめ
44. 46 【４】 【４】 時間積分の安定化と 時間積分の安定化と アーチファクトの除去 アーチファクトの除去 【４】 【４】 時間積分の安定化と 時間積分の安定化と アーチファクトの除去 アーチファクトの除去  Ryo Kikuuwe:“A Time-Integration Method for Stable Simulation of Extremely Deformable Hyperelastic Objects,” The Visual Computer, doi:10.1007/s00371-016-1225-0, 2016．  菊植，岡，山本，"極端な変形下でも安定な実時間非線形有限要素シミュレーション"，計算 工学講演会論文集，2012.  菊植:“陰解法ベース実時間シミュレーションのための線型方程式ソルバの選定"，第15回日 本バーチャルリアリティ学会大会論文集，2010.
45. 47 毎秒30回～1000回 構成則 時間積分 節点 位置 節点力 構成則と時間積分 構成則と時間積分 ユーザ からの入力  節点位置を更新する  陽解法  陰解法  節点に加わる力を求める  バネネットワークモデル  非線形有限要素法，など
46. 48  離散化・テーラー展開  未知数 は節点速度  行列 は剛性行列＋α  連立方程式ソルバ  1/30秒ほど毎に収束計算 時間積分：特に陰解法について 時間積分：特に陰解法について  陽解法よりは安定だが，しばしば「爆発」 時間積分 構成則
47. 49 もっとも単純な連立方程式ソルバ： もっとも単純な連立方程式ソルバ： 共役勾配法 共役勾配法  比較的高速  CG分野では好まれる  Baraff&Witkin [1998]  M ller&Gross [2004] ü  Kikuuwe et al. [2009]  の が正定対称で ないと「破綻」  非線形FEMにおいて，剛性行列 は正定とは限らない A
48. 50 剛性行列が正定な構成則を「作る」 剛性行列が正定な構成則を「作る」  Rotational Linear Methods  変位勾配テンソルを伸縮成分と 回転成分に分解（直交分解，QR 分解，特異値分解・・・）  線形FEMの剛性行列を四面体ご とに回転させる．  要素の裏返りを防げない．裏返っ たときの対応が一意ではない．  一度裏返ると復元しない． [Irving et al. 2006]  非線形FEM＋非正定ソルバの組み 合わせ（＝積分法を工夫する）という アプローチは，前例がなかった．
49. 51 非正定ソルバ 非正定ソルバ  の が正定対称行列でない場合． BiCG Lanczos, 1952 BiCG Lanczos, 1952 CGS Sonneveld, 1989 CGS Sonneveld, 1989 BiCGStab van der Vorst, 1992 BiCGStab van der Vorst, 1992  共役勾配法をやめて BiCGStabを実装すると ⇒比較的安定になったが， それでも発散 A [Zhou & Walker, 1994]
50. 52 数値計算分野 数値計算分野 vs vs 実時間シミュレーション 実時間シミュレーション  実時間シミュレーション  30～50ms毎に「それなり」の解が必要  収束の成否に関わらず計算打ち切り  収束の速さより，残差履歴の単調性が 重要 繰り返し回数 残差の対数 [Zhou & Walker, 1994]  数値計算分野  残差（　　　　　　 ）を速く収束 させることがポイント． [ANSYS Japanのウェブサイトより]
51. 53 解決法：最小残差スムージング（ 解決法：最小残差スムージング（MRS MRS） ） 初期値 [Zhou&Walker, 1994] 繰り返し計算 解 修正解  既存のソルバと組み合わせて，残差を単調減少させる．  CG／VR／実時間シミュレーション分野で応用例がない．  数値計算分野でもあまり使われていない？ 既存のソルバ （BiCGStabなど） 残差スムージング
52. 54 残差スムージング（ 残差スムージング（MRS MRS）の効果 ）の効果  BiCGStab＋MRS  BiCGStabのみ [菊植, VR学会, 2010]  発散はほぼ完全に防げた．
53. 55 さらに改良 さらに改良  BiCGStab+MRSのみ  ときおりパルス状の 異常な変形を示す．  BiCGStab+MRS +非線形性補正  異常変形がなくなった [菊植ら, 計算工学講演会,2012]
54. 56 QMR QMR（ （Quasi-Minimum Residual Quasi-Minimum Residual）法 ）法 MRS (Minimum Res. Sm.) Zhou & Walker,1994 MRS (Minimum Res. Sm.) Zhou & Walker,1994 QMRS (Quasi-Min. Res. Sm.) Zhou & Walker,1994 QMRS (Quasi-Min. Res. Sm.) Zhou & Walker,1994  BiCG + QMRS = QMR  QMRを使えばMRSは不要らしい BiCG Lanczos 1952 BiCG Lanczos 1952 CGS Sonneveld 1989 CGS Sonneveld 1989 第一世代の 非正定ソルバ 第二世代の 非正定ソルバ 残差スムーザ BiCGStab van der Vorst,1992 BiCGStab van der Vorst,1992 QMR Freund & Nachtigal 1991,1994 QMR Freund & Nachtigal 1991,1994 [Kikuuwe 2016, Visual Computers] [Zhou & Walker, 1994]
55. 57 経緯のまとめ 経緯のまとめ BiCGStab BiCGStab 非線形性 補償 非線形性 補償 QMR QMR VR学会, 2010 Visual Computers, 2016 （ジャーナル論文） CG CG BiCGStab BiCGStab 非線形性 補償 非線形性 補償 BiCGStab BiCGStab MRS MRS 高速化 FEM 高速化 FEM ACM-TOG, 2009 MRS MRS 計算工学講演会, 2012 挙動はそのままで 簡略化 かなり良くなった 安定になった 不安定 あまり改善せず 構成則 時間積分 しかし，たまに変な挙動
56. 58  要素数：2,968  頂点数：975  ステップ時間：20ms  Intel Core i7-X980 （6コア），3.33Ghz  OpenMPで並列化 （スレッド数5） 実装 実装 [菊植ら, 計算工学講演会,2012]
57. 59 比較 比較  QMR+非線形補正 vs BiCGStab  QMR+非線形補正 vs QMR  StVK-FEM vs Rotational Linear [Kikuuwe 2016, Visual Computers]
58. 60  実時間シミュレーション向け連立方程式ソルバ  忘れられていた（？）技術の再利用  時間積分の際の非線形性補正  極端な変形でもＯＫ 「時間積分の安定化」のまとめ 「時間積分の安定化」のまとめ
59. 61 【４】 【４】 おまけ：より現実に近い挙動を目指した おまけ：より現実に近い挙動を目指した シミュレーション研究 シミュレーション研究 【４】 【４】 おまけ：より現実に近い挙動を目指した おまけ：より現実に近い挙動を目指した シミュレーション研究 シミュレーション研究  X. Xiong et al.:“A Differential-Algebraic Multistate Friction Model,” Proc. SIMPAR, Springer, pp.77-88, 2012.  X. Xiong et al.:“A Differential-Algebraic Contact Model with Nonlinear Compliance,” Proc. DSCC, Paper 8642, 2012  R. Kikuuwe et al.:“Admittance and Impedance Representations of Friction Based on Implicit Euler Integration,” IEEE Trans. Rob., 22(6), pp.1176-1188, 2006.
60. 62 ご紹介した研究成果 ご紹介した研究成果 モデル化 実装 条 件 「ブレ」 「ブレ」  ざっくりした現象をざっくり現実的にシミュレーション モデル （微分方程式） モデル （微分方程式） 現実 現実 シミュレータ （アルゴリズム） シミュレータ （アルゴリズム） 起こりうる 現象の集合 シミュレーション 結果の集合
61. 63 もっと精密なシミュレーション もっと精密なシミュレーション 「ブレ」 「ブレ」  もっと工学的な研究もしています！！ 起こりうる 現象の集合 シミュレーション 結果の集合 モデル化 実装 条 件 モデル （微分方程式） モデル （微分方程式） 現実 現実 シミュレータ （アルゴリズム） シミュレータ （アルゴリズム）
62. 64 摩擦と接触：力と運動の関係が不連続 摩擦と接触：力と運動の関係が不連続  摩擦（クーロン摩擦）  接触（片側剛体拘束） 速度 摩擦 力 位置 接触力 速度 摩擦 力 接触 力 位置  符号関数：  ダイオード関数： （微分包含式； differential inclusions）  未開拓領域：「微分包含式」で表現されて， 実験データともよく合う摩擦・接触モデルの構築
63. 65 「微分包含式」ベースの力学モデル 「微分包含式」ベースの力学モデル  柔軟体の接触力モデル  滑り出す前のヒステリシスなど を再現する摩擦力モデル  転がり面接触摩擦モデル ⇒ タイヤの摩擦力の再現  Dr. Xiaogang Xiong  2014年，当研究室で 博士修了  現在，ハルビン工科 大学（深セン） 助教授 [Xiong, Kikuuwe & Yamamoto, 2014, ASME-JAM] [Xiong, Kikuuwe & Yamamoto, 2013, Tribology Letters]  現在論文執筆中 [菊植，三上，山本，SI2007] [フー・菊植 2016（発表予定）， 　JSME交通・物流部門大会] [東，菊植，SI2015]
64. 66 接触力モデル（柔軟体） 接触力モデル（柔軟体）  摩擦と接触の不連続性を考慮に入れた上で，過去 の実験データと整合する微分包含式を発見 Hunt-Crossleyモデル 提案モデル [Xiong et al., 2014, ASME-JAM] 過去論文の実験データ [Yamamoto et al., ICRA2013]
65. 67 摩擦力モデル 摩擦力モデル  過去の実験データと整合する波形を描く微分包含 式を発見 [Xiong, Kikuuwe & Yamamoto, 2013, Tribology Letters]
66. 68 タイヤの摩擦の力学モデル タイヤの摩擦の力学モデル  不連続性を含む偏微分方程式：「偏微分代数包含式」 微小部位の バネ・ダンパ力 微小部位の 摩擦力 微小部位の 対地速度 タイヤにかかる総力 微小部分の 変形速度 [菊植，三上，山本，SI2007]  地面とタイヤの間の相互作用・・とても複雑！  接触面内で混在する静止摩擦状態と動摩擦状態
67. 70 タイヤの摩擦の力学モデル：実装 タイヤの摩擦の力学モデル：実装  いくつかの非線形項を組み合わせて，経験式（Magic Formula） の出力値と近づくように改良中  スエギリ時のタイヤの挙動も再現  「停車してハンドルを切った状態から，徐行発進すると，ハンドルが ゆっくりスルスルと戻っていくような状況」を 物理演算に基づいてシミュレーション可能（になるはず） [フー・菊植，JSME交通・物流部門大会, 　 2016（発表予定）]
68. 71 全体のまとめ 全体のまとめ 全体のまとめ 全体のまとめ
69. 72 [1] シミュレーション研究に対する私のスタンス [2] 超弾性体の実時間シミュレーションについて [3] 構成則計算の高速化 [4] 時間積分の安定化とArtifactの除去 [5] おまけ：摩擦・接触力のモデル化について  ご清聴ありがとうございました． まとめ まとめ