アジマススラスタ首振角に制約を課す場合の動的最適推力配分に関して
正 員
正 員

梶 原 宏 之∗中 村 昌 彦∗∗∗学生員

上 田 泰 秀 ∗∗正 員

木 村元∗

On dynamically optimal thrust allocation under the azimuth thrusters’ angle constraints
byHiroyuki Kajiwara, Member
Yasuhide Ueda, Student Member
Masahiko Nakamura, Member Hajime Kimura, Member

Key Words: DPS (Dynamic Positioning System), Thrust Allocation，Closed-loop Simulation
2. 推力配分法

1. 緒 言本稿は，アジマススラスタを装備した浮体の推力配分問題について考察する。いま，浮体の運動としては surge，sway，yawのみを考え，望ましい運動を得るために，ある時刻 t において，surge 方向の力 X(t)，sway 方向の力 Y (t)，yaw 方向のモーメント N (t) が必要であるとする。これらをまとめて制御力と呼ぶことにする。この制御力は，各アジマススラスタの推力 Ti (t)（i = 1, · · · , n）の合成和として得られることから，実際には各アジマススラスタへの推力配分問題を解く必要がある。これは，数学的には，未知数の数が方程式の数より大きい場合の連立一次方程式の求解問題であるため，解は一意には定まらない。すなわ
ち，ある制御力をもたらす様々な推力配分の仕方がある。一番よく知られた方法は，推力の大きさの２乗和を最小とするもので，ラグランジュ法と呼ばれている。この従来法は実装が簡単である
が，著者らは，平水中の模型実験において，スラスタが頻繁に首振りを行うことを観察していた。
Fossen1) のグループは，アジマススラスタ首振角などの制約を考え，制約条件付きの線形計画問題または２次計画問題として定式化しているが，実時間計算できる実装法に至っていないようである。これに対して，大坪ら 2) は，従来法による配分結果に，当該連立一次方程式から得られる推力配分行列の零化空間基底の適当な線形結合を加えることで，アジマススラスタ首振角などの制約を満たすことができることを示した。また五百木ら 3) は，この方法を実時間で厳密に実装できる方法を提案した。しかしながら，この提案法を用いた平水中の模型実験においては，スラスタの首振角は所定の範囲内に収まるが，それでもなおその範囲内で頻繁に首振りを行うことが観察されていた。これは，各時刻では最適配分が行われているが，時間区間でみると首振り量の総和が抑制されていないためと考えられ，スラスタの動特性を考慮したいわゆる動的推力配分問題の定式化が必要であるとされていた。
以上の背景のもと，本稿は，静的推力配分問題を解いて得られた従来法と提案法を，閉ループ系シミュレーションを通して評価し，スラスタの振れ回り現象を再現し，その原因について考察し，対策について検討した結果をまとめることを目的とする。
∗九州大学大学院工学研究院海洋システム工学部門∗∗

九州大学大学院工学府海洋システム工学専攻∗∗∗

九州大学応用力学研究所

2.1

推力配分方程式

いま，n 個のアジマススラスタを装備した浮体を考え，i 番目の取り付け位置を (xi , yi )，時刻 t における推力 Ti (t) の surge(x)成分と sway(y) 成分をTix (t) = Ti (t) cos δi (t), Tiy (t) = Ti (t) sin δi (t)とすると，yaw モーメントは£¤·¸

Tix (t)−yi xiTiy (t)となる。したがって，推力の合成力は制御力に等しくなければならいので，次式が成り立つ。
T1x (t) + · · · + Tnx (t) = X(t)T1y (t) + · · · + Tny (t) = Y (t)N1x (t) + · · · + Nnx (t) = N (t)これの行列表現すると

T1x (t)  T (t)  
10 ···
10X(t)
 1y ..
 01 ···01 
 = Y (t) .
 Tnx (t) −y1 x1 · · · −yn xnN (t)
