JP2007011460A - 物体の変位シミュレーション方法及び物体の変位シミュレーション装置、並びに力覚提示装置 - Google Patents

Abstract

【課題】空間で拘束条件を欠く物体に対して掴み等の操作を行う場合に、物体の変形を、有限要素法を利用して算出するに当たり、短時間で行うことができる物体の変位シミュレーション方法を提供する。
【解決手段】
物体Ｂに発生する変形を有限要素モデルの全体の剛性方程式に基づいて解析し、シミュレーションを行う。この際、全体の剛性方程式における有限要素モデルの全ノードに関する変位を、自由ノードの未知変位と、指先ＵＳに接触する接触面から１つ抽出された拘束面に関する既知変位と、拘束面以外の残りの接触面に関し、等価荷重で変換可能な既知変位とに区分し、有限要素モデルの全体の剛性方程式から算出された剛性行列の逆行列を、ウッドバリーの公式が適用できるように予め変形する。未知変位をウッドバリーの公式に基づいて算出する。
【選択図】 図６

Description

本発明は、ＶＲ（バーチャルリアリティ）技術に応用できる物体の変位シミュレーション方法及び物体の変位シミュレーション装置並びに力覚提示装置に関するものである。

近年、ＶＲ（バーチャルリアリティ）技術の発達により、ディスプレイに視覚情報が提示された仮想空間に力覚を提示することが可能である。又、物体の弾性モデルを導入し、柔らかい物体（以下、仮想柔軟物体という）の力覚表現も可能となり、仮想環境のリアリティは向上している。従来から、仮想柔軟物体の力覚表現方法には、大別すると、バネマスモデルによる方法と、有限要素法による手法がある。バネマスモデルによる方法は仮想柔軟物体の変位、反力を少ない計算量で求めることができる（非特許文献１）。しかし、バネマスモデルによる方法は、仮想柔軟物体のヤング率、ポアソン比等の特性値をモデルに直接組み込むことができないため、仮想柔軟物体の変形は、現実と異なるものとなる。

一方、有限要素法による手法は、仮想柔軟物体の特性値をモデルに組み込むことが可能であるため、仮想柔軟物体の変形を正確に表現できる。しかしながら、有限要素法による手法は、後述するように、仮想柔軟物体の変形を求める計算量が膨大であるため、力覚提示処理のボトルネックとなっている。

しかし、最近は、有限要素法による仮想柔軟物体の力覚表示に関して、モデルの縮退化や、変位の「等価荷重」を計算する手法（非特許文献２）を用いることにより、短時間で変形の計算処理ができる技術が提案されている。非特許文献２が提案している、変位の「等価荷重」による手法は、初期の仮想柔軟物体の拘束条件を変化させないことを前提に、予め剛性行列の逆行列を求め、その逆行列を用いて、仮想指先と、仮想柔軟物体との接触箇所の変位を「等価荷重」に変換して、変形を算出するものである。前記非特許文献２の手法は、少ない計算量で、変形が、算出できることが特徴であるが、剛性行列の逆行列を保証するために、常に仮想柔軟物体は拘束条件を持つ、すなわち、仮想柔軟物体の一部が空間に拘束されている必要がある。なお、この非特許文献２の手法については、後述する。 非特許文献３は、この手法を利用して、仮想箸による仮想柔軟物体の変形を可能としている。

なお、非特許文献４は、後述するウッドバリー（Woodbury）の公式（後述する式（１６））が記載されている文献である。この非特許文献４は、後では、参考文献という。
宮崎愼也，安田孝美，横井茂樹，鳥脇純一郎，「仮想弾性物体の対話操作のためのモデル化と実現」，電子情報通信学会論文誌，１９９６年１１月，Vol,J79-A，No.11，pp.1919-1926 広田光一，金子豊久，「柔らかい仮想物体の力覚表現」，情報処理学会論文誌，１９９８年１２月，Vol.39，No.12，pp.3261-3268 北村喜文，銅子研一，来山真，岸野文郎，「仮想箸のための物体の変形と力覚表示の一検討」，日本バーチャルリアリティ学会論文誌，２００２年９月，Vol.7，No.3，pp.291-299 Gene H.Goulb, Charles F.Van Loan,"Matrix Computations THIRD EDITION",The Johns Hopkins University Press, 1996.11.1,pp.50-51

非特許文献２が提案している手法では、図１（ａ）に示すように、仮想指先１０で仮想柔軟物体１２を突くような形態であれば、仮想柔軟物体１２が斜線部１２ａで示す箇所で空間に拘束されていても問題はない。しかし、同じ条件で図１（ｂ）に示すように２本の仮想指先１０，１１で仮想柔軟物体１２を掴むような場合、斜線部１２ａが常に空間に拘束されるため、仮想柔軟物体１２は、移動することができず、正しい変形及び反力が得られない。

図１（ｃ）に示すように、仮想柔軟物体１２の仮想指先１０，１１に対する接触箇所１２ｂに拘束条件を設けることができれば、仮想柔軟物体１２の変形及び反力を正しく表現できる。しかし、仮想指先１０，１１の動きに伴い、仮想柔軟物体１２に対する接触箇所１２ｂは変動するため、拘束条件は常に変化することになる。拘束条件が変化することで、剛性行列は変化し、その逆行列を短時間で求めることは難しくなる問題がある。すなわち、剛性行列をＲ^ｎ×ｎ（なお、ｎ×ｎはｎ行ｎ列を表す）とした場合、逆行列演算算出に、通常Ｏ（ｎ^３）の計算量を要する。このことについては後述する。

ここで、従来の、物体変形演算方法について説明する。
まず、前提として、ノード数ｎで表現される線形な有限要素モデルを想定する。有限要素法では、対象物体（有限要素モデル）全体の変形と荷重の関係を全体の剛性行列により表現する。すなわち、この想定した有限要素モデルによる剛性方程式は、荷重Ｆ（∈Ｒ^ｎ）と変位ｕ（∈Ｒ^ｎ）と剛性行列Ｋ（∈Ｒ^ｎ×ｎ）を用いて式（１）で与えられる。

変位ｕを与えられたノード、すなわち、拘束条件をもつノードにおいては、変位ｕは既知であり、荷重Ｆは未知である。そこで、既知の変数と未知の変数とをまとめ、式（１）の剛性方程式を書き換えると、

になる。ここで、ｕ_ｋは既知変位であり，Ｆ_ｋは既知荷重である。又、ｕ_ｕは拘束されていないの自由ノードの未知変位とし，Ｆ_ｕは未知荷重とする。そして、式（２）を変形し、

とする。ここで、Ｋ_０，Ｋ_１はＫの部分行列とする。Ｅは単位行列である。式（３）において、Ｋ_０，Ｋ_１，Ｅを含む剛性行列をＡとする。この解は、式（４）のように与えられる。

ただし、拘束が与えられていない場合は、逆行列Ａ^−１が求まらないことから解は求まらない。又、この従来の変形演算方法では、拘束箇所が変化する場合、その逆行列Ａ^−１を短時間で求めることは、剛性行列をＲ^ｎ×ｎとした場合、逆行列演算算出に、通常Ｏ（ｎ^３）の計算量を要する。

非特許文献２の手法では仮想指先１０，１１と仮想柔軟物体１２との接触箇所の変位を等価な荷重（以下、この荷重のことを等価荷重Ｆ^＊という）に変換することにより、容易に物体変形を求めている。ここで、非特許文献２の手法を説明すると、式（４）に相当する式は式（５）で書き表すことができる。

ここで、Ｌは逆行列Ａ^−１の部分正方行列、ｕ_ｋ ^＊は等価荷重Ｆ^＊が働くノードの変位である。このノードは、拘束されていないが、仮想指先１０，１１と接触していることから、その変位は、既知となる。仮に、拘束ノードの変位をｕ_ｋ＝０とおいた場合、等価荷重Ｆ^＊は、式（６）となる。

この式（６）を式（５）に代入することにより、容易に拘束されていない自由ノードの未知変位ｕ_ｕを求めることができる。すなわち、物体の変形算出を行うことができる。しかし、非特許文献２の手法では、図１（ｃ）に示すように掴む操作に対しては、対応できない問題がある。すなわち、掴む動作においては、仮想柔軟物体１２は空間に拘束されておらず、仮想指先１０，１１に拘束されている。仮想指先１０，１１は仮想柔軟物体１２表面上を移動するため、行列Ａの要素は変化する。このため式（５）の逆行列Ａ^−１の要素は変化し、式（５）の逆行列Ａ^−１を短時間で求めることは難しく、この結果、変形を有限要素法を利用して算出するに当たり、短時間で行い難い問題がある。なお、本明細書で、空間で拘束条件を欠く物体とは、空間に拘束されていない状態（すなわち、非拘束状態）であって、空間に位置する物体が、前記仮想指先等の拘束物以外のものには拘束されていない状態をいう。

