基于遺傳算法的半潛式平臺(tái)動(dòng)力定位系統(tǒng)動(dòng)態(tài)約束可行域推力分配法

梁海志，喬?hào)|生，李蘆鈺，張紀(jì)剛 ，歐進(jìn)萍

（1.大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；2.青島理工大學(xué) 土木工程學(xué)院，山東 青島266033）

0 引 言

動(dòng)力定位系統(tǒng)最早在鉆井船等大型船舶上得到廣泛應(yīng)用，隨著海洋石油開(kāi)采業(yè)向深海進(jìn)軍，動(dòng)力定位系統(tǒng)以其不受作業(yè)水深影響等優(yōu)勢(shì)逐漸取代傳統(tǒng)錨泊系統(tǒng)被應(yīng)用到深海海洋平臺(tái)上。在風(fēng)浪流聯(lián)合作用下，動(dòng)力定位推力系統(tǒng)提供控制力使半潛式平臺(tái)位移響應(yīng)保持在安全范圍內(nèi)，考慮到推力器出力方式及運(yùn)行經(jīng)濟(jì)性，僅控制縱蕩、橫蕩和艏搖三個(gè)自由度的低頻慢漂運(yùn)動(dòng)，而推力系統(tǒng)由8個(gè)推力器組成以保證系統(tǒng)的冗余度，因此整個(gè)系統(tǒng)為過(guò)驅(qū)動(dòng)系統(tǒng)。在動(dòng)力定位模塊化設(shè)計(jì)過(guò)程中，如圖1所示，首先由控制算法得到三自由度控制力，然后通過(guò)分配算法將其分配到底層執(zhí)行機(jī)構(gòu)推力器上，最終實(shí)現(xiàn)定位目標(biāo)。

動(dòng)力定位系統(tǒng)的底層推力系統(tǒng)一般采用全回轉(zhuǎn)推力器、舵等，在推力器方位角已知的情況

下，S?rdalen[1]采用偽逆法對(duì)推力直接進(jìn)行分配，分配過(guò)程中存在的奇異構(gòu)型問(wèn)題，通過(guò)修改奇異值來(lái)解決。付海軍等[2]提出了基于自適應(yīng)遺傳算法的加權(quán)偽逆策略，改善推進(jìn)器角度飽和問(wèn)題。偽逆法的特點(diǎn)是快捷有效，但是無(wú)法給出最優(yōu)推力分配方案。針對(duì)分配過(guò)程中的奇異構(gòu)型問(wèn)題，徐勝文等[3]提出了一種高效地避免奇異結(jié)構(gòu)的簡(jiǎn)單方法。序列二次規(guī)劃法是推力分配策略中采用較為普遍的優(yōu)化方法，如文獻(xiàn)[4-6]等。推力分配過(guò)程中，影響動(dòng)力定位系統(tǒng)性能及其運(yùn)行經(jīng)濟(jì)性的主要因素為推力器間尾流相互干擾造成的推力損失，針對(duì)該問(wèn)題，文獻(xiàn)[4-8]提出了對(duì)推力器推力方向設(shè)置推力禁區(qū)的方法，使推力器間的相互影響降到最低，但由此得到的推力可行域?yàn)榉峭辜?，無(wú)法直接采用序列二次規(guī)劃法，處理方法較為復(fù)雜。

本文針對(duì)以上存在的問(wèn)題，基于文獻(xiàn)[5]建立了以推力變化率和轉(zhuǎn)角變化速率為控制參數(shù)，以耗能最小為目標(biāo)的優(yōu)化模型。在考慮推力器推力損失的基礎(chǔ)上通過(guò)動(dòng)態(tài)設(shè)置推力可行域，利用遺傳算法對(duì)控制力進(jìn)行分配，在提高運(yùn)行經(jīng)濟(jì)性的同時(shí)也避免了序列二次規(guī)劃法對(duì)可行域要求嚴(yán)格的問(wèn)題。

圖1 動(dòng)力定位系統(tǒng)工作流程Fig.1 Workflow diagram of DPS

1 推力分配模型

