應(yīng)用圖形處理器實現(xiàn)無功優(yōu)化并行計算

黃玉龍 劉明波

（華南理工大學(xué)電力學(xué)院 廣東省綠色能源技術(shù)重點實驗室 廣州 510640）

1 引言

非線性原對偶內(nèi)點法具有良好的收斂性和魯棒性[1-2]，為大型電力系統(tǒng)無功優(yōu)化的快速求解提供了可能。然而，隨著現(xiàn)代電網(wǎng)規(guī)模的迅速擴大，傳統(tǒng)在CPU上實現(xiàn)的無功優(yōu)化還無法滿足在線計算和實時控制的需求。為此，文獻[3]提出一種粗粒度的多機并行計算方法求解動態(tài)無功優(yōu)化問題，4機并行計算加速比最高只可以達到1.81。文獻[4-5]建立多分區(qū)無功優(yōu)化并行計算模型加速計算，118節(jié)點系統(tǒng)分3區(qū)并行計算加速比只有1.89[4]，而且存在1%左右的誤差。

近年來，圖形處理器（Graphic Processing Unit，GPU）不斷演變?yōu)楦叨炔⑿小⒍嗑€程多核處理單元，具有很高的處理速度和存儲帶寬。GPU越來越適用于圖形處理外的通用并行計算[6-8]，如代數(shù)運算[9]、最優(yōu)化計算[10]、電力系統(tǒng)暫態(tài)穩(wěn)定仿真[11]、電網(wǎng)可視化[12]和醫(yī)療成像[13]等。文獻[9]基于GPU求解線性方程組，獲得一定的加速效果，但還只能實現(xiàn)方程維數(shù)低于4096的單精度浮點運算。

為了進一步提高大電網(wǎng)無功優(yōu)化計算的速度，本文利用NVIDIA 公司GeForce GTX 285 GPU和臺式計算（3.0GHz Intel雙核CPU，2.0GB內(nèi)存）構(gòu)建并行計算環(huán)境，GPU和CPU之間接口采用PCIe 2.0x16總線，帶寬為8GB/s。通過兩者的協(xié)同配合來并行實現(xiàn)內(nèi)嵌離散懲罰的非線性原對偶內(nèi)點算法。將計算密集部分—線性修正方程的求解在圖形處理器上實現(xiàn)，其余部分在CPU上執(zhí)行。這樣可以充分利用GPU的高速浮點運算和高帶寬，極大地提高了整體計算速度，可以適應(yīng)大電網(wǎng)的計算需要。

2 內(nèi)嵌離散懲罰的非線性原對偶內(nèi)點法

無功優(yōu)化在數(shù)學(xué)上表現(xiàn)為典型的連續(xù)和離散變量混合的非線性優(yōu)化問題[1-2]，其模型可描述為

引入松弛變量，將不等式約束式（3）和式（4）轉(zhuǎn)變成等式約束：

式中，su，sl，sh，sw為由松弛變量組成的列向量，

引入對數(shù)壁壘函數(shù)消去松弛變量的非負(fù)性約束，對等式約束引入拉格朗日乘子，并引入正曲率二次罰函數(shù)來處理QC和ΤB這類離散變量[1-2]，從而形成如下的增廣拉格朗日函數(shù)：

式中，y，yu，yl，yh，yw為拉格朗日乘子向量，，且yu，yh≤0，yl，yw≥0；μ為壁壘參數(shù)，且μ≥0；vj為離散變量x1j的罰因子，x1jb為其鄰域中心。

根據(jù)Karush-Kuhn-Tucker最優(yōu)性條件可得滿足最優(yōu)性必要條件的非線性方程組，然后線性化得到如下的簡化修正方程[1-2]：

式中，系數(shù)矩陣中各元素表達式見式（12）～式（15）；Yu，Yl，Yh，Yw，Su，Sl，Sh，Sw分別為以yu，yl，yh，yw，su，sl，sh，sw的分量為對角元素的對角陣。

當(dāng)k=j=1時

當(dāng)k≠1且j≠1時

且