本発明の目的は、空間で拘束条件を欠く物体に対して掴み等の操作を行う場合に、物体の変形を、有限要素法を利用して算出するに当たり、従来よりも短時間で算出を行うことができる物体の変位シミュレーション方法及び変位シミュレーション装置並びに力覚提示装置を提供することにある。

上記目的を達成するために、請求項１は、空間で拘束条件を欠く物体を拘束する拘束物が前記物体に対して移動することによって発生する前記物体の変形を、前記物体に関する有限要素モデルについての全体の剛性方程式に基づいて演算するシミュレーションを行う変位シミュレーション方法において、前記全体の剛性方程式における前記有限要素モデルの全ノードに関する変位を、自由ノードの未知変位と、前記拘束物に接触する接触面から、一つ抽出された拘束面に関する既知変位と、前記拘束面以外の残りの接触面に関し、等価荷重で変換可能な既知変位とを算出し、前記有限要素モデルの全体の剛性方程式から算出される剛性行列の逆行列を、予め変形することにより算出した変形式に対して、ウッドバリーの公式を適用することにより、前記全ノードに関する変位のうち既知変位以外の変位である未知変位を算出することを特徴とするシミュレーションを行う変位シミュレーション方法を要旨とするものである。

請求項２の発明は、請求項１において、前記ウッドバリーの公式には、対角行列を含み、該対角行列は、接触面に関与している対角要素が１，それ以外が０で構成することを特徴とする。

請求項３の発明は、請求項１において、前記ウッドバリーの公式には、対角行列を含み、該対角行列は、対角要素を全て１とすることを特徴とする。
請求項４の発明は、請求項１乃至請求項３のうちいずれか１項において、前記未知変位を算出した後、前記拘束物への反力を算出することを特徴とする。

請求項５の発明は、空間で拘束条件を欠く物体を拘束する拘束物が前記物体に対して移動することによって発生する前記物体の変形を、前記物体に関する有限要素モデルについての全体の剛性方程式に基づいて演算するシミュレーションを行う変位シミュレーション装置において、入力された前記拘束物の位置情報と、前記有限要素モデルを構成するノードの位置情報に基づいて、前記全体の剛性方程式における前記有限要素モデルの全ノードに関する変位のうち、前記拘束物に接触する接触面から一つ抽出された拘束面に関する既知変位、及び、前記拘束面以外の残りの接触面に関し、等価荷重で変換可能な既知変位の両既知変位を算出する既知変位算出手段と、前記有限要素モデルの全体の剛性方程式から算出される剛性行列の逆行列を予め変形することにより算出した変形式に対して、ウッドバリーの公式適用することにより、前記全ノードに関する変位のうち既知変位以外の変位である未知変位を算出する未知変位算出手段を備えていることを特徴とするシミュレーションを行う変位シミュレーション装置を要旨とするものである。

請求項６の発明は、請求項５において、前記算出された未知変位に基づいて前記拘束物への反力を算出する反力算出手段を備えたことを特徴とする。
請求項７の発明は、空間で拘束条件を欠く物体を拘束する拘束物が前記物体に対して移動することによって発生する前記物体の変形を、前記物体に関する有限要素モデルについての全体の剛性方程式に基づいて演算するシミュレーションを行い、力覚を提示する力覚提示装置において、入力された前記拘束物の位置情報と、前記有限要素モデルを構成するノードの位置情報に基づいて、前記全体の剛性方程式における前記有限要素モデルの全ノードに関する変位のうち、前記拘束物に接触する接触面から一つ抽出された拘束面に関する既知変位、及び、前記拘束面以外の残りの接触面に関し、等価荷重で変換可能な既知変位の両既知変位を算出する既知変位算出手段と、前記有限要素モデルの全体の剛性方程式から算出される剛性行列の逆行列を予め変形することにより算出した変形式に対して、ウッドバリーの公式適用することにより、前記全ノードに関する変位のうち既知変位以外の変位である未知変位を算出する未知変位算出手段と、前記算出された未知変位に基づいて前記拘束物への反力を算出する反力算出手段と、前記算出された反力に基づいて力を発生する力生成手段を備えたことを特徴とする力覚提示装置を要旨とするものである。

（作用）
ここで、本発明で行われる、物体の変形を短時間で求めるための手法について説明する。

好適には、短時間で物体の変形を求めるために下記（１）〜（３）の方法を行うことにより、計算量をＯ（ｎ）に低減することができる。
（１） 物体の拘束箇所を複数の接触面から１つ抽出する。

（２） 全ての接触面を構成するノード全部に着目し、演算する。
（３） ウッドバリー（Woodbury）の公式（式（１６））を適用する。
なお、以下では、拘束面は、拘束物で拘束する拘束箇所のことをいう。図６においては、指先ＵＳが物体Ｂを掴み操作することによって拘束する拘束物として示されている。又、接触面は、拘束物と物体Ｂの接触点が存在する３角形パッチ面である（図６参照）。なお、接触面を構成する３つのノードには、同じ変位が加わるものとする。又、物体Ｂは、柔軟性を備えたものとしている。

（拘束面の抽出）
図６に示すように、拘束物（図６では指先ＵＳ）で物体Ｂを掴む場合、複数の接触面が存在する。そこで、拘束される面（すなわち、拘束面）を常に複数の接触面の中から１つ選択する。そして、拘束面には既知変位ｕ_ｋを、残りの接触面には等価荷重Ｆ^＊を与える。現時点では等価荷重Ｆ^＊の値は未知であるが、後述する方法で求める。

（接触面構成ノードの着目）
行列Ｗ（∈Ｒ^ｎ＋ｎ）は対角行列で、接触面に関与している対角要素は１，それ以外は０で構成されている。式（７）に示すように行列Ｗを導入しても、その解は変化しない（請求項２）。

この結果、式（５）は、以下の通りに書き直すことができる。

すなわち、接触面以外のノードは変形に関与しないことから、接触面のノードが関与するノードだけを求めれば変形は求めることができる。

（ウッドバリー（Woodbury）の公式の適用）
拘束物（図６では指先ＵＳ）を物体Ｂの拘束面に対して動かしたとき、拘束物の動きに伴い拘束面が変化するため、剛性行列Ａの要素も変化する。いま、剛性行列Ａを下記の通りに表す。

Ａ＝Ａ_０＋ΔＡと表す。
ここで、Ａ_０は、「初期処理」でそれ自身の逆行列Ａ_０ ^−１を保証するため、１つの拘束面を与えたときの行列である。又、ΔＡは、変化分である。剛性行列Ａの要素は、拘束が解かれるノードと新たに拘束されるノードが、以下のように変化する。

拘束面の変化を式（１０）のように表す。

Ｕ（∈Ｒ^ｎ＋ｐ）はｐ個の単位ベクトルで構成される行列である。又、Ｖ（∈Ｒ^ｎ＋ｐ）はΔＡの非ゼロベクトルで構成される行列である。なお、Ｔは転置行列を意味する。ｐは「初期処理」で拘束した面（拘束面）と、接触面から抽出した拘束面の２つの面を構成しているノードのベクトル次数である。たとえば、後述するフローチャート（図４参照）で、「初期処理」とは、前記フローチャートが実行される前に行われる処理であり、「接触面から抽出した拘束面」とは、前記フローチャートを実行する際に「接触面から抽出した拘束面」を指す。具体的に３次元空間を考えると、図２に示すように１つの面は３角形パッチで構成されているため、３つのノードＰ１，Ｐ２，Ｐ３を含み、且つ、１つのノードは３次元の位置座標をもつことにより、１つのノード自由度は３であるから、１面当たり９のベクトル次数を持つ。従って、ｐとして取り得る最大数はｐ＝１８である。

式（９）の変化は式（１０）を用いることにより、以下のとおり展開できる。

式（１１）より、式（８）は下記の通りとなる。

ウッドバリー（Woodbury）の公式（後述する式（１６））を用いて、（Ａ_０＋ＵＶ^Ｔ）^−１を以下の通りに算出する。

ここで、Ｂ＝Ａ_０ ^−１Ｗ，Ｃ＝Ａ_０ ^−１Ｕであり、これらは、乗算計算をしなくても要素の代入操作で求めることができる。すなわち、これらは、行列計算せずに、要素の配列の代入により、求めることができる。Ｂは零列ベクトルを含む疎行列ベクトルであり、非零ベクトルの列リストを保持することにより、行列乗算の計算を低減することができる。

なお、対角行列Ｗを無視して、すなわち、実質的に対角行列Ｗの対角要素が全て「１」とした場合であっても、変形を求めるに当たって従来よりも短時間で処理することは可能である（請求項３）。

