基于MATLAB的自主水下航行器經(jīng)典控制算法仿真分析

劉華棟，康寶臣，馮麗娜

（1.海軍裝備部，山西 太原 030027；2.山西汾西重工有限責(zé)任公司，山西 太原 030027）

0 引言

在海洋資源開(kāi)發(fā)和海軍軍事任務(wù)雙重需求的牽引下，水下無(wú)人航行器以及自主水下航行器成為當(dāng)前的熱門(mén)研究課題。自主水下航行器就是離開(kāi)發(fā)射平臺(tái)后仍具備自行航行能力至預(yù)定水域偵查、探測(cè)和打擊為一體的航行器。為進(jìn)一步提高布設(shè)的隱蔽性，研制自主水下航行器就顯得特別重要。自主水下航行器能在敵方防潛設(shè)施外啟動(dòng)航行，利用其自身的動(dòng)力裝置航行到預(yù)定海域?qū)持匾劭诨蛩蜻M(jìn)行偵查、探測(cè)和打擊。

自主水下航行器無(wú)論是在軍用市場(chǎng)還是民用市場(chǎng)，全球各國(guó)都在加緊研制與開(kāi)發(fā)采用最新科技成果的水下無(wú)人航行器。未來(lái)的水下無(wú)人航行器結(jié)構(gòu)更加簡(jiǎn)單，功能多樣，自適應(yīng)能力強(qiáng)，可游刃有余地執(zhí)行復(fù)雜的任務(wù)[1]。而各種控制算法，如神經(jīng)網(wǎng)絡(luò)控制、模糊控制、滑?？刂坪湍：?PID控制等智能控制算法均被應(yīng)用到了UUV的運(yùn)動(dòng)控制中[2-7]。文獻(xiàn)[2]采用了經(jīng)典的 PID控制算法，文獻(xiàn)[3]采用了以太網(wǎng)的分布式控制系統(tǒng)設(shè)計(jì)，文獻(xiàn)[4]采用了滑模理論控制算法，文獻(xiàn)[5]則將經(jīng)典 PID算法、自抗擾技術(shù)、模糊策略和滑?？刂扑惴ǚ謩e進(jìn)行了分析與設(shè)計(jì)，文獻(xiàn)[6]采用了神經(jīng)網(wǎng)絡(luò)控制算法，文獻(xiàn)[7]采用了魯棒滑?？刂扑惴?。

進(jìn)入21世紀(jì)，仿真技術(shù)得以發(fā)展的主要原因，是它所帶來(lái)的巨大社會(huì)經(jīng)濟(jì)效益。仿真主要應(yīng)用于航空、航天、航海、電力、化工以及其他工業(yè)過(guò)程控制等工程技術(shù)領(lǐng)域。在航空工業(yè)方面，采用仿真技術(shù)使大型客機(jī)的設(shè)計(jì)和研制周期縮短20%。利用計(jì)算機(jī)仿真，不僅節(jié)省大量燃料和經(jīng)費(fèi)，而且不受氣象條件和場(chǎng)地的限制。在航天工業(yè)方面，采用仿真實(shí)驗(yàn)代替試航試驗(yàn)可使實(shí)彈試驗(yàn)的次數(shù)減少80%。在電力工業(yè)方面采用仿真系統(tǒng)對(duì)核電站進(jìn)行調(diào)試、維護(hù)和排除故障，1年即可收回建造仿真系統(tǒng)的成本?，F(xiàn)代仿真技術(shù)不僅應(yīng)用于傳統(tǒng)的工程領(lǐng)域，而且日益廣泛地應(yīng)用于社會(huì)、經(jīng)濟(jì)、生物等領(lǐng)域，如交通控制、城市規(guī)劃、資源利用、環(huán)境污染防治、生產(chǎn)管理、市場(chǎng)預(yù)測(cè)、世界經(jīng)濟(jì)的分析和預(yù)測(cè)、人口控制等。因此，利用仿真技術(shù)來(lái)研究這些系統(tǒng)就具有更為重要的意義。

1 總體參數(shù)與流體動(dòng)力參數(shù)

自主水下航行器航行控制數(shù)學(xué)模型屬于灰箱系統(tǒng)，其運(yùn)動(dòng)規(guī)律遵循牛頓經(jīng)典力學(xué)定律，質(zhì)量守恒定律和能量守恒定律。數(shù)學(xué)模型中的流體動(dòng)力參數(shù)可通過(guò)水洞試驗(yàn)、風(fēng)動(dòng)試驗(yàn)或自主水下航行器流體動(dòng)力參數(shù)辨識(shí)得到，某航行器所需要的總體參數(shù)明細(xì)如表1所示。

表1 總體參數(shù)Table 1 General parameters