為使動(dòng)力定位系統(tǒng)能準(zhǔn)確執(zhí)行頂層控制策略得到的控制力，推力器推力和控制力之間需滿(mǎn)足如下關(guān)系

式中，τ包括縱蕩、橫蕩控制力和艏搖控制力矩，T為推力器推力，T=［T1…Tm］T，m為底層推力系統(tǒng)推力器個(gè)數(shù)，B=［B（ α1）…B（ αm）］， ［lxi,lyi］為第i個(gè)推力器在圖2所示隨體坐標(biāo)下的布置坐標(biāo)，αi為第i個(gè)推力器轉(zhuǎn)角，Ti為第i個(gè)推力器的推力值，如圖3所示。

推力分配策略在滿(mǎn)足控制力要求的前提下，更應(yīng)使整個(gè)系統(tǒng)燃油消耗最小以提高運(yùn)行經(jīng)濟(jì)性，系統(tǒng)耗能功率為

式中：KQ、KT分別為螺旋槳扭矩系數(shù)和推力系數(shù)，ρ為海水密度，D為螺旋槳直徑。針對(duì)推力分配系統(tǒng)，還需進(jìn)一步考慮推力方向變化速率，由此得到動(dòng)力定位系統(tǒng)推力分配策略的優(yōu)化目標(biāo)

圖2 推力器布置示意圖Fig.2 Propeller layout diagram

圖3 推力器推力及推力方向示意圖Fig.3 Propeller’s thrust and thrust direction

約束條件：

（3）式中，T∈Rm×1，Δα∈Rm×1為優(yōu)化參數(shù)；Ω 為推力方向變化權(quán)矩陣，Ω∈Rm×m。 （3）、（4）式建立了以推力T和推力方向變化Δα為變量的推力分配基本數(shù)學(xué)模型，但是對(duì)于推力器而言，受其自身動(dòng)力性能限制，推力器很難實(shí)時(shí)響應(yīng)由上述模型得到的推力T，若以推力變化ΔT作為優(yōu)化變量將更滿(mǎn)足客觀

現(xiàn)實(shí)，為此采用ΔT替換推力T對(duì)上述模型進(jìn)行改進(jìn)，對(duì)（2）式系統(tǒng)耗能功率進(jìn)行泰勒展開(kāi)，可得

式中：Ti0為t0時(shí)刻第i號(hào)推力器推力，將上式代入（3）式，則得到以推力變化ΔT和推力方向變化Δα為變量的優(yōu)化目標(biāo)方程

同時(shí)，約束條件式(4)中的等式約束方程相應(yīng)的線性化為

推力變化ΔT和轉(zhuǎn)角變化Δα約束條件為

由此得到如圖4所示的推力器運(yùn)行可行域，并可通過(guò)α=α0+Δα，T=T0+ΔT得到推力器出力方向和推力，最后還需考慮推力幅值約束條件

式中：Tmin、Tmax為推力器的推力上下限。

圖4 單個(gè)推力器可行域Fig.4 Schematic diagram of single propeller’s thrust feasible region

綜上所述，建立了以（6）式為優(yōu)化目標(biāo)方程，u= ［ΔT1×mΔα1×m］T為優(yōu)化變量，（7）、（8）式和（9）式為約束方程的優(yōu)化模型。 但是對(duì)于半潛式平臺(tái)動(dòng)力定位系統(tǒng)而言，從圖2中可看出，同一浮箱上推力器間安裝距離較近，受推力器間尾流的相互影響，勢(shì)必會(huì)造成推力器推力損失問(wèn)題，若不加以考慮，不僅會(huì)增加運(yùn)行成本，更會(huì)導(dǎo)致定位失位問(wèn)題，本文從約束推力方向入手，最大程度地降低因槳—槳干擾造成的推力損失問(wèn)題，提高系統(tǒng)可靠性及經(jīng)濟(jì)性。

2 推力方向動(dòng)態(tài)約束域設(shè)置

2.1 推力損失模型