ここで、式（１３）の乗算計算量を、行列の乗算順に示すと、下記の通りとなる。
・ Ｔ_１＝Ｖ^ＴＣ
１８×３ｎ×１８
・ Ｔ_２＝Ｖ^ＴＢ
１８×３ｎ×９（１＋m）
・Ｔ_３＝（Ｅ＋Ｔ_１）^−１Ｔ_２
１８×１８×９（１＋m）
・Ｔ_４＝ＣＴ_３
３ｎ×１８×９（１＋m）
ここで、３ｎは全ノードの自由度、ｍは接触面数、（１＋ｍ）は「初期処理」で仮に拘束条件を与えた面１つと接触面数ｍを指す。又、「１８」は、先に述べたpの最大値である。又、ｎ≫ｍより、逆行列（Ｅ＋Ｔ_１）^−１は計算量として考慮していない。

これより、式（１３）に係る計算量はＯ（ｎ）である。
なお、拘束面以外の接触面が変化する場合は、その接触面を構成する逆行列のノードだけを更新すればよく、計算量は上記計算量よりも小さくなる。

（等価荷重Ｆ^＊の計算）
拘束されていない自由ノードの未知変位ｕ_ｕを求めるためには、等価荷重Ｆ^＊を求める必要がある。式（１２）のｕ_ｋ ^＊とｕ_ｋが既知であることから、（Ａ_０＋ＵＶ^Ｔ）^−１Ｗの部分正方行列の逆行列を用いることにより、等価荷重Ｆ^＊を算出できる。すなわち、

であるから、以下の通りになる。

ここで、α，βは（Ａ_０＋ＵＶ^Ｔ）^−１Ｗの部分正方行列である。これより、未知変位ｕ_ｕを求めることができる。

（ウッドバリー（Woodbury）の公式）

なお、この欄のウッドバリー（Woodbury）の公式で使用している行列Ａ，Ｕ，Ｖは、既に説明しているＡ，Ｕ，Ｖの定義とは異なり、この欄のみに適用されるものと理解されたい。すなわち、この欄では、Ａは、ｎ×ｎ行列であり、Ｕ，Ｖはそれぞれｎ×ｐ行列（ただし、ｐ＜ｎ，通常はｐ≪ｎ）である。なお、参考文献では、式（１６）中、（Ｅ＋Ｖ^ＴＡ^−１Ｕ）^−１の部分は、（Ｉ＋Ｖ^ＴＡ^−１Ｕ）^−１となっているが、この参考文献での「Ｉ」は単位行列（＝Ｅ）を意味しているため、本明細書では、説明を整合するために上記式（１６）では「Ｉ」の代わりに「Ｅ」を使用している。

なお、接触面における、拘束物に対する反力の計算は、式（２）で算出される（請求項４）。
請求項５では、既知変位算出手段は、入力された前記拘束物の位置情報と、前記有限要素モデルを構成するノードの位置情報に基づき前記拘束物に接触する接触面から、１つ抽出された拘束面に関する既知変位と、前記拘束面以外の残りの接触面に関し、等価荷重で変換可能な既知変位を算出する。又、未知変位算出手段は、前記有限要素モデルの全体の剛性方程式から算出された剛性行列の逆行列が予めウッドバリーの公式に関連付けされ、前記剛性方程式から算出された変形式を使用して、前記未知変位を、前記ウッドバリーの公式に基づいて算出する。

さらに、請求項６では、前記算出された未知変位に基づいて、反力算出手段は前記拘束物への反力を算出する。
請求項７の力覚提示装置においては、請求項６に記載の反力算出手段が算出した拘束物への反力に基づいて、力生成手段は、力を発生する。このことにより、力覚を提示する。

請求項１乃至請求項６の発明によれば、空間で拘束条件を欠く物体に対して掴み等の操作を行う場合に、物体の変形を、有限要素法を利用して算出するに当たり、従来よりも短時間で物体の変形の算出を行うことができる。

又、請求項２の発明では、ウッドバリーの公式に対角行列を含み、該対角行列の接触面に関与している対角要素は１，それ以外は０で構成することにより、乗算計算量を少なくでき、物体の変形を、有限要素法を利用して算出するに当たり、従来よりも短時間で物体の変形の算出を行うことができる。

請求項３の発明では、ウッドバリーの公式に対角行列を含み、該対角行列の対角要素を全て１とすることにより、乗算計算量を少なくでき、物体の変形を、有限要素法を利用して算出するに当たり、従来よりも短時間で物体の変形の算出を行うことができる。

請求項４の発明では、短時間で未知変位が算出されるため、拘束物への反力の算出するまでの時間も短縮することができる。
請求項５の発明では、請求項１の方法を直接使用して、物体の変形を、有限要素法を利用して算出するに当たり、短時間で行うことができ、シミュレーションを短時間で行うことができる。

請求項６の発明では、短時間で未知変位が算出されるため、拘束物への反力の算出するまでの時間も短縮することができる変位シミュレーション装置を提供することができる。
請求項７の発明では、短時間で未知変位が算出され、拘束物への反力の算出するまでの時間も短縮することができるため、力覚提示するまでの時間も短縮することができる。

（第１実施形態）
以下、本発明を具体化した第１実施形態を図３〜１０を参照して説明する。
本実施形態の物体の変位シミュレーション装置（以下、単にシミュレーション装置１００という）は、図５に示すように所定の物体Ｂを定義し、この物体Ｂを掴む際に、該物体Ｂの変位等のシミュレーションが可能とされている。シミュレーション装置１００は、力覚提示装置を兼用しており、一対の力覚提示ユニット２０Ａ，２０Ｂと、一対のパーソナルコンピュータ（ＰＣ３０Ａ，３０Ｂ）と、ＣＲＴや液晶表示装置等のディスプレイ４０とから構成されている。

一対の力覚提示ユニット２０Ａ，２０Ｂは、いずれもロボットアーム２２、位置検出部２４、力生成部２６を備えている。各力覚提示ユニット２０Ａ，２０Ｂの位置検出部２４、力生成部２６は、それぞれＰＣ３０Ａ，３０Ｂに電気的に接続されている。力生成部２６は力生成手段に相当する。各力覚提示ユニット２０Ａ，２０Ｂは、オペレータに対するマン・マシン・インターフェイスとして機能し、オペレータから与えられる操作量を入力するとともに、オペレータに対して力を返す働きをするようにされている。本実施形態では、１つの力覚提示ユニット２０Ａ，２０Ｂは、それぞれ１本の指としての機能を果たす。本実施形態では、物体Ｂに対する掴み操作により、物体Ｂの変位等をシミュレートするシミュレーション装置として構成されており、２組の力覚提示ユニット２０Ａ，２０Ｂが用意されている。

ロボットアーム２２上の所定の単数（又は複数）の作用点は、ディスプレイ４０の画面上では、物体Ｂに対して、接触可能な操作点として設定されており、オペレータにより印加された力（荷重）に基づいて、該作用点は、三次元の任意方向に移動可能とされている。

各位置検出部２４は、ロボットアーム２２上の各作用点の三次元座標系（例えば、ｘｙｚの直交座標）上での位置を検出可能である。力覚提示ユニット２０Ａの位置検出部２４は、ＰＣ３０Ａに対して、検出したロボットアーム２２の作用点の位置を位置情報として出力する。又、ＰＣ３０Ａは、ＰＣ３０Ｂに対して相互通信可能にＬＡＮ等の通信線を介して接続されており、互いの各種情報の通信が可能である。例えば、ＰＣ３０Ｂには、ＰＣ３０Ａから、力覚提示ユニット２０Ａの位置検出部２４から入力した位置情報等が入力可能にされている。又、力覚提示ユニット２０Ｂの位置検出部２４は、ＰＣ３０Ｂに対して、検出したロボットアーム２２の作用点の位置を位置情報として出力する。

ディスプレイ４０上では、ＰＣ３０Ｂは、入力された力覚提示ユニット２０Ａ，２０Ｂからの各作用点の位置情報に基づいて、力覚提示ユニット２０Ａ，２０Ｂの各ロボットアーム２２に対応した仮想指先（以下、単に指先ＵＳという）が表示される。すなわち、該指先ＵＳの表面上の所定の点が前記操作点に対応する。前記指先ＵＳは、拘束物に相当する。

ＰＣ３０Ｂは、図示しない外部記憶装置を備えており、該外部記憶装置には、物体Ｂを定義するための基本情報及び、シミュレーションプログラムが格納されている。そして、ＰＣ３０Ｂは、該基本情報に基づいて、ディスプレイ４０上に、物体Ｂを表示制御可能である。この基本情報は、例えば、物体Ｂの形状情報（立方体、直方体、円錐等の立体形状を表す三次元情報）や、ポアソン比や、ヤング率等の情報等を含む。本実施形態の物体Bは、前記ポアソン比や、ヤング率により、柔軟性を有する。

