張寧宇，高 山，趙 欣

（東南大學(xué) 電氣工程學(xué)院，江蘇 南京 210096）

0 引言

半定規(guī)劃 SDP（SemiDefinite Programming）是線性規(guī)劃的一種推廣，它是在滿足約束“對(duì)稱矩陣的仿射組合半正定”的條件下使線性函數(shù)極大（極?。┗膯栴}。目前為止，已有學(xué)者對(duì)SDP在電力系統(tǒng)優(yōu)化調(diào)度中的應(yīng)用作了初步研究，文獻(xiàn)[1]使用SDP算法來求解中期水火電力系統(tǒng)調(diào)度問題，得到了較好的結(jié)果。文獻(xiàn)[2-3]針對(duì)機(jī)組組合模型中存在0／1整數(shù)變量的特點(diǎn)，引入一種整數(shù)約束條件，同時(shí)將目標(biāo)函數(shù)和各種約束條件轉(zhuǎn)換成變量的二次形式，最后采用對(duì)偶變尺度算法對(duì)得到的SDP模型進(jìn)行求解。SDP算法與拉格朗日松弛法[4]和傳統(tǒng)分支定界法[5]將機(jī)組給合（UC）模型分割成上下層循環(huán)迭代求解不同，它屬于一種直接求解算法，尤其是在考慮網(wǎng)絡(luò)安全約束情況下，可將所有約束化成相應(yīng)的半定矩陣形式后使用內(nèi)點(diǎn)法求解。

內(nèi)點(diǎn)SDP存在以下2個(gè)問題：所需內(nèi)存較大，由于每個(gè)約束條件均對(duì)應(yīng)一個(gè)半定矩陣，隨著機(jī)組組合模型變量和規(guī)模的增加，所需的內(nèi)存空間急劇增加；內(nèi)點(diǎn)法每次迭代中求解大規(guī)模、稀疏、半正定線性方程組需要消耗大量時(shí)間。針對(duì)第1個(gè)問題，可利用半定矩陣對(duì)稱、稀疏的特點(diǎn)，采用稀疏存儲(chǔ)技術(shù)減少所需內(nèi)存。第2個(gè)問題歸納為解Ax=b所示的大型稀疏線性方程組，解法可分為直接法和迭代法2類。直接法主要是通過對(duì)系數(shù)矩陣進(jìn)行變換，如Gauss消元、Cholesky分解、QR分解等，將原方程組化為三角或三對(duì)角等容易求解的形式，然后通過回代或追趕等方法得到方程組的解。該方法一般能得到準(zhǔn)確解，但由于固有的前推回代的特點(diǎn)，使得很難實(shí)現(xiàn)并行化，因而所需的內(nèi)存與計(jì)算量均很大。迭代法可分為古典迭代法和Krylov迭代法兩大類，其中，古典迭代法有 Jacobi迭代法、Gauss-Seidel迭代法、SOR等。目前很少用古典迭代法直接求解大規(guī)模稀疏線性方程組，常結(jié)合Krylov迭代法來使用。Krylov子空間方法的主要思想是為各迭代步遞歸構(gòu)造殘差向量，即第n步的殘差向量rn通過系數(shù)矩陣A的某個(gè)多項(xiàng)式與第1個(gè)殘差向量r1相乘得到。迭代多項(xiàng)式的選取應(yīng)使所構(gòu)造的殘差向量在某種內(nèi)積意義下相互正交，從而保證某種極小性，達(dá)到快速收斂的目的。由于迭代法的存儲(chǔ)開銷極大減少，同時(shí)每步迭代中只包括向量與矩陣的乘法和加法或者向量?jī)?nèi)積運(yùn)算，便于并行加速，因此成為稀疏線性方程組并行算法研究的熱點(diǎn)。

作為顯卡上計(jì)算核心的圖形處理器GPU（Graphic Processing Unit），是一種用于密集數(shù)據(jù)計(jì)算的多核并行處理器，計(jì)算單元數(shù)量要遠(yuǎn)超過CPU。已有部分學(xué)者對(duì)GPU在電力系統(tǒng)中的應(yīng)用進(jìn)行了研究，并取得了一定的成果。文獻(xiàn)[6]把GPU應(yīng)用到潮流計(jì)算中，取得了一定的加速比。文獻(xiàn)[7-9]利用GPU實(shí)現(xiàn)了電力系統(tǒng)暫態(tài)仿真并行計(jì)算，相比傳統(tǒng)的串行算法取得了良好的加速比。

