特許 7202567 多体系の最小力学パラメータ同定装置、方法及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2022-12-28

(45)【発行日】2023-01-12

(54)【発明の名称】多体系の最小力学パラメータ同定装置、方法及びプログラム

(51)【国際特許分類】

   G01M 1/10 20060101AFI20230104BHJP

   G01M 1/12 20060101ALI20230104BHJP

【ＦＩ】

G01M1/10

G01M1/12

【請求項の数】 9

(21)【出願番号】P 2019034164

(22)【出願日】2019-02-27

(65)【公開番号】P2020139784

(43)【公開日】2020-09-03

【審査請求日】2022-02-15

(73)【特許権者】

【識別番号】304021417

【氏名又は名称】国立大学法人東京工業大学

(74)【代理人】

【識別番号】100100011

【弁理士】

【氏名又は名称】五十嵐 省三

(72)【発明者】

【氏名】大熊 政明

(72)【発明者】

【氏名】本間 貴大

【審査官】森口 正治

(56)【参考文献】

【文献】特開２０１５－８９５８５（ＪＰ，Ａ）

【文献】国際公開第２００９／１４７８７５（ＷＯ，Ａ１）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０１Ｍ １／００－９９／００

(57)【特許請求の範囲】

【請求項1】

複数の物体が接続された多体系モデルを疑似的周辺自由境界条件で搭載し、前記多体系モデルの自由振動の固有角振動数を検出して複数姿勢での全体慣性特性を計測するための全体慣性特性計測ユニットと、
前記多体系モデルを疑似的周辺自由境界条件で搭載し、前記多体系モデルをリンク間内力で相対運動させて振動又は運動を計測するための相対運動計測ユニットと、
前記全体慣性特性計測ユニット及び前記相対運動計測ユニットに接続され、前記複数姿勢での前記全体慣性特性及び前記振動又は運動を連立させて前記各物体の最小力学パラメータを同定するための制御ユニットと
を具備する多体系の最小力学パラメータ同定装置。

【請求項2】

前記制御ユニットは前記複数姿勢での全体慣性特性に基づいて前記各物体の質量に関する第１の最小力学パラメータ、前記各物体の重心に関する第２の最小力学パラメータ及び前記各物体の慣性テンソルの慣性乗積成分に関する第３の最小力学パラメータを同定するようにした請求項１に記載の多体系の最小力学パラメータ同定装置。

【請求項3】

前記制御ユニットはさらに前記複数姿勢での全体慣性特性及び前記各物体の前記相対運動に基づいて前記各物体の慣性テンソルの慣性対角項成分に関する第４の最小力学パラメータを同定するようにした請求項２に記載の多体系の最小力学パラメータ同定装置。

【請求項4】

複数の物体が接続された多体系モデルを疑似的周辺自由境界条件で該多体系モデルの自由振動の固有角振動数を検出して複数姿勢での全体慣性特性を計測するための全体慣性特性計測工程と、
前記多体系モデルを疑似的周辺自由境界条件で搭載し、前記多体系モデルをリンク間内力で相対運動させて振動又は運動を計測するための相対運動計測工程と、
前記複数姿勢での前記全体慣性特性及び前記振動又は運動を連立させて前記各物体の最小力学パラメータを同定するための最小力学パラメータ同定工程と
を具備する多体系の最小力学パラメータ同定方法。

【請求項5】

前記最小力学パラメータ同定工程は、
前記複数姿勢での前記全体慣性特性に基づき前記各物体の質量に関する第１の最小力学パラメータを同定するための第１の最小力学パラメータ同定工程と、
前記複数姿勢での前記全体慣性特性に基づき前記各物体の重心に関する第２の最小力学パラメータを同定するための第２の最小力学パラメータ同定工程と、
前記複数姿勢での前記全体慣性特性に基づき前記各物体の慣性テンソルの慣性乗積成分に関する第３の最小力学パラメータを同定するための第３の最小力学パラメータ同定工程と
を具備する請求項４に記載の多体系の最小力学パラメータ同定方法。

【請求項6】

前記最小力学パラメータ同定工程は、さらに、前記複数姿勢での全体慣性特性及び前記各物体の前記相対運動に基づき前記各物体の慣性テンソルの慣性対角項成分に関する第４の最小力学パラメータを同定するための第４の最小力学パラメータ同定工程を具備する請求項５に記載の多体系の最小力学パラメータ同定方法。

【請求項7】