Ni (t) = −Tix (t)yi + Tiy (t)xi =Tny (t)となり，これを推力配分方程式と呼ぶ。
つぎに，アジマススラスタの首振角に制約をつけることを前提として，i 番目の取り付け角度を δi0 とし，これを基準とした首振角を θi (t) とすると

δi (t) = δi0 + θi (t)を得る。これから·¸ ·¸·¸cos δi (t)cos δi0 − sin δi0cos θi (t)=sin δi (t)sin δi0cos δi0sin θi (t)が成り立つことに注意してT̃ix (t) = Ti (t) cos δi (t), T̃iy (t) = Ti (t) sin δi (t)を定義すれば¸·¸ ·¸·T̃ix (t)Tix (t)cos δi0 − sin δi0=
T̃iy (t)Tiy (t)sin δi0cos δi0を得る。これを上の推力配分方程式に代入すれば，五百木ら 3)によって示された新しい推力配分方程式
ÃT̃ (t) = τ (t)(1)を得る。ここでT̃ (t) =£T̃1x (t)T̃1y (t)···T̃nx (t)T̃ny (t)¤T£¤Tτ (t) = X(t) Y (t) N (t)また，行列 Ã の第 2i −列は
1 列と第 2i ·10cos δi0Ã(: .[2i − 1, 2i]) =01
sin δi0−yi xi
2.2は，大坪ら 2) が近似解について，五百木ら 3) が厳密解について検討している。

− sin δi0cos δi0¸

3. 閉ループ系シミュレーション
3.11)
Fossenに従って，次の浮体運動方程式を考える。




cos ψ(t) − sin ψ(t) 0u(t)ẋ(t)
 ẏ(t)  =  sin ψ(t) cos ψ(t) 0   v(t) ψ̇(t)001r(t)推力配分に関する従来法と提案法

推力配分問題は，線形方程式 ÃT̃ (t) = τ (t) の求解問題に他ならない。その一般解は，Ã の特異値分解をÃ = U£|Σ1{z0¤£}|V1ΣV2{z¤Tm11 00}

VT

Σ1 = diag{σ1 , σ2 , σ3 }, (σ1 ≥ σ2 ≥ σ3 > 0)U T U = I3 , V T V = Inで表すとき，次式で与えられる。
T̃ (t) = V1 Σ1−1 U T τ (t) + V2 c

閉ループ系の構成

(2)ここで，c は 2n − 3 個の要素をもつ任意に選ばれたベクトルである。T̃ (t) の２ノルムはT22kT̃ (t)k2 = kΣ−11 U τ (t)k + kckのように計算できることから，２ノルム最小解は，c = 0 としてTT̃ ∗ (t) = V1 Σ−11 U τ (t)のように得られる。これは最適化問題½min kT̃ (t)k2
subject to ÃT̃ (t) = τ (t)を，ラグランジュの未定定数法を用いて解いた結果T̃ ∗ (t) = ÃT (ÃÃT )−1 τ (t)と一致する。すなわち，アジマススラスタの首振角に制約をつけない場合のラグランジュ法は，望ましい制御力 τ (t) に疑似逆行列 ÃT (ÃÃT )−1 をかけて得られる。
これに対して，次のアジマススラスタの首振角制約を考えた場合の最適化問題は，次のように表わされる。
2
 min kT̃ (t)ksubject to ÃT̃ (t) = τ (t)|θi (t)| ≤ θ̄i ≤ π2 (i = 1, · · · , n)となる。いま，i = 1, · · · , n に対して
|θi (t)| ≤ θ̄i ⇔ Tix (t) > 0, |Tiy (t)| ≤ Tix (t) tan θ̄iおよび，(2) 式の 2i − 1 要素と 2i 要素の表現式

0m22m32

0u̇(t)m23   v̇(t)  +m33ṙ(t)







