計測自動制御学会東北支部 340 回研究集会 (2022.12.23)
資料番号 340-5

協調牽引により溶接ケーブルを最適形状に維持する群ロボットシステム-第 6 報: 多体モデルの関節角記号微分に基づく
弾性ケーブルの 3 次元形状予測の高速化Cooperative Towing by Multi-Robot System That Maintains Welding Cable
in Optimized Shape-6th report: fast 3D shape estimation of elastic cable based on joint angle
symbolic diﬀerentiation of multibody model○小熊一矢*，岡田佳都*, 衞藤晴彦**, 坂井郁也**, 大野和則*, 多田隈建二郎*, 田所諭*
○ Kazuya Oguma*, Yoshito Okada*, Haruhiko Eto**, Yuya Sakai**, Kazunori Ohno*,Kenjiro Tadakuma*, Satoshi Tadokoro**東北大学, **住友重機械工業株式会社*Tohoku University, **Sumitomo Heavy Industries, Ltd.
キーワード : 三次元ケーブル形状 (3D cable shape), 最適化問題 (Optimization), 溶接ロボット (Welding robot)
連絡先 : 〒 980-0845 宮城県仙台市青葉区荒巻字青葉 468-1 東北大学 情報科学研究科応用情報科学専攻 田所研究室
小熊一矢 E-mail: oguma.kazuya@rm.is.tohoku.ac.jp

1.

緒言

船舶等の大型構造物の溶接においてロボットが導入され省人化と溶接時間短縮に貢献している 1) ．しかし，現状のロボットの多くは自走しない大型マニピュレータ 2) か自律走行機能を持たない小型台車 3) であるため．狭所や高所における長尺溶接に適さない．
そこで著者らは、複数台の車両による溶接トーチとケーブルの協調運搬により狭所・高所での長尺溶接に適用可能な群ロボットシステムを研究している (図1)．第 1 報から第 4 報 4) および 5) では，群ロボット

Fig. 1 Cable Towing with Multiple Robots

実機による水平面上のケーブル牽引を実現した.
第 5 報 6) では，水平面ではなく溶接現場で多く見する．従って，変形したケーブルが環境と接触しないられる斜面や曲面におけるケーブル牽引を想定し，車ケーブルの形状の予測は牽引車両のオンライン誘導

両により両端が固定把持されたケーブルの三次元形状アルゴリズムの一部となるため，十分に高速である必と拘束力を予測する手法を提案した．牽引中，ケーブ要がある．一方で前稿の形状予測には 300-400ms の時間を要しており 6) ，最大 10Hz での実行を想定する

ルは重力の影響を受けながら車両間で三次元的に変形よう，ケーブル形状の事前予測が必要であった．

–1–

誘導アルゴリズムに組み込むには不十分であった．
そこで本稿では，ケーブルの三次元形状予測の20ms 以下への高速化を目標とする．前稿の予測手法
では，形状予測を最適化問題に帰着させた上で，最適化中の勾配算出に数値微分を用いた．数値微分は計算に時間を要するため，予測にかかる時間を増大させる主要因となっていた．本稿では，ケーブルの幾何に基Fig. 2 The Calculation Model of 3D Cable Shape

づく記号微分への置換による高速化手法を提案する．
また計算時間を実測し，目標とした 20ms 以下への高速化を確認した．

aJi ,1 ，弾性係数を k とする．各リンクの質量を m とする．
また，リンク Li と LN −1 を端とする閉区間をセグ

2.

関連研究

メントと呼び，Si で表す．セグメント S0 はモデル全索状弾性体の力学的モデリングや形状推定に関する既存手法には，梁の大変形理論解析に基づく方法 7)や多体モデルに基づく方法 8) 9) がある．
7) は連続梁の大変形理論に基づいて両端固定時の弾性ベルトの二次元形状を導いている．8) は剛体リンクを回転関節で直列に接続した多体モデルとして板バネをモデル化した上で，重力は考慮せず形状推定を行っている．9) の対象は，片端が地面に固定され，自由端に加わる推力で浮上する索状ロボットである．多体に相当する．関節 Ji の回転に伴い，セグメント Siが Ji を中心に回転すると考える．
各関節，リンクおよびセグメントに座標系を設定す
る．各座標系の原点は，関節の場合は回転中心，リンクおよびセグメントの場合は重心とする．世界座標系は W とする．座標系 A から B への変換を，並進ベクトル B pA および回転クォータニオン B q A で表す．座標系 A から見た B の状態ベクトル x を A xBと書く．

体モデルを用いて浮上の安定性を検討している．
いずれも，重力が影響する三次元空間におけるケー3.2

ブルを考える本研究とは目的や対象が異なる．

ケーブル平衡形状算出のための最適化問題

溶接トーチの溶接中の移動速度は 0.01[m/s] と低速

3.