複数の物体が接続された多体系モデルを疑似的周辺自由境界条件で該多体系モデルの自由振動の固有角振動数を検出して複数姿勢での全体慣性特性を計測するための全体慣性特性計測手順と、
前記多体系モデルを疑似的周辺自由境界条件で搭載し、前記多体系モデルをリンク間内力で相対運動させて振動又は運動を計測するための相対運動計測手順を具備し、
前記複数姿勢での前記全体慣性特性及び前記振動又は運動を連立させて前記各物体の最小力学パラメータを同定するための最小力学パラメータ同定手順と
を具備する多体系の最小力学パラメータ同定プログラム。

【請求項8】

前記最小力学パラメータ同定手順は、
前記複数姿勢での前記全体慣性特性に基づき前記各物体の質量に関する第１の最小力学パラメータを同定するための第１の最小力学パラメータ同定手順と、
前記複数姿勢での前記全体慣性特性に基づき前記各物体の重心に関する第２の最小力学パラメータを同定するための第２の最小力学パラメータ同定手順と、
前記複数姿勢での前記全体慣性特性に基づき前記各物体の慣性テンソルの慣性乗積成分に関する第３の最小力学パラメータを同定するための第３の最小力学パラメータ同定手順と
を具備する請求項７に記載の多体系の最小力学パラメータ同定プログラム。

【請求項9】

前記最小力学パラメータ同定手順は、さらに、前記複数姿勢での全体慣性特性及び前記各物体の前記相対運動に基づき前記各物体の慣性テンソルの慣性対角項成分に関する第４の最小力学パラメータを同定するための第４の最小力学パラメータ同定手順を具備する請求項８に記載の多体系の最小力学パラメータ同定プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は多体系の最小力学パラメータ同定装置、方法及びプログラムに関する。

【背景技術】

【0002】

物体には動的な挙動の解析、評価に必要な１個のパラメータよりなる質量、３個のパラメータよりなる重心及び６個のパラメータよりなる慣性テンソルの合計１０個のパラメータよりなる固有の慣性物性が定義されている。尚、慣性テンソルは慣性モーメント（慣性対角項）及び慣性乗積よりなる。たとえば、重心の位置は自動車等の乗り物や、重機等において横転防止等につながる重要なパラメータである。また、自動車等の乗り物においては、全般的に，操縦のし易さ、快適性、運動性能等の評価に慣性対角項のみならず慣性乗積が大きく影響する。

【0003】

自由振動計測による従来の単体モデルの慣性特性同定ユニットは、図１に示すように、単体モデル１００を搭載するためのプラットフォーム１０１と、プラットフォーム１０１を自由振動可能に疑似的に支持する（疑似的周辺自由境界条件での支持という）支持部材（ばね／ワイヤ）１０２と、プラットフォーム１０１の振動を検知するためのセンサ１０３（たとえば加速度センサ）とよりなる（参照：特許文献１、２）。この場合、プラットフォーム１０１と単体モデル１００とは一体の剛体とみなされ、加振器、インパルスハンマ等で単体モデル１００とプラットフォーム１０１の一体としての剛体を加振し、弾性モードが生じない低周波域での周波数応答を得る。予め、プラットフォーム１０１の慣性特性を精度よく同定しておき、プラットフォーム１０１と単体モデル１００との一体の剛体特性からプラットフォーム１０１の分を引くことにより単体モデル１００の慣性特性を得ることができる。このようにして、ばね、ワイヤ等の支持部材１０２の影響をモデル化して取除き高精度の同定を可能にすると共に、自由振動計測からの周波数スペクトル分析による並進、回転に基づく６個の固有角振動数の値のみの把握で同定可能であるのでセンサの感度校正を必要とせず、計測誤差も小さい。尚、この単体モデルの慣性特性同定ユニットの発展型を本発明に係る多体系モデルの全体慣性特性計測ユニットとして用いる。

【0004】

他方、ロボット工学、人間工学等においては、複数の構成物体が対偶的に接続された多体系の動力学モデルを作成するために、構成物体（リンク）長等の幾何パラメータ及び各構成物体（リンク）の慣性特性が要求される。この場合、各構成物体のＣＡＤデータがある場合には、各慣性特性はＣＡＤモデルによる計算で一応は求めることができる。しかし、組立てられた状態での機械機構（たとえばロボット）や、当然切離すことができない状態での身体においては、各慣性特性の同定は実験的に計測するしかない。その際、幾何パラメータは静的な計測から容易に同定可能であるが、各構成物体（リンク）の慣性特性は多体系の動力学モデルに対して冗長であるので、構成物体（リンク）の運動と関節トルクまたは外力データからは同定できない。そこで、多体系の動力学モデルを表現するのに必要最小限の力学パラメータを最小力学パラメータと定義している。