為提高內(nèi)點(diǎn)SDP算法中大型稀疏線性方程組的計(jì)算效率，本文提出了一種基于GPU的Krylov子空間并行算法。該并行算法針對(duì)系數(shù)矩陣稀疏、半正定的特點(diǎn)，采用預(yù)條件處理的擬最小殘差法（QMR法），并以矩陣分塊技術(shù)為基礎(chǔ)，在CSR（Compressed Sparse Row）存儲(chǔ)格式下由GPU實(shí)現(xiàn)了稀疏矩陣的Incomplete Cholesky預(yù)處理方法和所有迭代計(jì)算。實(shí)驗(yàn)分析表明，它是一種求解稀疏、對(duì)稱、半正定線性方程組的有效方法。對(duì)10～100機(jī)24時(shí)段6個(gè)不同的算例進(jìn)行仿真，結(jié)果表明：本文算法是一種十分有效的求解機(jī)組組合問題的算法，所達(dá)到的加速效果隨著算例規(guī)模的增加而更加明顯。

1 機(jī)組組合問題的SDP模型

1.1 SDP原理

SDP是線性規(guī)劃的推廣，是一種非光滑凸優(yōu)化問題。將現(xiàn)代內(nèi)點(diǎn)法應(yīng)用于SDP后，可保證求解諸多凸優(yōu)化問題時(shí)在多項(xiàng)式時(shí)間內(nèi)收斂于最優(yōu)解[10]。

SDP模型的標(biāo)準(zhǔn)形式如下。

其中，Ai、C、Z 為 n×n 實(shí)對(duì)稱矩陣，A=[A1A2…Am]T；y 為 m 維向量；b=[b1b2…bm]T；“≥”表示左側(cè)矩陣為半正定矩陣，即該矩陣的特征根均大于等于0；“·”為跡運(yùn)算符，即。

比較成熟的SDP內(nèi)點(diǎn)法有原始對(duì)偶內(nèi)點(diǎn)法及對(duì)偶變尺度法等，本文所采用的不可行內(nèi)點(diǎn)法屬于前一類，可把一個(gè)不在可行域內(nèi)的解作為初始解進(jìn)行迭代運(yùn)算，簡(jiǎn)化了初始化計(jì)算。

1.2 不可行內(nèi)點(diǎn)法

求解式（1）所示的SDP原問題等價(jià)于求解其對(duì)偶問題的對(duì)數(shù)障礙問題，其一階KKT最優(yōu)條件為：

其中，X、Z≥0，用 XZ=μI代替式（2）中第 3 式（μ 為障礙因子，隨著的值為相應(yīng)的最優(yōu)值），并利用泰勒級(jí)數(shù)進(jìn)行展開，得到以下線性系統(tǒng)：

對(duì)式（3）進(jìn)行求解便可得到X、y和Z的搜索方向（ΔX，Δy，ΔZ）。但是直接求解得到的 ΔX 不是對(duì)稱矩陣，需引入對(duì)稱化算子對(duì)式（3）中的第3式進(jìn)行處理，算子如下[11]：

其中，V為需進(jìn)行對(duì)稱化處理的矩陣；P表示不同的搜索方向，本文取 P=Z-1／2。

得到：

其中，E=P-TZ?sP，F(xiàn)=PX?sP-T，?s為對(duì)稱克勞內(nèi)特積運(yùn)算；svec表示將對(duì)稱矩陣轉(zhuǎn)化為向量的運(yùn)算，smat為逆運(yùn)算svec運(yùn)算的逆運(yùn)算，表示將向量轉(zhuǎn)化為對(duì)稱矩陣。消去式（5）中的ΔX和ΔZ可得：

其中，G=AE-1FAT，h=rp-AE-1（Rc-FRd）。求解線性方程組式（6）得到Δy后，便可代入式（5）求出ΔX和ΔZ，然后通過線性搜索得到迭代步長(zhǎng)，當(dāng)對(duì)偶間隙滿足指定精度時(shí)中止迭代。

1.3 機(jī)組組合問題的SDP模型