車両間ケーブル形状の予測手法の提案であるため，トーチから延びるケーブルも低速で牽

本稿ではケーブルを多体モデルで表現した上で，最小ポテンシャルエネルギーの原理に基づく形状予測手引される．本稿では，ケーブル牽引は準静的であり，ケーブルは車両間に張られた状態で釣り合い形状をとると考える．そこで，前節で説明したケーブルの力学法を提案する．

モデルの平衡状態を考えることにより，ケーブルの形3.1状を算出する．

車両間ケーブルの力学モデル

モデルが安定な力学的平衡状態にあるとき，その2 台の車両の間に張ったケーブルを，本稿では図 2

ポテンシャルエネルギーは極小値となる．言い換
えれば，ケーブルのポテンシャルエネルギーを目のようにモデル化して考える．
モデルは，N 個の剛体リンク {L0 , ..., LN −1 } を的関数 F (θ J )，ケーブル両端の位置姿勢に関するN − 1 個の関節 {J1 , ..., JN −1 } で直列に接続して構制約条件を G(θ J ) = 0，ケーブル内の関節角 θ J =成する．各関節は曲げ方向の 2 自由度に回転でき(θJ1 ,0 , θJ1 ,1 , ..., θJN −1 ,0 , θJN −1 ,1 ) を最適化パラメータる．関節 Ji の各自由度について，関節角を θJi ,0 おとする最小化問題の解が，ケーブルの釣り合い形状とよび θJi ,1 ，回転軸を表す単位ベクトルを aJi ,0 およびなる．
このポテンシャルエネルギー最小化問題を以下の通–2–り定式化した．式 (1) が目的関数 F (θ J ) であり，初で書ける．

項が各リンクの重心位置に依存する位置エネルギーθ̇Ji ,k A aJi ,k ×A

∂ p Ljを，次項が各関節の関節角に依存する弾性エネルギー∂θJi ,kを表す．ただし，g は重力ベクトルである．式 (3) が== aJi ,k ×A

制約関数 G(θ J ) であり，ケーブル両端のリンクの相((Ap L j − A p Ji

θ̇Ji ,kAp L j − A p Ji))
(4)対位置姿勢と添え字 r で示した目標値との誤差を表

また，関節角 θJi ,k に対するリンクの姿勢 A q Lj のす．ただし，Vec(q) はクォータニオン q の虚部を 3偏微分は，その定義より A q Lj の時間微分と Ji の回次元ベクトルとして取り出す関数である．

転角速度の比となる．そのため，実部が w で虚部がv であるクォータニオンを (w, v) で表すと，式 (5) で最小化

F (θ J ) =N−1∑(−mW pLi (θ J ) · W g

i=0+N−1 (∑i=1

制約[G(θ J ) =書ける．

)

1k(θJ2i ,0 + θJ2i ,1 )2(1)∂ A q Lj)pLN −1 (θ J ) − L0 pLN −1 ,r )(Vec L0 q LN −1 (θ J )L0 q −1
LN −1 ,rL0∂θJi ,k(2)

1= 2(0, θ̇Ji ,k A aJi ,k )A q Ljθ̇Ji ,k1 A= (0, aJi ,k )A q Lj
2

式 (4) および (5) は，リンクのみならずセグメント]の位置姿勢に対しても成り立つ．

=0(3)

3.3.2

上述の通り定式化した最小化問題は，等式制約付き非線形最適化問題の標準的な形式をとっているため，一般の最適化ソルバを用いて解ける．

目的関数の記号微分

目的関数 (1) は，ケーブル全体に相当する質量 N mのセグメント S0 を用いて式 (6) のように書ける．

F (θ J ) = −N m pS0 (θ J )· g+W