一方、ＰＣ３０Ｂは、前記作用点の位置と、物体Ｂの位置に基いて、各指先ＵＳと物体Ｂとの接触状態を判定し、この接触状態に応じて拘束物としての各指先ＵＳに生じるべき反力（物体Ｂ側からオペレータ側へと返される力）の算出が可能である。このように、ＰＣ３０Ｂは各指先ＵＳの操作点に生じるべき反力を算出すると、この反力に基いて各操作点に加えるべき制御力を示す力制御ベクトルが算出される。そして、これら力制御ベクトルを示すデータは、力覚提示ユニット２０Ａの力生成部２６に対しては、ＰＣ３０Ａを介して付与し、力覚提示ユニット２０Ｂの力生成部２６に対しては直接付与される。

力覚提示ユニット２０Ａ，２０Ｂの力生成部２６は、付与された力制御ベクトルに基いてロボットアーム２２の作用点に加える力を発生させることが可能である。
（作用）
さて、上記のように構成されたシミュレーション装置１００の作用を、図４を参照して説明する。図４は、ＰＣ３０Ｂが前記シミュレーションプログラムに従って実行する物体Ｂの変位シミュレーションに関するフローチャートであり、一旦開始した後は、終了の割り込み指令がない限り、継続的にループ処理が行われる。なお、このフローチャートを実行する前の「初期処理」において、ヤング率や、ポアソン比や、物体Ｂのノード情報（三次元座標系の位置情報）に基づいて、全体の剛性行列Ｋは求められており、ＰＣ３０Ｂの図示しない外部記憶装置に格納されているものとする。又、ノードには、物体Ｂの重心に位置するノード、或いは、該重心にもっとも近いノードが重心ノードとされている。又、このフローチャートが実行される前の「初期処理」において、１つの拘束面を与えたときの行列としてＡ_０が、初期値として算出されているものとし、それ自身の逆行列Ａ_０ ^−１の算出が保証されている。そして、オペレータは、このフローチャートを実行するときにおいて、ロボットアーム２２を操作する。

（Ｓ１０）
ＰＣ３０Ｂは、各位置検出部２４による各ロボットアーム２２の作用点の位置情報に基づいて、指先位置の更新を行う。

（Ｓ２０）
ＰＣ３０ＢはＳ１０で更新した指先位置と物体Ｂとの衝突判定を指先ＵＳ毎に行う。すなわち、衝突判定では、指先ＵＳの位置情報と、変形していない状態の物体Ｂの各ノードの位置情報とで、衝突判定が行われる。ここで、変位していない状態の物体Ｂとは、該物体Ｂの移動を反映させた、変形していない物体Ｂをいう。

Ｓ２０で、衝突なしと判定した場合、このフローチャートを一旦終了する。Ｓ２０で、衝突ありと判定された場合、Ｓ３０に移行する。
（Ｓ３０）
次に、ＰＣ３０Ｂは指先ＵＳが接触した複数の接触面の中から、拘束面を１つ選択することにより、既知変位ｕ_ｋを求めるとともに、等価荷重Ｆ^＊が働くノードの変位ｕ_ｋ ^＊を求める。この場合の拘束面の選択は、任意でよいが、本シミュレーションプログラムでは、一方の指先ＵＳの所定の部位に対し接触した接触面を拘束面としている。このようにして、ＰＣ３０Ｂは拘束面及び接触面における既知変位（変位量）ｕ_ｋ ^＊とｕ_ｋを求めることにより、指先ＵＳと接触する物体Ｂの接触面、及び指先ＵＳの侵入量を求める。ここで、侵入量は、前述の変位ｕ_ｋ ^＊とｕ_ｋに相当する。

なお、変位ｕ_ｋ ^＊は前述したように等価荷重Ｆ^＊が働くノードの変位である。このノードは、拘束されていないが、指先ＵＳと接触していることから、その変位は、指先ＵＳの位置情報から算出する。

このように、ＰＣ３０Ｂは既知変位算出手段に相当する。
（Ｓ４０）
前回の処理サイクル時における接触面と、今回処理サイクルの接触面が変化しているか否かの判定を行う。すなわち、少なくとも１つの接触三角形パッチに関して接触の有無が変化している場合には、接触面が変化しているものと判定されて、Ｓ５０に移行する。

（Ｓ５０）
このステップでは、ＰＣ３０Ｂにより、「逆行列計算」を行う。ここでいう、「逆行列計算」とは、（Ａ_０＋ＵＶ^Ｔ）^−１Ｗを求める式（１３）の演算である。ここで、式（８）と、式（１２）とを見れば分かるように、剛性行列の逆行列Ａ^−１は、（Ａ_０＋ＵＶ^Ｔ）^−１に変形されている。そして、（Ａ_０＋ＵＶ^Ｔ）^−１は、ウッドバリーの公式である式（１６）の左辺と同じである。従って、（Ａ_０＋ＵＶ^Ｔ）^−１Ｗは、ウッドバリーの公式に基づいて算出される。

そして、この「逆行列計算」は、Ｓ４０で前記変化したと判定された接触面が、拘束面であった場合と、拘束面でなかった場合とで異なっている。すなわち、Ｓ４０で前記変化ありと判定された接触面が、拘束面であった場合、式（１３）の右辺を、最初から計算し直すことにより行われる。

一方、Ｓ４０で前記変化ありと判定された接触面が、拘束面でない場合は、前回の処理サイクル時で求めた（Ａ_０＋ＵＶ^Ｔ）^−１Ｗの一部を求め直すことにより行われる。このように、最初から演算し、或いは、一部を求め直して得られた（以後、これを「（Ａ_０＋ＵＶ^Ｔ）^−１Ｗの更新」という）、（Ａ_０＋ＵＶ^Ｔ）^−１Ｗは、前記外部記憶装置により記憶保持される。

なお、このステップにおいては、式（１３）の右辺の計算を行うに当たり、Ｂ＝Ａ_０ ^−１Ｗ，及びＣ＝Ａ_０ ^−１Ｕは、シミュレーションプログラム上において、要素の配列の代入で求める。Ｓ５０の計算の終了後、Ｓ６０に移行する。

（Ｓ６０）
このステップでは、「移動／変形／反力の算出」が、ＰＣ３０Ｂにより行われる。
「変形の算出」は、すなわち、拘束されていない自由ノードの未知変位ｕ_ｕの算出のことである。まず、ＰＣ３０Ｂは、式（１５）を使用して等価荷重Ｆ^＊を算出する。そして、ＰＣ３０Ｂは求めた等価荷重Ｆ^＊に基づいて、剛性方程式から変形された変形式である式（１２）を使用して未知変位ｕ_ｕの算出を行う。ここで、式（１２）は、前記剛性方程式から算出される剛性行列の逆行列が、ウッドバリーの公式が適用できるように予め変形されており、前記剛性方程式から算出されたものである。

なお、Ｓ４０において、「変化なし」と判定して、Ｓ６０に移行している場合、式（１２）の使用時に、ＰＣ３０Ｂは前回の処理サイクル時に記憶保持されている（Ａ_０＋ＵＶ^Ｔ）^−１Ｗを使用する。又、Ｓ５０からＳ６０に移行している場合、Ｓ５０において、ＰＣ３０Ｂは最初から演算し、或いは、「（Ａ_０＋ＵＶ^Ｔ）^−１Ｗの更新」された（Ａ_０＋ＵＶ^Ｔ）^−１Ｗを式（１２）の使用時に使用する。

又、接触面における「反力」の算出は、既知変位ｕ_ｋと、前記算出された未知変位ｕ_ｕと式（２）を使用して行われる。又、「移動の算出」は、重心ノードの変位を算出することである。重心ノードの変位は、物体Ｂの移動量に相当する。

このステップで算出された自由ノード（重心ノードを含む）の未知変位ｕ_ｕと、既知変位ｕ_ｋ及び変位ｕ_ｋ ^＊に基づいて、ＰＣ３０Ｂは、ディスプレイ４０の画面に対して画像制御を行って、変形状態及び移動状態の少なくともいずれか一方の状態となっている物体Ｂを表示する。

併せて、ＰＣ３０Ｂは、算出した反力に基づいて各操作点に加えるべき制御力を示す力制御ベクトルを算出し、これら力制御ベクトルを示すデータを、力覚提示ユニット２０Ａの力生成部２６に対しては、ＰＣ３０Ａを介して付与し、力覚提示ユニット２０Ｂの力生成部２６に対しては直接付与する。

この結果、力覚提示ユニット２０Ａ，２０Ｂの力生成部２６は、付与された力制御ベクトルに基いてロボットアーム２２の作用点に加える力を発生させ、力覚を操作者に提示する。このように、ＰＣ３０Ｂは、未知変位算出手段及び反力算出手段に相当する。Ｓ６０の処理が終了後、このフローチャートを一旦終了する。

