"加速または減速運動が行われる運動体において、加速または減速運動に伴って運動体に生じるn自由度nは自然数かつ時変の振動を抑制する振動抑制方法において、運動体に付加された第1の入力によって運動体に生じる第1の応答に、第1の応答の変位が0となる時刻において、第2の入力を運動体に付加して生じる第2の応答を重ね合わせて、互いに打ち消し合うインプットシェイピング法を適用する場合に、n自由度かつ時変の振動における複数のモード間のエネルギの移り変わりのうち、モード間の相互作用によるものを無視し、固有振動数および減衰比の変化によるものを考慮して、第1の入力に対する第2の入力の大きさを決定する、振動抑制方法。加速または減速運動が行われる運動体において、加速または減速運動に伴って運動体に生じるn自由度nは自然数かつ時変の振動を抑制する振動抑制方法において、運動体に付加された第1の入力によって運動体に生じる第1の応答に、第1の応答の変位が0となる時刻において、第2の入力を運動体に付加して生じる第2の応答を重ね合わせて、互いに打ち消し合うインプットシェイピング法を適用する場合に、n自由度かつ時変の振動における複数のモード間の相互作用によるエネルギの移り変わりを無視し、加速または減速運動を行っている間に運動体自体の姿勢が変化することにより生じるコリオリがした仕事を考慮して、第1の入力に対する第2の入力の大きさを決定する、振動抑制方法。加速または減速運動が行われる運動体において、加速または減速運動に伴って運動体に生じるn自由度nは自然数かつ時変の振動を抑制する振動抑制方法において、運動体に付加された第1の入力によって運動体に生じる第1の応答に、第1の応答の変位が0となる時刻において、第2の入力を運動体に付加して生じる第2の応答を重ね合わせて、互いに打ち消し合うインプットシェイピング法を適用する場合に、n自由度かつ時変の振動におけるモード間の相互作用によるエネルギの移り変わりをゼロと見なして、固有振動数および減衰比の時間変化によるモード間のエネルギの移り変わりを考慮して、第1の入力に対する第2の入力の大きさを決定する、振動抑制方法。各振動モードにおいて、第1の入力に対応する入力を打ち消す第2の入力に対応する入力の大きさを個別に算出して、振動モードごとに算出された入力を合成することで第2の入力の大きさを決定する、請求項1から3のいずれか1つに記載の振動抑制方法。運動体が多関節型ロボットアームであって、ロボットアームに生じるn自由度かつ時変の振動を抑制する、請求項1から4のいずれか1つに記載の振動抑制方法。ロボットアームは、様々な作業に対応可能とするために、一般に多関節および多自由度を有している。そのため、ロボットアームにおいて、一の関節が駆動された運動動作を行っている間に他の関節が駆動されると、運動動作を行っているアーム自体の姿勢が変化する。そのため、アームの固有振動数が時間によって変化することになる。このような場合にあっては時変の振動系となるため、時不変の振動系を対象として従来のインプットシェイピング法を単に適用するだけでは、振動抑制に対応することができない。このような時変の振動系において振動を抑制する方法としては、いくつかの提案がなされている例えば、特許文献1および非特許文献2参照。しかしながら、時変の振動抑制に対して、インプットシェイピング法を効果的に適用できる振動抑制方法が求められている。特に、自由度が2以上のロボットアームに対しては、複数の固有振動数が時々刻々と変化していくため、インプットシェイピング法を用いて多自由度かつ時変の振動を効果的に抑制することが求められている。従って、本発明の目的は、上記課題を解決することにあって、加速または減速運動に伴って運動体に生じるn自由度かつ時変の振動を、インプットシェイピング法を効果的に適用して抑制する振動抑制方法を提供することにある。上記目的を達成するために、本発明は以下のように構成する。本発明の第1態様によれば、加速または減速運動が行われる運動体において、加速または減速運動に伴って運動体に生じるn自由度nは自然数かつ時変の振動を抑制する振動抑制方法において、運動体に付加された第1の入力によって運動体に生じる第1の応答に、第1の応答の変位が0となる時刻において、第2の入力を運動体に付加して生じる第2の応答を重ね合わせて、互いに打ち消し合うインプットシェイピング法を適用する場合に、n自由度かつ時変の振動における複数のモード間のエネルギの移り変わりのうち、モード間の相互作用によるものを無視し、固有振動数および減衰比の変化によるものを考慮して、第1の入力に対する第2の入力の大きさを決定する、振動抑制方法を提供する。本発明の第2態様によれば、加速または減速運動が行われる運動体において、加速または減速運動に伴って運動体に生じるn自由度nは自然数かつ時変の振動を抑制する振動抑制方法において、運動体に付加された第1の入力によって運動体に生じる第1の応答に、第1の応答の変位が0となる時刻において、第2の入力を運動体に付加して生じる第2の応答を重ね合わせて、互いに打ち消し合うインプットシェイピング法を適用する場合に、n自由度かつ時変の振動における複数のモード間の相互作用によるエネルギの移り変わりを無視し、加速または減速運動を行っている間に運動体自体の姿勢が変化することにより生じるコリオリがした仕事を考慮して、第1の入力に対する第2の入力の大きさを決定する、振動抑制方法を提供する。本発明の第3態様によれば、加速または減速運動が行われる運動体において、加速または減速運動に伴って運動体に生じるn自由度nは自然数かつ時変の振動を抑制する振動抑制方法において、運動体に付加された第1の入力によって運動体に生じる第1の応答に、第1の応答の変位が0となる時刻において、第2の入力を運動体に付加して生じる第2の応答を重ね合わせて、互いに打ち消し合うインプットシェイピング法を適用する場合に、n自由度かつ時変の振動におけるモード間の相互作用によるエネルギの移り変わりをゼロと見なして、固有振動数および減衰比の時間変化によるモード間のエネルギの移り変わりを考慮して、第1の入力に対する第2の入力の大きさを決定する、振動抑制方法を提供する。本発明の第4態様によれば、各振動モードにおいて、第1の入力に対応する入力を打ち消す第2の入力に対応する入力の大きさを個別に算出して、振動モードごとに算出された入力を合成することで第2の入力の大きさを決定する、第1態様から第3態様のいずれか1つに記載の振動抑制方法を提供する。本発明の第5態様によれば、運動体が多関節型ロボットアームであって、ロボットアームに生じるn自由度かつ時変の振動を抑制する、第1態様から第4態様のいずれか1つに記載の振動抑制方法を提供する。本発明によれば、加速または減速運動に伴って運動体に生じる、n自由度かつ時変の振動を、インプットシェイピング法を用いて抑制する場合に、n自由度かつ時変の振動における複数のモード間のエネルギの移り変わりのうち、モード間の相互作用によるものを無視し、固有振動数および減衰比の変化によるものを考慮して、第2の入力を決定するという手法を採用している。そのため、第1の入力によって生じたn自由度かつ時変の振動第1の応答を抑制できる第2の入力を算出することが可能となる。よって、加速または減速運動に伴って運動体に生じるn自由度かつ時変の振動を、インプットシェイピング法適用して効果的に抑制する振動抑制方法を提供することができる。以下に、本発明にかかる実施の形態を図面に基づいて詳細に説明する。実施の形態本発明の実施の形態にかかる運動体の振動抑制方法について、運動体の一例であるロボットアーム10を用いて説明する。多関節型ロボットアームの運動方程式図1に示すように、ロボットアーム10は、多関節型ロボットアームであり、n個のリンクL1,,Lnと、それぞれのリンクLを動作可能に接続するn個の関節J1,,Jnとを備えているnは自然数。また、関節Jの部分を示す図2および図1に示すように、このロボットアーム10のそれぞれのリンクLおよび関節Jは、1つの平面内で動作可能となっている。図2に示すように、関節J1の角度をθ1、その指令角をθ01、ばね定数をk1とし、関節J1以外の任意の関節jの角度をθj、その指令角をθ0j、ばね定数をkjとする。上付き添え字Tは転置を表すとし、数1~数3のように定義する。数1~数3の定義と、ラグランジュの方法により、多関節型ロボットアーム10の運動方程式は、数4のように表される。ここで、Jθは慣性行列、h左辺第2項はコリオリの力と遠心力を表す項、gθは重力付加を表す項、CVは粘性減衰行列である。テイラー展開次に、非線形な式である数4をθ=θ0まわりの1次のテイラー展開により線形化を行うと、数5が得られる。なお、数5におけるφ、Ch、Km、Kh、Kgは、数6~数10の通りである。多自由度線形時変系のモード分解モード解析は、本来は線形時不変系を対象とする解析手法の一つである。ここでは、このモード分解を連続時間で記述される線形時変系に適用し、入力値整形に利用する。一般にn自由度線形時変系の運動方程式は、n次正方実行列の質量行列M、減衰行列C、剛性行列K、n行縦ベクトルの入力f、偏差角ベクトルφを用いて、数11のように表せる。多関節型ロボットアーム10の場合には、指令値θ0に対して参照指令値θrefを与えるものとして、数12~数15とすれば、数5は数11の形の多自由度線形時変系とみなすことができる。ただし、href、grefは数16、数17のように表される。ここで釣り合いの位置からの偏差角ψは数18の通りである。この時間に関する1階微分と2階微分はそれぞれ数19、数20のようであるので、釣り合いの位置からの偏差角ψに関する運動方程式は数21のようになり、数22のように表される。Enをn次単位行列とすると、数21の運動方程式に対応する状態方程式は、数23~数26にように表される。ここで、Aの時刻tにおける瞬間的な固有ベクトルと固有値について、常に共役なペアがn個存在する場合、つまり、全モードの応答が振動的であり、過減衰または臨界減衰でない場合についてのみ考える。Aの時刻tにおけるi番目の瞬間的な固有ベクトルをpi、対応する瞬間的な固有値をλiとし、固有ベクトルからなる行列Pと固有値を対角要素にもつ行列Λを、数27、数28のように定義するiは自然数。ただし、数29の関係式が成立するようにiを選ぶものとする。ここで、モード座標qを数30のように定義する。すると、状態方程式はモード座標について、数31のように表すことができる。ここで、数31に示される状態方程式における右辺第一項である数32にて表される行列のうち対角要素のみを取り出した行列Dを考え、この行列Dを数33のように表す。数32の行列において、非対角要素の絶対値が対角要素の絶対値に比べて十分に小さいことから、数33の行列Dを用いて、数31の状態方程式を数34のように近似できる。なお、この状態方程式の近似の考え方については後述する。数31の状態方程式では、数32の行列において、非対角要素に有意な値があるため、このままでは方程式を解くことができない。しかしながら、数31の状態方程式を数34の状態方程式にように近似することで、各モードの応答を独立な方程式を解くことで求めることが可能となる。これにより、i次モードのインパルス応答関数を数35のように導くことができる。なお、τは単位インパルスを与える時刻である。モードごとの入力値整形次に、ロボットアーム10に付加された第1の入力によってロボットアーム10に生じる第1の応答に、ロボットアーム10に第2の入力を付加して生じる第2の応答を重ね合わせて、互いに打ち消し合うような第2の入力の値を、振動モードごとに算出方法について説明する。ここで、インプットシェイピング法におけるインパルス応答の重ね合わせの考え方について、図3に示す。図3ではインパルス応答の変位振動量:縦軸と時間横軸との関係を示している。図3に示すように、ある時刻τに与えられる第1のインパルスP1第1の入力に対する第1の応答Hi第1の応答を打ち消すため、第2のインパルスP2第2の入力を与える時刻をτ数式上では、τの上にと表す。以降同じ。>τとし、第1のインパルスP1に対する第2のインパルスP2の大きさの比を実数α数式上では、αの上にと表す。以降同じ。とする。ここで、αとτは未知数であり、これらは以下のように求めることができる。まず、第1のインパルスP1により生じる第1の応答とインパルスP2により生じる第2の応答とが、t>τにおいて互いに打ち消しあうための条件は、数36にて表される。数35より、t>τでは、数37の関係式が成り立つため、数36は数38のようになる。一般にHi≠0であるので、数36を満たすための条件は、数39のようになる。この数39は、第2のインパルスP2の大きさを求める式である。τは未知であり、αは実数であるので、数40に示す関係を満たさなければならない。数40を用いて、時刻τで与えられた第1のインパルスP1に対する第1の応答の虚部が0になる時刻τを求めることができる。しかし、数40を満たすτは無数に存在する。そこで、数40を満たす最小のτをτとする。なお、τは力学系においては、第1のインパルスP1によって生じる第1の応答に対して、変位が最初に0になる時刻である。また、は、時不変系で言うところの半周期に相当する。時変系であれば、τは第1のインパルスP1を与える時刻τにより異なるので、τをτの実関数として、数41にように表す。さらに、αiを数42のように定義すれば、数43の関係式が成り立つ。次に、このインパルスの組み合わせを用いて、i次モードに対する任意の入力uiによる応答を打ち消すことができる追加の入力wiを求める。任意の入力uiはある時刻τSからτEの間にのみ任意の複素数をとるものとし、それ以外の時刻を乗じ、τについてτSからτEまで積分すると、τ>ηiτEのとき、数44が成り立つ。数44の左辺第一項は畳み込み積分の形になっており、入力uiに対する応答を表している。左辺第二項についても畳み込み積分の形にするため、積分変数をτからτに変換する。なお、ηiの逆関数をηi1、ηi1の1次導関数を数45のように表すと、数46、数47のようになる。これらより、数44は数48のように表される。ただし、数49~数51のような関係になる。数48の左辺第一項は、uiによる応答、左辺第二項はwiによる応答を表すことになる。さらに、これらの足し合わせがt>τEにおいて0になることを、数48は示している。つまり、uiwiをi次モードに対する第2の入力とすれば、このi次モードに生じた残留振動は抑えられることになる。なお、第1のインパルスP1が与えられた時刻τから第1の応答の虚部が最初に0になるまでの時間τ-τを時変系における半周期と呼ぶことにすれば、数51においてηi1は、時刻τから見て半周期前の時刻τを表すことになる。また、αiは時刻τに与えられた第1のインパルスP1に対する第1の応答が、半周期の間にどれだけ減衰するのかを表しており、uiは半周期前の入力の大きさを表すことになる。また、数52は、時変系における固有周期の変化に対する補正量である。ここで、時変系において連続的に第1のインパルスが入力される場合を考える。図4に示すように、第1の入力として、連続的にインパルスが入力された第1のインパルス列P3により生じる第1の応答を打ち消すように連続的にインパルスが入力された第2のインパルス列P4を考える。時変系の場合は、第1のインパルス列P3の時間間隔がdτであるのに対して、第2のインパルス列P4の時間間隔は数48における積分変数の変換により数53のように表される。そのため、この時間間隔を補正するために追加の入力に対しては、数52に示される時変系における固有周期の変化に対する補正量を乗じる必要がある。ここで、元々の入力uiに対してwiを足し合わせると、例えばこの入力を力とみれば、系に加わる力積が単純に増えるということになる。そこで、入力の正規化係数を数54で定義する。数48に任意の複素定数を乗じても右辺が0となることから、uiwiに任意の複素定数を乗じても2つの入力は打ち消し合うことになる。そこで、数54で定義される正規化係数を用いて、各モードにおける整形された入力値は、数55のように表すことができる。複数モードを対象とした入力値整形次に、複数モードを対象とした入力値整形の方法を示す。まず、入力fとして任意の参照入力frefを与えると、対応する各モードへの入力uiは、数56のように表される。この数56をもとに数55により整形された入力uiをモードごとに求め、任意の実定数κ1κnを用いて、数57のように表す。BをBの疑似逆行列とすると、数58で表されるfshapedが、数56で表される多自由度線形時変系において、全てのモードの残留振動を抑制する入力となる。κiuiがi次モードの残留振動を抑制する入力であり、さらにその他のモードを励起しないため、そのような入力の足し合わせであるfshapedは、結局、全てのモードについての残留振動を抑制することになる。加減速パターンの生成次に、ロボットアーム10の残留振動を抑制する加減速パターン数59を生成する方法について説明する。本手法により加減速パターンを整形した結果、対応する速度パターン数60、変位指令値数61は、数62、数63から大きく変わらないものとする。これにより、加減速パターンを整形した後もモード形状の時間推移は変わらないものとみなせる。すると、数15、数22より、数64で表される加減速パターン数59が残留振動を抑制する加減速パターンとして生成される。ただし、κ=κ1κnは動作終了時に、対応する速度パターン数60が0となるように、信頼領域Dogleg法によりκを変化させながら数56~数58、数64を繰り返し評価することにより決めることが望ましい。また、Γは対角実行列であり、整形後の指令値数61が参照指令値数63と同じ値で終わるようにその要素を決めることが望ましい。状態方程式の近似についてここで、上述した数31の状態方程式において、数32にて表される右辺第一項の行列のうちの非対角要素を無視してすなわち0として、対角要素のみを取り出した行列Dを用いて、数34のように近似した本発明の考え方について、さらに説明する。2関節ロボットアーム20を例として説明する。図5に示すように、2関節ロボット20では2種類の振動の形態が考えられる。図5aに示す振動の形態が1次モードであり、図5bに示す形態が2次モードである。このような振動の形態は、1次モードについては固有ベクトルp1によって決まり、2次モードについては固有ベクトルp2によって決まる。実際には、この2つの振動モードが合成されて混在してアームが振動し、2つの振動の形態が合成された動きとなる。時変系の場合は、数31において、行列数32の非対角要素が0でないため、q1の値がq2に影響を及ぼし、q2の値がq1に影響を及ぼすことになる。すなわち、時変系の場合には振動モード間での相互作用が発生することになる。これに対して、時不変系の場合には振動モード間での相互作用は発生しない。時変系におけるこの相互作用の大きさが行列数32の非対角要素の値の大きさと対応する。しかしながら、運動体がロボットアームの場合には、振動モード間での相互作用はあまり大きくない。具体的には、行列数32の対角要素モード毎の自らの時変性が示されている要素に比して、非対角要素モード間の相互作用の大きさが示されている要素が十分に小さい。そこで、この非対角要素を無視して0とみなし、数32の対角要素のみを取り出したDという行列数33を用いることで、数31の状態方程式を数34のように近似することで、モードごとの応答を独立な方程式を解くことで算出することを可能とするというのが、本発明の考え方である。なお、このような観点からは、本発明において振動モード間の相互作用とは、数31において、q1の値がq2に影響を及ぼし、q2の値がq1に影響を及ぼす作用を意味すると言うことができる。また、モード間のエネルギのやり取りの観点からも、この相互作用について説明できる。系に減衰が無く、外力も与えられていないとすれば、エネルギ保存則が成り立つ。1次モードのエネルギをE1、2次モードのエネルギをE2とすれば、E1E2は常に一定の値となる。時不変系であれば、さらにモードごとにエネルギ保存則が成り立つ。つまり、E1とE2がそれぞれ一定の値となる。これに対して、時変系では、モードごとにエネルギ保存則が成り立たず、モード間でエネルギのやり取りが発生する。例えば、E1が大きくなるときには、E2がその分だけ小さくなる。逆に、E1が小さくなるときには、E2はその分だけ大きくなる。ここで、このような時変系におけるモード間でのエネルギのやり取りには2つの要因がある。1つ目は、行列数32の対角要素、つまり固有振動数と減縮比が時間により変化することである。2つ目は、行列数32の非対角要素により異なるモード間にて相互に作用し合うことによるものである。本発明では、時変系において、1つ目の要因を考慮しながら、2つ目の要因を無視することにより、数34のような近似式を得るという考え方に基づいている。すなわち、本発明の振動抑制方法では、n自由度かつ時変の振動における複数のモード間のエネルギの移り変わりのうち、モード間の相互作用によるものを無視し、固有振動数および減衰比の変化によるものを考慮して、第1の入力により生じる第1の応答を打ち消すような第2の入力の大きさを決定する、という考え方に基づいている。次に、このような考え方に基づく、n自由度かつ時変の振動を抑制する振動抑制方法について、具体的な実施例を用いて説明する。実施例1:1自由度線形時変系まず、実施例1として、1自由度線形時変系について説明する。具体的には、上記実施の形態にて説明した手法を、図6に示す可変ばねを有する1自由度振動系に適用し、本手法の時変系に対する有効性を確かめる。特に固有周期の変化に対する補正が、高速にパラメータが変化する際に重要となることを本実施例1にて示す。運動方程式図6に示す可変ばねを有する1自由度振動系を考える。ばねの自然長は一定である。質量は1に正規化し、図示右向きを正として質量の位置をξとする。固有振動数をωnとすると、ばね定数はωn2となる。また、右向きを正として指令値をξ0とする。すると、この系の運動方程式は数65で表される。ここで、質量の変位と指令値の偏差z=ξ-ξ0について運動方程式を書きなおすと、数66で表される。なお、初期条件は数67にて表される条件とする。この偏差に関する運動方程式に対応する状態方程式は、数68のように表される。数68において、Aの時刻tにおける瞬間的な固有ベクトpと固有値λはそれぞれ数69、数70にて表される。モード座標への変換次に、モード座標への変換を行う。数71および数72とすれば、数73にて表されるモード座標を用いて、状態方程式は数74と書き換えられる。ここで、数75に表されるようにし、数74の非対角要素を無視すれば、数74は数76のように近似できる。各モードへの入力を数77に示すようにすると、モードごとの運動方程式は、数78にて表すことができる。ここで、数79とすれば、dの時間積分は数80と書けるので、インパルス応答関数は、数81と求まる。入力値整形次に、数40を満たす数82の最小値をτとすると、数83となる。これをτとτについて解き、それぞれτ=η,τ=η-1とする。また、α=-Hτとする。数51より、数84により励起される振動を打ち消す加速度指令値は、数85のように求めることができる。シミュレーションシミュレーションの例として、固有振動数が数86のようになる場合について考える。ここで、ω0はt=0における初期固有振動数、aは固有振動数の時間変化率を表す実定数である。すると、各種関数が、数87~数91のように、初等関数のみを使って求めることができる。参考のため、a=π,ω0=πの場合の数90を図7に示し、数88、数89および数91を図8に示す。ここで、初速度0の状態から数92の状態まで加速するための参照加速度指令値として、数93を与えるものとする。非対角要素を無視した数76の近似の妥当性を調べるため、数93を入力として与えたとき、数74による応答qと近似式である数76による応答qを、それぞれルンゲクッタ法の一種であるDormandPrince法により図9のように得た。図9に示されるように、2つの応答はよく似ており、数76の近似は妥当であると言える。参照加速度指令値数93が与えられれば、数85を用いると、初等関数のみを用いて明示的に振動を打ち消す加速度指令値数85を求めることができる。このようにして求めた振動を打ち消す加速度指令値を図10に示す。ここで、残留振動を抑えつつ、加速期間終了時t=ηに速度数94が1m/sとなるように、加速度指令値を数95とする。参照加速度指令値数96、整形後の加速度指令値数97、Choの方法により整形した加速度指令値を図11に示し、それぞれに対する応答を図12に示す。図12より明らかなように、本発明の手法で残留振動が大幅に抑制されていることが分かる。なお、Choの方法については参考文献ChoJKandParkYSVibrationreductioninflexiblesystemsusingatimevaryingimpulsesequenceRoboticaVol13Issue3pp303313に開示されている。特に着目すべき点は、本発明の手法では追加の入力に相当する部分がChoの方法に比べて大きい点である。これは、時間が経つにつれて固有周期が短くなっているためである。このような場合には参照入力と追加の入力が与える力積を等しくするために、追加の入力の方を大きくなるよう補正する必要がある。本発明の手法では、この固有周期の変化に対する補正を数52で示す項で行っている。一方で、図13に示すように、Choの方法により得られる加速度指令値は数52で示す項を1として整形した加速度指令値と一致することから、Choの方法では固有周期の変化に対する補正がされてないと言える。このことから、固有周期の変化に対する補正が時変系の残留振動抑制に必要であることが分かる。また、数75の積分が解析的に求まり、さらに、数83が解析的に解けるのであれば、本発明の手法では単純な演算のみで入力値整形を行うことができる。実際、固有振動数の変化が数86のように表せるため、初等関数のみを用いて入力値整形を行うことができた。このように、1自由度時変振動系を対象とする場合には、本発明の手法による入力値整形の演算量は少なくて済む。このため、多くの行列演算を必要とするChoの方法と比べて、本発明の手法は実機に実装するにあたって有用である。高速時変系に対する効果本発明の手法の特徴は、パラメータが高速に変化する場合にでも適用可能な点にある。先に説明した例を用いて、固有振動数の変化率aが大きい場合でも残留振動が低減可能であることを数値計算によるシミュレーションにより示す。まず、数98に対する応答の絶対値zのt>1における最大値をzrefmaxとする。また、数99に対する応答の絶対値zのt>ηにおける最大値をzmaxとする。すると、zmax/zrefmaxが入力値整形によりどれだけ残留振動が低減されたかを表す指標になる。そこで、aをω0/10から10ω0の間で変えながら運動方程式をDormandPrince法により解いた。どのようにzmax/zrefmaxが変化するのかを、Choの方法による結果とともに図14に示す。図14に示すように、本発明の手法では、固有振動数が初期固有振動数ω0に対して単位時間あたりで10倍も変化するような場合でも、残留振動が3以下にまで低減されていることが分かる。また、Choの方法と比べると、固有振動数が高速に変化する場合でも効果的に残留振動が抑制されていることが分かる。効果本発明の手法では、多自由度の複雑な系であっても、入力値整形に必要な関数を補間などの方法により求めることで計算時間の短縮を狙うことができる。αは減衰などによるエネルギ消散量に関連しており、インパルス応答が分かれば容易にαの関数形状を想像することができる。例えば、図7より、振幅が時間とともに増大していく場合には、対応するαは常に1より大きくなり、振幅の変化が小さくなるとαは1に近づいていく。また、η-1についても瞬間的な固有周期の時間変化率に関連している。図7に示したように固有周期が短くなる場合には、2つ目のインパルスの値は大きくなるようにせねばならず、η-1は常に1以上となる。また、固有周期をT数100とすると、その変化率は数101のように表され、時間の経過とともに固有周期の変化が小さくなることが分かる。このことから、時間の経過ともにη-1が1に収束していくことが予想され、実際、図8ではそのようになっている。このように、固有周期の変化をあらかじめ知ることができればη-1の形状を予測できる。すると、これらの関数は補間などの手法により得ることができ、設計者が補間に必要なだけの点数を選ぶことが可能になる。関数の形状さえ再現できるのであれば、不必要に時間間隔を細かく区切る必要が無くなる。このことから、本発明の手法では入力値を整形するのに必要な演算量作業量を低減することが可能であると言える。実施例2:水平2関節型ロボットアーム次に、実施例2として、水平2関節型ロボットアームへの適用例について説明する。図15に示す水平面内で動作する2関節ロボットアーム20について、残留振動を抑制するための加減速パターンを生成し、数値計算と実験により上述にて説明した入力値整形の有効性を検証した。この検証にて用いた実験装置実験モデルを図16に示す。この実験装置は、第1リンク11を水平面内にて回動させる第1関節21にACサーボモータRSF14A100を用い、第2リンク12を水平面内にて回動させる第2関節22にACサーボモータRSF14A50を用いた水平面内で動作する2関節ロボットアーム20である。これらのACサーボアクチュエータには波動歯車装置が用いられている。各関節のモータM1、M2は、それぞれACサーボドライバHA655により図17の表に示す設定で位置制御モードにより制御されており、外部指令装置からパルスを送ることで駆動する。第1関節21は1パルスあたり45×10-3度、第2関節は90×10-3度回転する。外部指令装置にはDS1104RDコントローラボードのDAC出力を用い、1秒あたり最大で25×104パルスを同時に2つのACサーボドライバに送ることができる。また、動作終了後のアーム先端の残留振動の大きさの測定には、レーザー変位計HLC105BHKを用いた。なお、この測定値はフィードバックすることなく、測定のためだけに用いた。運動方程式関節jの回転角度をθj,その指令値をθ0j,リンクlの質量をml,リンクlの関節j周りの慣性モーメントをJj,リンクlの長さをbl,関節jからリンクlの質量中心までの距離をlj,関節jのねじり剛性をkjとする。運動方程式は、数103~数111の条件に基づいて数102のように表される。数5~数10によるテイラー展開により、参照指令値をθref=θref1θref2Tとすると、水平2関節型ロボットアームの運動方程式は数21で表される。ただし、各パラメータは数12~数15より、数112~数115のようになる。ただし、数116~数122に示す関係にある。なお、水平多関節型ロボットアームの場合には重力項が無いため、数18よりψ=φである。パラメータ同定本発明の手法により加減速パターンを生成するために必要なパラメータの同定方法について述べる。まず、電子はかりによる測定と幾何形状から求まるパラメータを図18の表に示す。ここで、未知のパラメータはJ1、k1、k2、Cvである。これらは以下に説明するようにして同定した。モーダルパラメータの同定まず、全てのモータをサーボオンにした状態で、第1関節21に対する指令角をθ01=0,θ02=θr2として固定したときの自由応答を測定した。ここで、θr2は、数123のように表される。自由応答は、図19のように、加速度ピックアップをアーム先端から9mmのところに貼り付け、その裏側をインパクトハンマにより加振し、アーム先端の加速度を測定した。θ02=θr2のときにFFTアナライザにより得られたアクセレランスをGrとする。なお、サンプリング周波数を256Hz,計測時間を16秒,ωk=2πfk,fkを0,00625,0125,,128Hzとした。ここで、一般粘性減衰をもつ2自由度系のアクセレランスLr数124に対し、Grを1~100Hzの領域で曲線適合させるため、数125にて表される値が最小となるようパラメータUri,Vri,ωrni,ζriを信頼領域Reflective法により求めた。一例として、θ02=θ42=π4radのときに計測した伝達関数G4と、求まったパラメータをもとに求めた伝達関数L4を図20に示す。なお、θ02=π2,7π12radのときは、鉛直方向の1次固有振動数と水平方向の2次固有振動数が近く、インパクトハンマによる打撃によって鉛直方向の振動が少なからず励起されるとともに、加速度ピックアップが鉛直方向の加速度をも拾ってしまうため、曲線適合は一般粘性減衰をもつ3自由度系のアクセレランスLr数126を用いて行った。慣性モーメントとねじり剛性の同定次に、ωrniをi次モードの固有振動数とみなし、このパラメータをもとにJ1,k1,k2を以下のように同定した。数5において数127に示すようにしたときの偏差角φに関する運動方程式数128より実固有値解析によって求まる固有振動数をωrniとする。そして、数129に示す値が最小となるよう信頼領域Reflective法によりJ1,k1,k2を図21に示す表のように同定した。アクセレランスの曲線適合により求めた固有振動数ωrniと、同定したJ1,k1,k2を運動方程式数130に代入し、実固有値解析により求まる固有振動数ωniを図22に示す。なお、第1関節21に用いられているモータM1のねじり剛性は、ねじり角が小さい場合には47×103Nmradである。また、第2関節22では34×103Nmradである。同定したねじり剛性k1,k2がこれらの値よりも小さいのは、関節に用いられている軸の軸径が小さいことにより機械的剛性が低くなっている部分があることと、サーボ機構により系全体の剛性が下がっていることによる。粘性減衰行列の同定次に、粘性減衰行列Cvを求める。ζriをi次モードのモード減衰比とみなし、θ02=θr2のときの粘性減衰行列が、定数αr,βrを用いて数131というRayleigh減衰の形で与えられるとすると、定数αr,βrは、i次モードに関して数132という関係を満たす。これより、1次のモーダルパラメータωrn1,ζr1と、2次のモーダルパラメータωrn2,ζr2を用いて、定数αr,βrがそれぞれ数133、数134と求まる。数131、数133、数134よりCvrが求まるので、その平均をとることで数135のようにCvを決めた。座標系の定義アーム先端の面内振動を評価するため、図23のようにベースから見た座標系ΣB,アーム先端の指令位置から見た座標系ΣH,動作終了後のアーム先端の位置から見た座標系ΣEを定める。なお、θE=θ1Eθ2ETは動作終了時における関節の角度である。図23より、座標系ΣBにおけるアーム先端の位置は、数136、数137のようになる。座標系ΣHにおけるアーム先端の位置は、同次変換数138により、数139、数140となる。この座標系ΣHは動作中のアーム先端の振動を評価するために用いる。同様にして,座標系ΣEにおけるアーム先端の位置は、数141、数142となる。座標系ΣEにおけるアーム先端の位置は、レーザー変位計により直接測定できるもので、動作終了後のアーム先端の残留振動を評価するために用いる。突き出し動作次に、図24に示す突き出し動作について考える。動作開始時の各関節の角度θSは数143にて表され、動作終了時の関節の角度θEは数144にて表される。最大角加速度を4πrads2,最大角速度をπradsとし、2つの関節が同時に動作終了するようにすると、最短時間制御は図25に示すような指令値となる。これを参照指令値θrefとし、対応する加減速パターン数145,速度パターン数146についても図25に示す。ここで、i次モードの疑似固有振動数と疑似減衰比をそれぞれωi=di,ζi=-realdiωiとすると、動作中の疑似固有振動数と疑似減衰比は図26のように変化する。1次固有振動数と2次固有振動数は、動作開始時にはそれぞれ79Hz,183Hzであるが、動作終了時にはそれぞれ62Hz,280Hzとなる。近似の妥当性この突き出し動作において、上述にて説明した加減速パターンを導出する際に行った近似の妥当性を検証する。元々の非線形な運動方程式数4、剰余項を無視して線形化した運動方程式数21、非対角項を無視することによって導いたモードごとの非連成な運動方程式数34をそれぞれ数値計算により解いた。その結果得られる偏差角をそれぞれψ=ψ1ψ2T,ψL=ψL1ψL2T,ψM=ψM1ψM2Tとし、図27に示す。また、非線形な運動方程式数34の解ψからの誤差を数147、数148のように定義し、図28に示す。なお、数値計算にはDormandPrince法を用いた。図28より、線形化したことによる誤差が5以下であり、モード分解をしても誤差が5.5以下であることから、この突き出し動作における数34の近似は妥当であると言える。加減速パターンの生成加減速パターンを生成するにあたって、数149を解析的に求めることは困難であるため、差分を用いて数149を計算する。しかしながら、図25に示す加減速パターンは時刻tに関して不連続であるため、微小時間ΔtPを用いて数150とすると、Pが不連続になる時刻tにおいて、数151の各要素が他の時刻のときと比べて非常に大きな値をとり、生成される加減速パターンは残留振動を抑制しなくなる。そこで、数151を数152とする。ただし、med関数は引数として与えられた行列のうち、フロベニウスノルムが中央値をとる行列を返す関数である。このように数149を定義しなおすと、ΔtPを小さくしても、Pが不連続な点において数149の各成分が異常に大きな値をとることがなくなる。この数149を用いて、数64により得られる加速度パターン数59と対応する数60数61を図29に示す。ただし,ΔtP=0001sとした。整形前の図25と整形後の図29とを比べることにより、両者の速度パターンと指令値に大きな差異がないことが分かる。よって、この場合にはモード形状の時間推移は整形後も変わらないものとみなすことが妥当であると言える。比較のため、InputShaperの中で最も基本的なZVShaperにより数153を整形した加減速パターン数154を図30に示す。なお,ZVShaperのパラメータとして、同定した1次モードのモーダルパラメータの平均値を数155、数156のように用いた。つまり、この加減速パターンは1次モードの固有振動数が69Hz,減衰比が0037である線形時不変系の残留振動を抑制する入力に相当する。時変系である2関節型ロボットアームに適用しても、ある程度残留振動が抑制されることが期待されることから、本発明の手法の比較対象として以下のように取り上げる。ここで、動作を開始してから最後に指令値を送るまでの時間を動作時間とし、各加減速パターンにより生成される指令値の動作時間を比較する。参照指令値ではt=0834sで動作終了するのに対し、整形後の指令値ではt=0915sで動作終了する。この参照指令値に対する動作終了時刻の遅れは、動作終了時の1次モードの半周期πω1E=0081sに相当する。ここで、ωiEは動作終了時のi次モードの固有振動数である。本発明の手法では、一般に参照指令値の動作時間に対して1次固有周期の半分πω1Eだけ動作が遅れる。これは、モードごとに2パルスで入力値整形を行うため、最も固有周期の長い1次モードにおいて入力が終わるまでにかかる時間が動作時間となるからである。そのため、高次のモードを扱う場合でも1次固有周期の半分だけ動作が遅れることになる。これに対し、Choの方法などの2パルスのConvolved型のInputShaperでは、一般に、残留振動を抑制したいモードの固有周期の合計の半分数157だけ動作終了時刻が遅れる。そのため、本発明の手法の方が動作終了時刻を短縮できることが分かる。なお、ZVShaperでは平均的な固有振動数を用いて加減速パターンを整形するため、動作終了時刻はπωZVだけ遅れる。すなわち、この突き出し動作の場合はt=0907sで動作が終了する。結果と考察生成した指令値θshapedを与えた場合に、参照指令値θrefを与えた場合やZVShaperにより生成した指令値θZVを与えた場合と比べて残留振動がより抑制されることを数値計算と実験により検証した。なお、以下に示す数値計算による結果は、数4をDormandPrince法により解いて得られるθから計算したものである。まず、動作終了後の残留振動を評価するため、数値計算により求めた動作終了後のyEを図31に示し、実験においてレーザー変位計により測定したyEを図32に示す。図31、図32より、実験結果が数値計算の結果をよく再現していることが分かるとともに、本発明の手法がZVShaperと比べてより効果的に残留振動を抑制していることが分かる。ここで、動作開始から振動振幅が025mm以下に収まるまでの時間を位置決め時間とすると、実験結果より、θrefを与えた場合には位置決め時間が25秒以上、θZVを与えた場合には161秒となった。それらに対し、θshapedを与えた場合には、動作終了時には振動振幅が025mm以下となっており、位置決め時間は092秒となった。これより、本発明の手法を用いると位置決め時間が短縮され、作業効率を向上させられることが分かる。次に、動作中の振動を評価するため、数値計算により求めた動作中と動作終了後のxHとyHを図33に示す。図33より、θshapedを与えると動作中の振動も抑制されることが分かる。これは、時変性を考慮して入力値整形をしているためである。それに対し、θZVを与えた場合には、動作中である03~058秒付近において偶然にも振動が収まっているが、時変性を考慮せずに入力値整形をしているため、動作中のその他の期間と動作終了後において振動が生じる結果となっている。最後に、数値計算の結果より、座標系ΣBにおいてアーム先端が描く軌跡を図34に示す。θshapedを与えたときのアーム先端の軌道は、必ずしも参照指令値θrefを与えたときと同じ軌道を描くとは限らないことが図34より分かる。そのため、アーム先端の軌道に制約がある場合などには、参照指令値θrefを変えて加減速パターンを生成しなおす必要がある。この突き出し動作の例の場合には、アーム先端がyB>0の領域に大きく侵入していることが分かる。これは、図25、図29に示すように、第2関節の参照加速度数158に比べて、生成した第2関節の加速度パターン数159では動作初期において大きく加速し、第1関節よりも速く動作するためである。このため、例えばyH=001mの位置に障害物があるような場合には、参照指令値の第2関節の角速度を落とすなどしてから加減速パターンを生成しなおす必要がある。動作時間と残留振動非減衰1自由度線形時不変系において、固有周期の整数倍の時間だけ一定外力を与えると残留振動が生じないということから、図25のような一定加速度により加減速する加減速パターンを用いても、加速時間と減速時間を調整することにより残留振動をある程度まで抑えられることが予想される。同様に、時変性を考慮しないZVShaperを用いても加減速時間と系の固有周期の関係により残留振動が効果的に抑えられる場合があると考えられる。しかし、これらの方法では系の固有周期により動作時間を決めなければならないため、必ずしも位置決め時間を短縮することができるとは限らない。それに対し、本発明の手法により生成する加減速パターンは、時変性を考慮しているため、系の固有周期と加減速時間の関係によらず残留振動を抑制でき、位置決め時間を短縮できると考えられる。このことを示すため、動作時間と発生する残留振動の大きさの関係を数値計算と実験により調べた。まず、一定加速度により加減速する動作について考える。θdを数160とする。参照指令値として最大角速度をvmax,最大角加速度をAmaxとし、2つの関節の動作が同時に終了する最短時間制御を行う。θd<v2maxAmaxの場合には図35のような三角形速度パターンとなり、θd≧v2maxAmaxの場合には図36のような台形速度パターンとなる。このとき加速時間ta、動作時間tmは数161、数162のようになる。ここで、vmax=πradsとし、最大角加速度Amaxを与えると、ta,tm,および数163が決まる。この加減速パターンに対応する参照指令値θrefを与えたとき、数値計算により得られる動作終了後のアーム先端の最大変位maxyEを図37に示す。ここで、横軸は動作終了時の1次固有周期に対する加速時間の比数164である。図37より、加速時間が1次固有周期の整数倍に近いと発生する残留振動が小さくなることが分かる。しかし、2関節ロボットアームでは非線形性や時変性により、加速時間が1次固有周期の整数倍であっても少なからず残留振動が発生することが分かる。次に、数値計算により得られる動作終了後のアーム先端の最大変位maxyEを、横軸を動作時間tmとして図38に示す。また、ZVShaperを適用した場合の最大変位と、本発明の手法を適用した場合の最大変位も示す。なお、最高角速度vmaxをπradsとしているので、参照指令値の最高角加速度を大きくして最も制御時間が短くなるようにしても、一定加速度の加減速パターンでは、数165だけ動作に時間がかかる。同様にして、ZVShaperにより整形した加減速パターンでは、数166だけ動作に時間がかかり、本発明の手法により生成された加減速パターンでは数167だけ動作に時間がかかる。ここで、実験においてレーザー変位計により測定した最大変位を図39に示す。図38、図39より、動作時間が短い場合には本発明の手法がZVShaperと比較して動作時間によらず効果的に残留振動を抑制できることが分かる。図39より、動作時間によっては本発明の手法よりも発生する残留振動の大きさがZVShaperを用いた方が小さくなる場合があることが分かる。しかし、アーム先端の負荷が変わるなどして固有周期が変化すれば、残留振動が小さくなる動作時間は変化する。そのため、残留振動が生じない動作時間を見つけるのは困難である上、動作時間が限られてしまうため、実際にZVShaperを多関節型ロボットアームに適用して任意の動作の残留振動を抑制するのは難しい。これに対し、本発明の手法で加減速パターンを生成すると、残留振動を抑制しながら動作時間を短縮できることが分かる。なお、本発明の手法において、動作時間が短い場合、すなわち、角加速度が大きい場合には、角加速度が小さい場合と比べて残留振動の抑制効果が小さいことが図38より分かる。本発明の手法も含む各種InputShaping法では、完全に打ち消し合わせられない振動的な応答が残る場合に、加振力に相当する角加速度が大きいことにより残留振動が大きくなる。時変性を考慮したにも関わらず本発明の手法により振動的な応答が残るのは、加減速パターンの整形前後においてモード形状の時間推移が比較的大きく変わったためだと考えられる。Amax=618rads2のときの参照指令値の加減速パターンとそれに対応する速度パターン、指令値を図40に示す。また、この参照指令値をもとに本発明の手法で生成した加減速パターンとそれに対応する速度パターン、指令値を図41に示す。両者の速度パターンを比較すると、特に減速期間において速度パターンが異なることが分かる。しかし、加減速パターンを生成する際には加減速パターンの整形前後においてモード形状の時間推移が変わらないことを前提としているため、この前提をもとに生成した加減速パターンでは残留振動が完全に抑制できないということになる。ここで、図38と図39を比較することにより、動作時間が長い場合、すなわち、角加速度が小さい場合には、数値計算による結果と比べて実験ではあまり残留振動が抑制されていないことが分かる。これは、指令角の量子化の精度が粗いためである。動作終了前後における関節1に対する理想的な指令角θshaped1と、実際にサーボドライバに送信されるデジタル化された指令角を図42に示す。図42aには角速度が大きい場合の指令角を示し、図42bには角速度が小さい場合の指令角を示している。図42より、角加速度が小さい場合には、動作終了直前に非常にゆっくりと回転するように指令をするが、量子化の精度が粗いために、断続的に角度を指令することになる。このため、動作時間が長い場合に残留振動が発生してしまっていると考えられる。このことから、本発明の手法を実際にロボットアームに適用する場合には量子化の精度に注意することが望ましい。振り回し動作図43に示すアームを振り回す動作により、ZVShaperとの違いを明確にする。この動作は、図24に示す突き出し動作において第2関節22の回転方向を逆向きにしたものである。つまり、第2関節22への指令角θ02の符号を反転させただけであり、固有振動数の時間推移は等しくなる。しかし、1次モードに寄与するトルクが大きくなるため、その分だけ残留振動が大きくなり、時変性を考慮しない従来法との違いが顕著になる。ここで、最大角加速度を4πrads2,最大角速度をπradsとし,2つの関節が同時に動作終了するような最短時間制御を参照指令値θrefとする。これをもとに、ZVShaperでθZVを、本発明の手法でθshapedを生成し、数値計算と実験により残留振動の大きさを調べた。yEについて、数値計算の結果を図44に、実験結果を図45に示す。これより、時変性を考慮しないZVShaperと比べて、本発明の手法の方がより効果的に残留振動を抑制できることが分かる。最大角速度をπrads,最大角加速度をAmaxとし、2つの関節が同時に動作終了する最短時間制御を参照指令値θrefとする。Amaxを変化させて、この参照指令値をもとにZVShaperで生成した指令値θZVと、本発明の手法で生成した指令値θshapedを与えたとき、残留振動のアーム先端の最大変位を数値計算により調べた結果を図46に示し、実験により調べた結果を図47に示す。突き出し動作と比べると、1次モードに寄与するトルクが大きいために、実験では残留振動の振幅が全体的に大きくなってしまっている。しかし、時変性を考慮しないZVShaperでは動作時間により残留振動が抑制される場合とそうでない場合があるのに対し、本発明の手法では動作時間によらず残留振動が抑制されている。このことから、本発明の手法が多関節型ロボットアームに有効であると言える。折りたたみ動作上述の説明にて述べたとおり、本発明の手法では全てのモードの残留振動を抑制する。しかしながら、これまでに見てきた動作例では2次の振動モードがほとんど励起されなかった。そこで、2次モードの振動も抑制されることを示すため、図48に示すアームを折りたたむ動作により発生する残留振動を、実験によりレーザー変位計で測定した。測定した座標系ΣEにおけるアーム先端の位置を図49に示す。図49よりyE方向の変位には2次モードの成分が多く含まれていると考えられる。そこで、yEのパワースペクトルを図50に示す。図50より、2次モードについてもピークレベルが下がっており、本発明の手法により複数の振動モードが抑制されることが分かる。なお、図50より、残留振動を抑制し、偏差角が小さくなると、固有振動数が変化することが分かる。これは、波動歯車装置の非線形剛性や、サーボ機構によるものだと考えられる。本発明の手法では、このような剛性の非線形性を厳密に考慮することはできない。しかしながら、上述にて説明した3つの動作例による実験結果から、剛性が線形であるとしてアクセレランスなどをもとにパラメータ同定をすることである程度まで振動が抑えることができることが分かる。実施例3:垂直2関節型ロボットアーム次に、実施例3として、垂直2関節型ロボットアームへの適用した例について説明する。具体的には、図51に示す鉛直面内で動作する2関節ロボットアーム20について、残留振動を抑制するための加減速パターンを生成し、数値計算により重力がある場合でも本発明の手法が有効であることを検証した。運動方程式関節jの回転角度をθj,その指令値をθ0j,リンクlの質量をml,リンクlの質量中心まわりの慣性モーメントをJGl,リンクlの長さをbl,関節jからリンクlの質量中心までの距離をlj,関節jのねじり剛性をkjとする。運動方程式は数168のようになる。ただし、数169~数179の関係にある。また、g=98ms2は重力加速度であり、図51に示す方向に重力が働いているものとする。数5~数10によるテイラー展開により、参照指令値をθref=θref1θref2Tとすると、垂直2関節型ロボットアームの運動方程式は数21で表される。ただし、各パラメータは数11~数15より、数180~数183のようになる。なお、数184~数190の関係にある。ここで、減衰力は各関節の偏差角の変化率に比例して働くものとすると、粘性減衰行列Cvは、数193のように表される。なお、本実施例3では、図52の表に示す諸元を用いて以下の数値計算を行った。突き出し動作での加減速パターンの生成図24に示す突き出し動作について考える。重力による釣り合いの位置を考慮すれば動作開始時の指令値θ0S=θ01Sθ02STは、数194を満たすθ0Sとなる。ここで数195である。同様にして、動作終了時の指令角θ0E=θ01Eθ02ETは、数196を満たすθ0Eとなる。図24のように、θ1S=7π24=0916rad,θ2S=-7π12=-183rad,θ1E=θ2E=0radである場合には。θ01S=0926rad,θ02S=-183rad,θ01E=0016rad,θ02E=0004radとなる。このことを考慮し、最大角加速度を4πrads2,最大角速度をπradsとすると、最短時間制御は図53に示すような指令値となる。これを参照指令値θrefとし、対応する加減速パターン数197,速度パターン数198についても図53に示す。K1grefの時間微分を解析的に求めるのは困難であるので、数152と同様に数199、数200と定義する。数199、数200、および数152を用いて、数64により得られる加速度パターン数59と対応する数60、数61を図54に示す。ただし、ΔtP=0001s,ΔtG=105sとした。整形前の図53と整形後の図54とを比べることにより、両者の速度パターンと角度指令値は大きく違わないことが分かる。これにより、モード形状の時間推移は整形後も変わらないものとみなすことが妥当であると言える。結果と考察参照指令値θrefを与えた場合と比べて、整形した加減速パターン数59に対応する角度指令値θshapedを与えたときの方が、残留振動が抑制されることを数値計算により示した。θrefを与えた場合とθshapedを与えた場合について、それぞれDormandPrince法により数4を数値的に解いた。解いた結果得られるθから計算した、動作終了後のアーム先端の座標系ΣBにおけるy方向の変位を図55に示す。これより、本発明の手法が、重力を考慮しなければならない垂直2関節ロボットアームの残留振動抑制にも有効であることが分かる。また、数値計算の結果より、動作中の各関節の角速度θを図56に示し、釣り合いの位置からの偏差角ψ=ψ1ψ2Tを図57に示す。これらより、水平2関節型ロボットアームの場合と同様に、動作中の振動が抑制されているのが分かる。これらの数値計算の結果より、釣り合いの位置からの偏差角ψを数18と定義し、ψに関する運動方程式を数21のように書くことで、水平2関節型ロボットアームと同じように残留振動を抑制する加減速パターンを生成できることが分かる。本発明によれば、加速または減速運動に伴って運動体に生じるn自由度かつ時変の振動を、インプットシェイピング法を用いて抑制する場合に、n自由度かつ時変の振動における複数のモード間のエネルギの移り変わりのうち、モード間の相互作用によるものを無視し、固有振動数および減衰比の変化によるものを考慮して、第2の入力を決定するという手法を採用している。また、本発明では、n自由度かつ時変の振動におけるモード間の相互作用によるエネルギの移り変わりをゼロと見なして、固有振動数および減衰比の時間変化によるモード間のエネルギの移り変わりを考慮して、第1の入力に対する第2の入力の大きさを決定するという手法を採用していると言うこともできる。また、別の観点からは、本発明では、n自由度かつ時変の振動における複数のモード間の相互作用によるエネルギの移り変わりを無視し、加速または減速運動を行っている間に運動体自体の姿勢が変化することにより生じるコリオリがした仕事を考慮して、第1の入力に対する第2の入力の大きさを決定するという手法を採用していると言うこともできる。このような手法を採用していることにより、各振動モードにおいて、第1の入力に対応する入力を打ち消す第2の入力に対応する入力の大きさを個別に算出して、振動モードごとに算出された入力を合成することで第2の入力の大きさを決定することが可能となる。そのため、本発明によれば、第1の入力によって生じたn自由度かつ時変の振動第1の応答を抑制できる第2の入力を算出することが可能となる。よって、ロボットアーム、工作機械、およびクレーンなど、その動作中に姿勢変化を伴うような動作を行う運動体に生じる振動、すなわち、n自由度かつ時変の振動を、インプットシェイピング法を適用して効果的に抑制する振動抑制方法を提供することができる。なお、上記様々な実施形態のうちの任意の実施形態を適宜組み合わせることにより、それぞれの有する効果を奏するようにすることができる。"