基于GPU的大規模無人機編隊控制并行仿真方法*

李文光,王 強,曹 嚴

(北京理工大學宇航學院, 北京 100081)

0 引言

相比于單個無人機,多無人機作戰具有顯著的優勢,比如作戰半徑大、偵查范圍廣等,因此受到很多學者的關注,成為研究的熱點[1]。編隊控制算法是多無人機協同作戰的一項關鍵技術,基于一致性的協同編隊控制算法是眾多編隊控制算法的一種,也是一種典型的分布式算法,與傳統的集中式算法相比,具有通信與控制結構靈活、無人機規模不受限制等優點,因此在編隊控制問題中具有重要的應用[2]。

仿真為算法的可行性驗證提供了良好的手段。但是當仿真模型規模較大時,仿真耗時問題相對突出,其成為制約仿真發展的重要因素。為提高仿真效率,并行仿真憑借其高效、可重用等特點越來越多地受到青睞[3]。并行仿真的主要任務是將仿真目標分解為多個子目標并將其分布在不同處理機上同時仿真,從而提高仿真的整體效率。在大規模的無人機編隊仿真中,由于無人機數量龐大,同樣會面臨仿真耗時長的問題?；诖?文中研究了一種并行仿真方法,并充分利用GPU強大的計算能力和高效的并行性以提高仿真效率。

1 基于一致性的協同控制建模

無人機的運動學方程為：

(1)

式中：xi(t)為無人機的位置向量;vi(t)為無人機的速度向量；ui(t)為無人機的控制向量；N為無人機數量。

(2)

(3)

ui(t)=ui-form(t)+ui-vel(t)

(4)

根據個體控制律式,可得系統的控制向量u(t)：

(5)

式中：n為無人機模型維度,L為系統的Laplace矩陣,E與期望的編隊隊形相關,具體定義如下:

(6)

(7)

Li=[li1,li2,…,liN]

(8)

(9)

整理可得閉環系統的狀態方程：

(10)

2 復雜系統的并行仿真方法

并行仿真的第一步是將整個模型拆分為多個子模型以減小模型的規模。編隊控制模型是一個連續時間模型,其一般形式如式(11)所示：

(11)

而f=[f1,f2,…,fN]T為狀態轉移矩陣,x=[x1,x2,…,xN]T為狀態向量。

連續時間模型的分割實際上就是將式(11)中的N個等式方程劃分為M組,每組中的全部狀態變量及其狀態轉移方程即組成了新的子模型Si(i=1,2,…,M)。模型分割要考慮的首要問題是降低子模型間的耦合,因為耦合是造成時延誤差的關鍵因素。

構造N×N階雅可比矩陣,如式(12)：

(12)

圖1 矩陣變換最終結果示意圖

并行仿真的第二步是將已分割的子模型分布在多個仿真機上并行求解。由于子模型間存在耦合,需要以合適的通信步長完成子模型間的數據通信。通信步長的選取至關重要,是影響求解誤差的重要因素[4]。

3 GPU體系結構與CUDA

3.1 GPU及其內部體系結構

GPU(graphic processing unit)是一種擴展的計算設備,稱為協處理器。與CPU相比,GPU的優勢體現在其高效的并行性、強大的浮點計算能力等。

GPU的優勢得益于其內部體系結構,GPU有多種架構,比如NVIDIA公司的Tesla、Fermi、Kepler、Maxwell、Pascal等等,不同GPU架構的內部體系結構略有差別[5]。以Fermi為例,其內部體系結構簡化圖如圖2所示。

圖2 Fermi架構GPU內部體系結構

圖2包含了幾個GPU中非常重要的結構,分別是：

1)Core：流處理器,也稱為SP(streaming processor)。是GPU中最基本的計算單元,能單獨完成雙精度運算和32位整數運算。

2)寄存器：32位寄存器,是每個線程私有的,用于存儲局部變量。

3)共享存儲器：用于存儲共享變量。與寄存器相比,共享存儲器中的數據是一個線程塊中所有線程共享的,但訪問速度較慢。

4)線程調度器：用于線程調度。

5)SM：多核流處理器。一個GPU中包含一個或多個SM,一個SM通常包含線程調度器、共享存儲器、多個SP以及上萬個寄存器。

6)全局存儲器：在某種意義上等同于GPU的顯存,一個kernel函數中所有線程都能訪問,存儲空間大,但訪問速度最慢。

雖然每個SM中有數以萬計的SP,但并不是所有的SP都同時被調度。warp是SM中的線程調度單元,每個warp包含連續的32個線程,物理上占用32個SP用于計算。SM在任一時刻按照單指令多數據模式通過線程調度器執行一個warp中的所有線程,也就是說,同一個warp中的32個線程同時執行同一條指令,但分別處理各自的數據。

3.2 CUDA

CUDA是NVIDIA公司推出的編程模型,為用戶提供簡單的接口以實現CPU和GPU異構系統開發。

CUDA編程模型中引入主機端(CPU)和設備端(GPU)的概念,一個完整的CUDA程序包括主機端和設備端兩部分代碼。主機端代碼在CPU上執行；設備端代碼又稱為kernel函數,在GPU上執行。

