【手続補正書】
【提出日】平成16年7月29日(2004.7.29)
【手続補正１】
【補正対象書類名】明細書
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【書類名】                  明細書
【発明の名称】  マルチボディシステムとして取扱い可能な可動体の解析処理装置および制御方法
【特許請求の範囲】
【請求項１】  等式拘束で表わされるホロノミック拘束と非ホロノミック拘束からなるマルチボディシステムとして取扱い可能な可動体の解析処理装置において、ホロノミック等式拘束を非ホロノミック等式拘束の形式に変換して定式化し、また両者の拘束誤差を１個のスカラーのみで定義して統合し、マルチボディシステムの解析処理をすることを特徴とするマルチボディシステムとして取扱い可能な可動体の解析処理装置。
【請求項２】  請求項１において、拘束誤差をスカラー（Ｊ）で定義する時に、拘束誤差を行列（Ｓ）で重み付けすることを特徴とするマルチボディシステムとして取扱い可能な可動体の解析処理装置。
【請求項３】  請求項１において、拘束力を可変ゲイン誤差フィードバックにより補正することを特徴とするマルチボディシステムとして取扱い可能な可動体の解析処理装置。
【請求項４】  請求項１において、拘束の安定化にかかわるパラメータを運動の状態に適応させて自動的に変化させることを特徴とするマルチボディシステムとして取扱い可能な可動体の解析処理装置。
【請求項５】  ホロノミック等式拘束と非ホロノミック等式拘束および不等式拘束からなるマルチボディシステムとして取扱い可能な可動体の解析処理装置において、まず不等式拘束をホロノミック等式拘束または非ホロノミック等式拘束の形式に変換し、次にホロノミック等式拘束を非ホロノミック等式拘束の形式に変換して定式化し、また各拘束誤差を１個のスカラーのみで定義して統合することを特徴とするマルチボディシステムとして取扱い可能な可動体の解析処理装置。
【請求項６】  請求項１ないし請求項５のいずれかにおいて、可動体は、多関節体または柔軟体であることを特徴とするマルチボディシステムとして取扱い可能な可動体の解析処理装置。
【請求項７】  等式拘束で表わされるホロノミック拘束と非ホロノミック拘束からなるマルチボディシステムとして取扱い可能な可動体の制御方法において、ホロノミック等式拘束を非ホロノミック等式拘束の形式に変換して定式化し、また両者の拘束誤差を１個のスカラーのみで定義して統合し、マルチボディシステムを制御することを特徴とするマルチボディシステムとして取扱い可能な可動体の制御方法。
【請求項８】  請求項７において、拘束誤差をスカラー（Ｊ）で定義する時に、拘束誤差を行列（Ｓ）で重み付けすることを特徴とするマルチボディシステムとして取扱い可能な可動体の制御方法。
【請求項９】  請求項７において、拘束力を可変ゲイン誤差フィードバックにより補正することを特徴とするマルチボディシステムとして取扱い可能な可動体の制御方法。
【請求項１０】  請求項７において、拘束の安定化にかかわるパラメータを運動の状態に適応させて自動的に変化させることを特徴とするマルチボディシステムとして取扱い可能な可動体の制御方法。
【請求項１１】  ホロノミック等式拘束と非ホロノミック等式拘束および不等式拘束からなるマルチボディシステムとして取扱い可能な可動体の制御方法において、まず不等式拘束をホロノミック等式拘束または非ホロノミック等式拘束の形式に変換し、次にホロノミック等式拘束を非ホロノミック等式拘束の形式に変換して定式化し、また各拘束誤差を１個のスカラーのみで定義して統合することを特徴とするマルチボディシステムとして取扱い可能な可動体の制御方法。
【請求項１２】  請求項７ないし請求項１１のいずれかにおいて、可動体は、多関節体または柔軟体であることを特徴とするマルチボディシステムとして取扱い可能な可動体の制御方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、等式拘束で表わされるホロノミック拘束と非ホロノミック拘束からなるマルチボディシステムとして取扱い可能な可動体の解析処理装置および制御方法に関するものであり、たとえば複数の機械要素がジョイントで結合されたマルチボディシステムや離散化処理可能な運動方程式を持つ柔軟体をトータルに取り扱って設計あるいは制御を行うことを可能にする。
【０００２】
【従来の技術】
多くの機械は、複数の構成要素（ボディ）が幾何学的拘束（ジョイント）により結合されたマルチボディシステム（多体系）として捉える事ができる。典型的な多体系は、関節を有するリンク機構や、連続体を有限要素などの集合として離散化したモデル等である。  マルチボディシステムの設計や制御に際して、それらの構成要素の運動や変形を計算により予測することは不可欠であるが、この種のシステムでは、一般的なニュートンやラグランジュの力学原理を直接適用して動力学的な定式化により運動を予測するのは困難である。その主たる要因は以下のようなものである。
（Ａ）力学原理から直接得られる運動方程式は、システムの自由度の３乗に比例して計算量が増加するため、複雑なシステムでは実際上計算が困難になり、原理的に扱えなくなる。
（Ｂ）運動方程式が数値的に不安定になる場合が多く、そのままでは解が発散して、意味のある結果が得られない。
（Ｃ）幾何学的拘束が特異になると、拘束力を一意に計算できない。
（Ｄ）幾何学的拘束は、一般化速度を含まないホロノミック拘束とそれ以外の非ホロノミック拘束に分類される。ホロノミック拘束は、拘束条件が代数方程式や積分可能な微分方程式で表わされるものであり、独立な一般化座標を選択する事で比較的容易に運動方程式を導出できる場合が多い。一方、非ホロノミック拘束は、拘束条件が不等式や積分できない微分方程式で表わされるなど、普遍的な数学的記述法すら存在していないものである。その中で典型的なのは速度を含んだ等式として表される拘束であるが、ラッチ機構の様に不等式拘束となる場合も有る。非ホロノミック等式拘束は通常ラグランジュの未定乗数によって処理され、その結果ダイナミクスは微分−代数方程式で記述されることになるが、これは一般にスティフで数値的不安定に陥りやすい。伝統的な理論ではこれらのホロノミックと非ホロノミックの拘束を別個に取り扱うので、両者が混在する現実のシステムを効率よくモデル化できない。
【０００３】
このようなマルチボディシステムの背景のもとで、これまでに開発された主な解析手法と問題点を次に示す。
（1）  拘束の誤差ベクトルが指数減衰するように設定する方法（参考文献１，２）。−→上記（Ｂ）を改善するもので、他の問題は解決されない。
（2）  運動方程式を漸化的に解くことで、計算量をシステム自由度の１乗のオーダに減らす方法（参考文献６）。−→上記（Ａ）以外は改善されない。（Ａ）についても、拘束が複雑な場合には改善効果が十分に発揮されず、また並列計算ができないので、さらなる効率化には障害となることもある。
（3）  拘束の誤差ベクトルに「定数」を乗じて拘束力の未定乗数とするペナルティ法（参考文献３〜５）。−→上記（Ａ）、（Ｃ）は通常改善される。しかし（Ｂ）の問題に対しては、誤差の挙動が予め分からないので、定数の値が不適切であると発散する場合もあることから、完全な解決にはならない。
（4）上記ペナルティ法の「定数」を可変とし、拘束の誤差ノルムを指数減衰させる事で安定性を保証する本発明者により提案された可変ゲインペナルティ法（参考文献７）。−→上記の問題は概ね改善される。しかし、数値計算の精度と効率が不十分である。
（4）の可変ゲインペナルティ法は、拘束の安定性が保証され、誤差の目標値も確実に維持され、プログラミングが容易で数値計算コストが低く、並列計算に適している上、ホロノミック・非ホロノミックを問わず全ての等式拘束を統一的に扱える利点をもつ。しかし現実のシステムには、不等式で表現される幾何学的拘束を含むものも存在する。ラッチ機構の様に、一般化座標が或る範囲内の間は自由で、それを越えると拘束される場合がその典型である。可変ゲインペナルティ法を含む既存の定式化の殆どは、この種のシステムを直接的にモデル化する事ができない。最も簡単な方法として、不等式の境界で剛性が大きく変化する仮想的なバネなどで表現する事も可能であるが、これは一般に運動方程式をスティフにするため、計算コストや数値的安定性の観点から実用的とは言えない。不等式拘束の中でも境界で反射が起きる時は、その前後の不連続な速度変化が表面の材質等に大きく依存するので、動力学的にモデル化せざるを得ない。しかしラッチ機構の如く不等式境界に捕捉される場合は純粋な幾何学的拘束として扱う方が遥かに効率的であり、一般性も高い。
【０００４】
【発明が解決しようとする課題】
この様に、既存の手法は何らかの形で運動方程式を緩和するため、数値解がエネルギーや運動量等の観点から見て力学的にどの程度妥当なのか厳しく検証せねばならない。しかし関連する論文の多くは拘束の安定性ばかりに注目し、この点を疎かにしている感が有る。またアルゴリズムのパラメータと実際に得られる解の精度との関係が複雑で決定に経験を要求される。しかもほとんどの場合ホロノミック拘束と非ホロノミック拘束で定式化が異なるため、両者が混在するシステムを扱うのが困難である。従って実用上の制約も大きい。
【０００５】
本発明の課題は、ホロノミック拘束と非ホロノミック拘束を含むマルチボディシステムを統一的な運動方程式で取り扱える手法を実現することにある。
【０００６】
【課題を解決するための手段】
本発明は、非常に広範囲の剛体／柔軟マルチボディシステムに適用できる新たな一定式化手法を提供する。
【０００７】
本発明は、ホロノミック系要素と非ホロノミック系要素が混在するマルチボディシステムとして取扱い可能な可動体の解析処理装置および制御方法として実現され、以下のように構成される。
（１）  等式拘束で表わされるホロノミック拘束と非ホロノミック拘束からなるマルチボディシステムとして取扱い可能な可動体の解析処理装置および制御方法において、ホロノミック等式拘束を非ホロノミック等式拘束の形式に変換して定式化し、また両者の拘束誤差を１個のスカラーのみで定義して統合し、マルチボディシステムの解析処理をすることを特徴とするマルチボディシステムとして取扱い可能な可動体の解析処理装置および制御方法の構成。
（２）  前項（１）において、拘束誤差をスカラー（Ｊ）で定義する時に、拘束誤差を行列（Ｓ）で重み付けすることを特徴とするマルチボディシステムとして取扱い可能な可動体の解析処理装置および制御方法の構成。
（３）  前項（１）において、拘束力を可変ゲイン誤差フィードバックにより補正することを特徴とするマルチボディシステムとして取扱い可能な可動体の解析処理装置および制御方法の構成。
（４）  前項（１）において、拘束の安定化にかかわるパラメータを運動の状態に適応させて自動的に変化させることを特徴とするマルチボディシステムとして取扱い可能な可動体の解析処理装置および制御方法の構成。
（５）  ホロノミック等式拘束と非ホロノミック等式拘束および不等式拘束からなるマルチボディシステムとして取扱い可能な可動体の解析処理装置および制御方法において、まず不等式拘束をホロノミック等式拘束または非ホロノミック等式拘束の形式に変換し、次にホロノミック等式拘束を非ホロノミック等式拘束の形式に変換して定式化し、また各拘束誤差を１個のスカラーのみで定義して統合することを特徴とするマルチボディシステムとして取扱い可能な可動体の解析処理装置および制御方法の構成。
（６）  前項（１）ないし（５）のいずれかにおいて、可動体は、多関節体または柔軟体であることを特徴とするマルチボディシステムとして取扱い可能な可動体の解析処理装置および制御方法の構成。
【０００８】
図１は、上記した本発明の構成によりマルチボディシステムを定式化する処理の概略フローである。
【０００９】
図１において、ステップ（Ｓ−１）では、マルチボディシステムの構成要素中に、ラッチ機構のような不等式拘束が含まれているかどうかを調べ、不等式拘束が含まれている場合には、ステップ（Ｓ−２）で、不等式拘束をホロノミック等式拘束か非ホロノミック等式拘束に変換し、マルチボディシステムから不等式拘束形式の表現を取り除く。このとき拘束誤差は，スカラーで定義する。次にステップ（Ｓ−３）で、マルチボディシステム中にホロノミック等式拘束と非ホロノミック等式拘束が混在しているかどうかを調べ、混在している場合には、ステップ（Ｓ−４）で、ホロノミック等式拘束を非ホロノミック等式拘束の形式に変換し、拘束誤差をスカラーで定義する。このようにして、マルチボディシステム全体を、非ホロノミック等式拘束の表現で統一し、定式化する。
【００１０】
次に、本発明による定式化処理の理論的根拠について、以下に述べる。
１．等式拘束
１. １  運動方程式
ＮB 個のボディから成る多体システムを定式化することを考える。幾何学的拘束を付加する前は、システムを構成する個々のボディの運動方程式が以下の様に表せる（参考文献７）。
【００１１】
【数１】
2002288237.tif 000001
【００１２】
但しここで、
【００１３】
【数２】
2002288237.tif 000002
【００１４】
が第ｉボディ（Ｂi ）の一般化速度ベクトルであり、一般化座標の変化率
【００１５】
【数３】
2002288237.tif 000003
【００１６】
でも良いし、またこれに関して線形で相互に変換可能な関数であれば何でも構わない。
【００１７】
【数４】
2002288237.tif 000004
【００１８】
は対応する一般化力、
【００１９】
【数５】
2002288237.tif 000005
【００２０】
はそれ以外の項で、遠心力やコリオリ力等が含まれる。そして、これらを合わせた力の総計を
【００２１】
【数６】
2002288237.tif 000006
【００２２】
で表している。左辺のＭi は質量行列である。どの様な柔軟体であっても、一般化速度ベクトル
【００２３】
【数７】
2002288237.tif 000007
【００２４】
を選択し、ボディ内に在る任意の微小質量ｄｍの位置ベクトル
【００２５】
【数８】
2002288237.tif 000008
【００２６】
とこれに作用する力のベクトル
【００２７】
【数９】
2002288237.tif 000009
【００２８】
（共に慣性系表示とする）
が与えられたとすると、これらの関数は次式より計算できる。
【００２９】
【数１０】
2002288237.tif 000010
【００３０】
但し上式を含めここでは一貫して、ベクトル関数
【００３１】
【数１１】
2002288237.tif 000011
【００３２】
のベクトル
【００３３】
【数１２】
2002288237.tif 000012
【００３４】
に関するヤコビ行列
【００３５】
【数１３】
2002288237.tif 000013
【００３６】
を
【００３７】
【数１４】
2002288237.tif 000014
【００３８】
と略記し、また
【００３９】
【数１５】
2002288237.tif 000015
【００４０】
が零の時に
【００４１】
【数１６】
2002288237.tif 000016
【００４２】
が取る値を
【００４３】
【数１７】
2002288237.tif 000017
【００４４】
で表す。
【００４５】
幾何学的拘束を付加する前のシステム方程式は、以下の如く各ボディの方程式を並べるだけで簡単に得られる。
【００４６】
【数１８】
2002288237.tif 000018
【００４７】
システムに幾何学的拘束が与えられると、これを満足させるための拘束力が必要となるので、運動方程式が以下の様に変化する。
【００４８】
【数１９】
2002288237.tif 000019
【００４９】
上式の
【００５０】
【数２０】
2002288237.tif 000020
【００５１】
が拘束力である。
１．２  幾何学的拘束の緩和
ここでは、次式で与えられる等式拘束を考える。
【００５２】
【数２１】
2002288237.tif 000021
【００５３】
但し、
【００５４】
【数２２】
2002288237.tif 000022
【００５５】
は式（３ｃ）と同様に定義されたシステムの一般化座標ベクトルであり、
【００５６】
【数２３】
2002288237.tif 000023
【００５７】
は零ベクトルである。通常式(5a)と式(5b）は各々非ホロノミック等式拘束及びホロノミック拘束に対応する。但し後者も、式(5a)の形に容易に変換できる。例えば時間の逆数の次元を持つ適当な正の定数ωを用いて
【００５８】
【数２４】
2002288237.tif 000024
【００５９】
と置けば、右辺が一般化速度
【００６０】
【数２５】
2002288237.tif 000025
【００６１】
を含み、かつ
【００６２】
【数２６】
2002288237.tif 000026
【００６３】
で
【００６４】
【数２７】
2002288237.tif 000027
【００６５】
となるので、非ホロノミック等式拘束として扱うことができる。更に第２節で詳述するようにラッチ型の不等式拘束も式(5a)または式(5b)の形で表現できる。よって以降は式(5a)の非ホロノミック等式拘束だけを仮定して説明を進める。
【００６６】
式(5a)の左辺の
【００６７】
【数２８】
2002288237.tif 000028
【００６８】
は拘束の誤差ベクトルであって、理想的な状態では常に零であるが、これは数値的な困難を生じるので、微小な値は許すものとする。そのため、ここでは誤差ベクトルの重み付ノルム
【００６９】
【数２９】
2002288237.tif 000029
【００７０】
（以下誤差ノルムと呼ぶ）を適当な許容値に維持する方法を考える。但し上式のＳは、数値的特性を良好にするための正定置な対角行列（定数）である。その導出は１．４節で述べるが、これを以下スケーリング行列と呼ぶ。
【００７１】
さて、ノルムＪを収束させる方法は無数に考えられるが、最も簡単なのは、以下の微分方程式を設定することである。
【００７２】
【数３０】
2002288237.tif 000030
【００７３】
ここでωは式(6) にも登場した時間の逆数の次元を持つ定数である。これを拘束の減衰係数と名付ける。上式が満足されれば誤差ノルムが零に向かって指数的に減少するが、ωはその時定数の逆数に相当する。パラメータωは拘束の固有値に相当し、簡便には定数としても構わない。しかし数値積分のステップ幅Δｔにほぼ反比例するシステムの最大固有値より少し大きめの値を維持すると、計算の精度・コスト両面で有利である。従って、次式を使って適当な間隔でωを修正する事が推奨される。
【００７４】
ω＝β／Δｔ      （８ａ）
但しβは、積分法を選択すると決まる定数である。
１．３  拘束力の緩和
前節で緩和された幾何学的拘束において、誤差ベクトルのノルムが式(8) に従うという条件は与えられているが、ベクトルの各成分には特別な制約が無いため、このままでは式(4) の拘束力
【００７５】
【数３１】
2002288237.tif 000031
【００７６】
が一意に決まらない。そこで、以下に挙げる２つの条件を拘束力に付加する。その第１は殆どすべての研究者が仮定しているが、拘束力は緩和しても理想的な拘束力と同様に仮想変位と直交するというものである。すると拘束力は、適当な未定乗数
【００７７】
【数３２】
2002288237.tif 000032
【００７８】
を用いて以下の様に表せる。
【００７９】
【数３３】
2002288237.tif 000033
【００８０】
そして第２の条件として、上式における未定乗数
【００８１】
【数３４】
2002288237.tif 000034
【００８２】
が常に誤差ベクトルと次式の関係に在ると仮定する。
【００８３】
【数３５】
2002288237.tif 000035
【００８４】
上式の第一項は未定乗数
【００８５】
【数３６】
2002288237.tif 000036
【００８６】
の推定値であり、初期時刻で零とし、以降は直前の積分ステップにおける値を用いる。第二項はこれを修正する効果を持つが、ここでＳは式(7) で導入したスケーリング行列で、
【００８７】
【数３７】
2002288237.tif 000037
【００８８】
はスケーリングされた拘束誤差ベクトルである。またαは時変のスカラーであり、これをここではゲインと呼ぶ事にする。これを定数と置くのが通常のペナルティ法（参考論文３〜５）であるが、拘束を安定化するゲインの最小値は一般に変動するので、余程大きな値にしておかないと危険である。しかしこれは運動方程式をスティフにし、無用な計算コストを増大させる。またゲインを定数とすると、仮に安定であっても一般に拘束誤差は指数減衰せず、より複雑な履歴をたどるため、計算効率も提案する可変ゲインの場合より劣ると考えられる。
【００８９】
式(10)に式(11)を代入すると、上記の２条件を満たす拘束力は次の形で与えられる事が分かる。
【００９０】
【数３８】
2002288237.tif 000038
【００９１】
但しここで、以下の記号を定義した。
【００９２】
【数３９】
2002288237.tif 000039
【００９３】
前者は修正前の拘束力である。また後者は単位ゲイン当りに修正される拘束力に相当するので、これを拘束力修正方向ベクトルと名付けるものとする。式(13)の拘束力を式(4) の運動方程式に代入し、一般化加速度
【００９４】
【数４０】
2002288237.tif 000040
【００９５】
に関して解くと、
【００９６】
【数４１】
2002288237.tif 000041
【００９７】
の様になる。但しここで、
【００９８】
【数４２】
2002288237.tif 000042
【００９９】
なる記号を導入した。上式の
【０１００】
【数４３】
2002288237.tif 000043
【０１０１】
は拘束力の修正が無い場合の一般化加速度であり、
【０１０２】
【数４４】
2002288237.tif 000044
【０１０３】
は単位ゲイン当りの一般化加速度の修正量に相当する。
【０１０４】
ゲインαを適切に与えれば、誤差ノルムを式(8) に従って指数的に収束させられる事を次に示す。式(5a)左辺の誤差ベクトルを時間で微分すると、次式を得る。
【０１０５】
【数４５】
2002288237.tif 000045
【０１０６】
特にホロノミック拘束を式(6) で変換した場合、上式右辺の諸量はその定義に従うと
【０１０７】
【数４６】
2002288237.tif 000046
【０１０８】
で与えられる。式(15)を式(17)に代入すると、拘束誤差は以下の様なダイナミクスを持つ事が分かる。
【０１０９】
【数４７】
2002288237.tif 000047
【０１１０】
式（8 ）を満たすゲインαを得るため、式(7) を
【０１１１】
【数４８】
2002288237.tif 000048
【０１１２】
の如く時間微分する。これを式（8)に代入し分母を払うと、
【０１１３】
【数４９】
2002288237.tif 000049
【０１１４】
となる。上式左辺の第１項は式（19）および式（14）より
【０１１５】
【数５０】
2002288237.tif 000050
【０１１６】
と置けるので、これを代入して整理すると、ゲインαに関する１次方程式
【０１１７】
【数５１】
2002288237.tif 000051
【０１１８】
を得る。但しここで、以下のスカラ一量を導入した。
【０１１９】
【数５２】
2002288237.tif 000052
【０１２０】
式（２４ａ）の係数Ｐは、式（16b ）より
【０１２１】
【数５３】
2002288237.tif 000053
【０１２２】
であり、質量行列Ｍが正定値である事を考えると、拘束力方向ベクトル
【０１２３】
【数５４】
2002288237.tif 000054
【０１２４】
が零でない限り常に正となる事に注意されたい。これは拘束の正則性よりかなり緩い条件である。何故なら式（14）においてヤコビ行列
【０１２５】
【数５５】
2002288237.tif 000055
【０１２６】
のランク落ちが生じて拘束が特異になっても、
【０１２７】
【数５６】
2002288237.tif 000056
【０１２８】
が偶然その零空間内に在って
【０１２９】
【数５７】
2002288237.tif 000057
【０１３０】
となる事は実用上まず無いからである。よって本手法は特異な拘束に対して頑健である事が分かる。この特性は後述の数値例で実証する。式（14）より拘束誤差が全く無い場合には拘束力方向ベクトルも零となるため式（23）からゲインαを求める事はできないが、拘束が完璧に満たされている瞬間は拘束力も修正する必要がないので、ゲインも零としておけば良い。即ち、
【０１３１】
【数５８】
2002288237.tif 000058
【０１３２】
より簡単にゲインαが得られる。但し実際の計算においては、拘束誤差が零でなくても非常に小さく機械精度の限界に近いときは、上式右辺の分母Ｐが有効桁を失って意味の無い過大な拘束力が出る恐れが有るため、強制的にα＝０とするべきである。
【０１３３】
結局、順ダイナミクスに必要な一般化加速度
【０１３４】
【数５９】
2002288237.tif 000059
【０１３５】
の計算は、以下の手順（ａ）〜（ｅ）で行えば良い事が分かる。
（ａ）システムの状態変数と時刻ｔから、運動方程式（4)におけるシステムの質量行列Ｍ（式（3b）) 
と拘束力以外の力
【０１３６】
【数６０】
2002288237.tif 000060
【０１３７】
を計算し、これより拘束力の修正が無い場合の一般化加速度
【０１３８】
【数６１】
2002288237.tif 000061
【０１３９】
を求める。
（ｂ）拘束誤差ベクトル
【０１４０】
【数６２】
2002288237.tif 000062
【０１４１】
を計算し、そのノルムＪ（式(7) ）が十分小さい時は上記の
【０１４２】
【数６３】
2002288237.tif 000063
【０１４３】
をシステムの一般化加速度とする。そうでない場合は以下を実行する。
（ｃ）拘束誤差の速度オフセット
【０１４４】
【数６４】
2002288237.tif 000064
【０１４５】
、拘束力修正方向ベクトル
【０１４６】
【数６５】
2002288237.tif 000065
【０１４７】
、および対応する一般化加速度
【０１４８】
【数６６】
2002288237.tif 000066
【０１４９】
を計算する。
（ｄ）上記の諸量を用いて、式（24）の係数Ｐ，Ｐ0 を計算し、式（23）を解いてゲインαを計算する。
（ｅ）すると、式（15）の右辺より、システムの一般化加速度
【０１５０】
【数６７】
2002288237.tif 000067
【０１５１】
が直ちに求まる。
対応する式を見れば明らかであるが、このアルゴリズムは主要部分が固定数の内積と一次結合のみで構成されおり、これに関与するベクトルはシステムの一般化速度
【０１５２】
【数６８】
2002288237.tif 000068
【０１５３】
もしくは拘束誤差
【０１５４】
【数６９】
2002288237.tif 000069
【０１５５】
の次元である。行列演算が式（16）（拘束前の質量行列の逆行列とベクトルの積）及び式（14）（拘束のヤコビ行列の転置とベクトルの積）に現れるが、アルゴリズムは行列そのものを必要としない。しかもこれらの演算は式（3b) などを利用するとボディまたはジョイントごとに分割して実行できる。従って計算コストがシステムの自由度（ボディ数）のみならず拘束される自由度（ジョイント数）に関しても１乗のオーダーであって効率が良く、並列計算にも適している。また最終的に扱う行列の次元が極めて小さくなり、殆どの場合行列を用いず直接インラインに式を展開できる。実際、後述する汎用プログラムは行列というデータ型を用いず、ベクトルのみを使って極めて簡潔に実装されている。
１．４  誤差ベクトルのスケーリング
式（7)で導入したスケーリング行列Ｓについては正定値かつ対角という以外何も説明していなかったが、この選択が数値的特性を大きく左右する。確かにゲインを式（23）から計算すれば、スケーリング行列Ｓが正定値である限り誤差ノルムＪは式(8）に従って全く同じ指数減衰をするが、誤差ベクトル
【０１５６】
【数７０】
2002288237.tif 000070
【０１５７】
の各成分の挙動は大きく異なる。この点を明らかにするため、ゲインの方程式（23）左辺の係数Ｐに注目する。これは現在の誤差ベクトルがどの程度誤差ノルムに影響を与えるかを表すものと考えられるが、その定義式（24a）を式（16b）と式（14）を用いて以下の様に書き直してみる。
【０１５８】
【数７１】
2002288237.tif 000071
【０１５９】
で与えられ、幾何学的拘束を選択した時点で決まってしまうものである。ここでＳを単位行列としてスケーリングを行わない最も素朴なペナルティ法を考えて見ると、行列Ｄの最大・最小固有値が非常に異なるときには、感度の高い前者の固有ベクトルに対応する方向に集中して拘束力を出そうとする。そしてこれが限界に達すると他の固有ベクトル方向に急激に移行する、という過程を繰り返すので、たとえ誤差ノルムが滑らかに指数減衰しても拘束力の向きが激しく変動する事になる。スケーリング行列Ｓは、行列Ｄの条件数を小さくして上記の不具合を防止するために必要となる。
【０１６０】
スケーリングには行列Ｄの固有値が必要となるが、これを正確に求める事は非現実的である。しかし幸い式（27）より行列Ｄは常に対称で、しかも質量行列Ｍが正定値なので非負定値（拘束が特異でなければ正定値）
である事が保証されている。従って行列Ｄの固有値はその対角成分の大きさで見積もれる。この考察から、最も簡便には、スケーリング行列Ｓを行列Ｄの両側から乗じると対角成分がすべて同じ値になる様に設定すれば良い事が分かる。即ち
【０１６１】
【数７２】  
2002288237.tif 000072
【０１６２】
但しＤｉ，ｉは行列Ｄの第１対角成分である。また式（28a ）のＳmax は0 ＜Ｄｉ，ｉなるｓi内最大のものを表すが、これで割っているのは、スケーリング行列Ｓの成分の最大値が１となる様に正規化するためである。対角成分Ｄｉ，ｉは式（27）の右辺の両側から単位ベクトル
【０１６３】
【数７３】
2002288237.tif 000073
【０１６４】
を乗じる事で得られるが、これは
【０１６５】
【数７４】
2002288237.tif 000074
【０１６６】
の様に書く事ができるため、拘束力修正方向ベクトル（式(14) )、対応する加速度（式(16b) ）及び係数Ｐ（式(24a))のサブルーチンを使って簡単に計算する事ができる。
【０１６７】
スケーリング行列Ｓは少なくとも初期時刻で計算しておく必要があるが、積分ステップごとに評価するべきものではなく、多くの問題では定数としてもかまわない。但し行列Ｄが時変でしかも条件数が大きく変動する場合には、適当な間隔（通常は1000ステップに１回程度で十分）で行列Ｓを更新するのが望ましい。但し再スケーリングを行う場合には、その前後で誤差ノルムＪの値が不連続にならない様、行列Ｓに定数倍の補正を行う。
２．不等式拘束
２．１  速度を含まない場合
前節では等式拘束のみを扱った。しかし例えばラッチ機構の様に、一般化座標が或る範囲内の間は自由で、それを越えると拘束が生じる場合も存在する。この様な機構は、一般化座標に関する不等式拘束
【０１６８】
【数７５】              
2002288237.tif 000075
【０１６９】
でモデル化する事ができる  ( 拘束誤差は簡単のためしばらくの間スカラーとする) 。この種の幾何学的拘束を前節までの定式化で直接扱う事はできないが、これと等価な等式拘束に変換ずれば可能である。変換法は一意ではないが、最も簡単なのは不等式が満たされている間は強制的にδ＝0 と置く事であるが、この場合不等式の境界で速度や加速度が不連続になる恐れが有る。この不具合を回避する方法は、上式(31)左辺の関数δを
【０１７０】
【数７６】
2002288237.tif 000076
【０１７１】
の如く不等式が満たされている時0 となる様に修正し、式(5b)のホロノミック拘束として表現する事である。但し式(32b) のεは適当な微小量である。式(32a) よりδIN( δ) は元の拘束誤差δのｎ次関数であるが、後者の大きさがεの時に前者も同じ絶対値となる。次数ｎは適当な整数であるが、δ＝0でδINの加速度まで連続とするため通常は3 以上にする。この拘束は, 0≦δの間は明らかに何もしないが、δ＜0になると拘束力が発生してδが−εより小さくならないことを保証するので、不等式(31)が維持される。但し、ラッチ機構の場合は、―度不等式の境界に達するとこれに捕捉されるので、その後は例え０≦δとなってもδ＜0の式を使ってδINを計算しなければならない。なお、不等式の境界δ＝０では衝突現象 (速度不連続) が生じ得るが、式(8) よりその時定数は1/ωであり、この過程は減衰を伴うので―般に非弾性的である点が注意される。
【０１７２】
さて、等式拘束に変換した式(32a) を適用するためには、通常のホロノミック拘束と同様に拘束誤差の速度と加速度が必要である。そこで時間微分すると、以下を得る。
【０１７３】
【数７７】
2002288237.tif 000077
【０１７４】
式（31）の両辺は簡単のためスカラーとしたが、この様な不等式拘束が複数存在する場合も各々に対して式（32a ）の変換を行なえば良い。例えばｉ番目の拘束を不等式に変換したいなら、式（32a ）に従って拘束誤差べクトル
【０１７５】
【数７８】
2002288237.tif 000078
【０１７６】
を以下の
【０１７７】
【数７９】
2002288237.tif 000079
【０１７８】
に置き換える。
【０１７９】
【数８０】
2002288237.tif 000080
【０１８０】
同様にして、拘束誤差の速度及び加速度の式（33a),(33b ）も以下の如く簡単にべクトルへ拡張することができる。
【０１８１】
【数８１】
2002288237.tif 000081
【０１８２】
【数８２】
2002288237.tif 000082
【０１８３】
中の他の式を不等式にする場合も、ＡIN、ＢIN及び
【０１８４】
【数８３】
2002288237.tif 000083
【０１８５】
の対応する成分を定数の１や０から各々ａIN、ｂIN及びｃINへ次々に置き換えて行くだけで良い。
【０１８６】
等式拘束に変換した式（32a ）を適用するには通常のホロノミック拘束と同じく式（18）の右辺に現れる諸量が必要であるが、これらは式（35a ）と式（35b）より
【０１８７】
【数８４】
2002288237.tif 000084
【０１８８】
のように直ちに得られ、そしてこれと式（14）より、対応する拘束力修正方向べクトル
【０１８９】
【数８５】
2002288237.tif 000085
【０１９０】
とすればよいことが分かる。
２．２  速度を含む場合
前節と全く同じ手法を用いて、速度を含む不等式拘束
【０１９１】
【数８６】
2002288237.tif 000086
【０１９２】
を定式化することができる。（不等号はべクトルの成分毎に適用されるものとする。）即ち、式（34a ）と同様に行列関数ＡINを用いて
【０１９３】
【数８７】
2002288237.tif 000087
【０１９４】
の如く定義し、等式拘束
【０１９５】
【数８８】
2002288237.tif 000088
【０１９６】
として扱えば良い。もちろん式（34b ）の様に行列ＡINの対応する成分を１とすれば、等式拘束を混在させることも可能である。式（38）を時間微分すると、式（35ａ）に対応する式
【０１９７】
【数８９】
2002288237.tif 000089
【０１９８】
を得るが、これより変換後の等式拘束を特徴付ける諸量（式（17）参照）は、
【０１９９】
【数９０】
2002288237.tif 000090
【０２００】
となり、これと式（14）より、対応する拘束力修正方向べクトルは
【０２０１】
【数９１】
2002288237.tif 000091
【０２０２】
と置けば良いことが分かる。
【０２０３】
上記の手法を用いると、任意の等式拘束はその一部または全部をラッチ型の不等式拘束に変換することができる。
【０２０４】
【発明の実施の形態】
本発明は、基本的にはペナルティ法に基づいているが、拘束の安定性を保証し、かつ誤差の目標値を確実に維持させるため、ゲインを可変とする。本発明の手法は、ホロノミック、非ホロノミックを問わず、あらゆる等式拘束と、不等式拘束を統一的に扱うことができる。またシステム自由度のみならず拘束される自由度に関しても計算コストがＯ (Ｎ）であって効率が良く、しかもボディ・ジョイント毎に計算を分散できるので、並列計算に適している。更に、固定数のベクトル演算だけで順ダイナミクスが解けるため、行列を使用しないで極めて簡潔なコーディングを可能にする。そして本発明の原理に基づく汎用プログラムを開発し、いくつかのシステムのシミュレーションを通じて計算結果の妥当性を検証し、有用性を確認できた。以下、本発明の実施例について説明する。  
３．数値計算例
本発明の手法に基づく汎用柔軟多体動解析プログラムを開発し、その妥当性と有効性を幾つかのシステムについて検証した。典型的な解析例を以下に示す。
なお、簡単のため、これらの解析においてパラメータωは定数、式（１１）における未定数の推定値
【０２０５】
【数９２】
2002288237.tif 000092
【０２０６】
は零とし、また誤差ノルムＪは零でなく与えられた目標値εに収束させた。基本的な特性に変わりは無いが、これらの簡略化をせずに前述のアルゴリズムを完全に実装すると、解の精度と計算時間の両面から見た総合的な成績が更に数百倍向上する事も確認済みである。
３．１剛体３リンク自由振り子の動解析
図２に示す２次元剛体システムのシミュレーションを行う。システムは線密度ρ＝１ｋｇ／ｍの剛体棒Ｂ1 、Ｂ2 、Ｂ3 と、Ｂ2 上に固定された質点ｍ＝５ｋｇより成り、全体が一様重力（−ｙ方向）のみを受けている。剛体棒の長さは
【０２０７】
【数９３】
2002288237.tif 000093
【０２０８】
とし、図中の寸法ｄは０．０５ｍと置く。また、同図の如く、剛体棒Ｂ1 、Ｂ2 、Ｂ3 の各々に沿って質量中心を原点とした座標軸ｒx1、ｒx2、ｒx3が固定されている。剛体Ｂ1 、Ｂ3 は、各々回転ジョイントＪ1 、Ｊ4 で水平なベース上に拘束されており、更に回転ジョイントＪ2 、Ｊ3 によって剛体Ｂ2 と連結されている。従ってこのシステムのトポロジは閉じており、独立な一般化座標に関して記述される通常の定式化では解析が困難な問題である。
なお、図２においてジョイント角θ1 、θ2 、θ3 を
【０２０９】
【数９４】
2002288237.tif 000094
【０２１０】
と置くと剛体棒Ｂ2 が水平となり、しかも幾何学的拘束が満足される。この静止状態をｔ＝0として、重力によって誘起される振り子運動をｔ＝１．５secまで計算した結果を図３に示す。但し式（8)における拘束の減衰係数はω＝１０3 sec−１、目標誤差はε= １０−３、数値積分の許容打ち切り誤差はεint ＝１０−５とした。なお、目標誤差εは式(6) で定義される変換された誤差ベクトル
【０２１１】
【数９５】
2002288237.tif 000095
【０２１２】
に対するものであるから、実際のホロノミック拘束誤差
【０２１３】
【数９６】
2002288237.tif 000096
【０２１４】
のノルムに換算すると、ε／ω＝１０−５となる点に注意されたい。本システムにおいて誤差ベクトル
【０２１５】
【数９７】
2002288237.tif 000097
【０２１６】
は距離［ｍ］の次元を持つので、これは即ちジョイント結合部のずれを約１０μmに抑える事に相当する。
【０２１７】
図３（ａ）は、図２のジョイント角θ1 、θ4 である。周期的な振り子運動が生じている事が分かる。一見滑らかな応答であるが、システムはｔ＝0.395sec とｔ＝1.159sec の２箇所で特異姿勢を通過している。これは、ジョイント角速度の比
【０２１８】
【数９８】
2002288237.tif 000098
【０２１９】
をプロットした図３（ｂ）が不連続になる事から判断できる。両時刻共にジョイント各はθ1 ＝296.5deg 、θ4 ＝270.0deg である。即ち、リンクＢ3 が鉛直下方を向いた時に特異姿勢となる。Baumgarte の方法では特異姿勢になると、拘束力が計算できないし、他の多くの数値計算法でも解析が困難となる。しかしながら本シュミレーションにおいては拘束の誤差ノルムＪが目標値ε＝10−３へすみやかに収束し完全に維持されている（図３ (ｃ））。しかも図３（ｄ）よりエネルギー保存則も満たされており、（ΔＴ，ΔＵ，ΔＥは各々システムの運動エネルギー、位置エネルギーおよび全エネルギーの変化を表す）、力学的に妥当な解が安定に得られている事が分かる。これは提案したアルゴリズムが特異な拘束に対しても頑健である事を示している。図３（ｅ）はペナルティ法におけるゲインαの履歴である。式（１１）や式（１３）からも分かる様に、これは全体的な拘束力の大きさに対応しているが、同図を見ると、特異姿勢の通過時に極めて大きな値が生じている事が分かる。また、このように急峻なピークを持つ力を受ける運動を追跡するには当然積分刻みを小さくする必要が有るので、特異姿勢の付近は計算コストが余計にかかる。これは、シミュレーション開始時から計測したＣＰＵ時間ＴCPU を示した図３（ｆ）を見れば明らかである。また同図から、全体の計算時間が１０分３６秒である事も分かる。但しこれにはアニメーションとデータ出力の時間も含まれており、それを除く純粋な数値計算だけならば９分３８秒で済む。何れにしろ、十分実用的な解析ができているといえる。
３．２  柔軟骨組構造物の動解析
柔軟体の運動方程式は、適当な離散化を施すと殆どの場合、式（１）の形式に帰着される。本プログラムは剛体のみならず、そのような柔軟体をすべて扱う事ができる。そこで次に、図４に示す柔軟骨組宇宙構造物のシミュレーションを行った。本システムは20枚の正６角形パネルから成っており正５角形の穴が１２箇所に有る。炭素原子が中空球殻状に結合した分子を一般にフラーレン（fullerene ）と呼ぶが、同図の構造物はその典型であるサッカーボール型炭素分子Ｃ６０と同じ形状をしている。正６角形パネルは円管部材を辺と対角線に配したラーメン構造であり、隣接するパネル同士も剛に結合されている。円管部材はすべて自然長２ｍ、外径１ｃｍ、厚さ１ｍｍであり、密度ρ＝2.80×103 ｋｇ／ｍ3 、縦弾性係数Ｅ＝7.26×1010Ｐａの７０７５Ａ１−Ｔ６アルミ合金製とする。各円管の柔軟変形は、３次のエルミート補間多項式を用いて慣性系で離散化した。この場合システムの一般化座標ベクトル（拘束前）は2520次元、拘束式の数は1560となる。
【０２２０】
この種の球殻状構造は同一表面積の平面形状より剛性を高める事が容易で、オムニ性も有るため、太陽光発電や通信・観測などに有用である。そこで、ここでは本構造物の正６角形パネルを太陽電池とし、これを用いて電気推進を行う場合を想定する。スラスタは１基で、正６角形パネルの中央に置き（図４）、推力ｆ＝０．１Ｎとする。なお、骨組み構造自身の剛性に注目するため、パネルに張る太陽電池シートの質量や柔軟性等は無視し、また推力以外の外力は働かないものとする。
【０２２１】
本システムのシミュレーションをｔ＝3.0secまで行った。但し式（8 ）における拘束の減衰係数はω＝103 sec−１、目標誤差はε= 10−１、数値積分の許容打ち切り誤差はεint ＝10−７とした。得られた結果を図５に示す。図５（ａ）のΔｈは、本構造物の進行方向の高さｈ（より正確には、スラスタの在るノードと点対称位置のノードの進行方向距離）の変化を表す。静止状態にステップ状の推力が加わる事で、構造物全体がこの方向に約15ｍｍ縮み、約1.8Ｈｚの残留振動が生じている。図５（ｂ）と図５（ｃ）は、各々システムの力学的エネルギー変化ΔＥ及び運動エネルギー変化ΔＴの履歴である。前者はスラスタの行った仕事に相当する。後者にはΔｈに見られない周波数の振動も乗っているが、これは運動エネルギーがすべての柔軟モードの影響を受けるためである。図５（ｄ）から明らかな様に、巨大な自由度を持つ本システムにおいても拘束の目標精度が安定かつ完全に達成されている。図５（ｅ）はペナルティ法におけるゲインα、図５（ｆ）は前節と同一の計算機環境において計測したＣＰＵ時間ＴCPU を示す。同図より、強い拘束力が生じている時は計算負荷が比較的大きいことが見て取れる。なお、総計算時間はアニメーションとデータ出力も含めて５６３分であった。このままでは十分実用的とは言い難いが、これは飽くまで１台のＰＣを使用した場合である。前述した様に、提案した定式化はボディまたはジョイントごとにほぼ完全な並列化が可能である。本システムは２１０個のボディ（部材）で構成されているが、これは同じＰＣでも210台に分散させると、理想的には計算時間を約２．７分に短縮できる事を示唆している。実際、上記システムの10倍の自由度を有する巨大柔軟構造物の運動を9台のＰＣで並列計算したところ、所要時間が約1／9に短縮される事が確認された。
３．３  ラッチ−バネ機構付き柔軟２リンク系の動解析
前述した不等式拘束を含むマルチボディシステムの定式化手法を、先に開発した汎用プログラムに適用し、図６に示すラッチ−バネ機構付柔軟２リンク系のシュミレーションを行なった。システムは自然長２ｍ、外径1cm、厚さ1mmの円管Ｂ1 、Ｂ2 と両者を接続するジョイントＪ1 で構成されている。円管は密度ρ＝2.80×103 ｋｇ／ｍ3 、縦弾性係数Ｅ＝7.26×1010Ｐａの７０７５Ａｌ−Ｔ６アルミ合金製とする。各円管の柔軟変形は、３次のエルミート補完多項式を用いて慣性系で離散化した。ジョイントＪ1 は、球面ジョイントにラッチ機構（ラッチ角120deg）と捻りバネ（剛性ｋ＝0.1Nm／rad、自然状態の角度90deg）を取り付けたものである。システム全体は初期時刻において、慣性空間Ｏ−ｘｙｚの原点OにジョイントＪ1 を一致させ、ｙ軸に関して対称となる形で静止しており円管Ｂ1 、Ｂ2 はｘ軸に対して0.01degの初期角度をもっている。また外力は働かないものとする。計算は３次元で行なうが、実現する運動は２次元的（ｘｙ平面内）である。
【０２２２】
本システムのシミュレーションをｔ＝1.2secまで行なった結果を図７に示す。但し拘束の減衰係数はω＝104 sec−１：目標誤差はε＝10−３、数値積分の許容打切り誤差はεint＝10−10 とした。図７（ａ）は円管Ｂ2 先端の垂直位置ｙtip （図６参照）を示す。捻りバネによって円管が閉じ始めるが、約0.685秒後にラッチがかかってジョイントＪ1 における角度が固定され、その衝撃で残留振動が生じている。図７（ｂ）はシステムの力学的エネルギー変化△Ｅの履歴である。ラッチの瞬間にＪ1 回りの角速度が失われるのでシステムエネルギーは不連続に減少するが、その前後では保存しているので、得られた解は力学的に妥当であると考えられる。また図７（ｃ）より拘束の誤差ノルムＪも目標値ε＝10−３へ速やかに収束し完全に維持されているので、拘束の精度と安定性も実証されている。図７（ｄ）はペナルティ法におけるゲインαの履歴であるが、式（１１）や式（１３）からも明らかな様に、これは全体的な拘束力の大きさに対応している。同図から、ラッチの瞬間に衝撃力が生じている事が見て取れる。しかしこれによる計算負荷の増加は殆ど無い事が、ＣＰＵ時間ＴCPU をプロットした図７（ｅ）から分かる。もしこの機構をラッチの瞬間に剛性が大きく変化する仮想的なバネ等でモデル化したならば、運動方程式がスティフとなって非現実的な計算時間を要するであろう。これは２節で提案した不等式拘束の定式化の有効性を示している。また同図より、総計算時間は９分６秒である（ただしパソコン［Power   Macintosh   Ｇ３  ３５０ＭＨｚ］を使用し、アニメーションとデータ出力も行う場合）。
【０２２３】
このように、多様なマルチボディシステムに対して本発明による定式化を行なうことにより、システムモデルをコンピュータで取り扱い易い形に構築できるため、システム制御に適用する場合も、制御に必要なデータ処理量を少なくでき、制御を高速化できる。
【０２２４】
【発明の効果】
本発明は、極めて広範囲の剛体／柔軟マルチボディシステムに適用して定式化を行なうことができる。しかも幾何学的拘束の安定性は常に保証され、かつその目標精度は確実に維持される。従ってアルゴリズムのパラメータ選択が容易であり、経験が要求されないため実用的である。またシステム自由度のみならず、拘束される自由度に関しても計算量は単純に比例するだけなので、効率が良い。またボディ・ジョイント毎に計算を分散できるので、並列計算に適している。更に、固定数のベクトル演算だけで順ダイナミクスが解けるため、行列を使用しない極めて簡潔なコーディングを可能にする。このため、システムモデルの運動予測や制御を、従来方法に比べて少ないハード量でより高速に実行することを可能にする。
＜参考論文＞
１：J.Baumgarte,
“Stabilization of Constraints and Integrals of Motion in Dynamical Systems",
Computer Methods in Applied Mechanics and Engineering,
Vol.1,pp.1-16,1972. 
２：J.W.Baumgarte,
“A New Method of Stabilization for Holonomic Constraints",
Journal of Applied Mechanics, 
Vol.50,pp.869-870,December 1983.
３：K.C.Park and J.C.Chiou, 
“Stabilization of Computational Procedures for Constrained Dynamical  Systems",
Journal of Guidance,Control,and Dynamics, 
Vol.11,No.4,pp.365-370,July-August 1988.
４：K.C.Park and J.C.Chiou and J.D.Downer,
“Explicit-Implicit Staggered Procedure for Multibody Dynamics Analys
is",
Journal of Guidance,Control,and Dynamics, 
Vol.13,No.3,pp.562-570,May-June 1990. 
５：E.Bayo and J.G.de Jalon and M.A.Serna,
“A Modified Lagrangian Formulation for the Dynamic Analysis of Constrained Mechanical Systems", 
Computer Methods in Applied Mechanics and Engineering,
Vol.71,pp.183-195,1988. 
６：A.Jain and G.Rodriguez, 
Recursive Dynamics Algorithm for Multibody Systems with Prescribed Motion,
Journal of Guidance,Control,and Dynamics, 
Vol.16,No.5,pp.830-837,September-October,1993.
７：長田  隆
マルチボディ・システムの統一的定式化法（可変ゲインペナルティ法による数値解析アルゴリズム），機論，66-649，Ｃ(2000).2919-2925
【図面の簡単な説明】
【図１】
本発明によるマルチボディシステム定式化処理の概略フローである。
【図２】
本発明を適用して動解析を行なった２次元剛体システムの説明図である。
【図３】
２次元剛体システムの動解析結果のレスポンス波形のグラフである。
【図４】
本発明を適用して動解析を行なったＣ６０フラーレン構造の説明図である。
【図５】
Ｃ６０フラーレン構造の解析結果のレスポンス波形のグラフである。
【図６】
本発明を適用してシミュレーションしたラッチ−バネ機構付柔軟２リンク系の説明図である。
【図７】
ラッチ−バネ機構付柔軟２リンク系のシミュレーション結果のレスポンス波形のグラフである。