（試験例）
さて、以上のように構成された、シミュレーション装置１００の具体的な試験例について説明する。

（行列サイズによる計算量の比較）
ＰＣ３０Ｂとして、Windows2000（マイクロソフト社）を搭載したパーソナルコンピュータ（Xeon 1400MHz（インテル社），ＲＡＭ；５１２ＭＢ）を使用して、図１（ｃ）に示すように物体Ｂを２点で掴み操作する際、行列サイズに応じて、計算量がどの程度変化するかを試験した。その結果を図７に示す。

実施例１，２では、シミュレーション装置１００の力覚提示ユニット２０Ａ，２０Ｂは、PHANToM 1.5/6DOFとPHANToM DESKTOP（いずれもSensAble Technologies社製）からなる。実施例１は、前記（１）〜（３）の方法を全て行った前記実施形態に沿って行った場合である。実施例２は、前記（１）、（２）の方法を使用するとともに、かつ（３）のウッドバリー（Woodbury）の公式（式（１６））を使用するが、式（８）中、行列Ｗを無視、すなわち、行列Ｗの対角要素が全て「１」として計算したものである。

比較例１のシミュレーション装置のハード構成は、実施例１と同じとし、従来手法である式（４）により、算出した場合の例である。
物体Ｂは、６種類のノード数の異なるモデルを使用し、接触面数ｍ＝２とした。なお、モデルには縮体化を施した。その結果は、図７に示す。図７において、横軸は、行列サイズ（ｎ×ｎ）を示し、縦軸は時間（sec）である。同図に示すように、実施例１，及び実施例２は、行列サイズに応じて、処理時間は増えるが、比較例１と比較して、いずれの行列サイズにおいても短時間で処理できることが示されている。特に、実施例１では、いずれの行列サイズにおいても行列サイズに比例して、短時間で計算が可能であることが分かる。

（反力の検証）
次に、実施例１の構成を利用して、掴み操作の際の反力の検証を行った。ここで、反力の検証の対象として、図５に示すように一辺７ノードの六面体で構成された立方体を物体Ｂとした。この物体Ｂの行列サイズは、モデルの縮退化を施して、６５４行６５４列、ヤング率：50000［Pa］、ボアソン比：0.4とし、立方体の長さ：１［ｍ］とした。なお、予め拘束するＡ_０の逆行列Ａ_０ ^−１は初期処理で求めておいた。図５では左右から、物体Ｂを指先ＵＳにより挟んだときの変形を示している。なお、図５では、指先ＵＳは、いわゆる指先形状で示しておらず、シンボルとして表示されている。

この状態で侵入量｜ｕ_１｜を与えたときの反力ｆ１（同図中、右側の面に働く反力），ｆ２（同図中、左側の面に働く反力）を空間に拘束箇所を持つ場合（比較例２では物体Ｂの底面において、拘束されている。：比較例２）と、実施例１の２ケースで求めた。なお、初期処理において、逆行列Ａ_０ ^−１を求めるために、立方体の底の一部を空間に拘束した。

図８（ａ）は比較例２の場合において、ｘ方向の侵入量を示している。図８（ｂ）は、比較例２において、反力ｆ１，ｆ２を示している。図８（ａ），（ｂ）に示すように、反力ｆ１，ｆ２は侵入量に応じてその大きさは変化するが、｜ｆ１｜≠｜ｆ２｜であることが分かる。

図９（ａ）は実施例１の場合において、ｘ方向の侵入量を示している。図９（ｂ）は、実施例１において、反力ｆ１，ｆ２を示している。図９（ａ）に示すように、｜ｕ_１｜＝｜ｕ_２｜となった。又、図９（ｂ）に示すように、反力ｆ１，ｆ２は侵入量に応じてその大きさは変化し、｜ｆ１｜＝｜ｆ２｜であることが分かる。このように、実世界では、物体Ｂを掴み、安定であるときは、物体Ｂに働く合力は０であることから、実施例１は、現実の変形と比較し、差異がないことが分かる。

（拘束箇所の移動に伴う検証）
図１０（ａ）〜（ｄ）は、実施例１のシミュレーション装置において、物体Ｂを掴んだ状態で、指先ＵＳの拘束位置を上から下へ移動させたときの０．５［sec］毎の物体Ｂの変形状態を示している。なお、図１０では、指先ＵＳは、いわゆる指先形状で示しておらず、シンボルとして表示されている。

同図に示すように、指先ＵＳによる拘束位置の移動に伴い物体Ｂの変形も変化している。この際、拘束面が変化したときに係る（Ａ_０＋ＵＶ^Ｔ）^−１Ｗの算出時間は、０．０４［sec］、拘束面以外の接触面が変化したときに係る「（Ａ_０＋ＵＶ^Ｔ）^−１Ｗの更新」の計算時間は、０．０１［sec］であり、短時間で行われていた。なお、接触面が変化しないときは、計算は不要である。

又、オペレータが感じる反力に違和感がないことがら、実施例１は、本発明の効果を達することができることが分かる。
以上の実施例の効果を下記にまとめる。

・ 図７に示すように、パーソナルコンピュータ（Xeon 1400MHz（インテル社），ＲＡＭ；５１２ＭＢ）で演算した場合、比較例１（従来技術）では、行列サイズが３００行３００列の場合、１秒を超えており、リアルタイムに力覚を提示できない。しかし、実施例２では、行列サイズが３００行３００列の場合、０．１３秒であり、リアルタイムに力覚を提示できた。

・ 又、実施例１では、行列サイズが１２００行１２００列の場合であっても、図７に示すように、０．１１秒であり、リアルタイムに力覚を提示できた。
・ このように、実施例２は、行列サイズが、６００行６００列程度まで、実施例１は、さらに、大きな行列サイズまで、リアルタイム（実時間）処理が可能であり、状来よりも、拘束条件を欠く物体に対して把持動作する場合の変形算出及び力覚提示を短時間で行うことが可能となる。

なお、この先、電子計算機が発達することにより、演算速度がさらに向上すると思われるが、本実施例で例示したように、本発明を利用すれば、演算処理の向上に伴って、ノード数、すなわち、行列サイズを大幅に増加することが可能となる。

（第２実施形態）
第２実施形態は、変位シミュレーション装置及び力覚提示装置として、多指触覚インターフェイスに具体化したものである。この多指触覚インターフェイスの構成を、図１１〜図１３を参照して、説明する。図１１は多指触覚インターフェイスの概略を示す斜視図、図１２は本実施形態の触覚指の斜視図である。

多指触覚インターフェイス１１０は、ベース１１１、アーム部１１２、ハンド１１３を備えている。アーム部１１２は、第１アーム１１４及び第２アーム１１５からなる。第１アーム１１４は、ベース１１１に対して、３自由度を有する第１腕関節１１６を介して連結されている。第１腕関節１１６は図１１に示すように第１腕関節１１６の中心において互いに直交する軸ｘ０, y０,z０の周りでそれぞれ第１アーム１１４を回転する。この軸ｘ０,y０,z０は基準座標系の軸となる。又、第２アーム１１５は第１アーム１１４に対して１自由度を有する第２腕関節１１７（肘部）を介して連結されている。第２腕関節１１７は、第２アーム１１５を図１１に示すように第２腕関節１１７の中心線である軸ｍの周りで能動回転する。

ハンド１１３は触覚指ベース１２０と、触覚指ベース１２０に設けられた５本の触覚指１２１〜１２５を備えている。５本の触覚指１２１〜１２５は、人間の拇指〜小指の５本の指と相対可能に、列状に配置されている。すなわち、第１触覚指１２１が人間の拇指、第２触覚指１２２が人差し指、第３触覚指１２３が中指、第４触覚指１２４が薬指、第５触覚指１２５が小指と向かい合うように触覚指ベース１２０に配置されている。触覚指ベース１２０の基端側に設けられた取付部は、第２アーム１１５の先端部に対して手首関節１１８（図１３参照）を介して取付けされている。手首関節１１８は図１１に示すように手首関節１１８の中心において互いに直交する軸ｘｍ,yｍ,ｚｍの周りでそれぞれ触覚指ベース１２０を回転する。

この軸ｘｍ,yｍ,ｚｍは触覚指ベース座標系の軸となる。第１腕関節１１６、第２腕関節１１７及び手首関節１１８は、駆動源としての駆動モータをその自由度の数に一致する個数分備えており、各駆動モータを制御することにより、前述した軸の周りで回転する。

又、第１腕関節１１６、及び第２腕関節１１７及び手首関節１１８において、前述した各軸の周りで回転駆動する駆動モータには回転角度を検出する回転角度検出手段としてのロータリエンコーダＡＲＥ１〜ＡＲＥ７が設けられている(図１３参照)。これらのロータリエンコーダＡＲＥ１〜ＡＲＥ７により、前記各軸の周りの回転角度が検出される。なお、図１３では、説明の便宜上、ロータリエンコーダＡＲＥ１〜ＡＲＥ７をまとめて図示している。前記ロータリエンコーダＡＲＥ１〜ＡＲＥ７にて検出された回転角度は、それぞれアーム関節角度ｑｉ（ｉ＝１，…７）という。アーム関節角度ｑｉのｉは、アーム機構の自由度の数に対応している。