啟動kernel函數時需要指定線程網格劃分,之后CUDA運行時系統生成一個兩級結構的線程網格：一個網格由多個線程塊組成,每個線程塊又由多個線程組成。所有線程執行同一個kernel函數,每個線程在運行時在物理上占用一個SP和多個寄存器(寄存器占用數量跟kernel函數中代碼相關),因此,GPU設備對kernel函數的網格劃分是有數量限定的[6]。

此外,CUDA提供了一系列函數用于操作和管理GPU設備。比如數據傳輸函數cudaMemcpy (),用于完成CPU與GPU間數據的傳輸；線程管理函數cudaDeviceSynchronize(),用于CPU和GPU間的同步等等[5]。

4 仿真設計與驗證

仿真實驗所用計算機的CPU為Inter core i7,頻率3.4 GHz,內存8 GB；GPU為Fermi體系結構,SM數量為1,Core數量為48,寄存器數量為32 768個。

4.1 模型分割

首先對模型進行分割。系統的雅克比矩陣即為式(10)中矩陣A,計算可得系統Laplace矩陣如式(13)所示：

(13)

將式(13)代入矩陣A,可以看出,每架無人機的位置向量和速度向量之間存在耦合關系；不同無人機之間一定存在且只存在與相連兩架無人機的耦合。因此,同一個無人機模型的所有狀態變量應通過行列置換作為一個整體,不同無人機模型按位置順序排列。此時,在任一位置將模型分割,分割位置相連兩無人機模型間的耦合即為分開后兩模型間的耦合,子模型間的耦合最小。考慮到負載均衡和實驗設備的限制,將整個模型拆分為5個子模型,每個子模型包含連續的200架無人機模型。

每個子模型定義如下：

(14)

4.2 CUDA程序設計與優化

求解式(14)所示的子模型,設計的CUDA程序偽代碼如圖3所示。其中,kernel函數實現的功能是迭代一個仿真步長。子模型的求解采用Euler法,因此kernel函數的具體實現是計算ξi=ξi+h(Aiξi(t)+Bi+ui)。在該kernel函數中,網格劃分為1×200,每個線程求解不同無人機的位置狀態量和速度狀態量。

圖3 子模型求解CUDA偽代碼

GPU中,當一個warp中的線程訪問全局存儲器的連續地址時,這些線程的訪問請求會被合并成對連續單元的合并訪問,以提高數據讀取速度?；诖?考慮對上述程序優化。在kernel函數中,每個線程求解時需要獲取系數矩陣A和狀態向量ξ的數據。在一個warp中,32個線程同時從全局存儲器中獲取矩陣A的同列不同行的元素,如圖4中虛線框所示,此時每個線程同時從矩陣A中讀取的元素是不連續的。

圖4 優化前訪問矩陣元素示意圖

考慮將矩陣A轉置后再復制到GPU存儲器中。這時,kernel函數在計算矩陣相乘時如圖5所示,同一warp中不同線程訪問的是連續的元素。

將設計的CUDA程序分布到不同工作站上求解。不同工作站間的通信采用集中式通信：由另外一臺計算機作為服務器,用于接受源節點的數據并發送到目標節點,并且集中控制通信步長,同步子模型的仿真時鐘。仿真實驗分為三組：第一組完全在CPU下求解；第二組采用未優化的CUDA程序求解子模型；第三組采用優化后的CUDA程序求解子模型。三組仿真耗時情況如圖6所示。

圖5 優化后訪問矩陣元素示意圖

圖6 CPU、GPU仿真耗時對比圖

圖6體現出相比CPU,GPU在求解這一類問題上更加高效,而且隨著無人機規模的增大,其加速效果更加顯著。同時也體現出經過優化后的CUDA程序,其效率有了明顯的改善。

4.3 通信步長對誤差和加速比的影響

通信步長是影響誤差與仿真耗時的關鍵因素。為探究通信步長對誤差和仿真耗時的實際影響,并為實際仿真應用中通信步長的選取提供依據,設計了本節仿真實驗。實驗結果如圖7和圖8所示,縱坐標分別為相對誤差和加速比,橫坐標為通信步長與仿真步長的比值,即X=H/h；相對誤差是指GPU下并行仿真結果相對CPU下串行仿真結果的誤差,取所有狀態量相對誤差的均值作為觀測值；加速比是指GPU下并行仿真耗時與CPU下串行仿真耗時的比值。

圖7 X與相對誤差關系圖

由圖7可以看出,X和誤差成正比；由圖8可以看出,X對加速比的影響在X很小時相對較大,隨著X的增大,其對加速比的影響逐漸減小。結合圖7、圖8可知,當X取仿真步長的5～15倍時,具有較好的效果。

圖8 X與加速比關系圖

5 結論

從仿真結果可知,通信步長對求解誤差以及求解效率有顯著的影響,當通信步長取仿真步長的5～15倍時,相對誤差僅為0.01～0.02,驗證了該并行仿真方法的正確性；加速比達到15～20,驗證了方法的高效性。圖6體現出了GPU在大規模數值仿真中的優勢,并且實驗所用GPU性能較差,計算能力(compute capab-ility)僅為2.1,最新的GPU計算能力能夠達到7.1,因此GPU在大規模仿真中有著廣闊的應用前景。