基于改進粒子群優化算法的船舶縱向運動參數辨識

戴運桃 ， 趙希人 ， 劉利強

（哈爾濱工程大學 a.理學院；b.自動化學院，哈爾濱 150001）

1 引 言

水動力參數是船舶操縱設計的重要參數，但是有些參數難以用實驗方法獲得，而理論計算值和實驗值都與實航實驗有較大區別。因此，水動力參數的獲得一般都是應用系統辨識技術，從船舶的運動狀態觀測數據中辨識出來，從而直接建立船舶的水動力參數和運動狀態之間的數學模型。水動力參數辨識方法一般采用經典的辨識方法如極大似然辨識方法[1]和預報誤差法[2]等。但由于縱向運動狀態對阻尼力和力矩等參數不敏感，靈敏度非常低，而常規辨識算法對靈敏度較低的參數是非常難以求解到真值的。

粒子群優化算法（Particle Swarm Optimization，PSO）[4]是一種新的基于群體智能的優化算法，它仿效生物學中進化和遺傳的過程，同時考察多個候選解，淘汰劣質解，鼓勵發展優質解，逐步提高解群體的質量，從而逼近所研究問題的最優解，為解決優化問題提供了新思路[10]。本文針對常規辨識方法在水動力參數辨識中的局限性，探索了模擬生物智能的粒子群優化算法在船舶水動力參數辨識中的應用。

2 問題描述及數學模型

2.1 狀態方程

根據船舶水動力理論[8]，船舶縱向運動的簡化方程如公式（1）所示。

式中：z表示船體的升沉，θ表示船體的縱搖，I5表示縱搖慣性矩，Δ表示排水量，FR表示水平舵升力，XR為舵升力中心至船體重心的縱向距離，Fw3表示海浪升沉干擾力，Mw5表示海浪縱搖干擾力矩；a33表示附加質量；a35，a53表示質量矩（kg·m）；a55表示轉動慣量（kg·m2）；b33，b35表示阻尼力系數；c33，c35表示復原力系數；b53，b55表示阻尼力矩系數；c53，c55表示復原力矩系數。 a33，a35，a53，a55，b33，b35，b53，b55，c33，c35，c53，c55統稱為水動力參數。

令 x1=z，x3=θ，設狀態變量為

則得到系統狀態方程如下：

其中，A=E-1A*，B=E-1B*，C=E-1C*，

para=[a33，a35，a53，a55，b33，b35，b53，b55，c33，c35，c53，c55]為待辨識參數，W表示海浪擾動力和力矩。

2.2 觀測方程

選擇狀態x1、x3為測量狀態，則

式中：

Y為觀測向量，V為二維測量噪聲，通常情況下可以認為是白噪聲，它的方差陣按一級精度的傳感器可以取為Qvv=diag[ 20.3e-4 2.25e-6]。

2.3 離散方程

在進行計算機仿真前，首先需要對狀態方程和觀測方程進行離散化。設Ts為采樣周期，在連續的控制對象前有一個零階保持器，即u（t）=u（k ），kTs＜t＜（k+ 1） Ts。 忽略量化效應，并將連續的控制對象和它前面的零階保持器一起離散化，從而使整個系統變為純粹的離散系統。則離散化后狀態方程為：

其中，Φ=eATs，

3 粒子群優化算法

粒子群優化（Particle Swarm Optimization，PSO）算法是由Kennedy和Eberhart借鑒鳥類尋找食物的自然現象提出的一類基于種群的隨機全局優化技術[4-5]。在鳥群的飛行中，每只鳥在初始狀態下處于隨機位置，且朝各個方向隨機飛行，但隨著時間的推移，這些初始處于隨機狀態的鳥通過相互學習、相互跟蹤，自組織地聚集成一個個小的群落，并以相同的速度朝著相同的方向飛行，最終整個群體聚集到同一位置—食物源。