又、アーム部１１２、第１腕関節１１６、第２腕関節１１７及び手首関節１１８はアーム機構を構成している。従って、アーム機構は、７自由度を備えている。前記アーム機構により、前記触覚指ベース１２０は空間運動が可能にされている。

第１触覚指〜第５触覚指１２１〜１２５は基本的に同じ機能を有する部材にて構成されているため、以下では、図１２に示す第２触覚指１２２の構成について説明し、他の触覚指については説明を省略する。

第２触覚指１２２は第１リンク１３５、第２リンク１３６、第１指関節１２７、第２指関節１２８、受動球関節１２９、及び指固定部材１３０を備えている。第１リンク１３５は第１指関節１２７を介して触覚指ベース１２０に取付けされている。第１指関節１２７は２軸の関節であって、図１２に示すように第１指関節１２７の中心において互いに直交する第１軸１２７ａ及び第２軸１２７ｂの周りでそれぞれ第１リンク１３５を能動回転する。

すなわち、第１軸１２７ａの周りで第１指関節１２７にて第１リンク１３５が回転する場合は、隣接する他の触覚指に対し近接又は離間する方向（内外転方向）に向かう回転（内外転）となる。第２軸１２７ｂの周りで第１指関節１２７にて第１リンク１３５が回転する場合は、前後屈回転となる。なお、ここで、前とは、触覚指１２２が握り方向に屈曲する方向であり、後は、その１８０度反対方向である。

又、第２リンク１３６は第１リンク１３５に対して前後屈への１自由度を有する第２指関節１２８を介して連結されている。すなわち、第２指関節１２８は、第２リンク１３６を図１２に示すように第２指関節１２８の中心線である第３軸１２８ａの周りで能動回転する。

第１指関節１２７及び第２指関節１２８は、駆動源としての駆動モータ（図示しない）をその自由度の数に一致する個数分備えており、各駆動モータを制御することにより、前述した第１軸１２７ａ〜第３軸１２８ａの周りでそれぞれ回転する。

又、第１指関節１２７及び第２指関節１２８において、前述した第１軸１２７ａ〜第３軸１２８ａを駆動する駆動モータには回転角度を検出するロータリエンコーダＵＲＥ１〜ＵＲＥ３がそれぞれ設けられている。これらのロータリエンコーダＵＲＥ１〜ＵＲＥ３により、前記各軸の周りの回転角度が検出される。

前記ロータリエンコーダＵＲＥ１〜ＵＲＥ３にて検出された回転角度は、それぞれ触覚指関節角度ｐｉｊ（ｉ＝１，…５，ｊ＝１，…３）という。なお、ｐｉｊのｉは、触覚指の数に対応する。ｐｉｊのｊは、触覚指の自由度の数に対応し、すなわち、第１軸１２７ａ〜第３軸１２８ａという軸の数（軸数）に対応する。

受動球関節１２９は、円球状をなし、第２リンク１３６の先端に半球状に凹設された取付凹部１２６ａに対して着脱自在に、かつ３６０度のいずれの範囲においても回転自在に嵌合されている。受動球関節１２９は鉄等の強磁性体からなる。永久磁石１３１は、第２リンク１３６の取付凹部１２６ａに近接した部位において第２リンク１３６に内装されている。永久磁石１３１の磁力（吸着力）により、受動球関節１２９は常には取付凹部１２６ａ内に嵌合保持されている。しかし、受動球関節１２９に、前記磁力を超える引っ張り力が加わった際には、受動球関節１２９は、取付凹部１２６ａから離脱可能である。指固定部材１３０は受動球関節１２９の露出した頂部に一体に固定されており、先端の周面は指取付部１３２が凹設されている。指取付部１３２は、長さ方向の先端側半分は断面半円状に凹設され、基端側半分は、指の先端（指先）が、受動球関節１２９の頂部に近接する位置まで挿入可能に断面円形状に形成されている。

又、指取付部１３２に対応して、指固定部材１３０の先端の周面には、リング状に形成された有弾性の締付バンド１３３が固定されている。締付バンド１３３は合成ゴム等からなるバンドにて形成されている。締付バンド１３３と指取付部１３２との間には人の人差し指が差し込み可能な空間が形成され、人差し指を差し込みした際、締付バンド１３３が同指を安定して指取付部１３２に対して添わせる。

なお、第２触覚指１２２と同様に第１触覚指１２１、第３触覚指１２３〜第５触覚指１２５の指取付部１３２及び締付バンド１３３は、それぞれ人間の拇指、中指、薬指、小指を接続可能である。

第２リンク１３６の先端周部には３軸力センサＫ２が設けられている。３軸力センサＫ２は、第２触覚指１２２に働く、互いに直交する３軸方向の並進力を検出可能である。
なお、他の触覚指の３軸力センサについては、第１触覚指〜第５触覚指１２１〜１２５の数え順の番号に一致した符号を図１３に示すようにＫに付すものとする。なお、他の触覚指におけるロータリエンコーダＵＲＥ１〜ＵＲＥ３は、図１３では、説明の便宜上、それぞれ触覚指毎にまとめて図示している。

次に多指触覚インターフェイス１１０に関する装置の電気的構成について説明する。３次元位置姿勢計測装置１４１は、多指触覚インターフェイス１１０を扱う人（以下、オペレータＨという）の手首Ｒに装着された３次元位置姿勢センサ１４２と無線で通信可能にされている。３次元位置姿勢センサ１４２は、オペレータＨの手Ｈａの位置（３次元的な空間位置）、及びその姿勢を検出し、その検出信号を３次元位置姿勢計測装置１４１に無線で通信する。

なお、前記手Ｈａの位置とは、オペレータＨの３次元位置姿勢センサ１４２を装着した手Ｈａの座標系（以下、ハンド座標系（ｘｈ，ｙｈ，ｚｈ）という）の原点位置Ｏｈを示す３つのパラメータである。又、手Ｈａの姿勢は、軸ｘｈ，ｙｈ，ｚｈが軸ｘ０，ｙ０，ｚ０へそれぞれのなす角度からなる３つのオイラーパラメータで示され、ハンド座標系の姿勢と一致する。

３次元位置姿勢計測装置１４１は、３次元位置姿勢センサ１４２から送信された検出信号を受信すると、その信号を手Ｈａの位置及びその姿勢に関するハンド座標系の位置データ及び姿勢データとして信号処理し、制御装置１４０に出力する。

制御装置１４０は、ＣＰＵ１４０ａ（中央処理装置）、ＲＯＭ１４０ｂ、ＲＡＭ１４０ｃを備えている。ＲＯＭ１４０ｂには、多指触覚インターフェイス１１０を制御するための制御プログラムや、シミュレーションプログラムを始め各種制御プログラムが格納されている。ＲＡＭ１４０ｃは、ＣＰＵ１４０ａが演算処理する際の作業用メモリである。ＣＰＵ１４０ａは、既知変位算出手段、未知変位算出手段及び反力算出手段に相当する。

第２触覚指駆動装置１４５は、制御装置１４０に電気的に接続され、制御装置１４０からの第２触覚指制御信号に基づいて、第２触覚指１２２の第１軸１２７ａ，第２軸１２７ｂを備えた第１指関節１２７、及び第３軸１２８ａを備えた第２指関節１２８を駆動する。第２触覚指駆動装置１４５は、力生成手段に相当する。

第１触覚指１２１、第３触覚指１２３〜第５触覚指１２５においても、第２触覚指１２２の第２触覚指駆動装置１４５に相当する駆動装置を備えており、制御装置１４０からの第１、第３〜第５触覚指制御信号に基づいてそれぞれの第１指関節１２７及び第２指関節１２８を同様に駆動する。第１触覚指１２１、第３触覚指１２３〜第５触覚指１２５において、第２触覚指駆動装置１４５に相当する図示しない駆動装置は、力生成手段に相当する。

なお、図１３では、説明の便宜上、第２触覚指１２２の第２触覚指駆動装置１４５等についてのみ図示し、他の触覚指に関する、駆動装置、第１指関節１２７及び第２指関節１２８については省略している。

アーム駆動装置１４６は制御装置１４０に電気的に接続され、制御装置１４０からのアーム制御信号に基づいて、第１腕関節１１６、第２腕関節１１７及び手首関節１１８を駆動する。