修正方程（簡寫為AX=B，A是N×N矩陣，X和B是N維列向量，N=3n+p+q）維數(shù)N隨系統(tǒng)規(guī)模迅速增大。例如，某省級1133節(jié)點系統(tǒng)N為4960，某區(qū)域2212節(jié)點系統(tǒng)N為9637。在CPU上反復(fù)求解如此高維數(shù)的修正方程耗費大量計算時間。經(jīng)計算表明，在CPU上修正方程的求解時間占全部優(yōu)化計算時間的96%以上。由于修正方程求解比較適合于并行計算（計算密集、高度并行），所以本文將該部分在GPU上執(zhí)行，其余部分在CPU上執(zhí)行。又由于修正方程系數(shù)矩陣在迭代過程中是不斷變化的，所以對修正方程的求解不采用LU分解法，而采用高斯消去法。計算機中采用單精度浮點數(shù)比采用雙精度浮點數(shù)舍入誤差大，而列主元高斯消去法能夠減少這種誤差的影響，可提高計算精度，保證其數(shù)值穩(wěn)定性。故此，本文采用列主元高斯消去法求解修正方程，并且采用單精度和雙精度浮點數(shù)運算加以對照。

3 GPU的硬件結(jié)構(gòu)和并行計算原理

3.1 GPU特點和開發(fā)工具

NVIDIA公司的GeForce GTX 285 GPU有30個多處理器，時鐘1.476GHz；每個多處理器能夠同時處理32個warp，每32個并行線程稱之為1個warp。顯存帶寬為158.796GB/s（giga bytes per second）、容量為1GB。同時支持單精度和雙精度運算，理論最高運算速度為1062.72Gflops （giga floating-point operations per second）。而目前主流的雙核處理器只有不超過25Gflops的浮點運算速度。造成CPU和GPU性能之間的差異主要原因是，GPU是一種針對圖形處理的專用處理器，它的特點就是控制部件相對較少，芯片上主要部分都為計算器所占據(jù)，片內(nèi)Cache也很少，因此其專用性決定了它的計算速度是很快的。

GPU是專門為圖形運算而設(shè)計的，是一種高度并行化的流式處理器，對所有的像素進行并行操作。即，相同或類似的運算會在海量的數(shù)據(jù)上重復(fù)運行，這恰恰符合單指令多數(shù)據(jù)流的概念。因此，GPU在通用計算方面適合用于諸如科學(xué)運算、數(shù)據(jù)分析、線性代數(shù)、流體模擬等需要大量重復(fù)的數(shù)據(jù)集運算和密集的內(nèi)存存取的應(yīng)用程序[6]。相比之下，CPU本質(zhì)上是一個標(biāo)量計算模型，計算單元偏少，主要針對復(fù)雜控制和低延遲，而非高帶寬優(yōu)化。

2006年11月，NVIDIA公司發(fā)布了一個在GPU上的開發(fā)工具，即統(tǒng)一計算設(shè)備構(gòu)架（Compute Unified Device Architecture，CUDA）。CUDA是一個面向通用計算和圖形計算的解決方案，是C語言的一種擴展。CUDA程序構(gòu)架分為兩部分：主機和設(shè)備；主機指的是CPU，設(shè)備指的是GPU。主程序由CPU來執(zhí)行，而當(dāng)遇到數(shù)據(jù)并行處理部分，CUDA就會將程序編譯為GPU能執(zhí)行的代碼，傳送到GPU執(zhí)行。這個并行處理程序在CUDA里稱為“核”程序。當(dāng)調(diào)用“核”程序時，m個線程并行地執(zhí)行相同指令，加速比達到m。GPU作為運行 C 語言程序主機的協(xié)同處理器，“核”在 GPU 上執(zhí)行，而程序的其他部分在CPU上執(zhí)行[14-15]，即串行代碼在CPU上執(zhí)行，并行代碼在GPU上執(zhí)行。

3.2 GPU的硬件結(jié)構(gòu)