PSO算法在可行解空間中初始化產生一組粒子，每個粒子在搜索空間中以一定的速度飛行，并根據對個體和集體的飛行經驗的綜合分析來動態調整個體自身的飛行速度，同時向著自己以前經歷過的最好位置和其它微粒曾經經歷過的最好位置飛行，依靠這種個體間的信息交換來達到整個群體的共同演化。實現復雜空間全局最優解的搜索。與其他進化算法相比，PSO算法的主要特點是：概念簡單，實現容易，需要調整的參數相對較少，可以用于解決復雜非線性優化問題。下面以一個最小化問題的求解為例，給出PSO算法的流程。

步驟1：在整個搜索空間中隨機地初始化粒子群中各粒子的位置xi和速度vi，并計算每個粒子在當前位置處的適應度值fitnessi=fitness（ xi）。相應的初始化粒子自身的最優解（個體極值pbesti）和整個粒子群的當前最優解（全局極值 gbest），即 pbesti=fitnessi，gbest=min（fitness1,…， fitnessn），gbest越小越好。

步驟2：在每次迭代過程中，粒子通過跟蹤其自身的個體極值和全局極值，按公式（6）來進行自身速度和位置的更新。

式中：vi表示粒子的速度向量；xi代表粒子的位置向量；c1、c2為常數，稱為學習因子；r1、r2是[0，]1上均勻分布的隨機數；w為慣性權重。

步驟3：更新個體極值pbesti和全局極值gbest。

步驟4：如果已經滿足停止準則，如gbest達到某個閾值或者已經達到最大迭代次數，算法中止計算；否則跳轉到“步驟2”。

4 改進粒子群優化算法

4.1 算法基本思想

許多優化算法在尋優的過程中都會出現早熟問題，粒子群算法也不例外。在粒子群算法中，各粒子是根據自身的最優位置與全體粒子的最優位置來調整搜索方向的，在算法運行的初始階段，收斂速度比較快，但運行一段時間后，速度開始減慢甚至停滯，當所有粒子的速度幾乎為0時，粒子群喪失了進一步進化的能力，可以認為算法執行已經收斂。但在許多情況下，算法并沒有收斂到全局極值，甚至連局部極值也未達到。這是因為發生該現象時粒子群高度聚集，嚴重缺乏多樣性，粒子群會長時間或永遠跳不出聚集點[9]。

針對這一問題，本文設計了一種改進粒子群優化算法，在算法運行的每一迭代步中計算全局最優粒子的適應值globalBestValue（即全局最小適應值），若在設定的迭代步changeNum內，全局最小適應值globalBestValue一直都不變化，則代表各粒子已經開始高度聚集，算法在當前的位置和速度下已無法找到更優解。此時，利用高斯變異算子變異各粒子的位置，使粒子從新的位置開始繼續搜索。隨著粒子變異次數的增加，算法不斷縮小高斯變異算子的變異范圍，從而使粒子在逐漸縮小的搜索空間進行精細尋優，直到搜索到全局最優解。

4.2 高斯變異算子設計

選擇高斯密度函數作為各粒子的位置密度分布模型。如（7）式所示。

其中，xmini，d為適應值最小的粒子的位置，σd表示本次變異中的第d個參數的方差。

算法對粒子的位置密度分布函數進行采樣，得到各粒子的位置。采樣過程中，使用函數p（ xi）所對應的隨機數生成器進行，該生成器如公式（8）所示。

其中，randl為均勻分布在[0，1]之間的一個隨機變量。

隨著算法變異次數的增加，利用公式（9）不斷縮小粒子位置的變異范圍。

其中，ρ（ 0＜ρ＜ 1）為收縮因子，k為變異次數。 在算法開始運行階段，方差 σd較大，這樣可以使粒子在較大的范圍內變異，從而保持粒子的多樣性；隨著粒子位置不斷變異，方差σd越來越小，算法開始在更小范圍內進行精細搜索。

5 船舶縱向運動參數辨識算法設計