【0005】

第１の従来の多体系の最小力学パラメータ同定方法は、複数のリンクを関節によって接続した多体系において、各関節にトルクセンサを設け、関節トルクに対する運動を計測することによって各構成物体（リンク）の慣性特性の最小力学パラメータを同定する（参照：非特許文献１）。

【0006】

第２の従来の多体系の最小力学パラメータ同定方法は、運動とその運動に対する床反力を計測することによって各構成物体（リンク）の最小力学パラメータを同定する（参照：非特許文献２）。

【先行技術文献】

【特許文献】

【0007】

【文献】特開２０１２－１８０９２号公報

【文献】特表２０１４－５１５４９２号公報

【非特許文献】

【0008】

【文献】H. Kawasaki : ロボットアームのパラメータ同定, 計測と制御, 28 巻4 号p. 344-350, 1989.

【文献】K. Ayusawa, G. Venture, Y. Nakamura : Minimal-set of inertial parameters identification of legged robots using base-link dynamics, 日本ロボット学会誌, Vol. 27 No. 9, 2009.

【文献】R. Klopper : A measurement system for rigid body properties enabled by gravity-dependent suspension modeling, PhD.Thesis, 2009.

【文献】R. Klopper, M. Okuma : Experimental identification of rigid body inertia properties using low-frequency unbalance excitation, Proc. 50th ASME Structural dynamics conference, 2007.

【文献】R. Klopper, H. Akita, M. Okuma, S. Terada : An experimental identification method for rigid body properties enabled by gravity-dependent suspension modelling, Proceedings of joint inter-national conference on multibody system dynamics, 2010.

【発明の概要】

【発明が解決しようとする課題】

【0009】

しかしながら、上述の第１の従来の多体系の最小力学パラメータ同定方法は、トルクという力を高精度に計測する必要があるので、適用の制約が大きい。たとえば、スポーツ科学や工学における人の動きの解析、ヒューマノイドロボットの解析においては、関節すべてにトルクセンサを配置することが困難であり、また、人間を対象とする場合には関節トルクを計測することが困難であるので、適用の制約が大きいという課題がある。

【0010】

また、上述の第２の従来の多体系の最小力学パラメータ同定方法においても、床反力という力を計測する必要があるので、適用の制約が大きい。つまり、床反力としての抗力及び摩擦力の両方を高精度に計測することは困難であり、かつモーションキャプチャを用いた３次元運動の計測での最適な運動の選定が困難なため，多体系の３次元の慣性特性の同定への一般化は困難であるので、適用の制約が大きいという課題がある。

【課題を解決するための手段】

【0011】

上述の課題を解決するために、本発明に係る多体系の最小力学パラメータ同定装置は、複数の物体が接続された多体系モデルを疑似的周辺自由境界条件で搭載し、多体系モデルを搭載した系の自由振動の固有角振動数を検出する全体慣性特性計測ユニットと、その全体慣性特性計測ユニットに接続され、多体系モデルを疑似的周辺自由境界条件で搭載し，多体系モデルをリンク間内力で相対運動させてその振動または運動を計測するための相対運動計測ユニットと，全体慣性特性計測ユニット及び相対運動計測ユニットに接続され，複数姿勢での多体系モデルの全体慣性特性及びその振動または運動に基づく方程式を連立させて各物体の最小力学パラメータを同定するための制御ユニットとを具備するものである。なお，全体慣性特性計測ユニットと相対運動計測ユニットはハードウエアとしては同一にすることもできるし，別々のユニットとして準備することもできる。

【0012】

また、本発明に係る多体系の最小力学パラメータ同定方法は、複数の物体が接続された多体系モデルを疑似的周辺自由境界条件での全体慣性特性計測ユニットに搭載して複数姿勢での自由振動の固有角振動数を計測するための全体慣性特性計測工程と、多体系モデルを疑似的周辺自由境界条件で搭載してリンク間内力で相対運動させてその振動または運動を計測するための相対運動計測工程と、複数姿勢での多体系モデルの全体慣性特性及びその振動または運動に基づく方程式を連立させて各物体の最小力学パラメータを同定するための最小力学パラメータ同定工程とを具備するものである。

【0013】