CPU的硬件結(jié)構(gòu)如圖1所示。通過CUDA編譯時，GPU可以被視為能執(zhí)行高數(shù)量并行線程的計算設(shè)備，運行在主機上的并行數(shù)據(jù)和高密度計算應(yīng)用程序部分被卸載到這個設(shè)備上。設(shè)備被看作一組帶有片上共享內(nèi)存的多處理器，每個多處理器有8個標(biāo)量處理器、一個多線程指令單元和片上共享內(nèi)存。每個多處理器使用四種片上內(nèi)存：寄存器、共享內(nèi)存、只讀常量緩存和只讀紋理緩存。另外，每個多處理器還有2個特殊函數(shù)處理單元用來計算諸如正弦和余弦函數(shù)等特殊函數(shù)。主機和設(shè)備使用它們各自的DRAM：主機內(nèi)存和設(shè)備內(nèi)存（見圖1），并能相互復(fù)制數(shù)據(jù)。“核”程序只能在設(shè)備內(nèi)存上執(zhí)行。

圖1 GPU硬件結(jié)構(gòu)Fig.1 GPU hardware structure

每個多處理器使用單指令多線程架構(gòu)（Single-Instruction Multiple-Thread，SIMT）[16]：在任何時鐘周期內(nèi)，每個線程被映射到一個標(biāo)量處理器上，各自獨立地按照其指令地址和寄存器狀態(tài)執(zhí)行；SIMT按照每32個并行線程為單位（稱為1個warp）來創(chuàng)建和組織執(zhí)行各個線程。同一個warp中的線程按照相同的程序地址啟動執(zhí)行，之后卻可以按照不同的分支路徑獨立執(zhí)行。這樣，當(dāng)同一warp內(nèi)的32個線程執(zhí)行相同的路徑時，效率最高；而當(dāng)這32個線程遇到流控制指令（if，switch，do，for，while）發(fā)生分支而執(zhí)行不同路徑時，不同路徑線程將被序列化地執(zhí)行，效率下降。

3.3 線程調(diào)用層次結(jié)構(gòu)

CUDA將計算分層，如圖2所示。計算的最小單位是線程；線程塊是指線程的批處理，是一個多處理器上運行的最小單位。各線程塊包含相同數(shù)量的線程，它通過快速的共享內(nèi)存有效地分享數(shù)據(jù)，并且可以在指定的內(nèi)存訪問中同步線程地執(zhí)行。

圖2 線程塊柵格結(jié)構(gòu)Fig.2 Grid of thread blocks

執(zhí)行同一“核”程序的塊可以合成一批線程塊的柵格，因此通過一個“核”程序發(fā)送請求的線程總數(shù)是每個塊內(nèi)的線程數(shù)與塊數(shù)量的乘積。柵格是主機調(diào)用GPU的最小單位。通常調(diào)用一個“核”程序（名為kernel）的形式如下：

其中，Dg為指定柵格的維數(shù)和大小；Db為指定每個塊的維數(shù)和大小；A，B是“核”程序計算用參數(shù)。圖2中Dg為（8，8）二維柵格，共有64個塊；Db為（16，16）二維塊，每塊共有256個線程。這樣，“核”程序發(fā)送請求的線程總數(shù)為64×256個線程。

柵格中線程塊的數(shù)目應(yīng)當(dāng)大于多處理器的數(shù)量，使得每個多處理器都處于工作狀態(tài)。進而，為了使得在線程同步時處理器不會處于等待狀態(tài)，每個多處理器需要數(shù)倍數(shù)量的活動線程塊。同時，為了適應(yīng)設(shè)備升級的需要，應(yīng)當(dāng)為“核”程序配置數(shù)百個或數(shù)千個線程塊[17]。

塊內(nèi)線程數(shù)量應(yīng)當(dāng)是warp尺寸（即32）的倍數(shù)，避免在線程不足的warp上浪費計算能力，也有利于實現(xiàn)訪問內(nèi)存聯(lián)合，最好采用128～256[17]。

3.4 內(nèi)存體系結(jié)構(gòu)

CUDA 線程可在執(zhí)行過程中訪問多個存儲器空間的數(shù)據(jù)，如圖3所示。每個線程都有一個私有的本地內(nèi)存。每個線程塊都有一個共享內(nèi)存，該內(nèi)存對于塊內(nèi)的所有線程都是可見的，并且與塊具有相同的生命周期。所有塊內(nèi)的線程都可訪問相同的全局內(nèi)存。此外只讀常量緩存空間和只讀紋理緩存空間可由所有線程訪問。

圖3 內(nèi)存體系結(jié)構(gòu)Fig.3 Memory hierarchy