本文采用文獻(xiàn)[6]所示的機(jī)組組合模型，機(jī)組的啟動(dòng)費(fèi)用采用傳統(tǒng)的冷熱啟動(dòng)三段式模型，約束條件包括機(jī)組出力約束、最小啟停時(shí)間約束、機(jī)組爬坡約束、功率平衡約束和旋轉(zhuǎn)備用約束等。該模型本身并非凸優(yōu)化問題，原因是機(jī)組啟停變量為0／1整數(shù)變量，為此引入輔助變量Q2=1，同時(shí)將0／1整數(shù)變量約束用凸二次約束u2i，t-ui，tQ=0 代替，ui，t表示機(jī)組i在時(shí)段t的啟停狀態(tài)，這樣UC問題就可描述為一個(gè)凸優(yōu)化問題，進(jìn)而用內(nèi)點(diǎn)SDP法求解。關(guān)于機(jī)組組合SDP模型的具體介紹可見文獻(xiàn)[6]。

2 GPU

2.1 GPU結(jié)構(gòu)

近年來，隨著人們對(duì)計(jì)算機(jī)圖像顯示效果的要求越來越高，顯卡上核心處理器GPU的計(jì)算量和吞吐量也越來越大。與CPU注重邏輯控制不同，GPU主要負(fù)責(zé)大規(guī)模的密集型數(shù)據(jù)并行計(jì)算。

圖1 CPU和GPU的結(jié)構(gòu)設(shè)計(jì)Fig.1 Structure of CPU and GPU

如圖1所示，CPU采用復(fù)雜的控制邏輯，用指令來控制單線程可執(zhí)行程序的執(zhí)行，還采用大型緩存，既可以減少訪問復(fù)雜應(yīng)用程序的指令和數(shù)據(jù)時(shí)產(chǎn)生的延時(shí)，又能夠節(jié)約帶寬。而GPU包含了大量的計(jì)算單元，通過多線程技術(shù)來提升運(yùn)算速度與吞吐量。硬件充分利用因?yàn)榈却L問內(nèi)存而產(chǎn)生較長(zhǎng)延時(shí)的大量線程，減少了控制邏輯中需要執(zhí)行的線程，因此簡(jiǎn)化了邏輯控制和緩存單元。圖中Control表示邏輯控制單元，ALU表示運(yùn)算單元，Cache表示緩沖單元，DRAM表示內(nèi)存單元。

2.2 統(tǒng)一計(jì)算設(shè)備架構(gòu)編程

統(tǒng)一計(jì)算設(shè)備架構(gòu)CUDA（Compute Unified De-vice Architecture）是NVIDIA在2007年推出的一種用于GPU的編程模型。該編程模型不需要借助于圖形學(xué)應(yīng)用程序編程接口（API），采用比較容易掌握的類C語言進(jìn)行開發(fā)。

CUDA模型中的線程采用線程網(wǎng)格（grid）和線程塊（block）2級(jí)結(jié)構(gòu)。線程塊網(wǎng)絡(luò)和線程塊中可以分別有多個(gè)線程塊和線程（thread），且都有唯一的ID標(biāo)志。線程是最基本的計(jì)算單元，可以分配到單個(gè)數(shù)據(jù)的計(jì)算中。以向量A[n]和B[n]相加為例，如果該程序在CPU上執(zhí)行需要循環(huán)n次；當(dāng)采用CUDA編程時(shí)，可以分配n個(gè)線程，單個(gè)線程根據(jù)自身ID來執(zhí)行向量對(duì)應(yīng)位置上元素的和運(yùn)算，如A[i]+B[i]，由此實(shí)現(xiàn)了并行化計(jì)算。

3 基于GPU的預(yù)條件處理QMR法

原-對(duì)偶內(nèi)點(diǎn)法求解機(jī)組組合SDP模型時(shí)，求解式（6）所示的大型線性方程組為主要步驟，矩陣G的維數(shù)等于模型的約束條件數(shù)，例如10機(jī)24時(shí)段算例的約束條件數(shù)為1488，20機(jī)算例為2928，100機(jī)算例達(dá)到14448，且呈現(xiàn)高維數(shù)、稀疏、對(duì)稱、半正定的特點(diǎn)，使用直接法（如Cholesky分解）求解需要大量的存儲(chǔ)空間以及消耗很長(zhǎng)的時(shí)間，而Krylov子空間法作為20世紀(jì)90年代才完整提出的一類迭代法，具有存儲(chǔ)量少、計(jì)算量少且易于并行等優(yōu)點(diǎn)，非常適合于并行求解大型稀疏線性方程組，且結(jié)合預(yù)條件處理技術(shù)可獲得良好的收斂特性和較高的數(shù)值穩(wěn)定性，目前已是求解大型稀疏線性方程組的最主要方法。關(guān)于Krylov子空間法的基本原理詳見文獻(xiàn)[12]，這里不再贅述。

