一種基于gpu的多波前潮流計算方法和系統(tǒng)的制作方法
【專利摘要】本發(fā)明公開了一種基于GPU的多波前潮流計算方法，該方法包括：1)波前矩陣鏈異步并發(fā)執(zhí)行；2)CPU-GPU間任務(wù)的分配；3)GPU上矩陣乘法、減法算法的優(yōu)化。波前矩陣鏈采用異步并發(fā)執(zhí)行的模式處理，GPU的資源得到了充分的利用；根據(jù)矩陣的大小決定采用CPU還是GPU處理，使得單個矩陣的處理時間最少，結(jié)合波前矩陣鏈的異步并發(fā)執(zhí)行，能盡可能的減少CPU、GPU的空閑等待時間；在GPU上對矩陣乘法、減法進行優(yōu)化，并適時的利用shared memory，使程序的性能達到最優(yōu)。上述3種方法的結(jié)合，可以顯著地提升分解矩陣的性能。
【專利說明】一種基于GPU的多波前潮流計算方法和系統(tǒng)

【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于計算機高性能計算領(lǐng)域，更具體地，涉及一種基于GPU的多波前潮流 計算方法和系統(tǒng)，該方法能夠在GPU上利用多波前算法高效求解潮流計算。

【背景技術(shù)】
[0002] 如今，各行各業(yè)都已離不開電，在一個大城市中，一次意外的斷電事故將會造成巨 大的損失，因此，對電網(wǎng)能否安全穩(wěn)定的運行進行預(yù)測顯得尤為重要。電力系統(tǒng)潮流計算是 對電力系統(tǒng)進行安全分析的基礎(chǔ)，其實質(zhì)是求解一組多元非線性方程組。牛頓-拉夫遜法 因收斂速度快的特性，是求解潮流方程最常用的方法，它通過引入雅可比矩陣將非線性方 程組的求解轉(zhuǎn)變?yōu)橄鄳?yīng)的線性方程組的求解。由于雅可比矩陣非常的稀疏，當(dāng)電網(wǎng)的規(guī)模 非常大時，如何對雅可比矩陣進行快速的分解是一個具有挑戰(zhàn)性的研宄課題。
[0003] 目前，求解上述線性方程組主要有兩種方案：⑴稀疏三角分解技術(shù)：利用電力系 統(tǒng)中方程組的稀疏特性，盡量減少不必要的計算以提高求解的效率，但該方案不太易并行 化；（2)多波前方法：將大的稀疏的雅可比矩陣轉(zhuǎn)化為一系列小的稠密矩陣（波前矩陣）， 然后處理這些稠密矩陣，該方案可以很好的并行化。由于稀疏三角分解技術(shù)的原理較簡單， 如今還是稀疏三角分析技術(shù)應(yīng)用的比較多，因為不易并行化的特點，也就無法與GPU很好 的結(jié)合。而多波前方法因并行性好，可以很好的發(fā)揮GPU的計算能力，但目前研宄的還比較 少。


【發(fā)明內(nèi)容】