每個處理器有一組本地寄存器。當(dāng)存在從寄存器讀取數(shù)據(jù)依賴于對寄存器寫的結(jié)果時，或者當(dāng)半個warp的線程訪問寄存器出現(xiàn)沖突時，訪問寄存器需要延時[14]。除此之外，讀寫寄存器不花費額外時間。

由于共享內(nèi)存是片上的，所以它比本地和全局內(nèi)存空間更快；同時為了提高內(nèi)存帶寬，共享內(nèi)存空間被劃分為大小相等的內(nèi)存模塊，稱為bank[14,17]。當(dāng)訪問共享內(nèi)存的w個地址分別位于w個不同的bank時可以同時訪問，帶寬提高到原來的w倍；而當(dāng)內(nèi)存請求的多個地址位于同一個bank中，就存在bank沖突，訪問必須被序列化，帶寬降低。事實上，對于一個warp的前半或者后半16個線程訪問共享內(nèi)存不存在bank沖突時，訪問共享內(nèi)存和訪問寄存器一樣快。

本地和全局內(nèi)存作為設(shè)備內(nèi)存的讀寫區(qū)域沒有被緩存，相對片上內(nèi)存而言，具有較大的時間延遲、較低的帶寬。經(jīng)測試，本機主存到設(shè)備全局內(nèi)存的帶寬只有1.304GB/s、反向傳輸帶寬1.327GB/s、設(shè)備到設(shè)備全局內(nèi)存之間的傳輸帶寬只有0.3MB/s，遠遠低于其顯存帶寬，制約其性能的有效發(fā)揮。算法設(shè)計中應(yīng)當(dāng)盡量減少通過本地和全局內(nèi)存進行數(shù)據(jù)傳輸。另外，為了獲得較高的帶寬，訪問全局內(nèi)存要遵循特殊的訪問模式。半個warp的16個線程訪問地址連續(xù)的若干全局內(nèi)存空間，稱為聯(lián)合的全局內(nèi)存訪問模式。全局內(nèi)存訪問聯(lián)合時，一次或兩次讀、寫操作就可以完成訪問；反之，訪問被序列化，有效帶寬降低[17]。算法設(shè)計中應(yīng)當(dāng)首先保證全局內(nèi)存訪問被聯(lián)合。

由于紋理作為設(shè)備內(nèi)存的讀區(qū)域被緩存了，當(dāng)紋理緩存沒有被錯過時，一個紋理操作只消耗紋理緩存的一個讀操作，否則，消耗一個從設(shè)備內(nèi)存的讀操作。所以，通過紋理緩存讀取設(shè)備內(nèi)存相對于從全局或常量內(nèi)存讀取設(shè)備內(nèi)存有更高的帶寬[14-17]。另外，紋理拾取可以避免從全局內(nèi)存未聯(lián)合存取[17]，增加有效帶寬。

4 線性方程組單精度和雙精度算法

4.1 三種計算方案

本文在GPU上使用三種方案對線性方程組（AX=B）求解。方案一采用紋理緩存空間存取矩陣、向量，只支持單精度浮點數(shù)據(jù)格式；方案二、三采用global內(nèi)存空間存取矩陣、向量，方案二用單精度浮點數(shù)據(jù)格式，而方案三用雙精度浮點數(shù)據(jù)格式。

各方案采用填充技術(shù)存放矩陣A，在每行最后一個元素后面填充一些存儲空間使得每行有256的倍數(shù)個元素，從而半個warp的16個線程可以方便地訪問地址連續(xù)的若干全局內(nèi)存空間，使得global存儲訪問區(qū)域被聯(lián)合，提高訪問效率。圖4表示212×212維線性方程組第k+1次消去運算的情形。圖4中，A、B各元素分別用a、b加下標(biāo)表示，* 符號表示內(nèi)存填充區(qū)。本機最多能夠填充32 768個單精度浮點數(shù)、16 384個雙精度浮點數(shù)，從而本算法可以求解維數(shù)低于32 768單精度、維數(shù)低于163 84雙精度線性方程組。

圖4 解212×212線性方程時第k+1次消去運算Fig.4 The k+1th elimination for solving 212×212 linear equation