表1中：流體的密度ρ=999.866 kg/m3；航行器的體積V=1.1 m3；航行器的質(zhì)量m=1 300 kg；航行器的長(zhǎng)度L=7 m；航行器的最大橫截面積S=0.224 m2；航行器質(zhì)心距浮心的距離xc=-25 mm；航行器質(zhì)心下移量yc=-20 mm；航行器質(zhì)心側(cè)移量zc=0；航行器繞航行器坐標(biāo)系x軸轉(zhuǎn)動(dòng)的轉(zhuǎn)動(dòng)慣量Jxx=602kg·m；航行器繞航行器坐標(biāo)系軸y轉(zhuǎn)動(dòng)的轉(zhuǎn)動(dòng)慣量Jyy=6 1232kg·m；航行器繞航行器坐標(biāo)系軸z轉(zhuǎn)動(dòng)的轉(zhuǎn)動(dòng)慣量Jzz=6 083 kg·m2。

流體動(dòng)力參數(shù)明細(xì)如表2所示。

表2 流體動(dòng)力參數(shù)Table 2 Hydrodynamic parameters

2 縱向彈道數(shù)值模擬

2.1 縱向運(yùn)動(dòng)存在的條件

航行器在空間的運(yùn)動(dòng)可分解為縱向運(yùn)動(dòng)和側(cè)向運(yùn)動(dòng)，這是一種簡(jiǎn)化的研究方法。在滿(mǎn)足一定條件下，這種簡(jiǎn)化研究方法，在工程精度范圍內(nèi)是可以允許的。航行器的縱向運(yùn)動(dòng)為一平面運(yùn)動(dòng)，其運(yùn)動(dòng)平面與航行器的縱對(duì)稱(chēng)面重合，且為空間垂直平面（地面坐標(biāo)系平面）。航行器縱向運(yùn)動(dòng)單獨(dú)存在的條件是：

1）航行器對(duì)其縱對(duì)稱(chēng)面是對(duì)稱(chēng)的，在縱向運(yùn)動(dòng)中，不產(chǎn)生側(cè)向力、偏航力矩和橫滾力矩；

2）航行器控制系統(tǒng)是理想的，使航行器縱對(duì)稱(chēng)面保持在垂直面內(nèi)；

3）航行器航行的彈道為垂直面內(nèi)的平面彈道，如爬潛和尋深彈道，直航和垂直面導(dǎo)引彈道。

2.2 原始縱向運(yùn)動(dòng)方程組-數(shù)學(xué)模型

縱向運(yùn)動(dòng)是指航行器在鉛垂平面（地面坐標(biāo)系的ox0y0平面）內(nèi)的運(yùn)動(dòng)，并且航行器坐標(biāo)系的oxy平面與運(yùn)動(dòng)平面重合。縱向運(yùn)動(dòng)由航行器重心在oxy平面內(nèi)的運(yùn)動(dòng)和繞oz軸的轉(zhuǎn)動(dòng)運(yùn)動(dòng)組成。由于地面坐標(biāo)系的ox0軸方向可以任意選取，如果將ox0軸也選在運(yùn)動(dòng)平面內(nèi)，則地面系的ox0y0平面，航行器系的oxy平面，速度系的ox1y1平面，半速度系的ox1y1平面都是同一平面。這時(shí)有9個(gè)運(yùn)動(dòng)參數(shù)，即β、υz、ωx、ωy、φ、ψ、z0、Ψ、Φc均為0，余下的非零運(yùn)動(dòng)參數(shù)為α、θ、Θ、ωz、υx、υy、υ、x0、y0，這9個(gè)運(yùn)動(dòng)參數(shù)稱(chēng)為縱向運(yùn)動(dòng)參數(shù)，描述縱向運(yùn)動(dòng)參數(shù)變化規(guī)律的方程稱(chēng)為縱向運(yùn)動(dòng)方程。

假定航行器具有oxy及oxz共2個(gè)對(duì)稱(chēng)面，于是根據(jù)航行器空間運(yùn)動(dòng)方程組可以得到如下方程[8]：

由式（1）-（9）組成的方程組就是航行器原始的縱向運(yùn)動(dòng)方程組。如果橫舵角δe已知，則由9個(gè)方程可以求解9個(gè)運(yùn)動(dòng)參數(shù)，從而確定縱向運(yùn)動(dòng)和彈道。在固定舵角下的運(yùn)動(dòng)和彈道稱(chēng)為非操縱運(yùn)動(dòng)和彈道。航行器制導(dǎo)系統(tǒng)通過(guò)傳感器測(cè)量航行器運(yùn)動(dòng)參數(shù)，通過(guò)自導(dǎo)系統(tǒng)測(cè)量目標(biāo)相對(duì)運(yùn)動(dòng)參數(shù)，根據(jù)這些參數(shù)，制導(dǎo)系統(tǒng)給出操舵指令，形成閉環(huán)系統(tǒng)。制導(dǎo)系統(tǒng)的工作原理可以寫(xiě)成制導(dǎo)方程（控制方程和導(dǎo)引方程），在縱向運(yùn)動(dòng)中需要補(bǔ)充制導(dǎo)方程，與縱向運(yùn)動(dòng)方程聯(lián)立求解縱向運(yùn)動(dòng)參數(shù)和橫舵角δe，確定航行器運(yùn)動(dòng)及其彈道。