矩陣G具有部分特征根接近于0的特點(diǎn)，常用于求解對(duì)稱矩陣的CG（Conjugate Gradient）法可能會(huì)導(dǎo)致數(shù)值的不穩(wěn)定性，因此，本文采用預(yù)條件處理QMR法來求解大型稀疏線性方程組。

3.1 預(yù)條件處理的QMR法

設(shè)有n階線性方程組：

其中，J為n×n維的實(shí)對(duì)稱矩陣，k為n維向量，w為方程的解。

設(shè)w0為該方程組迭代計(jì)算的初值，算法流程如下。

a.初始化 w0，計(jì)算 r0=k-Jw0，v0=0，q0=M-1r0，ρ0=r0Tq0，n=1。

d.若 ρn-1＜ε，結(jié)束計(jì)算；否則，qn=un+βnqn-1，n=n+1，返回步驟 b。

由此可知，除去步驟c中un=M-1rn的計(jì)算，n次循環(huán)迭代的計(jì)算量主要包括n次矩陣向量乘法、3n次向量?jī)?nèi)積運(yùn)算和3n次向量加減法。上述運(yùn)算可通過并行稀疏矩陣向量乘法和并行向量運(yùn)算在GPU上實(shí)現(xiàn)并行運(yùn)算，故并行算法實(shí)現(xiàn)的難點(diǎn)在于使用GPU計(jì)算預(yù)處理矩陣M和實(shí)現(xiàn)un=M-1rn運(yùn)算。

3.2 基于矩陣分塊的Incomplete Cholesky預(yù)處理矩陣計(jì)算

預(yù)處理技術(shù)通過改變系數(shù)矩陣的條件數(shù)來達(dá)到加快迭代法收斂速度的目的[12]。目前為止，尚沒有一種適用于所有線性方程組的預(yù)處理技術(shù)，實(shí)際計(jì)算中應(yīng)根據(jù)系數(shù)矩陣的特點(diǎn)來設(shè)計(jì)相應(yīng)的預(yù)處理矩陣。文獻(xiàn)[14-15]根據(jù)Chebyshev多項(xiàng)式求出矩陣A的近似逆矩陣作為M-1，加快了收斂速度，但僅適用于強(qiáng)正則矩陣。文獻(xiàn)[16]基于最小化Frobenius范數(shù)‖I-AM‖的方法來設(shè)計(jì)預(yù)處理矩陣，取得了一定的效果。文獻(xiàn)[17]提出了一種對(duì)稱超松弛（SSOR）法預(yù)處理技術(shù)，通過矩陣相乘和相加運(yùn)算便可得到矩陣M。Incomplete Cholesky[18]對(duì)于正定型系數(shù)矩陣是一種有效的預(yù)處理方法，但在分解過程中前后行（列）元素之間是順序進(jìn)行的，不利于并行化計(jì)算。

基于矩陣分塊技術(shù)，本文提出一種Incomplete Cholesky并行預(yù)處理方法，在CSR稀疏矩陣存儲(chǔ)格式下使用GPU實(shí)現(xiàn)了預(yù)處理矩陣的并行計(jì)算，提高了計(jì)算效率。

在求解UC問題的半定規(guī)模過程中，矩陣G的結(jié)構(gòu)如圖2所示，考慮到對(duì)稱性，只對(duì)下三角矩陣部分進(jìn)行分析，陰影部分表示非0元素，對(duì)角線的元素可分為各個(gè)獨(dú)立的2×2矩陣，如Fi、Bi所示?？梢钥闯鼍仃嘑i之間完全獨(dú)立，可進(jìn)行并行Cholesky分解；分解完成后，剩余矩陣部分B中非0元素利用Fi的分解結(jié)果進(jìn)行Incomplete Cholesky分解，由于B中各行元素互相獨(dú)立，因此可實(shí)現(xiàn)并行化；計(jì)算完成后，可在矩陣B對(duì)角線上繼續(xù)尋找互相獨(dú)立的2×2矩陣進(jìn)行并行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é)束計(jì)算。