推力分配策略的首要任務(wù)是使各推力器協(xié)調(diào)作用以滿(mǎn)足控制力要求，整個(gè)系統(tǒng)中影響定位性能的除推力器自身動(dòng)力性能外，還有推力器之間的相互干擾造成推力損失，致使推力器實(shí)際推力達(dá)不到指令要求造成定位失位。

針對(duì)推力損失問(wèn)題，Ekstrom等[9]進(jìn)行了槳-槳干擾敞水試驗(yàn)，從試驗(yàn)結(jié)果中可看出尾流在相當(dāng)大的范圍內(nèi)會(huì)對(duì)下游螺旋槳產(chǎn)生影響。Nienhuis[10]給出了平板下槳-槳干擾試驗(yàn)結(jié)果。為了能在實(shí)際工程中對(duì)下游螺旋槳推力減額進(jìn)行定量估算，Dang等[11]分析了槳間距和推力方向?qū)ν屏p失的影響，并總結(jié)出了槳—槳干擾推力損失經(jīng)驗(yàn)公式，前后推力器同向布置，槳間距變化時(shí)推力損失經(jīng)驗(yàn)公式為

式中：t為推力減額因數(shù)，T0為敞水系柱推力，T為下游螺旋槳所產(chǎn)生的推力，x為兩螺旋槳間距離。進(jìn)一步，固定槳間距，推力器間相對(duì)推力方向變化時(shí)的推力損失經(jīng)驗(yàn)公式為

式中：Φ為兩推力器之間的相對(duì)夾角，t為x/D為定值、Φ=0時(shí)的推力減額系數(shù)。

2.2 推力可行域設(shè)置

在動(dòng)力定位系統(tǒng)中，若不對(duì)推力器推力方向加以限制，總會(huì)出現(xiàn)如圖5所示的情況，下游推力器完全處于上游推力器的尾流中，由推力損失模式經(jīng)驗(yàn)公式可知，此時(shí)下游推力器推力損失最嚴(yán)重。文獻(xiàn)[8]針對(duì)推力損失問(wèn)題總結(jié)了常用的3種處理方法：（1）接受推力損失，通過(guò)增加推力轉(zhuǎn)速來(lái)彌補(bǔ)推力損失；（2）設(shè)置推力禁區(qū)，如圖6（a）陰影部分所示，對(duì)推力方向不施加約束，當(dāng)推力方向在該禁區(qū)時(shí)設(shè)置推力為零，；（3）設(shè)置推力禁區(qū)，始終保持推力器出力方向處于禁區(qū)外。方法（1）最易實(shí)現(xiàn)，但會(huì)明顯增加燃油消耗；方法（2）中，當(dāng)1號(hào)推力器推力為零時(shí)（如圖6（a）所示），由于整體系統(tǒng)的推力突然減小從而導(dǎo)致短時(shí)失位問(wèn)題，并且若其余推力器已滿(mǎn)負(fù)荷工作，那么由于無(wú)法提供足夠的推力亦會(huì)導(dǎo)致定位失位；方法（3）中如果分配指令要求1號(hào)推力器的推力方向從a快速到達(dá)b，那么由于禁止角的存在，將使整個(gè)系統(tǒng)耗能無(wú)法在最短時(shí)間內(nèi)實(shí)現(xiàn)最優(yōu)。基于以上方法的缺點(diǎn)，本文依據(jù)推力損失模型，通過(guò)動(dòng)態(tài)設(shè)置下游推力器推力可行域來(lái)實(shí)現(xiàn)推力分配。

圖6 推力器推力可行域示意圖Fig.6 Schematic diagram of thrust feasible region

針對(duì)如圖2所示定位系統(tǒng)，本文將相臨的兩個(gè)推力器分為一組，如將1號(hào)和2號(hào)分為一組，首先，對(duì)每組中的各推力器推力方向按圖6（a）進(jìn)行分區(qū)，陰影區(qū)域?yàn)闃?槳干擾最嚴(yán)重范圍；其次，若上游推力器推力方向進(jìn)入陰影區(qū)域，啟動(dòng)推力方向約束條件以避免槳-槳干擾，下游推力器推力方向可行域根據(jù)主推力器的推力角進(jìn)行設(shè)置，下面針對(duì)第1組推力器介紹其可行域的具體設(shè)置方法。