ディスプレイ２００は、ＣＲＴや液晶表示装置等のディスプレイ４０とからなり、制御装置１４０により、表示制御可能である。
上記のように構成された多指触覚インターフェイス１１０を使用する場合、オペレータＨは、５本の指を、第１触覚指１２１〜第５触覚指１２５の各指固定部材１３０の指取付部１３２に対して締付バンド１３３にて保持する。そして、オペレータＨが手Ｈａを動かすと、各触覚指に力が作用し、その力は、各３軸力センサＫ１〜Ｋ５にて検出される。

又、オペレータＨの手Ｈａの位置及びその姿勢は、３次元位置姿勢センサ１４２が検出し、３次元位置姿勢計測装置１４１に検出信号を送信する。そして、３次元位置姿勢計測装置１４１では、その信号をハンド座標系の位置データ及び姿勢データとして信号処理し、制御装置１４０に出力するようになっている。又、アーム機構のロータリエンコーダＡＲＥ１〜ＡＲＥ７は、アーム機構の各関節における各軸の周りの回転角度を検出する。すなわち、ロータリエンコーダＡＲＥ１〜ＡＲＥ７はアーム関節角度ｑiを検出し、制御装置１４０に出力する。又、各触覚指１２１〜１２５のロータリエンコーダＵＲＥ１〜ＵＲＥ３は、各関節の各軸の周りの触覚指関節角度ｐijを検出し、制御装置１４０に出力するようになっている。

又、アーム機構のロータリエンコーダＡＲＥ１〜ＡＲＥ７は、アーム機構の各関節における各軸の周りの回転角度を検出する。すなわち、ロータリエンコーダＡＲＥ１〜ＡＲＥ７はアーム関節角度ｑiを検出し、制御装置１４０に出力するようになっている。

なお、ｑi中、ｉ＝１，…，７は、アーム機構の第１腕関節１１６、第２腕関節１１７及び手首関節１１８（アーム用能動関節）において、７自由度を付与する各軸に対応させたものである。又、各触覚指１２１〜１２５のロータリエンコーダＵＲＥ１〜ＵＲＥ３は、各関節の各軸の周りの触覚指関節角度ｐijを検出し、制御装置１４０に出力する。又、ｐij中、ｉは、ｉ＝１，…，５で、それぞれ第１触覚指１２１〜第５触覚指１２５に対応する。又、ｊは、ｊ＝１，…，３で、第１軸１２７ａ，第２軸１２７ｂ，第３軸１２８ａにそれぞれ対応する。

なお、この多指触覚インターフェイス１１０は、オペレータＨの手Ｈａを自由に動かせるように、制御装置１４０のＣＰＵ１４０ａはアーム部１１２及び各触覚指を駆動制御する。この制御について説明する。

まず、ＣＰＵ１４０ａは、触覚指ベース１２０の位置及びその姿勢を、オペレータＨの手Ｈａの位置、及びその姿勢と対向するように制御する。又、同時に、ＣＰＵ１４０ａはオペレータＨの指先で作られる仮想平面αに対して面対称となるように、アーム駆動装置１４６にアーム制御信号を出力してアーム部１１２を駆動制御する。

すなわち、ＣＰＵ１４０ａは、３次元位置姿勢計測装置１４１から、オペレータＨの手Ｈａの位置及びその姿勢に関するハンド座標系の位置データ及び姿勢データを入力すると、これらのデータを基準座標系を介して触覚指ベース座標系に変換する。

そして、この触覚指ベース座標系において、触覚指ベース１２０の位置及びその姿勢を手Ｈａの位置及びその姿勢と対向するように、オペレータＨの指先で作られる仮想平面αに対して面対称となるように、アーム駆動装置１４６にアーム制御信号を出力してアーム部１１２を駆動制御する。

前記アーム部１１２、第１腕関節１１６、第２腕関節１１７及び手首関節１１８からなるアーム機構は７自由度を有するため、この対向はアーム機構の動作可能範囲であれば常に実現できる。

ここで、仮想平面αについて説明する。仮想平面αは、オペレータＨの拇指、人差し指、中指の３本の指先で作られる平面を基本としている。なお、これら３本の指がほぼ直線状に並ぶときは、他の指を含め、最小２乗法で指先を乗せる平面を決定している。

又、手Ｈａの指先の位置は、第１触覚指１２１〜第５触覚指１２５の各指固定部材１３０の指取付部１３２に保持されていることから、各触覚指の先端位置（指取付部１３２の位置）と同じである。このため、ＣＰＵ１４０ａは仮想平面αを決定するために必要な各触覚指の先端位置（３次元位置）を、触覚指の触覚指関節角度ｐij、触覚指の長さ、アーム機構のアーム関節角度ｑi、アーム機構の長さに基づいて算出する。この各触覚指の先端位置（指取付部１３２の位置）の位置情報が、第１実施形態で説明した「作用点」の位置情報に相当する。

なお、各触覚指の触覚指関節角度ｐij及びアーム機構のアーム関節角度ｑiは、各触覚指のロータリエンコーダＵＲＥ１〜ＵＲＥ３の検出値及びアーム機構のロータリエンコーダＡＲＥ１〜ＡＲＥ７の検出値に基づいてそれぞれ算出する。又、各触覚指の長さ及びアーム機構の長さは、触覚指及びアーム機構を構成している各部材の長さであり、これらの値は予めＲＯＭ１４０ｂに格納されている。

さて、上記のように構成された多指触覚インターフェイス１１０を利用して、仮想空間での仮想指先と、仮想柔軟物体との接触が第１実施形態と同様に生じた場合、第１実施形態と同様のシミュレーションプログラムに従い、仮想柔軟物体の変形のための演算を行う。この処理は第１実施形態の図４のフローチャートと同様であるため、説明を省略するが、例えば、各触覚指の先端位置（すなわち、作用点の位置情報）が、図４のフローチャートにおいて、Ｓ１０で仮想指先の位置の更新に使用される。又、Ｓ２０では、各触覚指の先端位置（すなわち、作用点の位置情報）が仮想柔軟物体との衝突判定に使用される。

上記実施形態、実施例１、及び実施例２は下記の効果を奏する。
○ 前記実施形態では、全体の剛性方程式における物体Ｂ（有限要素モデル）の全ノードに関する変位を、自由ノードの未知変位ｕ_ｕと、指先ＵＳ（拘束物）に接触する接触面から、１つ抽出された拘束面に関する既知変位ｕ_ｋと、拘束面以外の残りの接触面に関し、等価荷重Ｆ^＊で変換可能な既知変位ｕ_ｋ ^＊とに区分するようにした。そして、物体Ｂの全体の剛性方程式から算出される剛性行列Ａの逆行列Ａ^−１を、ウッドバリーの公式が適用できるように予め変形した。さらに、前記未知変位ｕ_ｕを、ウッドバリーの公式に基づいて算出するようにした。このようにすることにより、空間で拘束条件を欠く物体の掴み操作を行う場合に、物体の変形を、有限要素法を利用して算出するに当たり、短時間で行うことができる。

○ 実施例１では、ウッドバリーの公式には、対角行列を含み、該対角行列は、接触面に関与している対角要素は１，それ以外は０で構成しているため、乗算計算量を少なくでき、物体の変形を、有限要素法を利用して算出するに当たり、短時間で行うことができる。

○ 実施例２では、ウッドバリーの公式には、対角行列を含み、該対角行列は、対角要素を全て１としているため、乗算計算量を少なくでき、物体の変形を、有限要素法を利用して算出するに当たり、短時間で行うことができる。

○ 前記実施形態では、未知変位を算出した後、指先ＵＳへの反力を算出することにより、短時間で未知変位が算出されるため、拘束物への反力の算出するまでの時間も短縮することができる。

○ 第１実施形態のシミュレーション装置１００は、ＰＣ３０Ｂは、入力された指先ＵＳの位置情報と、有限要素モデルを構成するノードの位置情報に基づいて、指先ＵＳに接触する接触面から、１つ抽出された拘束面に関する既知変位と、前記拘束面以外の残りの接触面に関し、等価荷重で変換可能な既知変位を算出するようにした。又、ＰＣ３０Ｂは、前記有限要素モデルの全体の剛性方程式から算出された剛性行列の逆行列がウッドバリーの公式が適用できるように予め変形され、前記剛性方程式から算出された変形式を使用して、前記未知変位を、前記ウッドバリーの公式に基づいて算出するようにした。この結果、物体の変形を、有限要素法を利用して算出するに当たり、短時間で行うことができ、シミュレーションを短時間で行うことができる。

○ 又、第１実施形態のシミュレーション装置１００は、未知変位ｕ_ｕを算出した後、指先ＵＳへの反力を算出するようにした。この結果、短時間で未知変位が算出されるため、拘束物への反力の算出するまでの時間も短縮することができる。

なお、本発明は前記実施形態に限定されるものではない。下記のように具体化してもよい。
○ 第２実施形態では５本指を有するハンド１１３を有するマンマシンインタフェイスとしたが、グローブ形のマンマシンインタフェイスとしてもよい。例えば、このようなものとして、図１４に示すマンマシンインタフェイス３００を挙げることができる。

