複雑系の制御・設計論グループ - Department of Precision Engineering

複雑系の制御・設計論グループ
ヒトの神経筋骨格モデルに基づく歩行生成
工学研究科航空宇宙工学専攻
青井 伸也
Abstract: Animals including human beings generate various motions by cooperatively manipulating their
complicated and redundant musculoskeletal systems controlled by the nervous system. To elucidate the
mechanisms, this study carried out numerical simulations of human bipedal locomotion based on the planar
musculoskeletal model and nervous model from the physiological knowledge. The skeletal system consists
of seven rigid links that represent the trunk including the head and arms, thighs, shanks, and feet. Each leg
has nine principal muscles modeled as massless wires. The nervous model focuses on the information processing based on the feedback structure. The numerical simulations established stable bipedal locomotion
that compares human locomotion.
Key words: human bipedal locomotion, neuromusculoskeletal model, CPG, numerical simulation
1. はじめに
ヒトや動物は，その冗長で複雑な筋骨格系を巧みに且つ協調的に動かすことで適応的な歩行を
実現する．この運動は，脳，脊髄の神経回路網によって支配されている．これら神経回路網による
運動の生成と制御のメカニズムを明らかにするために，神経回路モデルと筋骨格系の力学モデル
を統合した神経筋骨格モデルに基づいたシステム論的な研究が始められている．特に，この神経
筋骨格モデルに基づいて数値シミュレーションを行い，そこから得られる結果を，姿勢運動，床
反力，EMG データなどの運動学・生理学的な実測データと比較検討することで，脳・脊髄神経回
路網の構成，作動原理を明らかにするとともに，脳・脊髄系の神経回路網における情報処理過程
の運動制御的な役割を明らかにすることが目的である．
我々はこれまでにニホンザルを対象として，解剖学的に詳細な骨格モデルと，神経生理学的な
知見に基づく神経回路モデルを構築し，数値シミュレーションを用いて研究を進めてきた [4]．そ
して更に，解剖学的に詳細なデータに基づいて筋系のモデル化を行い，より精密なシミュレーショ
ンを行うことを目標としている．本研究では，まずは解剖学的に詳細なデータが得られているヒ
トの歩行運動を対象として始めている研究について報告する．特にここでは，ヒトの詳細な 2 次
元筋骨格モデルと神経生理学的な知見から神経制御系のフィードバックに基づく情報処理に着目
して神経回路モデルを構築し，数値シミュレーションを行った結果について報告する．
2. ヒトの筋骨格モデル
(b)
ヒトの骨格モデルとして，頭部・上肢を (a)
含んだ体幹 1 節，それぞれ，もも・すね・足
の 3 節からなる左右両脚の合計 7 節から構
Trunk
成される 2 次元直鎖型剛体リンク系を用い
Gluteus
る (図 1(a))．この骨格モデルを形成する関節
Iliopsoas
Maximus
位置の配置とリンクの形状，質量・慣性モー
Rectus
Biceps
メントなどのパラメータは文献 [3] のデータ
Femoris
Thigh
Femoris Long
を用いる．この骨格モデルは，身長 164 cm，
Biceps
Vastus
Femoris Short
体重 50 kg であり，床はバネ，ダンパーを用
Shank
Gastrocnemius
いてモデル化する．
Soleus
筋は，質量を持たないワイヤーとしてモ
Tibialis
Foot
Anterior
デル化し，図 1(b) のようにそれぞれの脚に
9 本ずつ配置する．各筋は，収縮部と並列弾
図 1: Musculoskeletal model
性・粘性部からなり，発生する筋張力は文
献 [3] に従うものとして以下のように書ける．
CE
PD
PE
Fm = F̄m
· k(ξm ) · h(ηm ) · ym + Fm
+ Fm
k(ξm ) = 0.32 + 0.71 exp[−1.112(ξm − 1)] sin[3.722(ξm − 0.656)]
h(ηm ) = 1 + tanh(3.0ηm )
PD
Fm
PE
Fm
=
=
cPmD L̇m
PE
km
{exp[15(Lm
(1)
− L̄m )] − 1}
CE は収縮部の最大筋張力，y は
ここで，添え字 m は筋の種類を示し，Fm は発生する筋張力，F̄m
m
P
D
P
E
筋の活性度，Fm ，Fm はそれぞれ並列粘性部，弾性部で発生する筋張力である．k(ξm )，h(ηm )
はそれぞれ長さ，速度関係を示し，ξm ，ηm はそれぞれ筋長 Lm ，収縮速度 L̇m を規格化したもので
ある．ここで，筋長はそれぞれの筋の自然長，収縮速度は全ての筋において 3.0 m/s で規格化して
いる．ただし，筋長は股関節，膝関節，足首関節がそれぞれ直立状態に対して 0.0◦ ，29.8◦ ，15.7◦
のとき自然長となるように定めている．また，それぞれの筋の関節に対するモーメントアームは
関節角に関係なく一定としている．
3. 神経回路モデル
歩行中の筋の収縮を支配する運動ニューロンは，各種
感覚器からの入力を受けるとともに，歩行中枢パターン
supraspinal center
発生器 (CPG) からの入力を受ける．この CPG は，介在
CPG
ニューロンを介した階層的な構造を持ち，リズム発生部
(RG) とパターン形成部 (PF) から構成されていると推測
Rhythm Generator
されている [1, 5]．RG は基本となるリズムを生成し，各
種感覚器からの入力に基づいて修正する．PF はこのリ
Pattern formation
ズムをもとに，フィードフォワード (FF) やフィードバッ
feedfoward
feedback
ク (FB) に基づく情報処理を介して，運動ニューロンへ
の入力の時空パターンを形成する．FF による情報処理
では，大脳皮質，小脳，歩行誘発野等からの入力に基づ
musculo-skeletal structure
き，FB による情報処理では，上位中枢からの入力と体
性感覚器からの入力に基づいてパターンが形成される．
PF は固定された回路ではなく，多くの介在ニューロンが 図 2: Schematic of CPG model (bold lines
RG の信号により再構成され，FB あるいは FF に基づく indicate focused pathways in this paper)
情報処理により，多様な時空パターンを形成していると
推測されている．
これらの生理学的な知見をもとに，歩行運動の CPG の機能ブロック図を図 2 に示す．ただし，
太線は今回焦点を当てて考察する経路である．PF においては，FB および FF に基づく情報処理と
して，それぞれ個別の数学モデルを対応させている．本研究では，この歩行運動の CPG の機能ブ
ロック図に対して具体的に数理モデルを構成する．CPG は各関節の伸筋，屈筋の組に対して存在
すると推測されている [1]．ここでは，これらを全てまとめて，左右の脚それぞれに対して 1 つの
CPG を設定する．RG の基本的な機能は基準となるリズムの生成であり，その振幅，位相が中枢
および末梢からの入力によって柔軟に変化する．この RG の数理モデルとして，ここでは単純な
非線形振動子を用い，その位相を φ とする．
フィードバックに基づく情報処理のモデルとして，Feldman [2] は α − γ 連関を含む脊髄反射系
に基礎をおいた λ モデルを提案している．この λ モデルでは，緊張性伸張反射の閾値 (λ) を調節す
ることで，筋長を制御している．本研究では，この λ モデルに基づく数理モデルを構成する．ま
ず，RG のモデルとして用いる振動子の位相 φ に基づいて，基準となる λm (φ) を構成する．本研究
[m]
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
-0.2
0
1
2
3
4
5
6 [m]
図 3: Stick diagram of stable bipedal locomotion
では，この λm (φ) は Winter [7] の歩行データをもとにして構成する．具体的には，歩行中のキネマ
ティクスデータからそれぞれの関節の運動を構成し，それぞれの筋の筋長を求め，これを λm (φ)
として用いる．ただし，本研究で用いた歩行キネマティクスのデューティー比は 0.61 であり，歩
行周期は 1.0 s である．そして，この構成した λm (φ) を用いて入力を形成する．その際，それぞれ
の筋の活性度 ym は次から形成される．
τ ẏm = −ym + κ(Lm − λm (φ)) + σ L̇m
(2)
ただし，ym < 0 のときは ym = 0 とする．ここで，τ は時定数，κ，σ はゲインパラメータであり，
全ての筋で同じ値を用いる．また上述したように，RG は各種感覚器からの入力を受け取り，リズ
ムを修正する．本研究では，足の接地感覚情報に着目し，遊脚運動時に足が接地して床反力を受
けると，対応する脚の振動子の位相を支持脚相の開始位置にリセットすることで修正する．
4. シミュレーション結果
Time constant τ [ms]
上記の筋骨格モデルと神経回路モデルに基
100
づいて，2 足歩行運動の数値シミュレーショ
ンを行った．その結果，転倒することなく定
80
with a phase reset
常で安定な歩行運動が実現された．その様子
60
を図 3 に示す．ただし，胴体と右脚は実線，
40
左脚は破線で表示している．
without a phase reset
また，足の接地感覚情報に基づいた位相リ
20
セットからリズムの調整を行う場合と，行わ
0
1000
2000
3000
4000
5000
ない場合についての比較を行った．具体的に
Gain parameter κ
は，λ モデルのゲインパラメータが一定とい
う条件の下，筋活性度の時定数を徐々に上げ 図 4: Comparison between with and without a phase
ていき，転倒したときの時定数を記録するこ reset
とで比較を行った．その結果を図 4 に示す．
ここで，ゲインパラメータに対して τ = 0 を示している結果は安定な歩行が全く実現されなかっ
たことを示している．位相リセットを用いてリズムの調整を行う方が、より時定数の大きいとこ
ろでも、転倒することなく歩行が持続されている様子が確認できる。すなわちこの結果から，足
の接地感覚情報に基づいてリズムの調整を行うことで，より安定な歩行が実現されることが確認
できる．
5. おわりに
ヒトや動物は，複雑な筋骨格系を巧みに協調させ，適応的な歩行を実現している．こうした動
物の優れた歩行生成知能の解明に向けて，我々は数理モデルを構築し，数値シミュレーションを
通して考察を行っている．本研究では，特にヒトの歩行運動を対象に，2 次元骨格モデル，筋モデ
ル，神経回路モデルを構築して数値シミュレーションを行い，安定な歩行を実現した．ここで神
経回路モデルは，特にフィードバックに基づく情報処理に着目して構築している．その際，この
数理モデルに基づく数値シミュレーションでは，ヒトの運動学データから再構成した関節角の運
動がほぼ常に実現されている状況で，直接筋駆動系で制御されない姿勢の運動，すなわち接地点
まわりの運動をいかにして生成することで安定な歩行を実現するかという問題になっている．こ
のシミュレーションの結果として明らかとされたことは，関節角としての内部状態の局所的な安
定性が満たされている状況の下，姿勢運動という外部状態の安定性を実現するためには，接地感
覚器を通して与えられる離散的なイベントとしての接地情報に基づいて，内部状態の位相関係な
いし周波数を調節することが重要であるということである．実際これを利用することで，RG の生
成するリズムの基本となる周波数を変化させ，歩行速度を変えたとしても，広い範囲で安定な歩
行を実現することもできた．また，実際生理学的な研究からも，このような接地感覚情報に基づ
く運動調整の重要性が示唆されている．
しかしながら，生理学的な研究から明らかとされているところでは，この神経回路モデルで用
いているような FB に基づく情報処理が，歩行中常に働いているとは限らない．むしろ，接地のよ
うな局所的な所で非常に強く働き，あとは FF に基づく情報処理で運動が生成されていると考えら
れる．特に，この FF に基づく情報処理として，内部モデルに基づいたモデル予測制御の枠組みが
広く知られており [6]，現在 FB に基づく情報処理系だけでなく，FF に基づく情報処理系を統合し
た神経回路モデルの作成に取り組んでいる．
参考文献
[1] R.E. Burke, A.M. Degtyarenko, and E.S. Simon: ”Patterns of locomotor drive to motoneurons and
last-order interneurons: Clues to the structure of the CPG”, J. Neurophysiol. 86(1): 447–462, 2001.
[2] A.G. Feldman: ”Once more on the equilibrium-point hypothesis (λ-model) for motor control”, J.
Mot. Behav. 18(1): 17–54, 1986.
[3] N. Ogihara and N. Yamazaki: ”Generation of human bipedal locomotion by a bio-mimetic neuromusculo-skeletal model”, Biol. Cybern., 84: 1–11, 2001.
[4] N. Ogihara, M. Nakatsukasa, Y. Sugimoto, S. Aoi, and K. Tsuchiya: ”Adaptive locomotion mechanisms inherent in the musculoskeletal structure”, Proc. SICE-ICASE Int. Conf.,, 2006.
[5] I.A. Rybak, K. Stecina, N.A. Shevtsova, and D.A. McCrea: ”Modelling spinal circuitry involved
in locomotor pattern generation: insights from the effects of afferent stimulation”, J. Physiol.
577(2):641–658, 2006.
[6] D.M. Wolpert and Z. Ghahramani: ”Computational principles of movement neuroscience”, Nat.
Neurosci. 3: 1212–1217, 2000.
[7] D.A. Winter: ”Biomechanics and motor control of human movement (3rd edition)”, New York:
John Wiley & Sons, 2004.
受動的歩行の安定解析と応用に関する研究
工学研究科航空宇宙工学専攻
杉本 靖博
Abstract: The purpose of this research is to analyze the stability of Passive Dynamic Walking (PDW)
using a linearized analytical Poincaré map. In particular, we focus on a bifurcation phenomenon in PDW.
Although the bifurcation of the walking period is one of the well-known features of PDW, it have not
been studied sufficiently so far. Using techniques similar to our previous research, we derive an analytical
Poincaré map for 2-period walking and discuss the stability of PDW with this map. In addition, we point
out that there is a similar interesting structure in this Poincaré map.
Key words: Passive Dynamic Walking, Poincaré map, stability analysis, bifurcation
1. はじめに
近年歩行ロボットの研究において McGeer により提唱された Passive Dynamic Walking(PDW: 受
動的歩行) が着目されている [1]．この受動的歩行には，いくつかの興味深い性質が知られており，
これらの特徴をうまく利用することでより高度なロボットの歩行を達成する一助になると考えら
れ，McGeer 以来，受動的歩行に関する様々な研究がなされてきた．
特に，受動的動歩行を周期運動と考え，そのリミットサイクルを用いた安定性の議論もなされ
ている [2, 3, 4]．しかしながら，対象とするロボットが非常に簡単なものであっても，状態のジャ
ンプを含むある種のハイブリットシステムとなるため，これまでの研究ではリミットサイクルや
それに伴うポアンカレマップが解析的な扱いをされることはほとんどなく，数値的な安定性の議
論に終始していた．その結果，受動的歩行が安定であることは示されてきたものの，何故安定に
なるのか，そのメカニズムに関しては明らかにされてはいなかった．最近になり，ようやく，対象
を非常に簡単なモデルと限定したり，線形化したロボットの運動方程式を用いることで，状態の
ジャンプを有するシステムとなる受動的歩行の解析的な安定解析を行なうことがなされるように
なり，ポアンカレマップの平衡点が何故安定になるのか，つまりは，
「何故受動的歩行が安定にな
るのか」を解明しようとする研究も報告されるようになってきた [5, 6]．
ところで，受動的歩行においては，斜面の傾斜角やロボットの物理パラメータを変化させると
歩行周期が二周期，四周期と多周期となっていく「分岐現象」という興味深い現象が起こること
も知られている [2, 3, 7]．この分岐現象は，単に数値シミュレーションだけでなく，実機にて実際
にその現象が観測されてもいる [7]．しかしながら，受動的歩行の安定化メカニズム同様，そうし
た現象が起きることはよく知られているものの，何故そうした現象が起きるのかについては明ら
かにされていない．
そこで本研究では，[6] で求められている解析的なポアンカレマップと同様の手法を用いること
で，二周期歩行を行っている際のポアンカレマップを導出し，それを用いた受動的歩行の安定解
析，特に分岐が起きる場合での受動的歩行の安定解析を行った．
2. 歩行ロボットモデル
図 1 に本研究で対象とする歩行ロボットのモデルを示す．
ロボットの運動を記述する基礎式として，以下の 3 つを用いる．詳細については [6] 参照．
まず，ロボットの運動方程式としては，Lagrange 法から導いた，
M (θ)θ̈ + N (θ, θ̇)θ̇ + g(θ, α) = u(t)
を，線形化して得られる
ẋ = Ax + b
を用いる．ただし，[θ, θ̇]T = x．
（今回は入力を考えないため，u(t) = [0, 0]T とした．
）
(1)
次に，遊脚着地時における支持脚と遊脚の速度の切り
替え方程式は，いくつかの仮定の元，角運動量保存則か
ら導くことができ，
"
#
Ta
0
+
x (t) =
· x− (t) := R(x− (t))
(2)
0 Tr Pa−1 Pb
とできる．
最後に，遊脚と地面の衝突は，ロボットの状態が状態空
間内のある幾何学的条件を表わす平面に到達した際に生
じると考えることができる．この条件は，定数行列 C =
[1 1 0 0] を用いて，
mb
l
rg
θp
−θw m
α
図 1: A compass model of walking robot
Cx− = 0
(3)
と書くことができる．
3. 衝突点に関するポアンカレマップ
3.1 一周期歩行でのポアンカレマップ 受動的歩行を行なっている際の衝突点及び歩行周期を
−
−
x−
∗ , τ∗ とし，k 回目の衝突時の状態偏差を ∆xk := xk − x∗ とおくと，k 回目の衝突時の状態偏差
∆xk を ∆xk+1 へと移すポアンカレマップ Pk :
∆xk+1 = Pk (∆xk )
(4)
¯
k¯
を定義することができる．すると，[6] で述べられているように，Pk の線形一次近似 Pk := ∂P
∂x ¯
は，
Pk = SeAτ∗ Rd , S = I −
v∗ C
, v∗ = Ax−
∗ + b.
Cv∗
x=x−
∗
(5)
と求めることができる，
3.2 二周期歩行を行っている際の衝突点に関するポアンカレマップ ポアンカレマップ (5) は，一
周期歩行を行っている際の衝突点 x−
∗ 及び歩行周期 τ∗ 近傍でのシステムの挙動を表しているため，
二周期歩行を行っている際のシステムの安定性を考察するには適切なものではないと考えられる．
そこで，一歩行周期での解析的なポアンカレマップ (5) を求めたときと同様に，二周期歩行を
行っている際のポアンカレマップ:
∆xk+1 = P̃k−1 (∆xk−1 )
(6)
を求め，その線形一時近似 P̃k−1 を用いて，二周期歩行をしている際の受動的歩行の安定性の考察
を行う．
二周期歩行を行っている際には，受動的歩行を行なっている際の衝突点及び歩行周期は二つ存
−
在するので，それらをそれぞれ x−
∗ , x∗∗ 及び τ∗ , τ∗∗ とする．さらに k − 1 回目及び k 回目の衝突時
−
−
−
の状態偏差を ∆xk := x−
k−1 − x∗ 及び ∆xk+1 := xk − x∗∗ とすると，k − 1 回目の衝突時の状態偏
差 ∆xk−1 を ∆xk+1 へと移すポアンカレマップ P̃k−1 の線形一時近似 P̃k−1 は以下の様に求めるこ
とが出来る．
¯
た だ し ，R∗
=
∂R(x) ¯
∂x ¯x=x
∆xk+1 = P̃ · ∆xk−1
v∗ C Aτ∗
v∗∗ C Aτ∗∗
P̃ = (I −
)e R∗∗ (I −
)e
R∗
Cv∗
Cv∗∗
¯
∂R(x) ¯
,
R
=
, v∗∗ = Ax−
∗∗
∗∗ + b,v∗
∂x ¯
−
−
∗
x=x∗∗
(7)
(8)
=
Ax−
∗ + b
2
1.8
そこで，二周期歩行を行う際のポアンカレマップ (8)
を用いて，傾斜角 α を変化させた際の安定性の変化を考
察する．α を 3.0→6.7[deg] に変化させた際の P̃k の固有
値の絶対値の最大値をプロットしたものを図 3.2 に示す．
これによると，α が大きくなるにつれ，max(|λi (P̃k )|) が
増加し，α=6.6 付近で 1 より大きくなり不安定化してい
α[deg]
ることが確認できる．一方で，[8] によると，式 (5) で表
される一周期歩行のポアンカレマップ Pk を用いた場合
図 2: max(|λi (P̃k )|)
では，傾斜角が α=5.3[deg] より大きくなったあたりで，
max(|λi (P̃k )|) は 1 より大きくなり Pk は不安定化する．つまり，α=5.3-6.6[deg] の範囲では一周期
歩行としては不安定であるが，二周期歩行としては安定となっている．言い換えれば，この傾斜
角のもとで継続した安定な歩行を実現するためには，一周期歩行ではなく二周期歩行を行う必要
があったといえる．
1.6
1.4
|λ i |
1.2
1
0.8
0.6
0.4
0.2
0
3
3.5
4
4.5
5
5.5
6
6.5
Remark 一周期歩行のポアンカレマップ (5) 及び二周期歩行のポアンカレマップ (8) を比較する
と，四周期歩行のポアンカレマップは次のように求められると考えられる．
P̂ = P (4) · P (3) · P (2) · P (1)
P (1) = (I −
(9)
v(1) C Aτ(4)
v(2) C Aτ(1)
R(1) , P (2) = (I −
R(2) , ...
)e
)e
Cv(1)
Cv(2)
同様に 8 周期や 16 周期歩行のポアンカレマップも求めることが出来ると予想できる．しかしな
がら，より高次になると，ポアンカレマップの線形化時における誤差がより大きくなってしまう
ため，PDW の安定性の議論を行う際，こうした線形化したポアンカレマップを用いるのは徐々に
不適当になっていくと思われる．
4. ポアンカレマップ内のフィードバック構造
次に，二周期歩行に対応するポアンカレマップ (8) に
内在している構造について考察を行う．
[6] にて述べられているように，一周期歩行のポアンカ
レマップ (5) 内にはフィードバック構造が存在している
ことが知られている．
そこで，二周期歩行のポアンカレマップ (8) について
も同様の置き換えをしてやると，
P̃k−1 = A⋆P − BP⋆ KP⋆
∆xk+1 = A⋆P ∆xk−1 + BP⋆ ∆τ ⋆
∆τ
⋆
=
−KP⋆ ∆xk−1
∆xk+1
+
−
BP∗
1
z
−
∆xk−1
A∗∗
P
BP∗∗
KP∗∗
KP∗
⇓
∆xk+1
−
(12)
∆xk
A∗P
(10)
(11)
1
z
+
1
z2
∆xk−1
∗∗ ∗
A∗P (A∗∗
P − BP KP )
∗∗ ∗
と書くことができ，図 3 の様に P̃ 内にもフィードバック
BP∗
KP∗∗ (A∗∗
P − BP KP )
構造が存在していることがいえる．式 (10)- (12) の各要
素の詳細は [8] 参照．
図 3: Feedback structure in P̃
ここで，二つのフィードバック構造間の関係を考える．
前節での安定性の議論より，傾斜角が変化した際，歩行が継続していくためには，ポアンカレマッ
プ内の構造が一周期のものから二周期のものへと変化する必要があった，つまり，ロボットは二
周期歩行をする必要があったと考えられる．ここで，二周期歩行は一周期歩行を内包する，つま
り一周期歩行として安定なら二周期歩行としても安定であることに注意する．
さらに，前節の Remark で述べたように，4
周期，8 周期，…のポアンカレマップも求める
ことができ，かつ，これらのポアンカレマップ
にも同様のフィードバック構造が存在している
と考えられる．一周期と二周期との関係と同様
に，各倍周期歩行間の内包関係を考えれば，以
下が言えると考えられる [9]．
1
∆x
1
z
+
BP
∆x
+A
−
−
∆Ǽ
k
k+
P
1
1
Z
∗
BP
∆x
k−
14
z
∗2)
k+
x
∆
∗∗ KP
∗∗ −
∗ (A P
−
AP
BP
+
∗∗
P
∗ ∗ (A
KP
KP
2
k−
+
A∗P)
∗∗ KP
−
BP
+
KP
+
BP
.... n-period
....
1-period
2-period 4-period
• PDW のシステム全体には，図. 4 のような階
層的なフィードバック構造が存在している． 図 4: The hierarchical implicit structure in PDW
• 傾斜角やロボットの物理パラメータに応じ robot system
て，歩行を継続するために最も適切な構造
が自動的に選択される．
• 選択された構造に対応する歩容（一周期歩
行，二周期歩行，…）が現れる．
逆に，こうした興味深い構造が存在していることにより，環境などに対する PDW の適応的な歩
行が実現していると考えられる．
5. おわりに
二周期歩行に対応する解析的なポアンカレマップを求め，受動的歩行に見られる興味深い現象
の一つである分岐現象の解析を行った．また，一周期歩行に対応するポアンカレマップと同様に
二周期歩行に対応するポアンカレマップに関しても，その内部にフィードバック構造が存在する
ことを示し，一周期歩行から二周期歩行に変化することで，フィードバック構造にどういった変
化が起こるかを示した．
参考文献
[1] T. McGeer : Passive Dynamic Walking, The Int. J. of Robotics Research, 9-2, (1990), 62-82.
[2] A. Goswami, B. Thuilot and B. Espau : Compass-Like Biped Robot Part I: Stability and Bifurcation
of Passive Ga its, Technical Report 2996 INRIA (1996).
[3] A. Goswami, B. Thuilot and B. Espiau : A Study of the Passive Gait of a Compass-Like Biped
Robot: Symmetry and Chaos, The Int. J. of Robotics Research, 17-12, (1998), 1282-1301.
[4] M. Garcia, A. Chatterjee, A. Ruina and M. Coleman : The Simplest Walking Model: Stability,
Complexity and Scaling, ASME J. of Biomechanical Engineering, 120-2, (1998), 281-288.
[5] 平田, 小亀：状態にジャンプを有する線形システムの周期運動 –Compass Walking のモデリン
グ, 安定解析, フィードバック制御, システム制御情報学会論文誌, 17-12, (2004), 553-560.
[6] 杉本, 大須賀：受動的動歩行の安定性に関する一考察–ポアンカレマップの構造解釈からのア
プローチ–, システム制御情報学会論文誌, 18-7, 255/260 (2005).
[7] 大須賀, 桐原：受動的歩行ロボット Quartet II の歩行解析と歩行実験, 日本ロボット学会誌, 18-5,
(2000), 737-742.
[8] 杉本, 大須賀：解析的なポアンカレマップを用いた受動的歩行における分岐現象の解析, 第 11
回ロボティクスシンポジア予稿集, 2006, 92-97.
[9] Y. Sugimoto, K. Osuka: ”Hierarchical Implicit Feedback Structure in Passive Dynamic Walking”,
Proc. of 2007 IEEE/RSJ International Conference on Intelligent Robots and Systems(IROS 2007).
Adaptive Output Regulation of Affine Nonlinear Systems
工学研究科航空宇宙工学専攻
Mai Bando
Abstract: In this paper the output regulation problem for multilinear systems with unknown parameters is
considered. Based on the Lyapunov stability theory, an adaptive controller which stabilize the system is
derived. Then sufficient conditions for the output regulation problem with full information to be solvable
are established. Simulation results are given to illustrate the theory.
Key words: Output Regulation, Adaptive Control, Nonlinear Control, Tracking, Distrubance Rejection
1.
Introduction
Consider a affine nonlinear system of the form
ẋ = A(x)x + B1 (x)w + B2 (x)u,
z = C1 x + D11 w + D12 u,
(1)
where x ∈ Rn is the state, u ∈ Rm is the control, and z ∈ Rq is the output to be regulated. The signal
w ∈ Rs denotes disturbances or reference signals generated by an anti-stable exosystem
ẇ = Sw.
(2)
The output regulation problem is to find a control law such that z(t) converges to zero as t → ∞ for any
initial conditions of the plant and the exosystem. When the state and the signal are available, it is called
the output regulation with full information and a control of the form u = F x + Gw is sought. If A, B1
and B2 are constant matrices, the following results are known [1, 2, 3].
Theorem 1. Suppose (A, B2 ) is stabilizable. Then the output regulation problem with full information
is solvable if and only if there exist two matrices Π and Γ which satisfy the regulator equation
AΠ − ΠS + B1 + B2 Γ = 0,
C1 Π + D11 + D12 Γ = 0.
(3)
Under this condition admissible controllers are given by
u(t) = F x(t) + (Γ − F Π)w(t),
where F is any matrix such that A + BF is exponentially stable.
Output regulation for the affine nonlinear system with special structure is considered in [4, 5]. In
this paper, we consider the output regulation problem for affine nonlinear systems and then consider the
adaptive output regulation. First we consider the output regulation for exactly known affine nonlienar
system (1). Next we consider affine nonlinear system described by multiple linear models with unknown
parameters which is special case of (1). We introduce an estimator of the system and adaptive laws which
decrease the estimation error of the state and parameters. We then introduce a feedback law based on the
algebraic Riccati equation. We give a sufficient condition for the local asymptotic stability of the closedloop system. Based on this result, we consider the local output regulation problem with full information
and give sufficient conditions for solvability in terms of the regulator equation. The paper is organized
as follows. Section 1 is introduction. Section 3 is concerned with output regulation of affine nonlinear
systems. As an example, output regulation of the van der Pol equation is considered and simulation
results are given. Finally in Section 4, adaptive output regulation of nonlinear systems described by
multiple linear models is considered. As an example, adaptive tracking of a step and sine function by a
nonlinear mass-spring system is considered, and simulation results are given.
2.
Output Regulation of Affine Nonlinear Systems
We shall consider the output regulation problem for (1) under the following condition:
Assumption 1. D12 = 0 and (Ai , B2i ) is controllable canonical form of the same structure each i and
C1 has the structure
£
¤
C1 = C 1 C 2 · · · C m ,
(4)


0
 ..

 .



i

 ∈ Rq×ni ,
C = 1
(5)

 .

 ..

0
0
where the i-th elements of the first column of C i are one.
2.1 Stabilization Consider first the regulation
problem of (1) with w = 0. Let Q be positive-definite.
√
Then (C, A(x)) is observable, where C = Q. Since (A(x), B2 (x)) is by Assumption 1 stabilizable for
all x, there exists a positive stabilizing solution X = X(x) of the algebraic Riccati equation
A(x)T X + XA(x) + Q − XB2 (x)B2 (x)T X = 0.
(6)
This is known as State-Dependent Riccati equation (SDRE) [6, 7]. We can show the regularity of the
solution X in x = (xi ).
Lemma 2. For all x, there exists a unique positive stabilizing solution X = X(x) of (6) which is
continuous and continuously differentiable with respect to x. Moreover, if x stays in a compact domain,
∂X
is bounded for any i.
∂xi
Now we introduce the control law
u = −B2T X(x)x,
(7)
and we shall examine the stability of the closed-loop system
ẋ = [A(x) − B2 (x)B2 (x)T X(x)]x.
(8)
Theorem 3. The equilibrium xe = 0 of the system (1) is locally asymptotically stable.
2.2 Output regulation Consider the output regulation problem associated with (1). By Assumption
1, there exists a solution (Π, Γ) of the regulator equation
A(x)Π − ΠS + B1 + B2 (x)Γ = 0,
C1 Π + D11 + D12 Γ = 0,
(9)
such that Π is independent of x [8]. We choose the controller
u = −B2T (x)X(x)x + (Γ(x) + B2T (x)X(x)Π)w,
(10)
where X(x) is the solution of the Riccati equation (6).
Theorem 4. Under Assumption 1, the local output regulation is fulfilled, i.e., lim z(t) = 0 for any small
t→∞
x0 and w0 .
We shall apply Theorem 4 to the van der Pol equation.
2.3 Example
2.5
x
u(t)=0
2
3
w(1)
x(1)
u(t)=0
2
1.5
1
0.5
x(2)
x(1),w(1)
1
0
-0.5
0
-1
-1.5
-1
-2
-2.5
-2.5 -2 -1.5 -1 -0.5
-2
0
5
10
15
t
20
25
30
1
1.5
2
2.5
Figure 2: The trajectories of the state in the phase
plane.
Figure 1: The trajectories of the state.
Example 2.3.1
0 0.5
x(1)
Consider the van der Pol equation
ξ¨ = 0.2(1 − ξ 2 )ξ˙ − ξ + u.
(11)
This system can be represented by the following affine nonlinear system:
ẋ = A(x)x + Bu,
where
·
A(x) =
0
1
−1 0.2(1 − x21 )
¸
·
, B=
(12)
0
1
¸
·
, x=
x1
x2
¸
·
=
ξ
ξ˙
¸
.
We design a state feedback controller such that x1 (t) → 0.3sin(t). In this case we set C1 = [1 0],
D11 = [−1 0] and take the following exosystem
·
¸
0 1
S=
.
−1 0
Then Π = I2 and Γ = [−1, −0.2(1 − x21 )]. The free system (11) has a periodic solution with amplitude
= 2. The simulation results with x(0) = [2 0]T , w(0) = [0 1]T are shown in Figures 1 and 2. The
dotted lines are trajectories of the free system with u = 0, and the solid lines are controlled trajectories.
3.
Nonlinear Systems Described by Multiple Linear Models
Consider the nonlinear system described by multiple linear models
ẋ =
r
X
λi (x)Ai x + B1 w +
i=1
r
X
λi (x)B2i u,
i=1
(13)
z = C1 x + D11 w + D12 u,
where the constant matrices Ai and B2i contain unknown parameters, and the other matrices are assumed
to be known. We assume Ai , B2i and C1 satisfy Assumption 1. Note that all pairs (Ai , B2i ) are in
r
r
X
X
the same controllable canonical form and hence (
λi (x)Ai ,
λi (x)B2i ) is also in the controllable
i=1
i=1
canonical form for any unknown parameters of (Ai , B2i ) and for any x. λi (x) are known functions of x
such that
r
X
λi (x) = 1, λi (x) ≥ 0.
(14)
i=1
Note that (13) is a special case of (1) if Ai and B2i are known. Hence Theorems 3 and 4 can be applied.
3.1 Adaptive regulation First we consider the stabilization problem. Following [9] we introduce an
estimator and adaptive laws of the form
x̂˙ = Am x̂ +
r
X
λi (x)(Âi − Am )x + B1 w +
i=1
r
X
λi (x)B̂2i u,
(15)
i=1
˙
Âi = Φ̇i = −λi (x)P exT ,
˙
B̂2i = Ψ̇2i = −λi (x)P euT ,
(16)
where Am is an n × n stable matrix, P is the solution of the following matrix equation
ATm P + P Am = −Q0
for some positive-definite matrix Q0 and
e = x̂ − x, Φi = Ai − A, Ψ2i = B̂2i − B2 .
Then the error equation is given by
ė = Am e +
r
X
λi (x)(Φi x + Ψ2i u).
(17)
i=1
If some elements of Ai and B2i are known, we can omit their adaptive laws in (16), but for notational
convenience we use (16).
Lemma 5. If x(t) and u(t) are bounded for all t ≥ 0, then (16) and (17) are globally stable.
P
Consider (6) with (A(x), B2 (x)) replaced by (Â(x, t), B̂2 (x, t)), where Â(x, t) = ri=1 λi (x)Âi (t),
P
B̂2 (x, t) = ri=1 λi (x)B̂2i (t). Let Q be positive-definite. Then for each t fixed, (6) becomes
r
r
r
r
X
X
X
X
(
λi (x)Âi )T X + X(
λi (x)Âi ) + Q − X(
λi (x)B̂2i )(
λi (x)B̂2i )T X = 0.
i=1
i=1
i=1
(18)
i=1
X = X(x, θ(t)) is continuously differentiable function of x and θ(t), where
θ = [vec Â1 , · · · , vec Âr , vec B̂21 , · · · , vec B̂2r ]T .
Now introduce the control law
u = −B̂2T X(x, θ(t))x̂
(19)
and consider the stability of the adaptive control system.
Theorem 6. For sufficiently small x(0), x̂(0), e(0) and θ(0), x, x̂, e and θ are bounded and lim e(t) =
0. Moreover, if w = 0, then lim x(t) = 0.
t→∞
t→∞
3.2 Adaptive output regulation Now we consider the adaptive output regulation problem associated
with (13), (15) and (16). In this case regulator equation is given by
Â(x, θ)Π − ΠS + B1 + B̂2 (x, θ)Γ = 0,
C1 Π + D11 + D12 Γ = 0.
(20)
By Assumption 1, there exists a solution (Π, Γ) of (20) such that Π is a constant matrix. We choose the
controller
u = −B̂2T X(x, θ)x̂ + (Γ(x, θ) + B̂2T X(x, θ)Π)w,
(21)
where X(x, θ) is the solution of the Riccati equation (18) corresponding to (13), (15) and (16). Then the
following result is obtained.
Theorem 7. For sufficiently small x(0), x̂(0), e(0) and θ(0), the adaptive local output regulation is
fulfilled i.e., lim z(t) = 0.
t→∞
Remark We can extend Theorem 6 and 7 to the system (13) where λi (x) are continuously differentiable functions of x but does not satisfy (14). For example, van der Pol equation (11) can be written
as the three linear models:
r
´
³
X
(22)
ẋ = A0 +
λi (x)Ai x + Bu,
i=1
where
·
¸
·
¸
·
¸
· ¸
0 1
0 1
0 1
0
A0 =
, A1 =
, A2 =
, B=
,
0 0
a1 0
0 a2
1
·
¸ · ¸
ξ
x1
x=
= ˙ , λ1 (x) = 1, λ2 (x) = 1 − x21 , a1 = −1, a2 = 0.2.
x2
ξ
By replacing  =
r
³
´
X
λ
(x)
Â
in
(15),
(18)
and
(20)
by
A
+
λ
(x)
Â
i
0
i
i , Theorem 6 and 7 hold
i=1 i
Pr
i=1
for this system.
3.3 Examples
Example 3.3.1 Consider the nonlinear mass-spring system
ξ¨ = −0.01ξ − 0.67ξ 3 + u,
(23)
The nonlinear term satisfies the following conditions for ξ ∈ [ − 1 1]:
−0.67ξ 6 −0.67ξ 3 6 0ξ,
0ξ 6 −0.67ξ 3 6 −0.67ξ,
ξ > 0,
ξ 6 0.
Hence it can be represented by the following two linear models:
ẋ =
2
X
λi (x)(Ai x + Bi u),
(24)
i=1
where
·
A1 =
0 1
a1 a2
¸
·
, A2 =
0 1
a3 a4
¸
·
,Bi =
0
1
¸
·
, x=
x1
x2
¸
·
=
ξ
ξ˙
¸
,
λ1 (x) = 1 − x21 , λ2 (x) = x21 , a1 = −0.01, a2 = 0, a3 = −0.68, a4 = 0.
Here a1 and a3 are regarded unknown and λi are given functions. The stabilization by (19) is considered.
The simulation result with x(0) = [0.8 0]T , [â1 (0), â3 (0)] = [0, 0] is shown in Figure 3.
Example 3.3.2 We design a state feedback controller for (23) such that x1 (t) → 0.3sin(t). In this case
we set C1 = [1 0], D11 = [−1 0] and take the following exosystem
·
¸
0 1
S=
.
−1 0
Then Π = I2 and Γ = [−1 + (x21 − 1)â1 − x21 â3 , (x21 − 1)â2 − x21 â4 ]. The simulation result with
x(0) = [0.8 0]T , w(0) = [0 0.3]T , [â1 (0), â3 (0)] = [0, 0] is shown in Figures 4.
Example 3.3.3 For the van der Pol equation (11), , we design a state feedback controller such that
x1 (t) → sin(t) using the above remark. Here a1 and a2 are regarded unknown and λi are given functions. The simulation results with x(0) = [2 0]T , w(0) = [0 1]T and [â1 (0), â2 (0)] = [0, 0] are shown
in Figures 5 and 6.
0.8
1
x
x
x
x
0.6
0.6
0.2
0.4
x
x(1),w(1)
0.4
0
0.2
-0.2
0
-0.4
-0.2
-0.6
0
5
10
15
t
20
25
x(1)
w(1)
z
0.8
30
-0.4
0
Figure 3: Stabilization.
5
10
15
t
1.5
1
0.5
x(2)
1
x(1),w(1)
30
x
u(t)=0
2
w(1)
x(1)
u(t)=0
2
25
Figure 4: Sine tracking.
2.5
3
20
0
0
-0.5
-1
-1.5
-1
-2
-2
-2.5
-2.5 -2 -1.5 -1 -0.5
0
5
10
15
t
20
Figure 5: Sine tracking.
25
30
0 0.5
x(1)
1
1.5
2
2.5
Figure 6: The trajectories of the state in the phase
plane.
References
[1] A. Saberi, A. A. Stoorvogel and P. Sannuti, Control of Linear Systems with Regulation and Input
Constraints, Springer, 2000.
[2] B. A. Francis, The linear multivariable regulator problems, SIAM J. Contr. Optimiz. 15 (1974),
486-505.
[3] B. A. Francis and W. M. Wonham, The internal model principle for linear multivariable regulators,
Appl. Math. Optim. 2 (1975), 170–194.
[4] , Isidori, A. and Byrnes, CI., Output regulation of nonlinear systems, Automatic Control, IEEE
Transactions on, 35, 131–140, 1990.
[5] , Serrani, A. and Isidori, A. and Marconi, L., Semi-global nonlinear output regulation with adaptive
internalmodel, Control, IEEE Transactions on, 46, 1178–1194, 2001.
[6] Cloutier, J.R. and D’Souza, C.N. and Mracek, C.P., Nonlinear regulation and nonlinear H-infinity
control via the state-dependent Riccati equation technique. I- Theory, II- Examples, International
Conference on Nonlinear Problems in Aviation and Aerospace, 1 st, Daytona Beach, FL, 117–141,
1996.
[7] Mracek, C.P. and Cloutier, J.R., Control designs for the nonlinear benchmark problem via the statedependent Riccati equation method, International Journal of Robust and Nonlinear Control 45,
401–433,1998 117–141, 1996.
[8] M. Bando and A. Ichikawa, Adaptive Output Regulation of Nonlinear Systems described by Multiple Linear Models, Proceedings of 9th IFAC Workshop on Adaptation and Learning in Control and
Signal Processing 2007 (ALCOSP’2007), 2007.
[9] K. S. Narendra, A. M. Annaswamy, Stable Adaptive Systems, Prentice Hall, 1989.
自律型ロボットのための自己位置推定に関する研究
工学研究科機械理工学専攻
金田 さやか, 中西 弘明
Abstract: Precise localization is important for an autonomous robot, such as a rover for small planetary
bodies including asteroids and comets. Conventional localization methods are not adequate on small body
surface because of the planet’s irregular shape and small mass. In this paper, we propose an effective localization method of a rover for a small planetary body. It uses two-way range measurement between a rover
and a mother spacecraft. The measurements are conducted repeatedly. The effectiveness of this method is
certified by numerical simulations assuming ITOKAWA-size-planet(with less than 1[km] radius):the accuracy of the rover position is 1[m] by 6 hours observation.
Key words: Localization, Two-Way Range, Radio Wave, Autonomous Robot
1. はじめに
ロボットに関する研究において自己位置同定は重要である．自己位置はロボットの位置制御に
不可欠であるのは自明である．近年，ロボットの環境適応性に関する関心が高いが，環境の全て
をセンシングすることはできないことから，自己の振る舞いを通じて環境を認識する必要がある．
つまり，ロボットの自己状態同定は環境適応性においても中心的役割を果たす．自己状態の中でも
位置は利用価値が高く，ロボットの自己位置同定法については現在も極めて関心が高い．様々な
研究がなされている．本稿では，小天体探査用ローバを取りあげ，その移動制御のための自己位
置同定システムについて述べる．
ローバによる小天体表面探査は世界的にもいまだ実現していないが, 月面上および火星表面上の
探査機はこれまでにも幾つか実現例があった. しかし，過去の惑星探査ローバは, 探査範囲がラン
ダの近傍であって, 惑星との相対的な距離で考えるとごく局所的な範囲の移動であった. 一方，小
天体探査ローバは直径数百メートルの端から端まで，全ての表面を移動することが予想される．さ
らに，ローバの移動目標の候補の中で最も理学的興味の高い表面試料採取時の衝突痕は，大きく
とも直径数メートルであろう．直径数百メートルの天体上で，数メートルサイズの目標地への誘
導を行うには，リアルタイムかつメートルオーダの位置同定精度が要求される．このような要求
を満たす既存の位置同定手法は存在しない．画像マッチングによる位置同定は，ローバに搭載す
るカメラの分解能と周回衛星が作製した地図との精度が見合わないことが予想される．太陽や星
の方向を用いた位置同定は重力方向を基準とするが，小天体では表面で計測できる重力方向が必
ずしも天体重心を通るとは限らないために精度が見込めない．また，地上の GPS のように複数衛
星を小天体周回軌道に乗せるということも現実的ではない．
本稿では，小天体探査ローバとその母船間の相対距離を測定することにより, 位置同定を行う手
法を提案する．提案手法の位置同定精度および速度について，数値シミュレーションによって検
証を行い，小天体上でも適用可能であることを示す．
2. 電波の伝播遅延時間測定による位置同定手法の概要
提案手法では，位置同定のために小天体探査ローバとその母船間の距離を，電波の伝播遅延時
間によって測定する．その際，図 1 に示すように，母船から電波を送信し，次にローバで電波を
反射，再び母船で電波を受信，往復の距離 (two-way range) を測定することとする．母船側で測定
値を処理し，ローバ位置を同定する．ローバは母船との交信で自己位置を知り，目的地へ移動す
る．その際，以下の仮定をおく．
• ローバは位置同定中は小天体に対して静止している．ローバの運動は，小天体の自転運動のみ
に依存する．小天体では重力が小さく，脱出速度が低くなることから，ローバ移動速度の大き
さには大きな制約がある．このためローバが小天体に対して静止しているという仮定は十分な
mother spacecraft
rotation
axis
planet surface
X(ti,e)
X(t)
X(ti,ref)
radio signal
rover
x(t)
X(ti,r)
mother
spacecraft
X(ti+1,e)
X(ti+1,ref)
X(ti+1,r)
rover
図 1: Definition of Coordinates of the Proposed
Model
x(ti,e)
x(ti,ref)
x(ti,r)
x(ti+1,e)
x(ti+1,ref)
x(ti+1,r)
orb
it
図 2: Two-Way Range Measurement of the Proposed Model
近似といってよい．また，加速度計やオドメトリをローバに搭載することにより，加速度や対
地速度などの情報が得られるが，これらの情報を利用するように提案手法を拡張することは容
易である．
• 小天体の自転運動は既知である．小天体の重心，自転軸および自転角速度は既知であるとした．
通常の天体探査では，ローバによる表面探査に先んじて，母船によって天体の詳細な地図や重
力場地図が作製される．このため，一般に探査する天体の自転運動はローバ投入時には既知で
あり．適切な仮定である．
• 母船の軌道パラメータは既知である．従来の技術により，母船軌道は太陽方向や星配置，天体
のランドマークを基準として決めることができる．よって，母船の軌道パラメータは既知とし
てよい．
• 増幅回路などにおける遅延時間の統計的性質は既知である. 提案手法では，母船から発信され
た電波をローバにおいて反射する．小天体表面における反射波と区別するために増幅などをか
ける必要がある．このために回路遅延が発生するが，遅延時間の平均や分散はあらかじめ測定
できるため，この仮定は適切である．
提案手法は，母船において電波の往復伝播遅延測定し，異なる時刻における複数の測定結果か
らローバの位置同定を行う．異なる時刻の測定結果を用いるために単一の母船で実行可能であり，
往復距離から位置同定を行うアルゴリズムは処理計算が比較的簡単である．また，ローバ側に高
精度の時計と位置同定のための処理計算の負担がないことも利点である．
3. 電波の往復伝搬時間に基づく自己位置推定
3.1 座標系の定義 図 1 のように天体重心を原点とした慣性空間固定座標系をとり，天体の自転
軸を z 軸とする．母船の軌道面と z 軸に直交する平面との交線を x 軸とする．時刻 t におけるロー
バ，母船の位置をそれぞれ x(t)，X(t) と表す．
3.2 ローバの運動方程式
仮定よりローバの運動は小天体の自転運動のみに依存する．
x(t) = Rz (ψ(t − t0 )) x(t0 )
ただし，ψ は小天体の自転角速度を表す． Rz は z 軸まわりの回転行列である．
(1)
3.3 母船の運動方程式 母船は周回軌道上に存在すると仮定する．その運動について, 軌道 6 要
素 (昇交点経度 Ω, 軌道傾斜角 I, 軌道長半径 A, 軌道離心率 ε, 近地点引数 ω, 平均近点角 M0 ) を用
いて以下のように記述できる [1]．平均運動 n = (Gm/A3 )1/2 を用いれば，任意の平均近点角 Mt
は Mt = n(t − t0 ) + M0 と表される．ケプラーの方程式
M (t) = E(t) − ε sin E(t)
を満たす E(t) より，軌道上における母船位置 X obt (t) は
⎞ ⎛
⎞
⎛
A(cos
E(t)
−
ε)
Xobt (t)
⎟ ⎜ √
⎟
⎜
⎝ Yobt (t) ⎠ = ⎝ A 1 − ε2 sin E(t) ⎠
Zobt (t)
0
と表される．これを天体重心中心，慣性座標系で表すと母船位置 X(t) は
⎞
⎛
⎛
⎞
Xobt (t)
X(t)
⎟
⎜
⎜
⎟
⎝ Y (t) ⎠ = Rz (−Ω)Rx(−I)Rz (−ω) ⎝ Yobt (t) ⎠
Zobt (t)
Z(t)
(2)
(3)
(4)
である．ここで，Rx は x 軸周りの回転行列を表す．
3.4 観測方程式 母船と天体表面上ローバとの間の two-way range を測定する．電波の送信時刻
を te ，ローバにおける反射時刻を tref ，母船における受信時刻を tr とすると，観測量 τ は図 2 よ
り幾何学的に
1
τ = (|X(te ) − x(tref )| + |X(tr ) − x(tref )|)
(5)
c
と表される．ただし，c は光速を表す．以下，
1
· ρ(x(te ))
c
1
1
τd = · ρd = |X(te ) − x(tref )|
c
c
1
1
τu = · ρu = |X(tr ) − x(tref )|
c
c
τ=
(6)
(7)
(8)
と表す．
3.5 拡張カルマンフィルタの利用 時間間隔 Δt 毎に観測を行い, 電波送信時刻 te におけるロー
バ位置 xe を推定する．観測方程式が非線形のため, 拡張カルマンフィルタを用いる [2]. i は i 番目
の観測（時刻 t0 + iΔt），ξ は観測雑音，w は伝播雑音，∗ は真値，ˆ は推定量を表す．
観測量:
1
τ ob (i) = ρ(x∗ (te , i)) + ξ(i)
c
(9)
更新 (UpDate)：
1
x̂(i|i) = x̂(i|i − 1) + K(i)[τob (i) − ρ(x̂(i|i − 1))]
c
P (i|i) = P (i|i − 1) − K(i)H(i)P (i|i − 1)
(10)
伝播 (Propagation)：
x̂(i + 1|i) = Φx̂(i|i) + w(i)
(11)
P (i + 1|i) = ΦP (i|i)Φ + Q(i)
(12)
T
ただし，
−1
K(i) = P (i|i − 1)H(i)T H(i)P (i|i − 1)H(i)T + r
r(i) = E[ξ(i)2 ]
Q(i) = E (w(i) − w̄(i)) (w(i) − w̄(i))T
である. P は推定誤差の共分散を表す 3×3 行列. H は
∂ρ
H=
∂x(te )
(13)
(14)
(15)
で定義される 1×3 行列である．Φ は x の運動を与え，式 (1) より Φ = Rz (ψΔt) である．
3.5.1 観測量 τ の推定 i 番目の観測時における推定量 x̂(i|i− 1) = x(te ) より，観測量 τ (i+ 1|i) =
τd (i + 1|i) + τu (i + 1|i) を推定する．まず τ̂d を求めるためには，
cτ̂d = |X(te ) − Rz (τ̂d )x(te )|
(16)
cτ̂u = |X(te + τd + τ̂u ) − Rz (τd )x(te )|
(17)
を解く．得られた τ̂d を用いて
より τ̂u を求める．ここで，X(te + τd + τ̂u ) はケプラー方程式 (2) を解く必要があるが，(2) は非線
形方程式であるので逐次法によって求めなければならない．
3.5.2 ローバ位置による観測量変化 カルマンフィルタによる推定式 (10) における H は，式 (15)
より各電波送信時刻におけるローバ位置 x(te ) が微小変化した場合の，i 番目の観測量 τ の変化で
ある．すなわち，H = ∂τ /∂x(te ) = ∂τ /∂xe = ∂τd /∂xe + ∂τu /∂xe を求める．式 (16)，(17) より
1
∂τd
=
(xe − Rz (−ψτd )X e )T
∂xe
Aτd
∂xref
1
∂X r ∂τd
∂τu
T
·
=
(X r − xref )
−
∂xe
Aτu
∂τ ∂xe
∂xe
(18)
(19)
ただし，
Aτd = c2 τd + (Xe xe + Ye ye )ψ sin(ψτd ) + (Xe ye − Ye xe )ψ cos(ψτd )
∂X r
Aτu = c2 τu − (X r − xref )T
∂τ
∂X obt
∂X r
= Rz (−Ω)Rx (−I)Rz (−ω)
∂τ
⎞ ∂τ
⎛
∂E
A(− sin E) ∂τ
∂X obt ⎜ √
⎟
= ⎝ A 1 − ε2 cos E ∂E
∂τ ⎠
∂τ
0
Gm
∂E
A3
=
∂τ
1 − ε cos E(τ )
⎞T
⎛
−xe ψ sin(ψτd ) + ye ψ cos(ψτd )
∂xref
∂τd ⎜
⎟
= Rz (ψτd ) +
⎝ −xe ψ cos(ψτd ) − ue ψ sin(ψτd ) ⎠
∂xe
∂xe
0
(20)
(21)
(22)
(23)
(24)
(25)
であり, xref および X r に関しては, 伝搬時間の推定値 τ̂d , τ̂u を用いた x̂ref = x(te + τ̂d ), X̂ r =
X(te + τ̂d + τ̂u ) とした.
100
Δx 2 + Δy 2 + Δz 2
Δx
15
Variance Value[m2]
Estimation Error[m]
20
Δy
Δx + Δy + Δz
2
2
Δz
2
10
Δz
5
Pxx
Pzz
60
Pzz
40
Pyy
Pxx
Δx
1
Pyy
80
20
Δy
0
0
2
3
4
5
6
0
0
1
2
3
4
5
6
Observation Time[h]
Observation Time[h]
図 3: (left)Position Estimation Error and Observation Time:While x and y element errors are correctly determined fast, the z element error remains. (right) Variance of Each Element and Observation
Time:Firstly the variance of x is minimized, the second is y element, and finaly z element reduced.
4. 数値シミュレーション
4.1 シミュレーションパラメータ 数値シミュレーションに用いた各パラメータ値を Table 1 に
まとめる．母船の昇交点経度, 近地点引数および平均均点角は 0[rad] とした．雑音は白色ガウス乱
数を仮定した．I は単位行列を表す．
4.2 位置同定精度の数値シミュレーション結果 図 3 に拡張カルマンフィルタを用いた数値シミュ
レーション結果の 1 例を示す. 上は位置推定誤差, 下は分散行列 P の各対角成分を表す. 100 回シ
ミュレーションした結果の平均値を示す.
横軸は観測開始からの経過時間 [h], 縦軸は真値との差 [m] を示す. 実線が真値からの誤差 (Δx2 +
Δy 2 + Δz 2 )1/2 , 鎖線が Δx, 破線が Δy, 点線が Δz の各成分の推定誤差を表す.
観測開始から不可視（ローバの視野から母船が補足不可能）状態に入るまでのおよそ 5 時間の
観測により, 1[m] の精度でローバ位置同定ができていることがわかる. また x, y 成分の同定精度は
高く, z 成分の位置同定精度が悪いことがわかる. また, 誤差分散の推移から, まず x 成分が同定さ
れ, 次に y 成分, 最後に z 成分が同定されることがわかる.
4.3 計算量軽減法 電波受信時刻における母船位置 X r に，実測値 X r,ob を用いて，位置同定の
処理計算量を軽減する．すなわち，式 (16) で τ̂d を推定し，
ρ̂d = cτ̂d
(26)
ρ̂u = |X(tr,ob ) − Rz (τ̂d )xe |
(27)
ρ̂ = ρ̂d + ρ̂u
(28)
により得られる
と観測値 ρob = cτob との差分から，推定値 x(te ) を修正する．
本来であれば，下りの伝送時間 τd および上りの伝送時間 τu の 2 つを推定し， 受信時刻におけ
る母船位置を推定する必要があるが，受信時刻における母船位置 X r に実測値を用いると，下り
の伝送時間 τd のみの推定でよいことになる． これにより，観測更新の計算処理負担はおよそ半分
に軽減される．
図 4 に同一条件において通常の計算量で自己位置推定を行った結果と計算量を軽減した方法で
推定を行った結果を示す．図より，提案手法により電波の往復伝搬時間により自己位置が推定で
きることが分かる．また，計算量を軽減したにも関わらず， 同程度の精度で位置同定が可能であ
ることもわかる．
Estimation Error[m]
20
Regular
Calculation
㩷
Mitigated
㩷 㩷 㩷 㩷 㩷 㩷 㩷 㩷Version
㩷㩷㩷㩷㩷㩷㩷㩷㩷㩷㩷
表 1: 数値シミュレーションに用いたパラメー
タ値
15
初期推定誤差
1.4 × 10−4 [rad/s]
(period : 12[h])
300[m]×300[m]×600[m]
2500[kg/m3 ]
1[rad]
3000[m]
0.2
平均 0, 分散 10−8
天頂方向から 80[deg]
T
10 10 10
初期共分散
100I
小天体の自転角速度
10
5
0
0
1
2
3
4
5
6
Observation Time[h]
図 4: Localization Efficiency without Estimation
of X r :Instead of estimate Xr , the mitigated version used the observation data Xr directly.
小天体の各軸直径
小天体の密度 [3]
母船の軌道傾斜角
母船の軌道長半径
母船の軌道離心率
観測雑音
ローバの視界
5. 結言
本稿ではロボットの移動制御に不可欠な位置同定手法について，特に小天体表面における探査
ローバに対しても有効な位置同定手法を提案した．提案手法はローバと母船との電波の伝播遅延
時間を繰り返し測定することで実現する．提案手法は以下の 3 点で現実的である．
1. 往復の伝播遅延を母船側で測定することで，ローバへの搭載機器の負担を軽減可能．
2. 繰り返し測定することで，母船 1 機で位置同定可能．
3. 天体の大きさを問わず，直径数百メートル級の小天体表面上における位置同定も可能．
このほか，今年度は本稿で示した位置同定手法の感度などについて解析や非線形フィルタの適
用などに関する研究や自律ロボットの位置同定手法により環境推定に対して適用する手法につい
ても研究をすすめたが，本報告書提出〆切時には学会未発表であるため本稿では割愛する [4, 5]．
[謝辞] ローバの自己位置同定のために母船とローバーの間の往復伝搬遅延時間測定を用いる
考えは筆頭著者の中谷一郎先生，および吉光徹雄先生の指導による修士論文研究 (東京大学) に基
づいている．本稿で示した位置同定アルゴリズムや数値計算については著者が博士課程進学後に
行ったものであるが，論文発表の際に両先生には貴重なご意見をいただいた．ここに記して感謝
する．
参考文献
[1] Guochang Xu, ”GPS : Theory, Algorithms and Applications”, Springer-Verlag, Berlin, Federal
Republic of Germany(DEU), 2003
[2] 中野道雄, ”パソコンで解くカルマンフィルタ”, 丸善, 平成 5 年
[3] 道上達広, ”小惑星 25143ITOKAWA のレゴリス厚さと岩塊（ブロック）の数密度の推定”, 日
本惑星科学会誌, Vol.13, No.2, 2004
[4] 中西，金田，椹木，”小型無人航空機のための地磁気偏角推定法”，2007 年 12 月 20 日∼22 日
計測自動制御学会システムインテグレーション部門講演会 2007 にて発表予定
[5] Sayaka Kanata, Hiroaki Nakanishi, Tetsuo Sawaragi, Tetsuo Yoshimitsu, Ichiro Nakatani, Radio
Wave Based Localization of a Rover for a Small Planetary Body, The 27th IASTED International
Conference on Modelling, Identification, and Control -MIC2008-, (accepted)
共有自律系における機能分担の設計に関する研究
工学研究科機械理工学専攻
堀口 由貴男
Abstract: This work investigated a principle of human-machine interaction designs to facilitate human
operator’s correct awareness of the operating state (i.e., the mode) of the automation in their joint activity
to control their common target plant. Based on the comprehension that a human being is a proactive agent
who would try to affect the external world to obtain new cues to foster its understandings, the relationship
between input operations by the operator and output behaviors by the automation was focused on. Those
individual mapping structures were termed ‘input-output relations’ here. It was estimated that different
control modes with similar input-output relations should be easily mixed up with one another since interactions mediated by those relations would intimate a same mode. This hypothesis was verified through
an experiment using a driving simulator that examined human drivers on their recognitions of the operating mode of ACC (Adaptive Cruise Control) while driving with it. Another experiment confirmed that a
human-machine interaction design to differentiate modes in their input-output relations could successfully
reduce mode confusion errors.
Key words: Human-machine systems, function allocation, user interface design
1.
はじめに
機械化は，ある目的を達成するために作業システムが実現すべき諸機能のうちの特定の機能を
本来の作業主体である人間から新たにシステムに導入されたコンピュータ制御（自動化）に委譲
することを意味する．その目的はシステムのパフォーマンス向上や人間作業者の負担軽減と安全
にある．一方で，遭遇し得る状況を限定できない現実環境では自動化単独で対処できない事態が
必ず存在するために，作業システムは必然的に人間と自動化による共同作業の形態を取らざるを
得ない．複数の意思決定主体による円滑な共同作業を実現するには，各主体が共同作業者の状態
を適切に把握して協調作業の土台となる状況認識を整合させることが重要である．人間と自動化
の共同作業系において，そのような「コミュニケーション」が成立しない状況は，“single-minded
automation” [1] や “automation surprises” [2] として指摘される．人間同士の共同作業においても共
同作業者へのアウェアネス（気づき）は，タスク環境へのアウェアネスとともに，円滑な共同作
業遂行の必須要件である．そして，G. Klein らが提起 [3] しているように，それは人間機械系設計
の研究領域における大きな研究課題を構成している．本研究では，人間操作者と自動化が作業シ
ステムの制御機能を分担する共有自律系に対して，人間操作者と自動化が整合した状況認識を共
有するためのインタラクション設計について，特に人間操作者の自動化の稼働状態に対するアウェ
アネスについて検討を行った．
作業対象が備える性質がすべて直接知覚可能な形式で表出しているとは限らないため，知覚者
はまず何らかの操作を加えることで積極的に対象に働きかけてそこに隠されている状態や構造を
揺さぶり出し，さらに作業場に人工的な秩序や構造を与えて代替の効果的な認知資源を作り出す
ことで，対象の把握にインタラクションをうまく活用している．すなわち，人間は積極的に外界
に対して働きかけることを通じて活動に必要な意味解釈をつくり上げる能動的な認知主体である．
このような理解に基づいて，行為と密接に結びつくフィードバック情報を人間の状況認識のため
の効果的な手がかりであると見なし，操作者の操作入力とそれに対する自動化の反応挙動の間の
対応構造に注目した．この対応構造を『入出力関係』と呼ぶことにする．
本研究ではこの入出力関係に注目し，人間−機械間のインタラクション設計が人間操作者の自動
化に対するアウェアネスに与える影響について被験者実験を用いて調査した．実験環境では，ACC
(Adaptive Cruise Control) 機能を自動化として実装したドライビングミュレータを利用し，その制
御モードを調査するアウェアネスの対象である自動化の稼働状態と見なした．そして，運転中の
ドライバのモード認識を問う実験を実施し，入出力関係とモード認識の相関について検討した [4]．
cruise-on
engaged
[when speed <
lower limit]
active
carfollowing
target car is out
detect car in front
constantspeed
brake
or
cancel
press gas pedal
release gas pedal
override
set
and
[when speed
>=
lower limit]
resume
canceled
off
not-engaged
armed
on
idle
図 2: ACC の制御モードと制御モード間の遷移
運転操作時の自動化動作モードの識別
2.
ACC による定速走行・車間距離制御機能を模擬
したドライビングシミュレータ（図 1）を構築し，
それを用いてドライバの運転中のモード誤認識を
調査する被験者実験を実施した．
2.1 ACC 機能の複雑さ ACC は，設定速度を正
確に維持して自車両を走行させるとともに，安全
な車間距離を保つように走行速度を適宜調節して
先行車に追従する機能を備える自動車運転支援の
ための自動化技術である．その振る舞いは多数の
制御モードと制御モード間の複雑な遷移によって
実現される．図 2 は，本研究のドライビングシミュ
図 1: ドライビングシミュレータ
レータ上で実装した ACC のマシンモデルを示し
ている．idle（未起動状態），armed（待機状態），
canceled（一時中断状態），override（オーバーライド状態），constant-speed（定速走行状態），carfollowing（追従走行状態）の 6 つの制御モードがあり，それらは共通する遷移イベントによって
いくつかのグループに階層的にまとめられる．図 2 において実線で描かれた有向リングはドライ
バが起動するモード遷移を表している．それらは，ステアリングコントローラ内側に取り付けら
れた 4 つのボタン（図 1 右下の拡大図）とアクセル／ブレーキペダルを使って操作される．一方，
破線で描かれたリンクは ACC の独自の状況判断によって自動的に実行されるモード遷移を表す．
両リンクを修飾するラベルは，モード遷移を引き起こすイベント条件を表している．
2.2 モード識別実験 ドライビングシミュレータ上で仮想市街地を走行する作業課題を被験者に
課し，運転操作中の自動化の動作モード認識を調査する実験を実施した．多数の車両が走行して
いる市街地の道路は信号機によって交通がコントロールされており，被験者はそのような環境内
を実験者に指示されたコースに沿って運転した．実験では，約 10 分間の市街地走行の間に約 1 分
間隔で不意の強制的な割り込みが入りシミュレーションが停止する．この割り込みが発生すると，
モニタ画面は全体が白く塗り潰され，3 秒間 ACC の動作状態に関する情報をディスプレイから入
手できない状態になる1 ．これは被験者が認識している ACC の動作モードを調査するもらうこと
1
運転中の ACC の動作モードはインストルメントパネル内のインジケータの表示状態から知ることができる．
armed
8.3
canceled
6.9
armed
0
12.2
canceled
0
0
constantspeed
2.83
2.83
2.83
carfollowing
2.83
2.83
2.83
0
override
2
2
2
2
idle
armed
canceled
constantspeed
0
1.0
6.3
carfollowing
0
0
1.0
12.0
override
0
4.1
4.3
4.3
idle
armed
canceled
0
constantcarspeed
following
表 1: モード誤認識発生率の分布
2
constantcarspeed
following
表 2: モードベクトル間の距離
を目的としたもので，実験者は被験者に動作モードについて瞬時の回答を求めた．延べ 35 人の被
験者から得た 303 個のサンプルデータのうち，被験者が実際とは異なる ACC の動作モードを回答
したものが 59 個あり，全体的なモード誤認識率は 19.5% であった．表 1 にモード誤認識の発生率
の分布を実験結果として示す．
2.3 システムの入出力関係とモード誤認識の相関 本研究の実験系においてドライバが弁別可能
な機械システム（車両−自動化系）の入出力関係は以下の 4 種類に整理される．
1
2
3
4
アクセルペダルを踏み込んだ時にアクセルの開度に応じて加速するかどうか
アクセルペダルから足を放した時に減速するかどうか
“resume/accelerate” ボタンを押下した時に目標走行速度の数値が増加するかどうか
“set/coast” ボタンを押下した時に目標走行速度の数値が減少するかどうか
被験者のモード判断の瞬間における入出力関
*
係の状態を属性として，さらに 1 つ前の動作
response to footing gas pedal = nothing
lead-car indication = nothing
モードと直前の遷移イベント2 ，および各イン
< constant-speed >
ジケータの表示状態を属性に加えて，被験者の
lead-car indication = displayed
< car-following >
モード認識がそれらの属性に関するどのような
response to footing gas pedal = accelerating car
car-icon indication = nothing
情報に依存していたのかを決定木抽出アルゴリ
< idle >
ズムである C4.5 [5] を用いて推定した．図 3 に
car-icon indication = red
< override >
抽出された決定木の一部を示す．被験者のモー
car-icon indication = yellow
< canceled >
ド認識が正確であればモード遷移履歴やインジ
car-icon indication = blue
ケータの表示状態だけで被験者の判断がモデル
set-speed indication = nothing
< armed >
化できるはずだが，実際のモード認識は不完全
set-speed indication = displayed
< canceled >
なために入出力関係に基づく決定が導出された
モデルには含まれている．すなわち，被験者の
モード認識は彼らの記憶やユーザインタフェー 図 3: C4.5 [5] によって得た被験者のモード判断
ス上の表示データだけでは説明できないことが を説明する決定木
示唆されている．
被験者のモード認識の入出力関係への依存についてさらに検討するために，モード間の入出力
1 の関係では「加速する」場合には T (1 0)，
「加速しな
関係の類似性を定量的に評価した．例えば
T
1 ∼
4 の入出力関係を独立した 2 次元ベクトルに変換して，
い」場合には (0 1) というように，
全入出力関係についてこれを連結し，2 値ベクトル（以下，モードベクトル）として制御モードを
コード化した．このコード化により入出力関係の観点でのモード間の非類似性はモードベクトル間
の距離で評価する．モードベクトル間のユークリッド距離を表 2 にまとめる．この表より，{idle,
armed, canceled} や {constant-speed, car-following} というモードの組合せはモードベクトル間の距
離が 0 であり，入出力関係が完全に一致することが分かる．そして，表 1 の結果と対照すること
2
これらはモード遷移履歴に関するドライバの記憶に対応する．
で，これらの組合せでは高い確率でモード誤認識（モードの混同）が発生していたことが確認で
きる．さらに，モードの組合せ 1 つを 1 事例として，モード間距離とモード誤認識発生確率の間
の相関計数を計算したところ，−0.83 という非常に高い負の相関を得た．これらの結果より，被
験者のモード認識が他のモードとの入出力関係の類似性に強く影響を受けていたと考察できる．
3.
入出力関係の差別化がモード認識に与える影響
本研究ではさらに，前節で明かになった入出力関係とモード認識の相関をうまく利用して人間操
作者の自動化に対するアウェアネスを向上させるインタラクション設計について検討した．具体的
には，識別が必要なモード間で入出力関係を差別化するようなインタラクション挙動の変化を自
動化に実装し，そのモード認識改善への効果について調査した．表 1 の実験結果において誤認識率
が非常に高かった {constant-speed, car-following} の混同を改善の対象とし，両モードベクトル間の
距離を広げるための「ステアリング操作時に感じる反力が硬いか柔らかいか」という新たな入出力
5 を追加した．constant-speed モードの時はステアリング操作に対する反力を強くし，それ以
関係
外のモードでは弱くすることでこの入出力関係の差別化を実現した．これにより，constant-speed
と car-following のモード間距離は 0 から 1 に増加する．そして，先の実験と同じ条件で実験を実
施し，延べ 13 人の被験者に対して 213 個のサンプルデータを得た．被験者が実際とは異なる ACC
の動作モードを回答したのは 29 個で，全体のモード誤認識率は 13.5%であった．そして，2 つの
モードの混同は 12.0% だったものが 2.5% と大きく低減できた．この結果より，システムの入出力
関係とドライバのモード認識の結び付きの強さを再確認するとともに，入出力関係に関してモード
間を差別化するインタラクション設計がドライバのモード認識の改善に寄与することを確認した．
4.
おわりに
本研究では，入出力関係の類似性が複雑な自動化の稼働状態に対するアウェアネスに影響するこ
とを，ACC 機能を実装したドライビングシミュレータを用いた被験者実験を通じて確認した．さ
らに，モード誤認識の危険を低減するための，モード間で入出力関係を差別化する人間−自動化
間のインタラクション設計の有効性を同様の被験者実験によって確認した．
人間と機械の共有自律系において，システムのモードの管理と把握はユーザの主たる作業では
ない．それは，複雑なモード遷移を内包する機械を相手にしたことで発生した副次的なタスクで
ある．それゆえユーザの本来の作業とモードの管理が彼の認知資源（特に，知覚モダリティ）に関
して競合する場合，たとえシステムの動作モードについての情報が表示装置に正しく提示されて
いたとしても，その種の情報は十分には活用されない可能性が高い．このようなタスク環境に対
して行為と密接に結びつく知覚情報を人間−機械間のイタラクションに埋め込むことは，人間−
機械間の状況認識整合のための効果的な支援になると考える．
参考文献
[1] Sheridan, T. B. Humans and Automation: System Design and Research Issues. Wiley Series in
Systems Engineering and Management. Wiley-Interscience, July 2002.
[2] Sarter, N. B., D. D. Woods, and C. E. Billings. Automation Surprises. Handbook of Human Factors
& Ergonomics (ed.Salvendy, G.). Wiley, 2nd edition, 1997.
[3] Klein, G., D. D. Woods, J. M. Bradshaw, R. R. Hoffman, and P. J. Feltovich. Ten Challenges for
Making Automation a “Team Player” in Joint Human-Agent Activity. IEEE Intelligent Systems,
Vol. 19, No. 6, pp. 91–95, 2004.
[4] Horiguchi, Y., R. Fukuju, and T. Sawaragi. Differentiation of Input-Output Relations to Facilitate
User’s Correct Awareness of Operating Mode of Automated Control System. Proceedings of 2007
IEEE International Conference on Systems, Man and Cybernetics, pp. 2570–2575, 2007.
[5] Quinlan, J. R. C4.5: Programs for Machine Learning. Morgan Kaufmann Publishers, 1993.
環境適応性に優れたヒトの手の機能の解明および
その機械的実現
工学研究科機械理工学専攻
横小路 泰義
Abstract: For the purpose of realizing the functions of the human hand on mechanical systems, fundamental understanding of the human hand is necessary through intensive observations of manipulation tasks by
humans. As an example of a skilled manipulation task, we target origami folding. We already succeeded in
developing a robot hand that can fold an origami work “Tadpole”. However, the robot often fails folding and
the success rate is very low (about 20%). This year, we aim to increase the success rate of “Valley-fold” by
using the statistical motion data directly taught by a human. To model the statistical motion by the human,
we use the Hidden Marcov Model (HMM). We then extract some critical motion segments where sensory
feedback seem to play an important role. We use the HMM to control the robot for such critical segments,
taking the current sensor signals as the inputs. We show that the success rate is improved by the proposed
method.
Key words: Multi-fingered hand, Origami folding, Dextrous manipulation, Flexible objects
1. はじめに
本研究課題では，環境に対して適応的に働くヒトの手
の機能の解明およびその機械的な実現を，機構，センシ
ᠲ૞ᜰ㧝
࿕ቯᜰ㧝
ング，アクチュエーションの側面から目指し，人の動作
を詳細に観察し解析を重ね，その結果に基づく設計論や
数理モデルを提案し，機械的な機能の実現を通してヒト
の手の機能の本質的理解を目指してきた．一昨年度から
「物体を操る手」にテーマを集中し，手を用いる複雑な動
作として「折り紙」という具体的作業を取り上げて，人
࿕ቯᜰ㧞
ᠲ૞ᜰ㧞
間による作業の観察から機能の機械的実現までを試みて
図 1: 折り紙ロボットの概観
きた．
一昨年度は図 1 に示すような各指 3 自由度からなる 4 本の指をもつハンドを設計・試作し，昨
年度は，図 2 に示す「おたまじゃくし」[1] という作品を一連で折ることに成功した [2]．しかし，
現状でこのロボットハンドは様々な問題点を抱えており，特に問題となるのがタスクの成功確率
の低さ (約 20%) である．現状ではハンドの指先の軌道は人間が事前に試行錯誤的に与えたもので
あり，作業中に視覚や力覚を用いたフィードバックによる軌道修正も行っていない．そのため紙
を折る際の失敗の原因のほとんどは，与えた軌道の不適切さや不確実な紙の挙動によるものと考
えられ，より適当な指先の軌道を与えると共に人間と同様にセンサフィードバックに基づいた器
用な操りを実現する必要がある．
柔軟物操作のような器用なタスクをロボット
に行わせる一つの方法として，人間の動作を教
示させることが考えられる．巧みな人間の運動
をロボットに移植させようとする研究はこれま
でにもいくつか行われてきた [3][4]．他の方法
図 2: おたまじゃくしの折り手順
として，強化学習のようにロボット自身に最適
な動作を学ばせる方法も考えられるが，折り紙のような複雑なタスクにおいては効率が悪く，う
まくいく動作が得られる保証はない．本研究で開発したロボットハンドは，機構自体は人間の手
の外観とは異なるが，人間が折り紙を折る時と同じように折るように設計したため，人間にとっ
ては比較的直感的に動作教示を行うことができる．よって本研究では，教示による動作改善の手
法をとることにする．
紙のように対象物が確率的な振る舞いをする場合，人間の動作はその挙動に合わせて試行毎に
変化する．従って教示も一回だけでは不十分であり，複数の試行を記録する必要がある．また，記
録した教示データの運動軌跡を単に再生するだけでは，毎回異なる紙の挙動に対処しきれないの
で，複数回行った教示データの統計的性質から人間が行った適切なセンサフィードバック則を抽
出する必要がある．
このような対象物の動きに対応して変化す
る人間の動作をロボットに模倣させる研究とし
て，大西ら [5] は患者の抱き上げタスクを対象
に，DP マッチングを用いて作成した人間の動
(a) 成功
(b) 失敗
作モデルから，運動の生成を行っている．人間
の運動をロボットに反映させることでタスクの
図 3: 谷折りの失敗例
成功確率が改善されたが，大西らの方法では空
間的な運動の修正方法を明示的に設計できなければならない．折り紙作業のような複雑なタスク
では人間の運動の修正方法は明示的ではないため，人間が与えた教示データから修正方法を抽出
する必要がある．
本研究では，人間の教示動作から，運動の軌跡だけでなく，人間が行うセンサフィードバックも
抽出し，ロボットに再現させる．その際に，実際の人間の作業ではセンサ情報を使って運動をし
ている部分と，そうではなくランダムに運動している部分とが考えられるため，まず初めに，セ
ンサ情報を使って運動している部分の絞込みを行う．そして，その部分のみに教示データから抽
出したセンサーフィードバックを適応する．このような方法を用い，タスクの成功確率を改善す
ることを目指す．
x [mm]
2. 折り紙ロボット
150
200
300
250
350
0
y [mm]
50
100
11
z [mm]
-100
2.1 対象作業 図 1 に示す折り紙ロボットは，主に
10
-80
折り紙を操る目的に用いる 2 本の操作指 (操作指 1 &
-60
2 ) と主に折り紙を固定する目的に用いる固定指 (固
12
-40
4
5
「おた
定指 1 & 2 ) の 4 本の指からなる．現状では，
1
7
6
-20
8
まじゃくし」を一連に折る作業の成功確率は 20% 程
9
2
13
3
0
度でしかない．主に失敗するところとしては図 3 に
20
示すように谷折り時に指が紙から外れるという失敗
40
がある．現在の指先の軌道は図 4 に示すように，ま
図 4: 指先軌道の例
ず離散的な指先の参照点をあらかじめ与え，それを
直線補間するように軌道を生成している．参照点の位置は試行錯誤的に決定され，適当な参照点
の数も明確ではないため，このような方法で与えられた軌道では不十分であり，また決められた
軌道を再生するのみでは試行毎の紙の挙動の変化には対応できない．
本研究では，一連の「おたまじゃくし」の折
りを対象とする前に，ます基本的な折り紙作業
である谷折り作業に特化して成功率の向上を目
指すことにする．谷折り作業は，おたまじゃく
しの中には 1/2 折り，1/4 折り，尻尾折りの計 3
Manipulation finger 1
Holding finger
回含まれるが，そのなかでも図 5 に示すような
Holding finger beneath the paper
Manipulation finger 2
最も基本となる 1/2 折りを取り上る．
図 5: 谷折り作業の手順
2.2 指先センサ 折り紙作業において重要な
センサ情報として考えられるものに視覚と指先の力触覚が考えられる．人間による作業において
視覚は最も重要な情報源となるが，ロボットの場合，折り紙の３次元的な変形を十分な精度でか
つリアルタイムで取得できるセンサをすることは容易ではない．力触覚においても，人間の持つ
機能と同等のセンサを実装することは容易ではないので，ここではまず紙を折る際の反力を検出
するセンサを実装することにする．
折り紙を折る際の指先にかかる力はかなり微小であり，その力はおおよそ 0.1∼10.0gf である．
市販されているセンサでは，このような微小な力を検出するのに適当なものが見当たらなかった
ため，ひずみゲージを 3 箇所に配置して , , ３軸方向の力の検出を可能とする力センサを自作
した．谷折り時に実際に紙を捕らえて折る役目をする指は操作指 2 であるため，このセンサは操
作指 2 のみに取り付た．
指先に対しての紙の傾きが変化すれば，３軸方向の力の成分が変化するため，このセンサでも
ある程度折り紙の状態を検出することは可能であると思われる．実際に，谷折り作業を行った際
の指先センサの履歴を見たところ，谷折りが成功した場合と失敗した場合とで 軸方向の力履歴
に顕著な違いが見られ，今回作成した力センサを用いてある程度紙の状態を認識できることが分
かった．
3. 人間による動作教示と HMM によるモデル化
Human hand
3.1 人間による直接教示 人間の教示データ
Sensor 2
を取得する方法として，本研究では図 6(a) に示
すように，人間がハンドの指先を直接持って，
Sensor 1
指先を動かし，谷折り作業を行う．教示データ
としては谷折り作業中の指先位置 ，速度 ，
前節の微小力センサ (Sensor1) による力 を記
録する．直接人間がハンドを持って教示を行う
(a) 直接教示の様子
(b) ２つの力覚センサ
ためには，ロボットハンドがバックドライバブ
図 6: 人間による直接教示
ルである必要があるが，このハンドは鉛直方向
にボールねじを用いているためこの方向はバックドライバブルではない．そこで紙の反力検出用
のセンサ (Sensor 1) とは別に教示用に市販の力センサ (Sensor 2) を図 6(b) のように取り付けること
で，人間が鉛直方向に動かそうとした際に加えた力に応じてロボットが動くようにした．
教示の際は図 6(a) のように 2 本の固定指は紙を押さえつけたまま固定されているとし，操作指
2 本をそれぞれ両手で持って教示データをとった．試行回数は 20 回とした．図 6(a) に実際に得ら
れた指先軌道の教示データを示す．図から分かるように，人間が教示した折りの軌道には試行毎
に揺らぎがあることが分かる．
z [mm]
-80
Human demonstration data
3.2 HMM による人間の教示動
Averaged trajectory
作のモデル化 教示データのモ
-60
デル化には隠れマルコフモデル
-40
(HMM) を用いる．HMM は音声
認識の分野で盛んに用いられて
-20
いるが，稲邑らによる見まね学
習の研究 [6] や板橋らによる人間
0
の技能の獲得を目指した研究 [7]
150
200
250
300
x [mm]
などのように，近年，ロボット
(a) 教示データと平均軌道 (b) HMM の各状態の平均位置
の運動や人の技能のモデル化に
も応用されている．このような
図 7: 教示データと HMM によるモデル化
時間的，空間的にばらつきのあ
る人間の動作のモデル化を行うのに HMM は適している．
HMM の学習のための出力時系列データ は取得した教示データから とし，出力の確率分布 には次に示すような単一ガウス分布を仮定
した．
(1)
HMM の学習にはよく用いられる手法である Baum-Welch Algorithm[8] を用いた．その際，学習
すべきパラメータ は初期状態確率 ，状態遷移確率 ，状態 のガウス分布の平均 とその分
散 である．図 7(b) に学習した結果の各状態における指先の平均位置を示す．このように HMM
を用いることで作業が適切な位置で分節化されていることが分かる．なお，状態数が適当かどう
かの判断には AIC(Akaike’s Information Criterion)[9] を用いた．また図 7(a) には，モデル化された
HMM の出力を平均して得られる平均軌道も示してある． HMM によって生成される平均軌道が，
教示データの平均として適当であることが分かる．
4. 教示データに基づいた動作生成
Variance X Correlation factor [(mm/s) 2]
Velocity variance [(mm/s)2 ]
4.1 センサフィードバックが必
B
A
要な区間の抽出 本研究ではタ
スクのセンサフィードバックを
行っていると思われる重要区間
を抽出するために，まず速度の
分散に注目する．前節で教示デー
State
State
タにより学習した HMM の状態
における速度の分散をプロット
(a) 速度の分散
(b) 速度の分散 力との相関
したものを図 8(a) に示す．ここ
図 8: 速度分散と力との相関による動作の分節化
で速度の分散は，次式のように
共分散行列 を対角化した対角項の和として与えている．
2000
1000
0
0
10
20
600
400
200
0
0
10
20
(2)
図 8(a) を見ると教示データには速度の分散の大きいところとそうでないところがあることが確認
できる．速度の分散が小さいところはほぼ確定された動作を行っているところと考えられ，あまり
動作に揺らぎのない部分である．一方，速度の分散が大きいところは，人間が力情報に頼らずに
指先を動かしている部分と，指先力に応じて速度を変えて制御を行っている部分とがあると考え
られる．今回の谷折り作業では途中に指先が紙に触れない部分が存在するため，その部分 (図 8(a)
の領域 A) は，教示動作において明らかに力センサによるフィードバック動作をしていない区間で
あるのであらかじめ省いておく．
さらに，速度分散の大きなところから力情報によるセンサフィードバックを行っている部分の
みを抽出するために，力と速度の相関を考える．力と速度の相関としては正準相関係数 [10] を用
いる．正準相関はベクトル ， の相関を求める際に， ， それぞれをあるベクトル ， に
射影したときの値 ， を考え，その値同士の相関が最も大きくなるように相関係数およ
びベクトル ， を求めるものである．
具体的には状態 のガウス分布の分散 を用いて以下の式の を求める．
(3)
(4)
式 3 より は の固有値であり， は対応する固有ベクトルである．固有値
の中で最も絶対値の大きいものに対応する相関係数 を第 1 正準相関係数，2 番目に大きいも
の を第 2 正準相関係数という．この第 1 正準相関係数を力との相関の重みとして速度の分散に
掛け合わせる．これにより速度分散の大きい部分からさらに力情報を用いて運動を行っている部
分の抽出を行う．
速度の分散に相関係数を掛け合わせたものを図 8(b) に示す．図で値の大きいところが人間が力
情報を用いて制御を行っていると予想されるところであり，折り作業において重要な部分であると
考えられる．図から速度分散と相関係数を掛け合わせた値が大きくなっている部分が 3 つあるこ
とが分かる．1 つ目の山 (状態 3 から 5 付近) は，紙を押さえてたわませている部分に相当し，2 つ
目の状態 10 付近の山は操作指 1 を紙の下にもぐりこませるために，操作指 2 で紙を引っ掛けなが
ら持ち上げる部分に相当する．さらに最後の部分は実際に操作指 2 が紙を折り始める部分に相当
する．今回は折りの失敗に直接関係していると考えられる 3 つ目のピークに相当する図 8(b) の領
域 B に注目し，この部分に人間による教示データから抽出されるフィードバック動作を適応する．
4.2 統計モデルに基づくセンサフィードバック動作の生成 3.2 節で作成したモデルを用いて，指
先の運動の生成を行う．まず，方法 1 として人間の平均の指先軌道を与える制御を行う．これは
稲邑らの手法 [6] を用いて，学習した HMM から平均の指先軌道を生成し，その軌道に指先を追従
させる方法である．図 7(a) に示した平均軌道は，この方法によって作成したものである．
次に方法 2 として，4.1 節で抽出した領域にセンサフィードバックを適応する．まず HMM をセ
ンサフィードバック区間として抽出した部分とそうでない部分とで分割し，後者の部分に関して
は，方法 1 同様，平均軌道を作成し，その軌道に追従させる．次に前者の部分に関しては，その
ときの状態 の確率分布から，現在の指先にかかる力 の条件付き分布の平均速度 を以下のよ
うに与える．
(5)
なお状態の遷移については，モデルから求められる平均的な状態遷移に従って進むものとする．
5. 実験
B
0
5.1 実験結果 前節で提案した運動生成法を用いて，実際の
15
A
5
谷折り作業が改善されるのか実験を行った．その際，谷折り
動作の失敗の判断基準として，最終的に折りが終了した状態
で紙から指が外れてしまった場合に加えて，図 9 に示すよう
図 9: 谷折りの失敗の定義
に端から 5 mm の領域 B に指先がある場合も，紙が外れそう
な不安定な状態と見なし失敗とた．逆に図 9 の領域 A に指先
がある場合を成功とした．従来の方法と今回行った 2 つの方法の比較を表 10 に示す．なおここで
従来の方法とは図 4 に示した方法である．
表より教示データの平均軌道を与えた方法 1
図 10: 実験結果
で，従来の方法から成功率がかなり改善されて
success fail total
いることが分かる．これは人間の動作から軌道
Original method
17
13
30
を生成することで，指先の軌道が従来のものよ
りも適当なものとなったためである．また作業
Averaged trajectory
27
3
30
時間においても，従来の方法では 1 回の谷折り
Sensory feedback
30
0
30
作業で 45 秒ほどかかっていたが， 15 秒程度で
作業を行うことができた．しかし，この方法でもまだ失敗することがあり，平均軌道だけでは紙
の挙動に対応しきれていないことを意味する．次に前節でのセンサフィードバック動作を入れた
方法 2 を行った場合では，30 回の試行で失敗は見られず，紙のゆらぎに対応したより安定な折り
を実現することが確認できた．図 11 に方法 2 を用いたときの試行毎の軌道の違いを示す．紙の反
力から軌道を変化させている様子がわかる．
Averaged trajectory
5.2 考察 前節で提案手法により作業の改善を行うことに成
Trajectories with
-60
sensory feedback
功したが，今回設定した成功の定義では，紙が折れていれば
成功と見なしており，次の作業である紙の 2 頂点合わせへの
-40
移行のしやすさは考慮されていない．人間が実際に紙を折る
際は，折りの見栄えや次の作業の行いやすさを考慮して折る
-20
ため，今後はこのような点も評価して折りの失敗を定義しな
おす必要がある．
0
もう一つの課題として，教示回数の多さがある．今回は谷
150
200
250
300
折り作業のみで 20 回の教示を行ったが，一連の作業において
も 20 回もの教示を行うのは，教示者にとって負担である．も
図 11: センサフィードバックを用
しも今回の谷折り作業の教示結果から，おたまじゃくしに含
いた際の指先軌道
まれる他の 2 回の谷折り作業にも適応できるような一般性の
ある制御則を抽出することができれば，他の 2 回の谷折り作業に対しての教示回数を減らすこと
ができる可能性がある．そのためにも統計的性質から抽出したフィードバック則の汎化が課題と
なる．
6. まとめ
本研究では，折り紙作業の中の基本となる谷折り作業において人間の教示データを HMM を用
いてモデル化し，そのモデルを用いてより適当な指先軌道の生成を行った．また教示データから
抽出した作業の重要部分に対し力情報を用いたフィードバック動作を適用する方法を示した．こ
れらの方法を用いることで実際に作業の成功確率を改善できることを実験により示した．
ロボットには目的の作業を行える機能を有することが必要であるだけでなく，目的とする作業
の動作を生成する機能も必要である．動作生成の方法としては人間による直接教示が最も簡単な
方法の一つであるが，本研究で得られた成果は挙動に揺らぎのある柔軟物体を扱うような複雑な
作業に対する直接教示のための方法の一つを示したといえる．フィードバックが必要な区間の抽
出方法として現状の方法が妥当であるかは今後さらに検討する必要があるが，このような作業の
中の重要区間を知ることは，人間の技能の理解にも繋がると考えられる．
参考文献
[1] 新宮文明, “おりがみくらぶ”, http://www.origami-club.com/
[2] K. Tanaka, Y. Kamotani, and Y. Yokokohji, “Origami Folding by a Robotic Hand,” Proc. International Conference on Intelligent Robots and Systems (IROS 2007), October 29-November 2, San Diego, CA, USA,
pp.2540–2547, 2007.
[3] Y. Kuniyoshi, M. Inaba, and H. Inoue: “Leaning by watching: Extracting Reusable Task Knowledge from
Visual Observation of Human Performance”, IEEE Transactions on Robotics and Automation, Vol.10, No.6,
pp.799-822, 1994.
[4] 平井真一, 野口英昭, “変形しやすい管状物体の挿入作業における人間の運動と作業識別則のロボットマ
ニピュレータへの移植”, 日本ロボット学会誌, Vo.15, No.8, pp.1172-1179, 1997.
[5] 大西正輝, 小田島正, 羅志偉, “環境と接するロボットの感覚運動統合による動作模倣”，電気学会論文
誌，Vol.125-C, No.6, pp.856-862, 2005.
[6] 稲邑哲也, 谷江博昭, 中村仁彦, “連続分布型隠れマルコフモデルを用いた時系列データのキーフレーム
抽出とその復元”, ロボティクス・メカトロニクス講演会, 2P1-3F-C6, 2003.
[7] 板橋界児, イエ世薫, 鈴木達也, 大熊繁, “隠れマルコフモデルの基づく人間の技能獲得”, 計測自動技術学
会論文集, Vol.34, No.8, pp.890-897, 1998.
[8] 中川聖一, “確率モデルによる音声認識”, 電子情報通信学会, 1988.
[9] 赤池弘次, “情報基準量 AIC とは何か—その意味と将来への展望”, 数理化学, pp.5-11, 1976.
[10] 浅野長一, 塩谷実, “多変量解析論”, 共立出版, 1967.
歩行面の振動と歩行リズムの相互作用に関する研究
（引き込み現象を実現する歩行モデルの作成に関する研究）
工学研究科機械理工学専攻
宇津野秀夫，後藤一孝
Abstract: Human locomotion system has ‘pull-in effect’, in which if the ground swings with specific
frequencies in one’s frontal plane, a human’s gait rhythm changes into the ground swinging
frequency from one’s natural gait rhythm. However the cause of the effect is still not clear. This paper
notices human’s posture controlling system, and refers to the cause of pull-in effect with suggesting a
dynamical human gait model which has two legs and posture controller. Firstly the features of pull-in
effect are clarified experimentally and secondly the posture controller which works for stabilizing
frontal body behavior is studied. Then a dynamical human gait model is established from the results.
The model consists of two beams and a mass representing legs and center of gravity of one’s body. It
basically behaves passively, but only the posture is controlled actively. The model shows the features
of pull-in effect very well, and looking at the behavior of the model, it is explained that the posture
controller causes pull-in effect.
Key words: Synchronization, human walk, vibration, Human Gait Model, Pull-In Effect, Posture Control
1. はじめに
人間の歩行リズム生成機構には，人体の前額面（体の正面）に対して特定の振動数域で横振動する歩行面
上を歩行する時，歩行リズムが歩行面の振動数に同調するという引き込み現象が生じることが報告 1)されている．
そこで歩行を模擬した力学モデルを作成し，そのモデルを用いて強制変位加振台上で引き込み現象を実現し，
引き込み現象の発生を説明する研究が行われた．しかし引き込み現象の性質を物理的に解明することに重き
をおいたため，人間に比べて引き込まれる振動数域が広いことに加え，人間の横方向加振力を定量的に評価
できていない問題点が残されている．
本研究ではこの力学モデルをベースに，人間と同等の横方向加振力と引き込み振動数域をもつモデルに
拡張し，強制変位台における引き込み現象の定量的な解析を可能にした．
2. コンパス形歩行モデル
2.1 コンパス型モデル 力学モデルを作成する上で，McGeerの提案する受動歩行モデル 2)を用いる．受動
歩行モデルとは，人間の足を模した二脚構造を有するモデルが，各関節に能動的なトルクを加える制御なしに，
倒れこみ動作により自立的に歩行運動を生成する二足歩行モデルであり，その挙動は実際の人間の歩行をよ
く再現する．McGeer の受動歩行モデルでは矢軸面内の挙動と前額面内の挙動が独立にモデル化されている
ので，本稿では前額面内の歩行面の振動が歩行に与える影響について考察するため前額面内の歩行運動を
模擬する力学モデルのみを用いる．図 1(a)に示すような人間の前額面内歩行運動を模擬するため，図 1(b)に
示す二脚を有するコンパス型受動歩行モデルを用いる．図 1(b)に示すように，錘の質量，脚部長さ，脚の開き
角をそれぞれ m, l , ϕ 0 とし，脚部は質量を持たないものとする．ここでコンパスモデルの脚部長さとは人間に
おける足裏から重心（へそのあたり）までの長さに相当する．
重心位置
m = 60kg
m
（へそのあたり）
左足
l
θ
ϕ0 = 8o
(< 0)
l = 0.9m
右足
θ
(> 0)
θ =0
振動台
振動台
0
(a) 歩行時の人間の模式図
x0
(b) 歩行時の人間のモデル化
図 1 前額面内の歩行モデル
(a) 左足を軸とした
回転運動
(b) 両足着地
図 2 歩行モデルの挙動
(c) 右足を軸とした
回転運動
このモデルの運動は 2 種類の挙動に分けられる．一つは図 2(a)のように左足を支点とする倒立振子としての挙
動，他方は図 2(c)のように右足を支点とする倒立振子としての挙動である．強制加振を加える台の変位を x 0 と
表す．また，倒立振子の回転運動は ϕ で観測するのが一般的であるが，コンパスモデルの挙動は人間の重心
の横揺れと対応させていることに考慮し，中心軸の鉛直線からのずれ θ で計測する．ここで ϕ と θ の関係は次
のようになる．パラメータは実際の人間に基づき表 1 のように定めた．
表 1 定めたパラメータ
右足支点
左足支点
( > 0)
θ = −ϕ0 + ϕ ( < 0 )
θ = ϕ0 + ϕ
質量 m
(1)
脚部長さ l
60 kg
(2)
脚の開き角
0.8 m
ϕ0
8°
2.2 台上のコンパスモデル 減衰を考慮したコンパスモデルの運動方程式を考える．図 2(a)(c)のように片
足を軸とした回転運動は，運動する台の上の倒立振子の挙動とまったく同じである．よって倒立振子の挙動を
扱いやすい ϕ に関する運動方程式を用い，式(1),(2)を用いてコンパスモデルの挙動θ に変換することとする．
図 3 より
振子
J ϕ&& = Vl sin ϕ − Hl cos ϕ − cϕ&
(3)
d2
m 2 ( x0 + l sin ϕ ) = H
dt
d2
m 2 ( l cos ϕ ) = V − mg
dt
y
(4)
(5)
V
※ 重心まわりの慣性モーメント J = 0
H
台車
Mx&&0 = U − H
となり，式(3)～(6)から， H とV を消去して次式を得る．
ϕ&& −
ϕ
x
&&
x
g
sin ϕ + 0 cos ϕ + Cϕ& = 0
l
l
(C = c
2
)
( M + m ) &&x0 + mlϕ&& cos ϕ − mlϕ& sin ϕ = U
l
H
V
Mg
(6)
ml
mg
(7)
U
図 3 倒立振子
(8)
2.3 コンパスモデルの撃力と減衰 減衰のあるコンパスモデルの挙動のうち，右足を支点に回転する
場合を考える（図 4）
．その挙動θ は図 5 のようになる．時刻 t1 に右足を支点とした回転運動がはじまり，
時刻 t2 に最高点に達し，時刻 t3 に左足が接地し右足支点の運動が終了する．最高点の角変位をθ max とす
る．図 5 から読み取れるように，減衰の影響で，回転速度の絶対値 θ& が時刻 t1 に比べて，時刻 t3 では遅
くなっていることがわかる．そのままでは，左足支持による運動後，点線波形（
）のように再び最高点
θ
θ max
時刻 t1
t2
t3
0 t
1
t2
t3
時刻 [s]
−θ max
図 4 右足支持の様子
図 5 減衰をもつコンパスモデルの振れ角
−θ max に達することができず，いずれ収束し周期運動が維持できない．そこで両足が接地する時刻 t3 で周期
運動を維持できるような撃力を加える．この撃力は人間の左足による歩行面の蹴り上げに対応する．表 2 に足
踏み歩行時の減衰係数測定結果を示す．個人差もあるが同じ被験者でもばらつきが見られたので，目安とし
て μ が 1.1～1.2 の間におさまるように減衰係数を決定する．
表 2 減衰係数の測定実験結果
被験者 A
被験者 B
被験者 C
被験者 D
平均
μ
1.17
1.36
1.11
1.17
1.03
2.4 人間の横方向加振力 人間が歩行する際，足裏には床反力が生じ，それは進行方向・上下方向・左
右方向にわけて観測することができる．この床反力は通常，体重すべてがかかった時の力（＝体重×重
力 [N]）を 100%として表される．3 方向の中でも左右方向の床反力が，モデル作成に必要な横方向加振
力である．人間の横方向床反力は上下方向加振力に比べて個人差や歩き方に大きく左右されるということは
想像に難くない．実際に論文 3,4)によっても結果がさまざまである．図 6 は横方向床反力の概形である． t1 が片
方の足が着地した瞬間で， t3 がその足が地面を離れる瞬間である．おおよその評価として，歩行モデルによっ
て生じる横方向床反力を平均 10％前後に収める用にパラメータを決定していく．不安定な床面を歩く際は少し
足幅を広げて歩くため一般的な歩行より少し高めに目標とする床反力を設定した．
2.5 作成モデルまとめ 図 7 に作成したモデルの挙動をまとめた．表 3 に本章で調整したパラメータをまと
め，それによって自動的にきまる通常歩行時におけるモデルの性質を表 4 に示した．通常歩行時における実
際の人間の性質をよく模擬していると言える．
%
60 kg
15
0.8 m
10
5
8°
t1
t3
時刻 [s]
-5
0図 6
あがっていくにしたがって、足裏の
支点を外側に変化させていく
衝突後、蹴り上げる力でいつも
同じ初速度であがり始める
-1
人間の片足から生ずる横方向床反力（概
図 7 作成モデルのまとめ
表 3 調整したパラメータ
質量 m
60.0 kg
脚部長さ l
脚の開き角
ϕ0
最大角変位
θ max
減衰係数 C
o
o
1.20 /s
4.00
8.00
表 4 自動的に定まる通常歩行時の性質
μ
初期角速度
歩行リズム f p
横方向床反力 R
0.800 m
29.1
o
s
0.866 Hz
1.14
11.0%
4. 自由振動台上での引き込み現象
4.1 自由振動台上を 1 人の人間モデルが歩行した場合 最大角変位や支点移動機構の差異によって引
き込み周波数域での位相差の特徴にわずかな差をもつモデル A（固有歩調 0.84Hz）に，それぞれ固有振
動数 0.8Hz,0.87Hz の 2 種類の台上を歩行させた．
固有歩調より低い固有振動数をもつ台上に歩かせた場
合(図 8)，自身の横方向加振力により台の振幅が増幅し，それによって台と歩行モデルの振動数が想定
通り同相で同調し，台の固有振動数付近に引き込まれていることが読み取れる．振動数の収束値が台の
固有振動数から少しずれているのは，人が歩行することによる横方向加振力の影響である．また，高い
固有振動数をもつ台上を歩かせた場合(図 9)，振幅があまり増幅することがなく，台と歩行モデルの振
動数は想定していた逆相ではなく同相で同調し，
台の固有振動数よりは自身の固有歩調付近に収束した．
共に同相に同調した原因は，強制変位台上での引き込み周波数域で位相差がやや同相に偏っていたため
と考えられる．
4
4
2
2
0
0
-2
-2
-4
0
5
10
15
20
(1) 角変位の過渡応答
25
30
-4
0
0.02
0.01
0.01
0
0
-0.01
-0.01
10
5
15
20
(2) 台変位の過渡応答
25
時刻 [s]
30
-0.02
0
180
90
90
0
0
-90
-90
5
10
15
20
(3) 位相差の過渡応答
25
時刻 [s]
30
-180
0
0.9
0.9
0.85
台
0.8
5
5
0.8
歩行モデル
0.6
0
5
10
15
20
25
10
15
20
10
15
20
30
時刻 [s]
25
歩行モデル
30
時刻 [s]
25
(3) 位相差の
過渡応答
1
0.7
15
(2) 台変位の
過渡応答
180
-180
0
10
(1) 角変位の
過渡応答
0.02
-0.02
0
5
時刻 [s]
30
時刻[s]
台
0.75
20
(4) 振動数の過渡応答
25
時刻 [s]
30
A(0.84Hz)の歩行
図8 固有振動数0.8Hz 台上でのモデル
0.7
0
5
10
15
20
(4) 振動数の
過渡応答
25
30
時刻[s]
0.87Hz 台上でのモデル
A(0.84Hz)の歩行
図9 固有振動数
5. 結 言
本研究では静止した歩行面での歩行を模擬した歩行モデルを作成し，それを用いて引き込み現象の発生
機構の物理的な説明，さらに実際の橋の上での引き込み現象を説明することを目的とした．はじめに，前額面
の歩行を模擬したコンパス型受動歩行モデルを改良することで人間の通常歩行時に与える横方向加振力を模
擬する歩行モデルを作成した．次に作成した歩行モデルを強制変位する台上で歩行させるシミュレーションを
行い，人の歩行の引き込み現象の特徴を再現することを確認した．その結果，強制変位台上での引き込み現
象の定量的な解析が可能となった．最後に，横揺れする橋上を歩行する状況の再現を目指し，自由振動する
台上でモデルを歩行させるシミュレーションも行った．その際，人間と台の重心位置の位相差が台の揺れに与
える影響を説明し，作成した歩行モデルが自由振動台上を歩行した際の収束値に関する特徴を示した．
文献
(1) 松久寛他，歩道橋の横揺れによる歩行リズムの同期現象， D&D, 2003, CD-ROM 論文集, No.458
(2) T. McGeer, Passive Dynamic Walking, Int. J of Rob. Res., Vol.9, No2, pp, 62-82, 1990.
(3) 高見正利・他, 床反力計による健常者歩行の研究 特に年齢および性別による違いについて,「リハ bi
リテーション医学」誌,第 24 巻,93～101 頁,1987.
(4) Adrienne E. Hunt etc, Inter-segment foot motion and ground reaction forces over the stance
phase of walking, Clinical Biomechanics 16, 592~600, 2001.
環境負荷履歴が再生組織の潤滑機能発現・形成に及ぼす影響
(相対すべり環境と静置環境の相互作用に関して)
工学研究科機械理工学専攻
山本 浩司
Abstract: The physiological and physical characteristics of articular cartilage are constructed and maintained
through complex joint movements, involving in fluctuations of compression, rotation and friction. In the
previous reports, it was shown that the static incubation period after stirring-chamber culture for regenerated
cartilage plays an important role for a maturation of tribological properties in the tissue. However, detail effects
of post-static incubation have not clarified yet. In this study, effects of loading and static incubation period on the
maturation of tribological performances were examined. From the results of this study, it was suggested that the
expression of tribological functions during cartilage regeneration was associated with the initial tribological
stimulations and that the effect was gradually matured in the increase of static incubation period.
Key words: Regenerated cartilage, Tribological property, Loading history, Bio-environment designing
1. はじめに
関節軟骨の再生医療において，生体外で培養した組織を生体内の患部に移植する際，どのような状態
にある組織を移植すべきかということに関する評価指標や，その状態に至るよう培養する設計プロセス
の概念は未だ確立されているとは言い難い．現在提唱されている様々な培養方法は，遺伝子発現量やタ
ンパク質産生量の増加を目指したものであり，量的なレベルで生体関節軟骨に近づけることを目標とし
ている場合が多い．しかしそれらを実際に生体内移植し評価した例は少なく，生体内における効果は未
知である．再生段階や，着目するものによって量的な関係が重要な指標となる場合もあるが，生体の機
能は動的な環境変化を介して，その時間的・空間的履歴を内包していることを念頭に置く必要がある．
例えば，軟骨の主要機能である潤滑特性は，わずかな物理的培養環境の違いによって組織量と必ずしも
相関しなくなることを報告してきた[1]．また，仮に生体関節軟骨とまったく同じ再生軟骨を培養できた
場合でもレシピエント組織との接合が問題になると考えられ，生体組織に近づけることが必ずしも最適
なアプローチであるとは言い切れない．以前の研究において，組織量や構造全体を反映した潤滑機能そ
のものを評価指標とし，相対すべり運動下培養で機能改善が見られた組織は，生体内に移植した際，再
生・修復過程においても機能的に作用することを我々は報告してきた[2]．また昨年度の研究において，
相対すべり運動後に静置培養期間を設けることによって培養軟骨の耐摩耗特性が向上するという結果を
得た[1]．本年度は，より詳細に相対すべり運動下培養と静置培養の負荷期間を変化させ，in vitro 培養軟
(a)
(b)
(c)
(d)
図1 Histological observations of sagittal sections on the chondrocyte-seeded fibroin-sponge scaffolds after
3-month allografting. (a) Chondrocyte containing fibroin-sponge scaffold with 2-week dynamic in stirring
chamber group. (b) Chondrocyte-containing fibroin-sponge scaffold with 2-week static incubation group. (c)
Scaffolds insertion group. (d) Defect only (non-insertion group). Toluidine blue staining of the implanted
scaffold (Scale bar = 1 mm) [2].
骨における潤滑機能の発現・形成に及ぼす影響を検討した．
2. 再生軟骨組織の培養および機能評価法
2.1 軟骨細胞の単離および播種 4週齢の日本白色家兎の大腿骨，脛骨，上腕骨より軟骨片を採取し，
コラゲナーゼ処理によって軟骨細胞を単離した．単離した軟骨細胞を一度継代し，フィブロインスポン
ジ(直径：8 mm，厚さ：1 mm)にマイクロシリンジを用いて5.0×105 cells播種した．
2.2 相対すべり運動負荷培養 以前の研究で報告した相対すべり運動負荷培養装置を用いた[1,2]．図 2
に示すように，細胞播種面を底面に向け(ディッシュ上面に相対する方向)，各組織に対して約 30 mg の
ステンレスメッシュをおもりとして設置した．回転速度は 5 rpm(組織中心のすべり速度：12.0 mm/s)に設
定し，相対すべり運動を負荷する群には 24 時間/day 連続した回転を行った．
2.3 環境負荷履歴の設定 図 3 に各群の培養プロセスを示す[3]．播種後 3 日間は細胞の培養担体への接
着期間(Ad)として静置培養を行った．その後，相対すべり運動負荷培養装置に移し，合計 2 週間の間で相対
すべり運動下培養(D)と静置培養(S)期間の負荷履歴を変化させ，合計 5 群の比較を行った(D0S11，D1S10，D4S7，
D7S4，D11S0)．例えば D1S10 群は 3 日間の Ad の後に 1 日の相対すべり運動を加え，さらに 10 日間の静置培養
を継続したことを意味する．また，D0S11，D11S0 はそれぞれ静置培養のみ，相対すべり運動下培養のみを負荷
した単一培養条件群である．各培養において CO25%，37°C であり，静置培養には 24well 培養プレートを用
いた．
2.4 潤滑特性評価 潤滑特性は摩擦時間に対する動摩擦係数の変化を各群の培養2週間目の組織に対し
て測定した．摩擦実験には，往復動摩擦試験機を用い[3]，ステンレス平板(Ra=0.06 µm)に対して，生理食
塩中にて摩擦速度 0.8 mm/s，垂直荷重 29 mN の条件で測定した．なお荷重および摩擦力は，ワイヤー放電に
よって加工された平行平板ばねにより設定・計測し，荷重および摩擦速度はステッピングモータを用いて制御
した．
2.5 生化学的評価 再生軟骨組織の成熟度を評価する指標の一つとして，軟骨組織の細胞外基質産生量
が用いられている．各試料において，再生軟骨組織中の glycosaminoglycan をコンドロイチン硫酸量を測
定することによって評価した．
3. 実験結果
図 4 に示すように，期間に関係なく，相対すべり運動を負荷した後に静置培養へ移行した群(D1S10，D4S7，
D7S4)は他群に比べて摩擦初期の摩擦係数が低い値を示す傾向にある．一方，相対すべり運動を継続して負
荷した D11S0 群は，摩擦の初期段階から静置培養群よりも高い摩擦係数を示した．長期の摩擦特性を比較する
と，D4S7，D7S4 群が他群に比べて低い摩擦係数を維持しているが，D1S10 群に関しては徐々に上昇する傾向が
図 2 Schematic diagrams of the culture method
using the lab-made stirring chamber.
図 3 Design of loading history. Ad: Cell adhesion
period. Di: Culture with the stirring chamber. Sj: Static
incubation with 24-well plate. (i,j: culture period
[day]).
図 4 Time-dependent transitions of friction coefficient
of regenerated cartilage. (n=3, Mean±SD.)
図 5 Synthesis of extracellular matrix evaluated by
glycosaminoglycan content. (n=3, Mean±SD., *:
p<0.05)
あった．また，単一の環境で培養し続けた D0S11 および D11S0 群は，20 分の摩擦でフィブロインスポンジのみの
摩擦係数と同等の値を示すまで上昇した．図 5 は組織内コンドロイチン硫酸量の比較である．D1S10 群が他群
に対して有意に高い値を示した．しかし，摩擦特性に顕著な差が現れた D4S7，D7S4，D11S0 群における有意差
は見られなかった．
4. 考察
生体関節軟骨の潤滑においては，関節液や組織表層に存在する糖タンパク質およびリン脂質が重要な役
割を担っていることが報告されている[4,5]．特に近年では lubricin や superficial zone protein(SZP)として知られ
る proteoglycan4(PRG4)に関連した物質がヒアルロン酸と協調して低摩擦機構を達成することなどが示唆されて
いる[6]．また PRG4 の発現には，動的なせん断刺激負荷の関連性も報告されていることから[7]，相対すべり運
動群は静置群より潤滑に寄与する物質が表層部に産生されていた可能性が高い．しかし，組織形成過程に継
続的な相対すべり運動を負荷すると，初期・長期の摩擦特性が悪化する結果となった．このことから，本培養系
では，すべり運動による組織表層における物理的な摩耗と，組織の摩擦低減や耐摩耗改善に寄与する物質産
生および構造構築への効果が混在していると考えられる．相対すべり運動を 1 日しか負荷しなかった D1S10 群
においても初期の潤滑機能改善が見られたことから，細胞レベルへの効果は運動負荷の比較的早い段階で
達成されたと考えられる．しかし耐摩耗特性の改善にはある程度の期間相対すべり運動を負荷することが必要
であると言える．
図 7 は摩擦特性に改善の見られた D4 群に着目して静置培養期間と摩擦特性の変化を調べた結果である．
図 6 Time points when frictional experiments were
performed during static incubation of D4 group.
図 7 Time-dependent transitions of friction coefficient
of regenerated cartilage in D4 group. (n=3,
Mean+or-SD.)
(図 6 は摩擦測定のタイムポイントを示す．) 図 7 に示すように，相対すべり運動を負荷し終わった直後(D4S0)
は負荷前(D0S0)と変わらない摩擦特性を示すが，静置培養の日数に応じて摩擦特性が変化していくことがわか
る．特に初期の摩擦係数変化にその傾向が強く現れている．
これらの結果から，軟骨細胞が組織化するプロセスにおける潤滑機能形成を考えると，初期のある程度の刺
激(量，期間)およびその効果に伴って生じる状態(遺伝子発現，タンパク質産生など)の成熟期間が重要である
と考えられよう．
5. 結言
すべり環境の負荷履歴が軟骨再生過程の潤滑特性に及ぼす影響を調べた結果，比較的早期，短期間の
相対すべり運動負荷によって軟骨基質産生が増加するが，摩擦および耐摩耗構造の構築のためには，負荷
の期間とその効果が発現する培養環境に対するより詳細な制御が必要であった．今後潤滑に寄与する遺伝
子などの発現履歴を調べ，移植に適した組織の培養環境設計に繋げる予定である．
参考文献
[1] 山本浩司，環境負荷履歴が再生組織の潤滑機能発現・形成に及ぼす影響，21 世紀 COE プログラム 平成
18 年度年次活動報告書，(2007)，369-372.
[2] Chueh S., Tomita N., Yamamoto K., Harada Y., Nakajima M., Terao T., Tamada Y., Transplantation of allogenic
chondrocytes cultured in fibroin sponge and stirring chamber to promote cartilage regeneration, Tissue Eng., 13,
(2007), 483-492.
[3] 山本浩司，高家理気，川上雅弘，玉田靖，富田直秀，相対すべり運動の負荷履歴が軟骨再生過程におけ
る潤滑機能形成に及ぼす影響，日本トライボロジー学会予稿集，(2007)，521-522.
[4] Pickard J.E., Fisher J., Ingham E., Egan J., Investigation into the effects of proteins and lipids on the frictional
properties of articular cartilage, Biomaterials, 19, (1998), 1807-1812.
[5] Jay GD., Cha CJ., The effect of phospholipase digestion upon the boundary lubricating ability of synovial fluid, J
Rheumatol., 26(11), (1999), 2454-2457.
[6] Schmidt TA., Gastelum NS., Nguyen QT., Schumacher BL., Sah RL., Boundary lubrication of articular
cartilage: role of synovial fluid constituents, Arthritis Rheum., 56(3), (2007), 882-891.
[7] Nugent GE., Aneloski NM., Schmidt TA., Schumacher BL., Voegtline MS., Sah RL., Dynamic shear stimulation
of bovine cartilage biosynthesis of proteoglycan 4, Arthritis Rheum., 54(6), (2006), 1888-1896.
動的機能を有する機械構造システムの創成設計法の開発
工学研究科航空宇宙工学専攻 西脇 眞二
Abstract: In vibration optimization problems, eigen-frequencies are usually maximized in the optimization since
resonance phenomena in a mechanical structure must be avoided, and maximizing eigen-frequencies can provide a
high probability of dynamic stability. However, vibrating mechanical structures can provide additional useful
dynamic functions or performance if desired eigen-frequencies and eigen-mode shapes in the structures can be
implemented. In this research, we propose a new structual optimization method for designing vibrating structures
that targets desired eigen-frequencies, using level set-based structural optimizations. Several numerical examples
are presented to confirm that the method presented here can provide optimized vibrating structures applicable to
the design of mechanical resonators and actuators.
Key words: Level Set-based Structural Optimization, Vibrating Structures, Eigen-Frequency, Eigen-Mode
1. はじめに
一般の機械構造は，破壊や共振現象を回避し，システムの安定性を確保するため，剛性の最大化や固
有振動数の最大化を目標として設計されることが多い．しかし，機械構造に動的な柔軟性や，特定の固
有振動数と固有モードを有する動的特性を適切に付加することが可能となれば，全く新しい機能を有す
る複雑な機械システムを設計することができる．本研究では，これらの新しい機械構造システムを創成
する方法をトポロジー最適化[1]の考え方に基づき構築した．そして，現在までに，それらの方法により，
特定の振動特性をもつ機械式レゾネータ[2]の創成設計法を開発した．さらに，これらの最適化の考え方
をマルチフィジックス現象に展開し，ピエゾ電気材料を利用したアクチュエータ[3]，ピエゾ抵抗効果を
利用したセンサー[4]，マイクロ波導波管の断面形状設計[5]，誘電体共振アンテナ[6]の設計法，フォノニ
ックバンドキャップ構造の形状創成設計法[7]を開発してきた．しかし，トポロジー最適化は，動的な目
標関数の設定が複雑になると，目標関数の多峰性によりグレースケールを多く含む物理的に意味のない
最適構造を創出する欠点をもつ．他方，グレースケールなどの数値的問題を抜本的に解決しながら構造
の形状と形態の変更を可能とする新しい形状最適化の方法として，レベルセット法に基づく構造最適化
[8,9]が提案されている．ここでは，特定の振動特性をもつ機械式レゾネータの設計を目的に，振動特性
の適正化を目指した構造最適化法を，レベルセット法よる形状表現に基づき構築する．
2. レベルセット法による形状表現
物体により占められている領域 Ω+ (以下，物体領域と記述する)と，物体により占められていない領域
Ω − (以下，空洞領域と記述する) の混在が許容される固定設計領域 D において，物体領域の形状を最
適化する問題について考える．レベルセット法では，物体境界からの符号付き距離関数を，レベルセッ
ト関数 φ ( x ) として導入し，物体領域 Ω+ 及びその外形形状 ∂Ω を，レベルセット関数を用いて表現する．
すなわち，次式に示すように，物体領域 Ω+ において正，物体の境界 ∂Ω ではゼロ，空洞領域 Ω− において
負となるレベルセット関数を用いて物体の形状表現を行う．
⎧φ ( x ) > 0 ⇔ x ∈ Ω +
⎪
⎨φ ( x ) > 0 ⇔ x ∈ ∂Ω
−
⎪
⎩φ ( x ) < 0 ⇔ x ∈ Ω
(1)
さらに，次式に示すように，物体領域で 1，空洞領域で 0 の値をとるヘビサイド関数 H (φ ( x ) ) を導入す
る．
⎧⎪0 if φ ( x ) < 0
H (φ ( x ) ) = ⎨
⎪⎩1 if φ ( x ) ≥ 0
(2)
そして，最適化過程における物体領域の変動を，レベルセット関数の変動により表現する．今，レベル
セット関数の変動を移流とみなし，レベルセット関数の移流速度をV ( x ) とすると，次の移流方程式を得
る．
∂φ
+ V ( x ) ⋅ ∇φ ( x ) = 0
(3)
∂t
ここで，移流速度V ( x ) のレベルセット関数の等位面に対する法線成分を VN ( x ) とすると，
VN ( x ) = V ( x ) ⋅
∇φ ( x )
(4)
∇φ ( x )
であるから，式(3)に代入して，次式を得る．
∂φ
+ VN ( x ) ∇φ ( x ) = 0
(5)
∂t
上式の VN ( x ) を，物体領域の変動を決定する外形形状感度に設定すれば，上式を解くことにより最適構
造を得ることができる．さらに，レベルセット関数が符号付距離関数であれば，∇φ ( x ) = 1 であるから，
上式は次式のように簡略化される．
∂φ
+ VN ( x ) = 0
(6)
∂t
3. 最適化問題の定式化
固定設計領域 D において，制約汎関数 Gφ の制約のもと，目的汎関数 Fφ を最小化する次式の最適化問
題を考える．
Minimize Fφ
(7)
Subject to Gφ − Gmax ≤ 0
(8)
ここで， G max は Gφ の上限値である．さらに，目的汎関数の変分 δFφ と制約汎関数の変分$は，レベル
セット関数の任意の変分 δφ を用いて，次式のように表されるものとする．
δFφ = ∫ f (φ (x ))δφdΩ
(9)
δGφ = ∫ g (φ (x ))δφdΩ
(10)
Fφ = Fφ + λ (Gφ − Gmax )
(11)
D
D
ここで， f (φ (x )) と g (φ (x )) は，それぞれ目的汎関数と制約汎関数の変分を記述するためのφ (x ) の関数
である．次に，上の最適化問題の KKT 条件を導く．まず，ラグランジェアン Fφ をラグランジュ乗数 λ
を用いて，次式のように表す．
上式を用いて KKT 条件を導けば，次式となる
Fφ = 0, λ ( Gφ − Gmax ) = 0, λ ≥ 0, Gφ − Gmax ≤ 0
(12)
上式の Fφ = 0 を展開すれば
(
)
(
)
f φ ( x (t ) , t ) + λ g φ ( x (t ) , t ) = 0
(13)
が得られ，この式を満足するレベルセット関数が最適解となるが，これを直接求めるのは困難である．
そこで，仮想的な時間 t を導入し，次式に示す偏微分方程式を解くことにより最適解を求める．
∂φ (x(t ), t )
+ f (φ (x(t ), t )) + λg (φ (x(t ), t )) = 0
∂t
(14)
式(6)と式(14)を比較すれば，本最適化問題においては，VN ( x ) は次式となることがわかる．
(
)
(
VN ( x ) = f φ ( x ( t ) , t ) + λ g φ ( x ( t ) , t )
)
(15)
さらに，もし式(8)の不等式制約が非活性れあれば, λ は 0 となり，活性であれば
れる．
λ=−
∫ f (φ (x(t ), t ))g (φ (x(t ), t ))dΩ
∫ {g (φ (x(t ), t ))} dΩ
D
λ は次式により得ら
(16)
2
D
このようにラグランジュ乗数を決定し，式(14)に基づきレベルセット関数を初期値から更新していき，
変動が収束したときのレベルセット関数が，最適解となる．
4. 目的関数と制約条件の設定
境界 Γu を完全拘束した固定設計領域 D において，線形弾性体が自由振動している場合を考える．こ
のときの最低次から k 番目の固有振動数ω k ，固有振動モードを u k と表せば，次式に示す固有方程式が
成立する．
−ωk 2 ∫ ρ u k ⋅ u k H (φ ( x ) ) d Ω + ∫ ε ( u k ) : D : ε ( u k ) H (φ ( x ) ) d Ω = 0
D
D
(17)
ここで，ε はひずみテンソルであり，微小線形ひずみを仮定している． D と ρ は，それぞれ，線形弾性
体の弾性定数テンソルと質量密度である．また，H (φ (x )) はヘビサイド関数であり，次式に示すように，
物体領域において 1，空洞領域において 0 をとる．
⎧⎪0 if φ ( x ) < 0
H (φ ( x ) ) = ⎨
⎪⎩1 if φ ( x ) ≥ 0
(18)
ρu k ⋅ u k H (φ (x ))dΩ = 1
(19)
さらに，通常固有振動モードは，
∫
D
となるように正規化するので，k 番目の固有振動数ω k は
ωk 2 = ∫ ε ( u k ) : D : ε ( u k ) H (φ ( x ) ) d Ω
(20)
D
として与えられる．
まず，より動的に安定な構造物を得るために，最低次の固有振動数の最大化を考える．このとき目的汎関数
Fφ は次式となる．
Fφ = −ω1
2
(21)
なお，この場合，目的汎関数の最大化を最小化に置き換えるため， ω12 に負符号を付加し，目的汎関数
としている．この場合，式(9)の f (φ (x )) は次式で与えられる．
f (φ ( x ) ) = − {ε ( u1 ) : D : ε ( u1 ) − ω12 ρ u1 ⋅ u1}
dH
dφ
(22)
次に，特定の振動数において共振現象を示す機械構造物を得るために，s 番目の固有振動数 ωs を目標
とする値ω ojb, s に一致させる問題を考える．このとき，目的汎関数 Fφ は次式のように定式化される．
(
Fφ = ω s − ω obj , s
2
)
2 2
(23)
この場合，式(9)の f (φ (x )) は次式で与えられる．
(
f (φ (x )) = 2 ω s − ω obj , s
2
2
){ε(u ) : D : ε(u ) − ω
s
s
2
s
ρu s ⋅ u s
}dH
dφ
(24)
他方，制約汎関数 Gφ $は，上の二つの場合について物体領域の体積とする．すなわち，
Gφ = ∫ H (φ (x ))dΩ
(25)
D
この場合，式(10)の g (φ ( x ) ) は次式で与えられる．
g (φ (x )) =
dH
dφ
(26)
5. 数値例
幾つかの平面応力状態を想定した二次元の振動特性最適化問題を解くことにより，本研究で提案する
最適化法の有用性を検証する．
図１に設計領域と境界条件を示す．
設計領域の左辺と右辺を完全拘束し，
中央部分に 0.16kg の質点を設定する．図２に初期構造を示す．図中灰色の領域が物体領域，白色の領域
が空洞領域である．設計領域を構成する要素には双一次近似アイソパラメトリック四角形平面応力要素
5
3
を用いる．また，図の物体領域は，ヤング率 2.1× 10 MPa，ポアソン比 0.3，質量密度 7.8 ×10 kg/m3
−5
の線形弾性体で構成されているものとする．さらに，物体領域の体積の上限値 Gmax を1.5 × 10 m2 と
し，最適化を図る．なお，初期形状における最低次の固有振動数は 1189Hz，二番目の固有振動数は 4077Hz，
三番目の固有振動数は 14052Hz であった．
図１ 設計領域と境界条件
図２ 初期構造
最初に，最低次の固有振動数を最大化する場合について，最適構造を求める．図３に最適構造を示す．
図３(a)は，トポロジカルデリバティブを用いて形態変更を行った場合で，このときの最低次の固有振動
数は 2055Hz となった．図３(b)は，形態変更を行わなかった場合で，このときの最低次の固有振動数は
1432Hz となった．図３(a)より，最適構造は，設計領域に二つの孔をもち，y 方向に幅もつ構造が創出し
ていることがわかる．この場合の最低次の固有振動モードは，設計領域の中心に腹をもつ y 方向の一次
曲げモードであるが，この構造により腹部分の剛性が向上し，最低次の固有振動数は高くなっていると
考えらる．また，図３(b)より，形態変更を行わなかった場合には，形態変更を行った場合より腹部分の
剛性が向上しないため，最低次の固有振動数は，形態変更を行った場合ほど高くならなかったと考えら
れる．
(a) 形態変更を行った場合
(b) 形態変更を行わなかった場合
図３ 最低次の固有振動数を最大化したときの最適構造
次に，設定した特定の固有振動数をもつ構造を求める．図４に，最低次の固有振動数を特定の値に設
定した場合の最適構造を示す．図４(a), (b), (c), (d), (e)は，最低次の固有振動数をそれぞれ 300Hz，600Hz，
900Hz，1500Hz，1800Hz に設定した場合である．なお，この場合に，最適構造の最低次の固有振動数は
それぞれ，300Hz，600Hz，901Hz，1496Hz，1799Hz となった．
(a) 最低次の固有振動数の設定値 = 300Hz
(b) 最低次の固有振動数の設定値 = 600Hz
(c) 最低次の固有振動数の設定値 = 900Hz
(d) 最低次の固有振動数の設定値 = 1500Hz
(e) 最低次の固有振動数の設定値 = 1800Hz
図 4 最低次の固有振動数を目標値に設定した場合の最適構造
これらより，どの固有振動数の場合においても，物理的に意味のある明瞭な構造が得られていること
がわかる．そして，最低次の固有振動数を 300Hz，600Hz，900Hz に設定した場合には，体積制約は非
活性となっているうえ，形態の変更はない．これに対して，最低次の固有振動数を 1500Hz，1800Hz に
設定した場合，体積制約は活性となっているうえ，構造物の形態は変更されている．すなわち，最低次
の固有振動数について，比較的低い値に設定した場合には，固有振動数の変更に体積の増加を必要とし
ないが，形態変更を行わなかった場合の最大値である 1432Hz を越える設定では，形態変更と体積の増
加を必要とすることがわかる．
6. まとめ
本研究では，動的機能を有する機械構造システムの創成設計法の開発する一環として，特定の振動特
性をもつ機械式レゾネータの設計を目的に，レベルセット法にも基づく構造最適化の方法により，振動特性の
適正化を目指した構造最適化法を構築した．併せて，簡単な数値例により，本研究で提案する方法の有用性
を検証した．その結果，提案する方法により，指定の固有振動数をもつ明瞭な最適構造が得られることがわか
た．
参考文献
[1] Bendsøe, M. P. and Kikuchi, N., Generating optimal topologies in structural design using a homogenization
method, Comput. Meth. Appl. Mech. Engrg., 71, (1988), 197-224.
[2] Maeda, Y., Nishiwaki, S., Izui, K, Yoshimura, M., Matsui, K, and Terada, K., Structural Topology Optimization
of Vibrating Structures with Specified Eigen-frequencies and Eigen-mode Shapes, Int. J. Numer. Meth. Eng., 67,
(2006), 597-628.
[3] Carbonari, R.C., Silva, E. C. N., and Nishiwaki, S., Design of piezoelectric multiactuated microtools using
topology optimization, J. Smart Struct. Mat., 14, (2005), 1431-1447.
[4] Rubio, W. M., Silva, E. C. N., and Nishiwaki, S., Piezoresistive sensor design using topology optimization,
Structural and Multidisciplinary Optimization, to appear.
[5] 木下慎也，西脇眞二，泉井一浩，吉村允孝，野村壮史，佐藤和夫，平山浩一，周波数特性を設計目標とす
る電磁波共振器のトポロジー最適化，日本計算工学会論文集，2007, No. 20070025.
[6] Nomura, T., Sato, K., Taguchi, K., Kashiwa, T., and Nishiwaki S., Structural topology optimization for the
design of broadband dielectric resonator antennas using the finite element time domain technique, Int. J. Numer.
Meth. Eng., 71, (2007), 1261-1296.
[7] Mogami, K., Nishiwaki, S., Izui, K., Yoshimura, M., and Nomura, T., Structural Optimization Using Discrete
Elements for the Design of Structures Implementing Particular Wave Propagation Functions, 7th World Congress on
Structural and Multidisciplinary Optimization(WCSMO-7), May 21-25,2007, Seoul, Korea, 2033-2042.
[8] Wang, M. Y., Wang, X. and Guo, D., A level set method for structural topology optimization, Comput. Meth.
Appl. Mech. Engrg.., 192, (2003), 227-246.
[9] Allaire, G. and Jouve, F., A level-det method for vibration and multiple loads structural optimization, Comput.
Meth. Appl. Mech. Engrg., 194, (2005)., 3269-3290.
人間機械共働系としての生産システムにおける改善の研究
工学研究科機械理工学専攻
水山 元
Abstract: It is studied how autonomous small groups running kaizen activities in a complex manufacturing system
should be managed, so as to obtain not only their ordinal intangible effects but also coherent overall effects on the
tangible measures like costs, quality, lead-times, etc. from a number of kaizen tasks conducted by the distributed
kaizen teams. First, possible options are formulated for the organizational structure of the small groups, the budget
system for the kaizen tasks, and the performance evaluation system of the groups. Then, communications among
the kaizen teams are modeled as algorithmic negotiation processes. Finally, an agent-based simulation model is
developed to investigate how these institutional design elements on kaizen team management affect the whole
performance of the kaizen process. Numerical examples of the simulation are also provided to show how the
approach works.
Key words: Agent-based simulation, Kaizen activities, Autonomous small groups, Institutional design
1. はじめに
生産システムは，複数の人間と機械設備が協働する系であり，競争力の維持・向上，収益性の向上，
環境変化への対応などのために，常に変革を迫られている．そこでは，人間は，客体としての生産シス
テムの構成要素であると同時に，それに変革を加える主体でもある．一方，それを取り巻く市場環境は，
当該生産システムの他，顧客，競合他社などの複数の意思決定主体から構成される相互作用系である．
すなわち，生産システムは，それ自身の内部環境も，それを要素として含む外部環境も，ともに多主体
系であるという特徴を持っている．
本研究では，
これらの特徴を考慮して，
「変革の主体としての人間に，
対象生産システムの問題点や改善すべき方向を示すインタフェースの解明と設計」
，
「内部環境の多主体
性を考慮した生産システムとその変革のモデル化」
，
「外部環境の多主体性を考慮した改善・革新の評価
フレームワークの構築」などの課題に取り組んでいる．
今日の改善活動には，作業者の士気や能力を向上させる，現場の課題を表面化しやすくする，改善の
PDCA サイクルを高速化する，などのインタンジブルな効果に加えて，さらに，原価低減，品質向上，
納期短縮などのタンジブルな効果を生み出していること，そして，それらの効果が，単に全方位的なも
のではなく，当該企業の競争戦略に合致したものであること，が求められている．暗黙知を反映させた
ローカルな意思決定をベースにして自律分散的に進められる改善小集団活動において，インタンジブル
な効果だけでなく，そうしたタンジブルな面でも，総体として首尾一貫した効果が得られるようにする
ためには，その制度設計・環境設計を適切に整えることが重要である．以下では，そうした観点から，
改善小集団の組織構成，改善予算制度，業績評価制度などを検討するための新たなアプローチを開発し
た成果の概要を報告する．
2. 改善小集団活動の制度設計問題
改善とは，
「ロス，ムダなどと一般に呼ばれる，現場の課題（以下では，まとめてロスと呼ぶ）に対し
て，適切な対策を施し，それを解消あるいは低減すること」と捉えることができる．生産システムに生
じ得るロスの分類については，例えば，文献[1-3]などにまとめられており，また，企業がその業態に応
じて個別に定義していることも多い．任意のロスを解消あるいは低減した場合，結果として生じ得るタ
ンジブルな効果としては，
「余分な資源消費が抑制されることによる原価の低減，
在庫が減少することに
よる資本金利の低減，などの金銭面での効果」
，
「不良率が低下することなどによる品質面での効果」
，
「リードタイムが短縮されること，安定することなどによる納期面での効果」などが挙げられる．本研
究では，これらの多面的な効果は金銭尺度にまとめて一元化されているものと仮定し，その一元化され
た尺度を（広義の意味で）
生産技術部隊 E
エージェントシミュレーションモデル
「製造コスト」と呼ぶこと
各小集団の意思決定，
改善提案
予算・技術サポート
小集団間のコミュニケ
にする．そして，その製造
ーション等をモデル化
コストに基づいて，各改善
小集団がそれぞれ疑似プロ
現場小集団 g
現場小集団 g’
現場小集団 g”
フィットセンタとして評価
されている状況を考える．
ロス発生状況
改善活動計画
各小集団は，他小集団と
生産システムシミュレーションモデル
のコミュニケーションと，
改善策実施による効果，
自小集団の意思決定とを重
各種ロスの発生状況等
を評価
ねながら，それぞれ自小集
団の利益増大を目指して改
図 1 提案フレームワークの概要
善活動を進めていく．問題
は，そうした自律分散的な
改善活動の総体が，結果として，前述の広義の製造コストを効果的に減少させていくように，環境を適
切に整えることである．具体的には，例えば，次のような項目を検討する必要がある：
(1) 生産システムの現場への分割の仕方：これによって，任意の現場小集団が把握可能なロスが異なり，
それらと他の現場のロスとの関係なども異なる．
(2) 広義の製造コストの定義とその各現場への配賦基準：これによって，同じ改善でも，その効果額，
効果の及ぶ範囲などが異なる．
(3) 技術者を，それぞれいずれかの現場小集団に専属とするか，あるいは，まとめて，どの現場も支援
できる独立した小集団を形成するとするか．
(4) 技術者が独立型の場合，予算を現場ごとの個別枠とするか，全体枠として技術者小集団が配分する
か．また，現場小集団から技術者小集団への報酬をどのように定義するか．
(5) 改善による製造コスト低減額を現場の標準（振替価格）に反映させるまでの期間はどの程度がよい
か．
3. エージェントシミュレーションの展開
図 1 は，改善活動の制度設計について検討するための，提案フレームワークの全体像を示したもので
ある．これは，
「各改善小集団の意思決定，小集団間のコミュニケーションなどをモデル化するエージェ
ントシミュレーションの部分」と「改善策実施による効果，各種ロスの発生状況などを評価する生産シ
ステムシミュレーションの部分」とから構成される．ここでは，それらのうち，上部モジュールのエー
ジェントシミュレーションの基本設計の概要について述べる（詳細は文献[4]を参照）
．
提案法では，個々の現場小集団 g（ G ）をそれぞれ一つの意思決定エージェントとしてモデル化し
た．任意の現場エージェント g は，改善案 i（ L ）の好ましさを，毎 t 期，ある尺度で比較しており，
~
その選好尺度V gtiG の値の大きい改善案から優先的に開始しようと行動する．また，技術者が独立型の場
合は，技術者小集団も一つの意思決定エージェントであると考える．技術者エージェントも個々の改善
案 i に対して独自の選好尺度VtiE をもっており，その尺度の大きい改善案から優先的に開始しようと行動
する．なお，各意思決定エージェントの，任意の改善案に対する選好尺度は，それによって当該エージ
ェントにもたらされる総利益額を制約リソース（技術者の工数あるいは予算）の消費量で割ったもので
ある．
技術者の組織構成が専属型の場合は，技術者エージェントは存在しない．また，技術者の工数と予算
の双方のリソースが，すでに現場エージェントの間に分配されているため，現場間でリソースの取り合
いも生じない．これらのことから，この場合，各期に開始する改善を決定する問題は，現場ごとに独立
した 2 次元ナップサック問題となる．一方，技術者の組織構成が独立型の場合は，技術者エージェント
が技術者の工数を（予算が全体枠ならば，さらに予算配分の決定権も）所持しているため，現場エージ
ェントが各期に新しい改善を開始するためには，そ
~
V gtiG
れに必要な工数（および予算）を技術者エージェン
トから提供してもらう必要がある．このとき，技術
この改善案が
選択される。
者エージェントは技術者サイドの選好尺度VtiE の高
いものから優先的に工数（や予算）を割り当てたい
y5
と考え，各現場エージェント g は，自現場の改善案
~
で，現場サイドの選好尺度V gtiG の高いものから優先
的に工数（や予算）を割り当てて欲しいと考える．
それらの優先順位は一般には一致するとは限らない
ので，期に一度，各エージェントが集まった会議な
どを通じて，工数（と予算）の配分について交渉す
VtiE
x5
るような状況を考える．そうした交渉は，通常，そ
れぞれが順に少しずつ妥協していって，最終的に実
図 2 現場と技術者の交渉
行可能解に到達するようなプロセスになる．ここで
は，それを，次のような，簡単な交渉アルゴリズム
によってモデル化した：
~
Step 1: 当該期に各現場で開始する改善の集合を S gtIni   g  G に初期化し，開始する改善の候補を
  L  StRun ，アルゴリズムのパラメータを m = n = 1 とおく．また， 集合  に含まれる全要素 i