實(shí)際計(jì)算中，矩陣G以CSR格式存儲(chǔ)，分塊計(jì)算中只需遍歷非0元，有效減少了計(jì)算時(shí)間；此外，采用原-對(duì)偶內(nèi)點(diǎn)法計(jì)算時(shí)，經(jīng)過少數(shù)幾次迭代后，矩陣G的非0元的位置便恒定下來，因此，矩陣分塊信息不需每次都進(jìn)行更新。以10機(jī)系統(tǒng)為例，Incomplete Cholesky預(yù)處理矩陣計(jì)算的循環(huán)次數(shù)從分塊前的1488次降至分塊后的96次，具體的計(jì)算消耗時(shí)間及加速比將在4.2中給出。

3.3 基于GPU的預(yù)處理QMR并行算法

CUDA并行計(jì)算kernel以線程網(wǎng)格（grid）的形式組織，每個(gè)線程網(wǎng)格由若干個(gè)線程塊（block）組成，而每個(gè)線程塊又由若干個(gè)線程（thread）組成，實(shí)質(zhì)上，kernel是以block為單位執(zhí)行的。

使用GPU實(shí)現(xiàn)基于矩陣分塊的Incomplete Cholesky預(yù)處理矩陣計(jì)算如圖3所示，給每個(gè)相互獨(dú)立對(duì)角線上的 2×2 矩陣（F1、F2、F3、F4）分配相應(yīng)線程 塊 （block0、block1、block2、block3） 進(jìn) 行 Cholesky分解；分解完成后，給矩陣B分配線程塊（block4、block5、block6）對(duì)相應(yīng)行的非0元素進(jìn)行Incomplete Cholesky分解。其中，由于矩陣 F1、F2、F3、F4中計(jì)算量較小，線程塊 block0、block1、block2、block3 中線程數(shù)取 32；block4、block5、block6 線程設(shè)計(jì)為二維形式，每一行線程對(duì)應(yīng)于矩陣B中某一行的計(jì)算。

圖3 GPU中的任務(wù)劃分Fig.3 Mission decomposition in GPU

圖4 CPU和GPU的主要計(jì)算內(nèi)容Fig.4 Main calculation contents of CPU and GPU

圖4所示為QMR并行算法在CPU和GPU中的任務(wù)劃分，將可以實(shí)現(xiàn)并行化的計(jì)算放在GPU中，而邏輯控制功能和簡(jiǎn)單的代數(shù)運(yùn)算由CPU實(shí)現(xiàn)，既充分利用了GPU適合密集型并行運(yùn)算的優(yōu)點(diǎn)，又發(fā)揮了CPU具有邏輯控制功能的優(yōu)勢(shì)。值得注意的是，GPU任務(wù)中線性方程組求解是指QMR算法中的un=M-1rn運(yùn)算，求解時(shí)存在前推和回代2個(gè)過程，前推過程可以利用已得到的矩陣G的正向分塊信息實(shí)現(xiàn)并行化，而回代過程需要對(duì)矩陣G進(jìn)行一次反向的分塊計(jì)算并根據(jù)得到的信息并行計(jì)算即可。

4 算例分析

4.1 開發(fā)環(huán)境及算例介紹

本文以Microsoft Visual Studio 2008和CUDA平臺(tái)為開發(fā)環(huán)境，采用C語言編寫了基于GPU的Incomplete Cholesky預(yù)處理QMR并行算法，并在同一平臺(tái)上編寫了基于CPU的Incomplete Cholesky串行分解和Cholesky直接法求解線性方程組的程序。硬件平臺(tái)：CPU為Intel Core 2，主頻為2.4 GHz，內(nèi)存為10 GB；GPU為NVIDIA Tesla C2050，顯存頻率為 1147 MHz。

4.2 Incomplete Cholesky并行預(yù)處理矩陣的計(jì)算