本文設計的參數辨識算法的思想是將粒子的位置作為待辨識的水動力參數。首先，用給定的水動力參數求得一組理論觀測值；然后，在算法運行過程中，利用各粒子的位置（辨識得到的水動力參數）求得實際觀測值，并使用適應值評價函數對本組水動力參數作出評價。這一評價即為粒子當前位置的適應度值。各粒子使用個體最優適應度值和全局最優適應度值進行速度更新和位置更新，算法迭代運行，從而引導各粒子向著適應度值高的方向前進。

適應值評價函數是評估算法辨識得到的水動力參數好壞的重要依據。如公式（10）所示，在算法中我們使用理論觀測值和實際觀測值之間對應數據的誤差平方和的均值作為適應值評價函數。函數值越小，代表該組辨識得到的水動力參數越好。

其中，Yreali是根據實船試驗得到的水動力參數計算所得的理論觀測值，Yi是根據辨識得到的水動力參數計算得到的實際觀測值，dataNum是觀測次數。算法的流程圖及運算步驟如下：

步驟1：初始化，設定學習因子c1、c2和慣性權重w，在允許范圍內，對各粒子的初始位置和速度進行隨機選取，利用給定的水動力參數和船舶縱向運動方程解算理論觀測值；

步驟2：使用每個粒子對應的位置信息（水動力參數）和船舶縱向運動方程解算實際觀測值，并使用適應值評價函數得到各粒子當前位置的適應度值；

步驟3：對各粒子比較其當前位置適應度值和個體極值pbest，如果更好，則用當前位置適應度值來更新pbest；

步驟4：用每個粒子的個體極值pbest與全局極值gbest比較，若更好，則更新gbest；

步驟5：判斷是否滿足終止條件，若滿足，則輸出全局極值gbest和其所對應的位置；否則，轉到“步驟 6”。

步驟6：判斷全局最優適應值是否有連續changeNum次未改變，如果是，則按照公式（8）重新生成粒子位置，否則，轉入“步驟2”。

6 仿真分析

圖1 參數辨識算法流程圖Fig.1 Flowchart of parameter identification algorithm

針對上述設計的參數辨識算法，本文使用Visual C++和MATLAB作為工具進行了計算機仿真。由文獻[3]可知，c33,c35,c53,c55為常數，因此，在辨識過程中，將c33,c35,c53,c55當作已知參數，然后對a33,a35,a53,a55，b33,b35,b53,b55進行辨識。在仿真實驗中，二維測量噪聲V取均值為零的高斯白噪聲，輸入舵角是幅值±10 為偽隨機數，理論觀測值使用文獻[8]中海情3 級航速 18 節，航向分別為 0°，45°，90°，135°，180°的五組水動力參數計算得到。

本次實驗采用基本的PSO算法與本文改進的PSO算法進行比較，兩算法基本參數設置相同，各參數如下：粒子數60，c1=c2=1.496 2，慣性權重隨迭代次數減少，即從初始的0.9線性遞減到0.4，各參數的搜索范圍為±200%，算法迭代次數為1 000次。在改進算法特有的參數中，收縮因子ρ=0.9，changeNum=20。 針對海情 3 級航速 18 節，航向分別為 0°，45°，90°，135°，180°的五組水動力參數的辨識問題，算法收斂的標準是所有的水動力參數的辨識誤差均收斂到±10%以內。

兩算法對于每個航向各運行100次，算法收斂速度和計算成功率如下表1所示。表1中各參數解釋如下：meanN表示所有達到收斂要求的迭代次數的均值；minN表示達到收斂要求的最小迭代次數；而maxN則是達到收斂要求的最大迭代次數；succeed是100次搜索中達到收斂要求的次數的比例。

表1 PSO與改進PSO算法在水動力參數辨識問題的對比Tab.1 Comparison between PSO and improved PSO in hydrodynamic parameter identification

從表1中我們可以看到，改進的PSO算法在求解船舶縱向運動水動力參數辨識問題時穩定性遠高于基本的PSO算法，而且平均迭代次數也小于PSO算法的迭代次數。