さらに、本発明に係る多体系の最小力学パラメータ同定プログラムは、複数の物体が接続された多体系モデルを疑似的周辺自由境界条件で該多体系モデルの固有角振動数を検出して複数姿勢での全体慣性特性を計測するための全体慣性特性計測手順と、多体系モデルを疑似的周辺自由境界条件で搭載してリンク間内力で相対運動させてその振動又は運動を計測するための相対運動計測手順と、全体慣性特性及びその振動又は運動に基づく方程式を連立させて各物体の最小力学パラメータを同定するための最小力学パラメータ同定手順とを具備するものである。

【発明の効果】

【0014】

本発明によれば、力の計測を行うことなく、多体系の最小力学パラメータを同定できる。

【図面の簡単な説明】

【0015】

【図1】従来の単体モデルの慣性特性同定ユニットの例を示し、（Ａ）は正面図、（Ｂ）は斜視図である。

【図2】本発明に係る多体系の最小力学パラメータ固定装置の実施の形態を示すブロック図である。

【図3】図２の多体系モデルの例を示す図である。

【図4】図２の全体慣性特性計測ユニットの市販製品の例を示し、（Ａ）は正面図、（Ｂ）は斜視図である。

【図5】図２の相対運動計測ユニットの例を示し、（Ａ）は正面図、（Ｂ）は斜視図である。

【図6】図２の制御ユニットの動作を説明するためのフローチャートである。

【図7】図６の姿勢変更ステップを説明するための多体系モデルを示す斜視図である。

【図8】図６の姿勢変更ステップにおける多体系モデルの複数姿勢たとえば３通りの姿勢を示す斜視図である。

【図9】図６の全体慣性特性計測ステップ６０２の詳細なフローチャートである。

【図10】図６の相対運動計測ステップ６０４を説明する図であって、（Ａ）は計測状態を示す斜視図、（Ｂ）は計測データを図式的に示したものである。

【図11】図８の各姿勢で全体慣性特性計測結果（７回計測しての平均値）を示す表である。

【図12】図１１の計測結果に基づく重心位置に関する最小力学パラメータの同定結果およびその精度を示す図である。

【図13】図１０（Ｂ）の計測データと図１１の計測結果に基づくベースリンクの慣性テンソルに関する最小力学パラメータの同定結果およびその精度を示す図である。

【図14】図１０（Ｂ）の計測データと図１１の計測結果に基づくリンク１の慣性テンソルに関する最小力学パラメータの同定結果の精度評価を説明するための図である。

【図15】本発明の他の実施の形態を説明するための図である。

【図16】本発明が適用される分野の例を説明する図である。

【発明を実施するための形態】

【0016】

尚、以下の数１式等におけるベクトル、行列の太字、斜体の数字、文字は明細書中においては通常の数字、文字とした。

【0017】

図２は本発明に係る多体系の最小力学パラメータ同定装置を示すブロック図である。

【0018】

図２においては、多体系モデル１０の複数姿勢での全体慣性特性を計測するための全体慣性特性計測ユニット２０及び多体系モデル１０の各リンク間の相対運動を計測するための相対運動計測ユニット３０を設け、さらに、全体慣性特性計測ユニット２０及び相対運動計測ユニット３０を制御するための制御ユニット４０（たとえばマイクロコンピュータ）を設ける。なお，ハードウエア構成としては，全体慣性特性計測ユニット２０と相対運動計測ユニット３０を統一化して設計することもできる．

【0019】

図３は図２の多体系モデル１０の例を示す図である。

【0020】

図３は、Ｎ個のリンクｊよりなる木構造をなし、リンク０はベース（基底）リンクであるが、どのリンクも基底リンクになり得る。尚、人間の体幹と四肢の構造モデルの場合、体幹リンクをベースリンクとする。

【0021】

図４は図２の全体慣性特性計測ユニット２０の一例を示し、（Ａ）は正面図、（Ｂ）は斜視図である。図４の全体慣性特性計測ユニット２０は図１の従来の単体モデルの慣性特性同定ユニットを発展させたものである。

【0022】

図４においては、多体系モデル１０を搭載するためのプラットフォーム２１、プラットフォーム２１を自由振動可能に、つまり、疑似周辺自由境界条件で支持する支持部材（ばね）２２、及びプラットフォーム２１の振動を検知するためのひずみセンサ２３を設ける。ひずみセンサ２３は制御ユニット４０に接続される。この場合も、プラットフォーム２１と多体系モデル１０とは一体の剛体とみなされ、加振器、インパルスハンマ等でその一体剛体を加振し、弾性モードが生じない低周波域での周波数応答を得る。予め、プラットフォーム２１の慣性特性を精度よく同定しておき、プラットフォーム２１と多体系モデル１０との一体の慣性特性からプラットフォーム２１の分を引くことにより多体系モデル１０の慣性特性を得ることができる。このようにして、ばね、ワイヤ等の支持部材２２の影響をモデル化して取除き高精度の同定を可能にすると共に、自由振動計測からの周波数スペクトル分析による並進、回転に基づく６個の固有角振動数の値のみで同定可能であるので、多体系モデル１０の慣性特性の同定結果の誤差も小さい。