3.3(5)WN−1 (∑i=1

最適化高速化のための記号微分の導出

1k(θJ2i ,0 + θJ2i ,1 )2(6)前節で定式化したポテンシャルエネルギー最小化問

式 (6) を関節角 θJi ,j で偏微分し，目的関数の記号題を解くにあたっては，制約条件 (3) の下で目的関数微分 (7) を得る．初項の位置エネルギーの微分にあ(1) の値を小さくする方向に最適化パラメータ θ J をたっては，関節の回転に伴い質量 (N − i)m のセグメ更新していく．θ J の更新方向の決定には，目的関数

ント Si のみが関節を中心に回転することに注意したおよび制約関数の勾配が必要となる．前稿では勾配の上で，式 (4) を用いた．

算出に数値微分を用いていたが，本稿では高速化のた(∂F (θ J )= − (N − i)m W aJi ,j (θ J )θJi ,j)(× W pSi (θ J ) − W pJi (θ J ) · W g + kθJi ,jめに記号微分を用いることとした．本節では目的関数および制約関数の記号微分を導出した後，計算量オー(7)ダーの考え方に基づき高速化効果を検討する．

3.3.1

3.3.3

数学的準備

制約関数 (3) を関節角 θJi ,j で偏微分し，制約関数

関節 Ji の回転に伴い，リンク Lj が Ji を中心に回転する状況を考える．このとき，関節角 θJi ,k に対するリンクの位置 A pLj の偏微分は，その定義より Lj

制約関数の記号微分の記号微分 (8) を得る．上の位置成分の微分にあたっては式 (4) を，下の姿勢成分の微分にあたっては式の並進速度と Ji の回転角速度の比となるため，式 (4)–3–)

(5) を用いた．
Table 1 PC speciﬁcation
CPUIntel Core i9-9880H

L0aJi ,j (θ J )) 
 (LL0 × 0 pL
(θ)
−p(θ)
JJJiN −1∂G(θ J ) =θJi ,j(1L0Vec 2 (0, aJi ,j (θ J ))
L0q LN −1 (θ J )L0 q −1
LN −1 ,r(8)3.3.4

Memory32GB

Table 2 Experimental Parameters
Parameters
Cable End Pose
Items

Value

x[m]

y[m]

z[m]

Cable length [m]Cable Radius [m]

2.0
0.02

1.0
1.4

0.0
0.5

0.0
0.5

Elastic Modulus [Nm/rad]

0.75

1.8

Cable Weight [kg]1.0Model Resolution [m/joint]

0.1

計算量オーダーに基づく高速化効果の検討

前稿の最適化計算に用いた数値微分と本稿で用いる記号微分の計算量オーダーの比較を通じて，記号微分の導入による計算高速化の効果を検討する．
まず数値微分の計算量オーダーを考える．数値微算できる．事前計算は全ての導関数を合わせて O(N )分においては，変数の値をわずかに変化させ，変化さであり，本計算も変数 θ J 全てに対して繰り返すとやせる前後の関数の値の差分をとることで，変数に対すはり O(N ) となる．従って，事前計算と本計算による関数の微分値を得る．従って，1 つの変数に対するり全変数に対する微分値を得るための計算量オーダー

数値微分の計算量オーダーは，関数のそれと等しい．

は，結局 O(N ) となる．

変数が多次元である場合は，それぞれの変数に対し同様の計算を行い，全変数に対する微分値を得る．

制約関数の導関数 (8) も，式中の aJi ,j ,pLN −1 ,pJi ,q LN −1の事前計算と，O(1) の本計算に分割できる．従って，目的関数 (1) は N 個のリンクの位置エネルギーとN − 1 個の関節の弾性エネルギーの総和であるため，全変数に対する微分値を得るための計算量オーダーは，目的関数の場合と同じく O(N ) となる．

その計算量オーダーはランダウの O 記法を用いて以上の検討より，数値微分から記号微分への置き換

O(N ) となる．従って，1 つの変数 θJi ,j に対する目えにより，最適化計算内の勾配算出の計算量オーダー的関数の数値微分の計算量オーダーも O(N ) である．

が O(N 2 ) から O(N ) に改善する．また，記号微分にN に比例する次元数の変数 θ J 全てに対する微分値をよる高速化の効果は，モデルのリンク数 N が大きい2得るための計算量オーダーは，結局 O(N ) となる．

ほど指数的に大きくなる．

制約関数 (3) はケーブル両端のリンクの相対位置姿勢の誤差であり，N 個のリンクの位置姿勢を順に計4.

評価4.1

実験手法算する必要があるため，その計算量オーダーは目的関数と同じく O(N ) となる．従って，変数 θ J 全てに対する微分値を得るための計算量オーダーも，目的関数2の場合と同じく O(N ) となる．

両端固定ケーブルの形状を予測する最適化問題の計算時間を，前稿の数値微分と本稿の記号微分を用いた

次に記号微分の計算量オーダーを考える．記号微場合で比較し，高速化の効果を定量的に検証した．

分においては，関数を解析的に微分して得た導関数に実験に用いた PC の仕様を表 1 に示す．実験パラ変数を代入し，微分値を計算する．従って，1 つの変メータを表 2 に示す．ケーブル両端の相対位置を変数に対する記号微分の計算量オーダーは，その変数にえ，最適化問題の求解にかかる時間を評価した．相対対する導関数のそれと等しい．変数が多次元である場

位置は x が 3 種類，y が 2 種類，z が 2 種類の 12 条

合は，それぞれの変数に対し同様の計算を行い，全変件とした．各条件に対し求解を 1000 回繰り返し，計数に対する微分値を得る．

算時間の平均と標準偏差を計測した．

目的関数の導関数 (7) は，式中の aJi ,j ,pSi ,pJi を事前計算することで，O(1) の計算量オーダーで本計–4–

Table 3 Experimental Result (Calculation Time)
Cable End Pose
Calculation Time
New[ms]Previous6) [ms]x[m] y[m] z[m]Ave.
SD Ave.
SD
1.0

0.0

0.0

7.70

2.34

388

9.18

1.0

0.0

0.5

11.78

3.58

346

16.20

1.0

0.5

0.0

10.82

3.50

345

9.29

1.0

0.5

0.5

6.05

2.69

299

5.32

1.4

0.0

0.0

4.71

1.98

323

6.48

1.40.00.55.35

2.36319

5.31

1.4

0.5

0.0

5.94

2.37295

4.94

1.4

0.5

0.5

7.04

2.81446

6.05

1.8

0.0

0.0

4.90

2.07313

4.45

1.8
1.8

0.0
0.5

0.5
0.0

6.96
7.25

2.83
2.67

352
411

8.88
7.67

1.8

0.5

0.5

12.25

3.03

404

11.26

4.2

実験結果

参考文献

ケーブルの形状推定に要した時間を表 3 に示す．
前稿と比較すると，平均の計算時間が約 50 分の 1 に高速化され，4-12[ms] 程度での予測が可能となった．
また標準偏差は 2-3[ms] 程度となり，計算時間のばらつきが少なく高い確率で平均前後の計算速度が見込めることがわかった．これらの結果から，記号微分の導入により両端固定ケーブルの形状予測が高速化され，目標とした 20[ms] の計算時間を達成できたことがわかった．

5.

結言

本稿では仮定したケーブルの両端の姿勢に基づくケーブルの三次元形状の予測を，記号微分の採用により高速化する手法を提案した．提案手法の計算時間を実測し，牽引車両のオンライン誘導アルゴリズムに含めることが可能となる 20[ms] 以内で形状予測が可能と確認した．
次報以降では，ケーブル形状と拘束力の予測に基づいて，斜面や壁面上での群ロボットによる最適なケーブル牽引動作を生成する手法を検討する．

–5–1) Donghun Lee: “Robotics in the shipbuilding industry”, Robotics and Computer-Integrated Manufacturing, Volume 30, Issue 5, pp.442-450 (2014)
2) Donghun Lee, Sungcheul Lee, Namkuk Ku,Chaemook Lim, KyuYeul Lee, Taewan Kim, Jongwon Kim: “Development and application of anovel rail runner mechanism for double hull structures of ships”, Proceedings of IEEE international
conference on robotics and automation(ICRA),
pp.3985-3991 (2008)
3) Xueqin Lü, Ke Zhang, Yixiong Wu: “The scam
position detection and tracking for the mobilewelding robot”, The International Journal ofAdvanced Manufacturing Technology 88.5-8, pp.
2201-2210 (2017)4) 小熊一矢, 岡田佳都, 衞藤晴彦, 坂井郁也, 大野和則,
多田隈建二郎, 田所諭: “協調牽引により溶接ケーブルを最適形状に維持する群ロボットシステム (第 4 報:
車間距離の維持と目標車両位置の連続性の考慮による狭所における無衝突牽引動作の実現)”, 第 23 回 計測
自動制御学会システムインテグレーション部門講演会(SI2022), (2022)
5) Ryota Suzuki, Yoshito Okada, Yoshiki Yokota,
Tatsuyoshi Saijo, Haruhiko Eto, Yuya Sakai,Kenichi Murano, Kazunori Ohno, Kenjiro
Tadakuma, Satoshi Tadokoro:“CooperativeTowing by Multi-Robot System That Maintains Welding Cable in Optimized Shape”,
IEEE Robotics and Automation Letters doi:
10.1109/LRA.2022.3183529., (2022)
6) 小熊一矢, 岡田佳都, 衞藤晴彦, 坂井郁也, 大野和則,
多田隈建二郎, 田所諭: “協調牽引により溶接ケーブルを最適形状に維持する群ロボットシステム (第 5 報：
最小ポテンシャルエネルギーの原理に基づく車両間ケーブルの三次元形状と両端拘束力の実時間予測)”,
第 23 回 計測自動制御学会システムインテグレーショ
ン部門講演会 (SI2022), (2022)
7) 小澤悠，渡辺将広，多田隈建二郎，田所諭: “弾性履
帯の適応変形により高走破性を実現する単輪クローラ機構-非線形梁に基づく大変形解析-”，第 40 回日本ロボット学会学術講演会 (RSJ), (2022)
8) Hiromi Mochiyama: “Model Validation of Discretized Spatial Closed Elastica” IEEE/RSJ International Conference on Intelligent Robots andSystems (IROS), (2016)9) 安部裕一, 山内悠, 石井昭裕, 昆陽雅司, 多田隈建二郎,田所諭: “2 軸能動空気噴射モデルを持つ索状体の浮上制御方法の提案”, 第 18 回システムインテグレーショ
ン部門講演会 (SI2017), (2017)