表2給出了辨識得到的船舶縱向運動水動力參數仿真結果。其中，水動力參數的理論值是基于切片理論通過船模在水池做大量實驗并經測試計算得到，而辨識值是通過改進的PSO算法辨識得到，N表示收斂到次數是±10%以內所需的迭代次數，相對誤差的計算是：相對誤差=（理論值-辨識值）/理論值×100%。

表2 船舶縱向運動參數辨識結果Tab.2 Parameter identification result of ship vertical motion

圖2給出了3級海情下，航速18Kn，航向45°時算法運行過程中各水動力參數的變化曲線。

圖2中橫坐標表示迭代次數，縱坐標為各水動力參數在不同迭代次數時的辨識值。直線表示各水動力參數的理論值，曲線為各參數在不同迭代次數時的值。從圖中可以看到，粒子出現震蕩收斂現象，在最初始迭代時震蕩較大，而后緩慢收斂。這是因為在初始階段，變異的粒子的位置是在較大方差下給出的，保證了粒子的多樣性，增加了跳出局部最優的概率。在隨后的過程中方差逐漸減小，粒子位置緩慢趨于理論值。

圖3和圖4分別為升沉曲線圖和縱搖曲線圖。橫坐標為觀測的次數，縱坐標分別為升沉值（單位為米）和縱搖值（單位為度）。“*”為理論上的觀測值，曲線表示通過辨識的參數得到的模型輸出值。從圖上可以看出，理論觀測值與模型輸出值基本吻合，說明辨識結果能較準確地描述船舶的動力學特性。

圖2 海情3級，航速18Kn，航向45°的船舶縱向運動參數收斂曲線圖Fig.2 Parameter identification of the ship vertical motion for 3 sea condition,18Kn of ship speed,45°of course

圖3 升沉的辨識曲線Fig.3 Identification of heave

圖4 縱搖的辨識曲線Fig.4 Identification of pitch

7 結 論

船舶水動力參數辨識是獲取船舶水動力參數的一種重要手段。本文使用改進的粒子群優化算法設計了一種船舶縱向運動水動力參數辨識算法，介紹了參數辨識算法的思想，給出了算法流程和適應值函數，并進行了計算機仿真。仿真結果表明，改進的PSO算法在求解船舶縱向運動水動力參數辨識問題時穩定性好，而且辨識精度高，能夠滿足實際需要。

[1]丁文鏡,羅仁凡等.受控航行體水動參數的極大似然辨識[J].清華大學學報（自然科學版）,1992,2(32):89-95.

[2]丁文鏡,羅仁凡等.辨識航行器水動力參數的預報誤差法[J].艦船科學技術,1994(5):22-27.

[3]Dong R R,Joseph Katz.On the structure of bow waves on a ship model[J].Journal of Fluid Mechanics,1997,346:77-115.

[4]Kennedy J,Eberhart R.Particle Swarm Optimization[C]//IEEE Int’l Confon Neural Networks.Perth,Australia,1995:1942-1948.

[5]Eberhart R,Kennedy J.A new optimizer using particle swarm theory[C]//In:Proc of the Sixth International Symposium on Micro Machine and Human Science.Nagoya,Japan,1995:39-43.

[6]Shi Y,Erberhart R C.A modified particle swarm optimizer[J].In:Proc of the IEEE CEC,1998:69-73.

[7]Liu Bo,Wang Ling,et al.Improved particle swarm optimization combined with chaos[J].In:Chaos,Solitons and Fractals,2005,25:1261-1271.

[8]陳虹麗等.船舶縱向運動水動力參數模型的建立及分析[J].儀器儀表學報,2006,27(6):1120-1121.

[9]李 莉,李洪奇.基于混合粒子群算法的高維復雜函數求解[J].計算機應用,2007,27(7):1754-1756.

[10]陳緯琪,顏 開等.水下航行體水動力參數智能辨識方法研究[J].船舶力學,2007,11(1):40-46.