線形パラメータ変動制御技術の海洋工学分野
への応用に関する研究
九州大学大学院工学府海洋システム工学専攻
大坪 和久
2006 年 2 月
目次
第 1 章 緒論
1.1 研究背景 . . . . . . . . . . . . . . . . . . . . .
1.1.1 風車の回転数制御による高効率化 . . .
1.1.2 ライザー管リエントリー作業の自動化
1.2 本論文の構成と概要 . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
第 2 章 LMI ベース LPV 制御系設計法について
2.1 はじめに . . . . . . . . . . . . . . . . . . . . . . .
2.2 LMI の概要 . . . . . . . . . . . . . . . . . . . . .
2.2.1 凸性と凸計画問題 . . . . . . . . . . . . . .
2.2.2 LMI の定義 . . . . . . . . . . . . . . . . .
2.3 アナリシス LMI . . . . . . . . . . . . . . . . . . .
2.3.1 安定条件 . . . . . . . . . . . . . . . . . . .
2.3.2 H∞ 性能条件 . . . . . . . . . . . . . . . .
2.3.3 極配置制約条件 . . . . . . . . . . . . . . .
2.4 シンセシス LMI . . . . . . . . . . . . . . . . . . .
2.4.1 状態フィードバック . . . . . . . . . . . .
2.4.2 出力フィードバック . . . . . . . . . . . .
2.5 ゲインスケジューリング制御の考え方 . . . . . . .
2.6 LPV システムの構造上の特徴と性能解析 . . . . .
2.6.1 構造上の特徴 . . . . . . . . . . . . . . . .
2.6.2 性能解析 . . . . . . . . . . . . . . . . . . .
2.7 LPV モデルに基づくゲインスケジューリング制御
2.8 まとめ . . . . . . . . . . . . . . . . . . . . . . . .
第3章
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
固定ピッチ型水平軸風車の可変速制御による高効率化
はじめに . . . . . . . . . . . . . . . . . . . . . . . . .
風力発電機の運転モード . . . . . . . . . . . . . . . .
風力発電機の状態方程式の導出 . . . . . . . . . . . .
3.3.1 運動の数学モデル . . . . . . . . . . . . . . . .
3.3.2 動作点周りでの線形化 . . . . . . . . . . . . .
最大発電効率を得るための最適回転数制御 . . . . . .
定格出力維持のための最適回転数制御 . . . . . . . . .
風力発電機に対する LPV 制御系設計 . . . . . . . . .
数値シミュレーション . . . . . . . . . . . . . . . . .
風車シミュレータを用いた検証実験 . . . . . . . . . .
3.8.1 風車シミュレータの考え方と装置の概要 . . .
3.8.2 検証実験 . . . . . . . . . . . . . . . . . . . . .
まとめ . . . . . . . . . . . . . . . . . . . . . . . . . .
1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
7
7
7
8
8
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
12
12
13
13
14
14
15
15
15
16
16
18
23
25
25
26
28
29
.
.
.
.
.
.
.
.
.
.
.
.
.
36
36
37
38
38
39
39
40
42
43
44
44
44
45
第 4 章 大水深ライザー管の動特性を考慮したリエントリー作業の自動化
4.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 大水深ライザー管の状態方程式の導出 . . . . . . . . . . . . . . . . .
4.2.1 運動の数学モデル . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 モード展開による運動近似化 . . . . . . . . . . . . . . . . .
4.3 大水深ライザー管のリエントリー制御 . . . . . . . . . . . . . . . . .
4.3.1 大水深ライザー管の LPV モデルの導出 . . . . . . . . . . . .
4.3.2 固定パラメータ開ループ系モデルの動特性解析 . . . . . . . .
4.3.3 大水深ライザー管の LPV コントローラの導出 . . . . . . . .
4.4 数値シミュレーション . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 水槽実験による検証 . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.1 実験の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.2 実験内容とその妥当性について . . . . . . . . . . . . . . . .
4.5.3 実験結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6 まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
第 5 章 結言
付 録A
A.1
A.2
A.3
A.4
A.5
A.6
A.7
55
55
56
56
57
57
57
58
59
59
60
60
60
61
61
74
Nomoto モデルに基づく船舶の回頭角制御について
はじめに . . . . . . . . . . . . . . . . . . . . . . . . . .
船舶の非線形操縦運動方程式 . . . . . . . . . . . . . . .
Nomoto モデルの動特性とステップ応答解析 . . . . . .
船舶の PID 制御則と GSPID 制御則の問題点 . . . . . .
有次元 Nomoto モデルに基づく船舶の LPV 制御系設計
数値シミュレーション . . . . . . . . . . . . . . . . . .
まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
付 録 B 大水深ライザー管の運動方程式の導出
B.1 運動エネルギー (T ) . . . . . . . . . . . . . . . . . . . . . . .
B.2 ポテンシャルエネルギー (U ) . . . . . . . . . . . . . . . . . .
B.2.1 位置エネルギー (Ugrav ) . . . . . . . . . . . . . . . .
B.2.2 軸方向の変形 (伸び) によるひずみエネルギー (Uaxial )
B.2.3 曲げによるひずみエネルギー (Ubend ) . . . . . . . . .
B.2.4 張力によるひずみエネルギー (Utens ) . . . . . . . . .
B.3 仮想仕事 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.3.1 流体力による仮想仕事 (Wf luid ) . . . . . . . . . . . .
B.3.2 上端の荷重による仮想仕事 (Wupper ) . . . . . . . . .
B.4 各種エネルギーの変分の導出 . . . . . . . . . . . . . . . . . .
B.4.1 運動エネルギーの変分 (δT ) . . . . . . . . . . . . . .
B.4.2 ポテンシャルエネルギーの変分 (δU ) . . . . . . . . .
B.5 Hamilton の原理による運動方程式の導出 . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
76
76
76
77
79
80
81
81
.
.
.
.
.
.
.
.
.
.
.
.
.
89
89
90
90
90
90
91
91
91
93
93
93
94
95
97
参考文献
100
謝辞
2
図目次
1.1
1.2
1.3
1.4
Wind farm of California in the United States of America
Offshore type wind farm of Middelgrunden in Denmark
Drilling vessel ’CHIKYU’ (JAMSTEC) . . . . . . . . . .
Problem of riser reentry operation . . . . . . . . . . . .
.
.
.
.
10
10
11
11
2.1
2.2
Convex set and Nonconvex set (Left side: Convex set, Right side: Nonconvex set) . . .
Convex function and Nonconvex function (Left side: Convex function, Right side:
Nonconvex function) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
α-stability region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Circle-stability region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Conic-sector-stability region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Extended LMI Region from single LMI region . . . . . . . . . . . . . . . . . . . . . . .
Gain-scheduled control system structure . . . . . . . . . . . . . . . . . . . . . . . . . .
Constrainted parameter region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Stable region of parameter varying system . . . . . . . . . . . . . . . . . . . . . . . . .
Single model based robust control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Multi model based robust control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Gain-scheduled control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Class of linear systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Polytopic type LPV system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
LFT type LPV system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
30
Operational modes of wind turbine . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
”Wind Lens” developed by Research Institute for Applied Mechanics in Kyushu University . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Wind turbine system model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Power coefficient CP and torque coefficient CT curves of Wind Lens . . . . . . . . . .
Target operational point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
dCT
Relationship between ρπr
2JT V dλ and wind speed V in mode 3 . . . . . . . . . . . . . .
Approxmated torque coefficient CT curve of Wind Lens . . . . . . . . . . . . . . . . .
LPV system varying with rotational speed ω and wind speed V in mode 3 . . . . . . .
Description of the LPV system using two scheduling parameters . . . . . . . . . . . .
Interconection for the mixed sensivity problem . . . . . . . . . . . . . . . . . . . . . .
Comparison of control performance between conventional ω 2 control and LPV control
in simulational results ((a):mode 2 region, (b):mode 3 region) . . . . . . . . . . . . . .
Wind Turbine Simulator (photo 1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Wind Turbine Simulator (photo 2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Wind Turbine Simulator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Caluculation flow of aerodynamical torque in the Wind Turbine Simulator . . . . . . .
Experimental results in mode 2 region using conventional ω 2 controller . . . . . . . . .
47
2.3
2.4
2.5
2.6
2.7
2.8
2.9
2.10
2.11
2.12
2.13
2.14
2.15
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10
3.11
3.12
3.13
3.14
3.15
3.16
3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
30
31
31
31
32
32
33
33
34
34
34
35
35
35
47
47
48
48
48
49
49
49
50
50
51
51
52
52
53
3.17 Experimental results in mode 2 region using LPV controller . . . . . . . . . . . . . . .
3.18 Experimental results in mode 3 region using conventional ω 2 controller . . . . . . . . .
3.19 Experimental results in mode 3 region using LPV controller . . . . . . . . . . . . . . .
53
54
54
4.1
4.2
Coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Simulation results of marine riser’s free motion (solid line: with hydrodynamical damping (Cd = 1.17), dotted line: without hydrodynamical damping (Cd = 0.0)) . . . . . .
Bode plots of open-loop marine riser system (solid line: |ṙ| = 0, dash line: |ṙ| = 1.5) .
Pole placement of open-loop marine riser system (×: |ṙ| = 0, ○: |ṙ| = 1.5) . . . . . .
Control system structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Concept of controller system design in reentry operation of riser . . . . . . . . . . . . .
Bode plots of closed-loop marine riser system (solid line: |ṙ| = 0, dash line: |ṙ| = 1.5) .
Pole placement of closed-loop marine riser system (×: |ṙ| = 0, ○: |ṙ| = 1.5) . . . . . .
Comparison between MmBR control and LPV control in numerical simulation ((a):top
point position, (b):connected riser’s angle, (c):bottom point position, (d):controller
output) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Simulation results of the marine riser’s sway motion in the reentry operation ((a):MmBR
control result, (b):LPV control result) . . . . . . . . . . . . . . . . . . . . . . . . . . .
Time histories of 1∼8 mode amplitudes in numerical simulation (solid line: LPV control
result, dotted line: MmBR control result) . . . . . . . . . . . . . . . . . . . . . . . . .
Riser model in Research Institute for Applied Mechanics, Kyushu University . . . . .
Parallel mechanism to control the upper end of the riser model in Research Institute
for Applied Mechanics, Kyushu University . . . . . . . . . . . . . . . . . . . . . . . . .
Arrangement of measuring marker . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Overview of CCD camera frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Calibration marker’s arrangement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Camera coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Measurement flow of the motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Experiment of reentry operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Adequacy of this experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Experimental results of riser’s sway motion in the reentry operation . . . . . . . . . .
Experimental results of XY-plane motion in the reentry operation ((a): MmBR control,
(b): LPV control) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Time histories of bottom point motion in control direction . . . . . . . . . . . . . . . .
Time histories of top end forces Fx and Fy in the reentry operation ((a): MmBR
control, (b): LPV control) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
Coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Block diagram of ship control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Step response of ship δc = 5 [deg] . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Step response of ship δc = 10 [deg] . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Bode diagram of Nomoto’s model varying with ship forward speed . . . . . . . . . . .
Concept of dimensional Nomoto’s model . . . . . . . . . . . . . . . . . . . . . . . . . .
Concept of controller system design based on dimensional Nomoto’s model . . . . . . .
Time histories of ship steering motion in the dimensional Nomoto’s and the nonlinear model simulation ((a):dimensional Nomoto’s model, (b):nonlinear model, speed
variation: case 1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
83
83
83
84
84
84
4.3
4.4
4.5
4.6
4.7
4.8
4.9
4.10
4.11
4.12
4.13
4.14
4.15
4.16
4.17
4.18
4.19
4.20
4.21
4.22
4.23
4.24
A.1
A.2
A.3
A.4
A.5
A.6
A.7
A.8
4
63
64
64
64
65
65
65
66
66
67
68
68
69
69
70
70
71
71
72
72
73
73
73
85
A.9 Time histories of ship steering motion in the dimensional Nomoto’s and the nonlinear model simulation ((a):dimensional Nomoto’s model, (b):nonlinear model, speed
variation: case 2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.10 Time histories of ship steering motion in the dimensional Nomoto’s and the nonlinear model simulation ((a):dimensional Nomoto’s model, (b):nonlinear model, speed
variation: case 3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.11 Speed variation pattern . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.12 Control system structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.13 Triangle region under the varying ship speed . . . . . . . . . . . . . . . . . . . . . . .
A.14 Time histories of ship steering motion in the dimensional Nomoto’s and the nonlinear model simulation ((a):dimensional Nomoto’s model, (b):nonlinear model, speed
variation: case 1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.15 Time histories of ship steering motion in the dimensional Nomoto’s and the nonlinear model simulation ((a):dimensional Nomoto’s model, (b):nonlinear model, speed
variation: case 2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.16 Time histories of ship steering motion in the dimensional Nomoto’s and the nonlinear model simulation ((a):dimensional Nomoto’s model, (b):nonlinear model, speed
variation: case 3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
85
86
86
86
87
87
88
88
表目次
3.1
3.2
3.3
Specifications of Wind Lens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Controller performance in numerical simualtion . . . . . . . . . . . . . . . . . . . . . .
Control performance index in experiment using Wind Turbine Simulator . . . . . . . .
38
43
44
4.1
4.2
Controller performance in numerical simulation . . . . . . . . . . . . . . . . . . . . . .
Characteristics of the riser model in Research Institute for Applied Mechanics, Kyushu
University . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59
A.1 Main data of Mariner Class Vessel . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
77
6
60
第1章
緒論
1.1
研究背景
産業分野では、非線形な時変系のように動特性が大きく変化する制御対象に対する有効な手段として、
ゲインスケジューリング (Gain Scheduling, 以後, GS と略記) 制御が以前から用いられてきた。GS 制御
とは非線形な時変系に対して、コントローラをいくつか前もって用意し、あるパラメータ (スケジューリ
ング変数) の変動に対して、コントローラを随時変えていく制御方式であり、応用例の報告は数多く行
われている。しかしながら、この方法の一つの問題点は、閉ループ系の時変系としての安定性を理論的
に保証できない点にある。1990 年代から急速に発展した線形行列不等式 (Linear Matrix Inequality, 以
後, LMI と略記) ベース制御系設計法の登場により、時変系の安定性を理論的に保証することが可能に
なった。具体的には、制御対象の運動方程式から線形パラメータ変動 (Linear Parameter Varying, 以
後, LPV と略記) モデルを導出し、LMI を使ってゲインスケジューリングコントローラを設計する。本
論文では、この制御方式を LPV 制御と呼ぶことにする。この LPV 制御は、理論的な研究が主に行われ
てきており、近年において、その技術を応用する段階に入りつつある [1][2]。
海洋工学分野における制御技術の特徴は、制御対象が流体力学と絡み合った複雑なシステムであるこ
とである。この点が他分野における制御問題と大きく異なる点である。それらの取り扱いを適切に行う
ことにより、対象のより良い運動制御性能を実現することが可能である。もちろん従来から、海洋工学
分野においても、適応制御などの Lyapunov 関数ベースの制御系設計法に関する研究が数多く行われて
いる。しかし、それらの制御手法は、設計を解析的に行っていくために、複雑な数学や安定性の知識を
必要とし、流力などがシステム全体に対する影響などは不明瞭なままである。一方で、LPV 制御は線形
制御理論をベースとするために、難解な知識を改めて必要とはせず、成熟された線形理論をそのまま利
用することができる。また、設計の際には LMI を利用し数値的に処理するために、深い制御工学の知
識を有さない技術者でも、MATLAB などの制御系 CAD がインストールされた PC 環境が手元にあれ
ば、それらを援用することにより、適切なコントローラを容易に設計することが可能である。したがっ
て、海洋工学分野の諸問題に対して、LPV 制御技術の有効性や問題点を明確化し、その設計技術を確
立させていくことは非常に意義あることと考えられる。そこで本論文では、以下に示す 2 つの海洋工学
分野における問題に取り組むことにする。
1.1.1
風車の回転数制御による高効率化
石油、天然ガス、石炭などの化石燃料の枯渇問題が目前に迫ってきている中で、近年、新しい代替エ
ネルギーとして自然エネルギーが非常に注目されるようになり、世界各国で開発や導入が急速に進めら
れている。自然エネルギーとは、風力、太陽光、水力、地熱、海洋温度差などを指し示すことが多く、
その名の通り、自然なエネルギー源であるために、地球上に無尽蔵に存在するため枯渇しない。また、
環境に対する負荷が非常に小さく、二酸化炭素の排出量が非常に少ないことが大きな利点である。自然
エネルギーの中でも、風力エネルギーは非常に魅力的なエネルギーである。世界各地域の風力発電のポ
テンシャル量の総計は、498, 000TWh/年 と日本の電力供給量 91.960TWh/年 の 5400 倍もの可能性を
秘めていること、他の自然エネルギーと比較しても発電コストが相対的に安いことや大容量であること
7
が有利な点である [3]。したがって、今後、急速に風力エネルギー利用が普及していくものと期待され
る。Fig.1.1 には、カリフォルニア州に設置されたウィンドファームを示す。
すでに我が国においても、新エネルギー・産業技術総合開発機構 (New Energy and Industrial Technology Development Organization, 以後, NEDO と略記) を中心として、風力発電に関する研究開発が
活発に行われている。また最近では、Fig.1.2 に示すように風力発電設備の設置場所が陸上から洋上へと
シフトする傾向にある。その理由として、陸地と比較して洋上には設置可能スペースが多く存在するこ
と、ロータ翼が回転することにより発生する騒音問題などを考慮する必要がないこと、我が国の造船技
術は超大型浮体式構造物 (メガフロート) などの他国にはない独創的な技術を有すること、そして何よ
りも、洋上には風力発電に適した強い風が吹いているために、安定した発電が期待できることなどが挙
げられる。そのような安定した風速が得られる環境下においては、わずかな発電効率向上も、長期的に
は大きな効果が期待できるため、風力発電の高効率化に関する研究がますます重要になってきている。
制御技術の観点から風力発電の効率改善のためには、2 つの問題を解決する必要がある。第 1 に、定
格風速以下の風速領域において、風エネルギーから最大のパワーを獲得すること。第 2 に、定格風速以
上の風速領域においては、風のエネルギーの流入を調節しながら定格出力に維持することである。これ
らの問題を解決する 1 つの方法として、ロータの回転数を風速に応じた最適値に制御する方法がある。
しかし、風速そのものが予測不可能なことや、空力特性などに起因する動力学の非線形性などから、現
在のところ、この問題を解決する技術は確立されていない。ゆえに、今後の風力エネルギー利用促進の
ためにも、高効率化を実現できる制御技術を提案することが、大いに求められている。
1.1.2
ライザー管リエントリー作業の自動化
現在、地球環境問題、地震発生メカニズムの解明を目指して、水深 3000 m より深い深海の海底掘削や
資源採取を目指した統合国際深海掘削計画 (Integrated Ocean Drilling Program, 以後, IODP と略記) が
進められている。この計画の前に行われていた国際深海掘削計画 (Ocean Drilling Program, 以後, ODP
と略記) には、世界 22ヶ国が参加し、恐竜絶滅の原因とされる小惑星の衝突を実証するなど輝かしい成
果を挙げてきた。しかし、当時の技術では、海底下 2000m の掘削が限度であり、マントルや巨大地震
発生域などの深部ターゲットへ到達することが不可能だった。そこで、我が国では科学技術に貢献すべ
く、海洋研究開発機構 (Japan Agency for Marine-Earth Science and Technology、以後、JAMSTEC
と略記) が中心となって、海底下 7000 m までの海底掘削技術を要する Fig.1.3 に示すような地球深部探
査船「ちきゅう」(CHIKYU) を建造し、2005 年 7 月に運航を開始させた。
地球深部探査船「ちきゅう」に搭載されているライザー管は、最長 7000 m 以上にもなるため、相対
的な剛性が低下し、まるで釣り糸のような振る舞いを示す。また当然のごとく、全体が海中に没してい
るために流体力の影響を強く受ける。ゆえに、そのライザー管の複雑な挙動を予測することはもちろん
のこと、制御することは非常に困難な技術である。しかしながら、外乱要因を多く持つ海上での作業は
緊急退避などの必要に迫られるため、そのような状況下では、一時的にライザー管を接続している防
噴装置頂部から切り離し、回避後は早急に再接続する必要がある。それをリエントリー作業と呼ぶ。そ
の際、Fig.1.4 に示すように、掘削作業者は ROV(Remotely Operated Vehicle) に備え付けられたカメ
ラで状況を観察し、長年の経験で培った勘に頼りながら船舶を位置制御し、先端部を目標点まで制御す
る。しかしながら、現状では、このリエントリー作業に多大な労力を要しており、さらに、掘削者の技
術継承問題などを考慮すると、高精度ライザー管先端位置制御によるリエントリー作業の自動化は早急
に求められている必要不可欠な技術である。
1.2
本論文の構成と概要
本論文は、線形パラメータ変動制御技術を「固定ピッチ型水平軸風車の可変速制御による高効率化」
と「大水深ライザー管の動特性を考慮したリエントリー作業の自動化」の 2 つの海洋工学分野の問題に
8
適用し、数値シミュレーションと実験によりその有効性を検討したものである。
本論文は以下の 5 章により構成されている。
第 1 章の「緒論」では、本論文の研究背景、論文構成と概要について述べる。
第 2 章の「LMI ベース LPV 制御の考え方」では、まず、LMI の概要について簡単に述べ、幾つかの
制御仕様を LMI により定式化する。次に、その定式化された LMI を使用して、極配置制約付きの H∞
制御系設計法について、状態フィードバック、出力フィードバックの場合について説明する。そして、
ゲインスケジューリング制御の考え方、及び、その具体的な制御手法について説明する。その中で特に、
LPV 制御方式に着目し、LPV システムの構造上の特徴と性能解析について具体的に説明をする。最後
に、LMI ベース LPV 制御系設計法について述べ、例題として船舶に対する応用例を示す。
第 3 章の「固定ピッチ型水平軸風車の可変速制御による高効率化」では、まず、風力発電の現状と課
題、及び、従来研究について述べる。次に、風力発電機の運転モード、対象とする風力発電機の数学モ
デルを説明し、動作点周りで線形化することにより状態方程式を導出する。そして、2 つの運転モード
に着目し、それぞれのモードに対応する LPV モデルを求め、LPV コントローラの導出を行う。コント
ローラの性能検証のために数値シミュレーション、及び、本論文で提案した風車シミュレータと呼ばれ
る風力発電機開発支援装置による検証実験により有効性を確認し、最後に研究のまとめと今後の課題を
述べる。
第 4 章の「大水深ライザー管の動特性を考慮したリエントリー作業の自動化」では、まず、ライザー
管のリエントリー作業の現状とその必要性、及び、従来研究について述べる。次に、大水深ライザー管
の運動方程式を Hamilton の原理から導出し、モード展開法と微小変位の仮定を施すことにより LPV モ
デルを求める。そして、その LPV モデルの動特性解析を行った後、LPV コントローラを導出する。性
能検証のために、数値シミュレーション、及び、九州大学応用力学研究所で実施した水槽実験により有
効性を確認し、最後に研究のまとめと今後の課題を述べる。
第 5 章の「結言」では、各章で得られた研究成果をまとめて本論文の結論とする。
9
Fig. 1.1: Wind farm of California in the United States of America
(URL: http://www.mhi.co.jp/nsmw/mwt/jp/moj89j.html)
Fig. 1.2: Offshore type wind farm of Middelgrunden in Denmark
(URL: http://www.yannarthusbertrand.com/yann2/affichage.php?pais=Danemark)
10
Fig. 1.3: Drilling vessel ’CHIKYU’ (JAMSTEC)
(URL: http://www.nishinippon.co.jp/)
Riser pipe
Target
Image
Information
Operation
Fig. 1.4: Problem of riser reentry operation
(URL: http://www.ep.sci.hokudai.ac.jp/)
11
第2章
LMI ベース LPV 制御系設計法について
本章では、本論文の基礎となる LMI ベース LPV 制御方式について述べる。本章の構成は以下の通り
である。はじめに、LMI の概要について簡単に述べ、極配置や H∞ 性能条件などの幾つかの制御仕様
を LMI により定式化する。次に、その定式化された LMI を使用して、極配置制約付きの H∞ 制御系設
計法について、状態フィードバック、出力フィードバックの場合について説明する。そして、ゲインス
ケジューリング制御の考え方、及び、その具体的な制御方式について説明する。その中で、特に LPV
モデルに基づいてゲインスケジューリングコントローラを設計する LPV 制御に着目し、LPV システム
の構造上の特徴と性能解析について具体的に説明をする。最後に、LMI ベース LPV 制御系設計法につ
いて述べ、例題として船舶の回頭角制御に対する応用例を示す。
2.1
はじめに
GS 制御とは非線形な時変系に対して、コントローラをいくつか前もって用意し、スケジューリング
変数の変動に対して、コントローラを随時変えていく制御方式である。しかし、単純にコントローラを
いくつか用意しておくだけで良いというものでなく、スケジューリング変数 θ(t) の変化速度の大きさに
よっては GS 制御では対応てきないことがあるので、モデリングの仕方などに多大な労力を要すること
になる。GS 制御について、一般的に、変化速度 θ̇(t) が小さいほど、各点 θ(t) での固定パラメータ線形
システム P (θ) に対して、局所的に保証された制御性能が θ(t) を可変とした非線形 GS 制御系でも成立
しやすい。同じコントローラを用いても、θ̇(t) が大きいシステムに対してはノミナルシステムの周りの
安定領域が狭くなり、小さいシステムはその反対である。従って、θ̇(t) が大きいシステムでは安定領域
の大きさを保つために各点での局所的な線形制御系設計でより厳しい安定性の仕様が要求されることに
なる。このように GS 制御系設計では、θ̇(t) の大きさを考慮して局所的な線形コントローラを設計する
ことが、全体の非線形 GS 制御系の性能を達成するために重要である。
産業分野では、動特性が大きく変化する制御対象に対する有効な手段として、GS 制御が以前から用
いられてきた。この方法については、数多くの応用例の報告がなされている。しかしながら、その方法
の一つの問題点は、閉ループ系の時変系としての安定性を理論的に保証できないために、確固たる設
計論が無いことである。ゆえに設計者は、問題に応じて試行錯誤により設計法を工夫してきた。そんな
中、1990 年代から急速に発展した LMI ベース制御系設計法の登場により、時変系の安定性を理論的に
保証することが可能になった。具体的には、制御対象の運動方程式から LPV モデルを導出し、LMI を
使ってゲインスケジューリングコントローラを設計する LPV 制御は、航空機や産業ロボット分野など
において、近年、数多くの報告が行われてきている [4][5][6]。しかしながら、著者が知る限りでは、海
洋工学分野に応用した例はあまり見られない [7][8]。
そこで本章では、海洋工学分野の問題に対して LMI ベース LPV 制御を適用するために、その周辺の
基礎知識について整理することが目的である。
12
2.2
LMI の概要
制御系設計において、複数の設計仕様を考慮したい場合がある。その際、その設計仕様を表現する拘
束条件が等式であれば、条件式が増すにつれて、満足するコントローラを設計することは非常に困難
になり、時には可解でなくなる事さえもありうる。しかし、その拘束条件を不等式に緩和することによ
り、その可解性は劇的に向上する。LMI ベース制御系設計とは、課せられる拘束条件を線形の行列不等
式で表現してコントローラを数値的に設計するという方法論である。LMI は以前から線形制御理論な
どではよく登場するありふれた考え方であったが、近年、それを解く計算アルゴリズムが改良されたこ
とや、それを解くソルバーを有する制御系 CAD が整備されてきたことなどから、大きな注目を浴びる
ようになってきた [9]。LMI とは、その名が示す通り、各要素が変数に線形に依存する行列に関する不
等式である。そして、なによりも、この不等式により記述された凸計画問題は必ず解けるという点が大
きな特徴である。つまり、凸計画法アルゴリズムを使って、数値的に解を求めることが可能である。し
たがって、制御系の解析と設計問題において、問題を LMI 問題へと書き直すことが可能であるならば、
効率的な解析、及び、設計が可能となる。このような事情から、最近では、LMI ベースの制御系設計が
非常に主流となりつつある。
2.2.1
凸性と凸計画問題
まず始めに、LMI を理解する上で必要な凸性と凸計画問題について、その考え方を説明する [10][11]。
定義 2.1 凸集合
集合 S は、以下の条件を満たすとき凸集合であるという。
x∈S, x̂∈S → αx + (1 − α)x̂ ∈ S, ∀ 0 ≤ α ≤ 1
(2.1)
■
上記の定義を理解するために、集合 S が R2 の時について考える。R は実数空間を意味する。要素ベ
クトル x ∈ S, x̂ ∈ S に対して、αx + (1 − α)x̂ は、その 2 点を結ぶ直線上で x と x̂ を α : 1 − α に内分
する点を表す。したがって、α が 0 から 1 まで変化させるとき、その軌跡は 2 点を結ぶ線分となる。つ
まり、任意の要素 x∈S、及び x̂∈S に対してそれらを結ぶ線分が集合 S に含まれるとき、S は凸集合であ
る。その考え方を Fig.2.1 に示す。集合 S1 、及び S2 が凸ならば、その交わり S1 ∩ S2 もまた凸集合であ
る。以上の考え方を関数にまで拡張すると、次のような凸関数を定義をすることが可能である。
定義 2.2 凸関数
関数 Φ : S → R は、以下の条件を満たすとき凸関数であるという。
Φ(x̄) ≤ αΦ(x) + (1 − α)Φ(x̂), ∀ 0 ≤ α ≤ 1, x∈S, x̂∈S
(2.2)
x̄ := αx + (1 − α)x̂
(2.3)
ただし、
■
この考え方を理解するため、関数の引数をスカラー実数に限定して考える。このとき x は、0 ≤ α ≤ 1
により横軸上 x, x̂ の間を移動する。このとき、αΦ(x) + (1 − α)Φ(x̂) は α を媒介変数として x, x̂ を結
ぶ線分を意味する。その際、x の値に関わらず、その線分が関数 Φ(x) よりも上にある時、その関数を
凸関数と定義する。その考え方を Fig.2.2 に示す。
以上で述べてきた凸集合と凸関数を基にすると、次のような凸計画問題を考えることができる。
13
凸計画問題 : 与えられた凸関数 Φ(x) と凸集合 S に対して、次の最適値 γ ∗ とそれを与える x∗ を求めよ。
γ ∗ = min Φ(x)
(2.4)
x∈S
この凸計画問題の大きな特徴は、数値計算などを利用すれば必ず解けるということである。また、一
般的に局所最適解などを見つける計算アルゴリズムは容易に作ることが出来るが、凸計画問題にそのア
ルゴリズムを利用すると、局所最適解が必ず大域的最適解になる。もし、凸関数 Φ(x)、凸集合 S のい
ずれか、及び両方がその凸性を失うと、局所最適解を求めてもその大域性を満足することは出来ないの
で、その中で大域的最適解を見つけることは非常に困難になる。
2.2.2
LMI の定義
LMI の一般形は次のような不等式で記述される。
F (x),F0 + x1 F1 + x2 F2 + · · · + xm Fm > 0
(2.5)
ここで、Fi (i = 0, 1, · · · , m) は与えられた対称行列、x,[x1 , · · · , xm ]T は変数ベクトルである。また、不
等号 > は行列関数 F (x) が正定、つまり最小固有値が正であることを意味している。この LMI の大き
な特徴は、凸性を有するために、効率良く数値的に大域的最適解を求めることが可能である点にある。
すなわち、F (x) > 0 は x に対して凸の制約条件を課すものであり、
F , {x : F (x) > 0}
(2.6)
は凸集合となる。この性質を理解するために次のように考える。x ∈ F と x̂ ∈ F を任意の要素とする。
すると、
F (x) > 0, F (x̂) > 0
(2.7)
が成立する。次に、α ∈ [0, 1] を任意に固定すると、
F (αx + (1 − α)x̂) = F0 +
m
X
(αxi + (1 − α)x̂i )Fi
i=1
= α(F0 +
m
X
xi Fi ) + (1 − α)(F0 +
i=1
m
X
x̂i Fi )
i=1
= αF (x) + (1 − α)F (x̂) > 0
(2.8)
となり、αx + (1 − α)x̂ ∈ F となり、定義 2.1 より確かに F は凸集合である。
2.3
アナリシス LMI
ここでは、次のような p 出力 m 入力を持つ n 次の安定な連続時間線形時不変システムを考える。
½
ẋ(t) = Ax(t) + Bu(t)
y(t) = Cx(t) + Du(t)
この線形システムに対して、主な制御仕様の LMI 定式化を以下に示す [12][13]。
14
(2.9)
2.3.1
安定条件
(2.9) 式の入力 u(t) をゼロとした線形自律系の安定条件は、次の Lyapunov 関数
V (x) = x(t)T P x(t),
P ∈ Rn×n
(2.10)
が存在し、その時間微分が
V̇ (x) < 0
(2.11)
P A + AT P < 0
(2.12)
となることである。これを LMI で表現すると
となる P > 0 が存在することである。また、Q = P −1 > 0 とおくと (2.12) 式の双対表現として
AQ + QAT < 0
(2.13)
を得ることが出来る。
2.3.2
H∞ 性能条件
(2.9) 式のシステムが H∞ 性能を満足するための必要十分条件は以下の定理である。この定理を有界
実補題 (Bounded real lemma)、または KYP 補題 (Kalman-Yacubovich-Popov lemma) と呼ぶ。
定理 2.1 有界実補題
(2.9) 式のシステムが安定かつ、H∞ ノルムが γ 以下であるための必要十分条件は、次の LMI 条件を
満足する正定行列 P > 0 が存在することである。


AT P + XP P B C T


(2.14)
BT P
−γ 2 I DT  < 0

C
D
−I
この双対形式は、下記条件を満足する正定行列 Q > 0 が存在することである。


QAT + AQ QB T C


BQ
−γ 2 I D  < 0

CT
DT −I
(2.15)
■
2.3.3
極配置制約条件
極配置はシステムの時間応答を決定づける重要な制御仕様である。また、コントローラが不必要に、
サンプリング周期を超える程の早いダイナミクスを持つことを回避するという観点からも、非常に重要
である。そこで、ここでは設計者が閉ループ系の極を任意に配置できるようにするための前準備とし
て、極配置と LMI の関係について説明する。まず最初に LMI 領域について定義する。
定義 2.3 LMI 領域
複素平面上の集合 D は、以下の式のような対称行列 α = [αkl ] ∈ Rm×m と行列 β = [βkl ] ∈ Rm×m が
存在するときに、LMI 領域と呼ぶ。
D = {z ∈ C : fD (z) < 0}
(2.16)
fD = α + zβ + z̄β T = [αkl + βkl z + βlk z̄]1≤k,l≤m
(2.17)
ここで、
15
■
システムのすべての極が領域 D の中に含まれる時、そのシステムを D-安定 (D-stable) であるという。
また、そのときの必要十分条件として、次の定理が与えられる。
定理 2.2 D-安定条件
(2.9) 式のシステムが D-安定であるための必要十分条件は、以下の式を満足する正定行列 P > 0 が存
在することである。
MD (A, P ) < 0
(2.18)
ここで、
MD := α⊗P + β⊗(AT P ) + β T ⊗(P A) = [λkl P + µkl AT P + µlk P A]1≤k,l≤m
(2.19)
を示しており、⊗ はクロネッカー積を意味する。
■
では、定理 2.2 を使って、幾つかの代表的な極配置領域と LMI との関係を示す。アルファ安定領域は
Fig.2.3、円安定領域は Fig.2.4、コニックセクタ安定領域は Fig.2.5 である。
アルファ安定領域
2αP + P A + AT P < 0
α > 0, ∃P > 0 :
⇐⇒
円安定領域
"
r > 0, ∃P > 0 :
"
⇐⇒
(2.20)
2αQ + AQ + QAT < 0
∃Q > 0 :
∃Q > 0 :
−rP
AT P
PA
−rP
−rQ AQ
QAT −rQ
#
<0
(2.21)
#
<0
コニックセクタ安定領域
π
0 ≤ θ ≤ , ∃P > 0 :
2
"
"
⇐⇒
2.4
2.4.1
∃Q > 0 :
sin θ(P A + AT P ) cos θ(P A − AT P )
− cos θ(P A − AT P ) sin θ(P A + AT P )
QAT )
QAT )
sin θ(AQ +
cos θ(AQ −
− cos θ(AQ − QAT ) sin θ(AQ + QAT )
#
<0
(2.22)
#
<0
シンセシス LMI
状態フィードバック
次の連続時間線形時不変系の一般化プラントについて考える。
½
ẋG (t) = AG xG (t) + BG1 w(t) + BG2 u(t)
z(t) = CG1 xG (t) + DG11 w(t) + DG12 u(t)
(2.23)
ここで、AG ∈ Rn×n , DG11 ∈ Rp1 ×m1 である。このシステムに状態フィードバック
u(t) = F xG (t)
を考えると、閉ループ系は
½
ẋG (t) = Acl xG (t) + BG1 w(t)
z(t) = Ccl xG (t) + DG11 w(t)
16
(2.24)
(2.25)
となる。ただし、Acl = AG + BG2 F, Ccl = CG1 + DG12 F を表す。このシステムが H∞ 性能を満足す
るには、(2.14) 式の有界実補題を満足する正定行列 P > 0 が存在することと等価であるので、(2.25) 式
を (2.14) 式に代入すると次のようになる。

 

Acl P + P ATcl BG1 P CclT
(AG + BG2 F )P + P (AG + BG2 F )T BG1 P (CG1 + DG12 F )T

 

T
T
T
T
BG1
−γ DG11
BG1
−γI
DG11

=

T
T
Ccl P
DG11 −γI
(CG1 + DG12 F ) P
DG11
−γI


T
T + P F T DT
AG P + BG2 F P + P ATG + P F T BG2
BG1 P CG1
G12


T
T
(2.26)
=
BG1
−γI
DG11
<0
CG1 P + DG12 F P
DG11
−γI
ここで、変数は F, P の 2 つであり、F P, P F T 項などの非線形項が存在するために、これは LMI では
ない。そこで、1 つの解決策として、変数変換法と呼ばれる方法を紹介する。新しい変数 Y を
Y = F P, Y ∈ Rm×n
(2.27)
と定義して代入すると、次のような LMI を導出することが可能である。
定理 2.3 有界実補題 (状態フィードバック)
(2.23) 式のシステムが安定、かつ H∞ 性能を満足するための必要十分条件は、以下の LMI 条件を満
足する行列 P > 0 と Y が存在することである。


T
T + Y T DT
AG P + BG2 Y + P ATG + Y T BG2
BG1 P CG1
G12


T
T
(2.28)
BG1
−γI
DG11

<0
CG1 P + DG12 Y
DG11
−γI
■
同様にして極配置制約にもこのアプローチを適用する。はじめに、(2.23) 式のシステムの極を Fig.2.3
のアルファ安定領域にあることを保証するためには、
Acl P + P ATcl + 2αP
= (AG + BG1 F )P + P (AG + BG1 F )T + 2αP
T
= AG P + BG1 |{z}
F P +P ATG + P
F T} BG1
+ 2αP < 0
| {z
Y
Y
(2.29)
T
したがって、整理すると、次のように LMI 定式化される。
系 2.1 アルファ安定保証領域 (状態フィードバック)
(2.23) 式のシステムの極を Fig.2.3 の領域にあることを保証するためには、以下の LMI 条件を満足す
る行列 P > 0 と Y が存在することである。
T
AG P + BG1 Y + P ATG + Y T BG1
+ 2αP < 0
(2.30)
■
(2.23) 式のシステムの極を Fig.2.4 の円安定領域にあることを保証するためには、
"
#
"
#
−rP Acl P
−rP
(AG + BG F )P
=
P T ATcl −rP
P (AG + BG F )T
−rP


−rP
AG P + BG |{z}
FP


Y

= 
T
T
 P AT + P F B T

−rP
G
G
| {z }
YT
17
(2.31)
したがって、整理すると、次のように LMI 定式化される。
系 2.2 円安定保証領域 (状態フィードバック)
(2.23) 式のシステムの極を Fig.2.4 のコニックセクタ安定領域にあることを保証するためには、以下
の LMI 条件を満足する行列 P > 0 と Y が存在することである。
"
−rP
AG P + BG Y
T
−rP
P ATG + Y T BG
#
<0
(2.32)
■
"
(2.23) 式のシステムの極を Fig.2.5 のコニックセクタ安定領域にあることを保証するためには、
#
"
#
sin θ(Acl P + P ATcl ) Block(2,1)T
sin θ{(AG + BG F )P + P (AG + BG F )T } Block(2,1)T
=
cos θ(P ATcl − Acl P ) Block(1,1)
cos θ{P (AG + BG F )T − (AG + BG F )P } Block(1,1)


´
³
T
T
T
sin θ AG P + BG |{z}
F P +P ATG + P
F
B
Block(2,1)
| {z } G




Y
YT
³
´
= 
(2.33)
<0
 cos θ P AT + P F T B T − AG P − BG F P

Block(1,1)
G
G
| {z }
|{z}
Y
YT
したがって、整理すると、次のように LMI 定式化される。
系 2.3 コニックセクタ安定保証領域 (状態フィードバック)
(2.23) 式のシステムの極を Fig.2.5 の領域にあることを保証するためには、以下の LMI 条件を満足す
る行列 P > 0 と Y が存在することである。
´
³
T
sin θ AG P + BG Y + P ATG + Y T BG
³
´

T −A P −B Y
cos θ P ATG + Y T BG
G
G

Block(2,1)T
Block(1,1)

<0
(2.34)
■
上記で述べてきた H∞ 性能と極配置の制御仕様を同時に満足させる多目的制御問題において、LMI
ベース制御系設計法は非常に容易であり、設定したい LMI をただ複数個並べてやるだけで良い。ただ
し、正定行列 P > 0 と Y は設定した LMI において共通解でなければならない。上記の極配置制約付き
H∞ 状態フィードバック制御問題は、次のような問題の中で γ を最小化する最適化問題として定式化さ
れる。
極配置制約付き状態フィードバック制御問題
Minimize γ Find P > 0, Y subject to (2.28),(2.30),(2.32),(2.34)
■
上記の問題を解くことにより、次のようにして極配置付き状態フィードバックゲインが決定される。
F = Y P −1
2.4.2
(2.35)
出力フィードバック
次の連続時間線形時不変系の一般化プラントについて考える。
( ẋG (t) = AG xG (t) + BG1 w(t) + BG2 u(t)
z(t) = CG1 xG (t) + DG11 w(t) + DG12 u(t)
yG (t) = CG2 xG (t) + DG21 w(t) + DG22 w(t)
18
(2.36)
ここで、AG ∈ Rn×n , DG11 ∈ Rp1 ×m1 , DG22 ∈ Rp2 ×m2 であるとする。また、簡単のために DG22 = 0
と仮定する。このシステムに
u(t) = F yG (t)
(2.37)
をフィードバックするために、次のような動的な出力フィードバックコントローラを考える。
½
ẋK (t) = AK xK (t) + BK y(t)
u(t) = CK xK (t) + DK y(t)
(2.38)
ただし、AK ∈ Rk×k , DK ∈ Rm2 ×p2 である。この時、閉ループ系は次のように記述される。
#
# "
"
#"
#
"
xG (t)
BG1 + BG2 DK DG21
ẋG (t)
AG + BG2 DK CG2 BG2 CK
+
w(t)
=
BK CG2
xK (t)
xK (t)
BK CG2
AK
|
{z
}
{z
}
|
Acl
B
cl
"
#
h
i x (t)
G
+ (DG11 + DG12 DK DG21 ) w(t)
z(t)
=
CG1 + DG12 DK CG2 DG12 CK
|
{z
}
|
{z
} xK (t)
D
cl
C
cl
(2.39)
上記の記述の簡単化のために、次のような表現を導入する。
"
#
"
#
h
i
AG 0
BG1
; B0 =
A0 =
; C0 = CG1 0
;
p1×(n+k)
0 0k
0
(n+k)×(n+k)
(n+k)×m1
"
#
"
#
0 BG2
0
Ik
(2.40)
B =
; C=
;
Ik
0
CG2 0
(n+k)×(k+m2 )
(k+p2)×(k+n)
"
#
h
i
0
D12 =
; D21 =
;
0 DG12
p1×(k+m2)
DG21
(k+p2)×m1
ここで簡単化のために、(2.38) 式のコントローラも、次のように 1 つの変数 ΩK ∈ R(k+m2 )×(k+p2 ) で表
現することにする。
#
"
AK BK
ΩK =
(2.41)
CK DK
その結果、(2.39) 式の Acl ∈ R(n+k)×(n+k) , Bcl ∈ R(n+k)×m1 , Ccl ∈ Rp1×(n+k) , Dcl ∈ Rp1×m1 は、
(2.41) 式を使ったアファイン表現をすることが可能である。
Acl = A0 + BΩK C,
Bcl = B0 + BΩK D21
Ccl = C0 + D12 ΩK C, Dcl = D0 + D12 ΩK D21
(2.42)
極制約付き H∞ 出力フィードバック制御問題を解くためには、次の 2 つの不等式を満足する正定行列 P
と、コントローラ行列 ΩK を求める必要がある。
[λkl P + µkl Acl P + µlk P ATcl ]1≤k,l≤m < 0


Acl P + P ATcl Bcl P CclT

T
T <0
Bcl
−γI Dcl


Ccl P
Dcl −γI
(2.43)
(2.44)
ここで、(2.43) 式と (2.44) 式は、いずれも LMI ではなく、状態フィードバックの時の (2.27) 式の変数
変換を適用することも出来ない。そこで、1 つの解決策として次のような方法を適用する。
19
正定対称行列 P とその逆行列 P −1 を
"
R
P =
MT
M
U
#
"
, P
−1
=
S
NT
N
V
#
(2.45)
とする。ここで、R ∈ Rn×n , S ∈ Rn×n は対称行列である。P と P −1 は逆行列の関係より、P P −1 =
I, P −1 P = I が成立するので、
"
# " #
"
#
" #
R
I
R
I
P −1 P = I ⇒ P −1
=
⇒
=P
(2.46)
T
T
M
0
M
0
同様にして、
"
PP
−1
=I
⇒
(2.46) 式と (2.47) 式より、
S
NT
P
#
"
=
I
0
#
"
⇒
I
0
#
"
=P
S
NT
#
# "
#
R I
I S
P
=
MT 0
0 NT
|
{z
} |
{z
}
(2.47)
"
Π1
(2.48)
Π2
ここで、P は正定対称行列なので、次のように表すことが可能である。
−T T
P = Π2 Π−1
1 = Π 1 Π2
(2.49)
したがって、その逆行列と (2.48) 式のような関係を有する正定対称行列においては、その行列そのもの
は、R, S, M, N のみによって表現することが可能であり、U, V には依存しない。つまり、言い換え
れば、U, V も R, S, M, N によって表現されることを意味する。
まず、閉ループ系の内部漸近安定性を保証する条件は、(2.12) 式の Lyapunov 不等式より
ATcl P + P Acl < 0
(2.50)
を満足する正定対称行列 P > 0 が存在することである。そこで、(2.49) 式を (2.50) 式の中に代入すると
次のようになる。
−1
−T T
T
∃Π2 Π−1
(2.51)
1 > 0 : Acl Π2 Π1 + Π1 Π2 Acl < 0
| {z } | {z }
P
左から
ΠT1 、右から
P
Π1 を掛けると
−T T
T T
T
ΠT1 (ATcl Π2 Π−1
1 + Π1 Π2 Acl )Π1 = Π1 Acl Π2 + Π2 Acl Π1
(2.52)
となるので、次の LMI と等価である。
T T
T
T
∃Π2 Π−1
1 > 0 : (Π2 Acl Π1 ) + Π2 Acl Π1 < 0
ここで、
"
ΠT2 Acl Π1
=
AG R + BCK
AK
A + BDK C
SA + BK C
(2.53)
#
(2.54)
である。ただし、
AK
BK
CK
DK
=
=
=
=
N AK M T + N BK CR + SBCK M T + S(A + BDK C)R
N BK + SBDK
CK M T + DK CR
DK
(2.55)
である。結果として、出力フィードバック問題における Lyapunov 安定定理は次のように定式化される。
20
定理 2.4 Lyapunov 安定定理 (出力フィードバック)
(2.36) 式のシステムが Lyapunov 安定性を満足するための必要十分条件は、次の LMI 条件を満足す
る R, S, AK , BK , CK , DK が存在することである。
1.
"
RAT AR + CTK B T + BCK
AK + AT + C T DTK B T
2.
"
ATK + A + BDK C
AT S + SA + C T BTK + BK C
R I
I S
#
<0
(2.56)
#
>0
(2.57)
■
R, S, AK , BK , CK , DK が求まれば、対応する出力フィードバンクは
I − SR = N M T
(2.58)
を満足する N, M を求め
AK
BK
CK
DK
=
=
=
=
N −1 (AK − S(A + BDK C)R − N BK CR − SBCK M T )M −T
N −1 (BK − SBDK )
(CK − DK CR)M −T
DK
(2.59)
のようにして、出力フィードバックコントローラ AK , BK , CK , DK を構成することが可能である。
次に、(2.36) 式のシステムが H∞ 性能条件を満足する、つまり、
sup
w∈L2
kzk2
<γ
kwk2
(2.60)
を満たすための条件は、次の有界実補題として表される。
定理 2.5 有界実補題 (出力フィードバック)
(2.36) 式のシステムが、H∞ 性能条件を満足するための必要十分条件は、以下の LMI 条件を満足す
る R, S, AK , BK , CK , DK が存在することである。
1.





RAT + AR + CTK B2T + B2 CK
T BT
AK + AT + C2T DK
2
T D BT
B1T + D21
K 2
C1 R + D12 CK
B1 + B2 DK D21
SB1 + BK D21
−γ 2 I
D11 + D12 DK D21
2.
ATK + A + B2 DK C2
AT S + SA + C2T BTK + BK C2
T BT
B1T S + D21
K
C1 + D12 DK C2

T
RC1T + CTK D12
T DT 
C1T + C2T DK
12 
<0
T
T
T
T 

D11 + D21 DK D12
−I
"
R I
I S
(2.61)
#
>0
(2.62)
■
極配置も同様に、出力フィードバック問題における D-安定条件が次のように定式化される。
21
定理 2.6 D-安定条件 (出力フィードバック)
(2.36) 式のシステムが D-安定であるための必要十分条件は、以下の LMI 条件を満足する R, S, AK ,
BK , CK , DK が存在することである。
1.
"
Φ⊗
R I
I S
#
"
+ Θ⊗
+
#
AR + BCK A + BDK C
AK
SA + BK C
"
#T
AR
+
BC
A
+
BD
C
K
K
ΘT ⊗
<0
AK
SA + BK C
2.
"
(2.63)
#
R I
I S
>0
(2.64)
ここで、⊗ はクロネッカー積を意味する。
■
上記の定理 2.6 を具体的な極配置領域に適用すると次のようになる。
系 2.4 アルファ安定保証領域 (出力フィードバック)
(2.36) 式のシステムの極を Fig.2.3 の領域にあることを保証するためには、次の LMI 条件
1.
"
2α
R I
I S
#
"
+
RAT + AR + CTK B T + BCK
T BT
AK + AT + C T DK
2.
"
R I
I S
AK + A + BDK C
T
A S + SA + C T BTK + BK C
#
<0
(2.65)
#
>0
(2.66)
を満足する R, S, AK , BK , CK , DK が存在すれば良い。
■
系 2.5 円安定保証領域 (出力フィードバック)
(2.36) 式のシステムの極を Fig.2.4 の領域にあることを保証するためには、次の LMI 条件
1.

"
−r


 "
 AR + BC
K

AK
2.
R I
I S
#
"
A + BDK C
SA + BK C
"
#
AR + BCK
AK
"
−r
R I
I S
A + BDK C
SA + BK C
#
R I
I S
# 


<0


(2.67)
#
>0
(2.68)
を満足する R, S, AK , BK , CK , DK が存在すれば良い。
■
22
系 2.6 コニックセクタ安定保証領域 (出力フィードバック)
(2.36) 式のシステムの極を Fig.2.5 の領域にあることを保証するためには、次の LMI 条件
1.
"
# "
#
sin θ cos θ
AR + BCK A + BDK C
⊗
+
cos θ sin θ
AK
SA + BK C
"
# "
#
sin θ − cos θ
RAT + CTK B T
ATK
⊗
<0
T B T AT S + C T B T
cos θ sin θ
AT + C T DK
K
2.
"
R I
I S
(2.69)
#
>0
(2.70)
を満足する R, S, AK , BK , CK , DK が存在すれば良い。
■
状態フィードバックと同様にして、極配置制約付き H∞ 出力フィードバック制御問題は、次のような
問題の中で γ を最小化する最適化問題として定式化される。
極配置制約付き出力フィードバック制御問題
Minimize γ Find R, S, AK , BK , CK , DK subject to (2.61), (2.62), (2.65), (2.67), (2.69)
■
上記の問題を解くことにより、(2.58), (2.59) 式を用いることで出力フィードバックコントローラ ΩK =
(AK , BK , CK , DK ) を構成することが可能である。
2.5
ゲインスケジューリング制御の考え方
スケジューリング変数に応じて動特性が変化する物理対象に対して、その変数の変化速度が安定領域
に強く影響を及ぼすことを理解するために、次のような数値例を考えてみる [13][14]。
"
#
−1 + αθ12 1 − αθ1 θ2
ẋ(t) =
x(t)
−1 − αθ1 θ2 −1 + αθ22
|
{z
}
(2.71)
A(θ)
ただし、パラメータ θi , (i = 1, 2) は Fig.2.8 に示すような θ12 + θ22 ≤1 範囲にあるものとする。
最初に簡単化のために、パラメータ θi , (i = 1, 2) を固定して線形時不変システムとして考える。その
ときのシステムの安定性は、行列 A(θ) の固有値を調べることにより確認することができるので、
det(sI − A(θ)) = s2 + (2 − α(θ12 + θ22 ))s + 2 − (θ12 + θ22 )
(2.72)
この行列式の解が負であるためには、α < 2 を満足する必要がある。つまり、線形時不変システムが安
定であるための必要十分条件は α < 2 である。しかし、Fig.2.8 の制約範囲内において、パラメータを
θ1 (t) = cos ωt, θ2 (t) = sin ωt と変動させ線形時変システムとして考えると、このときの解は次のよう
に与えられることが分かっている。
"
#
"
#
cos ωt sin ωt
−1 + α 1 − ω
x(t) = Φ(t, 0)x(0), Φ(t, 0) =
exp(Bt), B =
(2.73)
− sin ωt cos ωt
−1 + ω
−1
23
このシステムが安定であるための必要十分条件は、B の固有値の実部が負であるので、
α < 2,
1 − α + (1 − ω)2 > 0
(2.74)
である。Fig.2.9 にパラメータが変動する際の安定領域を示す。パラメータが固定されたときと変動す
るときとでは、安定領域が異なることはもちろん、その変化スピードによっても安定領域が変化するこ
とが分かる。つまり、パラメータを固定したシステムに対して設計したコントローラが得られたとして
も、そのコントローラをパラメータが変化するシステムに対して適応することは出来ない。また、パラ
メータの変化速度によっても安定領域は異なるために、Fig.2.7 に示すように、変化速度に応じたコン
トローラ設計が必要である。
近年、盛んに研究されてきたロバスト制御は、応用面においても十分な成果をあげている。ここでは、
そのロバスト制御と GS 制御の違いを整理しておく。ロバスト制御は、制御対象の持つ不確かさや非線
形性に対しても十分な制御性能が実現される制御手法である。つまり、Fig.2.10 に示すように、1 つの
固定コントローラにより制御対象の変化に対応する。また、対象が複雑になると、Fig.2.11 に示すよう
に、予め制御対象の変動を記述しておき、それに対して同時安定化を図る様な 1 つの固定コントローラ
を設計する。しかしながら、いずれにおいても、大きな、そして急速な制御対象の変動効果に対応する
ことには限界があるだけでなく、Fig.2.11 の方法についてはどうしても保守的なコントローラを設計す
ることになる。一方で GS 制御は、Fig.2.12 に示すように、対象の大きな変動や非線形性はコントロー
ラのスケジューリングで対応して、細かい小さな変動や不確かさはスケジューリング変数の変化速度を
考慮しつつローカルにロバスト設計でカバーすることにより、良い制御性能を実現することが可能であ
る。また、この方法は、基本的には線形制御理論の考え方を局所的に採用していくために、成熟した理
論や解析ツールを援用することができる。相手の変化に応じて、自分自身が変化していくという考え方
は、学習制御などによくある考え方であるが、GS 制御は生物の適応メカニズムほど柔軟な役割を果す
ことは出来ないが、失敗による繰り返し学習が許されない状況においては、理論的に性能保証を与えな
がら運動制御が出来るという意味で大きな意味を持つ [14]。
現在、GS 制御系設計法には2つの代表的な方法がある。1つはコントローラ補間法、もう1つが線
形パラメータ変動法である。ここでは、簡単に、2 つの方法について整理しておく。
コントローラ補間法
θ の属する領域を Θ とする。θ でパラメトライズされた線形システムの集合 G(θ) が制御対象のモデ
リングによって得られたとする。これに続く設計のプロセスとして、Θ 上の離散点 θ 毎に対してコント
ローラを設計し、それらを補間していく方法がある。それをコントローラ補間法という。この方法を整
理すると、次のようになる。
1. 離散的な点 pi ∈ Θ, i = 1, . . . , N 上の線形システム G(pi ) に対して、線形制御理論によりコント
ローラ K(pi ) を設計する。
2. K(pi ) を Θ 上で補間し、GS 制御則を構成する。
この方法は、実用面において実績を上げてきた方法である。しかし、離散点での制御対象の線形モデル
G(pi ) の情報しか必要としないが、系統的に設計を行うには、どのような pi を選択するか、その pi で
どの程度厳しい仕様で局所的な線形制御系設計を行うかを事前に決定することは難しい。結果的にどれ
くらいの変化速度まで全体の GS 制御系が許容できるかは、コントローラを定めた後に解析される。
線形パラメータ変動法
もう一つの方法が LPV 法である。この方法は、θ に連続的に依存した線形パラメータ変動モデル G(θ)
を用いて、コントローラも θ に連続的に直接設計してしまう方法である。LPV 法の大まかな特徴を述
べると以下のようになる。
24
1. Θ の離散的な点ではなく、すべての点に対してコントローラ設計
2. 制御対象 ( 近似でない ) の LPV モデル表現を導くことはいつでも可能というわけではなく、個別
の工夫が必要
3. θ の変化速度に依存することなく、時変系としての安定性を理論的に保証
この方法は近年盛んに研究されている方法であり、理論面においても十分な成果が挙げられている。本
章では、以後、理論的に体系化された LPV 法に着目して、詳細に説明を加えていく。
2.6
LPV システムの構造上の特徴と性能解析
ここで LPV 法について述べる前に、LPV システムについて説明する。LPV システムとは、線形時変
システム (Linear Time Varying, 以後, LTV と略記) のある特種なクラスであり、システム変動を表すパ
ラメータ θ(t) によって係数行列が決定される線形システムである。(2.75) 式にそれを示す。LPV、LTV
と線形時不変システム (Linear Time Invariant, 以後, LTI と略記) の 3 つのシステムの関係を Fig.2.13
に示す [8]。
½
ẋ(t) = A(θ(t))x(t) + B(θ(t))u(t)
(2.75)
y(t) = C(θ(t))x(t) + D(θ(t))u(t)
θ(t) はこのシステムの係数行列を定める変数である。また簡単のため、ここでは、A(θ), B(θ), C(θ), D(θ),
は有界で θ(t) に対して連続な行列関数とする。3 つのシステムの間には次のような関係がある。LPV シ
ステムの変動パラメータ θ(t) の軌跡を特定なパラメータ軌跡に固定する、つまり、θ(t) := θ∗ (t) とする
と、そのシステムは LTV システムとなる。さらに、そのパラメータ軌跡の時間をある時間 t0 に固定す
る、つまり、t = t0 とすると、LTV システムは LTI システムとなる。一方で、LPV システムの変動パ
ラメータを固定すると、つまり、θ(t) := θ0 とすると、LPV システムは LTI システムとなるという関
係を持つ。大きな特徴として、LTI、LTV システムはオフラインでも十分にシステムの特徴を把握でき
るのに対して、LPV システムはオンラインで変動パラメータを計測可能なときにのみ、システムの特
徴を把握することが可能である。
LPV システムを考えるにおいて、θ(t), θ̇(t) は適当な有界領域 Θ, ν に含まれているとすることが多
い。特に、θ の各要素 θk に対して、
"
# "
"
#
#
k0
X
A(θ) B(θ)
A0 B0
Ak Bk
θk
(2.76)
=
+
C(θ) D(θ)
C0 D0
C
D
k
k
k=1
の形で表されるとき、アファイン型 LPV システムと呼ばれ、モデリングを考える上で扱いやすい表現
である。また、このタイプのモデルは、解析、設計の両面から非常に取り扱い易いために、一般的には、
アファイン型 LPV システムの θ の範囲が有界な多角形のときポリトピック型に、ノルムで押さえられ
ている時には線形分数変換 (Linear Fractional Transformantion, 以後, LFT と略記) 型に変換して取り
扱うことが多い。これらについて簡単に説明する。
2.6.1
構造上の特徴
ポリトピック型 LPV システム
パラメータ θ(t) が次のような端点のポリトープ集合 Θ として表されるとする。
Θ = {θvi : i = 1, 2, · · · , N }
25
(2.77)
Θ に属する任意の時間のパラメータ θ(t) は、端点のパラメータ値を使って次のように分解することが
可能である。
N
N
X
X
ρ(t) =
αi (t)ρvi ,
αi (t) = 1, αi ≥ 0
(2.78)
i=1
i=1
(2.76) 式で表されるアファイン型 LPV システムは、次のようなポリトピック型 LPV システムに変換す
ることが可能である。Fig.2.14 にポリトピック型 LPV システムの構造を示す。
#
"
#
"
N
X
Ai Bi
A(θ(t)) B(θ(t))
αi
(2.79)
=
Ci Di
C(θ(t)) D(θ(t))
i=1
ここで、Ai , Bi , Ci , Di は各端点における値を意味する。
LFT 型 LPV システム
(2.76) 式で表されるアファイン型 LPV システムは、次のような LFT 型 LPV システムに変換するこ
とも可能である。
"
# "
# "
#
h
i
A(θ(t)) B(θ(t))
A B
Bθ
∆(I − Dθθ ∆)−1 Cθ Dθ
(2.80)
=
+
Dθ
C(θ(t)) D(θ(t))
C D
ここで、∆ はブロック対角構造をしたパラメータ行列で
∆ = diag{ρI, · · · , ρL I}
(2.81)
という構造をしている。Fig.2.15 に LFT 型 LPV システムの構造を示す。
2.6.2
性能解析
LPV システムの安定性や性能の解析は、LTV システム
½
ẋ(t) = A(t)x(t) + B(t)u(t)
y(t) = C(t)x(t) + D(t)u(t)
(2.82)
に対する結果が基本となる。ここで、A(t), B(t), C(t), D(t) は有界連続な時間関数であると仮定する。
また、LPV システムの性質には、θ の変化速度 θ̇ も関係するので、θ̇ の含まれる集合を ν と定義する。
ν としては各成分の θi の変化速度の上下限 ±νi を用いて、
ν = {θ̇| |θ̇i | ≤ νi }
(2.83)
と考える。LTV システムの安定性について、次を示すことが可能である。
定理 2.7 LTV システムの安定性
ここで、(2.82) 式の LTV システムにおいて u(t) = 0 の場合を考える。この時、有界で連続微分可能
な時間関数 P (t) が存在して
P (t) > 0, Ṗ (t) + A(t)T P (t) + P (t)A(t) < 0
が成り立てば、このシステムは大域的漸近安定である。
26
(2.84)
■
(2.84) 式は、Lyapunov 方程式の Lyapunov 関数 P (t) が時間の関数となっているため、Ṗ (t) が残ったも
のに対応する。この式は、次のように書き換えることが出来る。
定理 2.8 LPV システムの安定性
(2.75) 式の LPV システムにおいて u(t) = 0 の場合を考える。この時、有界で連続微分可能な P (θ) が
存在して、次の線形微分行列不等式
P (θ) > 0, r
X
θ̇i
i=1
∂P (θ)
+ AT (θ)P (θ) + P (θ)A(θ) < 0, ∀θ ∈ Θ, ∀θ̇ ∈ ν
∂θi
(2.85)
が成り立てば、このシステムは大域的漸近安定である。
■
ここで、(2.85) 式は Ṗ (θ) = θ̇(∂P /∂θ) であることから (2.84) 式と等価である。(2.85) を満たす定数行
列 P が存在するならば、微分項が消えて
P > 0, AT (θ)P + P A(θ) < 0, ∀θ ∈ Θ
(2.86)
となる。このときのシステムは二次安定である。二次安定の場合は、スケジューリング変数の任意の変
化速度 θ̇ に対して、漸近安定となるが、実際には物理的な理由から θ̇ に上下限値があることが多いので、
不必要に速い変化速度に対する安定性を保証することを意味する。
LPV システムは線形の時変システムであるので極の概念は無いが、次のように Lyapunov 関数を通
して状態変数の減衰率の上限を知ることができる。
定理 2.9 LPV システムの減衰率
(2.75) 式の LPV システムについて考える。正の実数 α, β, γ に対して、有界連続微分可能な P (θ) が
存在して
αI ≤ P (θ) ≤ βI,
r
X
i=1
θ̇i
∂P (θ)
+ AT (θ)P (θ) + P (θ)A(θ) < −γI, ∀θ ∈ Θ, ∀θ̇ ∈ ν
∂θi
が成立するとき、(2.82) で u(t) = 0、初期値 x(t0 ) = x0 に対して解 x(t) は次式を満たす。
³ β ´1/2
¢
γ ¡
t − t0 }, t0 ≤t
exp{−
kx(t)k ≤ kx0 k
α
2α
(2.87)
(2.88)
■
システムの入力信号 u(t)、出力信号 y(t) に対して、
sup
u6=0
kyk2
kuk2
を L2 ゲインと呼ぶ。ただし、k·k2 は次のように定義する。
µZ ∞
¶1/2
T
kxk2 =
x (t)x(t)dt
(2.89)
(2.90)
0
LTI システムに対して、L2 ゲインは H∞ ノルムと一致して、ロバスト制御で大きな役割を果す。(2.82)
式の LTV システムに対して、次の定理を言うことが可能である。
27
定理 2.10 LTV システムの L2 ゲイン性能
(2.82) 式の LTV システムについて考える。有界連続微分可能な時間関数 P (t) が存在して、


Ṗ (t) + AT (t)P (t) + P (t)A(t) P (t)B(t) C T (t)


P (t) > 0, 
B T (t)P (t)
−γI
DT (t)  < 0, ∀θ ∈ Θ, ∀θ̇ ∈ ν
C(t)
D(t)
(2.91)
−γI
が成立すれば、このシステムは γ 未満の L2 ゲイン性能を持つ。
■
これは、(2.14), (2.15) 式の LTI システムの有界実補題に Ṗ (t) が付加されただけのものである。これを
用いて、安定性の場合と同様にして、以下の結果を得ることが出来る。
定理 2.11 LPV システムの L2 ゲイン性能
(2.75) 式の LPV システムについて考える。この時、有界連続微分可能な P (θ) が存在して、次の線形
微分行列不等式
 X

r
∂P (θ)
T
T (θ)
+
A
(θ)P
(θ)
+
P
(θ)A(θ)
P
(θ)B(θ)
C
θ̇
i


∂θi


P (θ) > 0,  i=1
 < 0, ∀θ ∈ Θ, ∀θ̇ ∈ ν
T
T

B (θ)P (θ)
−γI
D (θ) 
C(θ)
D(θ)
−γI
(2.92)
が成立すれば、このシステムは γ 未満の L2 ゲイン性能を持つ。
■
L2 ゲイン性能についても二次安定性と同様に、条件としては保守的になるものの、設計と解析の簡単
さから X(θ) を定数行列とする方法がしばしば用いられる。
2.7
LPV モデルに基づくゲインスケジューリング制御
ここでは、LPV モデルに基づくゲインスケジューリング制御 (以後、LPV 制御) 系設計について説明
する。また、以後に述べる設計法は、計算量緩和のために定数解を用いた方法である。まずはじめに、
次のような LPV システムの一般化プラントを考える。
( ẋG (t) = AG (θ)xG (t) + BG1 (θ)w(t) + BG2 (θ)u(t)
z(t) = CG1 (θ)xG (t) + DG11 (θ)w(t) + DG12 (θ)u(t)
yG (t) = CG2 (θ)xG (t) + DG21 (θ)w(t) + DG22 (θ)w(t)
(2.93)
ここで、AG ∈ Rn×n , DG11 ∈ Rp1 ×m1 , DG22 ∈ Rp2 ×m2 であり、それぞれの係数行列は、有界連続な
θ の関数である。また、次のような仮定が成立するものとする。
1. スケジューリング変数は、θ ∈ Θ と θ̇ ∈ ν の条件を満足し、その中で (AG (θ), BG2 (θ)) が可安定、
(CG2 (θ), AG (θ)) は可検出である。
2. DG22 (θ) = 0
2 番目の条件は、簡単化のために導入したものであり、決して一般性を失うことはない。この LPV シ
ステムに対して、
u = F (θ)yG (t)
(2.94)
28
をフィードバックするために、次のような動的な出力フィードバックコントローラを考える。
½
ẋK (t) = AK (θ)xK (t) + BK (θ)y(t)
u(t) = CK (θ)xK (t) + DK (θ)y(t)
(2.95)
コントローラの次元は対象のプラントと同じ次元であるとする。LPV 制御のポイントは、LPV システ
ムの変化に応じてコントローラをスケジュールさせる点にある。そのためには、予め、対象に応じたコ
ントローラを設計しておく必要がある。また、その際、スケジューリング変数の変化率による安定性の
領域の変化にもうまく対応できなければならない。この要求を満足するためには、各端点において構成
した LMI を同時に解くことにより解決することできる。つまり、極配置制約付き出力フィードバック
LPV 制御問題は、次のように定式化される。
極配置制約付き出力フィードバック LPV 制御問題
(i)
(i)
(i)
(i)
Minimize γ Find R, S, AK , BK , CK , DK
subject to (2.61), (2.62), (2.65), (2.67), (2.69), (i = 1, · · · , N )
■
ここで、(2.61), (2.62), (2.65), (2.67), (2.69) 式は各端点における LMI である。また、上付添字 i は端
点におけるコントローラであり、N は端点の数を意味する。この端点ごとのコントローラを実装するた
めには、次のように補間すれば良い。
"
#
"
#
N
N
(i)
(i)
X
X
AK (θ) BK (θ)
AK BK
=
pi (θ(t))
,
pi (θ(t)) = 1
(2.96)
(i)
(i)
CK (θ) DK (θ)
CK DK
i=1
i=1
ただし、pi (θ(t)) はスケジューリング変数や端点の数などに応じて変化するが、例えば、スケジューリ
ング変数が 1 個、端点が 2 つの場合は、LPV コントローラは次のようになる。
#
"
"
"
#
#
(1)
(1)
(2)
(2)
θ(t) − θ2 AK BK
θ1 − θ(t) AK BK
AK (θ) BK (θ)
=
+
(2.97)
(1)
(1)
(2)
(2)
θ1 − θ2
θ1 − θ2
CK (θ) DK (θ)
CK DK
CK DK
もちろんではあるが、次の等式が成立する。
θ(t) − θ2 θ1 − θ(t)
+
=1
θ1 − θ2
θ1 − θ2
(2.98)
上記と同様の考え方の下で、極配置制約付き状態フィードバック LPV 制御問題も考えることが可能
である。参考までに、船舶の回頭角制御に状態フィードバック LPV 制御を適用した例を付録に示す。
2.8
まとめ
本章では、海洋工学分野の問題に対して LMI ベース LPV 制御を応用するために、その周辺の基礎知
識について整理した。以後の章で、具体的にこの方法論を応用することにする。
29
x^
x
x^
x
Fig. 2.1: Convex set and Nonconvex set (Left side: Convex set, Right side: Nonconvex set)
φ (x)
φ (x)
x
x^
x
x^
Fig. 2.2: Convex function and Nonconvex function (Left side: Convex function, Right side: Nonconvex function)
30
Imag
α
Real
Fig. 2.3: α-stability region
Imag
r
-r
r
Real
-r
Fig. 2.4: Circle-stability region
Imag
θ
Real
Fig. 2.5: Conic-sector-stability region
31
Imag
Imag
Imag
r
α
θ
-r
r
Real
Real
Real
-r
Imag
r
r
α
θ
Real
-r
Fig. 2.6: Extended LMI Region from single LMI region
w
z
G (θ)
u
θ
y
K (θ)
Fig. 2.7: Gain-scheduled control system structure
32
θ2
1
-1
1
θ1
-1
Fig. 2.8: Constrainted parameter region
α
2
α = (ω - 1) + 1
unstable region
2
stable region
ω
Fig. 2.9: Stable region of parameter varying system
33
Plant
Controller
Fig. 2.10: Single model based robust control
Plant 1
Contr.1
Plant 2
Fixed Controller
Contr.2
Fig. 2.11: Multi model based robust control
Plant 1
Plant (θ)
Scheduled
Contr.1
Plant 2
Scheduled
Contr.(θ)
Contr.2
Fig. 2.12: Gain-scheduled control
34
LPV
freeze the parameter
select a trajectory
θ (t):= θ 0
LTI
θ (t):= θ*(t)
freeze the time
t=t0
LTV
Fig. 2.13: Class of linear systems
θ(t)
G (s)
y
u
Fig. 2.14: Polytopic type LPV system
∆
w
z
y
G (s)
Fig. 2.15: LFT type LPV system
35
u
第3章
固定ピッチ型水平軸風車の可変速制御による高効率化
本章では、近年、クリーンエネルギーとして期待されつつある風力発電の高効率化に向けた可変速制
御について検討する。本章の構成は次の通りである。はじめに、問題背景について説明した後、本章で
対象とする固定ピッチ型水平軸風車の概要とその運転モードについて述べる。次に、その風力発電機の
数学モデルについて説明し、スケジューリング変数に依存した状態方程式を導出した後、各運転モード
に対する LPV モデルを導出し、LPV コントローラを設計する。そして、性能検証のために数値シミュ
レーションを用いて有効性を確認する。最後に、風車シミュレータという風力発電機開発支援装置を提
案し、その装置を利用して LPV 制御の検証実験の結果を示したのち、本論文のまとめと今後の課題に
ついて述べる [15][16][17][18]。
3.1
はじめに
近年、クリーンなエネルギーの一つとして風力発電が世界的に注目されている。中でも、デンマーク、
ドイツなどでは、すでに、電力需要の 10% 以上を風力発電により賄っており、EU 諸国全体では 2010
年に向けて電力需要の 22% を供給するように目標を掲げている。欧米では 2010 ∼ 2020 年にかけて、
電力の 5 ∼ 10% を風力発電により賄う時代が到来するとも言われている [19]。このような大きな電力需
要を目標に掲げることができるようになった主要因は、風車の大型化、空力性能の向上と、1990 年代
初期から研究開発が盛んに行われた制御技術の向上にある [20]。また最近では、風力発電設備の設置場
所が陸上から洋上へとシフトする傾向にある。洋上は陸上よりも一定風速以上の強い風が期待できるた
め、安定した電力供給の実現が可能であり、そのため、わずかな発電効率向上も、長期的には大きな効
果が期待できる。したがって、風力発電の制御技術の研究がますます重要になってきている。
制御技術の観点から風力発電の効率改善のためには、2 つの問題を解決する必要がある。第 1 に、定
格風速以下の風速領域において、風エネルギーから最大のパワーを獲得すること。第 2 に、定格風速以
上の風速領域においては、風のエネルギーの流入を調節しながら定格出力に維持することである。これ
らの問題を解決する方法として、ロータの回転数を風速に応じた最適値に制御する方法がある。風力発
電機の最適回転数制御の問題に関して、R.Datta と V.T.Ranganathan [21] は、従来から採用されてき
た方法について明らかにしている。その方法は、風力発電機の負荷トルク Tdtarget を
½
Tdtarget =
µ
¶ ¾
1 3
1
5
· ω2
ρπr CPmax
2
λopt
(3.1)
とするものである。ここで、ω はロータの回転速度、ρ は空気密度、r はロータの半径、CPmax は出力
係数の最大値を表し、λopt はその時の最適周速比の値を表す。この方法は、発電機の回転数の情報のみ
で必要な発電機トルクが決定され、非常に簡便な方法である。しかし、(3.1) 式は、風力発電機の静的
な釣合いの関係式から導かれたものであるために、風力発電機のダイナミクスが考慮されていない。さ
らに、空気密度は設置場所の気候等に大きく左右される事を考慮すると、この方法は高い発電効率を期
待することができない。
この問題を解決するため、涌井らは PID 制御を用いて、その制御方式の有用性について明らかにして
いる [22]。しかし、空力トルクが強い非線形性を有するため、彼らの方法は広い風速領域に適応するこ
36
とは困難である。そこで、W.E.Leithead らは、その非線形性に着目し、スケジューリング制御を適用し
ている [23]。同様にして、H.D.Battista らは、ロータ-シャフト-発電機システムに対して、Hamiltonian
から Controlled Lyapunov 関数を構成して非線形制御則を導出している [24]。しかし、以上の文献にお
いては、定格風速以上の風速領域における制御問題については言及されていない。
そこで本章では、風力発電機の数学モデルが風速変動に依存する非線形性を有することに着目し、LPV
制御方式を、最大発電効率を獲得するための回転数制御、及び定格出力を維持するための回転数制御の
2つの制御問題に適用し、数値シミュレーションと風車シミュレータと呼ばれる風力発電機開発支援装
置を利用した検証実験により、その有効性と問題点を明らかにすることを目的とする。
3.2
風力発電機の運転モード
一般的な風力発電機の運転モードは、Fig.3.1 に示すように分類される [19]。以下に、それぞれの運転
モードについて簡単に説明する。
モード 1. V < Vci
風速 V がカットイン風速 Vci 以下の領域である。風車が停止状態から起動する場合には、機械的
な摩擦などの抵抗よりも大きなトルクが必要となる。しかし、風速 V がカットイン風速 Vci 以下
の場合にはロータが始動トルクを発生することができないので、風は吹いているのに回転しない
という状況が起こる。したがって、この風速領域では発電を行うことができない。
モード 2. Vci < V < Vrated
風速 V がカットイン風速 Vci から定格風速 Vrated までの領域である。この領域では、風からの獲
得エネルギーを最大にするため、風速 V に応じた最適値にロータの回転数を制御することによ
り、効率的な発電が行われる。
モード 3. Vrated < V < Vco
風速 V が定格風速 Vrated 以上、カットアウト風速 Vco 以下の領域である。定格風速 Vrated 以上
の風速領域において、エネルギーを最大効率で獲得することは可能であるが、定格出力 Prated を
超えるような電力が獲得されたとしても、発電機が過負荷状態となってしまう。よって、定格風
速 Vrated からカットアウト風速 Vco までの領域では、ロータの回転数を制御することにより、定
格出力 Prated に維持することが必要になる。
モード 4. Vco ≤ V
風速 V が大きくなりカットアウト風速 Vco を超える領域である。この領域では、回転数や発電機
の発生トルクが大きくなり過ぎて、発電システム全体が危険になるために、強制的に回転を停止
させる必要がある。
本章では具体的な風力発電機として、Fig.3.2 に示す九州大学応用力学研究所で開発された風レンズ
風車と呼ばれる鍔付きデフューザ型の小型風車を取り上げる [25]。風レンズ風車の主要目を Table 3.1
に示す。また、本章では特に、モード 2 とモード 3 の 2 つの運転モードにおける風力発電機の回転数制
御について検討することとする。この風レンズ風車の洋上大型風車への展開は、風レンズ風車の持つデ
フューザの大型化が、現在の技術では不可能であるため、まだ考えられてはないが、本章は、風レンズ
風車のみに焦点を絞った研究ではないために、以後、述べることは、洋上で扱われるような大型風力発
電機にも適用可能である。
37
3.3
3.3.1
風力発電機の状態方程式の導出
運動の数学モデル
風力発電機のモデルは Fig.3.3 に示すように、ロータ等の機械部と直流発電機の電気部から構成され
る。本章では、シャフトや増速機などの伝達系は考えずに直結とし、また軸の捩れやシステム全体の機
械的なロス、電気的なロスはないとする。風力発電機に付加された電流コントローラが指令電流に応
じて抵抗値を速やかに調節することにより、電流が変化を起こし、ロータの回転数を制御する仕組みに
なっている。
Fig.3.3 に示した風力発電機の運動方程式は、
(
JT ω̇ = QA − KT i
(3.2)
Li̇ + Ri + e = KE ω
で与えられる。ここで、JT は慣性モーメント、L はインダクタンス、R は抵抗、KE は逆起電力定数、
KT はトルク定数を表す。また、e と i はそれぞれ、バッテリの電圧、電流、ω はロータの回転数、QA
はロータに発生する空力トルクで
1
QA = ρπr3 CT (λ)V 2
(3.3)
2
と表される。ここで V は風速、CT はトルク係数で、
CT =
QA
1
3 2
2 ρπr V
(3.4)
と定義される無次元値である。これは、風エネルギーから得られるトルクの割合を表す。同様にして、
出力係数 CP は、
QA ω
CP = 1
(3.5)
2 3
2 ρπr V
と表され、風エネルギーから取り出すことのできるパワーの割合を表す。出力係数 CP とトルク係数 CT
の間には、CP = λCT という関係があり、これらの係数はすべて、周速比
λ=
rω
V
(3.6)
を変数とした Fig.3.4 のような非線形関数として表される。ここで、出力係数 CP の最大値を CP max 、
その時の周速比 λ を λopt とすると、風力発電機は (λopt , CP max ) で最大発電可能となる。本章で対象と
する風レンズ風車は、(λopt , CP max ) = (4.51, 1.37) である。
Table 3.1: Specifications of Wind Lens
Rotor diameter
Blade number
Rated power
Cut-in wind speed
Rated wind speed
Cut-out wind speed
Rated rotational speed
Rated torque
r
N
Prated
Vci
Vrated
Vco
ωrated
Trated
600 mm
3
300 W
1 m/s
11 m/s
22 m/s
1580 rpm
1.91 N m
一般に、電気的時定数は機械的時定数よりも小さいので、(3.2) 式の第 2 式で、i̇ = 0 として近似的に
1
1
JT ω̇ = ρπr3 CT (λ)V 2 − KT (KE ω − e)
2
R
38
(3.7)
と表すことが可能である。したがって、抵抗を調整する制御問題と考えなければならないが、前述した
ように、本研究で対象とする風レンズは電流コントローラをもち、瞬時に指令電流を実現することが可
能であるので、本研究では、
1
JT ω̇ = ρπr3 CT (λ)V 2 − KT i
(3.8)
2
を風力発電機の数学モデルとして取り扱うことにする。
3.3.2
動作点周りでの線形化
はじめに、(3.8) 式の非線形方程式を動作点周りで一次近似することにより線形状態方程式を導出す
るため、動作点について確認する。動作点は静的な釣り合い状態であるので、
1
0 = ρπr3 CT (λ∗ )V 2 − KT i∗
2
(3.9)
となる。右肩添字の ∗ は動作点における値を意味する。よって、動作点における制御入力 i∗ は
i∗ =
1
ρπr3 CT (λ∗ )V 2
2KT
(3.10)
と表される。また、動作点における周速比 λ∗ を用いて、風速 V に応じた回転数 ω ∗ は (3.6) 式より、
ω∗ =
V ∗
λ
r
(3.11)
と表される。
次に (3.8) 式は、風速が一定の場合には ω̇ = f (ω, i) という関係で表されるので、これを (3.11) 式に示
した動作点周りで Taylor 展開し整理すると、次のような線形状態方程式を導くことができる。
¯
¯
∂f ¯¯
∂f ¯¯
d
∗
∗
(ω − ω ) +
(i − i∗ )
(3.12)
(ω − ω ) =
dt
∂ω ¯∗
∂i ¯∗
ただし、
f (ω, i) =
であり、
1
KT
ρπr3 CT (λ)V 2 −
i
2JT
JT
¯
¯
¯
∂f ¯¯
∂f ∂λ ¯¯
ρπr4 dCT ¯¯
=
=
,
V
∂ω ¯∗ ∂λ ∂ω ¯∗
2JT
dλ ¯∗
¯
∂f ¯¯
KT
=−
∂i ¯∗
JT
とする。ここで、導出された線形状態方程式 (3.12) 式 は、係数が風速 V に依存しているために、風
速 V に応じてシステムが変動することが分かる。
3.4
最大発電効率を得るための最適回転数制御
モード 2 における制御目的は、カットイン風速 Vci から定格風速 Vrated の間の領域において、風エネ
ルギーの獲得を最大にすることである。よって、Fig.3.5 に示すように、動作点としては風力発電機の
出力係数 CP が最大となる点を動作点 (λ∗ = λopt ) とする。
この時、風速 V に対する最適回転数 ω ∗ は、周速比の関係より以下のように表現される。
ω∗ =
λopt
V
r
(3.13)
ゆえに、制御目的は次のように表される。
½
λopt
V
lim ω −
t→∞
r
39
¾
=0
(3.14)
モード 2 では、動作点における周速比は常に動作点 (λ∗ = λopt ) で一定であるので、(3.12) 式は風速 V
に線形に依存したモデルであることが分かる。さらに、風速 V の変動域がカットイン風速 Vci から定格
風速 Vrated までの間であることを考慮すると、線形補間タイプのモデルとして、次のような LPV モデ
ルが導出される。
µ
¶
Vrated − V
V − Vci
KT
ζ̇ =
A1 +
A2 ζ −
u
(3.15)
Vrated − Vci
Vrated − Vci
JT
|
{z
} | {z }
Ap (V )
ただし、
3.5
Bp
ζ = ω − ω ∗ , u = i − i∗
ρπr4 dCT ¯¯
dCT ¯¯
ρπr4
A1 =
Vci
Vrated
, A2 =
¯
¯
2JT
dλ λopt
2JT
dλ λopt
定格出力維持のための最適回転数制御
モード 3 における制御目的は、発電機の過負荷を防ぐために、風力発電機の出力 P = QA ω を定格出
力 Prated に維持することであるので、獲得する風エネルギーと定格出力の間に、次のような関係
1
Prated = ρπr2 CP (λ)V 3
2
(3.16)
が満足される必要がある。ゆえに、動作点の周速比は出力係数 CP の逆関数を用いて、次のように求め
られる。
¶
µ
2Prated
λ∗ = CP−1
(3.17)
ρπr2 V 3
(3.17) 式の出力係数の逆関数は、Fig.3.5 に示すように、最適周速比 (λ = λopt ) 以外では解を 2 個持つ
多価関数である。これは、ロータの回転数を上げる、または回転数を下げて運転することによって、定
格出力を維持することが可能であることを意味する。安全性の面から後者の方法を採用すると、(3.17)
式は一意に表わされる。(3.17) 式を利用し、風速 V に対する最適回転数 ω ∗ は
µ
¶
V −1 2Prated
∗
ω = CP
(3.18)
r
ρπr2 V 3
と表される。ゆえに、制御目的は次のように表される。
µ
½
¶¾
V −1 2Prated
=0
lim ω − CP
t→∞
r
ρπr2 V 3
(3.19)
モード 3 においてはモード 2 とは異なり、動作点の周速比 λ∗ は (3.17) 式に従い変動するために、(3.12) 式
の係数 Ap (V ) の偏微分項 dCT /dλ も風速 V に応じて複雑に変動する。ゆえに、係数が風速に対して Fig.3.6
のように変動するために、線形補間モデルを導出することができない。
そこで本章では 1 つの解決策として、Fig.3.7 に示すように、空力特性 (トルク係数) 曲線を 2 次関数
近似する。
CT (λ) = c2 λ2 + c1 λ + c0
(3.20)
(3.20) 式を (3.8) 式に代入して線形化すると、次に示すような風速 V と回転数 ω の 2 つの変数をスケ
ジューリングパラメータに持つ LPV モデルが導出される。
ẋ = Ap (ω, V )x + Bp u
ただし、
Ap (ω, V ) =
µ
¶
1
1
KT
ρπr5 c2 ω + ρπr4 c1 V , Bp = −
JT
2
JT
40
(3.21)
モード 3 においては、風速 V と最適回転数 ω ∗ は (3.18) 式に従い、Fig.3.8 に示すような関係で表され
るので、本章では、それを包含するような三角形領域を選択し、その端点 Ak , (k = 1, 2, 3) を用いて次
のような線形補間タイプの LPV モデルを導出する。
µ
¶
KT
ẋ = Ap (ω, V )x+ −
u
JT
(3.22)
ただし、
Ap (ω, V ) = (1 − α)A1 + α(1 − β)A2 + αβA3
µ
¶
1
1
5
4
ρπr c2 ωi + ρπr c1 Vi , (i = 1, 2, 3) ,
Ai =
JT
2
α=
1
(1 − α)ω1 + αω2 − ω
(1 − α)V1 + αV2 − V
(p1 V − p2 ω + p3 ), β =
=
p4
α(ω2 − ω3 )
α(V2 − V3 )

p1



 p
2

p
 3


p4
= ω2 − ω3
= V2 − V3
= −V1 (ω2 − ω3 ) + ω1 (V2 − V3 )
= ω2 (V3 − V1 ) + ω3 (V1 − V2 ) + ω1 (V2 − V3 )
ここで (3.22) 式について説明する。まずはじめに Fig.3.9 に示すような三角形領域を考える。次に、
線分 V2 V3 に対して平行な線分 V12 V13 を考える。ただし、点 V12 と点 V13 はそれぞれ線分 V1 V2 、V1 V3
上にあるものとする。この時、点 V12 = (x12 , y12 ) が線分 V1 V2 を α : (1 − α) に内分するものとすると、
点 V13 = (x13 , y13 ) も同様に、線分 V1 V3 を α : (1 − α) に内分するために、次のような関係式が満足さ
れる。
(
x12 = (1 − α)x1 + αx2
(3.23)
y12 = (1 − α)y1 + αy2
(
x13 = (1 − α)x1 + αx3
y13 = (1 − α)y1 + αy3
(3.24)
y12 − y13
y13 x12 − y12 x13
x+
x12 − x13
x12 − x13
(3.25)
線分 V12 V13 上の点 V (x, y) は、
y=
という関係が導かれるので、(3.25) 式に (3.23)、(3.24) 式を代入し、整理すると
·½
α(y2 − y3 )
1
y=
x+
+ (1 − α2 )x1 y1 + α(1 − α)y1 x2 + α(1 − α)x1 y3
α(x2 − x3 )
α(x2 − x3 )
¾¸
¾ ½
2
2
2
+ α x2 y3 + (1 − α )x1 y1 + α(1 − α)y1 x3 + α(1 − α)x1 y2 + α x3 y2
(3.26)
と表すことが出来る。また、α 6= 0 と仮定すると、
½
¾
½
2
x2 (y3 − y1 ) + x3 (y1 − y2 ) + x1 (y2 − y3 ) α + y1 (x2 − x3 ) + x1 (y3 − y2 )
¾
+ (y2 − y3 )x − (x2 − x3 )y α = 0
41
(3.27)
α=
(x2 − x3 )y − (y2 − y3 )x − y1 (x2 − x3 ) + x1 (y2 − y3 )
x2 (y3 − y1 ) + x3 (y1 − y2 ) + x1 (y2 − y3 )
となる。さらに、点 V が線分 V12 V13 を β : (1 − β) に内分するものとすると、
(
x = (1 − β)x12 + βx13
y = (1 − β)y12 + βy13
という関係が満足される。ただし、

x12 − x
(1 − α)x1 + αx2 − x


=

 x12 − x13
α(x2 − x3 )
β=


y −y
(1 − α)y1 + αy2 − y

 12
=
y12 − y13
α(y2 − y3 )
である。(3.30) 式を整理すると
" #
" #
" #
" #
x
x1
x2
x3
= (1 − α)
+ α(1 − β)
+ αβ
y
y1
y2
y3
(3.28)
(3.29)
(3.30)
(3.31)
となる。仮定を満足しない α = 0 の場合について考えると、(3.31) 式から [x, y]T = [x1 , y1 ]T となって
おり、正しく表現されることが分かる。ゆえに、2 つのスケジューリングパラメータを持つシステムは、
1 つの平面内で変動するために、(3.22) 式のような線形補間タイプの LPV システムに変換することが
できる。
3.6
風力発電機に対する LPV 制御系設計
本章では LPV コントローラを設計する際、混合感度問題をベースとして制御系設計問題を考慮する
ために、Fig.3.10 に示す制御システム構造を用いる。ここで、WI と WD はそれぞれ積分重み、微分重
み、ωc は回転数指令信号、z は制御評価信号、y は観測出力を意味する。上記の制御システム構造につ
いて、拡大系を構成すると次のようになる。








ζ̇
x˙I
z1
z2
xI


 
 
 
=
 
 
 
Ap (V )
−WI
01×1
WD Ap (V )
01×1
01×1
01×1
I1×1
01×1
I1×1
01×1
Bp
WI
01×1
01×1 01×1
01×1 WD Bp
01×1 01×1



ζ


  xI 


 ω 
c 


u
(3.32)
上記のモデルに対して、以下に示す制御仕様を満足する極配置制約付き LPV 出力フィードバックコン
トローラを設計する。
• 閉ループ系は内部安定である。
• ωc から z までの H∞ ノルムが γ 以下である。
すべての制御仕様は LMI として実現されるので、カットイン風速 Vci の時と定格風速 Vrated の時のそ
れぞれ端点での数学モデルに対して LMI を作成し、すべてを同時に γ を最小化する最適化問題として
解くことにより、LPV コントローラを設計する。設計法の詳細については、2 章を参考にしてもらい
たい。
42
3.7
数値シミュレーション
本章では LPV 制御の有効性を検証するために、数値シミュレーションを行った。風速データとして
は、九州大学応用力学研究所のフィールド試験で計測された実風速データから、運転モードの風速領域
に合致したものを選択した。モード 2 の制御目的は、風力発電機を最適周速比 (λopt = 4.51) で運転する
ことであるので、周速比の時系列データの結果を Fig.3.11-(a) に、一方で、モード 3 の制御目的は、風
力発電機の出力を定格出力 (Prated = 300) に押さえ込むことであるので、出力の時系列データの結果を
Fig.3.11-(b) に示す。比較のために、(3.1) 式で示した回転数の 2 乗に比例する制御 (以後、従来 ω 2 制御
と表記) の結果も示す [21]。ここで、数値シミュレーション結果の定量的な評価のために、次のような
分散評価指標を定義する。

N

1 X

2

(xk − 4.51)2
in mode 2
σ
=

 4.51 N
k=1
(3.33)
N
X

1

2
2


(xk − 300.00) in mode 3
 σ300 = N
k=1
(3.33) 式の第 1 式は最適周速比 4.51 からの、第 2 式は定格出力 300 からの分散値を表す。さらに、モー
ド 2 における制御性能評価のために、次のような評価指標を導入する [19]。
Rt
1
2
3
2 ρπr CP (λ)V dt
2
3
0 2 ρπr CP max (λ)V dt
Ep (t) = R t 01
(3.34)
これは、時刻 t までに、流体力学的に獲得可能な風の総エネルギーに対して、どの程度の風エネルギー
が実際に獲得されているかを表す指標である。それらの評価指標の結果を Table3.2 に示す。
モード 2 の制御においては、Fig.3.11-(a) から、LPV 制御は従来 ω 2 制御に対して最適周速比近傍での
運転を実現しており、Tabel 3.2 の結果からも、最適周速比からの分散が、従来 ω 2 制御が 0.439 に対し、
LPV 制御は 0.083 と、LPV 制御は従来 ω 2 制御に対して、約 1/5 の分散値となっている。さらに、獲得
エネルギー EP (300) の比較で見ても、従来 ω 2 制御は 0.926 に対し、LPV 制御は 0.979 と約 6 % の改善
が見られる。また、モード 3 の制御についても、Fig.3.11-(b) から、LPV 制御は従来 ω 2 制御と比較して
出力変動が大きく低減している。Table 3.2 の定格出力 300 からの分散を見ると、従来 ω 2 制御は 252.839
であるのに対して、LPV 制御は 65.162 と、約 1/4 に分散が抑えられている。ここで、Fig.3.11-(b) に
は、出力値が 300 からオフセットを残して制御されているように観察されるが、これは、関数近似に
よる数値誤差の影響と考えられる。
Table 3.2: Controller performance in numerical simualtion
mode 2 region
Conventional ω 2 Control
LPV Control
2
σ4.51
0.439
0.083
EP (300)
0.926
0.979
mode 3 region
Conventional ω 2 Control
LPV Control
2
σ300
252.839
65.162
すべての運転モードの数値シミュレーションの結果において、高周波数の回転数変化に追従できてい
ない点は、制御コントローラの帯域幅の問題や、回転数を加速させるためには、電流を最小にして負荷
トルクを小さくし、回転数が大きくなることを待つしかないという風車の回転数制御問題の特徴に起因
するものと考えられる。
43
3.8
3.8.1
風車シミュレータを用いた検証実験
風車シミュレータの考え方と装置の概要
風力発電機の開発において、そのシステムが実際の風環境の中で効率的な運転を行うかどうかをテス
トすることは非常に重要である。その最も簡単な方法として、風洞試験による性能評価が挙げられるが、
基本的には風洞で実現可能な風速は定常風であるために、乱れやガストを含む自然風とは大きな違いが
ある。一方、より現実的な方法として、フィールド試験も考えられるが、自然風は季節ごとに特徴的な
様相を持つために、短期間の間にそのシステムの性能を評価することは非常に困難である。そこで、こ
れらの問題を解決するために、本章では風車シミュレータ (Wind Turbine Simulator) と呼ばれる風力
発電機制御システム開発支援装置を提案する。その装置全体の写真を Fig.3.12 と Fig.3.13、風車シミュ
レータの全体の構成を Fig.3.14 に示す。
風車シミュレータには、大きなロータ翼などは取り付けられていない。そこで、風車の持つ特有の空
力特性を実現するために、Fig.3.15 のようなプロセスを PC 上で実行し、それを AC サーボモータに指
令するものとする。本来、指令トルクが正確に、そして迅速に発生するようにフィードバック制御をす
る必要があるが、その場合、閉ループ系の極が虚軸上に配置させてしまうことを考慮して、本研究では
フィードフォワード制御のみで指令トルクを実現することとした。これにより、大きなロータを接続す
ることなく、実験室内において仮想風車を作り出し、様々な風速環境下での制御システムの性能評価を
容易に行うことが可能である。この開発支援装置の他の利点として、
• 風洞試験で実現困難な突風 (ガスト) モデルを実現可能
• 実装した際の、大きな機械的、電気的な破損を回避
• 風車の開発スピードを大幅にアップし、自然エネルギー利用を促進する
などが挙げられる。
3.8.2
検証実験
ここでは、風車シミュレータを用いた LPV 制御の性能検証実験の結果を示す。最初にモード 2 の実験
について説明する。対象とする風レンズ風車において、モード 2 の風速領域は 1(m/sec) から 11(m/sec)
であるが、摩擦の関係から速やかな回転を開始するような適当なトルクを発生する風速、また、大きな
風速時においては、非常に危険な回転数になる事を考慮して、ここでは適切なモード 2 の風速領域を、事
前の予備実験から 6(m/sec) から 9(m/sec) と設定することにした。風速としては、応用力学研究所で行
われたフィールド試験の計測データから適切な風速データを採用した。実験結果を、Fig.3.16, Fig.3.17
と Table 3.3 に示す。Fig.3.16-(a) を見ると、最適回転数と実験値の間に大きな違いがあることが分か
る。その結果、Fig.3.16-(b) から周速比も最適値からの分散が大きくなっており、Ep (300) も 0.60 と約
4 割のエネルギーを失ってしまっている。一方で、LPV 制御は Fig.3.17-(a) からも、最適回転数に非常
に良く追従されており、Fig.3.17-(b) から周速比も最適値近傍で運動されていることが分かる。効率を
見ると、Ep (300) = 0.98 とエネルギー効率は非常に高く、従来 ω 2 制御の約 1.5 倍もの高い効率が確認
される。
Table 3.3: Control performance index in experiment using Wind Turbine Simulator
mode 2 region
Conventional ω 2 Control
LPV Control
2
σ4.51
2.150
0.180
Ep (300)
0.600
0.980
44
mode 3 region
Conventional ω 2 Control
LPV Control
2
σ300
36241.780
50.324
次に、モード 3 の実験について説明する。風レンズ風車のモード 2 の風速領域は、11(m/sec) から
22(m/sec) であるが、回転数が上がりすぎることを考慮して、適切なモード 3 の風速領域を 13(m/sec)
から 18(n/sec) と設定した。風速としては、モード 2 と同様にして、応用力学研究所のフィールド試験
の計測データの中から適切なデータを選択した。実験結果を、Fig.3.18, Fig.3.19 と Table 3.3 に示す。
モード 3 においても、Fig.3.18-(a) と Fig.3.18-(b) の結果を見ると、従来 ω 2 制御は最適回転数、最適周
2 = 36241.78 といずれも大きな値である。この
速比に追従できておらず、その結果、出力の分散も σ300
大きな原因としては、物理定数の不確かさや、モデルに反映されない動摩擦力などの影響により、トル
クバランスが大きく崩れ、ω 2 制御ではこの問題を解決することは出来ないことが考えられる。これは、
K.E.Johnson らも文献の中で同様な指摘をしている [26]。また、モード 3 は不安定な強非線形性を有す
ることも原因の一つであろう。一方で、Fig.3.19-(a) と Fig.3.19-(b) の結果から、LPV 制御は最適回転
数、最適周速比に十分良く追従できている。ゆえに、出力変動を見ると、300W 近傍に押さえ込まれ、
2 = 50.324 と従来 ω 2 制御のおよそ 1/600 となっており、効率的な運転が実現されている様子が確認
σ300
される。
また、本風車シミュレータを用いた実験においては別途参考までに、固定ゲイン制御方式の検討も
行っている。モード 2 においては線形パラメータ変動制御方式とあまり大きな相違は見られなかった一
方で、モード 3 においては風速が大きくなる中で回転数をあえて減速させていくために、非常に不安定
となり、固定ゲイン制御方式では動特性の変化に対応出来ずに不安定化 (過回転状態) してしまうこと
が確認されている。
3.9
まとめ
本章では、小型の固定ピッチ型水平軸風車である風レンズ風車に対して、LPV 制御の有効性を確認
するために、数値シミュレーション、及び、風車シミュレータと呼ばれる風力発電機開発支援装置を利
用した検証実験を行った。そこで下記のような結論が得られた。
• モード 2 においては動作点が一定であるために、風速をスケジューリング変数に持つ LPV モデル
を導出することが可能であるが、モード 3 においては動作点が風速に応じて変化するために LPV
モデルを導出することができない。そこで、風車の大まかな空力特性を 2 次関数近似により捉え
ることにより、回転数と風速の 2 つのスケジューリング変数に持つ LPV モデルを導出することが
可能である。
• 従来から採用されている ω 2 制御による方法は、最適回転数追従制御が出来ずに、最適周速比
(λ = 4.51, P = 300 W) での近傍での運転をすることができない。その結果、モード 2 において
は風エネルギー取得の大きなロス、モード 3 においては、定格出力値からの大きな出力変動を発
生させてしまう。この大きな原因としては、物理定数の不確かさや、モデルに反映されない動摩
擦力などの影響により、トルクバランスが大きく崩れ、ω 2 制御ではこの問題を解決することは出
来ないことが考えられる。これは、文献 [26] においても、同様なコメントがされている。
• LPV 制御は、モード 2、モード 3 の両モードにおいて、非常によく最適回転数制御を実現するこ
とができ、最適周速比近傍での運転が実現可能である。
本研究の目的は、洋上風力発電機などで使用されるようなメガワットクラスの大型風力発電機に応用
するための基礎研究であるために、依然として多くの課題が残されている。例えば、
• 一般的に、大型風力発電機はロータ翼と回転数の 2 つを制御することにより効率化が行われてい
る。本研究は、回転数制御についてのみを検討しているに過ぎないために、ロータ翼についても
考慮する必要がある。
45
• 本章で提案した風車シミュレータは、確かに空力性能については実現することが可能であるが、
装置の関係からロータの慣性を考慮することができていない。実際の大型風車は慣性力が大きく
回転数制御に影響を及ぼすために、その点についても考慮しなければならない。 • 本章において応用した LPV 制御方式を、フィールド試験により検証し、問題点と有効性について
明確にする必要がある。
などが挙げられる。
46
Power [w]
Prated
Mode 1
Mode 3
Mode 2
Vci
0
Mode 4
Vrated
Vco
Wind Speed [m/sec]
Fig. 3.1: Operational modes of wind turbine
Fig. 3.2: ”Wind Lens” developed by Research Institute for Applied Mechanics in Kyushu University
ω
L
R
V
i
e
KE ω
Battery
JT
Fig. 3.3: Wind turbine system model
47
1.5
( λ opt ,Cpmax )
1
Cp-Curve
0.5
CT-Curve
0
1
2
3
4
5
6
Tip Speed Ratio λ
7
8
Fig. 3.4: Power coefficient CP and torque coefficient CT curves of Wind Lens
1.5
2P rated
Operational Point
in Mode 2
Cp ( ρπ r2 V3
1
(
-1
Operational Point
in Mode 3
0.5
0
1
2
3
4
5
6
7
Tip Speed Ratio λ
8
Fig. 3.5: Target operational point
2.5
ρπ r4 dC T
V
2J T dλ
2
Ap
Ap (V ) =
1.5
A co
1
0.5
0
A
rated
11
13
15
19
17
Wind Speed [m/sec]
Fig. 3.6: Relationship between
ρπr 4 dCT
2JT V dλ
48
21
and wind speed V in mode 3
0.3
0.2
0.1
Solid Line : CT - Curve
Dash Line : Approx. CT - Curve
0
2
4
Tip Speed Ratio
8
6
λ
Fig. 3.7: Approxmated torque coefficient CT curve of Wind Lens
Wind Speed V [m/sec]
A1 ( ω1 , V1)
20
16
( ω3 , V 3)
A3
A2
12
( ω2 , V 2)
80
100
120
140
160
Rotational Speed ω [rad/sec]
Fig. 3.8: LPV system varying with rotational speed ω and wind speed V in mode 3
V 1 (x 1 , y 1 )
α
α
V13 (x 13, y 13)
1-α
1-β
β
V ( x, y )
V 3 (x 3 , y 3 )
V12 (x 12, y 12)
1-α
V2 (x 2, y 2)
Fig. 3.9: Description of the LPV system using two scheduling parameters
49
1
V
λ*
r
ωc
WI
s
-
xI
i
contr.
Wind
Turbine
WD s
z1
z2
y
Wind Speed Data
Wind Speed (m/s)
Wind Speed (m/s)
Fig. 3.10: Interconection for the mixed sensivity problem
8
4
0
100
200
18
16
14
0
300
100
200
300
Time (sec)
Time (sec)
(b) Wind speed profile
(a) Wind speed profile
8
Rated Power
6
Power (W)
Tip Speed Ratio
λ
400
4
2
0
0
Optimal Tip Speed Ratio
100
200
300
200
0
300
100
200
300
Time (sec)
Time (sec)
(b) Conventional ω 2 control
(a) Conventional ω 2 control
8
Rated Power
6
Power (W)
Tip Speed Ratio
λ
400
4
Optimal Tip Speed Ratio
2
0
0
100
200
300
300
200
0
100
200
300
Time (sec)
Time (sec)
(a) Linear parameter varying control (b) Linear parameter varying control
Fig. 3.11: Comparison of control performance between conventional ω 2 control and LPV control in
simulational results ((a):mode 2 region, (b):mode 3 region)
50
Fig. 3.12: Wind Turbine Simulator (photo 1)
Fig. 3.13: Wind Turbine Simulator (photo 2)
51
Wind Data
Battery
Electrical Power
CT-Curve
Generator Controller
(Variable Resistor)
Electrical Current
Command
PC
Aero Torque
Command
Measured Rotational
Speed
Electrical Power
24
Generator Torque & Rotational AC-Servo Motor
Speed Detector
Fig. 3.14: Wind Turbine Simulator
Field Test Measured Wind Speed
V
Calculation of Tip Speed Ratio
λ
Decision of Torque Coefficient
from Torque Curve
CT
Calculation of Aerodynamical Torque
T=
1
rp r 3 V 2 CT
2
Command to AC ServoMotor
Rotational Speed is changed
Fig. 3.15: Caluculation flow of aerodynamical torque in the Wind Turbine Simulator
52
Tip speed ratio
Rotational Speed [rad/sec]
8
150
Optimum Speed
100
50
Experiment
0
0
20
40
6
2
0
0
60
Optimum Tip speed ratio
4
Experiment
20
100
Experiment
0
40
60
Wind speed profile [m/sec]
Output power [W]
200
20
60
(b) Tip Speed Ratio
(a) Rotational Speed (rad/sec)
0
40
Time [sec]
Time [sec]
15
10
Wind Speed
5
0
0
40
20
Time [sec]
60
Time [sec]
(c) Output Power (W)
(d) Wind Speed Profile (m/sec)
8
150
Tip speed ratio
Rotational Speed [rad/sec]
Fig. 3.16: Experimental results in mode 2 region using conventional ω 2 controller
Optimum speed
100
20
Optimum Tip speed ratio
4
Experiment
2
Experiment
50
0
6
40
0
60
20
Experiment
100
0
40
60
Wind speed profile [m/sec]
Output power [W]
200
20
60
(b) Tip Speed Ratio
(a) Rotational Speed (rad/sec)
0
40
Time [sec]
Time [sec]
15
10
Wind Speed
5
0
0
40
20
Time [sec]
Time [sec]
(c) Output Power (W)
(d) Wind Speed Profile (m/sec)
Fig. 3.17: Experimental results in mode 2 region using LPV controller
53
60
Tip Speed Ratio
Rotational Speed [rad/sec]
4
200
Optimum Speed
100
Experiment
0
0
100
200
Time [sec]
Optimum Tip speed ratio
2
Experiment
0
0
300
100
200
Time [sec]
300
(b) Tip Speed Ratio
(a) Rotational Speed (rad/sec)
Wind Speed [m/sec]
Output Power [W]
400
300 W
200
Experiment
0
0
100
200
20
Wind Speed
15
0
300
100
200
300
Time [sec]
Time [sec]
(d) Wind Speed Profile (m/sec)
(c) Output Power (W)
4
200
Tip Speed Ratio
Rotational Speed [rad/sec]
Fig. 3.18: Experimental results in mode 3 region using conventional ω 2 controller
Optimum Speed
100
Experiment
0
0
100
200
Time [sec]
Optimum Tip speed ratio
2
Experiment
0
0
300
100
200
Time [sec]
300
(b) Tip Speed Ratio
(a) Rotational Speed (rad/sec)
Experiment
Wind Speed [m/sec]
Output Power [W]
400
300 W
200
0
0
100
200
300
20
Wind Speed
15
0
100
200
Time [sec]
Time [sec]
(d) Wind Speed Profile (m/sec)
(c) Output Power (W)
Fig. 3.19: Experimental results in mode 3 region using LPV controller
54
300
第4章
大水深ライザー管の動特性を考慮したリエントリー作
業の自動化
本章では、大水深ライザー管のリエントリー作業の自動化について検討する。本章の構成は以下の通
りである。はじめに、Hamilton の原理から導出した大水深ライザー管の運動方程式について説明する。
次に、その運動方程式は無限自由度を有しており、非常に取り扱い難いために、モード展開近似を行い
有限次元化を行う。その有限次元モデルから微小変位の仮定を行い、ライザー管上端浮体部の速度をス
ケジューリング変数に持つライザー管の LPV モデルを求める。そのモデルに対しての流体抗力の及ぼ
す影響を調べた後、LPV コントローラを導出する。そして、性能検証のために数値シミュレーション
を用いて効果を確認する。最後に、九州大学応用力学研究所深海機器力学実験水槽で実施した水槽実験
の結果について示したのち、本論文のまとめと今後の課題について述べる [27][28][29][30][31]。
4.1
はじめに
地球規模の環境変動、地震発生のメカニズムの解明を進めるために、大水深ライザー管を用いた地
球深部掘削プロジェクト IODP が JAMSTEC を中心に進められており、その大きな役割を担う地球深
部探査船「ちきゅう」がいよいよ就航した (Fig.1.3)。この船に搭載されている長大弾性管 (通称、ライ
ザー管) の複雑な非線形運動を制御することは非常に困難であるものの、外乱要因を多く持つ海上での
作業は緊急退避などの必要に迫られるため、管を接続された防噴装置頂部から切り離し、回避後は早急
に再接続する必要がある。それをリエントリー作業と呼ぶ。その際、Fig.1.4 に示すように掘削作業者
は、ROV に備え付けられたカメラで状況を観察し、長年の経験で培った勘に頼りながら船舶を位置制
御し、先端部を目標点まで制御する。しかしながら、現状では、このリエントリー作業に多大な労力を
要しており、掘削者の作業時間低減化のための高精度ライザー管先端位置制御は早急に求められている
必要不可欠な技術である。
ライザー管に関する研究は盛んに行われており、近年においては、数値流体力学などの発展もあり、渦
励振などの流体力学的視点から管の動特性を考えようとする研究に重点が置かれるようになった [32][33]。
しかしながら、実際に運用上必要なライザー管の制御技術に関する研究は、鈴木らの研究 [34] が行われ
て以来、ほとんど発表されていない。最近では、制御理論に重点を置いた研究として、M.P.Fard の受
動性に基づく非線形制御に関する研究などがあるが、非線形無限次元系の安定性を保証する Lyapunov
関数の構築法などについて大きな成果が得られているが、まだ、水中の非線形運動を取り扱える実用的
なレベルにまでには到達していない [35]。
そこで本章では、大水深ライザー管の動特性が速度変動により大きく変化する点に着目し、LPV 制
御の有効性を数値シミュレーション、及び、水槽実験により確認することが目的である。
55
4.2
4.2.1
大水深ライザー管の状態方程式の導出
運動の数学モデル
本章では小寺山、鈴木らの方法に従い、ハミルトンの原理を用いて大水深ライザー管の運動方程式を
導出する [32][33][34][36]。はじめに Fig.4.1 に取り扱うライザー管の座標系を定義する。µ, η, θ はそれ
ぞれ、ライザー管の伸びと撓み、浮体部に対するライザー管剛体モードの回転角度、r は浮体部の水平
座標を表す。ここでは、運動方程式の導出の簡単化のために、次のような仮定を導入する。
• 断面形状は長さ方向に一様で、面内変形しない
• オイラーベルヌーイ梁として取り扱う
• ライザー管の回転方向の速度との相対速度を考慮した Morison 式から流体抗力を与える
• 入力としては、ライザー管上端浮体部の水平方向の推進力のみ
• 浮体部の Heave 運動は無いものとする
その結果、本章で取り扱う大水深ライザー管の運動方程式は次のように記述される [37]。なお、導出法
の詳細については付録を参考にしていただきたい。
水平方向:
½
¾
½
Z l·
2
(m̃ + mA ) cos θ η̈ + (z + µ)θ̈ + r̈ cos θ + 2µ̇θ̇ − η θ̇ + m̃ sin θ µ̈ − η θ̈ + r̈ sin θ − 2η̇ θ̇
0
¯
¯½
¾
¾¸
¯
¯
1
2
¯
¯
−(z + µ)θ̇ + ρo Cd Do cos θ¯η̇ + (z + µ)θ̇ + ṙ cos θ¯ η̇ + (z + µ)θ̇ + ṙ cos θ dz = u
(4.1)
2
回転方向:
½
¾
½
¾
Z l·
2
2
(m̃ + mA )(z + µ) η̈ + (z + µ)θ̈ + r̈ cos θ + 2µ̇θ̇ − η θ̇ − m̃η µ̈ − η θ̈ + r̈ sin θ − 2η̇ θ̇ − (z + µ)θ̇
0
¯
¯½
¾¸
¯
¯
1
+m̃g{η cos θ+(z+µ) sin θ}+ ρo Cd Do (z+µ)¯¯η̇+(z+µ)θ̇+ṙ cos θ¯¯ η̇+(z+µ)θ̇+ṙ cos θ dz = 0 (4.2)
2
たわみ方向:
¯
½
¾
¯
1
2
00 00
0 0
(m̃ + mA ) η̈ + (z + µ)θ̈ + r̈ cos θ + 2µ̇θ̇ − η θ̇ + (EIη ) − (Te η ) + m̃g sin θ + ρ0 Cd Do ¯¯η̇
2
¯½
¾
¯
+(z + µ)θ̇ + ṙ cos θ¯¯ η̇ + (z + µ)θ̇ + ṙ cos θ = 0
(4.3)
伸び方向:
½
¾
2
m̃ µ̈ − η θ̈ + r̈ sin θ − 2η̈ θ̇ − (z + µ)θ̇ − EAt µ00 − m̃g cos θ = 0
(4.4)
境界条件:
µ = η = 0 (at z = 0), EAt µ0 = 0 (at z = l), EIη 00 = 0 (at z = 0, l),
(−EIη 000 + Te η 0 )η = 0 (at z = l)
(4.5)
ここで、m̃ はライザー管の単位長さあたりの質量、mA はその付加質量、l はライザー管の全長、ρ0 は
流体密度、Cd は抗力係数、Do は弾性管の外径、g は重力加速度、E はヤング率、I は断面 2 次モーメン
ト、At はライザー管の断面積、Te はライザー管に働く張力、そして、u は浮体部に働く水平方向の推
0
進力を表している。ẋ (ドット) は時間微分、x (プライム) は z による空間微分を意味する。鈴木らは、
ライザー管の水深方向に幾つかのスラスターを配置し、それらを同時に使って制御を行う方法を提案し
ているが [34]、本章では浮体部に取り付けられたスラスターのみでライザー管を制御するものとする。
またその際には、簡単化のために潮流や波浪などの環境外乱などの影響は考えない。
56
4.2.2
モード展開による運動近似化
このライザー管の運動方程式を直接取り扱い、数値計算することは非常に困難であるために、本章で
は (4.1)∼(4.5) 式の運動方程式をモード分解近似して有限次元近似を行う。境界条件を満足するように、
(4.6) 式のような直交関数を用いて級数展開を行う。

N
X

1


Un (t) sin(n − )πz

 µ(z, t) =
2
n=1
(4.6)
N
X




Wn (t) sin(nπz)
 η(z, t) =
n=1
ここで、n は振動モードの次数を意味する。なお、上記のモード関数は無次元化された式に対するもの
である。本章では、8 次モードまでの数値シミュレーション環境を構築した。Fig.4.2 には、流体抗力の
影響効果を観察するために、ある定常状態からの自由応答の数値シミュレーション結果 ((a):ライザー
管上端部位置、(b):ライザー管全体傾斜角、(c):1 次モードの振動幅、(d):ライザー管最下端部位置) を示
す。すべての変数はライザー管の全長 l で無次元化して表示している。抗力係数 Cd をゼロとして流体
抗力の効果を取り除いたもの (破線: without hydrodynamical damping, Cd = 0.0) は、ダンピングの効
果を有さないために、流体抗力の効果があるもの (実線: with hydrodynamical damping, Cd = 1.17) に
比べて、振動が収まることなく全く異なった挙動をすることが分かる。一般的に、抗力係数 Cd は、Kc
数 (Keulegan-Carpenter 数) と渦励振の影響により大きく変化することが報告されている [38][39][40]。
よって、この流体抗力の変化が運動制御に悪影響を及ぼさないよう、ロバストなコントローラを設計す
ることが非常に重要である。
4.3
4.3.1
大水深ライザー管のリエントリー制御
大水深ライザー管の LPV モデルの導出
大水深ライザー管のリエントリー作業においては、鈴木らが指摘しているように [34]、ライザー管の
高次振動モードが発生するようなことは考えにくく、仮に発生したとしても流体抗力により速やかに消
滅していくものと考えられる。よって、低次の振動モードのみを考慮して制御系設計を行っても、十分
な制御性能を保証する低次のコントローラを設計することが可能である。そこで、本章ではライザー管
の運動の 1 次振動モードのみの数学モデルに基づき制御系設計を行うものとする。本章で取り扱う運動
方程式を行列表現すると (4.7) 式となる。ここで、x1 = [r θ W1 ]T , x2 = [ṙ ω Ẇ1 ]T であり、ω は
剛体モードの回転角速度を意味する。この時、θ ¿ 1 (sin θ ' θ, cos θ ' 1), ω ¿ 1, W1 ¿ 1 の微小
変位の仮定をすると、(4.8) 式のような状態方程式が導出される。|ṙ| の変化とともに A 行列が線形パラ
メータ変動することが分かる。
"
I3
03×3

03×3
M
#"
ẋ1
ẋ2
#
"
+
cos2 θ
03×3 −I3
K
C
1 + mA


 1 + mA
2 sin θ
M =
cos θ −
W1

2
π


2(1 + mA )
cos θ
π

0
 0
K=

0
#"
x1
x2
#
"
+
03×1
I3
1 + mA
2 sin θ
cos θ −
W1
2
π
1 + mA 1 2
+ W1
3
2
1 + mA
π

0
0

0
0


π 2 α2
α2 π 4 EI
+
g
0
2m̃l4
4l
57
#
"
G=
03×1
F
#
2(1 + mA )
cos θ
π
1 + mA
π
1 + mA
2
(4.7)











4
− sin θω


π

 l


ρ0 Cd Do |ṙ| 
+
W1 ω

 2m̃


0
1
2(1 + mA )
0 − sin θω −
cos θW1 ω

2
π

4mA 2

C= 0
W1 ω
−

3π

1 + mA
0
−
W1 ω
2

0

1
2
α2 g 
 sin θ + W1 cos θ
G=
2
π
l 

2
sin θ
π
ẋ = Ax + Bu ,


A=
03×3
A21
−
 , A21
ζ
|ṙ| I3
1 + mA






B=





4.3.2

I3
 0


=
 0


0
03×1
4α2 (π 2 − 6)
m̃l2 (1 + mA )(π 2 − 8)
−
−
6α2
m̃l2 (1 + mA )
4α2 π
m̃l2 (1 + mA )(π 2 − 8)
cos3 θ
1
cos2 θ
2
2
cos θ
π
1
cos2 θ
2
1
cos θ
3
1
cos θ
π










,


h
F =
u 0 0
iT
,
l2
α= 2
π
r
m̃
EI
y = Cx , x = [r θ W1 ṙ ω Ẇ1 ]T
µ
2
cos θ
π
1
cos θ
π
1
cos θ
2
¶
8
α2 g
2
π − 8 l(1 + mA )
6α2 g
−
l(1 + mA )
4πα2 g
−
l(1 + mA )(π 2 − 8)
3+
(4.8)
4π
12α2 g
ξ+
2
(1 + mA )(π − 8)
lπ(1 + mA )
24α2 g
−
lπ(1 + mA )
2π 2
−
ξ
(1 + mA )(π 2 − 8)





α2 π 4 EI
π 2 α2



ξ
=
+
g


2m̃l4
4l
 , C = [ I3 , 03×3 ],




 ζ = l ρ0 Cd Do


2m̃

固定パラメータ開ループ系モデルの動特性解析
この大水深ライザー管の LPV モデルに対しての流体抗力の影響を見るため、|ṙ| = 0 と |ṙ| = 1.5 の
場合について、|r/u|、|θ/u|、|W1 /u| のボード線図、極配置を調査した。ただし、|ṙ| = 1.5 は幾つかの
数値シミュレーションの結果から採用している。Fig.4.3 と Fig.4.4 にその結果を示す。ボード線図 (実
線: |ṙ| = 0、破線: |ṙ| = 1.5) からは、|ṙ| = 0 の時に |r/u| のゲインが低振動数領域で大きく、5Hz 近傍
で共振と反共振のピークが存在することが分かる。一方で、極配置 (×印: |ṙ| = 0、○印: |ṙ| = 1.5) か
らは、|ṙ| = 1.5 の時には、極が虚軸側と安定側 (左実軸上) に分かれているのに対して、|ṙ| = 0 の時に
は、すべての極が虚軸上に存在することが分かる。したがって、ライザー管のリエントリー作業に対す
る制御系設計の指針としては、|ṙ| = 0 の時の共振ピークを低減化し、物理パラメータの不確かさなど
の影響により、極が不安定側に容易に移行することのないように、虚軸上に張り付いた極を出来るだけ
安定側に極配置することが必要である。
58








4.3.3
大水深ライザー管の LPV コントローラの導出
本章では、Fig.4.5 に示すような制御システム構造の下で制御系設計を行うこととする。設計の考え
方を Fig.4.6 に示す。この一般化プラントは次の (4.9) 式のように表される。


 

ẋ
B
B
x
A
  


 z  =  Wd CA Wd CB Wd CB   w 
y
C
03×1
03×1
u
(4.9)
である。ここで、Wd は微分重み関数であり、z は評価信号、w は環境からの外乱入力、Wd は微分重み
を意味する。(4.9) 式に対して次のような制御仕様を満足する、極配置制約付き LPV 出力フィードバッ
クコントローラを設計する。
• 閉ループ系は内部安定である。
• w から z までの H∞ ノルムが γ 以下である。
すべての制御仕様は LMI として実現されるので、|ṙ| = 0 の時と |ṙ| = 1.5 の時の端点での数学モデル
に対して LMI を作成し、すべてを同時に γ を最小化する最適化問題として解くことにより、2 つの端点
コントローラを求め、それらのコントローラをスケジューリングさせることにより、LPV コントロー
ラを導出する [13]。設計法の詳細については 2 章を参考にしていただきたい。
前節と同様にして、|ṙ| = 0 と |ṙ| = 1.5 の場合について、閉ループ系の |r/w|、|θ/w|、|W1 /w| のボー
ド線図、極配置を調べる。Fig.4.7 と Fig.4.8 にそれを示す。ボード線図 (実線: |ṙ| = 0、破線: |ṙ| = 1.5)
からは、開ループ系で見られた 5Hz 近傍の共振ピークが低減化されており、極配置 (×印: |ṙ| = 0、○
印: |ṙ| = 1.5) からは、虚軸上にあった極が安定側にシフトしている様子が分かる。しかしながら、原点
近傍の極はほとんど動かすことが出来ずに、虚軸近傍に極配置されてしまっている。これは構造物の制
振問題で直面する回避できない現象であると思われる。
4.4
数値シミュレーション
本章では制御性能を評価するために数値シミュレーションを行った。Fig.4.9 にその結果 ((a):ライザー
管上端部位置、(b):ライザー管全体傾斜角、(c):ライザー管最下端部位置、(d):制御入力) を示す。すべ
ての変数は無次元値として表示している。また、性能比較のために、複数モデルに基づくロバスト制御
(Multi-model Based Robust、以後、MmBR 制御) の結果も同時に示す。MmBR 制御の結果は目標値
近傍で振動現象が発生している。この現象をより詳細に観察するために、Fig.4.10 にライザー管の挙動
を 2 次元的に表記したもの、Fig.4.11 に各モードの振幅の時系列の結果を示す。LPV 制御は目標点近
傍でも揺れる事なく制振されている一方、MmBR 制御は高次の振動モードまでもが発生していること
が分かる。これは、制御工学の観点から考えると、コントローラが |ṙ| = 0 のシステムに対応できずに、
その結果、虚軸近傍の極を不安定側に移行させ、高次の振動モードを誘起させたスピルオーバ現象と考
えられる。流体抗力の存在のために、大変形破壊を誘起するような大きな振動現象に発展することはな
いとは思われるが、疲労強度や作業効率などの観点から考えると、ぜひとも避けなければならない現象
である。
Table 4.1: Controller performance in numerical simulation
MmBR Control
LPV Control
Maximum Top Angle (deg)
3.783
1.777
59
Maximum Controller Output (N)
19.247
12.433
Table4.1 に MmBR 制御と LPV 制御のライザー管の上部傾斜角と制御入力の比較したものを示す。上
部傾斜角とは全体傾斜角とたわみによる傾斜の和 (本章では、1 次モードのたわみのみを考慮) である。
ライザー管の上部傾斜角の許容範囲は、ライザー中を通るドリルパイプとライザーが接触して破損する
可能性があるために、一般的に 3 度 (約 0.05 rad) 程度であると言われている [41]。LPV 制御は最も上
部傾斜角を押さえ込んでおり、MmBR 制御の上部傾斜角の約半分になっている。最大制御入力も過大
にならず、十分に可能な範囲に押さえ込まれていることを考慮すると、LPV 制御は MmBR 制御と比較
して、大きな上部傾斜角を必要とせずに迅速な制御が可能であると考えられる。
4.5
4.5.1
水槽実験による検証
実験の概要
数値シミュレーションの結果を踏まえ、九州大学応用力学研究所深海機器力学実験水槽 (65×5×7m)
において検証実験を行った。使用したライザー管モデルの基本的な物理定数等を Table 4.2 に示す。管の
材質はポリエチレンとテフロンの合成樹脂である。Fig.4.12 は、管を水面上から見たものである。この
管上端部を Fig.4.13 に示すようなパラレルリンクを用いて位置制御を行うことにより、ライザー管のリ
エントリー作業を模擬することとした。管全体の挙動は、計 14 台の CCD カメラを用いて、Fig.4.14 に
示すように配置されたマーカーの挙動を計測する。CCD カメラは、2 つの翼形状のフレームに 7 台ずつ
(CA1∼CA7,CB1∼CB7) 内包されており、Fig.4.15 に示すように配置されている。事前に Fig.4.16 に示
すようなマーカーを使ってキャリブレーションを行い、画像データを市販の解析ソフト HALCON を用
いて画像解析を行い、Fig.4.17 のように座標変換をして世界座標系で数値化することにした。Fig.4.18
に計測の全体的な流れを示す。
Table 4.2: Characteristics of the riser model in Research Institute for Applied Mechanics, Kyushu
University
Model length
Outer diameter
Inner diameter
4.5.2
6.5 (m)
22.5 (mm)
12.7 (mm)
Mass per length
Young’s Modulus
Bottom weight in water
0.4 (kg/m)
8.847 (MPa)
3.489 (N)
実験内容とその妥当性について
本章で行ったリエントリー実験は、Fig.4.19 に示すように、管を x 軸方向に 0.3 m 離れた地点を目標
点として制御するものである。その際、次の要求が満足されなければならない。
1. 迅速に目標点に到達し、振動発生を防止する。
2. 過大な制御入力を要求しない。
3. 上端部傾斜角を許容範囲内 (3 deg) に押さえる。
一つ目は作業効率、二つ目は船舶等に装備されたスラスタなどの性能限界、三つ目は、ライザー中を通
るドリルパイプとライザーが接触して破損する可能性からくる安全面からの要求である。本章で行った
実験において、上端部傾斜角を計測することができなかったため、本章においては、1, 2 の要求に着目
して実験を行った。
本章で行った実験は、画像計測処理のリアルタイム性がないために、予め、管上端部の動きを数値シ
ミュレーションにより計算し、その動きをパラレルリンクで実現するという実験手法を採用した。した
がって、本実験では、開ループ実験になってしまい正確な閉ループ系の挙動を得ることができない。本
60
実験の意味を考慮するために、伝達関数表現したものを Fig.4.20 に示す。この関係を用いると、
ye =
Gs
ys
Ge
(4.10)
という関係が得られる。したがって、近似的に Gs ≈ Ge であれば ys ≈ ye となり、本実験においても
シミュレーションと同様に、近似的に閉ループ系の挙動を得ることが可能であると考えられる。もちろ
ん、Gs ≈ Ge の仮定はすべての周波数領域で成立するものではないが、実用的な範囲内では十分妥当な
近似であることが確認されている [32][33][36]。
4.5.3
実験結果
XZ 平面上での管の挙動を Fig.4.21 に、管の XY 平面における挙動を z = 500, 3500, 6500 mm の
地点毎に表したものを Fig.4.22 に、管先端部の X 軸 (制御) 方向の時系列結果を Fig.4.23 にそれぞれ示
す。なお、すべては結果は無次元値として表記している。Fig.4.24 には、管上端部に働く X, Y 軸方向
の力の時系列の結果を示す。また、すべての結果において、比較検討のために MmBR 制御の結果も同
時に示す。なお、Fig,4.21 は、マーカーのデータ間をスプライン補間し、ライザー全体の形状を表現し
ているが、実験を行う上で幾つかのトラブルが発生したために、本来 7 個あったマーカーのデータが、
実際には 4 個しか使用することができなかったために、若干、形状が複雑なものとなっている点に注意
していただきたい。
Fig.4.21 の結果から、MmBR 制御は LPV 制御よりも目標点近傍での大きな振動が確認される。こ
れは、数値シミュレーションにおいても確認されていたスピルオーバ現象による振動であると考えられ
る。実験においては、数値シミュレーションで考慮できなかった低次流体減衰効果などにより、大きな
振動は発生しないと考えていたが、先端部は予想以上に大きく揺れることが、Fig.4.23 の結果からも同
時に確認された。また、MmBR 制御によるリエントリー時の管形状は、LPV 制御のときに比べて大き
く変形することも Fig.4.21 の結果から確認される。
Fig.4.22 の結果からは、MmBR 制御は LPV 制御と比較して y 軸方向、つまり制御方向と垂直方向へ
の大きなふらつきが顕著に確認される。MmBR 制御はリエントリー開始時においても固定されたコン
トローラであるために、比較的に大きな速度でリエントリーを開始する。その大きな速度に反応して、
MmBR 制御は渦励振を誘起させてしまったものと考えられる。一方で、LPV 制御はリエントリー開始
時においては、速度 0 に対応したコントローラが働くことになるので、緩やかな速度で出発する。それ
ゆえ、MmBR 制御で見られたような渦励振を誘起しないために、リエントリー時に垂直方向のふらつ
きが見られなかったものと考えられる。
Fig.4.24 の結果は管先端部に働く力であるが、言い換えれば、船舶などの浮体部に要求される力であ
るため、必要な制御入力と解釈することもできる。そう考えると、LPV 制御はコントローラ可変のた
め、目標点に近づくにつれ徐々に入力が減少される一方、MmBR 制御は常に大きな制御入力を要求し
ており、消費エネルギーの観点から考えると大きな問題である。
4.6
まとめ
本章では、ライザー管に働く流体抗力が速度変動により大きく変化することに着目し、LPV 制御に
よるリエントリー作業の自動化に関する研究を行った。その結果、下記のような結論が得られた。
• MmBR 制御による方法は、ライザー管の上端速度が零になる領域、つまり、目標点近傍において
大きな振動現象を発生させる可能性がある。また、過渡期においても許容範囲以上の上部傾斜角
を引き起こしてしまい、リエントリー作業には適した方法とは考えにくい。
• LPV 制御による方法は、速度変動により変化する動特性の変化に対してコントローラも変化させ
61
るために、MmBR 制御による方法で見られたような振動現象や大きな上端部傾斜角を誘起させな
いことが、数値シミュレーションと水槽実験により確認された。
• 本研究で対象としたライザー管モデルは実機との対応については、残念ながらあまり考慮されて
いないが、浮体部の速度に応じて動特性が変化することは実機においても同じである。したがっ
て、本研究で検討した線形パラメータ変動制御方式は十分な有効性を持つ可能性を有するものと
考えられる。
本章の内容は、大水深ライザー管のリエントリー作業に対する LPV 制御の可能性を示したに過ぎず、
依然として数多くの問題が残されている。ゆえに、今後の課題として、
• 本実験においても十分な妥当性のあると思われる実験を行ったが、あくまでも開ループ系の実験
であるために、今後、厳密な閉ループ系の実験により検証を行う必要がある。
• ライザー管の運動は 2 次元運動でしかないために、まず、3 次元運動の動力学について明らかに
し、それに対して、本制御方式の有効性と問題点を検討する必要がある。
• 波浪や潮流などの環境外乱は依然として考慮されていないので、それらに対して十分な検討を行
わなければならない。
などが挙げられる。
62
x
Floating Body
(r,0)
u
X
θ
Flexible Riser
( µ,η)
z
Z
Fig. 4.1: Coordinate systems
Riser angle [rad]
Normalized top point
r
θ
0.04
0.1
0
with hydrodynamical damp.
without hydrodynamical damp.
0
0
1
0.04
2
Normalized time
with hydrodynamical damp.
without hydrodynamcal damp.
0
1
2
Normalized time
(a) Top point r
(b) Riser’s angle θ
Normalized bottom point
Normalized amplitude
W1 0.01
0.1
0
0
with hydrodynamical damp.
without hydrodynamical damp.
−0.01
0
1
0
2
with hydrodynamical damp.
without hydrodynamical damp.
1
Normalized time
(c) Amplitude W1
2
Normalized time
(d) Bottom point
Fig. 4.2: Simulation results of marine riser’s free motion (solid line: with hydrodynamical damping
(Cd = 1.17), dotted line: without hydrodynamical damping (Cd = 0.0))
63
|r/u|
|q/u|
200
100
0
-100
-200
200
0
-200
|W1/u|
-400
200
0
-200
-400
-3
10
10
-2
-1
0
10
1
10
2
10
3
10
10
Frequency (Hz)
Fig. 4.3: Bode plots of open-loop marine riser system (solid line: |ṙ| = 0, dash line: |ṙ| = 1.5)
30
20
Imag
10
0
-10
-20
-30
-150
-100
-50
0
Real
Fig. 4.4: Pole placement of open-loop marine riser system (×: |ṙ| = 0, ○: |ṙ| = 1.5)
w
Riser
System
contr.
-
y
z
WD s
Fig. 4.5: Control system structure
64
Design Space
Physical Space
Infinite
Dimensional
System
finite
Dimensional
System
Mode
Expansion
Linearization
Riser ‘s
LPV Model
Design Spec.
Validation
by Exp.
Design
Validation
by Sim.
LPV
Controller
Fig. 4.6: Concept of controller system design in reentry operation of riser
|r/w|
200
100
0
-100
-200
200
|q/w|
0
-200
|W1/w|
-400
200
0
-200
-400 -5
10
-4
10
-3
10
-2
-1
10
10
0
10
1
10
2
3
10
10
Frequency (Hz)
Fig. 4.7: Bode plots of closed-loop marine riser system (solid line: |ṙ| = 0, dash line: |ṙ| = 1.5)
30
20
Imag
10
0
-10
-20
-30
-150
-100
-50
0
Real
Fig. 4.8: Pole placement of closed-loop marine riser system (×: |ṙ| = 0, ○: |ṙ| = 1.5)
65
Riser angle
Normalized top point
r
θ 0.08
LPV Controller
MmBR Controller
0.1
LPV Controller
MmBR Controller
Hardware Limitted Angle
0.04
0
0
0
1
2
3
4
5
0
1
2
3
4
Normalized time
5
Normalized time
(a) Top point position
(b) Connected riser’s angle
Normalized bottom point
Normalized controller output
u
LPV Controller
MmBR Controller
0.1
LPV Controller
MmBR Controller
20
0
0
0
1
2
3
4
−20
5
0
1
2
3
4
5
Normalized time
Normalized time
(c) bottom point
(d) Controller output
Fig. 4.9: Comparison between MmBR control and LPV control in numerical simulation ((a):top
point position, (b):connected riser’s angle, (c):bottom point position, (d):controller output)
0
Water Depth
Water Depth
0
−0.5
−1
−0.5
−1
Start Point
−0.1
Target Point
−0.05
Sway Displacement
Start Point
0
−0.1
(a) MmBR Control
Target Point
−0.05
Sway Displacement
0
(b) LPV Control
Fig. 4.10: Simulation results of the marine riser’s sway motion in the reentry operation ((a):MmBR
control result, (b):LPV control result)
66
0.02
Normalized amplitude W 1
Normalized amplitude W 2
0.01
0.005
0
0
0.01
0
1
2
3
4
5
Normalized time
0.005
0
0
0
1
2
3
4
5
Normalized time
0
1
(c) 3rd mode amplitude W3
0.001
0.0005
0
0
0
1
2
3
4
5
Normalized time
0.0005
4
5
Normalized time
0
1
2
3
4
5
Normalized time
(f) 6th mode amplitude W6
Normalized amplitude W 7
Normalized amplitude W 8
0.0005
0.0005
0
0
0.0005
3
Normalized amplitude W 6
(e) 5th mode amplitude W5
0.001
2
(d) 4th mode amplitude W4
Normalized amplitude W 5
0.001
0.001
4
5
Normalized time
0.002
0
0.002
3
Normalized amplitude W 4
Normalized amplitude W 3
0.002
0.002
2
(b) 2nd mode amplitude W2
(a) 1st mode amplitude W1
0.004
1
0
1
2
3
4
5
Normalized time
0.0005
(g) 7th mode amplitude W7
0
1
2
3
4
5
Normalized time
(h) 8th mode amplitude W8
Fig. 4.11: Time histories of 1∼8 mode amplitudes in numerical simulation (solid line: LPV control
result, dotted line: MmBR control result)
67
RIser model
CCD camera frame
Fig. 4.12: Riser model in Research Institute for Applied Mechanics, Kyushu University
Fig. 4.13: Parallel mechanism to control the upper end of the riser model in Research Institute for
Applied Mechanics, Kyushu University
68
Upper End
0.20 m
0.30 m
0.10 m
1.00 m
1.00 m
1.00 m
6.60 m
1.00 m
1.00 m
1.00 m
Fig. 4.14: Arrangement of measuring marker
-Y
Camera Frame
Load Cell
-X
C1B
C1A
C2B
C2A
C3B
C3A
C4B
C4A
C5B
C5A
Z
C6B
C6A
C7B
C7A
Bottom Weight Flexible Riser
Fig. 4.15: Overview of CCD camera frame
69
y
200
5
6
x
9
7
8
200
200
200
1
4 100 100 3
200
100
13 100
CameraA
200
10
C1A
12
11
2
C2A
CameraB
C1B
C2B
Z
Caliblation Marker
Image by Camera1A
Image by Camera1B
Fig. 4.16: Calibration marker’s arrangement
Camera Coordinate System
X
yc
zc
Y
xc
u
Z
P
World Coordinate System
v
Image Coordinate System
Fig. 4.17: Camera coordinate systems
70
Parallel Mechanism
PC
Generate command
to control riser's position
Load Cell
CCD Camera
Universal
Joint
Video Data
Riser
Motion
Video Signal
-> Mpg file ->Tif file
Motion Analysis
using Halcon
Flexible Riser Model
PC
Fig. 4.18: Measurement flow of the motion
Start
Finish
0.30 [m]
Start Point
Target Point
Fig. 4.19: Experiment of reentry operation
71
Transfer Function in Simulation
rc
K
u
Transfer Function in Experiment
ys
Gs
u
K
rc
u=
1 + Gs K
rc
Ge
ye
K
rc
1 + Gs K
Ge
K
1 + Gs K
Ge
ys
Gs K
=
rc 1 + G s K
ye
ye
ye
Ge K
=
rc 1 + G s K
Fig. 4.20: Adequacy of this experiment
0
Water Depth Z/Z 0
Water Depth Z/Z 0
0
0.5
1
Start Point
1
1
Target Point
0.5
Sway Displacement X/X0
0.5
Start Point
1
0
Target Point
0.5
Sway Displacement X/X0
(a) MmBR control
(b) LPV control
Fig. 4.21: Experimental results of riser’s sway motion in the reentry operation
72
0
0.2
0.2
Start point
0
−0.1
−1
0
−0.1
Target point
−0.2
Start point
0.1
Y/X0
Y/X0
0.1
−0.5
X/X0
Target point
−0.2
0
−1
(a) z = 500 mm
0.2
Start point
0
−0.1
−1
0
Target point
−0.1
Target point
−0.2
Start point
0.1
Y/X0
0.1
Y/X0
0
(b) z = 500 mm
0.2
−0.5
X/X0
−0.2
0
−1
(a) z = 3500 mm
0.2
Start point
Y/X0
0.1
0
−0.1
−0.5
X/X0
0
−0.1
Target point
−1
0
(b) z = 3500 mm
Start point
−0.2
−0.5
X/X0
0.2
0.1
Y/X0
−0.5
X/X0
Target point
−0.2
0
−1
(a) z = 6500 mm
−0.5
X/X0
0
(c) z = 6500 mm
Fig. 4.22: Experimental results of XY-plane motion in the reentry operation ((a): MmBR control,
(b): LPV control)
0
X/X0
X/X0
0
−1
−1
0
25
50
75
0
25
50
Time (sec)
75
Time (sec)
(a) MmBR control
(b) LPV control
Fig. 4.23: Time histories of bottom point motion in control direction
1
FX (N)
FX (N)
1
0
−1
0
25
50
0
−1
75
0
25
Time (sec)
(a) Fx
50
75
1
FY (N)
FY (N)
75
(b) Fx
1
0
−1
50
Time (sec)
0
25
50
0
−1
75
Time (sec)
0
25
Time (sec)
(a) Fy
(b) Fy
Fig. 4.24: Time histories of top end forces Fx and Fy in the reentry operation ((a): MmBR control,
(b): LPV control)
73
第5章
結言
本論文は、線形パラメータ変動制御技術を「固定ピッチ型水平軸風車の可変速制御による高効率化」
と「大水深ライザー管の動特性を考慮したリエントリー作業の自動化」の 2 つの海洋工学分野の問題に
適用し、数値シミュレーションと実験によりその有効性を検討したものである。
本論文により得られた主な成果を以下に示す。
1. 風力発電機の運転モードのうち、モード 2 とモード 3 の 2 つのモードに着目し、それぞれモード
において効率的な運転をするための回転数制御問題を定式化した。
2. それぞれのモードについて、風力発電機の数学モデルを動作点周りで線形化を行い、風速 V をス
ケジューリング変数にもつ線形パラメータ変動モデルを導出した。モード 3 においては、出力係
数の変微分項 dCT /dλ の関係から、容易に線形パラメータ変動モデルと導出ことができないため、
本論文では出力係数 CT の二次関数近似を行い、風速 V と回転数 ω の 2 つのスケジューリング変
数を持つ 2 変数タイプ線形パラメータ変動モデルを提案した。
3. 導出した LPV コントローラの性能検証のために、数値シミュレーションによる検討を行ったと
ころ、モード 2 では、最適周速比 4.51 からの分散は従来的に用いられている ω 2 制御と比較して
1/5 程度になっており、出力においては 6% 程度の改善が確認された。モード 3 においては、出力
300W からの分散は従来 ω 2 制御の 1/4 程度となることが分かった。ゆえに、数値シミュレーショ
ンにより、両モードにおいての LPV 制御の優位性が確認された。
4. フィールド試験や風洞試験の限界などを考慮して、風力発電機開発支援装置「風車シミュレータ」
を提案し、その空力性能の実現のための計算フローを示した。
5. 風車シミュレータを用いた検証実験を行ったところ、モード 2 においては、LPV 制御の出力は
0.980、従来 ω 2 制御は 0.60 であり約 1.5 倍もの出力を LPV 制御は獲得することが確認された。一
方で、モード 3 においては、LPV 制御の 300W からの分散は 50.324、従来 ω 2 制御は 36241.780
と約 1/600 の改善が確認された。ゆえに、風車シミュレータを用いた実験において、両モードに
おいて LPV 制御の優位性が確認された。
6. 大水深ライザー管の運動方程式を Hamilton の原理から導出し、モード展開による有次元化、平
衡点周りでの線形化を行うことにより、管上端部の速度をスケジューリング変数に持つ線形パラ
メータ変動モデルを導出した。
7. 導出した線形パラメータ変動モデルの極配置、ボード線図を調査することにより、修正モリソン
式で表現した流体抗力により、管上端部の速度が大きくなるにつれてシステムは安定化し、逆に
速度が小さくなるにつれてシステムの安定性は失われることを明らかにした。
8. 導出した LPV コントローラの性能検証のために、数値シミュレーションによる検討を行ったとこ
ろ、MmBR 制御は上端部最大傾斜角が 3.783 (deg)、最大コントローラ出力が 19.247 (N) である
一方、LPV 制御は、上端部最大傾斜角 1.777 (deg)、最大コントローラ出力 12.433 (N) と大きく
74
改善されることが分かった。また、MmBR 制御では、目標点近傍において振動的な現象が確認さ
れ、LPV 制御ではそのような現象を誘起させないことも確認された。
9. 九州大学応用力学研究所において水槽実験を行ったところ、実験装置の関係から厳密な閉ループ系
実験を行うことはできなかったが、十分妥当性のある開ループ系実験により、LPV 制御は MmBR
制御と比較しても優位な制御手法であることを確認した。
75
付 録A
Nomoto モデルに基づく船舶の回頭角制御について
A.1
はじめに
船舶の操縦性は、船速や積載荷物などの条件によって変化する。船速が速ければ舵の効果は大きく、逆
に船速が遅ければ小さくなる。ゆえに、操船の現場では、船速や積載荷重などの条件による操縦性の変化
に応じて、PID 補償器のゲイン再設定が一般的に行われている。そのような観点から、船舶の運動制御
に関する研究分野では、適応制御についての研究が主に行われてきた。例えば、最近では D.C.Donha に
よる H∞ 適応制御による方法 [42]、S.Nejim によるゲインセルフチューニング法による方法 [43] などが
挙げられる。
本章では、船舶の運動を表す基本的なモデルとして、野本により提案された 1 次モデル [44][45](以後、
Nomoto モデルと略記) に着目し、船速により大きく変化する動特性をコントローラのスケジューリン
グ、また、Nomoto モデルと非線形運動方程式の間のギャップをロバスト性によりカバーするという考
えの下、船舶に対する LPV 制御の可能性について検討した [46]。
A.2
船舶の非線形操縦運動方程式
本章で対象とする船舶は、文献 [45] の Mariner Class Vessel である。Table A.1 にその主要目を示す。
Fig.A.1 に示す座標系において、前後 (surge)、横流れ (sway)、旋回 (yaw) 運動の 3 自由度の船舶の操
縦運動を取り扱う。その運動は、次式に示す無次元化操縦運動方程式により示される。ただし、∆ はノ
ミナル成分からの摂動成分を意味する。

0
0
0

= ∆X 0
 (m − Xu̇ )∆u̇
0
0
0
0
0
0
0
(A.1)
(m − Yv̇ )∆v̇ + (m xG − Yṙ )∆ṙ = ∆Y 0

 (m0 − x0 )∆v̇ 0 + (I 0 − N 0 )∆ṙ0
0
= ∆N
ṙ
G
Z
∆X 0 , ∆Y 0 , ∆N 0 は流体力の項であり、次のように記述される。流力微係数については、文献 [45] を
参考にしていただきたい。以後、表記の簡単化のために、摂動成分を意味する ∆ は省略する。
0
0
0
0
0
0
0
∆X 0 = Xu0 ∆u0 + Xuu
∆u02 + Xuuu
∆u03 + Xvv
∆v 02 + Xrr
∆r02 + Xrv
∆r0 ∆v 0 + Xδδ
∆δ 02 + Xuδδ
∆u0 ∆δ 02
0
0
+ Xvδ
∆v 0 ∆δ 0 + Xuvδ
∆u0 ∆v 0 ∆δ 0
0
0
0
0
∆Y 0 = Yv0 ∆v 0 + Yr0 ∆r0 + Yvvv
∆v 03 + Yvvr
∆v 02 ∆r0 + Yvu
∆v 0 ∆u0 + Yru
∆r0 ∆u0 + Yδ0 ∆δ 0 + Yδδδ ∆δ 03
0
0
0
+ Yuδ ∆u0 ∆δ 0 + Yuuδ
∆u02 ∆δ 0 + Yvδδ
∆v 0 δ 02 + Yvvδ
∆v 02 ∆δ 0 + (Y 00 + Y 00u ∆u0 + Y 00uu ∆u02 )
0
0
0
0
0
∆N 0 = Nv0 ∆v 0 + Nr0 ∆r0 + Nvvv
∆v 03 + Nvvr
∆v 02 r0 + Nvu
∆v 0 ∆u0 + Nru
∆r0 ∆u0 + Nδ0 ∆δ 0 + Nδδδ
∆δ 03
0
0
0
0
+ Nuδ
∆u0 ∆δ 0 + Nuuδ
∆u02 ∆δ 0 + Nvδδ
∆v 0 ∆δ 02 + Nvvδ
∆v 02 ∆δ 0 + (N 00 + N 00u ∆u0 + N 00uu ∆u02 )
ここで、無次元化は次に従うものとする。
u, v 0 rL
X, Y
xG
, u0 , v 0 =
, r =
, X 0, Y 0 = 1 2 2
L
U
U
2 ρU L
N
m
I
z
N 0 = 1 2 3 , m0 = 1 3 , Iz0 = 1 5
2 ρU L
2 ρL
2 ρL
x0G =
76
L, U, ρ
u0 , v 0
r 0 , m0
IZ0
X 0, Y 0
N0
δ0
:
:
:
:
:
:
:
船長、前進速度、及び、流体の密度
船舶の前後速度、横流れ速度の無次元値
回頭角速度、船体質量の無次元値
z 軸まわりの慣性モーメントの無次元値
作用する x, y 軸方向の流体力の無次元値
重心まわりの回頭モーメントの無次元値
舵角の無次元値
Table A.1: Main data of Mariner Class Vessel
Length overall (Loa )
Maximum beam (B)
Design draft (T )
Design displacement (∇)
Design speed (U ∗ )
171.80
23.17
8.23
18541
7.72
(m)
(m)
(m)
(m3 )
(m/s)
舵の性能限界として、舵角と舵角変化率には (A.2) 式のような制約があるために、Fig.A.2 に示すよ
うな構成の下で制御問題を取り扱うことにする。
|δ| ≤ 10 (deg),
A.3
|δ̇| ≤ 5 (deg/s)
(A.2)
Nomoto モデルの動特性とステップ応答解析
船舶の操縦性などを評価する上で、Nomoto モデルは非常に重要である [44][45]。また、そのモデル
の簡潔さから、Nomoto モデルに基づいた制御系設計は、低次元コントローラを導出することを可能に
する。そこで本章では、(A.1) 式の非線形操縦運動方程式に対応する Nomoto モデルを確認し、そのモ
デルに基づいた制御系設計を行う。
まず、はじめに、(A.1) 式の Nomoto モデルを導出する。前後方向の運動変化は非常に小さい (u ' 0)、
及び、流体力の 2 次以上の高次の項を無視すると、次のような線形モデルが導出される。
#
#"
"
m0 x0G − Yṙ0
v̇ 0
m0 − Yv̇0
=
Iz0 − Nṙ0
m0 x0G − Nv̇0
ṙ0
|
{z
} | {z }
#"
#
" E
# ẋ "
#
"
0
0
0
0
Yv̇ Yṙ
v
Y
0
Yδ0
+
δ0 +
(A.3)
N 00
Nv̇0 Nṙ0
r0
Nδ0 |{z}
u
|
{z
} | {z } | {z }
| {z }
x
A
ここで、舵角 δ から回頭角速度
する。
r0
B
d
までの伝達関数 ĝ(s) を考えるために、出力方程式を次のように定義
#
0
v
y = r0 = 0 1
0
| {z } r
| {z }
C
h
i
"
(A.4)
x
このとき、(A.1) 式に対応する Nomoto モデルは次のように記述される。
ṙ0 (t) = −
1 0
K0 0
r
(t)
+
δ (t)
T0
T0
(A.5)
ただし、定常ゲイン K 0 と時定数 T 0 は (A.6) 式と (A.7) 式のように表される。
K0 =
Nv0 Yδ0 − Yv0 Nδ0
Yv0 Nr0 − Yr0 Nv0
77
(A.6)
T0 =
−(m0 − Yv̇0 )Nr0 − (Iz0 − Nṙ0 )Yv0 + (m0 x0G − Yṙ0 )Nv0 + (m0 x0G − Nv̇ )Yr0
Yv0 Nr0 − Yr0 Nv0
(m0 − Yv̇0 )Nδ0 − (m0 x0G − Nv̇0 )Yδ0
−
Nv0 Yδ0 − Yv0 Nδ0
(A.7)
なお、本章で対象とした Mariner Class Vessel の定常ゲインと時定数は、K = 0.1851 と T = 163.559 (s)
である。時間領域における Nomoto モデルの特徴を把握するために、ステップ応答の様子を Fig.A.3 と
Fig.A.4 に示す。Fig.A.3 は舵を 5 (deg)、Fig.A.4 は 10 (deg) にしたときの応答の様子を表している。比
較のために、(A.1) 式の非線形操縦運動方程式のステップ応答の結果も同時に示す。舵角を大きく取る
ほど、その応答特性は、Nomoto モデルと非線形操縦運動方程式とでは異なっていく様子が分かる。つ
まり、舵を大きくとるほど、非線形項の影響を強く受けることを意味している。しかしながら、Nomoto
モデルは非線形操縦運動方程式の運動の特徴を十分に表しているものと考えられる。
T 0 , K 0 は無次元値であるので有次元化すると、
U
L 0
T , K = K0
U
L
(A.8)
L 0
U∗ 0
∗
T
,
K
=
K
U∗
L
(A.9)
T =
また、U を設計速度 U ∗ とすると、
T∗ =
ここで、(A.8) 式と (A.9) 式から L を消去すると
T (U ) =
U∗ ∗
U
T , K(U ) = ∗ K ∗
U
U
(A.10)
となる。∗ は設計速度に応じた値であることを意味する。つまり、船速 U が大きくなるほど時定数は小
さくなり、定常ゲインは大きくなることが分かる。(A.10) 式を (A.5) 式に代入すると、次のようなモデ
ルが導出される。
µ ¶2 ∗
U
K
U 1
δ(t)
(A.11)
ṙ(t) = − ∗ ∗ r(t) +
∗
U T
U
T∗
本来、Nomoto モデルは、船速一定時における動作点周りにおける船舶の動特性を表すモデルである。
しかしながら、船舶のように緩やかに速度が変化するシステムにおいて、各船速時において Nomoto モ
デルが瞬時瞬時に成立すると考えることが可能である。つまり、船舶がゆるやかに変化していく下では、
軌跡上に各船速に対応する Nomoto モデルが離散的に存在すると近似的に考えることができる。その考
え方を Fig.A.6 に示す。それを連続の繋がりとして表現したものが (A.11) 式であると考えられる。そこ
で、本章では、(A.11) 式を有次元 Nomoto モデルと呼ぶことにする。これを状態方程式として行列表
記すると次のようになる。
d
dt
"
r(t)
Ψ(t)
#
"
#"
# "
#
−U/(U ∗ T ∗ ) 0
r(t)
(U/U ∗ )2 K ∗ /T ∗
=
+
δ(t)
|{z}
1
0
Ψ(t)
0
|
{z
} | {z } |
{z
} u
x
A(U )
(A.12)
B(U 2 )
ここで、ψ は回頭角を表す。船速 U の有次元 Nomoto モデルに及ぼす影響を確認するために、Fig.A.5
にボード線図を示す。実線は U = U ∗ 、破線は U = 1.5U ∗ 、点線は U = 0.5U ∗ のときの特性である。U
の変化に応じて、特性が大きく変化することが分かる。
有次元 Nomoto モデルの特徴は、船速の 2 乗で舵の効果が変化するという特性を陽に反映している点
である。したがって、この特徴をコントローラに反映させるために、船舶の制御系設計時においては、
有次元 Nomoto モデルに基づいた制御系設計を行うことが一つの手段であると考えられる。この有次元
Nomoto モデルに基づく制御系設計の概念について、Fig.A.7 に示す。
78
A.4
船舶の PID 制御則と GSPID 制御則の問題点
T.I.Fossen は船舶の PID 制御則について次のように述べている [45]。
Z t
δ(t) = Kp (ψc − ψ(t)) − Kd r(t) + Ki
(ψc − ψ(τ ))dτ
(A.13)
0
ただし、
T ωn2
2T ζωn − 1
ωn
, Kd =
, Ki =
Kp
(A.14)
K
K
10
ここで、ζ と ωn は閉ループ系の減衰係数と固有振動数であるが、一般的に次の条件を満足する中で設
計者が適切な値を決定する。
q
p
1
< ωn 1 − 2ζ 2 + 4ζ 4 − 4ζ 2 + 2 < ωδ
(A.15)
T
Kp =
右辺の ωδ はアクチュエータの帯域幅を意味する。また、減衰係数については、0.8 ≤ ζ ≤ 1.0 の間で選
択される。
一方で、Fig.A.5 に示すように、船速によって定常ゲインや時定数が変化する船舶に対して、PID 制
御則を適用しようとした場合、制御性能の劣化が予想される。そこで、(A.10) 式を (A.14) 式に代入する
と次のように、PID 制御則を拡張させたものとして、ゲインスケジュールド PID(Gain-scheduled、以
後、GSPID と略記) 制御則を構成することが可能である。この方法も、T.I.Fossen の文献の中で述べら
れている [45]。
Z t
δ(t) = Kp (U )(ψc − ψ(t)) − Kd (U )r(t) + Ki (U ) (ψc − ψ(τ ))dτ
(A.16)
0
ただし、

µ ∗¶ ∗
T 2
U


ω
Kp (U ) =


∗ n

U |K{z

}



kp



µ ∗ ¶2 ∗
µ ∗¶


U
U
T
1

Kd (U ) =
−
2ζω
n
∗
∗
U
K
U
K
{z
}
|
|{z}



kd1
kd2


µ ∗ ¶2 ∗ 3



U
T ωn


Ki (U ) =

∗

U


|K {z10}

(A.17)
ki
PID 制御則、GSPID 制御則の効果を確認するために、(A.11) 式の有次元 Nomoto モデル、(A.1) 式
の非線形操縦運動方程式に対しての数値シミュレーションによる検証を行った。有次元 Nomoto モデル
に対しての検証を行った理由は、線形システムと非線形システムに対する制御性能の特徴を把握するた
めである。シミュレーションを実施する際には、障害物回避などの際に要求されるような、減速しなが
ら回頭角制御をするような状況を想定し、速度変動パターンを Fig.A.11 のように与えた。また、PID
ゲイン選定の際には、船速の変化がないときの Nomoto モデルに対して、望ましい応答が得られるゲイ
ン (kp = 2.2096, kd1 = 62.4980, kd2 = 5.4039, ki = 0.0110) を選択した。case 1, 2, 3 に対応する数値シ
ミュレーション結果を Fig.A.8∼A.10 に示す。なお、この時の目標回頭角は ψc = 10 (deg) と設定して
いる。また、参考までに、速度変動がない場合の結果も同時に示している。有次元 Nomoto モデルに対
する一連の結果を見ると、速度変動パターンが case1 から case3 へと緩やかになるにつれて、振動的で
あった回頭角は、速やかに収束するようになる様子が確認される。一方、非線形操縦運動方程式に対す
る結果は、速度変動が緩やかになっても、その制御性能が良くなる様子は見られない。また、非線形操
縦運動方程式においては船速変化がなくても、PID 制御則、GSPID 制御則は有効な方法ではないこと
も確認された。よって、非線形操縦運動方程式に対しては、船速変化の仕方に関係無く、PID 制御則、
GSPID 制御則には限界があるということが言える。
79
A.5
有次元 Nomoto モデルに基づく船舶の LPV 制御系設計
本章では、(A.12) 式の有次元 Nomoto モデルに対して、状態フィードバック型の LPV コントローラ
を設計する。その詳細について述べる前に、線形行列不等式 (Linear Matrix Inequality, LMI) ベースの
LPV 制御系設計法について説明するため、次の一般化プラントについて考える。
(
ẋ(t) = A(θ)x(t) + B1 (θ)w(t) + B2 (θ)u(t)
(A.18)
z(t) = C1 (θ)xG (t) + D11 (θ)w(t) + D12 (θ)u(t)
ここで、A(θ) ∈ Rn1 ×n1 , D11 (θ) ∈ Rp1 ×m1 , D12 (θ) ∈ Rp2 ×m2 であり、θ ∈ Rq はスケジューリング変
数である。この一般化プラントに対して、θ に応じて変化する極配置制約付きの H∞ 状態フィードバック
u = −F (θ)x
(A.19)
を求める。
まず、閉ループ系が H∞ 性能を満足するには、有界実補題と等価な次の LMI を満足する正定行列
P > 0 と Yi が存在することと等価である [13]。

T
B1i
Ai P − B2i Y + P ATi − YiT B2i

T
B1i
−γI

C1i P − D12i Yi
D11i

T
T
T
P C1i − Yi D12i

T
(A.20)
D11i
 < 0, (i = 1, 2, · · · , n)
−γI
下付添字 i は、スケジューリング変数の端点での値であることを意味する。同様にして極配置条件は、
コニックセクタ安定条件、アルファ安定条件、円安定条件の三つの LMI を満足する共通解を求めるこ
とにより、Fig.2.6 に示すような扇形の極配置領域を指定することが可能である [13]。この制御仕様を課
すことにより、サンプリング周期以上の速いダイナミクスを有するコントローラを回避する。(A.21) 式
はアルファ安定条件、(A.22) 式は円安定条件、(A.23) 式はコニックセクタ安定条件を表す LMI である。
T
Ai P − B1i Yi + P ATi − YiT B1i
+ 2αP < 0
"
#
−rP
Ai P − Bi Yi
<0
P ATi − YiT BiT
−rP
³
´
³
´ 
cos θ Ai P − Bi Yi − P ATi + YiT BiT
sin θ Ai P − Bi Yi + P ATi − YiT BiT
³
´
³
´ <0

− cos θ Ai P − Bi Yi − P ATi + YiT BiT
sin θ Ai P − Bi Yi + P ATi − YiT BiT
(A.21)
(A.22)

(A.23)
スケジューリング変数の変動領域を包含する端点それぞれに (A.20),(A.21),(A.22),(A.23) 式の LMI を
構成し、それらを同時に γ を最小化する最適化問題を解き、正定行列 P > 0 と Yi を求めることにより、
次のように状態フィードバックゲインを求める。
Fi = Yi P −1 ,
(i = 1, 2, · · · , n)
(A.24)
その結果、次のように端点コントローラに重み pi (θ) を掛け、Fig.2.12 に示すように θ と共にスケジュー
ルさせることにより LPV コントローラを設計することができる。
F (θ) =
n
X
pi (θ)Fi , ただし
i=1
n
X
i=1
80
pi (θ) = 1
(A.25)
本章では、有次元 Nomoto モデルに対して LPV コントローラを設計する際、混合感度問題の枠組み
で考えるために、Fig.A.12 に示すような制御システム構造の下で制御系設計を行うこととする。なお、
フィルターのダイナミクスとして gf (s) = 1/(2s + 1) を選択した。この一般化プラントは次の (A.26) 式
のように表される。Wd は微分重みであり、p = [rc , ψc ]T は参照信号、z は評価信号、下付添字 f はフィ
ルターを意味する。



B(U 2 )


ẋ
A(U )
02×2 02×1 

 

=
Af
01×2 Bf  
 ẋf   01×2

z
Wd A(U ) Wd B(U 2 ) 02×2 02×1
x
xf
p
u





(A.26)
(A.26) 式に対して、(A.20),(A.21),(A.22),(A.23) 式の LMI を各端点毎に作成して、それらを同時に解
くことにより端点コントローラを求めた。その際、端点の選択として、Fig.A.13 に示すような 3 つの点
を考えることにする。本制御系設計において、端点 1 は (5.79, 33.52)、端点 2 は (9.65, 93.12)、端点 3 は
(8.11, 58.71) を選択した。また、端点コントローラを船速 U に応じてスケジュールさせる際の (A.25)
式の重み pi (U ) の詳細については、3 章の設計法を参考にしていただきたい。
A.6
数値シミュレーション
本章では、前節まで述べてきた有次元 Nomoto モデルに基づく LPV 制御の可能性を確認するために、
数値シミュレーションによる検証を行った。計算条件は前節のものと同じである。また、PID 制御則、
GSPID 制御則の結果と比較するために、前節で示した結果と重複するが示している。また、基本的に、
状態フィードバック型 LPV コントローラは、次数などの制御構造が PID 制御則や GSPID 制御則と異
なるため、本章では、公平な比較を行うために、整定時間が約 200 (s) になるように、LPV コントロー
ラを設計した。アルファ安定条件を 0.012、円安定条件を 0.8、コニックセクタ安定条件を π/3 (rad) と
して設計した。case 1, 2, 3 に対応する数値シミュレーション結果を Fig.A.14∼A.16 に示す。
PID 制御則や GSPID 制御則においては、急激な速度変動があるときは、非線形操縦運動方程式はも
ちろんのこと、有次元 Nomoto モデルに対しても十分な制御性能を実現することができていない一方、
LPV 制御の結果は、PID 制御則や GSPID 制御則で整定するまでの間に観察されたような大きな揺れ
を発生させることなく、速やかな収束を実現している様子が確認される。また、LPV 制御の性能は速
度変動パターンにほとんど影響されない。sin 関数などの速度変動の場合についても検討したが、同様
の結果が得られた。これは、端点モデルに対して同時安定化を行い、時変系に対する安定性を LPV コ
ントローラが保証していることが、この結果の要因ではないかと考えている。
A.7
まとめ
本章では、極配置制約付き状態フィードバック LPV 制御を理解するために、船舶の回頭角制御に応
用した。その際、船舶の運動を記述する Nomoto モデルが船速の 2 乗に大きく依存することに着目し、
有次元 Nomoto モデルに基づいた LPV 制御系設計法を提案した。その結果、下記のような結論が得ら
れた。
• PID 制御則や GSPID 制御則は、速度変動下においては、非線形操縦運動方程式に対してはもち
ろんのこと、有次元 Nomoto モデルに対しても振動的な過渡現象を発生させる傾向にある。有次
元 Nomoto モデルに対しては、速度変動が急激になるほどその傾向は顕著になるが、非線形操縦
運動方程式に対しては、速度変動がない場合においても、その振動現象が発生する場合がある。
• 有次元 Nomoto モデルに基づく LPV 制御による方法は、対象とする非線形操縦運動方程式と
Nomoto モデル間に大きなギャップがあるにも関わらず、PID 制御則で見られたような過渡期の
81
振動現象を発生させない有効な制御手法である。また、その制御性能は速度変動パターンにほと
んど影響されない。
本研究は、有次元 Nomoto モデルに基づく船舶の LPV 制御についての基礎的研究段階であるために、
まだ、一つの可能性を示しているにすぎない。したがって、今後の課題として、
• 風波などの環境外乱に対する本制御方式のロバスト性について検討する。
• バルクキャリアやタンカーなどの大型船舶などの進路不安定性を有する船舶に対してのオートパ
イロット時の本制御方式の有効性について検討する。
などが挙げられる。
82
X
Body-fixed Axis
ψ
surge
x
u
U
yaw
r
G
v
δ
y
sway
Y
Fig. A.1: Coordinate systems
Controller
δ*
δ
Ship motion
Rudder
Saturation
Rudder
Rate limitter
Fig. A.2: Block diagram of ship control
yaw angle (deg)
100
Nonlinear model
50
Nomoto model
0
0
100
200
time (s)
Fig. A.3: Step response of ship δc = 5 [deg]
yaw angle (deg)
100
Nomoto model
50
Nonlinear model
0
0
100
200
time (s)
Fig. A.4: Step response of ship δc = 10 [deg]
83
Magnitude (dB)
-10
U=1.5U*
-20
U=U*
-30
-40
-50
U=0.5U*
-60
-70
Phase (deg)
U=1.5U*
U=U*
U=0.5U*
10
-4
10
-3
10
-2
10
-1
10
0
10
1
Frequency (rad/s)
Fig. A.5: Bode diagram of Nomoto’s model varying with ship forward speed
Nomoto’s Model
2
U 3 1 ( ) U3 K* ( )
r t +
δ t
U* T*
U* T *
( )
r(t) =
Ship Trajectory
G
U=U3
Nomoto’s Model
( )
G
r(t) =
2
U 1 1 ( ) U1 K* ( )
δ t
r t +
U* T*
U* T *
U = U2
G
U = U1
Nomoto’s Model
r(t) =
2
U 2 1 ( ) U2 K* ( )
r t +
δ t
U* T*
U* T *
( )
Fig. A.6: Concept of dimensional Nomoto’s model
Physical Space
Nonlinear
Dynamics
Design Space
Extraction
Reduction
Design Spec.
Validation
Dimentional
Nomoto
Model
Design
Controller
Fig. A.7: Concept of controller system design based on dimensional Nomoto’s model
84
15
10
PID
GSPID
PID (No Speed Variation)
5
0
0
200
yaw angle (deg)
yaw angle (deg)
15
10
PID
GSPID
PID (No Speed Variation)
5
0
400
0
time(s)
0
PID
GSPID
PID (No Speed Variation)
200
(b) Yaw angle (deg)
rudder angle (deg)
rudder angle (deg)
10
0
400
time(s)
(a) Yaw angle (deg)
−10
200
10
0
PID
GSPID
PID (No Speed Variation)
−10
400
0
time(s)
200
400
time(s)
(a) Rudder angle (deg)
(b) Rudder angle (deg)
Fig. A.8: Time histories of ship steering motion in the dimensional Nomoto’s and the nonlinear
model simulation ((a):dimensional Nomoto’s model, (b):nonlinear model, speed variation: case 1)
15
10
PID
GSPID
PID (No Speed Variation)
5
0
0
200
yaw angle (deg)
yaw angle (deg)
15
10
PID
GSPID
PID (No Speed Variation)
5
0
400
0
time(s)
0
PID
GSPID
PID (No Speed Variation)
200
(b) Yaw angle (deg)
rudder angle (deg)
rudder angle (deg)
10
0
400
time(s)
(a) Yaw angle (deg)
−10
200
400
10
0
PID
GSPID
PID (No Speed Variation)
−10
0
time(s)
200
400
time(s)
(a) Rudder angle (deg)
(b) Rudder angle (deg)
Fig. A.9: Time histories of ship steering motion in the dimensional Nomoto’s and the nonlinear
model simulation ((a):dimensional Nomoto’s model, (b):nonlinear model, speed variation: case 2)
85
15
yaw angle (deg)
yaw angle (deg)
15
10
PID
GSPID
PID (No Speed Variation)
5
0
0
200
10
PID
GSPID
PID (No Speed Variation)
5
0
400
0
200
time(s)
(b) Yaw angle (deg)
rudder angle (deg)
rudder angle (deg)
(a) Yaw angle (deg)
10
0
PID
GSPID
PID (No Speed Variation)
−10
0
200
400
time(s)
10
0
PID
GSPID
PID (No Speed Variation)
−10
400
0
200
time(s)
400
time(s)
(a) Rudder angle (deg)
(b) Rudder angle (deg)
Design Speed (m/s)
Fig. A.10: Time histories of ship steering motion in the dimensional Nomoto’s and the nonlinear
model simulation ((a):dimensional Nomoto’s model, (b):nonlinear model, speed variation: case 3)
10
Case 1
Case 2
Case 3
5
0
0
200
400
time(s)
Fig. A.11: Speed variation pattern
Augmented system
WD s
p
Nomoto
Model
Filter
xf
u
contr.
y
Fig. A.12: Control system structure
86
z
Vertex2
2
Umax
U
2
Vertex3
2
Umin
Vertex1
Umin
Umax
U
Fig. A.13: Triangle region under the varying ship speed
15
yaw angle (deg)
yaw angle (deg)
15
10
PID
GSPID
LPV
5
0
0
200
10
PID
GSPID
LPV
5
0
400
0
time(s)
0
PID
GSPID
200
LPV
rudder angle (deg)
rudder angle (deg)
(b) Yaw angle (deg)
10
0
400
time(s)
(a) Yaw angle (deg)
−10
200
400
10
0
PID
GSPID
−10
0
time(s)
200
LPV
400
time(s)
(a) Rudder angle (deg)
(b) Rudder angle (deg)
Fig. A.14: Time histories of ship steering motion in the dimensional Nomoto’s and the nonlinear
model simulation ((a):dimensional Nomoto’s model, (b):nonlinear model, speed variation: case 1)
87
15
yaw angle (deg)
yaw angle (deg)
15
10
PID
GSPID
LPV
5
0
0
200
10
PID
GSPID
LPV
5
0
400
0
200
time(s)
(b) Yaw angle (deg)
10
0
PID
GSPID
0
200
LPV
rudder angle (deg)
rudder angle (deg)
(a) Yaw angle (deg)
−10
400
time(s)
10
0
PID
GSPID
−10
400
0
time(s)
200
LPV
400
time(s)
(a) Rudder angle (deg)
(b) Rudder angle (deg)
Fig. A.15: Time histories of ship steering motion in the dimensional Nomoto’s and the nonlinear
model simulation ((a):dimensional Nomoto’s model, (b):nonlinear model, speed variation: case 2)
15
yaw angle (deg)
yaw angle (deg)
15
10
PID
GSPID
LPV
5
0
0
200
10
PID
GSPID
LPV
5
0
400
0
time(s)
0
PID
GSPID
200
LPV
rudder angle (deg)
rudder angle (deg)
(b) Yaw angle (deg)
10
0
400
time(s)
(a) Yaw angle (deg)
−10
200
400
10
0
PID
GSPID
−10
0
time(s)
200
LPV
400
time(s)
(a) Rudder angle (deg)
(b) Rudder angle (deg)
Fig. A.16: Time histories of ship steering motion in the dimensional Nomoto’s and the nonlinear
model simulation ((a):dimensional Nomoto’s model, (b):nonlinear model, speed variation: case 3)
88
付 録B
大水深ライザー管の運動方程式の導出
本章では大水深ライザー管の運動方程式を導出する。はじめに、Fig.4.1 に示すような 2 つの座標系
を定義する。1 つは空間固定座標系であり、X, Z 軸を持つ。また、空間固定座標系において [r(t), 0] を
原点とし、これから θ 回転した座標系を管固定座標系とする。管固定座標系における管の変形後の位置
は [η(z, t), z + µ(z, t)] で表される。η(z, t), µ(z, t) は管の撓みと伸びを意味する。[r(t), 0] は管上端部位
置を表す。空間固定座標系における管の中心線の位置は以下のように表される。
(
X = η cos θ + (z + µ) sin θ + r
(B.1)
Z = −η sin θ + (z + µ) cos θ
速度成分は (B.1) 式を用いて以下のように表される。
(
Ẋ = η̇ cos θ + µ̇ sin θ + {−η sin θ + (z + µ) cos θ}θ̇ + ṙ
Ż
=
−η̇ sin θ + µ̇ cos θ − {η cos θ + (z + µ) sin θ}θ̇
(B.2)
同様にして、(B.1) 式を用いて加速度成分は

Ẍ = η̈ cos θ + µ̈ sin θ + {−η sin θ + (z + µ) cos θ}θ̈ + r̈





+2(−η̇ sin θ + µ̇ cos θ)θ̇ − {η cos θ + (z + µ) sin θ}θ̇2

Z̈ = −η̈ sin θ + µ̈ cos θ − {η cos θ + (z + µ) sin θ}θ̈




−2(η̇ cos θ + µ̇ sin θ)θ̇ + {η sin θ − (z + µ) cos θ}θ̇2
(B.3)
と表される。次に、Hamilton の原理から大水深ライザー管の運動方程式を導出するために、運動エネ
ルギー T 、ポテンシャルエネルギー U 、仕事エネルギー W をそれぞれ求める。
B.1
運動エネルギー (T )
ここでは運動エネルギー T を求める。管の同一断面上の点は同じ速度を持つと仮定すると、全長 l の
管の運動エネルギーは以下のように表される。
T
=
=
Z
Z
´
1 l ³ 2
1
2
V dm =
m̃ Ẋ + Ż 2 dz
2 m
2 0
Z lh
1
m̃
η̇ 2 + µ̇2 + θ̇2 {η 2 + (z + µ)2 } + ṙ02 − 2θ̇µ̇η + 2θ̇η̇(z + µ)
2
0
i
+ 2ṙ0 {η̇ cos θ + µ̇ sin θ + (−η sin θ + (z + µ) cos θ)θ̇} dz
89
(B.4)
ただし、
m̃
ρt
ρo
ρi
At
Ao
Ai
B.2
:
:
:
:
:
:
:
管の単位長さあたりの水中質量 (= ρt At + ρi Ai − ρo Ao )
管の密度
外部流体の密度 (通常は海水)
内部流体の密度
管の断面積 (= Ao − Ai )
外径までの面積
内径までの面積
ポテンシャルエネルギー (U )
ここではポテンシャルエネルギー U を求める。ポテンシャルエネルギーとしては、位置エネルギー Ugrav 、
伸びによる歪みエネルギー Uaxial 、曲げによる歪みエネルギー Ubend 、張力による歪みエネルギー Utens
が挙げられる。以後、それらについて説明する。
B.2.1
位置エネルギー (Ugrav )
管の密度を一定と考えると、位置エネルギーは以下のように表される。g は重力加速度を表す。
Z
Ugrav = −
B.2.2
Z
l
l
m̃gZdz = −
m̃g{−η sin θ + (z + µ) cos θ}dz
0
(B.5)
0
軸方向の変形 (伸び) によるひずみエネルギー (Uaxial )
長さ dz の要素に軸方向力 P が働いた時の変形を dµ とすると、応力-ひずみ関係は次のようになる。
σ=
P
,
At
ε=
dµ
= µ0 ,
dz
σ = Eε
(B.6)
ここで、E はヤング率である。(B.6) 式より、この要素の軸方向伸縮によるひずみエネルギー蓄積は以
下のように表される。
1
1 dµ
1
2
dUaxial = P dµ = P dz = EAt µ0 dz
(B.7)
2
2 dz
2
これを長さ方向に積分すると、軸方向変形による全体のひずみエネルギーが得られる。
Z
Uaxial =
B.2.3
0
l
dUaxial
1
=
2
Z
l
0
2
EAt µ0 dz
(B.8)
曲げによるひずみエネルギー (Ubend )
曲げによる歪みエネルギーは以下のように表される。
Ubend =
1
2
Z
l
M dψ =
0
1
2
Z
l
M
0
dψ
dz
dz
(B.9)
ここで考える管は線形弾性であり Hooke の法則に従うと仮定する。この時、曲げモーメントは以下の
方程式に従う。
M
1
(B.10)
κ= =−
ρ
EI
ここで、κ は曲率、I は断面2次モーメントである。
90
(B.10) 式中の曲率は
κ=
1
dψ
=
ρ
ds
(B.11)
と表される。また、回転角 ψ は
dη
tan ψ ∼
(B.12)
=ψ=
dz
と得られる。ψ は十分に小さいので cos ψ ∼
= 1 と近似できて、曲げモーメントによる変形のひずみエネ
ルギー (B.9) 式は (B.10)、(B.11) 及び (B.12) 式を用いて以下のように表される。
Ubend =
B.2.4
1
2
Z
0
l
1
EIκdψ ∼
=
2
Z
µ
l
EI
0
dψ
dz
¶2
dz =
1
2
Z
µ
l
EI
0
d2 η
dz 2
¶2
dz =
1
2
Z
l
2
EIη 00 dz
(B.13)
0
張力によるひずみエネルギー (Utens )
この内部張力 Te による横方向変形のひずみエネルギーは次のように表される。
Z
1 l
Utens =
Te sin ψdr
2 0
(B.14)
回転角 ψ は十分小さいので
dη ∼ dη ∼ 0
sin ψ ∼
(B.15)
=
=η
=
ds
dz
と仮定できる。従って、(B.14) 式は (B.15) 式を用いて以下のように表される。
Z
Z
1 l ¡ 0¢ 0
1 l
2
Utens =
Te r r dz =
Te η 0 dz
(B.16)
2 0
2 0
³
z´
張力は t と z の関数であるが簡単のために、Te = m̃gl 1 −
+ WB とする。WB は下端荷重を表す。
l
B.3
仮想仕事
ここでは、仮想仕事を求める。仮想仕事成分としては、流体力による仮想仕事 (Wf luid )、上端部に働
く荷重による仮想仕事 (Wupper ) が挙げられる。以後、それらについて説明する。
B.3.1
流体力による仮想仕事 (Wf luid )
流体力 Q は z 軸に垂直に働くと考え、z 軸方向には働かないものとする。ある断面に働く流体力の
X, Z 軸方向成分を QX , QZ とする。Q = [QX , QZ ]T は次のように与える。
¯ ¯
1
¯ ¯
Q(z, t) = −ρo Ao CA Ẍ n − ρo Cd Do ¯Ẋ n ¯ Ẋ n
(B.17)
2
ただし、
CA
Cd
D0
an
:
:
:
:
付加質量係数
抗力係数
円柱の直径
あるベクトル a の管の中心線に対する法線ベクトル
である。(B.17) 式において、流体抗力の推定には Morison 式を採用した。空間固定座標における、流体
力による仕事は以下の式によって求められる。
Z l
Z l
{QX δX + QZ δZ}dz
Q · δXdz =
(B.18)
δWf luid =
0
0
91
任意の法線ベクトル an は
"
#
"
#
aX
sin θ
n
a =
− (aX sin θ + aZ cos θ)
aZ
cos θ
"
# "
#
#"
aX
sin2 θ
sin θ cos θ
aX
=
−
aZ
sin θ cos θ
cos2 θ
aZ
"
#"
#
cos2 θ
− sin θ cos θ
aX
=
2
− sin θ cos θ
sin θ
aZ
#
"
# "
aX cos2 θ − aZ sin θ cos θ
cos θ(aX cos θ − aZ sin θ)
=
=
−aX sin θ cos θ + aZ sin2 θ
− sin θ(aX cos θ − aZ sin θ)
(B.19)
のように得られるので、ゆえに、管の運動速度の法線ベクトル Ẋ n は以下のように表される。
"
#
"
#
³
´
n
o
cos
θ
cos
θ
Ẋ n = Ẋ cos θ − Ż sin θ
= η̇ + (z + µ)θ̇ + ṙ cos θ
(B.20)
− sin θ
− sin θ
n = η̇ + (z + µ)θ̇ + ṙ cos θ とする。
ここで、Vrel
また、同様にして、加速度の法線ベクトルは以下のように表される。
#
#
"
"
n
o
³
´
cos
θ
cos
θ
= η̈ + (z + µ)θ̈ + 2µ̇θ̇ − η θ̇2 + Ẍ0 cos θ
Ẍ n = Ẍ cos θ − Z̈ sin θ
(B.21)
− sin θ
− sin θ
したがって、(B.20) 式と (B.21) 式を用いて、流体力ベクトル (B.17) 式は以下のように表される。

n
o
2 + r̈ cos θ cos θ

Q
(z,
t)
=
−ρ
A
C
η̈
+
(z
+
µ)
θ̈
+
2
µ̇
θ̇
−
η
θ̇
 X
o o A


n
o


1
n | η̇ + (z + µ)θ̇ + ṙ cos θ cos θ = −Q cos θ − Q cos θ = −Q cos θ

− 2 ρo Cd Do |Vrel
0
A
D
n
o
2

QZ (z, t) = ρo Ao CA η̈ + (z + µ)θ̈ + 2µ̇θ̇ − η θ̇ + r̈ cos θ sin θ



n
o


n | η̇ + (z + µ)θ̇ + ṙ cos θ sin θ = Q sin θ + Q sin θ = Q sin θ

+ 21 ρo Cd Do |Vrel
0
A
D
(B.22)
ここで、mA は単位長さあたりの付加質量 (= ρo Ao CA ) であり、
n
o

2 + r̈ cos θ

Q
=
m
η̈
+
(z
+
µ)
θ̈
+
2
µ̇
θ̇
−
η
θ̇
A
A


n
o
1
n | η̇ + (z + µ)θ̇ + ṙ cos θ
ρ
C
D
|V
Q
=
o
o
D
d
rel
2



Q0 = QA + QD
(B.23)
である。X, Z の変分 δX, δZ は (B.1) 式を用いて、
(
δX = δ{η cos θ + (z + µ) sin θ + r} = δη cos θ + δµ sin θ + {−η sin θ + (z + µ) cos θ}δθ + δr
δZ = δ{−η sin θ + (z + µ) cos θ} = −δη sin θ + δµ cos θ − δθ{η cos θ + (z + µ) sin θ}
(B.24)
と表されるので、(B.22) 式と (B.24) 式を用いると、流体力 (B.18) 式による仮想仕事は以下のように表
92
される。
Z
δWf luid =
l
(QX δX + QZ δZ) dz
0
Z
l
[(QX cos θ − QZ sin θ)δη + (QX sin θ + QZ cos θ)δµ
=
0
Z
=
0
l
+ {QX (−η sin θ + (z + µ) cos θ) − QZ (η cos θ + (z + µ) sin θ)}δθ + QX δr] dz
£
Q0 −(cos2 θ + sin2 θ)δη + (− sin θ cos θ + cos θ sin θ)δµ
Z
= −
0
B.3.2
+ {− cos θ(−η sin θ + (z + µ) cos θ) − sin θ(η cos θ + (z + µ) sin θ)}δθ − cos θδr] dz
l
Q0 {δη + (z + µ)δθ + cos θδr} dz
(B.25)
上端の荷重による仮想仕事 (Wupper )
管上端部から r 方向に力 u が働く際、仮想仕事は以下のように書くことができる。
Z
δW2 =
0
l
δD (z)uδXdz
(B.26)
ここで δD (z) は Dirac デルタ関数とする。最終的に、外力 u による仮想仕事は、式 (B.26) を用いて次
のように得られる。
δW2 = uδr
B.4
(B.27)
各種エネルギーの変分の導出
ここでは、運動エネルギー (T )、ポテンシャルエネルギー (U ) の変分成分を導出する。
B.4.1
運動エネルギーの変分 (δT )
運動エネルギーの変分は以下のように表される。
Z
Z
1
2
δT =
δV dm =
(Ẋδ Ẋ + Żδ Ż)dm
2 m
m
(B.28)
ここで、δK1 = Ẋδ Ẋ + Żδ Ż とする。また、速度成分の変分は以下のように得られる。


δ Ẋ = δ[η̇ cos θ + µ̇ sin θ + {−η sin θ + (z + µ) cos θ}θ̇ + ṙ]






= cos θδ η̇ + sin θδ µ̇ + δ ṙ + {−η sin θ + (z + µ) cos θ}δ θ̇ − θ̇ sin θδη







+θ̇ cos θδµ + [−η̇ sin θ + µ̇ cos θ − {η cos θ + (z + µ) sin θ}θ̇]δθ


δ Ż












(B.29)
= δ[−η̇ sin θ + µ̇ cos θ − {η cos θ + (z + µ) sin θ}θ̇]
= − sin θδ η̇ + cos θδ µ̇ − {η cos θ + (z + µ) sin θ}δ θ̇
−θ̇ cos θδη − θ̇ sin θδµ − [η̇ cos θ + µ̇ sin θ − {η sin θ − (z + µ) cos θ}θ̇]δθ
93
従って、
δK1 = ṙδ ṙ + (Ẋ cos θ − Ż sin θ) δ η̇ + (Ẋ sin θ + Ż cos θ) δ µ̇
|
{z
}
|
{z
}
Aη
Aµ
+ [Ẋ{−η sin θ + (z + µ) cos θ} − Ż{η cos θ + (z + µ) sin θ}] δ θ̇
|
{z
}
Aθ
+ (−Ẋ sin θ − Ż cos θ)θ̇ δη + (Ẋ cos θ − Ż sin θ)θ̇ δµ
|
{z
}
|
{z
}
Bη
Bµ
+ [Ẋ{−η̇ sin θ + µ̇ cos θ − (η cos θ + (z + µ) sin θ)θ̇} + Ż{−η̇ cos θ − µ̇ sin θ + (η sin θ − (z + µ) cos θ)θ̇}] δθ
|
{z
}
Bθ
これを空間と時間の両方で積分すると
Z t1Z l
Z l ¯
¯t1
¯
¯
m̃δK1 dzdt =
m̃ ¯Aη δη + Aµ δµ + Ẋδr + Aθ δθ¯ dz
0
0 0
0
Z t1Z l
+
m̃{(Bη − Ȧη )δη + (Bµ − Ȧµ )δµ + (Bθ − Ȧθ )δθ − Ẍδr}dzdt
0
ここで、
(B.30)
(B.31)
0


Ȧη = Ẍ cos θ − Z̈ sin θ − θ̇(Ẋ sin θ + Ż cos θ)






Ȧµ = Ẍ sin θ + Z̈ cos θ + θ̇(Ẋ cos θ − Ż sin θ)



Ȧθ = Ẍ{−η sin θ + (z + µ) cos θ} − Z̈{η cos θ + (z + µ) sin θ}





+Ẋ{−η̇ sin θ + µ̇ cos θ − (η cos θ + (z + µ) sin θ)θ̇}





+Ż{−η̇ cos θ − µ̇ sin θ + (η sin θ − (z + µ) cos θ)θ̇}
(B.32)
これらを用いると

Bη − Ȧη = −(Ẍ cos θ − Z̈ sin θ)




n
o


2

=
−
η̈
+
(z
+
µ)
θ̈
+
r̈
cos
θ
+
2
µ̇
θ̇
−
η
θ̇








Bµ − Ȧµ = −(Ẍ sin θ + Z̈ cos θ)



n
o




= − µ̈ − η θ̈ + r̈ sin θ − 2η̇ θ̇ − (z + µ)θ̇2




Bθ − Ȧθ = −[Ẍ{−η sin θ + (z + µ) cos θ} − Z̈{η cos θ + (z + µ) sin θ}]




= η(Ẍ sin θ + Z̈ cos θ) − (z + µ)(Ẍ cos θ − Z̈ sin θ)




n
o


2

= η µ̈ − η θ̈ + r̈ sin θ − 2η̇ θ̇ − (z + µ)θ̇




n
o



2

−(z
+
µ)
η̈
+
(z
+
µ)
θ̈
+
r̈
cos
θ
+
2
µ̇
θ̇
−
η
θ̇




h
i



= − {η 2 + (z + µ)2 }θ̈ + (z + µ)η̈ − η µ̈ − {η sin θ − (z + µ) cos θ}r̈ + 2{(z + µ)µ̇θ̇ + η η̇ θ̇}
(B.33)
B.4.2
ポテンシャルエネルギーの変分 (δU )
(B.5), (B.7), (B.13), (B.16) 式を用いて、ポテンシャルエネルギーの変分は
δU
= δUaxial + δUbend + δUtens + δUgrav
Z l
{EAt µ0 δµ0 + EIη 00 δη 00 + Te η 0 δη 0 − m̃gδZ}dz
=
0
94
(B.34)
と表される。これを部分積分すると以下のようになる。
δU
¯l
¯
= ¯EAt µ0 δµ + EIη 00 δη 0 + {−(EIη 00 )0 + Te η 0 }δη ¯0
Z l
− [{−(EIη 00 )00 + (Te η 0 )0 − m̃g sin θ}δη + (EAt µ00 + m̃g cos θ)δµ
0
−m̃g{η cos θ + (z + µ) sin θ}δθ}δη]dz
B.5
(B.35)
Hamilton の原理による運動方程式の導出
上記までの結果を整理すると、Hamilton の原理より次のような恒等式が導出される。
Z t1
Z t1
[δL + δW ]dt =
[δT − δU + δW ]dt
0
Z
0
¯
¯t1
¯
¯
=
m̃ ¯Aη δη + Aµ δµ + Aθ δθ + Ẋδr¯ dz
0
0
Z t1
¯
¯
¯EAt µ0 δµ + EIη 00 δη 0 + {−(EIη 00 )0 + Te η 0 }δη ¯l dt
−
0
0
Z t1Z l hn
o
−
−m(Bη − Ȧη ) + (EIη 00 )00 − (Te η 0 )0 + m̃g sin θ + Q0 δη
0 0
n
o
+ −m(Bµ − Ȧµ ) − EAt µ00 − m̃g cos θ δµ
n
o i
+ −m(Bθ − Ȧθ ) + m̃g(η cos θ + (z + µ) sin θ) + (z + µ)Q0 δθ dzdt
¾
Z t1 ½Z l
−
(m̃Ẍ + cos θQ0 )dz + u δrdt
0
0
Z t1
¯
¯
¯EAt µ0 δµ + EIη 00 δη 0 + {−(EIη 00 )0 + Te η 0 }δη ¯l dt
= −
0
0
Z t1Z l hn
´
³
2
−
(m̃ + mA ) η̈ + (z + µ)θ̈ + Ẍ0 cos θ + 2µ̇θ̇ − η θ̇
0 0
ª
+(EIη 00 )00 − (Te η 0 )0 + m̃g sin θ + QD δη
´
o
n ³
+ m̃ µ̈ − η θ̈ + r̈ sin θ − 2η̇ θ̇ − (z + µ)θ̇2 − EAt µ00 − m̃g cos θ δµ
n ³
+ m̃ (η 2 + (z + µ)2 )θ̈ + (z + µ)η̈ − η µ̈ − (η sin θ − (z + µ) cos θ)r̈
´
+2((z + µ)µ̇θ̇ + η η̇ θ̇)
l
+ m̃g(η cos θ + (z + µ) sin θ) + (z + µ)QD
´o i
³
+mA (z + µ) η̈ + (z + µ)θ̈ + r̈ cos θ + 2µ̇θ̇ − η θ̇2 δθ dzdt
Z t1 ·Z l n ³
−
m̃ η̈ cos θ + µ̈ sin θ + (−η sin θ + (z + µ) cos θ)θ̈ + r̈
0
0
´
+2(−η̇ sin θ + µ̇ cos θ)θ̇ − (η cos θ + (z + µ) sin θ)θ̇2
³
´
o
i
+mA cos θ η̈ + (z + µ)θ̈ + r̈ cos θ + 2µ̇θ̇ − η θ̇2 + cos θQD dz − u δrdt = 0
(B.36)
ゆえに、大水深ライザー管の運動方程式は次のように記述される。
95
水平方向
Z lh
n
o
(m̃ + mA ) cos θ η̈ + (z + µ)θ̈ + r̈ cos θ + 2µ̇θ̇ − η θ̇2
0
n
o
+m̃ sin θ µ̈ − η θ̈ + r̈ sin θ − 2η̇ θ̇ − (z + µ)θ̇2
¯
¯n
o¸
1
¯
¯
+ ρo Cd Do cos θ ¯η̇ + (z + µ)θ̇ + ṙ cos θ¯ η̇ + (z + µ)θ̇ + ṙ cos θ dz = u
2
(B.37)
回転方向
Z lh
n
o
(m̃ + mA )(z + µ) η̈ + (z + µ)θ̈ + r̈ cos θ + 2µ̇θ̇ − η θ̇2
0
n
o
−m̃η µ̈ − η θ̈ + r̈ sin θ − 2η̇ θ̇ − (z + µ)θ̇2 + m̃g {η cos θ + (z + µ) sin θ}
¯
¯n
o¸
1
¯
¯
+ ρo Cd Do (z + µ) ¯η̇ + (z + µ)θ̇ + ṙ cos θ¯ η̇ + (z + µ)θ̇ + ṙ cos θ dz = 0
2
(B.38)
たわみ方向
n
o
(m̃ + mA ) η̈ + (z + µ)θ̈ + r̈ cos θ + 2µ̇θ̇ − η θ̇2 + (EIη 00 )00 − (Te η 0 )0 + m̃g sin θ
¯
¯n
o
1
¯
¯
+ ρo Cd Do ¯η̇ + (z + µ)θ̇ + ṙ cos θ¯ η̇ + (z + µ)θ̇ + ṙ cos θ = 0
2
(B.39)
伸び方向
o
n
m̃ µ̈ − η θ̈ + r̈ sin θ − 2η̇ θ̇ − (z + µ)θ̇2 − EAt µ00 − m̃g cos θ = 0
境界条件
η = µ = 0 (at z = 0)
EAt µ0 = 0
EIη 00
(−EIη 000
=0
(at z = 0, l)
+ Te
(B.40)
(at z = l)
η0) η
= 0 (at z = l)
上記の境界条件において、1 番目の式は上端が単純支持されている事による変位無しの条件、2 番目は
下端が伸び方向に関しての自由端条件。3 番目は両端においてモーメントが 0 を表し、最後は下端を自
由端と考えた場合、せん断力 0 となることを表している。z 軸を両端を結ぶように採った場合、最後の
境界条件は η(l, t) = 0 となる。
96
参考文献
[1] H.Kajiwara, R.Gao, K.Ohtsubo, W.Kotereyama, and M Nakamura. LPV control simulation of
an underwater vehicle using virtual reality technology. In Proceedings of International Conference
on Marine Simulation and Ship Manueverability, Aug 2003.
[2] R.Gao, K.Ohtsubo, and H.Kajiwara. LPV design for a space vehicle attitude control benchmark
problem. In Proceedings of SICE Annual Conference, Aug 2003.
[3] 小玉成人. 風力発電機の出力変動抑制に関する研究: 八戸工業大学博士論文, 2002.
[4] H.Kajiwara, P.Apkarian, and P.Gahinet. LPV techniques for control of an inverted pendulum.
IEEE Control Systems Magazine, Vol. 19, No. 1, pp. 44–54, Feb 1999.
[5] J.S.Shamma and J.R.Clotier. Gain-scheduled missile autopilot design using linear parameter
varying transformations. Journal of Guidance, Control and Dynamics, Vol. 16, No. 2, pp. 256–
263, 1993.
[6] R.Gao and H.Kajiwara. Attitude tracking control using inner/outer loop LPV techniques. In
Proceedings of 4th Asian-Pacific Conference on Aerospace Technology and Science, Nov 2002.
[7] R.Gao, E.Kondo, H.Kajiwara, W.Koterayama, and M.Nakamura. Depth control of an underwater vehicle using linear parameter-varying techniques. International Journal of Offshore and
Polar Engineering, Vol. 13, No. 1, pp. 52–59, Mar 2003.
[8] R.Gao. Linear Parameter-Varying Control of an Underwater Vehicle. PhD thesis, Kyushu
University, 2002.
[9] P.Gahinet, A.Nemirovski, A.J.Laub, and M.Chilali.
MATLAB-. The Math Works, Inc., 1995.
LMI Control Toolbox -For Use with
[10] 岩崎徹也. LMI と制御. 昭晃堂, 1997.
[11] 藤森篤. ロバスト制御. コロナ社, 2001.
[12] M.Chilali and P.Gahinet. H∞ design with pole placement constraints: An LMI approach. IEEE
Transactions on Automatic Control, Vol. 41, No. 3, pp. 358–367, Mar 1996.
[13] P.Apkarian, G.Becker, P.Gahinet, and H.Kajiwara. LMI techniques in control engineering from
theory to practice. In Workshop Notes CDC, Dec 1996.
[14] SICE セミナーテキスト. ゲインスケジューリング制御の基礎と応用. (社) 計測自動制御学会, 2002.
[15] 大坪和久, 本田大作, 梶原宏之. 固定ピッチプロペラ型風力発電機の可変速 LPV 制御. 日本船舶海
洋工学会論文集, 第 1 号, pp. 1–8, 2005.
[16] 大坪和久, 本田大作, 梶原宏之. 風力発電機の最適回転数制御. 日本造船学会講演会論文集, 第 3 号,
pp. 173–174, 2004.
97
[17] 大坪和久, 梶原宏之, 桜井晃, 大屋裕二, 鳥谷隆. 風レンズ風車に対する LPV 制御系設計. 風力エネ
ルギー利用シンポジウム, 2004.
[18] K.Ohtsubo and H.Kajiwara. LPV rotational speed control of wind turbine using measured wind
speed. In Proceedings of OCEANS’04,MTS/IEEE/TECHNO-OCEANS’04, Nov 2004.
[19] 牛山泉. 風車工学入門-基礎理論から風力発電技術まで-. 森北出版株式会社, 2002.
[20] 加藤和彦. NEDO における風力発電技術の開発現状と課題. 風力エネルギー利用シンポジウム,
2002.
[21] R.Datta and V.T.Ranganathan. A method of tracking the peak power points for a variable speed
wind energy conversion system. IEEE Transactions on Energy Conversion, Vol. 18, No. 1, pp.
163–168, 2003.
[22] 涌井徹也, 橋詰匠, 斉藤久男, 長尾利夫. 自立電源用風力発電システムの変速制御運転に関する研究.
風力エネルギー利用シンポジウム, 2002.
[23] W.E.Leithead, D.J.Leith, F.Hardan, and H.Markou. Global gain-scheduling control for variable
speed wind turbines. In Proceedings of European Wind Energy Conference EWEC, 1999.
[24] H.D.Battista, P.F.Puleston, R.J.Mantz, and C.F.Christiansen. Energy-based approach to the
output feedback control of wind energy systems. International Journal of Control, Vol. 76,
No. 3, pp. 299–308, 2003.
[25] 大屋祐二, 鳥谷隆, 桜井晃, 井上雅弘. 風レンズ効果 (風エネルギーの集中) による風力発電の高出力
化. 風力エネルギー利用シンポジウム, 2001.
[26] K.E.Johnson, L.J.Fingersh, M.J.Balas, and L.Y.Pao. Methods for increasing region 2 power
capture on a variable-speed wind turbine. Transactions of the ASME, Vol. 126, pp. 1092–1100,
2004.
[27] 大坪和久, 五百木陵行, 梶原宏之. 速度変動に伴う流体抗力の影響を考慮した大水深ライザー管の
ゲインスケジューリング制御. 日本船舶海洋工学会論文集, 第 2 号, pp. 1–8, 2006.
[28] 大坪和久, 千賀英敬, 眞鍋崇寛, 小寺山亘, 梶原宏之. 大水深ライザー管のゲインスケジューリング
制御によるリエントリー実験について. 日本船舶海洋工学会論文集, 第 2 号, pp. 49–55, 2006.
[29] 大坪和久, 五百木陵行, 梶原宏之. 速度変動に伴う流体抗力の影響を考慮した大水深ライザー管の
ゲインスケジューリング制御. 西部造船会第 110 回例会論文梗概, pp. 37–42, 2005.
[30] 大坪和久, 千賀英敬, 眞鍋崇寛, 小寺山亘, 梶原宏之. 大水深ライザー管のゲインスケジューリング
制御の水槽実験について. 日本船舶海洋工学会講演会論文集, pp. 379–380, 2005.
[31] 五百木陵行, 大坪和久, 梶原宏之. IDA-PBC 法によるライザー管のリエントリー問題に対する非線
形制御. 西部造船会第 110 回例会論文梗概, pp. 43–48, 2005.
[32] W.Koterayama and H.Senga. Vortex induced vibration of a long riser oscillating regularly and
irregularly at its upper end. In Proceedings of International Symposium on Technology of Ultra
Deep Ocean Engineering, Feb 2005.
[33] H.Senga and W.Koterayama. An experimental and numerical study on vortex-induced vibrations
of a hanging flexible riser with its top in irregular motion. International Journal of Offshore and
Polar Engineering, Vol. 15, No. 4, pp. 274–281, Mar 2005.
98
[34] 鈴木英之, 吉田宏一郎, 南東浩, 村井基彦, 宇佐美陽生, 石田成幹. アクティブ制御による大水深ライ
ザーのリエントリー問題に関する研究. 日本造船学会論文集, 第 174 号, pp. 865–874, 1993.
[35] M.P.Fard. Modeling and Control of Mechanical Flexible Systems. PhD thesis, Norwegian University of Science and Technology, 2001.
[36] Y.P.Hong. A Study on Dynamics of a Flexible Marine Riser. PhD thesis, Kyushu University,
2004.
[37] 五百木陵行. 海洋で用いる長大弾性管の受動性に基づく非線形制御系設計に関する研究:九州大学
大学院工学府海洋システム工学専攻修士論文, 2005.
[38] R.D.Blevins. Flow-Induced Vibration. Van Nostrand Reinhold Company, 1990.
[39] T.Sarpkaya and M.Isaacson. Mechanics of Wave Forces on Offshore Structures. Van Nostrand
Reinhold Company, 1981.
[40] W.Koterayama. Wave forces acting on a vertical circular cylinder with a constant forward
velocity. Ocean Engineering, Vol. 11, No. 4, pp. 363–379, 1984.
[41] 鈴木英之, 渡辺啓介. 大水深ライザーの3次元傾斜角制御法の開発及び実験による検証. 日本造船
学会論文集, 第 188 号, pp. 335–342, 2000.
[42] D.C.Donha, D.S.Desanj, M.R.Katebi, and M.J.Grimble. H∞ adaptive controllers for auto-pilot
applications. International Journal of Adaptive Ccontrol and Signal Processing, Vol. 12, No. 8,
pp. 623–648, 2003.
[43] S.Nejim. Design of limited authority adaptive ship steering autopilots. International Journal of
Adaptive Ccontrol and Signal Processing, Vol. 14, No. 4, pp. 381–391, 2000.
[44] K.Nomoto, T.Taguchi, K.Honda, and S.Hirano. On the steering qualities of ships, technical
report. International Shipbuilding Progress, Vol. 4, , 1957.
[45] T.I.Fossen. Guidance and Control of Ocean Vehicles. Wiley, UK, 1994.
[46] 梶原宏之, 大坪和久. 速度変動を伴う Nomoto モデルに対する LPV 制御系設計について. 日本造船
学会講演会論文集, 第 2 号, pp. 149–150, 2003.
99
謝辞
本論文をまとめるにあたり、終始ご指導ご鞭撻を頂きました九州大学大学院工学研究院海洋システム
工学部門教授 梶原宏之先生には心より御礼申し上げます。
九州大学大学院工学研究院航空宇宙工学部門教授 桜井晃先生、九州大学大学院工学研究院海洋シス
テム工学部門教授 吉川孝男先生、九州大学応用力学研究所海洋大気力学部門助教授 中村昌彦先生には
本論文をまとめるにあたり、貴重なご助言、ご指導および暖かい激励を頂きました。ここに改めて心か
ら御礼申し上げます。
第 3 章の「固定ピッチ型水平軸風車の高効率化に関する研究」では、九州大学風レンズ研究会の皆様、
特に、九州大学応用力学研究所海洋大気力学部門教授 大屋裕二先生、同部門助教授 鳥谷隆先生、佐世
保工業高等専門学校校長 井上雅弘先生には、風車研究に関する数多くのご助言を頂きました。ここに
改めて心から御礼申し上げます。
第 4 章の「ライザー管リエントリー作業の自動化に関する研究」では、九州大学応用力学研究所海洋
大気力学部門教授 小寺山亘先生には、適切なご助言と実験を行う機会を与えて下さいました。また、九
州大学応用力学研究所技術員 稲田勝氏、及び九州大学大学院総合理工学府大気海洋環境システム学専
攻海中機器制御研究室の学生の皆様には、実験補助などの作業を手伝って頂きました。ここに改めて心
から御礼申し上げます。
九州大学大学院工学研究院海洋システム工学部門助教授 木村元先生、同部門技術員 木下順二氏、佐
賀大学海洋エネルギー研究センター助教授 豊田和隆先生、ならびに海洋システム工学部門教職員の方々
には多くのご支援を頂きました。ここに改めて心から御礼申し上げます。
九州大学大学院海洋システム工学専攻システム計画学研究室の卒業生、及び在校生の皆様、特に高鋭
氏 (現在: 株式会社トヨタテクノサービス)、本田大作氏 (現在: トヨタ自動車株式会社)、中島一誠氏 (現
在: トヨタ自動車株式会社)、五百木陵行氏 (現在: 博士後期課程 1 年)、眞鍋崇寛氏 (現在: 修士課程 2 年)
とは、寝食を共にしながら議論することにより、研究を深めていくことが出来ました。心より御礼申し
上げます。
最後に、今まで経済的、精神的に応援してくれた家族に心から感謝したいと思います。
100
ダウンロード

主論文 - 海上技術安全研究所