~
を，図 2 のように，横軸 x に技術者サイドの選好尺度VtiE ，縦軸 y に現場サイドの選好尺度V gtiG をと
った xy 平面にプロットする（複数の現場で把握可能なロス i については，それらの各現場に対応す
る複数の点をプロットする）
．そして，プロットされた点の x 座標（ y 座標）を大きいものから順
に，x1, x2, …（ y1, y2, …）と名付ける．
Step 0: 第 t 期に残されている技術者の工数あるいは予算に新たな改善を開始する余裕がなくなった場
合は終了．
~
Step 1: 集合  にVtiE  xm かつVgtiG  yn を満たすロス i があれば，それを選択し，     {i} とおく．
~
~
Step 2: 技術者の工数および予算の双方で実行可能であれば， S gtIni  S gtIni  {i} とする．また，必要に応
じて，集合  を更新する．
Step 3: n > m であれば，m = m+1，そうでなければ，n = n+1 として Step2 へ戻る．
4. 数値実験
簡単な数値例で，提案法の基本的な機能を確認するとともに，制度設計が全体的な改善効果に与える
影響の一端を実際に考察する．数値例では，現場は三つで，各現場に 100 個ずつのロスが存在し，第 1
期には，それらの程度はすべて hi  1 であるとした．生産量は各期とも 10000 個で，固定費と変動費は
それぞれ次式で与えられるとした：
300
c V  c 0V   hi  c iV
(1)
i 1
300
c F  c 0F  10000   hi  c iF
(2)
i 1
ただし，個々の ciV ，ciF の値は，[0, 1]の一様乱数で与えた．また，任意の改善 i は，ロス i の程度 hi を半
減させるものとし，簡単のため，改善に必要な期間はすべて 1 期で，投資額は製造コストと比べて無視
できる程度であるとした．さらに，技術者は 3 人で，改善予算は十分にあるとした．
ロス間の連関の影響をみるために，改善 i によってロス i 以外のロス j に及ぶ影響について，次の三つ
の場合を考えた：
Case1: ロス間に連関は無い．任意の改善 i はロス i にのみ影響を与える．
Case2: ロス間に中程度の連関がある（改善 i が別のロス j に影響を与える組合せは 3000 組）
．
Case3: ロス間に大程度の連関がある（改善 i が別のロス j に影響を与える組合せは 9000 組）
．
ただし，ロス間の連関によって間接的な影響が及ぶ組
合せはランダムに選択し，それによって程度が変化す
る倍率は[0.5, 1.5]からの一様乱数で与えた．また，制
度設計の基本パターンとしては，次の三つを取り上げ
た：
Pattern1: 技術者が専属型の場合
Pattern2: 技術者が従属型で，成果報酬型報酬の場合
Pattern3: 技術者が従属型で，時間差型報酬の場合
なお，改善予算は十分あると仮定したので，予算の全
体枠と個別枠の区別はしていない．また，どの改善の
必要工数も等しいと仮定したので，報酬が工数比例型
の場合は考えていない．
製造コストのうち，ロスの影響を受ける部分のみを
取り出すと次式となる：
表 1 改善効果の平均と信頼区間
Case
1
2
3
Pattern
平均
95%信頼区間
1
0.38
[0.37, 0.40]
2
0.38
[0.37, 0.40]
3
0.38
[0.37, 0.40]
1
0.36
[0.27, 0.44]
2
0.36
[0.29, 0.42]
3
0.26
[0.22, 0.29]
1
0.54
[0.17, 0.91]
2
0.48
[0.14, 0.82]
3
0.16
[0.10, 0.22]
300
C H  10000   (hi  c iV  hi  c iF )
(3)
i 1
ここでは，第 1 期の CH の値に対する 100 期後の CH の値の比率を，改善活動の総体としての効果の指標
と考える．表 1 は，各組合せで 10 回ずつシミュレーションを実行して得られた改善効果の平均と 95%
信頼区間を示したものである．表より，ロス間に連関がない場合は，制度設計による改善効果の違いは
ほとんどみられないが，連関が強くなるに従って，違いが大きくなっていることがわかる．ロス間の連
関が大きい場合に時間差型のパターンが良好な成績を収めているのは，技術者エージェントがロス間の
連関による間接的な効果を考慮できているからである．この影響を定量的に評価できているところに提
案法の一つの意義がある．なお，この結果から，単純に，時間差型を採用すれば問題が解決されると結
論づけることは早計である．ロス間の連関が大きい場合に，時間差型のみが良好な成績を収めているこ
とは，時間差型の技術者の選好尺度と，現場の選好尺度との間の矛盾が大きかったことを示している．
このようなことが続くと，現場の意欲がそがれる恐れがある．むしろ，ロス間の連関が大きい場合にも，
現場が選好尺度を適切に評価できるような，制度設計上の工夫が必要である．例えば，生産システムの
現場への分割の仕方，製造コストの現場への配賦基準，などの面で工夫が考えられる．ただし，具体的
な施策は生産システムの形態に依存すると思われる．したがって，提案したエージェントシミュレーシ
ョンモデル（上部モジュール）に，検討対象の生産システムに対応した生産システムシミュレーション
モデル（下部モジュール）を組み合わせることで，生産システムの形態に応じた詳細な検討を進めるこ
とが今後の課題となる．
5. まとめ
本報告では，自律分散的な小集団によって進められるボトムアップ的な改善活動のための制度設計を
システマティックに検討するための手段として，エージェントシミュレーションによるアプローチを新
たに提案するとともに，簡単な数値例でその基本的な機能を確認した．
参考文献
[1]
[2]
[3]
[4]
大野耐一，トヨタ生産方式，(1978)，ダイヤモンド社．
鈴木徳太郎，生産革新のための新 TPM 展開プログラム 装置工業編，(1992)，JIPM．
中嶋清一，白瀬国夫，生産革新のための新 TPM 展開プログラム 加工組立編，(1992)，JIPM．
水山 元，長尾 昂，自律分散的小集団による改善活動のための制度設計―エージェントシミュレー
ションによるアプローチ―，(2007)，日本設備管理学会誌，to appear．
マルチエージェント系における競争と協調の意思決定
工学研究科機械理工学専攻
田 雅杰
Abstract: A multi-agent system (MAS) is usually used to solve the cooperated problems in many fields. Although
we know that a MAS without competition of autonomous agents is not an efficient system, in the case only limited
resources are available the strong competitive intention of agents can also cause a conflict which makes a MAS fall
into some deadlock situations. Thus, a certain extent of cooperation among autonomous agents is also essential for
avoiding such a conflict. However, it is difficult to decide when the agent needs to make a decision of competition
or cooperation and to what extent of cooperation the agent needs to have. In this research, an idea is obtained from
the Drum-Buffer-Rope (DBR) method described in Theory of Constraints (TOC) for determining a behavior of
competition or cooperation for an agent. An approach called Three-Zone Buffer Management based on DBR is
proposed for making decisions of the opposing behaviors of competition and cooperation. Then, the approach is
applied to many engineering systems such as logistics and production systems.
Key words: Multi-agent system (MAS), Cooperation, Competition, Cooperation, Three-zone Buffer Management
1
はじめに
マルチエージェントシステム(MAS)の考え方はさまざまな分野で応用されている．例えば，ヒューマンコン
ピュータインタラクション，情報通信システム，経済システム，サプライチェーンマネジメント(SCM)システ
ムなどである．多く応用されたのは自律エージェント同士の協調交渉における協調問題解決である[2,5,6,7]．
MAS にある自律エージェントは自らの目標とシステム全体の目標を持っている．自らの目標に達成するため
競争行為を行うが，システム全体の目標を達成するために協調行為を行う．現実のシステムでは，リソース
が有限であるので，エージェントは競争行為のみであれば，衝突を発生しやすく，システムがデッドロック
な状態に陥る可能性がある．システムにとって，競争と協調の両者が不可欠であり，競争行為はシステム発
展の推進力であり，協調行為はシステム求心力として捉えることができる．従って，各エージェントが競争
と協調を有効的に決める必要があるが，どのような基準で「競争」と「協調」を決めるべきか，またどの程
度までの協調をよしとするのかの決定は困難であるが重要なものとなる．
本研究では，上記の問題を解決するために，知能を有する自律的な魚が複数集合することで構成される魚
群の MAS を取り上げて協調性の分析を行った [3,6]．更に制約理論（The Thoery of Constraints, TOC）[1,4]
で提案されている Drum-Buffer-Rope (DBR) を用いて，MAS を構成する自律エージェントの競争と協調の意思
決定のための方法を提案する．これは Three-zone Buffer Management (TZBM)というアプローチである．すな
わち，DBR を制約となるエージェントの競争と協調の切り替える時期を示す方法として用い，一方 TZBM は，
エージェントが制約になる前に競争と協調を切り替えるべきタイミングと協調の程度を定める方法として提
案する．提案された TZBM アプローチを用いて，現実の大規模複雑な問題を解決することを望ましい．
2
MAS の協調性と適応性について
2.1 魚エージェントの行動ルール
魚群にはリーダが存在せず，各魚は自らの目，側線などの器官をセンサーとし，他の魚の行動と協調しな
がら自らの行動を決め，群れで行動する．このような個体から構成された魚群は MAS と言える[6]．様々な水
槽に泳げている魚群の行動から次の行動モデルを提案した[3]．
ここで，m は魚の質量であり，N は 魚群にある魚の数である．水槽に泳ぐ魚エージェントが三つの力が作
用されて，それぞれを式（3）
，
（4）および (7) に示す． F i1 は魚エージェントの推進力であり， F i2 と F i3 は
他のエージェントとのインタラクションおよび環境（水槽の壁）とのインタラクションを表す．
Tilapia
Repulsive and
heterogeneous
Bitterling
Cooperative and homogeneous
図 1 協調性の異なる MAS の行動結果
表 1 MAS のデッドロックと分割の比率
MAS
B-MAS
T-MAS
M=10
M=19
M=10
M=19
Rd
8%
50%
12%
26%
Rb
0%
0%
14%
0%
bi (r ij ) （式(5)）はインタラクションの誘引力と反発力を表し， ci (r ij ) (式(6))は成群力を示し， r ij は魚エ
ージェント間の距離である． f +wil (式(8))と f −wil (式(9))は環境からの反発力と誘引力を表す．パラメータ M は
インタラクションを行うエージェントの数であり，MAS の協調性と競争性に影響を与えるパラメータとも言
える．
モデルに含まれるパラメータは，水槽実験データを用いて求めた．実験データとして 2 種類の魚群，Bitterling
と Tilapia を取り上げ，そのパラメータを表１と表２に示す．それぞれの魚エージェントは同じルールで行動
しているが，パラメータの値が個体間の差を表す．表 1 より，各個体のパラメータの差が小さく，より均質
であることが分かる．表 2 より，各個体のパラメータの差が大きく，より異質であることが分かる．更に，
2
の大きさにより約 1.5 倍であ
Bitterling の個体間の反発力のパラメータ k1bi の大きさが誘引力のパラメータ k bi
2
ることと比べて，Tilapia の個体間の反発力のパラメータ k1bi の大きさが誘引力のパラメータ k bi
の大きさによ
り約 10 倍であることが分かった．
論文[3,6]に調べたように Bitterling は同質協調の MAS(B-MAS)と呼ばれてお
り，一方 Tilapia 魚群は異質競争の MAS (T-MAS) と呼ばれている．
2.2 魚群 MAS の協調性と競争性により適応性の変化
魚群は一つの特徴のある自律的な MAS である[6]．本研究では，魚の種類として，Bitterling と Tilapia (それ
ぞれ 20 尾の魚)の２種類の魚群を取り上げ，各群の環境変化に対する協調性と適応性について調べた．状来の
研究において，Bitterling は同質協調の MAS(B-MAS)と呼ばれており，一方 Tilapia 魚群は異質競争の MAS
(T-MAS) と呼ばれている [6]．シミュレーションの環境中に小さい箱型のトラップを設定し，その中に入った
魚群に対して突然の環境変動を与える．数学モデル上で，MAS の協調性と競争性のバランスを調整できるパ
ラメータ M の値を変化させ，各値のもとで初期場所を変化させ，それぞれ 20 回のシミュレーションを行った
結果を図１に示す．横軸は協調性の強さを，縦軸は環境変化に対する適応性の強さ(トラップから脱出した時
間の平均値)を表している．その結果から B-MAS の M＝10 の場合，元も最も強い協調性の MAS に競争性を
増加させ，協調性を減少させるのが，最も高い適応性が得られる．一方，T-MAS の M＝19 の場合，元も最も
弱い協調性の MAS に協調性を増加させ，競争性を減少させるのが，最も高い適応性が得られることがわかる．
2.3 魚群 MAS の協調性とデッドロックの関係
性質の異なる MAS に対して，デッドロックの発生率を調べた．ここで，小さいトラップに入った魚群が動
かずに，速度の大きさが０で，速度の方向が左右同じ角度で変動する状態を 300 秒以上続けるであれば，そ
の MAS は動かずに，左右方向に変動している状態から 300 秒以上続ければ，この魚群がデッドロック状態に
陥ったと定義する．B-MAS と T-MAS それぞれ協調性のパラメータを M=10(競争性を強化)と M=19（協調性
を強化）と設定し，４つの MAS モデル B-MAS (M=10)，B-MAS (M=19)，T-MAS (M=10)および T-MAS (M=19)
を作成した．それぞれに対して初期値を変更して，50 回のシミュレーションを行って，デッドロックの発生
数と群れが分割された数を集計し，比率を計算した．その結果を表 3 に示す．デッドロックの発生率の最も
高い MAS は競争性が低く協調性が一番強い MAS で，群れ分割の発生率の最も高い MAS は協調性が低く競
争性が一番強い MAS であることがわかる．従ってこのことから，MAS の協調性が強ければ強いほどすべて
いいといえず，環境変化に対して，協調性の強い同質の MAS は環境変化に対して適応し難いことが分かった．
それは突然変化した環境で一つの纏まった意見を求めるのはデッドロックを発生しやすく，様々な意見を出
して，正しい意見を求めるのは有効である．
3
Three-zone Buffer Management (TZBM)を用いた競争と協調の意思決定
3.1 Drum-Buffer-Rope (DBR) の考え方
Drum: (decide the pace of the
一番早い人 troop)
一番遅い人
制約になるエージェントを改善し，システム全体の改善
を求めるのが DBR の考え方である[1,4]．MAS を行進の隊
列に喩えるならば，図２に示すように，歩行速度の一番遅
いエージェントと一番はやいエージェントの間で，弛みを
Buffer: (a suitable distance
許容されたロープで繋げられており，各エージェントは先
Rope: (limit the leading scout
between the slowest and the
walking not so fast)
leading scouts)
行者を追い抜かさないという条件のもとで，各々の歩行速
度を調整し，隊列全体が少しでも多くの距離を進むことを
図２ Drum-Buffer-Rope の考え方
目的とする MAS と見なすことができる．早いエージェント
は幾ら早くに歩行可能であっても，ロープの長さを超えて速度をあげることはできない．即ち，ロープの弛
みの余裕があるときには，すべてのエージェントが競争することで各々の速度を上げることが許されるが，
ロープの余裕がなくなったときには，歩行速度の速いエージェントは一番遅いエージェントに対して譲歩す
ることが求められる．従って，DBR は MAS のエージェント同士の競争と協調の切り替える時期を示してい
ると考えられる．
図 3 加工・調達におけるスケジューリング問題
図３は大手製造会社のサプライチェーンマネジメント（SCM）における加工・調達に関するスケジューリ
ング問題である．図 3 の①は加工生産のスケジューリング問題であり，生産スケジューリング問題を MAS と
して多くの研究を行っている．図 3 の②は加工工場から加工された部品を梱包，パレット積載，コンテナ積
載，船積みなどの作業工程を行って，組立工場に配送するスケジューリング問題である．図 3 の③は各工程
間および工程内での作業を平準化させ，各工程の作業時間を短縮し，残業時間を削減する問題である．図 3
の①，②および③はそれぞれ特徴のある MAS と考えられる．①の生産加工問題，或いは調達問題に DBR の
応用ができる．例えば，①の生産加工の問題に対して，一番早い工程と一番遅い工程の間で，時間のロープ
で繋げられており，時間の余裕がある場合，それぞれの工程エージェントが自らのために競争行為を行うが，
時間の余裕がない場合，一番遅い工程はドラムを持って，生産のリズムを設定し，他の工程はそのリズムに
従って協調行為を行う．時間の余裕がなくなる時点で，非制約の工程はシステム全体がデッドロック状態に
陥ることを防ぐために協調行為を行う．
3.2 MAS の競争と協調の意思決定の一つのアプローチ―TZBM
制約になってから改善するのでは遅く，エージェント
Red-zone
Yellow-zone
Green-zone
にとっては，どの程度で協調するのかも重要な決定項目
agent 1, …, agent K agent 1, …, agent M agent 1, …, agent N
となる．本研究ではエージェントの競争と協調の意思決
定 を 有 効 的 に 行 う ア プ ロ ー チ Three-zone Buffer
Management (TZBM)を提案し，DBR の実現方法とも言え
る．図４に示すように，各エージェントが制約条件を満
たす状況により Green-zone, Yellow-zone および Red-zone
に Three-zone に配分する．すなわち，Green-zone に配分
(b) Red-zone が減少した
分類されたエージェントはオリジナル制約条件を満たし， (a) Red-zone が増加したた
ため、協調を許す範囲以
め、協調を許す範囲を超
Yellow-zone に分類されたエージェントは緩和された制約
内である。
えた。
条件を満たし，Red-zone に分類されたエージェントは緩
図 4 Three-zone Buffer Management
和された制約条件も満たさず，制約違反なエージェントである．従って，システムの管理者は Red-zone に注
目し，制約になった状態から脱出するために他の非制約なエージェントとの協調交渉を行う．しかし，協調
交渉しているエージェントはどの程度で協調すべきなのかを決める必要がある．
エージェント同士の協調交渉の目的は制約になったエージェントを Red-zone から脱出し，MAS 全体の目標
に達成することである．従って，協調交渉の結果は図 4 の(a)になれば，制約になったエージェントが少し改
善されでも，非制約のエージェントが Green-zone から Red-zone に変更し，協調許す程度を超えたので，こう
いった協調交渉を成立させない．一方，図４の(b)になれば，非制約のエージェントが Green-zone から
Yellow-zone に変更したが，
制約になったエージェントがRed-zone からYellow-zone に変更したことで，
Red-zone
の数を減少したので，システム全体にとって，制約違反のエージェントが減少された．従って，こういった
協調交渉を成立させる．
TZBM は DBR の競争と情報の切り替える時期を示すだけでなく，協調のエージェントにどの程度で譲歩す
るかも判断できる．従って，TZBM は MAS の競争と協調の基準を与えるアプローチであると言える．このよ
うな TZBM を，サプライチェーンマネジメント(SCM)，生産スケジューリング，在庫管理，物流管理などに
応用展開している．
3.3 競争と協調の意思決定における応用事例
上記の図３に示す②のコンテナ積載工程での作業平準化を一つの事例として TZBM アプローチの応用を説
明する．仮にコンテナ積載工程の作業開始時刻と完了時刻を既に決められている．限られた時間空間で決め
られたコンテナ積載作業を行う．一日８時間を４作業区分，２時間ずつに分ける．各作業区分での作業時間
Tdi とし，一本のコンテナの平均作業時間が一人で Ot とする．各作業区分に作業員の数を平均に配分するが，
人員の能力により作業スピードが異なる．各作業区分の作業時間をオーバーできる緩和パラメータ
(0 ≤ ε t ≤ ε t ,max ) を設定し， ε t ,max は緩和パラメータの最大値であり，オーバー時間が ε t ,max を超える場合，
εt 制約になったエージェントとし，システム全体に影響を与えるものである．各コンテナをそれぞれの作業区
分に指定され，指定状態を A，B，C および F 中の一つに属する．状態 A は指定日作業区分以外が変更不可能
であり，状態 B は指定日作業区分以外に午前或いは午後の 2 つの作業区分内で変更可能であり，状態 C は指
定日作業区分以外に前日午後と当日午前の４つの作業区分内で変更可能であり，状態 F は全工程内で自由に
変更可能とする．
各作業区分は自律なエージェントとし，配置された人員で配分されたコンテナ積載作業を最小時間で完成
されることが目標とする．ここで，人員のコストが同じと設定し，人員や作業時間を増えれば，コストを増
えると考えている．コンテナ積載工程は一つの MAS であり，それぞれの作業区分の最大完了時刻を最小にな
ることを目標とする．しかし，それぞれの作業区分に従属する人員の能力や配分された作業が異なるので，
作業の完了時刻が大きく変わる可能性がある．各エージェントの競争行為により，その作業区分の作業時間
を短縮でき，コスト削減ができると考えられる．しかし，作業量の多い作業区分は指定された作業時間作業
内で完成できないので，他の作業区分に影響を与えてしまい，結局，コンテナ積載作業全体の目的に達成で
きなくなる．従って，各エージェントの競争と協調の意思および協調の程度を有効的に決める必要がある．
N di は第 d 日の第 i 個の作業区分に配分されたコンテナの数を表し，N ok は各作業区分が指定された作業時
間 で積載できるコンテナ数である．現場で各作業区分に指定されたコンテナの配分状況は平準化の初期値と
する． Rdi = (N di − N ok ) × Ot で各作業区分のオーバーする時間を計算できる．では，Green-zone, Yellow-zone
および Red-zone を次のように定義する．
Green-zone： Rdi = (N di − N ok ) × Ot ≤ 0 であれば，該当エージェントは Green-zone に従属する；
Yellow-zone： 0 < R di = (N di − N ok ) × Ot ≤ ε t であれば，該当エージェントは Yellow-zone に従属する；
Red-zone： Rdi = (N di − N ok ) × Ot > ε t であれば，該当エージェントは Red-zone に従属する；
Red-zone にあるすべてのエージェントは現時点の作業時間をオーバーする許容範囲を超えているものであ
り，制約になる前に調整して，解消すれば，システ
ム全体の目標を達成することに影響を与えらない．
Red-zone にオーバー時間の一番大きいエージェント
は MAS の制約になるものであり，このエージェン
トを解消しなければ，コンテナ積載全工程が納期遅
れになる．TZBM アプローチを用いて，エージェン
ト間の競争と協調の切り替えるべきタイミングと協
調の程度を定める方法を提案する．競争と協調の切
り替えるべきタイミングとは，作業オーバー時間が
図 5 協調交渉における Three zone サイズの変化
ε t より小さければ，競争可能の範囲が，ε t を超えれ
ば，協調をしなければ，制約違反になる．更に，非制約なエージェントの協調を許す基準とは，制約になっ
たエージェントを改善するのは前提条件で，かつ非制約なエージェントは Red-zone に入らない限り協調すれ
ばよいと判断する．
TZBM を用いてコンテナ積載作業平準化の問題へ応用した結果を図５に示す．図５により，最初の
n ∈ [0,150] あたり，制約緩和なしで，エージェント間の協調交渉により，制約になったエージェントの数が
24 個から１３個に減少した．即ち作業コストが増えない条件で，エージェント間の協調交渉だけで改善され
た．制約になったエージェントが競争行為で，非制約のエージェントは協調行為を行った．非制約のエージ
ェントは現時点での制約条件を違反しないことが協調を許す基準である．ただし，n > 150 とき，制約違反の
エージェントが存在する状況で，協調を許す基準を守る協調交渉が成立できなくなった．このとき，制約条
件を緩和させて，即ち作業時間を少しオーバーして，再び協調交渉を行う．これにより緩和されたエージェ
ントを増加しながら，制約違反のエージェントがすべて改善された．即ち，作業を要求された時間に完成で
きるエージェントは作業時間の最大値をオーバーしたエージェントに許す範囲で譲歩させて，MAS の全体の
最大の作業時間を短縮され，残業時間を減少させた．最終的にコンテナ積載工程のコストを削減できたと言
える．
参考文献
[1] E.M. Goldratt and J. Cox: The Goal : A Process of Ongoing Improvement, The North River Press, (1984).
[21] S. Kraus, ``Negotiation and cooperation in multi-agent environments'', Artificial Intelligence, Vol.94, pp. 79-97
(1997).
[3] Y. Tian, N. Sannomiya and H. Nakamine: A simulation Study on Adaptability to Environmental Variations Based on
Ecological Systems, Proc. of the American Control Conference, pp. 1759-1763 (1999).
[4] D. Smith: The Measurement Nightmare: How the Theory of Constraints Can Resolve Conflicting Strategies, Policies
and Measures, St. Lucie Press, (2000).
[5] Y. Tian, Y. Liu, K. Shimohara and T. Sawaragi: Cooperative Negotiation Strategy in Multi-agent System, Proc. of the
42nd IEEE Conference on Decision and Control, pp. 2549-2554, (2003).
[6] Y. Tian, N. Sannomiya, H. Inoue and K. Shimohara: Cooperation of Multi-agent System and Its Composition, Proc. of
the 16th World Congress of International Federation of Automatic Control, pp. Mo-M21-TO_3 pp.1-6 (2005).
[7] Y, Liu, Y. Tian and T. Sawaragi: Applying DBR and MAS to Solving Container Loading Problem, Journal of
Japan Society for Fuzzy Theory and Intelligent Informatics, Vol.18 No.6, pp.837-848 (2006).
機械加工における長いスパンにおけるプロセス制御に関する研究
工学研究科マイクロエンジニアリング専攻 茨木創一
Abstract: Numerous researches have been reported in the literature on the cutting force control in end milling
processes. There have been, however, very few practical applications actually employed in commercial products.
In this research, as a practically feasible and effective cutting force control scheme, we propose a long-term cutting
force control strategy with a particular focus on the regulation of the tool life. First, the cutting force is monitored
only at every “check point” set on the tool path. During actual cutting, the feedrate profile is updated at every check
point, such that the cutting force is regulated along the target profile. The target profile is given such that the
desired cutting distance can be machined before reaching the end of tool life. The effectiveness of the present
approach is experimentally investigated by cutting experiments of hardened steel.
Key words: end milling, cutting force control, tool life, intelligent machining.
1. はじめに
金型加工や部品加工などの機械加工において，必要な加工精度を確保し，かつ可能な限り高能率な加
工を実現するためには，工程設計や加工条件の設定，加工状態のモニタリングに基づく加工条件の適応
的変更など，熟練技術者の経験と技術が必要不可欠であるとされている部分は多い．特に近年，工作機
械の高速・高加速度化や，工具のコーティング技術の発達，難削材の加工需要の増大など，機械加工に
関する技術の進歩は早く，それらを最大限に活用して高い能率で加工を行うことは，熟練技術者にとっ
ても難しいケースがしばしば見られる．工作機械や工具，ワークに対する十分な経験がない場合には，
安全性のために必要以上に保守的な加工条件で加工を行わざるを得ず，加工能率が犠牲になる．
工作機械や工具，ワークに対する十分な経験がない場合にも，工作機械や工具の性能を最大限活用し
安全かつ高能率な加工を行うためには，加工条件の設定や加工状態のモニタリングを自律的に行うこと
ができ，オペレータによる意思決定を支援するシステムが必要不可欠であることは，以前から多く指摘
されてきた．この観点から，機械加工プロセスの知的制御に関する研究は 1970 年代頃から非常に多く報
告されている．それらのなかで主なアプローチのひとつは，加工中に工具に作用する切削抵抗を，加工
状態の安全性を表す重要なパラメータと捉え，その制御という形で加工条件の自律的最適化を図るアプ
ローチである．しかし，多くの研究努力にもかかわらず，現在までのところ，そのような切削抵抗の制
御システムが実用化された例は極めて少ない．
その主な原因としては，① 切削抵抗をモニタリングするために，信頼性の高いセンサを導入しなくは
ならない，
② 切削抵抗の制御自体が達成できたとしても，
コストに見合う実際的な成果が期待できない，
ことが挙げられる．結果的に，これまでに提案されたアプローチの多くは，適応制御理論の機械加工プ
ロセスへの実験室レベルでの適用試験，という域を出ることができなかった．
本研究では，工作機械の高速・高加速度化，工具のコーティング技術の発達，難削材の加工需要の増
大といった背景を踏まえ，現実的な課題に対応した，非熟練者であっても目標の加工能率を達成できる
ための加工条件の意思決定支援システムを構築することを目的とする．特に本研究で着目した点は，工
具の摩耗の進行など，長いスパンにおける加工プロセスの変化を観察し，それに応じて加工条件を最適
設定するためには，熟練オペレータの経験と技術が必要不可欠であるとされている点である．本研究で
は，工具摩耗の進行などより長いスパンにおける加工プロセスの変化を陽にモニタリングし，その制御
という形で，オペレータの意思決定を支援するシステムを，実用的な形で構築することを目的とする．
2. 制御法の概要
Feedrate scheduling
based on cutting force
prediction model
Feedrate scheduling
based on cutting force
prediction model
エンドミル加工における切削抵抗制御に関
Machining
Machining
する過去の研究においては，多くの場合，操
作しやすい加工条件として，送り速度が制御
変数として用いられている．切削抵抗が目標
Monitoring of cutting force
at every check point
値より大きい領域では，工具に過多の負荷が
作用し工具損傷が生じるのを避けるため，送
り速度を低減する．切削抵抗が目標値より小
Estimation of tool wear
prgoress
さい領域では，加工能率を上げるために送り
速度を増大させる．過去の研究の殆どにおい
Update of feedrate profile
て，切削抵抗制御の基本的な目的は以上の通
りである．
過去の研究において，切削抵抗の制御手法
(a) Conventional feedrate
(b) Proposed approach
は以下の二つに大別できる．
scheduling approach
1つ目は切削抵抗をリアルタイムで検出し，
図 1 提案する制御アプローチの概念図
フィードバック制御で最適な加工条件に制御
する方法である．フィードバック制御手法を
適用する際，実用的にもっとも大きな課題は
切削抵抗のモニタリング手法の信頼性である
と考えられる．圧電素子を用いた動力計など，
(a) On a straight path
(b) On a circular path
加工中の切削抵抗を十分な精度で測定できる
センサは普及しているが，コスト面など実用
的な適用には課題が多い．外部センサを用い
ず，サーボモータや主軸モータの電機子電流
から，切削抵抗を推定する手法はコストがか
からないが，それに基づいて常にフィードバ
図 2 エンドミル切削における工具－ワーク干渉
ック制御を行うためには信頼性が十分でない
のが一般的である．
2つ目は加工プロセスをモデル化し，モデルに基づいて切削抵抗の予測値を目標値に制御できるように
送り速度を予め調整したNCプログラムを作成するという，フィードフォワード制御的な手法である．こ
の方法は比較的シンプルであるが，実際の加工プロセスのモニタリングを行わない．このため，加工プ
ロセスの変化に対応することが不可能である．
加工の進行に伴って加工プロセスの変化をもたらす代表的な要因として，工具摩耗の進行が挙げられ
る．一般に，工具磨耗の進行に伴い，切削抵抗は上昇する．切削抵抗制御に関する過去の研究のなかで
は，工具摩耗の影響など，長期的な加工プロセスの変化を陽に考慮した制御手法を提案したものは殆ど
ない．本研究では，制御対象を主に長期的な加工プロセスの変化とし，工具磨耗の進行を自律的に制御
し，求められる加工能率を安全に達成することを支援する制御システムを構築することを目的とする．
具体的には，工具経路に切削抵抗のモニタリングをするための「チェックポイント」を設定し，チェ
ックポイントでのみ限定的に切削抵抗のモニタリングを行い，それから工具摩耗の進行度の推定を行い，
切削距離の目標値が得られるよう送り速度を調整する切削抵抗の制御手法を提案する．具体的には，チ
ェックポイント間で，常時モニタリングを必要としない加工プロセスのモデル化によるフィードフォワ
ード制御を基本とし，加工プロセスの変化に対応できるよう，チェックポイントでモニタリングされた
切削抵抗からモデルの更新によるフィードバック制御を組合せ，より長いスパンに対応する切削抵抗の
制御を提案する．図１は，提案する手法の概念図を，従来のフィードフォワード制御手法との比較で示
したものである．
600
Cutting force (normal direction) N
600
C1
550
550
C2
500
500
C3
450
400
C4
400
350
350
C2
300
250
200
200
0
200
400
600
800
C3
300
250
150
C1
Tool wear index N
S1
C4
450
1000
150
S1
0
Cutting distance m
図 3 工具寿命にいたるまでの切削抵抗の測定値
200
400
600
800
1000
Cutting distance m
図 4 図 3 の切削抵抗から換算した工具磨耗指標
3. 工具磨耗の進行度の推定とその制御
3.1 切削抵抗の予測モデル 前節で述べた通り，本研究では，加工条件から切削抵抗を予測するための
モデルを用いるフィードフォワード制御手法を基本としている．本節では，本研究で用いる切削抵抗予
測モデルについて述べる．
本研究では垣野ら[1]が提案した切削抵抗予測モデルを用いる．このモデルは，変形前の最大切りくず
厚さtm と切削円弧長L をパラメータとして，以下の2次の応答曲面モデルを用いて切削抵抗を推定する．
パラメータtm 及びLの定義については図2を参照のこと．
Y = β 0 + β 1 X 1 + β 2 X 2 + β 11 X 12 + β 22 X 22 + β 12 X 1 X 2
(1)
ただし，X1 と X2 はそれぞれ tm と L を正規化したものであり，Y は切削抵抗の予測値を示す．未知の 6
個の係数βi は回帰分析により実験により同定する．
3.2 工具磨耗の進行度の推定 一般に，切削抵抗と工具摩耗の進行度には相関があることはよく知られ
ている．ここでは，加工条件や加工形状が変化している条件下でも，加工中に測定した切削抵抗から，
工具摩耗の進行度を推定するための方法を提案する．
チェックポイント毎に切削抵抗のモニタリングを行い，その切削抵抗を用いて「現時点で，仮に標準
条件によって加工した場合の切削抵抗値」を推定する．ここで，標準条件とは，3.1 節に示した切削抵抗
予測モデルにおける tm と L の中央値に対応する直線加工の加工条件と定義する．この場合，標準条件に
おいては X1=X2=0 であるため，この値は切削抵抗予測モデル式（1）の定数項βi に等しいと考えられる．
すなわち，チェックポイントごとに測定された切削抵抗からβi を更新し，それを持って工具摩耗の進行
度を推定するための指標とする．今後，このように同定されたβi を工具磨耗指標と呼ぶ．
加工条件の違いによらない工具摩耗の進行度の推定方法の有効性を確かめるために，基礎的な切削実
験を行った．ワーク材質は高硬度鋼SKD61（HRC53）で，工具は（Al，Ti）Nコーティッド超硬エンド
ミル（φ6mm，6枚刃）を用いた．異なる5つの条件（S1, C1～C4）のもと，直線切削を工具寿命に達する
まで続ける．チェックポイント毎に切削抵抗を工具動力計で測定した結果を図3に示す．
当然のことながら，各条件によって加工中の切削抵抗に大きな差があった．また，工具寿命に達した
ときの切削抵抗も各条件によって違いがあった．このため，単純に切削抵抗を見るだけでは，工具寿命
を判断することは困難だとわかる．
次に，上記と同じ結果を用いて，工具摩耗指標βi を計算した結果を図4に示す．いずれの条件でも，
工具摩耗指標の増加は，比較的近くなり，いずれの場合も350～450N程度に達したときに工具寿命とな
っていることがわかる．これより，工具寿命にいたる切削抵抗のレベルは加工条件によって大きく異な
るが，工具摩耗指標を導入することによって，
加工条件が異なっても共通して工具摩耗の進
行度の判定をより容易に行えるといえる．
420
400
cutting force N
380
predicted tool life
under constant feedrate
β0 under constant
o
desired tool life
o
feedrate
360
3.3 工具寿命の制御 切削抵抗を常時モニ
target curve
340
タリングし，
それに基づく制御を行うことは，
320
β0 under feedrate regulation
特に低コストな（実用的な）モニタリング手
300
法を適用する場合，信頼性に問題がある．し
280
かし，工具経路上で，定常的な加工が続き，
260
比較的切削抵抗のモニタリングがしやすい区
間に「チェックポイント」を設置し，そこで
240
0
200
400
600
800
1000
cutting length m
のみ切削抵抗のモニタリングを行うことは，
はるかに容易である．本研究では，チェック
図 5 工具寿命にいたるまでの切削抵抗の測定値
ポイント毎に切削抵抗のモニタリングを行い，
まず工具摩耗の進行度を推定する．工具摩耗
の進行度から，このままの条件で加工を続けた場合の工具寿命を推定し，さらに，チェックポイント毎
に送り速度のオーバライドを変更することで，工具寿命の目標値を得る制御手法を提案する．
構築した手法を用いて，与えられた切削抵抗上昇線に沿って切削抵抗を制御できるか確認することを
目的として，基礎的な切削実験を行った．
実験に使用したワーク，工具，工作機械は 3.2 節の実験で使用したものと同一である．最初に，一定
の送り速度のまま，工具寿命に達するまで直線切削を続けたところ，図 5 に示すように工具寿命指標は
上昇し，切削距離約 650m で工具寿命にいたった．次に，同様の実験を提案した制御手法を適用して行
う．切削距離 30m ごとに設けられたチェックポイントにおいて，切削抵抗を動力計を用いて測定し，工
具寿命指標を計算，それが目標軌跡に沿って上昇するように送り速度の制御の更新を行う．ここで，目
標軌跡は工具寿命が約 1000m まで延長することを目的に設定した．
図 5 に示されているように，送り速度のチェックポイント毎の更新により，工具寿命指標は目標軌跡
に沿って良好に制御することができた．その結果，工具寿命は約 900m となり，目標とした工具寿命が
ほぼ達成でき，制御を行わない場合と比較して約 250m 延長された．
4. おわりに
機械加工におけるオペレータの意思決定を支援するシステムとして，本研究では，加工プロセスの監視と送
り速度の制御による，長いスパンのおける加工プロセスの変化を制御する手法に取り組んだ．これまでに多く
研究報告されている切削抵抗制御の手法と比較して，本研究で提案した手法の特徴は以下の通りにまとめら
れる．
(1) 切削抵抗を常時モニタリングすることは要求せず、「チェックポイント」毎にのみ行いため、より｢安価な｣モ
ニタリング手法が採用できると考えられる。
(2) 工具磨耗の進行に伴う切削抵抗の上昇を制御することに着目し，所定の工具寿命を安全に得ることがで
きることを目的とした制御手法を提案した．
提案した手法の実際的な加工システムへの適用，より複雑で一般的な工具パスへの適用について，今後引
き続き検討を行う．
参考文献
[1] Kakino, Y., Ohtsuka, H., Nakagawa, H., and Hirogaki, T., Sasaki, M., 2000, A Study on Endmilling of Hardened
Steel (1st Report) - Simplified Prediction Model for Cutting Forces and Control for Constant Cutting Forces Using
this Model -, Journal of the Japan Society for Precision Engineering, Vol.66, No.5, pp.730-734 (in Japanese)