マンマシンインタフェイス３００は、グローブ体３１０、コンピュータ、ディスプレイ装置（ともに図示しない）を備えている。グローブ体３１０の側部には、各指に対してワイヤＹ１〜Ｙ４を取り付け支持する取付部３２０〜３５０が設けられている。各取付部３２０〜３５０の背面側（図１４において、上面側）には、ワイヤＹ１〜Ｙ４をガイド支持する支持部材３６０ａ〜３６０ｄが設けられている。ワイヤＹ１〜Ｙ４は、それぞれ親指，人指し指、中指、薬指用のものである。ワイヤＹ１〜Ｙ４の各先端は、各取付部３２０〜３５０の先端に取付けされている。ワイヤＹ１〜Ｙ４の各基端は、図示はしないが、ワイヤ毎に設けられた力生成部に連結されている。そして、該力生成部が、ワイヤＹ１〜Ｙ４を引っ張ることにより、取付部３２０〜３５０の先端に填められた人の指先を背面側から引っ張るようにされている。前記力生成部は、力生成手段に相当する。又、力生成部は、たとえば、第２実施形態の制御装置１４０に相当する図示しない前記コンピュータに接続されている。該コンピュータは、図示しないディスプレイ装置に接続されている。

又、グローブ体３１０には、手首に対応して位置するように３次元位置姿勢センサ（図示しない）が設けられ、オペレータの手の位置及びその姿勢は、３次元位置姿勢センサが検出可能である。

又、取付部３２０〜３５０には、図示はしないが、各指の３次元位置を検出する３次元位置センサが備えられ、前記コンピュータに電気的に接続されている。
３次元位置姿勢センサ、３次元位置センサの検出信号に基づいて、手の位置、姿勢、及び指の３次元位置が算出可能である。そして、前記コンピュータを、既知変位算出手段、未知変位算出手段及び反力算出手段として機能させることにより、マンマシンインタフェイス３００を変位シミュレーション装置や、力覚提示装置として具体化することが可能である。

上記のように、本実施形態と、第２実施形態の構成の相違により、機械的作用が下記のように異なっている。すなわち、前記第２実施形態がオペレータの指先に対して反力を付与する場合、指先の頭を押圧するようにされているのに対して、本実施形態では、オペレータの指先の背面側をワイヤＹ１〜Ｙ４にて引っ張るようにしているのが異なっている。

○ 第１実施形態では、シミュレーション装置１００に一対の力覚提示ユニット２０Ａ，２０Ｂを設けたが、一対に限定するものではなく、３つ以上の力覚提示装置を設けてもよい。

○ 又、３つ以上の力覚提示装置を設ける場合、各力覚提示装置に対して、それぞれＰＣを接続するようにしてもよい。
○ 第１実施形態では、シミュレーション装置１００は、一対のＰＣ３０Ａ，３０Ｂを設けたが、１つのＰＣで構成してもよい。

○ 前記各実施形態においては、仮想柔軟物体に対して、拘束物としての指による掴みについて説明したが、本発明は、掴み操作に限定されるものではなく、仮想柔軟物体を、複数の拘束物が挟むように操作されるものについても適用できる。

（ａ）〜（ｃ）は物体に対する指先の接触状態の説明図。 ３角形パッチの説明図。 シミュレーション装置の電気的構成を示すブロック図。 シミュレーション装置が処理するシミュレーションのフローチャート。 物体に対する指先の接触状態の説明図。 拘束面の説明図。 行列サイズに応じて、計算量がどの程度変化するかの計算量を示すグラフ。 （ａ）は比較例２のｘ方向の侵入量を示す時間変位グラフ、（ｂ）は、比較例２の、反力ｆ１，ｆ２の時間変位グラフ。 （ａ）は実施例１のｘ方向の侵入量を示す時間変位グラフ、（ｂ）は、実施例１の、反力ｆ１，ｆ２の時間変位グラフ。 （ａ）〜（ｄ）は、実施例１のシミュレーション装置において、物体Ｂを掴んだ状態で、指先ＵＳの拘束位置を上から下へ移動させたときの０．５［sec］毎の物体Ｂの変形状態の説明図。 他の実施形態におけるシミュレーション装置及び力覚提示装置としての多指触覚インターフェイスの概略を示す斜視図。 同じく力覚提示装置の第２触覚指の概略図。 同じく力覚提示装置の電気的構成を示すブロック図。 他の実施形態の力覚提示装置の概略斜視図。

符号の説明

２６…力生成部（力生成手段）
３０Ｂ…ＰＣ（既知変位算出手段、未知変位算出手段及び反力算出手段）
１００…シミュレーション装置（力覚提示装置）
１４０ａ…ＣＰＵ（既知変位算出手段、未知変位算出手段及び反力算出手段）
１４５…第２触覚指駆動装置（力生成手段）
ＵＳ…指先（拘束物）
Ｂ…物体

Claims (7)

1. 空間で拘束条件を欠く物体を拘束する拘束物が前記物体に対して移動することによって発生する前記物体の変形を、前記物体に関する有限要素モデルについての全体の剛性方程式に基づいて演算するシミュレーションを行う変位シミュレーション方法において、
   前記全体の剛性方程式における前記有限要素モデルの全ノードに関する変位を、自由ノードの未知変位と、前記拘束物に接触する接触面から、一つ抽出された拘束面に関する既知変位と、前記拘束面以外の残りの接触面に関し、等価荷重で変換可能な既知変位とを算出し、
   前記有限要素モデルの全体の剛性方程式から算出される剛性行列の逆行列を、予め変形することにより算出した変形式に対して、ウッドバリーの公式適用することにより、前記全ノードに関する変位のうち既知変位以外の変位である未知変位を算出することを特徴とする物体の変位シミュレーション方法。
2. 前記ウッドバリーの公式には、対角行列を含み、該対角行列は、接触面に関与している対角要素が１，それ以外が０で構成することを特徴とする請求項１に記載の物体の変位シミュレーション方法。
3. 前記ウッドバリーの公式には、対角行列を含み、該対角行列は、対角要素を全て１とすることを特徴とする請求項１に記載の物体の変位シミュレーション方法。
4. 前記未知変位を算出した後、前記拘束物への反力を算出することを特徴とする請求項１乃至請求項３のうちいずれか１項に記載の物体の変位シミュレーション方法。
5. 空間で拘束条件を欠く物体を拘束する拘束物が前記物体に対して移動することによって発生する前記物体の変形を、前記物体に関する有限要素モデルについての全体の剛性方程式に基づいて演算するシミュレーションを行う物体の変位シミュレーション装置において、
   入力された前記拘束物の位置情報と、前記有限要素モデルを構成するノードの位置情報に基づいて、前記全体の剛性方程式における前記有限要素モデルの全ノードに関する変位のうち、前記拘束物に接触する接触面から一つ抽出された拘束面に関する既知変位、及び、前記拘束面以外の残りの接触面に関し、等価荷重で変換可能な既知変位の両既知変位を算出する既知変位算出手段と、
   前記有限要素モデルの全体の剛性方程式から算出される剛性行列の逆行列を予め変形することにより算出した変形式に対して、ウッドバリーの公式適用することにより、前記全ノードに関する変位のうち既知変位以外の変位である未知変位を算出する未知変位算出手段を備えていることを特徴とする物体の変位シミュレーション装置。
6. 前記算出された未知変位に基づいて前記拘束物への反力を算出する反力算出手段を備えたことを特徴とする請求項５に記載の物体の変位シミュレーション装置。
7. 空間で拘束条件を欠く物体を拘束する拘束物が前記物体に対して移動することによって発生する前記物体の変形を、前記物体に関する有限要素モデルについての全体の剛性方程式に基づいて演算するシミュレーションを行い、力覚を提示する力覚提示装置において、
   入力された前記拘束物の位置情報と、前記有限要素モデルを構成するノードの位置情報に基づいて、前記全体の剛性方程式における前記有限要素モデルの全ノードに関する変位のうち、前記拘束物に接触する接触面から一つ抽出された拘束面に関する既知変位、及び、前記拘束面以外の残りの接触面に関し、等価荷重で変換可能な既知変位の両既知変位を算出する既知変位算出手段と、
   前記有限要素モデルの全体の剛性方程式から算出される剛性行列の逆行列を予め変形することにより算出した変形式に対して、ウッドバリーの公式適用することにより、前記全ノードに関する変位のうち既知変位以外の変位である未知変位を算出する未知変位算出手段と、
   前記算出された未知変位に基づいて前記拘束物への反力を算出する反力算出手段と、
   前記算出された反力に基づいて力を発生する力生成手段を備えたことを特徴とする力覚提示装置。

Images