衡量并行算法優(yōu)劣的指標(biāo)為所需存儲(chǔ)量和加速比。對(duì)文中所述的Incomplete Cholesky預(yù)處理QMR而言，GPU并行算法所需的存儲(chǔ)量明顯小于CPU串行算法，因此只對(duì)算法的加速比進(jìn)行分析。表1列出了對(duì)不同規(guī)模的矩陣進(jìn)行Incomplete Cholesky分解時(shí)，CPU串行算法和GPU并行算法消耗的時(shí)間?？梢钥闯?，矩陣維數(shù)為1488和2928時(shí)，Incom-plete Cholesky分解的串行算法消耗的時(shí)間要小于并行算法，這是由于計(jì)算量不是太大且數(shù)據(jù)在GPU內(nèi)存和CPU內(nèi)存間的通信需要消耗一定的時(shí)間造成的，隨著矩陣維數(shù)從2928增加到14448，系數(shù)矩陣的分塊信息沒有發(fā)生改變，而每次循環(huán)計(jì)算時(shí)增加的計(jì)算量可以通過增加線程塊的方式實(shí)現(xiàn)并行化，因此 GPU并行算法消耗的時(shí)間逐漸小于CPU算法，維數(shù)為14448時(shí)可取得1.420的并行加速比。

表1 不同規(guī)模矩陣的Incomplete Cholesky并行分解Tab.1 Incomplete Cholesky parallel decomposition for matrixes of different sizes

4.3 基于GPU的QMR預(yù)處理并行算法

QMR算法中的參數(shù)ε設(shè)為10-10，表2所示為Cholesky直接法和基于GPU的QMR預(yù)處理并行算法對(duì)不同規(guī)模線性方程組求解所需的時(shí)間。

表2 不同規(guī)模的線性方程組求解Tab.2 Calculation for linear equation sets of different sizes

與Incomplete Cholesky并行算法類似，當(dāng)維數(shù)為1488時(shí)，直接法的運(yùn)行時(shí)間要小于并行算法，隨著維數(shù)的增加，并行算法的速度優(yōu)勢(shì)逐漸體現(xiàn)出來，最終可取得7.45的并行加速比。在不同規(guī)模算例情況下，本文并行算法在經(jīng)過幾次迭代后可取得滿足計(jì)算精度的結(jié)果，表明了算法的正確性。表3為N=8688時(shí)，QMR并行算法運(yùn)行時(shí)間的分布。表中解線性方程表示使用前推和回代求解方程組；向量運(yùn)算包括向量間、向量與常數(shù)以及向量?jī)?nèi)積等計(jì)算；其他部分包括CPU和GPU之間的數(shù)據(jù)傳遞、CPU的邏輯判斷程序和CPU上代數(shù)計(jì)算等。雖然經(jīng)過矩陣分塊技術(shù)可以減少Incomplete Cholesky分解和使用前推回代法求解un=M-1rn的循環(huán)計(jì)算次數(shù)，然而其本身所具有的串行性質(zhì)，使得其求解時(shí)間占總運(yùn)行時(shí)間的絕大部分，而稀疏矩陣與向量乘法和向量相關(guān)運(yùn)行運(yùn)算在GPU可實(shí)現(xiàn)良好的并行化。

表3 QMR并行算法時(shí)間分布Tab.3 Time consumption distribution ofparallel QMR method

4.4 基于GPU的內(nèi)點(diǎn)法求解機(jī)組組合問題的SDP模型

為了驗(yàn)證并行算法的有效性，采用 10、20、40、60、80和100機(jī)24時(shí)段6個(gè)測(cè)試系統(tǒng)，20～100機(jī)算例通過對(duì)10機(jī)算例的擴(kuò)展得到，其發(fā)電機(jī)參數(shù)及各時(shí)段負(fù)荷數(shù)據(jù)見文獻(xiàn)[19]，所有機(jī)組均考慮爬坡約束，旋轉(zhuǎn)備用取系統(tǒng)總負(fù)荷的10%。

用SDP求解機(jī)組組合問題是一種直接求解法，不需要多重循環(huán)，也不同求解一系列子問題，可同時(shí)考慮所有的約束條件，由表4可見，通過30次左右循環(huán)計(jì)算，便可得到解結(jié)果。當(dāng)算例為10機(jī)系統(tǒng)時(shí)，GPU并行內(nèi)點(diǎn)算法的運(yùn)行時(shí)間要略大于直接法，隨著算例規(guī)模的增大，求解線性方程組式（6）的計(jì)算量逐漸增加，GPU并行算法的并行效率得到了體現(xiàn)，加速比從10機(jī)系統(tǒng)的0.95變大到100機(jī)系統(tǒng)的2.61。但是，由于SDP內(nèi)點(diǎn)法除求解線性方程組外的其他計(jì)算也需要消耗一定的時(shí)間，因此取得的加速比要小于表2中的加速比。