【0023】

図５は図２の相対運動計測ユニット３０の一例を示し、（Ａ）は正面図、（Ｂ）は斜視図である。

【0024】

図５においては、多体系モデル１０を搭載するためのプラットフォーム３１、プラットフォーム３１を水平方向の微小変位運動に関して疑似周辺自由境界条件で支持する支持部材（ワイヤ）３２、及び各リンクの相対運動を検知するためのモーションキャプチャ３３を設ける。モーションキャプチャ３３は制御ユニット４０に接続される。

【0025】

図３の多体系モデル１０の運動方程式は数１式のごとくなる。

【数1】

ここで、ｉ＝Ｏ、Ｃとおいて、
Ｈ_i,j（ｊ＝１、２）は慣性行列
ｑ_ｏはベースリンクの一般化座標（６自由度）
ｑ_Ｃは関節角度ベクトル
ｂ_ｉコリオリ力、遠心力、重力を表すベクトル項
τは関節トルクを表すベクトル
Ｎ_Ｃは多体系に働く外力の総数
Ｆ_ｋは多体系に働くｋ番目の外力ベクトル
Ｊ_ｉｋはｋ番目の外力を一般化力へ変換する行列

【0026】

数１式の運動方程式は、数２式に示すごとく、変形でき、質量、重心、慣性テンソルの力学パラメータに対して線形な関係となる。

【数2】

φは数３式に示すＮ個のリンクｊの力学パラメータを並べたベクトルである。

【数3】

φ_jは数４式に示すリンクｊの１０個の力学パラメータによって構成されるベクトルである。

【数4】

ここで、ｍ_ｊはリンクｊの質量
Ｉ_ｊ，ｘｘ，Ｉ_ｊ，ｙｙ，Ｉ_ｊ，ｚｚ，Ｉ_ｊ，ｘｙ，Ｉ_ｊ，ｘｙ，Ｉ_ｊ，ｙｚはリンクｊの座標原点回りの慣性テンソルＩ_ｊの各成分
ｇ_ｊ，ｘ，ｇ_ｊ，ｙ，ｇ_ｊ，ｚはリンクｊの座標原点からの重心位置ベクトルｇ_ｊの各成分
Ｙｉはφ_ｊにかかるリグレッサ行列と呼ばれ、ｑ_ｊ、ｑ_ｉの１階時間微分及びｑ_ｉの２階時間微分に基づいて構成される関数行列

【0027】

数２式からは力学パラメータのすべてを同定できない。各リンクｊの慣性特性は多体系モデル１１に対して冗長であるためである。そこで、最小力学パラメータからなるベクトルφ_Ｂを導入することにより数２式は数５式に変形できる。

【数5】

ここで、Ｙ_ｉＢは最小力学パラメータベクトルφ_Ｂのリグレッサ行列

【0028】

本発明においては、数５式の上段のトルクτを含まない運動方程式を用い、リンク長、対偶の位置等の幾何パラメータ、各リンクｊの運動と計測ユニットのばね、ワイヤの計測台取り付け点の位置と運動を計測することで得られるデータを代入して最小力学パラメータを求める方程式を作成する。そのため、多体系モデルの関節トルクを計測する必要がない。さらに、ばねやワイヤに発生する力は、ばね、ワイヤの力学的パラメータ（ばね定数や長さ）を予め精密にモデル化しているので、位置の計測から計算でき、床反力といった外力の計測を必要としない。この方程式の解法に、リンクｊ間の相対運動をさせたときの基本姿勢状態で多体系モデル１０全体の慣性特性の同定結果を組合せる。

【0029】

次に、図２の制御ユニット４０の動作を図６を参照して説明する。

【0030】

始めに、ステップ６０１にて、多体系モデル１０を全体慣性特性計測ユニット２０のプラットフォーム２１に搭載し、多体系モデル１０の姿勢を変更する。姿勢は変更後に固定される。尚、多体系モデル１０を、図７に示すごとく、２リンク系とすれば、図８の（Ａ）、（Ｂ）、（Ｃ）に示す姿勢を想定する。