各方案用一個線程塊（圖4中索引號為Bx）負(fù)責(zé)矩陣一行元素的計算，使線程塊數(shù)量最大化，提高處理器的計算效率。塊內(nèi)線程數(shù)量采用256，在維數(shù)大于256維的線性方程組中用循環(huán)訪問矩陣A。每個線程（圖4中索引號為Tx）對應(yīng)計算矩陣的一個元素。這樣，求解212×212維線性方程組“核”程序線程塊總數(shù)為212，每個塊內(nèi)線程數(shù)為256。

“核”程序通過一個名為“紋理擷取”的設(shè)備函數(shù)讀取紋理內(nèi)存。“紋理擷取”函數(shù)的第一個參數(shù)指定一個叫“紋理參照”的對象，該函數(shù)通過紋理坐標(biāo)來擷取綁定到“紋理參照”的內(nèi)存。“紋理參照”有一些屬性，其一是紋理尺寸，它定義是否可以通過一元紋理坐標(biāo)指定紋理按一維數(shù)組尋址，或者通過二元紋理坐標(biāo)指定紋理按二維數(shù)組尋址，亦或通過三元紋理坐標(biāo)指定紋理按三維數(shù)組尋址；另外，定義輸入輸出數(shù)據(jù)類型和數(shù)據(jù)讀取模式等。在方案一中，系數(shù)矩陣A用單元素二元紋理坐標(biāo)尋址，向量B用單元素一元紋理坐標(biāo)尋址。由于GPU對二元紋理坐標(biāo)尋址進行了優(yōu)化，這樣就可以充分利用顯存的高帶寬[17]，同時利用二元紋理坐標(biāo)下的紋理元素和矩陣元素間的一一對應(yīng)關(guān)系，便于檢索。

由于在同一個“核”程序中，當(dāng)對綁定到紋理的global內(nèi)存地址修改后進行紋理擷取，將取出不可預(yù)知的數(shù)據(jù)[17]。這意味著，只有對已經(jīng)在前面“核”程序或者數(shù)據(jù)復(fù)制中更新過的global內(nèi)存地址進行紋理擷取才準(zhǔn)確。因此在方案一中，采用在global內(nèi)存中分配兩個區(qū)域存儲矩陣A及其中間結(jié)果。交替將計算結(jié)果輸出到其中一個區(qū)域。而其余向量不需要在同一個“核”程序中進行更新后擷取，只用一個區(qū)域即可。在方案二和方案三中，因為global內(nèi)存是可以在同一個“核”程序中讀寫，所以當(dāng)采用global存儲時，不必使用兩個區(qū)域，減少數(shù)據(jù)傳遞量。

為了利用共享內(nèi)存和寄存器的快速存取特性，所有“核”程序采用下面的編程模式：

（1）從設(shè)備內(nèi)存加載數(shù)據(jù)到共享內(nèi)存或寄存器。

（2）在共享內(nèi)存和寄存器上處理數(shù)據(jù)。

（3）與塊內(nèi)的其他線程同步，以便每條線程可以安全地讀取由不同線程寫入的共享內(nèi)存的存儲單元。

（4）把結(jié)果寫回到設(shè)備內(nèi)存中。

4.2 計算流程

方案一的算法流程如下：

（1）將A和B中各元素復(fù)制到GPU全局內(nèi)存；

（2）循環(huán)fork=0 toN-2：

①將源區(qū)域綁定到紋理；

②取出矩陣第k+1列向量，找出其中第k+1到N個元素（即圖4中akk到a(211)k）中絕對值最大的元素；

③在源區(qū)域進行行交換，對B進行元素交換；

④消去運算，輸出到目標(biāo)區(qū)域；

⑤將元素a（k，k）從目標(biāo)區(qū)域輸出到源區(qū)域；

⑥更換目標(biāo)區(qū)域和源區(qū)域；

（3）結(jié)束循環(huán)；

（4）進行標(biāo)準(zhǔn)化計算，即bk/akk（fork=0 toN-1），將結(jié)果輸出到主機，結(jié)束計算。

方案二、三只用一個global內(nèi)存區(qū)域存取矩陣A，其余的計算流程與方案一相同。