如圖6（b）所示，當(dāng)1號(hào)推力器位于上游時(shí)，其推力、推力角可行域僅需滿(mǎn)足推力器物理性能即可，其可行域如圖4所示，圖中Δα±為轉(zhuǎn)角變化范圍，ΔT±為推力變化范圍；對(duì)于2號(hào)推力器，除需滿(mǎn)足物理性能外，還應(yīng)通過(guò)設(shè)置兩推力器之間的夾角以避免受1號(hào)推力器尾流影響而造成的推力損失，因此當(dāng)1號(hào)推力器進(jìn)入陰影區(qū)域時(shí)，將2號(hào)推力器可行域設(shè)置為圖6（b）中所示，其約束方程為

綜合（8）、（9）、（12）式及（13）式，推力器可行域歸納為

（14）式建立了第一組推力器的約束條件，其余組可依此設(shè)計(jì)。若同組推力器未發(fā)生推力干擾問(wèn)題，則（14b）式約束條件將不被考慮。（14）式根據(jù)各個(gè)推力器的實(shí)時(shí)推力方向來(lái)決定其可行域，從而實(shí)現(xiàn)了動(dòng)態(tài)約束。

3 數(shù)值模擬分析

本文以某半潛式平臺(tái)動(dòng)力定位系統(tǒng)為研究對(duì)象，該平臺(tái)總長(zhǎng)114 m，總寬90 m，總高110 m，最大甲板可變載荷9 000 t，平臺(tái)作業(yè)吃水19.0 m，作業(yè)水深為1 500 m。波浪載荷采用Pierson-Moskowitz譜，有義波高為4.0 m，跨零周期為7.9 s，入射方向由45°到60°，利用模糊PID控制算法[15]得到控制力和控制力矩時(shí)程，如圖7所示，0-1 500 s波浪力入射方向?yàn)?5°，1 500-3 000 s入射方向?yàn)?0°。選擇W-B4-70導(dǎo)管螺旋槳作為推力器，各個(gè)螺旋槳額定功率為4 600 kW，螺旋槳的主要參數(shù)采用文獻(xiàn)[16]所給值，分別為：系柱狀態(tài)下推力系數(shù)KT0=0.445，扭矩系數(shù)KQ0=0.06，直徑D=4 m，Tmax=540 kN，Tmin=0 kN，螺旋槳轉(zhuǎn)速變化率為±0.2 s-2，則推力變化率范圍為 ΔT=±4.67 kN/s，轉(zhuǎn)角變化率范圍為 Δα=±1°/s。本文在對(duì)螺旋槳進(jìn)行布置時(shí)選擇了錯(cuò)位布置，如圖2所示，錯(cuò)位布置的主要目的是為了加大螺旋槳之間的距離，降低槳-槳干擾，表1為螺旋槳布置坐標(biāo)值。角度權(quán)矩陣Ω=10×I8×8。采用Matlab提供的遺傳算法工具箱，染色體用雙精度向量編碼，種群數(shù)為20，遺傳代數(shù)100代；適應(yīng)度函數(shù)為（6）式；采用隨機(jī)均勻選擇算子；選用Scattered法交叉，交叉概率0.2；變異過(guò)程采用默認(rèn)函數(shù)。

表1 螺旋槳坐標(biāo)值Tab.1 Propeller coordinates

圖7 控制力、控制力矩時(shí)程曲線Fig.7 Time-histories of DP control force

如圖6（a）所示，1號(hào)和2號(hào)推力器中心連線與x軸的夾角δ=55.5°。根據(jù)文獻(xiàn)[13]可知，在系柱狀態(tài)下離上游推力器約3倍直徑時(shí)，尾流速度分布為沿軸線最大的單峰狀態(tài)，尾流影響范圍直徑約為2D，基于此本文取θ=36°，即相臨推力器間理想的最小出力夾角為18°，此時(shí)下游推力器推力減額因數(shù)由（11）式計(jì)算得到為 0.882 5。1號(hào)推力器的禁止角為[37.5°,73.5°]，同理，2號(hào)推力器的禁止角為[217.5°,253.5°]，其余三組推力器均依此設(shè)置。