【0031】

次に、ステップ６０２にて、固定された姿勢の多体系モデル１０の単体としての全体慣性特性を計測する。

【0032】

次に、ステップ６０３にて、姿勢変更での全体慣性特性計測が終了か否かを判別する。つまり、図８の（Ａ）、（Ｂ）、（Ｃ）に示す姿勢での計測がすべて終了まで、ステップ６０１、６０２を繰返し、姿勢変更が終了であれば、ステップ６０４に進む。

【0033】

ステップ６０４では、ステップ６０２にて得られた複数姿勢での全体慣性特性に基づいて各リンクの質量に関する最小力学パラメータ、各リンクの重心に関する最小力学パラメータ及び各リンクの慣性テンソルの慣性乗積成分に関する最小力学パラメータを同定する。

【0034】

次に、ステップ６０５では、多体系モデル１０を図５の相対運動計測ユニット３０のプラットフォーム３１に搭載して内力によるリンクｊ間の相対運動計測を行う。

【0035】

次に、ステップ６０６では、ステップ６０２にて得られた複数姿勢での全体慣性特性及びステップ６０５にて得られたリンク間の相対運動計測データに基づいて構成される方程式（数５式の上段部分）より各リンクｊの慣性テンソルの対角項成分に関する最小力学パラメータを同定する。

【0036】

そして、ステップ６０７に図６のルーチンは終了する。

【0037】

図９は図６の全体慣性特性計測ステップ６０２の詳細なルーチンである。

【0038】

始めに、ステップ９０１にて、プラットフォーム２１及び多体系モデル１０の一体での単体を加振してその自由振動応答から固有角振動数を複数求めておく。なお、自由振動のさせ方で観測できる固有角振動数の数は異なる場合があるので、できるだけ６個またはそれにできるだけ近い多くの数が観測できるように振動させることが加振時の注意点である。

【0039】

次に、ステップ９０２にて、未知数である慣性特性の初期値を決定する。

【0040】

次に、ステップ９０３にて、プラットフォーム２１及び多体系モデル１０の一体での単体としての質量行列Ｍをその慣性特性値から作成する。尚、質量行列は質量、重心及び慣性テンソルによって構成される６行×６列の行列である。

【0041】

次に、ステップ９０４にて、剛性行列Ｋを作成する。固定支持位置とプラットフォーム２１間に取り付けられたばねの両端位置とばね定数から剛性行列Ｋを作成することができる。尚、剛性行列Ｋは可動部（多体系モデル１０及びプラットフォーム２１）が静的釣り合い状態から微小の並進及び回転変位したときの可動部に加わる力を表す６行×６列の行列である。

【0042】

ステップ９０５では、プラットフォーム２１及び多体系モデル１０の一体での単体としての静的釣り合い状態を数値解析法により計算する。尚、静的釣り合い状態はステップ９０４で求めた剛性行列Ｋとステップ９０３で決定されている単体の質量行列Ｍおよび重力との静的つり合い方程式に基づいて解くことができる。

【0043】

次に、ステップ９０６にて、ステップ９０３にて得られた質量行列Ｍ及びステップ９０５にて得られた剛性行列Ｋからステップ９０５で求めた静的つり合い位置における微小自由振動の理論的な固有値解析の方程式を構築して解くことで、その固有値の平方根を６個の理論固有角振動数ωとして計算する。つまり、
det(K -ω²M)=0
から理論的固有角振動数ωを計算する。但し、Iは単位行列である。

【0044】

次に、ステップ９０７にて、ステップ９０１にて得られた実験固有角振動数とステップ９０６にて得られた理論固有角振動数との差の自乗ノルムを計算する。

【0045】

次に、ステップ９０８では、ステップ９０７にて得られた自乗ノルムがステップ９０３から９０８までの反復計算の前回の値より小さくなっているかどうかを判定し、小さくなっており、まださらに小さくなる収束途中であると判断された場合には、収束方向への感度解析を行い、それに基づいて慣性特性を適切に少し変更してステップ９０３～９０７を繰り返す。自乗ノルムが十分小さく収束したと判断した場合には、その時点での質量行列の値から全体慣性特性を同定する。

【0046】

そして、ステップ９０９にて図９のルーチンは終了する。

【0047】