設(shè)步驟①中第k+1～N個元素用向量o表示。首先進行第一次比較，將向量o的前半部分和后半部分元素進行一一對應(yīng)比較，將絕對值較大的元素放入前半部分，后半部分保持不變；然后，將o的前半部分繼續(xù)分為前半、后半兩部分，按上面方法處理；繼續(xù)這一過程直到尋找出最大值。當(dāng)向量大小不是2的整數(shù)倍時，在后面補入零元素，使其成為2的整數(shù)倍。這樣對于N維向量只需要進行ceil(log2N)次比較就能找出最大值（ceil表示向上取整），而在CPU上需要進行N-1次比較，極大地提高了計算效率，N越大，效率越高。

步驟④將前代和回代過程合并在一起計算。即第k+1次迭代時，將矩陣A的右下部分和右上部分、向量B前部和后部同時進行消去計算（見圖4），共需更新（N-1）×（N-k）個元素。

4.3 計算效率

采用上述三種計算方案獲得的加速比（加速比為算法在CPU上執(zhí)行時間與在GPU上執(zhí)行時間的比值）如圖5所示。加速比隨著線性方程組維數(shù)的增大而不斷增大，達到一定維數(shù)時趨于飽和。單精度計算在方程組維數(shù)到達6000×6000時飽和；雙精度計算在方程組維數(shù)到達3000×3000時飽和。方案一最大加速比57.46，方案二最大加速比43.56，方案三最大加速比10.58。可見，方案二相比方案一計算速度下降不少，這是由于 global存儲速度明顯低于texture緩沖存取速度，方案三計算速度最低，這是由于GPU原本是為圖形處理而設(shè)計的，其單精度浮點運算速度是最快的。故此，雖然當(dāng)GPU計算能力達到1.3[18]時，CUDA支持雙精度浮點計算，在精度滿足要求時最好還是采用單精度浮點運算。

圖5 三種方案加速比對照Fig.5 Speedup comparison of three schems

從圖6中可以看到，由于當(dāng)同一warp中的線程發(fā)生分流時計算效率下降，在方案一中對一般線性方程組進行計算時，考慮零元素后會使加速比明顯下降。

圖6 考慮矩陣A零元素后對加速比的影響Fig.6 Effect on speedup after considering zero element in matrix A

分別采用Matlab CPU版本（簡稱Matlab）、C語言CPU版本（簡稱CPU版本）和方案一C語言GPU版本（簡稱GPU版本）求解線性方程組，并在圖7中給出了計算時間對比。當(dāng)方程維數(shù)大于256×256時，GPU版本的計算時間大于CPU版本的；當(dāng)方程維數(shù)大于580×580時，GPU版本的計算時間大于Matlab的。雖然Matlab采用相當(dāng)好的優(yōu)化計算，GPU版本對高維線性方程組的求解速度還是遠遠高于Matlab的。圖8給出了采用紋理擷取的GPU版本相對于Matlab的加速比，其曲線呈波浪形上升，在方程維數(shù)達到7680×7680以上時趨于飽和，最高加速比可達3.79。加速比會出現(xiàn)局部下降現(xiàn)象，這是由于GPU全局內(nèi)存的存取特性決定的。當(dāng)一個warp的前半或者后半16個線程存取聯(lián)合的全局內(nèi)存時，存取速度是最快的。GTX 285圖形處理器的紋理聯(lián)合區(qū)域和全局內(nèi)存聯(lián)合區(qū)域大小為256B。當(dāng)方程維數(shù)適合聯(lián)合內(nèi)存區(qū)域大小時，程序獲得最好的計算速度，反之，速度下降。

圖7 線性方程組計算效率比較Fig.7 Comparison of computational efficiency for solving linear equation

圖8 采用紋理擷取的GPU版本相對于Matlab的加速比Fig.8 Speedup of GPU version by texture fetching ralative to Matlab

文獻[9]采用高斯消去法和LU分解法基于GPU求解線性方程組，獲得一定的加速效果，但只能進行方程維數(shù)小于4096的單精度浮點運算。Matlab在本機只能求解維數(shù)低于8859的線性方程組，而本文算法則可以求解維數(shù)低于32 768單精度、維數(shù)低于16 384雙精度線性方程組，對大電網(wǎng)計算有較好的適應(yīng)性。

5 無功優(yōu)化計算性能比較