d11
00X(t)u(t) 0 d22 d23   v(t)  =  Y (t) N (t)0d32 d33r(t)これに対する DPS(Dynamic Positioning System) の設計は非線形制御理論に基づくべきであるが，ここでは推力配分法を閉ループ系シミュレーションを通して評価することが目的であり，DPS の制御性能を厳しく達成するには及ばない。そこで，簡単のため，surge,sway,yaw の各運動に対する PID 制御を用いる。
すなわち、目標値との偏差
ex (t) = xc (t) − x(t)ey (t) = yc (t) − y(t)eψ (t) = ψc (t) − ψ(t)に基づく，次のＰＩＤ制御を考える。
Rtdex (t) − k13 0 ex (τ )dτX(t) = −k11 ex (t) − k12 dt
RtdY (t) = −k21 ey (t) − k22 dt
ey (t) − k23 0 ey (τ )dτ
RtdN (t) = −k31 eψ (t) − k32 dt eψ (t) − k33 0 eψ (τ )dτ
当該浮体に設置されるアジマススラスタの数は３個とする。i番目のスラスタの推力 Ti (t) は，スラスタ軸の回転数 ni (t) ≥ 0との間に，次式が成り立つものとする。
Ti (t) = ki ni (t)|ni (t)| = ki n2i (t) (i = 1, 2, 3)
また，スラスタ軸の回転数 ni (t) と，前章で定義した首振角 θi (t)は，次の動特性をもつとする。
ṅi (t) = − T1ni ni (t) + T1ni uni (t) (i = 1, 2, 3)

·

¸·

¸·

¸ ·

¸

θ̇i (t)01
θi (t)
0=+uωi (t)2θ̈i (t)−ωi −2ζi ωi
θ̇i (t)ωi2
 uT τ (t) 1(i = 1, 2, 3)¸·¸· (1)σ(1)(1)
Tix (t)v2i−1,1 v2i−1,2 v2i−1,3  uT2 τ1(t) 
以下では，パラメータ値と目標値は，問題の特徴を再現できれ=


(1)
(1)
(1)
σ2v2i,1v2i,2v2i,3Tiy (t)uTばよいので，次の仮想的なものを用いた。
3 τ (t)σ3 m11 = m22 = m33 = 1, m12 = m21 = 0· (2)¸ c1 (t)(2)
v2i−1,1 · · · v2i−1,2n−3  .
 d11 = d22 = d33 = 0, d12 = d21 = 0+ ..
 k = 1, k = √1 , k = 0 (i = 1, 2, 3)(2)(2)
i1i2i3v2i,1· · · v2i,2n−32
c2n−3 (t)ki = 10 (i = 1, 2, 3)に留意する。このとき，大坪ら 2) は，上の最適化問題を解く代Tni = 0.1, ζi = √12 , ωi = 0.1 (i = 1, 2, 3)√
わりに，次の不等式制約下の２次計画問題½, θ̄1 = π3(x1 , y1 ) = (− √23 , 0.5), δ10 = 2π
23minc(t) kc(t)k3, θ̄2 = π3(x2 , y2 ) = (− 2 , −0.5), δ20 = − 2π3
subject to Bc(t) ≤ d(t)π(x3 , y3 ) = (1, 0), δ30 = 0, θ̄3 = 2を解いて，c(t) を定め，(2) 式によって推力配分を行うことを提(5 ≤ t ≤ 10)
 0.2(t − 5)
案した。ここで，B と d(t) の表現式については割愛する。
xc (t) =1(10 ≤ t ≤ 40) , = 0 (else)
この提案法は，ある時刻 t において，市販のソルバ (たとえば−0.2(t − 40) + 1 (40 ≤ t ≤ 45)
MATLAB の関数 lsqlin) を用いれば，容易にオフライン計算(25 ≤ t ≤ 30)
 0.2(t − 25)できるが，実時間制御におけるサンプリング周期内に解けること