[0004] 針對現(xiàn)有技術(shù)的以上缺陷或改進需求，本發(fā)明提供了一種基于GPU的多波前潮流 計算方法和系統(tǒng)，其目的在于解決現(xiàn)有方法中存在的并行性程度不高及GPU資源沒充分利 用的技術(shù)問題。
[0005] 為實現(xiàn)上述目的，按照本發(fā)明的一個方面，提供了一種基于GPU的多波前潮流計 算方法，包括以下步驟：
[0006] (I)CPU端讀取潮流計算的原始數(shù)據(jù)，對其進行處理得到線性方程組的系數(shù)矩陣J 和常數(shù)項b，并構(gòu)建波前矩陣鏈的拓撲關(guān)系圖TG;
[0007] (2)在CPU和GPU上分配計算資源處理線性方程組J·ΔX=b:在CPU上對波前 矩陣作預(yù)處理；在GPU上作矩陣-矩陣乘法、減法的數(shù)值運算；
[0008] (3)根據(jù)TG為每個CPU端線程分配一個待處理的波前矩陣鏈，之后每個CPU端線 程開始處理其對應(yīng)波前矩陣鏈中的波前矩陣；(4)-(9)表示第cthread-i個CPU線程處理 波前矩陣鏈的過程，CPU中各個線程并行執(zhí)行（4)-(9)步；
[0009] (4)cthread-i號線程判斷波前矩陣鏈中的波前矩陣是否都已被處理，如果沒有則 轉(zhuǎn)入步驟（6);否則轉(zhuǎn)入步驟（5);
[0010] (5)cthread-i號線程判斷TG中的所有波前矩陣鏈是否都已被處理，如果是則轉(zhuǎn) 入步驟（10);否則根據(jù)TG判斷目前是否有可以被處理的波前矩陣鏈j，如果有則將波前矩 陣鏈j分配給cthread-i，然后進入步驟（6)，否則再次轉(zhuǎn)（5)，進行循環(huán)查詢；
[0011] (6)cthread-i號線程從為其分配的波前矩陣鏈中取一個波前矩陣F處理，將F分 為4塊小矩陣：13、4、0、(3，其維度分別為1^*1^、1^*11、111*1^、111*11，1^、111、11為任意值，？的分解過程 如下所5，
[0012]

【權(quán)利要求】
1. 一種基于GPU的多波前潮流計算方法，其特征在于，所述方法包括以下步驟： (1) CPU端讀取潮流計算的原始數(shù)據(jù)，對其進行處理得到線性方程組的系數(shù)矩陣J和常 數(shù)項b，并構(gòu)建波前矩陣鏈的拓撲關(guān)系圖TG; (2) 在CPU和GPU上分配計算資源處理線性方程組J?AX=b:在CPU上對波前矩陣 作預(yù)處理；在GPU上作矩陣-矩陣乘法、減法的數(shù)值運算； (3) 根據(jù)TG為每個CPU端線程分配一個待處理的波前矩陣鏈，之后每個CPU端線程開 始處理其對應(yīng)波前矩陣鏈中的波前矩陣；(4)-(9)表示第cthread-i個CPU線程處理波前 矩陣鏈的過程，CPU中各個線程并行執(zhí)行（4)-(9)步； (4) cthread-i號線程判斷波前矩陣鏈中的波前矩陣是否都已被處理，如果沒有則轉(zhuǎn)入 步驟（6);否則轉(zhuǎn)入步驟（5); (5)Cthread-i號線程判斷TG中的所有波前矩陣鏈是否都已被處理，如果是則轉(zhuǎn)入步 驟（10);否則根據(jù)TG判斷目前是否有可以被處理的波前矩陣鏈j，如果有則將波前矩陣鏈 j分配給cthread-i，然后進入步驟（6)，否則再次轉(zhuǎn)（5)，進行循環(huán)查詢； (6) cthread-i號線程從為其分配的波前矩陣鏈中取一個波前矩陣F處理，將F分為4 塊小矩陣：13、4、0、(]，其維度分別為1^*1^、1^*11、111*1^、111*11，1^、111、11為任意值，？的分解過程如下 所示：
(7) 確定計算量為Cal=m*k*n,在cthread-i號線程中判斷Cal與閾值Threshold的 大小關(guān)系；若Cal〈 =Threshold,則轉(zhuǎn)入步驟（8)，否則轉(zhuǎn)入步驟（9); (8) cthread-i號線程執(zhí)行矩陣L2、U2的乘法運算，利用得到的結(jié)果與矩陣C做減法運 算，再將結(jié)果存儲在矩陣C中，之后轉(zhuǎn)入步驟（4); (9) cthread-i號線程將矩陣L2、U2、C存放在SQ的第sq-i個存儲單元中，然后GPU中 的第gblock-i個線程塊從第sq-i個存儲單元中取數(shù)據(jù)執(zhí)行，處理完成后將結(jié)果傳回給CPU 端的第cthread-i個線程；之后轉(zhuǎn)入步驟（4); (10) 終止GPU程序的運行； (11) 利用分解得到的矩陣L、U進行前代回代計算求得AX，然后利用AX去更新X，之 后判斷AX中絕對值最大的|AX」是否小于10,如果小于則轉(zhuǎn)入步驟（12);否則轉(zhuǎn)入步 驟⑶； (12) 根據(jù)所求得的X去更新潮流計算中的功率變量。
2. 根據(jù)權(quán)利要求1所述的多波前潮流計算方法，其特征在于，所述步驟（1)包括以下子 步驟： (1-1)讀取潮流計算的原始數(shù)據(jù)，根據(jù)電路理論中的基爾霍夫第一定律得到非線性方 程組I=YV，I、V、Y分別為電流、電壓、導(dǎo)納矩陣； (1-2)利用牛頓-拉夫遜法通過求導(dǎo)將非線性方程組轉(zhuǎn)化為線性方程組，轉(zhuǎn)化的公式
(1-3)對J作行列變換，之后根據(jù)多波前算法理論對J進行符號分解，得到消去樹T、每 個波前矩陣的信息、波前矩陣鏈的結(jié)構(gòu)及波前矩陣之間的依賴關(guān)系； (1-4)根據(jù)波前矩陣鏈的結(jié)構(gòu)、波前矩陣之間的依賴關(guān)系得到波前矩陣鏈的拓撲關(guān)系 圖TG。
3. 根據(jù)權(quán)利要求1或2所述的多波前潮流計算方法，其特征在于，所述步驟（2)具體包 括： (2-1)在CPU端分配n個線程，每個線程處理一個波前矩陣鏈； (2-2)在GPU端分配n個線程塊，并在GPU顯存的全局內(nèi)存上分配大小為n*unit的存 儲隊列SQ，其中每段的長度為unit，共有n段，n段是連續(xù)的； (2-3) CPU中的第cthread-i個線程與GPU中的第gblock-i個線程塊、存儲隊列中的第sq_i個存儲單元--對應(yīng)，其中 1〈 =cthread-i、gblock-i、sq_i〈=n; (2-4)利用GPU顯存上的紋理內(nèi)存處理SQ的同步操作。
4. 根據(jù)權(quán)利要求1或2所述的多波前潮流計算方法，其特征在于，所述步驟（6)具體包 括： (6-1)利用LU分解法將矩陣B分解為上三角矩陣L1，其維度為k*k，和下三角矩陣U1， 其維度為k*k; (6-2)利用（6-1)得到的1^、％，對矩陣D、A進行矩陣-矩陣乘法運算，得到矩陣L2，其 維度為m*k，以及U2，其維度為k*n; (6-3)利用（6-2)得到的L2、U2，通過矩陣-矩陣的乘法與減法，對矩陣C更新，再將結(jié) 果存儲在C中，C=C-L2U2。
5. 根據(jù)權(quán)利要求1或2所述的多波前潮流計算方法，其特征在于，所述步驟（7)中的閾 值Threshold具體為： 設(shè)CPU處理矩陣所花費的時間為Tepu，GPU處理矩陣所花費的時間為，下式表示了T^、U具體的值
&=mnk，N2=mk+nk+2mn，1表示做矩陣運算操作的次數(shù)，1表示數(shù)據(jù)傳輸?shù)拇笮。籥agpu分別為在CPU、GPU上執(zhí)行矩陣-矩陣乘法、減法操作的平均性能；0為CPU-GPU 間傳輸矩陣的平均帶寬
6. 根據(jù)權(quán)利要求4所述的多波前潮流計算方法，其特征在于，所述步驟（9)包括以下子 步驟： (9-l)cthread-i號線程判斷SQ的sq-i個存儲單元是否為空，如不為空則再次轉(zhuǎn) (9-1)，進行循環(huán)查詢；否則轉(zhuǎn)（9-2); (9-2)cthread_i號線程將L2、U2、C拷貝到SQ的第sq-i個存儲單元中，之后通知GPU的第gblock-i個線程塊有數(shù)據(jù)需要處理； (9-3)當(dāng)GPU中第gblock-i個線程塊的線程處于空閑狀態(tài)時，判斷SQ的第sq-i個存 儲單元是否為空，如果為空，則轉(zhuǎn)（9-3)，進行循環(huán)查詢；否則轉(zhuǎn)（9-4); (9-4)判斷Cal的大小，如果Cal>=Gthreshold(閾值），則將第sq-i個存儲單元中 的數(shù)據(jù)轉(zhuǎn)存到GPU的共享內(nèi)存中，然后執(zhí)行矩陣乘法、減法運算；否則直接對第sq-i個存儲 單元的數(shù)據(jù)執(zhí)行矩陣乘法、減法運算；處理完后，將結(jié)果傳回給CPU端的cthread-i號線程。
7. -種基于GPU的多波前潮流計算系統(tǒng)，其特征在于，包括： 第一模塊，用于在CPU端讀取潮流計算的原始數(shù)據(jù)，對其進行處理得到線性方程組的 系數(shù)矩陣J和常數(shù)項b，并構(gòu)建波前矩陣鏈的拓撲關(guān)系圖TG; 第二模塊，用于在CPU和GPU上分配計算資源處理線性方程組J?A X=b :在CPU上 對波前矩陣作預(yù)處理；在GPU上作矩陣-矩陣乘法、減法的數(shù)值運算； 第三模塊，用于根據(jù)TG為每個CPU端線程分配一個待處理的波前矩陣鏈，之后每個CPU端線程開始處理其對應(yīng)波前矩陣鏈中的波前矩陣； 第四模塊，用于使cthread-i號線程判斷波前矩陣鏈中的波前矩陣是否都已被處理， 如果沒有則轉(zhuǎn)入第六模塊；否則轉(zhuǎn)入第五模塊； 第五模塊，用于使cthread-i號線程判斷TG中的所有波前矩陣鏈是否都已被處理，如 果是則轉(zhuǎn)入第十模塊；否則根據(jù)TG判斷目前是否有可以被處理的波前矩陣鏈j，如果有則 將波前矩陣鏈j分配給cthread-i，然后進入第六模塊，否則再次轉(zhuǎn)第五模塊，進行循環(huán)查 詢； 第六模塊，用于使cthread-i號線程從為其分配的波前矩陣鏈中取一個波前矩陣F處 理，將F分為4塊小矩陣：B、A、D、C，其維度分別為k*k、k*n、m*k、m*n，k、m、n為任意值，F(xiàn) 的分解過程如下所示：
第七模塊，用于確定計算量為Cal=m*k*n,在cthread-i號線程中判斷Cal與閾值Threshold的大小關(guān)系；若Cal〈 =Threshold,則轉(zhuǎn)入第八模塊，否則轉(zhuǎn)入第九模塊； 第八模塊，用于使cthread-i號線程執(zhí)行矩陣L2、U2的乘法運算，利用得到的結(jié)果與矩 陣C做減法運算，再將結(jié)果存儲在矩陣C中，之后轉(zhuǎn)入第四模塊； 第九模塊，用于使cthread-i號線程將矩陣L2、U2、C存放在SQ的第sq-i個存儲單元 中，然后GPU中的第gblock-i個線程塊從第sq-i個存儲單元中取數(shù)據(jù)執(zhí)行，處理完成后將 結(jié)果傳回給CPU端的第cthread-i個線程；之后轉(zhuǎn)入第四模塊； 第十模塊，用于終止GPU程序的運行； 第十一模塊，用于利用分解得到的矩陣L、U進行前代回代計算求得AX，然后利用AX去更新X，之后判斷AX中絕對值最大的|AX」是否小于10,如果小于則轉(zhuǎn)入第十二模 塊；否則轉(zhuǎn)入第三模塊； 第十二模塊，用于根據(jù)所求得的X去更新潮流計算中的功率變量。
【文檔編號】G06F9/50GK104484234SQ201410670758
【公開日】2015年4月1日 申請日期:2014年11月21日 優(yōu)先權(quán)日:2014年11月21日
【發(fā)明者】徐得超, 陳勇, 王偉, 江涵, 鄭然
申請人:中國電力科學(xué)研究院, 華中科技大學(xué)