分別用上述三種方案進行無功優(yōu)化計算，以IEEE 118節(jié)點和538、1133和2212節(jié)點系統(tǒng)為例進行了測試，系統(tǒng)基本數(shù)據(jù)見表1。

表1 測試系統(tǒng)基本數(shù)據(jù)Tab.1 Basic dates of test systems

5.1 收斂特性

補償間隙 ζ 的變化反映了優(yōu)化過程中滿足不等式約束的情況。補償間隙 ζ 計算公式如下：

對上式進行誤差分析：

式中，e表示參數(shù)的誤差；λ1、λ2分別為向量x1、x2各元素的誤差系數(shù)。忽略各變量間誤差和誤差系數(shù)的不同，進行簡化處理：

式中，ex、λ為向量x1、x2近似誤差和近似誤差系數(shù)。由式（22）可見補償間隙誤差e(ζ)隨系統(tǒng)規(guī)模的增大而線性增大。為了消去系統(tǒng)規(guī)模對補償間隙的影響，定義修正補償間隙ζλ如下：

采用修正補償間隙ζλ做收斂判據(jù)，那么在ζλ＜ε1時輸出最優(yōu)解，結(jié)束計算。ε1為收斂精度。經(jīng)過對4個測試系統(tǒng)的計算，采用修正補償間隙做收斂判據(jù)能夠較好地適應(yīng)各種規(guī)模電網(wǎng)無功優(yōu)化計算的需要，使得收斂精度統(tǒng)一。本文收斂精度ε1取2×10-6，最大迭代次數(shù)K設(shè)為50。圖9所示為538和2212節(jié)點系統(tǒng)迭代中ζλ變化曲線，單精度和雙精度計算都可以很快地收斂，收斂曲線基本一致。

圖9 迭代中修正補償間隙的變化軌跡Fig.9 Trajectories of corrected complementary gap with iterations

5.2 加速比

按照Amdahl定律，將順序程序并行化后可以獲得的理論加速比β[18]

式中，φ是不可并行化程序所占比例；σ是并行處理器數(shù)量。考慮到CPU和GPU頻率差異的影響，上式應(yīng)當(dāng)修改為

式中，γ為CPU和GPU頻率比值。

4個系統(tǒng)的計算結(jié)果見表2。由表2可見，采用單精度時計算速度是采用雙精度時的4倍左右；最高加速比出現(xiàn)在單精度計算2212節(jié)點系統(tǒng)時，高達28.706，用時42min。但是，最高理論加速比仍然是最高加速比的4倍左右，主要原因有兩點：第一，GPU設(shè)備到設(shè)備全局內(nèi)存之間的傳輸帶寬只有0.3MB/s，遠低于顯存帶寬8GB/s，降低了程序的計算效率；第二，采用GPU編程不能夠很好地利用修正方程系數(shù)矩陣的高度稀疏性，以單精度為例進行對照見表3。由于CPU計算可以有效地利用修正方程的稀疏性，計算修正方程的加速比為計算一般線性方程加速比的一半。

表2 四個系統(tǒng)的無功優(yōu)化計算結(jié)果Tab.2 Results of reactive optimization in four power systems

表3 系數(shù)矩陣稀疏性對求解修正方程加速比的影響Tab.3 Effect of sparsity of coefficient matrix on speedup in solving the correction equation

6 結(jié)論

本文應(yīng)用GPU來并行實現(xiàn)內(nèi)嵌離散懲罰的非線性原對偶內(nèi)點算法，并應(yīng)用于四個測試系統(tǒng)的無功優(yōu)化計算，主要有以下優(yōu)點：

（1）充分利用了GPU的高速浮點運算和高帶寬，將線性修正方程的求解在GPU上實現(xiàn)，其余部分在CPU上執(zhí)行，極大地提高了整體計算速度。相對純CPU雙精度計算，單精度計算加速比達到近30，雙精度計算最高加速比超過6，遠高于多機并行計算速度。

（2）在GPU上可實現(xiàn)求解維數(shù)低于32768單精度、維數(shù)低于16384雙精度線性方程組，能夠適應(yīng)大電網(wǎng)計算需要。

