基于改進粒子群優(yōu)化算法的船舶縱向運動參數(shù)辨識

戴運桃 ， 趙希人 ， 劉利強

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

1 引 言

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

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

2 問題描述及數(shù)學(xué)模型

2.1 狀態(tài)方程

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

式中：z表示船體的升沉，θ表示船體的縱搖，I5表示縱搖慣性矩，Δ表示排水量，F(xiàn)R表示水平舵升力，XR為舵升力中心至船體重心的縱向距離，F(xiàn)w3表示海浪升沉干擾力，Mw5表示海浪縱搖干擾力矩；a33表示附加質(zhì)量；a35，a53表示質(zhì)量矩（kg·m）；a55表示轉(zhuǎn)動慣量（kg·m2）；b33，b35表示阻尼力系數(shù)；c33，c35表示復(fù)原力系數(shù)；b53，b55表示阻尼力矩系數(shù)；c53，c55表示復(fù)原力矩系數(shù)。 a33，a35，a53，a55，b33，b35，b53，b55，c33，c35，c53，c55統(tǒng)稱為水動力參數(shù)。

令 x1=z，x3=θ，設(shè)狀態(tài)變量為

則得到系統(tǒng)狀態(tài)方程如下：

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

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

2.2 觀測方程

選擇狀態(tài)x1、x3為測量狀態(tài)，則

式中：

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

2.3 離散方程

在進行計算機仿真前，首先需要對狀態(tài)方程和觀測方程進行離散化。設(shè)Ts為采樣周期，在連續(xù)的控制對象前有一個零階保持器，即u（t）=u（k ），kTs＜t＜（k+ 1） Ts。 忽略量化效應(yīng)，并將連續(xù)的控制對象和它前面的零階保持器一起離散化，從而使整個系統(tǒng)變?yōu)榧兇獾碾x散系統(tǒng)。則離散化后狀態(tài)方程為：

其中，Φ=eATs，

3 粒子群優(yōu)化算法

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

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

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

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

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

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

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

4 改進粒子群優(yōu)化算法

4.1 算法基本思想

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

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

4.2 高斯變異算子設(shè)計

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

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

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

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

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

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

5 船舶縱向運動參數(shù)辨識算法設(shè)計

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

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

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

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

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

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

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

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

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

6 仿真分析

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

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

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

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

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

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

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

表2 船舶縱向運動參數(shù)辨識結(jié)果Tab.2 Parameter identification result of ship vertical motion

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

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

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

圖2 海情3級，航速18Kn，航向45°的船舶縱向運動參數(shù)收斂曲線圖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 結(jié) 論

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

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

[2]丁文鏡,羅仁凡等.辨識航行器水動力參數(shù)的預(yù)報誤差法[J].艦船科學(xué)技術(shù),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]陳虹麗等.船舶縱向運動水動力參數(shù)模型的建立及分析[J].儀器儀表學(xué)報,2006,27(6):1120-1121.

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

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