圖8 推力方向時(shí)程曲線Fig.8 Time-histories of thrust directions

圖9 推力時(shí)程曲線Fig.9 Thrust time-histories

圖11 推力變化du時(shí)程曲線Fig.11 Time-histories of thrust changes du

圖8和圖9分別為ψ=18°時(shí)的推力角時(shí)程曲線和推力時(shí)程曲線。圖8（a）為第一組推力器的推力角時(shí)程曲線，圖中斜線所示區(qū)域?yàn)?號(hào)推力器禁止角，1號(hào)推力器進(jìn)入該區(qū)時(shí)成為主推力器，此時(shí)2號(hào)推力器在依據(jù)1號(hào)推力器出力角設(shè)置的可行域內(nèi)出力。在設(shè)置次推力器可行域時(shí)，充分考慮了推力器的物理性能，即其受轉(zhuǎn)角轉(zhuǎn)動(dòng)速率的限制，因此在短時(shí)間內(nèi)無(wú)法完全實(shí)現(xiàn)避免推力器間的推力損失問(wèn)題，從圖8（a）的700-950 s區(qū)間，以及圖8（b）和8（c）中均可看出。當(dāng)主推力器在禁止區(qū)作用時(shí)間較長(zhǎng)時(shí)，推力器的轉(zhuǎn)角得到充分轉(zhuǎn)動(dòng)，從圖8（a）的1 800 s到2 550 s間可看出兩推力器存在明顯的夾角并在約1 900 s處達(dá)到理想狀態(tài)，此時(shí)下游推力器推力損失降低到了最小。當(dāng)推力器均未處于推力禁止角時(shí)，兩者出力角基本一致。圖9（a）為第一組推力器推力時(shí)程曲線，由于2號(hào)推力器出力角受限，此時(shí)1號(hào)推力器出力較大。從圖中可看出，第三組情況與第一組基本相同。圖8（b）中，在約850 s以后推力角大于2π，主要原因是推力器出力角范圍未受限制，其出力角以2π為周期。為了方便比較，對(duì)3號(hào)推力器的推力角進(jìn)行取余，使其推力角落在[-2π,2π]內(nèi)，修正后可看出，當(dāng)該組推力器均未處于禁止角時(shí)出力方向基本一致。

為了保證動(dòng)力定位系統(tǒng)的安全性，每個(gè)推力器推力需留出20%的安全裕度，從圖9中可看出推力均未超過(guò)432 kN。圖10和圖11分別為推力角變化da和推力變化du的時(shí)程曲線，從圖中可以看出推力變化未超出[-4.67 kN/s 4.67 kN/s]范圍，轉(zhuǎn)角變化率亦滿(mǎn)足[-1°/s 1°/s]。圖10和圖11說(shuō)明該分配方法滿(mǎn)足推力器的物理性能。

4 結(jié) 論

本文在充分考慮了推力器物理性能的基礎(chǔ)上，以推力器推力變化率和推力方向變化率為控制參數(shù)，基于耗能最小建立了最優(yōu)分配模型，得到了以下結(jié)論：

（1）將遺傳算法應(yīng)用到推力分配策略中，可以處理復(fù)雜的邊界條件，避免了采用序列二次規(guī)劃法對(duì)可行域要求嚴(yán)格的問(wèn)題；

（2）充分考慮了推力損失模型，依據(jù)推力器的出力角可充分估計(jì)推力器間的相互影響，通過(guò)設(shè)置禁止角可有效地降低損失，提高推力器效率；

（3）以推力變化和推力角變化為優(yōu)化參數(shù)充分考慮了推力器的性能，避免由于推力器性能受限而無(wú)法實(shí)現(xiàn)指令要求造成的失位問(wèn)題。通過(guò)數(shù)值模擬結(jié)果可以看出該方法在推力分配問(wèn)題上是行之有效的。