yc (t) =1(30≤ t ≤ 40) , = 0 (else)を保証することは困難である。提案法のオンライン計算について−0.2(t − 40) + 1 (40 ≤ t ≤ 45) π − 30 (t − 15)(15 ≤ t ≤ 20)
ψc (t) =(20 ≤ t ≤ 40) , = 0 (else)(40 ≤ t ≤ 45)
閉ループ系シミュレーションには，MATLAB/Simulink を用いてをた。制御信号の流れは，「偏差の計算⇒制御力の計算⇒推力配分の計算⇒スラスタ軸の回転数および動特性の計算⇒合成力の計算⇒浮体の動特性の計算」の順番となる。ここで，(2) 式の推力配分計算は，Fig1による。上側のパスが (2) 式の第１項（従来法）を，下側のパスが (2) 式の第２項（提案法による付加項）を表している。また，下側のパスのブロック内部を Fig2に示す。
これは，制御力の範囲
X ≤ X(t) ≤ X, Y ≤ Y (t) ≤ Y , N ≤ N (t) ≤ Nを仮定し，これらを適当にメッシュ分割し，あらかじめ最適な cベクトルをオフライン計算しておき，これを３次元補間するものである（解の非線形性が観察される）。
− π6π(t − 40) − π630Case2: Fig.5 と Fig.9 に，上記閉ループ系の推力配分計算に提案法を適用した場合のシミュレーション結果を，case 1 と対比して示す。回転数をそれほど増大させることなく，首振角の制約が実現されていることが観察され，提案法が有効に働いていることが分かる。
ところが，この提案法を用いた平水中のある模型実験においては，スラスタの首振角は所定の範囲内に収まるが，それでもなおその範囲内で頻繁に首振りを行うことが観察されていた。この状
況を再現するために，センサ出力に正弦波信号を加え，その振幅と周期をいろいろ変えてみた。

Case3: Fig.6 と Fig.10 に，case 2 と同じ状況下で，浮体応答 x(t)，y(t)，ψ(t) に, それぞれ正弦波信号 0.005 sin 2πt，0.005 sin 2πt，0.0025 sin 2πt を加えた場合のシミュレーション結果を示す。明らかに制御力が正弦波信号の影響を受けており，微弱な推力であるが，頻繁な首振りを要求している。ただし，センサ出力に加えた正弦波信号は周期 1 秒であり，その物理的意味については検討の余地がある。
Case4: Fig.7 と Fig.11 に，case 3 と同じ状況下で，Fig.3を用いて，センサ出力に加えた正弦波信号の影響を受けた制御力をフィルタリングし，またデットゾーンを挿入する対策を講じた
場合のシミュレーション結果を示す。頻繁な首振りはある程度の緩和されているが，DPS の性能劣化が観察される。

Fig. 1 Thrust Allocation Block

Fig. 3 Filtering of Sensor Noise
4. 結 言本稿では，アジマススラスタ首振角に制約を課す場合の静的推力配分法を，閉ループ系シミュレーションを通して評価し，頻繁Fig. 2 Interpolation Block for the Optimal cな首振りが生じる場合があることを確認した。これを回避するた
めには，動的推力配分問題の定式化が必要であるが（ボトムアップ的アプローチ），同時に制御力を発生する DPS の設計問題に3.2

シミュレーションによる問題分析

Case1: Fig.4 と Fig.8 に，上記閉ループ系の推力配分計算に従来法を適用した場合のシミュレーション結果を示す。まず Fig.4 の左列は，目標値 xc (t)，yc (t)，ψc (t) および浮体応
答 x(t)，y(t)，ψ(t) を重ねて示している。PID 制御の制御性能が甘いため，少し振動気味であるが，浮体応答は目標値によく追従している。このときの制御力 X(t)，Y (t)，N (t) と合成力 Tx1 (t) + Tx2 (t) + Tx3 (t)，Ty1 (t) + Ty2 (t) + Ty3 (t)，N1 (t) + N2 (t) + N3 (t) を，Fig.4 の右列に重ねて示している。
つぎに，Fig.8 の左列に回転数 n1 (t)，n2 (t)，n3 (t)，右図に首振角 θ1 (t)，θ2 (t)，θ3 (t) を示している。首振角に制約を課していないため，スラスタは大きな回転を余儀なくされている。