[1] 程瑩,劉明波.含離散控制變量的大規(guī)模電力系統(tǒng)無功優(yōu)化[J].中國電機工程學(xué)報,2002,22(5): 54-60.Cheng Ying,Liu Mingbo.Reactive-power optimization of large-scale power systems with discrete control variables[J].Proceedings of the CSEE,2002,22(5):54-60.

[2] Liu M B,Tso S K,Cheng Y.An extended nonlinear primal-dual interior-point algorithm for reactive-power optimization of large-scale power systems with discrete control variables[J].IEEE Transactions on Power Systems,2002,17(4):982-991.

[3] 繆楠林,劉明波,趙維興.電力系統(tǒng)動態(tài)無功優(yōu)化并行算法及其實現(xiàn)[J].電工技術(shù)學(xué)報,2009,24(2):150-157.Miao Nanlin,Liu Mingbo,Zhao Weixin.Parallel algorithm of dynamic reactive power optimization and its implementation[J].Transactions of China Electrotechnical Society,2009,24(2): 150-157.

[4] 劉寶英,楊仁剛.采用輔助問題原理的多分區(qū)并行無功優(yōu)化算法[J].中國電機工程學(xué)報,2009,29(7):47-51.Liu Baoying,Yang Rengang.Multi-subarea parallel reactive power optimization based on APP[J].Proceedings of the CSEE,2009,29(7): 47-51.

[5] 程新功,厲吉文,曹立霞,等.基于電網(wǎng)分區(qū)的多目標(biāo)分布式并行無功優(yōu)化研究[J].中國電機工程學(xué)報,2003,23(10): 109-113.Cheng Xingong,Li Jiwen,Cao Lixia,et al.Multi-objective distributed parallel reactive power optimization based on subarea division of the power systems[J].Proceedings of the CSEE,2003,23(10):109-113.

[6] 吳恩華.圖形處理器用于通用計算的技術(shù)、現(xiàn)狀及其挑戰(zhàn)[J].軟件學(xué)報,2004,15(10): 1493-1503.Wu Enhua.State of the art and future challenge on general purpose computation by graphics processing unit[J].Journal of Software,2004,15(10):1493-1503.

[7] 吳恩華,柳有權(quán).基于圖形處理器(GPU)的通用計算 [J].計算機輔助設(shè)計與圖形學(xué)學(xué)報,2004,16(5):601-612.Wu Enhu,Liu Youquan.General purpose computation on GPU[J].Journal of Computer-Aided Design & Computer Graphics,2004,16(5): 601-612.

[8] John D Owens,Mike Houston,et al.GPU computing[J].Proceedings of the IEEE,2008,96(5):879-899.

[9] Nico Galoppo,Naga K Govindaraju.LU-GPU:efficient algorithms for solving dense linear systems on graphics hardware[C].Proceedings of the ACM/IEEE SC 2005 Conference,2005: 3-15.

[10] 姜珊珊.基于GPU的修正單純形方法的實現(xiàn)[D].吉林: 吉林大學(xué),2008.

[11] Vahid Jalili Marandi,Venkata Dinavahi.Large-scale transient stability simulation on graphics processing units[C].Power & Energy Society General Meeting,Calgary,2009: 1-6.

[12] Joseph Euzebe Tate,Thomas J Overbye.Contouring for power systems using graphical processing units[C].Proceedings of the 41st Annual Hawaii International Conference on System Sciences,Waikoloa,HI,USA,2008: 1351-1360.

[13] Pan Yongsheng,Ross Whitaker.Feasibility of GPU-assisted iterative image reconstruction for mobile C-arm CT[C].Proceedings of the SPIE-The International Society for Optical Engineering,Lake Buena Vista,FL,USA ,2009: 72585J (9 pp.).

[14] Nvidia.CUDA Programming Guide (Version 2.3)[M].July 2009.

[15] Ziyad S Hakura,Anoop Gupta.The design and analysis of a cache architecture for texture mapping[C].Proceedings of the 24th Annual International Symposium on Computer Architecture,1997: 108-120.

[16] Nvidia.CUDA Programming Guide (Version 2.1)[M].Dec 2008.

[17] Nvidia.CUDA C Programming Best Practices Guide[M].July 2009.

[18] Blythe D.Rise of the graphics processor[J].Proceedings of the IEEE,2008,96(5): 761-778.