張寧宇,高 山,趙 欣
(東南大學(xué) 電氣工程學(xué)院,江蘇 南京 210096)
半定規(guī)劃 SDP(SemiDefinite Programming)是線性規(guī)劃的一種推廣,它是在滿足約束“對稱矩陣的仿射組合半正定”的條件下使線性函數(shù)極大(極?。┗膯栴}。目前為止,已有學(xué)者對SDP在電力系統(tǒng)優(yōu)化調(diào)度中的應(yīng)用作了初步研究,文獻[1]使用SDP算法來求解中期水火電力系統(tǒng)調(diào)度問題,得到了較好的結(jié)果。文獻[2-3]針對機組組合模型中存在0/1整數(shù)變量的特點,引入一種整數(shù)約束條件,同時將目標(biāo)函數(shù)和各種約束條件轉(zhuǎn)換成變量的二次形式,最后采用對偶變尺度算法對得到的SDP模型進行求解。SDP算法與拉格朗日松弛法[4]和傳統(tǒng)分支定界法[5]將機組給合(UC)模型分割成上下層循環(huán)迭代求解不同,它屬于一種直接求解算法,尤其是在考慮網(wǎng)絡(luò)安全約束情況下,可將所有約束化成相應(yīng)的半定矩陣形式后使用內(nèi)點法求解。
內(nèi)點SDP存在以下2個問題:所需內(nèi)存較大,由于每個約束條件均對應(yīng)一個半定矩陣,隨著機組組合模型變量和規(guī)模的增加,所需的內(nèi)存空間急劇增加;內(nèi)點法每次迭代中求解大規(guī)模、稀疏、半正定線性方程組需要消耗大量時間。針對第1個問題,可利用半定矩陣對稱、稀疏的特點,采用稀疏存儲技術(shù)減少所需內(nèi)存。第2個問題歸納為解Ax=b所示的大型稀疏線性方程組,解法可分為直接法和迭代法2類。直接法主要是通過對系數(shù)矩陣進行變換,如Gauss消元、Cholesky分解、QR分解等,將原方程組化為三角或三對角等容易求解的形式,然后通過回代或追趕等方法得到方程組的解。該方法一般能得到準(zhǔn)確解,但由于固有的前推回代的特點,使得很難實現(xiàn)并行化,因而所需的內(nèi)存與計算量均很大。迭代法可分為古典迭代法和Krylov迭代法兩大類,其中,古典迭代法有 Jacobi迭代法、Gauss-Seidel迭代法、SOR等。目前很少用古典迭代法直接求解大規(guī)模稀疏線性方程組,常結(jié)合Krylov迭代法來使用。Krylov子空間方法的主要思想是為各迭代步遞歸構(gòu)造殘差向量,即第n步的殘差向量rn通過系數(shù)矩陣A的某個多項式與第1個殘差向量r1相乘得到。迭代多項式的選取應(yīng)使所構(gòu)造的殘差向量在某種內(nèi)積意義下相互正交,從而保證某種極小性,達到快速收斂的目的。由于迭代法的存儲開銷極大減少,同時每步迭代中只包括向量與矩陣的乘法和加法或者向量內(nèi)積運算,便于并行加速,因此成為稀疏線性方程組并行算法研究的熱點。
作為顯卡上計算核心的圖形處理器GPU(Graphic Processing Unit),是一種用于密集數(shù)據(jù)計算的多核并行處理器,計算單元數(shù)量要遠超過CPU。已有部分學(xué)者對GPU在電力系統(tǒng)中的應(yīng)用進行了研究,并取得了一定的成果。文獻[6]把GPU應(yīng)用到潮流計算中,取得了一定的加速比。文獻[7-9]利用GPU實現(xiàn)了電力系統(tǒng)暫態(tài)仿真并行計算,相比傳統(tǒng)的串行算法取得了良好的加速比。
為提高內(nèi)點SDP算法中大型稀疏線性方程組的計算效率,本文提出了一種基于GPU的Krylov子空間并行算法。該并行算法針對系數(shù)矩陣稀疏、半正定的特點,采用預(yù)條件處理的擬最小殘差法(QMR法),并以矩陣分塊技術(shù)為基礎(chǔ),在CSR(Compressed Sparse Row)存儲格式下由GPU實現(xiàn)了稀疏矩陣的Incomplete Cholesky預(yù)處理方法和所有迭代計算。實驗分析表明,它是一種求解稀疏、對稱、半正定線性方程組的有效方法。對10~100機24時段6個不同的算例進行仿真,結(jié)果表明:本文算法是一種十分有效的求解機組組合問題的算法,所達到的加速效果隨著算例規(guī)模的增加而更加明顯。
SDP是線性規(guī)劃的推廣,是一種非光滑凸優(yōu)化問題。將現(xiàn)代內(nèi)點法應(yīng)用于SDP后,可保證求解諸多凸優(yōu)化問題時在多項式時間內(nèi)收斂于最優(yōu)解[10]。
SDP模型的標(biāo)準(zhǔn)形式如下。
其中,Ai、C、Z 為 n×n 實對稱矩陣,A=[A1A2…Am]T;y 為 m 維向量;b=[b1b2…bm]T;“≥”表示左側(cè)矩陣為半正定矩陣,即該矩陣的特征根均大于等于0;“·”為跡運算符,即。
比較成熟的SDP內(nèi)點法有原始對偶內(nèi)點法及對偶變尺度法等,本文所采用的不可行內(nèi)點法屬于前一類,可把一個不在可行域內(nèi)的解作為初始解進行迭代運算,簡化了初始化計算。
求解式(1)所示的SDP原問題等價于求解其對偶問題的對數(shù)障礙問題,其一階KKT最優(yōu)條件為:
其中,X、Z≥0,用 XZ=μI代替式(2)中第 3 式(μ 為障礙因子,隨著的值為相應(yīng)的最優(yōu)值),并利用泰勒級數(shù)進行展開,得到以下線性系統(tǒng):
對式(3)進行求解便可得到X、y和Z的搜索方向(ΔX,Δy,ΔZ)。但是直接求解得到的 ΔX 不是對稱矩陣,需引入對稱化算子對式(3)中的第3式進行處理,算子如下[11]:
其中,V為需進行對稱化處理的矩陣;P表示不同的搜索方向,本文取 P=Z-1/2。
得到:
其中,E=P-TZ?sP,F(xiàn)=PX?sP-T,?s為對稱克勞內(nèi)特積運算;svec表示將對稱矩陣轉(zhuǎn)化為向量的運算,smat為逆運算svec運算的逆運算,表示將向量轉(zhuǎn)化為對稱矩陣。消去式(5)中的ΔX和ΔZ可得:
其中,G=AE-1FAT,h=rp-AE-1(Rc-FRd)。求解線性方程組式(6)得到Δy后,便可代入式(5)求出ΔX和ΔZ,然后通過線性搜索得到迭代步長,當(dāng)對偶間隙滿足指定精度時中止迭代。
本文采用文獻[6]所示的機組組合模型,機組的啟動費用采用傳統(tǒng)的冷熱啟動三段式模型,約束條件包括機組出力約束、最小啟停時間約束、機組爬坡約束、功率平衡約束和旋轉(zhuǎn)備用約束等。該模型本身并非凸優(yōu)化問題,原因是機組啟停變量為0/1整數(shù)變量,為此引入輔助變量Q2=1,同時將0/1整數(shù)變量約束用凸二次約束u2i,t-ui,tQ=0 代替,ui,t表示機組i在時段t的啟停狀態(tài),這樣UC問題就可描述為一個凸優(yōu)化問題,進而用內(nèi)點SDP法求解。關(guān)于機組組合SDP模型的具體介紹可見文獻[6]。
近年來,隨著人們對計算機圖像顯示效果的要求越來越高,顯卡上核心處理器GPU的計算量和吞吐量也越來越大。與CPU注重邏輯控制不同,GPU主要負責(zé)大規(guī)模的密集型數(shù)據(jù)并行計算。
圖1 CPU和GPU的結(jié)構(gòu)設(shè)計Fig.1 Structure of CPU and GPU
如圖1所示,CPU采用復(fù)雜的控制邏輯,用指令來控制單線程可執(zhí)行程序的執(zhí)行,還采用大型緩存,既可以減少訪問復(fù)雜應(yīng)用程序的指令和數(shù)據(jù)時產(chǎn)生的延時,又能夠節(jié)約帶寬。而GPU包含了大量的計算單元,通過多線程技術(shù)來提升運算速度與吞吐量。硬件充分利用因為等待訪問內(nèi)存而產(chǎn)生較長延時的大量線程,減少了控制邏輯中需要執(zhí)行的線程,因此簡化了邏輯控制和緩存單元。圖中Control表示邏輯控制單元,ALU表示運算單元,Cache表示緩沖單元,DRAM表示內(nèi)存單元。
統(tǒng)一計算設(shè)備架構(gòu)CUDA(Compute Unified De-vice Architecture)是NVIDIA在2007年推出的一種用于GPU的編程模型。該編程模型不需要借助于圖形學(xué)應(yīng)用程序編程接口(API),采用比較容易掌握的類C語言進行開發(fā)。
CUDA模型中的線程采用線程網(wǎng)格(grid)和線程塊(block)2級結(jié)構(gòu)。線程塊網(wǎng)絡(luò)和線程塊中可以分別有多個線程塊和線程(thread),且都有唯一的ID標(biāo)志。線程是最基本的計算單元,可以分配到單個數(shù)據(jù)的計算中。以向量A[n]和B[n]相加為例,如果該程序在CPU上執(zhí)行需要循環(huán)n次;當(dāng)采用CUDA編程時,可以分配n個線程,單個線程根據(jù)自身ID來執(zhí)行向量對應(yīng)位置上元素的和運算,如A[i]+B[i],由此實現(xiàn)了并行化計算。
原-對偶內(nèi)點法求解機組組合SDP模型時,求解式(6)所示的大型線性方程組為主要步驟,矩陣G的維數(shù)等于模型的約束條件數(shù),例如10機24時段算例的約束條件數(shù)為1488,20機算例為2928,100機算例達到14448,且呈現(xiàn)高維數(shù)、稀疏、對稱、半正定的特點,使用直接法(如Cholesky分解)求解需要大量的存儲空間以及消耗很長的時間,而Krylov子空間法作為20世紀(jì)90年代才完整提出的一類迭代法,具有存儲量少、計算量少且易于并行等優(yōu)點,非常適合于并行求解大型稀疏線性方程組,且結(jié)合預(yù)條件處理技術(shù)可獲得良好的收斂特性和較高的數(shù)值穩(wěn)定性,目前已是求解大型稀疏線性方程組的最主要方法。關(guān)于Krylov子空間法的基本原理詳見文獻[12],這里不再贅述。
矩陣G具有部分特征根接近于0的特點,常用于求解對稱矩陣的CG(Conjugate Gradient)法可能會導(dǎo)致數(shù)值的不穩(wěn)定性,因此,本文采用預(yù)條件處理QMR法來求解大型稀疏線性方程組。
設(shè)有n階線性方程組:
其中,J為n×n維的實對稱矩陣,k為n維向量,w為方程的解。
設(shè)w0為該方程組迭代計算的初值,算法流程如下。
a.初始化 w0,計算 r0=k-Jw0,v0=0,q0=M-1r0,ρ0=r0Tq0,n=1。
d.若 ρn-1<ε,結(jié)束計算;否則,qn=un+βnqn-1,n=n+1,返回步驟 b。
由此可知,除去步驟c中un=M-1rn的計算,n次循環(huán)迭代的計算量主要包括n次矩陣向量乘法、3n次向量內(nèi)積運算和3n次向量加減法。上述運算可通過并行稀疏矩陣向量乘法和并行向量運算在GPU上實現(xiàn)并行運算,故并行算法實現(xiàn)的難點在于使用GPU計算預(yù)處理矩陣M和實現(xiàn)un=M-1rn運算。
預(yù)處理技術(shù)通過改變系數(shù)矩陣的條件數(shù)來達到加快迭代法收斂速度的目的[12]。目前為止,尚沒有一種適用于所有線性方程組的預(yù)處理技術(shù),實際計算中應(yīng)根據(jù)系數(shù)矩陣的特點來設(shè)計相應(yīng)的預(yù)處理矩陣。文獻[14-15]根據(jù)Chebyshev多項式求出矩陣A的近似逆矩陣作為M-1,加快了收斂速度,但僅適用于強正則矩陣。文獻[16]基于最小化Frobenius范數(shù)‖I-AM‖的方法來設(shè)計預(yù)處理矩陣,取得了一定的效果。文獻[17]提出了一種對稱超松弛(SSOR)法預(yù)處理技術(shù),通過矩陣相乘和相加運算便可得到矩陣M。Incomplete Cholesky[18]對于正定型系數(shù)矩陣是一種有效的預(yù)處理方法,但在分解過程中前后行(列)元素之間是順序進行的,不利于并行化計算。
基于矩陣分塊技術(shù),本文提出一種Incomplete Cholesky并行預(yù)處理方法,在CSR稀疏矩陣存儲格式下使用GPU實現(xiàn)了預(yù)處理矩陣的并行計算,提高了計算效率。
在求解UC問題的半定規(guī)模過程中,矩陣G的結(jié)構(gòu)如圖2所示,考慮到對稱性,只對下三角矩陣部分進行分析,陰影部分表示非0元素,對角線的元素可分為各個獨立的2×2矩陣,如Fi、Bi所示。可以看出矩陣Fi之間完全獨立,可進行并行Cholesky分解;分解完成后,剩余矩陣部分B中非0元素利用Fi的分解結(jié)果進行Incomplete Cholesky分解,由于B中各行元素互相獨立,因此可實現(xiàn)并行化;計算完成后,可在矩陣B對角線上繼續(xù)尋找互相獨立的2×2矩陣進行并行Cholesky分解,如此循環(huán)直至得到Incom-plete Cholesky預(yù)處理矩陣M。
圖2 系數(shù)矩陣GFig.2 Coefficient matrix G
矩陣G的分塊流程具體如下。
a.存放分塊信息的向量 Pos,flag=1,Pos[flag]=1,矩陣維數(shù)為 N,i取 2,j取 1。
b.j=Pos[flag]。
c.如果 G[i,j]為非 0 元且 i-j>2,flag=flag+1,Pos[flag]=i,轉(zhuǎn)步驟 e;否則轉(zhuǎn)步驟 d。
d.如果 j<i-1,j=j+1;否則,轉(zhuǎn)步驟 e。
e.如果 i<N,i=i+1,返回步驟 b;否則,結(jié)束計算。
實際計算中,矩陣G以CSR格式存儲,分塊計算中只需遍歷非0元,有效減少了計算時間;此外,采用原-對偶內(nèi)點法計算時,經(jīng)過少數(shù)幾次迭代后,矩陣G的非0元的位置便恒定下來,因此,矩陣分塊信息不需每次都進行更新。以10機系統(tǒng)為例,Incomplete Cholesky預(yù)處理矩陣計算的循環(huán)次數(shù)從分塊前的1488次降至分塊后的96次,具體的計算消耗時間及加速比將在4.2中給出。
CUDA并行計算kernel以線程網(wǎng)格(grid)的形式組織,每個線程網(wǎng)格由若干個線程塊(block)組成,而每個線程塊又由若干個線程(thread)組成,實質(zhì)上,kernel是以block為單位執(zhí)行的。
使用GPU實現(xiàn)基于矩陣分塊的Incomplete Cholesky預(yù)處理矩陣計算如圖3所示,給每個相互獨立對角線上的 2×2 矩陣(F1、F2、F3、F4)分配相應(yīng)線程 塊 (block0、block1、block2、block3) 進 行 Cholesky分解;分解完成后,給矩陣B分配線程塊(block4、block5、block6)對相應(yīng)行的非0元素進行Incomplete Cholesky分解。其中,由于矩陣 F1、F2、F3、F4中計算量較小,線程塊 block0、block1、block2、block3 中線程數(shù)取 32;block4、block5、block6 線程設(shè)計為二維形式,每一行線程對應(yīng)于矩陣B中某一行的計算。
圖3 GPU中的任務(wù)劃分Fig.3 Mission decomposition in GPU
圖4 CPU和GPU的主要計算內(nèi)容Fig.4 Main calculation contents of CPU and GPU
圖4所示為QMR并行算法在CPU和GPU中的任務(wù)劃分,將可以實現(xiàn)并行化的計算放在GPU中,而邏輯控制功能和簡單的代數(shù)運算由CPU實現(xiàn),既充分利用了GPU適合密集型并行運算的優(yōu)點,又發(fā)揮了CPU具有邏輯控制功能的優(yōu)勢。值得注意的是,GPU任務(wù)中線性方程組求解是指QMR算法中的un=M-1rn運算,求解時存在前推和回代2個過程,前推過程可以利用已得到的矩陣G的正向分塊信息實現(xiàn)并行化,而回代過程需要對矩陣G進行一次反向的分塊計算并根據(jù)得到的信息并行計算即可。
本文以Microsoft Visual Studio 2008和CUDA平臺為開發(fā)環(huán)境,采用C語言編寫了基于GPU的Incomplete Cholesky預(yù)處理QMR并行算法,并在同一平臺上編寫了基于CPU的Incomplete Cholesky串行分解和Cholesky直接法求解線性方程組的程序。硬件平臺:CPU為Intel Core 2,主頻為2.4 GHz,內(nèi)存為10 GB;GPU為NVIDIA Tesla C2050,顯存頻率為 1147 MHz。
衡量并行算法優(yōu)劣的指標(biāo)為所需存儲量和加速比。對文中所述的Incomplete Cholesky預(yù)處理QMR而言,GPU并行算法所需的存儲量明顯小于CPU串行算法,因此只對算法的加速比進行分析。表1列出了對不同規(guī)模的矩陣進行Incomplete Cholesky分解時,CPU串行算法和GPU并行算法消耗的時間。可以看出,矩陣維數(shù)為1488和2928時,Incom-plete Cholesky分解的串行算法消耗的時間要小于并行算法,這是由于計算量不是太大且數(shù)據(jù)在GPU內(nèi)存和CPU內(nèi)存間的通信需要消耗一定的時間造成的,隨著矩陣維數(shù)從2928增加到14448,系數(shù)矩陣的分塊信息沒有發(fā)生改變,而每次循環(huán)計算時增加的計算量可以通過增加線程塊的方式實現(xiàn)并行化,因此 GPU并行算法消耗的時間逐漸小于CPU算法,維數(shù)為14448時可取得1.420的并行加速比。
表1 不同規(guī)模矩陣的Incomplete Cholesky并行分解Tab.1 Incomplete Cholesky parallel decomposition for matrixes of different sizes
QMR算法中的參數(shù)ε設(shè)為10-10,表2所示為Cholesky直接法和基于GPU的QMR預(yù)處理并行算法對不同規(guī)模線性方程組求解所需的時間。
表2 不同規(guī)模的線性方程組求解Tab.2 Calculation for linear equation sets of different sizes
與Incomplete Cholesky并行算法類似,當(dāng)維數(shù)為1488時,直接法的運行時間要小于并行算法,隨著維數(shù)的增加,并行算法的速度優(yōu)勢逐漸體現(xiàn)出來,最終可取得7.45的并行加速比。在不同規(guī)模算例情況下,本文并行算法在經(jīng)過幾次迭代后可取得滿足計算精度的結(jié)果,表明了算法的正確性。表3為N=8688時,QMR并行算法運行時間的分布。表中解線性方程表示使用前推和回代求解方程組;向量運算包括向量間、向量與常數(shù)以及向量內(nèi)積等計算;其他部分包括CPU和GPU之間的數(shù)據(jù)傳遞、CPU的邏輯判斷程序和CPU上代數(shù)計算等。雖然經(jīng)過矩陣分塊技術(shù)可以減少Incomplete Cholesky分解和使用前推回代法求解un=M-1rn的循環(huán)計算次數(shù),然而其本身所具有的串行性質(zhì),使得其求解時間占總運行時間的絕大部分,而稀疏矩陣與向量乘法和向量相關(guān)運行運算在GPU可實現(xiàn)良好的并行化。
表3 QMR并行算法時間分布Tab.3 Time consumption distribution ofparallel QMR method
為了驗證并行算法的有效性,采用 10、20、40、60、80和100機24時段6個測試系統(tǒng),20~100機算例通過對10機算例的擴展得到,其發(fā)電機參數(shù)及各時段負荷數(shù)據(jù)見文獻[19],所有機組均考慮爬坡約束,旋轉(zhuǎn)備用取系統(tǒng)總負荷的10%。
用SDP求解機組組合問題是一種直接求解法,不需要多重循環(huán),也不同求解一系列子問題,可同時考慮所有的約束條件,由表4可見,通過30次左右循環(huán)計算,便可得到解結(jié)果。當(dāng)算例為10機系統(tǒng)時,GPU并行內(nèi)點算法的運行時間要略大于直接法,隨著算例規(guī)模的增大,求解線性方程組式(6)的計算量逐漸增加,GPU并行算法的并行效率得到了體現(xiàn),加速比從10機系統(tǒng)的0.95變大到100機系統(tǒng)的2.61。但是,由于SDP內(nèi)點法除求解線性方程組外的其他計算也需要消耗一定的時間,因此取得的加速比要小于表2中的加速比。
表4 不同算例的內(nèi)點法計算時間Tab.4 Time consumption of interior point method for different cases
建立機組組合模型時,機組啟動費用采用時間的指數(shù)函數(shù)或者冷熱啟動三段式費用更加符合實際情況,然而卻難以轉(zhuǎn)換成SDP所要求的凸優(yōu)化形式,因此,在建立UC問題的SDP模型時啟動費用采用固定值方式,在通過內(nèi)點法求解得到機組的最優(yōu)啟停狀態(tài)并進行經(jīng)濟調(diào)度時,目標(biāo)函數(shù)中再加入指數(shù)或者三段式的啟動費用,這樣可以保證計算結(jié)果的準(zhǔn)確性。表5所示為采用本文算法對6個算例求解得到的運行成本。表6為不考慮爬坡約束時本文算法與其他算法計算結(jié)果的對比,表中,LR、GA、EP、LRGA、GAUC分別指拉格朗日松弛法、遺傳算法、進化算法、拉格朗日-遺傳混合算法、基于分類的遺傳算法。
可以看出,本文算法可以很好地處理爬坡約束,得到的運行費用略大于不考慮爬坡約束的情況。從表6可以看出,SDP的并行內(nèi)點法每次計算只需固定的迭代次數(shù),隨機性較小,可獲得穩(wěn)定的計算結(jié)果,得到的費用優(yōu)于部分算法,并可得到較好的近似最優(yōu)解。
表5 不同系統(tǒng)的運行成本Tab.5 Operating cost of different systems
表6 各種算法計算結(jié)果的比較Tab.6 Comparison of calculating results among different algorithms
通過上面的分析,可得到以下結(jié)論。
a.基于矩陣分塊技術(shù)的Incomplete Cholesky并行預(yù)處理矩陣計算,可獲得一定的加速比。由于分解本身具有串行的性質(zhì),可取得的加速比有限,同時,在矩陣維數(shù)較小且計算量不大情況下,計算時間反而略大于串行算法。
b.基于GPU的QMR預(yù)處理并行算法是一種有效的大型稀疏線性方程組求解方法,可通過少數(shù)幾次迭代獲得解結(jié)果。Incomplete Cholesky預(yù)處理矩陣和前推回代的并行計算時間占解方程總時間的絕大部分,但該預(yù)處理法很好地處理系數(shù)矩陣半正定的特點,且在維數(shù)較大時可獲得良好的加速比。
c.SDP的并行內(nèi)點法可提高算法的計算效率,但如1.1節(jié)所示,除式(6)所示的大型稀疏線性方程組求解外,包括其他大量的數(shù)據(jù)處理以及矩陣和向量間的基本運算,仍需要消耗一定的時間,故算法整體的加速比小于QMR預(yù)處理并行算法求解線性方程組。在建立機組組合的SDP模型時,文中所采用的啟動費用的處理方式可以很好地解決啟動費用難以化成凸規(guī)劃形式的問題,提高了結(jié)果的準(zhǔn)確性。
SDP作為一種很有前景的規(guī)劃方法,存儲量大和計算時間長一直是制約其更廣泛應(yīng)用的瓶頸,文中所述的SDP并行內(nèi)點法可有效地減少程序的運行時間,并通過10~100機6個系統(tǒng)的仿真表明了算法的有效性,這為SDP在UC問題中的進一步應(yīng)用提供了可能,也為求解其他組合優(yōu)化問題帶來了新的思路和方法。