ついても考察する必要がある（トップダウン的アプローチ）といえる。
参 考 文 献1) T.I.Fossen: Marine Control Systems, pp.268-300,
pp.417-449, Marine Cybernetics, 2002
2) 大坪・梶原：アジマススラスタ首振角に制約がある場合の推力配分法, 日本船舶海洋工学会論文集，第 6 号，pp.177-182,
20073) 五百木・梶原：区分的線形補間による首振角制限付きアジマススラスタの最適推力配分法, 日本船舶海洋工学会論文集，
第 6 号，pp.183-190, 2007

1.5

1.5

110.5

0.5

00-0.50-0.50204060xc(t) & x(t)102040-20060020

110.50.500-0.5
0-0.5
060

yc(t) & y(t)

4060

40

4060

200

500-20060

20

theta1(t)1020

0n1(t)

1.5

40

00X(t) & Sum of Txi(t)

1.5

20

200

5020Y(t) & Sum of Tyi(t)40600n2(t)

20

4060

theta2(t)

0.5

1005

200
0-0.5
00-200

0-20
-40
020

4060

psic(t) & psi(t)

20

40

60

0

1.5

1.5

110.5

0.5

00-0.50-0.5
04060

xc(t) & x(t)

20

40

60

110.5

00-0.5
0-0.5
020

4060

yc(t) & y(t)4060204060theta3(t)10200

500-200
0204060

0204060

theta1(t)10200n1(t)1.50.560

Fig. 8 Revolutions and Angles in Case1X(t) & Sum of Txi(t)1.540n3(t)Fig. 4 Conventional Method Simulation (Case1)

2020N(t) & Sum of Ni(t)200500-200
020Y(t) & Sum of Tyi(t)40600n2(t)204060theta2(t)

0.5

100

500-200

2000-20-400204060-0.50psic(t) & psi(t)20406001.54060

0n3(t)

Fig. 5 Proposed Method Simulation (Case2)

204060theta3(t)Fig. 9 Revolutions and Angles in Case2

1.5

110.50.50100-0.5020N(t) & Sum of Ni(t)204060-0.50xc(t) & x(t)204060200500-2000X(t) & Sum of Txi(t)2040600n1(t)
1.5

204060theta1(t)

1.5

1

1

0.5

0.5

00-0.50-0.5020

4060yc(t) & y(t)20406010200

500-2000Y(t) & Sum of Tyi(t)2040600n2(t)204060theta2(t)0.51002000-20-40020

4060-0.50psic(t) & psi(t)204060500-200
0N(t) & Sum of Ni(t)2040600n3(t)

Fig. 6 Proposed Method Simulation with a Sinusoidal Signal (Case3)
1.5

1

1

0.5

0.5

00-0.5
0

20

4060

-0.5
0xc(t) & x(t)
1.5

1

1

0.5

0.5

00
20

4060

-0.5
0yc(t) & y(t)200

500-200404060204060theta1(t)10200

500-2002040

60

0n2(t)

60

204060theta2(t)10Y(t) & Sum of Tyi(t)0n1(t)

0202060

X(t) & Sum of Txi(t)

1.5

-0.5
04060

Fig. 10 Revolutions and Angles in Case3

02040

theta3(t)

101.5

20200

0.5
0500-200

0-200-4002040

psic(t) & psi(t)

60

-0.5
02040

60

2040n3(t)

60

02040

theta3(t)N(t) & Sum of Ni(t)Fig. 7 Proposed Method Simulation with DeadZone (Case4)

Fig. 11 Revolutions and Angles in Case460