黃玉龍 劉明波

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

1 引言

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

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

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

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

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

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

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

引入對(duì)數(shù)壁壘函數(shù)消去松弛變量的非負(fù)性約束，對(duì)等式約束引入拉格朗日乘子，并引入正曲率二次罰函數(shù)來(lái)處理QC和ΤB這類(lèi)離散變量[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)性條件可得滿(mǎn)足最優(yōu)性必要條件的非線性方程組，然后線性化得到如下的簡(jiǎn)化修正方程[1-2]：

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

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

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

且

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

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

3.1 GPU特點(diǎn)和開(kāi)發(fā)工具

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.1 三種計(jì)算方案

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

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

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

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

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

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

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

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

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

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

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

4.2 計(jì)算流程

方案一的算法流程如下：

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

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

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

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

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

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

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

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

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

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

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

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

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

4.3 計(jì)算效率

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

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

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

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

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

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

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

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

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

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

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

5.1 收斂特性

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

對(duì)上式進(jìn)行誤差分析：

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

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

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

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

5.2 加速比

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

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

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

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

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

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

6 結(jié)論

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

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

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

[1] 程瑩,劉明波.含離散控制變量的大規(guī)模電力系統(tǒng)無(wú)功優(yōu)化[J].中國(guó)電機(jī)工程學(xué)報(bào),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)動(dòng)態(tài)無(wú)功優(yōu)化并行算法及其實(shí)現(xiàn)[J].電工技術(shù)學(xué)報(bào),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] 劉寶英,楊仁剛.采用輔助問(wèn)題原理的多分區(qū)并行無(wú)功優(yōu)化算法[J].中國(guó)電機(jī)工程學(xué)報(bào),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)分布式并行無(wú)功優(yōu)化研究[J].中國(guó)電機(jī)工程學(xué)報(bào),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] 吳恩華.圖形處理器用于通用計(jì)算的技術(shù)、現(xiàn)狀及其挑戰(zhàn)[J].軟件學(xué)報(bào),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ì)算 [J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),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的修正單純形方法的實(shí)現(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.