図９のステップ９０５においては、剛性行列Ｋは釣り合い状態からの微小な剛体運動におけるばねやワイヤの影響をモデル化することによって導出する。ばねやワイヤによる影響要因は以下の３つと考えることができる。
1. ばねの復元力（ばねの軸線方向の復元力）
2. 振り子運動のように，釣り合い位置からずれると重力によって元に戻ろうとする幾何剛性
3. 釣り合い位置からモーメントの均衡が崩れた場合の復元モーメント（ばねの軸線に垂直な方向の復元力）
質量行列Ｍや剛性行列Ｋの詳しい導出は非特許文献３、４、５を参照されたい。このように、あらかじめ精密に特性値を把握したばねやワイヤを用いてその影響を精密にモデル化することで高精度な同定を可能にした。この同定法を用い複数姿勢で全体慣性特性を同定すると、図６のステップ６０４において、各リンクと全体の重心位置に関する数６式と慣性テンソルに関する数７式を計測姿勢に対応させて複数得ることができる。

【数6】

【数7】

ここで、下添字l（エル）はl（エル）回目の計測を示す
m_s、g_s、I_sは多体系全体の質量、ベースリンク座標原点からの重心位置ベクトルおよび慣性テンソル
p_jはベースリンク原点からリンクj座標原点へ向かうベクトル（j≧1）
R_jはベースリンク原点からリンクj座標原点への姿勢行列（j≧1）
[ω×] は数８式のような演算子

【数8】

【0048】

数６式、数７式から重心位置に関する最小力学パラメータ、慣性テンソルの慣性乗積に関する対応する部分に関する最小力学パラメータ、慣性テンソルの対角項成分の最小力学パラメータに関する関係式を得ることができる。たとえば図７に示す球対偶を持つ２リンク系でベースリンクの座標原点を球対偶の中心に設定した多体系においては、同定したい最小力学パラメータは、数９式、各リンクと全体の重心位置に関する式は、数１０式、各リンクと全体の慣性テンソルに関する式は、数１１式で表される。

【数9】

【数10】

【数11】

【0049】

数１０式の連立方程式を解くことにより，m_jg_jを得ることができ、数１１式の連立方程式と回転行列の性質から、I_j,xy，I_j,xz，I_j,yz，I_j,xx - I_j,yy，I_j,yy - I_j,zz，I_0,xx + I_1,xx，I_0,yy+ I_1,yy，I_0,zz + I_1,zz を得ることができる。そのため、図６のステップ６０５におけるある平面でのリンクの相対運動といった比較的単純な動作から、図６のステップ６０６にて残りの慣性テンソルに関する最小力学パラメータを得ることができ、すべての最小力学パラメータが同定可能である。

【0050】

次に、図６の内力によるリンク間の相対運動計測ステップ６０５について詳細に説明する。

【0051】

残りの慣性テンソルの対角項成分に関する最小力学パラメータを得るために、図１０に示すごとく、計測対象の多体系リンクにばねやワイヤを接続し、その反発力を内力として使ってリンク間に相対運動を起こさせ、その運動の計測を行う。モーションキャプチャ３３（図５）で計測を行う場合、あらかじめ各リンク上に取り付けた各マーカの位置からベースリンク０の並進ベクトル及び回転行列と各リンクの回転行列の時刻歴変化を求め、その並進ベクトル及び回転行列の時刻歴変化から速度、加速度、角速度、角加速度ベクトルを求める。並進ベクトル及び回転行列の各成分の時刻歴変化を各成分に対して最小二乗法的に最適な近似となる多項式に近似し、その多項式を微分することで速度、加速度、角速度、角加速度ベクトルの時刻歴変化を求め数値解析の安定性と結果の精度を高める。この計測法ではベースリンク０に疑似的周辺自由の境界条件をつくるために相対運動計測ユニットに設置されている計測系のばねやワイヤが発生する反力が外力となるため、剛性行列Ｋを用いて数５式のトルクの含まれない上段の式は数１２式のように表すことができる。

【数12】

【0052】

数１２式は並進3自由度、回転3自由度の計6自由度分の式となる。並進に関する式は全体の質量及び重心位置に関する最小力学パラメータから構成される。従って、すでに全体慣性特性計測から全体の質量と重心位置に関する最小力学パラメータは同定できているので、必要としない。数１２式の回転に関する式の内１つ以上の式と全体慣性特性同定で求めた最小力学パラメータの一部と関係式を用いることで慣性テンソルの慣性対角項成分に関する最小力学パラメータが得られ、すべての最小力学パラメータを同定することができる。たとえば、図７の場合は、数１２式は数１３式となる。

【数13】

ただし、

【数14】

【数15】