表4 不同算例的內(nèi)點(diǎn)法計(jì)算時(shí)間Tab.4 Time consumption of interior point method for different cases

建立機(jī)組組合模型時(shí)，機(jī)組啟動(dòng)費(fèi)用采用時(shí)間的指數(shù)函數(shù)或者冷熱啟動(dòng)三段式費(fèi)用更加符合實(shí)際情況，然而卻難以轉(zhuǎn)換成SDP所要求的凸優(yōu)化形式，因此，在建立UC問題的SDP模型時(shí)啟動(dòng)費(fèi)用采用固定值方式，在通過內(nèi)點(diǎn)法求解得到機(jī)組的最優(yōu)啟停狀態(tài)并進(jìn)行經(jīng)濟(jì)調(diào)度時(shí)，目標(biāo)函數(shù)中再加入指數(shù)或者三段式的啟動(dòng)費(fèi)用，這樣可以保證計(jì)算結(jié)果的準(zhǔn)確性。表5所示為采用本文算法對(duì)6個(gè)算例求解得到的運(yùn)行成本。表6為不考慮爬坡約束時(shí)本文算法與其他算法計(jì)算結(jié)果的對(duì)比，表中，LR、GA、EP、LRGA、GAUC分別指拉格朗日松弛法、遺傳算法、進(jìn)化算法、拉格朗日-遺傳混合算法、基于分類的遺傳算法。

可以看出，本文算法可以很好地處理爬坡約束，得到的運(yùn)行費(fèi)用略大于不考慮爬坡約束的情況。從表6可以看出，SDP的并行內(nèi)點(diǎn)法每次計(jì)算只需固定的迭代次數(shù)，隨機(jī)性較小，可獲得穩(wěn)定的計(jì)算結(jié)果，得到的費(fèi)用優(yōu)于部分算法，并可得到較好的近似最優(yōu)解。

表5 不同系統(tǒng)的運(yùn)行成本Tab.5 Operating cost of different systems

表6 各種算法計(jì)算結(jié)果的比較Tab.6 Comparison of calculating results among different algorithms

通過上面的分析，可得到以下結(jié)論。

a.基于矩陣分塊技術(shù)的Incomplete Cholesky并行預(yù)處理矩陣計(jì)算，可獲得一定的加速比。由于分解本身具有串行的性質(zhì)，可取得的加速比有限，同時(shí)，在矩陣維數(shù)較小且計(jì)算量不大情況下，計(jì)算時(shí)間反而略大于串行算法。

b.基于GPU的QMR預(yù)處理并行算法是一種有效的大型稀疏線性方程組求解方法，可通過少數(shù)幾次迭代獲得解結(jié)果。Incomplete Cholesky預(yù)處理矩陣和前推回代的并行計(jì)算時(shí)間占解方程總時(shí)間的絕大部分，但該預(yù)處理法很好地處理系數(shù)矩陣半正定的特點(diǎn)，且在維數(shù)較大時(shí)可獲得良好的加速比。

c.SDP的并行內(nèi)點(diǎn)法可提高算法的計(jì)算效率，但如1.1節(jié)所示，除式（6）所示的大型稀疏線性方程組求解外，包括其他大量的數(shù)據(jù)處理以及矩陣和向量間的基本運(yùn)算，仍需要消耗一定的時(shí)間，故算法整體的加速比小于QMR預(yù)處理并行算法求解線性方程組。在建立機(jī)組組合的SDP模型時(shí)，文中所采用的啟動(dòng)費(fèi)用的處理方式可以很好地解決啟動(dòng)費(fèi)用難以化成凸規(guī)劃形式的問題，提高了結(jié)果的準(zhǔn)確性。

5 結(jié)論

SDP作為一種很有前景的規(guī)劃方法，存儲(chǔ)量大和計(jì)算時(shí)間長(zhǎng)一直是制約其更廣泛應(yīng)用的瓶頸，文中所述的SDP并行內(nèi)點(diǎn)法可有效地減少程序的運(yùn)行時(shí)間，并通過10～100機(jī)6個(gè)系統(tǒng)的仿真表明了算法的有效性，這為SDP在UC問題中的進(jìn)一步應(yīng)用提供了可能，也為求解其他組合優(yōu)化問題帶來了新的思路和方法。