2.3 縱向運(yùn)動(dòng)數(shù)值模擬軟件程序編制

縱向運(yùn)動(dòng)程序算法采用Rung_Kutta（四階）方法，由MATLAB實(shí)現(xiàn)。

ODE解算器時(shí)，需要將微分方程組轉(zhuǎn)化成MATLAB可接受的標(biāo)準(zhǔn)格式，如下：

ODEs由多個(gè)高階或隱式微分方程組成，將其變換成一階顯式常微分方程組。本文的航行器縱向運(yùn)動(dòng)方程組是由9個(gè)隱式微分方程構(gòu)成的ODEs，可轉(zhuǎn)換成一階顯示微分方程，然后求解出右函數(shù)。

四階龍格-庫(kù)塔法（The Fourth-Order Runge－KuttaMethod）形式如下[9]：

航行器縱向控制方程采用如下形式[10]：

3 仿真結(jié)果

本文共進(jìn)行了4種工況下的縱向運(yùn)動(dòng)仿真，包括穩(wěn)態(tài)航速8 m/s定深5 m時(shí)的縱向運(yùn)動(dòng)仿真、穩(wěn)態(tài)航速8 m/s定深50 m時(shí)的縱向運(yùn)動(dòng)仿真、穩(wěn)態(tài)航速8 m/s定深5 m無(wú)自由角時(shí)的縱向運(yùn)動(dòng)仿真和穩(wěn)態(tài)航速8 m/s定深50 m無(wú)自由角時(shí)的縱向運(yùn)動(dòng)仿真；由于篇幅限制，只列出穩(wěn)態(tài)航速8 m/s定深50 m時(shí)的縱向運(yùn)動(dòng)仿真結(jié)果如下。

圖1 深度曲線(xiàn)Fig.1 Depth curve

圖2 彈道傾角曲線(xiàn)Fig.2 Trajectory inclination curve

圖3 攻角曲線(xiàn)Fig.3 Angle of attack curve

圖4 彈道曲線(xiàn)Fig.4 Trajectory curve

圖5 俯仰角曲線(xiàn)Fig.5 Pitch angle curve

圖6 速度曲線(xiàn)Fig.6 Velocity curve

圖7 水平舵角曲線(xiàn)Fig.7 Horizontal rudder angle curve

圖8 俯仰角速度曲線(xiàn)Fig.8 Pitch angular velocity curve

通過(guò)多次仿真計(jì)算和試參數(shù)，最終確定了Ky、Kθ和δfree，并對(duì)有無(wú)自由角進(jìn)行了對(duì)比，得出了自由角的取值合適時(shí)，縱向運(yùn)動(dòng)無(wú)超調(diào)，但穩(wěn)定時(shí)間稍長(zhǎng)，自由角偏離理想值時(shí)，縱向控制系統(tǒng)會(huì)產(chǎn)生一定的振蕩和超調(diào)。

4 結(jié)束語(yǔ)

以上4種算例的仿真結(jié)果均給出了航行器在縱平面內(nèi)的深度曲線(xiàn)、彈道曲線(xiàn)、攻角曲線(xiàn)、速度曲線(xiàn)、彈道傾角曲線(xiàn)、水平舵角曲線(xiàn)、俯仰角曲線(xiàn)和俯仰角速度曲線(xiàn)。從仿真結(jié)果可以看出：由于有重力、浮力和其它因素的影響，自主水下航行器的水平舵角不為0，而是有一個(gè)很小的度數(shù)，從而產(chǎn)生升力與負(fù)浮力等作用力相平衡。此外，自由角的取值合適時(shí)，縱向運(yùn)動(dòng)無(wú)超調(diào)，但穩(wěn)定時(shí)間稍長(zhǎng)，自由角偏離理想值時(shí)，縱向控制系統(tǒng)會(huì)產(chǎn)生一定的振蕩和超調(diào)。調(diào)整好自由角的大小和傳遞系數(shù)就會(huì)很好地抹除此控制算法下系統(tǒng)產(chǎn)生的振蕩和系統(tǒng)調(diào)整所耗費(fèi)的時(shí)間。此仿真算例得出了合理的結(jié)果，已經(jīng)驗(yàn)證和支持了某自主水下航行器的湖上試驗(yàn)，為某自主水下航行器項(xiàng)目實(shí)航演示樣機(jī)的研制提供理論支持，具有一定的工程參考價(jià)值。