ω_jは絶対座標系でのリンクjの角速度ベクトル
R₀は絶対座標原点からベースリンク座標原点への姿勢行列

【0053】

このようにして、数１３式、数１４式、数１５式に全体慣性特性計測での数１０式、数１１式から得た、m_jg_j，I_j,xy，I_j,xz，I_j,yz，I_j,xx － I_j,yy，I_j,yy － I_j,zz，I_0,xx + I_1,xx，I_0,yy + I_1,yy，I_0,zz + I_1,zz を代入することで、慣性テンソルの対角項成分に関する最小力学パラメータを得ることができる。

【0054】

図１１は図８の各姿勢Ａ、Ｂ、Ｃで７回計測し平均化して得られた全体慣性特性結果と回転行列の成分を示す。図１０（Ｂ）の計測データと図１１の計測結果に基づいて、図１２はベースリンク０及びリンク１の重心位置に関する最小力学パラメータの実験同定結果を示し、図１３はベースリンク０の慣性テンソルに関する実験同定結果を示し、図１４はリンク１の慣性テンソルに関する実験同定結果を示す。

【0055】

図１２に示すごとく、質量と重心位置で構成されるパラメータm_jg_jの実験同定結果の最大誤差成分は±０．００２８ｋｇｍであり、ベースリンク０の単位質量当り±０．０００４ｍである。これはシミュレーションで理論的に計測データに含まれると考えられるランダムノイズレベルを考慮して１万回同定したときの平均二乗誤差平方根（ＲＭＳＥ）で規定される誤差範囲±ＲＭＳＥ内である。

【0056】

図１３、図１４に示すごとく、慣性テンソルの実験同定結果の誤差は±３．１％であり、前記と同様に、これはシミュレーションで１万回同定したときの平均二乗誤差平方根ＲＭＳＥで規定される誤差範囲±ＲＭＳＥ内である。

【0057】

尚、図６、図９に示すフローチャートはプログラムとして制御ユニット４０のメモリに格納される。

【0058】

また、上述の実施の形態においては、２リンクよりなる多体系を例としたが、本発明は３以上のリンクよりなる多体系にも適用し得る。たとえば、多体系として人間個人の最小力学パラメータ同定の場合、図１５の（Ａ）に示すごとく、複数姿勢に基づく全体慣性特性計測を行い、図１５の（Ｂ）に示すごとく、各リンク間の相対運動計測を行う。

【0059】

さらに、スポーツ工学の分野において、今まで人間個人の最小力学パラメータが同定できなかったが、本発明により人間個人の最小力学パラメータが同定可能となれば、図１６の（Ａ）に示すようなスポーツ選手ごとの関節への負担（関節トルクなど）を算出することができ、怪我を事前に防ぐための評価値を構築することが可能となる。また、その最小力学パラメータを持つ人はどのような動きをした方がよいというようにその人にフィットした運動の提案が可能となる。

【0060】

さらにまた、自動車産業においては、ドライバーひとりひとりに合わせた自動車設計が期待されている。そのためには、人種別、成人男性、女性、年齢層別の子供などの区分、または究極としては個々人に近い図１６の（Ｂ）に示すダミー人形を使用して試験を行う有用性が指摘されている。生きている人間個人の最小力学パラメータが要求されるが今までは同定できないとされていた。本発明を適応することで個々人の最小力学パラメータが同定でき、対応したダミー人形の開発などを可能にする。

【0061】

さらに、自動車工学においては、自動車全体の慣性特性だけでなく、自動車を分解せずに、ボディーの慣性特性やエンジン（またはパワートレイン（Power-train））の慣性特性をも同定できれば理想的である。図１６の（Ｃ）に示すごとく、このような弾性結合で複数のコンポーネントが結合した機械構造物を「多体構造系」とし、本発明は、この種の対象物の主要コンポーネントの剛体特性を非分解で同定可能とする。

【0062】

また、本発明は上述の実施の形態の自明の範囲のいかなる変更も適用し得る。

【産業上の利用可能性】

【0063】

本発明は、上述のごとく、スポーツ工学、自動車工学等の分野に利用できる。

【符号の説明】

【0064】

１０：多体系モデル
０：ベースリンク
１、２、…：リンク
２０：全体慣性特性計測ユニット
２１：プラットフォーム
２２：ばね
２３：ひずみケージ
３０：相対運動計測ユニット
３１：プラットフォーム
３２：ワイヤ
３３：モーションキャプチャ
４０：制御ユニット
１００：単体モデル
１０１：プラットフォーム
１０２：支持部材
１０３：加速度センサ

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】