配電網(wǎng)的回路阻抗法潮流計算與仿真可視化

王文舉，李光耀

（同濟大學 電子與信息工程學院，上海 201804）

潮流計算是電力網(wǎng)絡運行分析的基礎，可為電力系統(tǒng)運行人員進行電力系統(tǒng)網(wǎng)絡重構(gòu)、故障處理、無功優(yōu)化和狀態(tài)估計提供參考依據(jù)，在電網(wǎng)調(diào)度、運行分析、操作模擬和設計規(guī)劃中都發(fā)揮著重要的作用［1－2］.

因配電網(wǎng)在結(jié)構(gòu)上具有閉環(huán)結(jié)構(gòu)、開環(huán)運行的特性，穩(wěn)態(tài)運行時多呈輻射狀，在發(fā)生故障或倒換負荷時則呈短時環(huán)網(wǎng)運行結(jié)構(gòu)；且具有在線路參數(shù)上R／X（電阻／阻抗）較大，網(wǎng)絡的PQ（P 為有功率，Q為無功率）節(jié)點多、PV（V為電壓幅值）節(jié)點較少等特點，使得傳統(tǒng)的、面向環(huán)狀結(jié)構(gòu)高壓輸電網(wǎng)進行潮流計算分析的牛頓拉夫遜算法［3］、快速解耦［4］等方法出現(xiàn)病態(tài)收斂，不再有效.為此，許多學者從上世紀90年代開始對配電網(wǎng)潮流計算進行了深入研究：一方面針對配電網(wǎng)結(jié)構(gòu)特性對牛頓法［5］、快速解耦法［6］改進，但造成求解復雜化，失去了原有算法的收斂性和穩(wěn)定性；另一方面陸續(xù)提出前推回代類算法［7－8］、回路阻抗法［9］等新配電網(wǎng)潮流計算方法.

前推回代類算法具有線性收斂、魯棒性好、求解速度快、內(nèi)存占用少等優(yōu)異特性，但應用于環(huán)狀或含有PV節(jié)點的配電網(wǎng)潮流計算分析時，需進行特殊處理，從而增加了求解難度，其通用性、收斂性不強.

與上述算法相比，Goswami［9］等人提出的回路阻抗法，具有網(wǎng)孔處理能力強、收斂性好、穩(wěn)定性高、適用范圍廣等優(yōu)點.但在計算機求解過程中存在下述問題：① 為利用稀疏矩陣技術簡化回路阻抗矩陣的LU（lower－upper）分解，文獻［9－10］將配電網(wǎng)樹狀結(jié)構(gòu)轉(zhuǎn)化成標準二叉樹形式優(yōu)化節(jié)點編號，加大了算法復雜度，限定了中間節(jié)點狀態(tài)的求解順序，改變了配電網(wǎng)原有的網(wǎng)絡結(jié)構(gòu)，不便于網(wǎng)絡重構(gòu)；文獻［11］提出以鄰接表與二叉樹結(jié)合的方式優(yōu)化節(jié)點編號，便于網(wǎng)絡重構(gòu)，但也無法進行有環(huán)配電網(wǎng)阻抗矩陣的自動化生成.②CPU對阻抗矩陣進行LU分解時，文獻［12－13］采用一定的技巧，對對角線元素和相鄰元素進行了優(yōu)化，但因阻抗矩陣為對稱滿秩矩陣，在計算過程中矩陣元素仍將占用大量內(nèi)存，計算速度會隨著節(jié)點數(shù)目的大量增加而大幅下降.

此外，電網(wǎng)潮流計算仿真結(jié)果通常以文本或二維圖表的形式給出.雖然近年來已出現(xiàn)基于物理拓撲結(jié)構(gòu)的二維餅圖、三維柱狀圖的表征形式［14－15］，但一般用點狀或線狀等抽象符號表達電力設備，無法客觀、真實地反映與刻畫復雜的三維空間對象本身結(jié)構(gòu)和相互之間的關聯(lián)關系，也缺乏動態(tài)處理和時空分析的能力，存在很大的局限性.

針對上述工作的不足，筆者基于配電網(wǎng)回路阻抗法原理及配電網(wǎng)接線特點，提出了一種配電網(wǎng)潮流計算和仿真信息可視化的新方法，以快速、有效地模擬仿真配電網(wǎng)運行狀況：采用十字鏈表對配電網(wǎng)物理拓撲結(jié)構(gòu)進行存儲，定義反映各節(jié)點之間電流流向關系的關聯(lián)矩陣，在兩者的基礎上，設計了阻抗矩陣生成算法，用于有環(huán)、無環(huán)配電網(wǎng)阻抗矩陣中各元素的自動化計算，可以任意順序計算各支路電流、內(nèi)節(jié)點電壓；提出了基于right－looking LU分解法的并行高斯消去算法，利用GPU（graphic processing unit）對回路阻抗矩陣快速精確求解；結(jié)合GIS和虛擬現(xiàn)實技術，實現(xiàn)了仿真運算結(jié)果信息在三維GIS（geographic information system）系統(tǒng)中的可視化展示.

1 回路阻抗法求解原理

配電網(wǎng)的基本單元是饋線，每條饋線以輻射型網(wǎng)絡連接若干臺配電變壓器.如果將節(jié)點負荷用恒定阻抗表示，且不考慮配電線對地電容，并將中間變壓器等值化簡（其數(shù)學模型網(wǎng)絡結(jié)構(gòu)如圖1所示），則從饋線根節(jié)點起到每1個負荷節(jié)點都形成1條回路.以每條回路電流為所求解變量，根據(jù)基爾霍夫電壓定律，具有m個節(jié)點、n個負荷節(jié)點的配電網(wǎng)回路電流方程組可寫為

圖1 配電網(wǎng)數(shù)學模型Fig.1 Mathematic model

式中：VS為根節(jié)點電壓；Ii為流經(jīng)負荷節(jié)點vi的第i條回路電流；Zi，i為vi與VS之間的支路阻抗和加上負荷節(jié)點vi自身的負荷阻抗Zi，稱為負荷節(jié)點i的自阻抗；Zi，j為負荷節(jié)點vi和vj到VS的共同支路阻抗相加之和，稱為vi和vj的互阻抗［9］.

當配電網(wǎng)只有一個vi且在該處形成閉環(huán)結(jié)構(gòu)時，其環(huán)路如圖2a所示.設ip1，ip2為獨立的2條回路（p1，p2）的電流，Zcom為2條回路的共同支路阻抗，Zload為vi的負荷阻抗，根據(jù)KVL定律（基爾霍夫電壓定律，Kinhhoffs Voltage Law），該閉環(huán)結(jié)構(gòu)電路回路方程可寫為

此時，需將該環(huán)在vi處解開分成2條獨立的支路（圖2b）求解.整理上式，寫成以ip1，ip2為求解變量的方程組如下：

圖2 弱環(huán)節(jié)點電流計算Fig.2 Equivalent loop of meshed network

當配電網(wǎng)有多負荷節(jié)點在負荷節(jié)點vi處形成弱環(huán)時，由式（1），（3）可推導出下式：

式中：Zvi，vr＋Zvi可簡記為Zvi，vr；Zvr，vi＋Zvi可簡記為Zvr，vi，且兩者值相同；vr為vi破環(huán)后對應的虛擬節(jié)點.這樣，回路阻抗方程組增加1階.同理，若配電網(wǎng)在多個負荷節(jié)點形成多個環(huán)網(wǎng)，均可寫成式（4）形式，但實質(zhì)上已轉(zhuǎn)化為式（1）形式求解.

2 潮流計算的自動化、快速求解

2.1 阻抗矩陣的自動化生成

回路阻抗法潮流計算自動化求解的核心是阻抗矩陣自動化生成，實現(xiàn)的前提是找到一個靈活的數(shù)據(jù)結(jié)構(gòu)，能對配電網(wǎng)的物理拓撲進行存儲，易于配電網(wǎng)物理狀況的判斷和結(jié)構(gòu)調(diào)整.

2.1.1 十字鏈表的物理拓撲結(jié)構(gòu)存儲

配電網(wǎng)在正常運行情況下呈輻射狀結(jié)構(gòu)；在故障處理倒負荷操作或網(wǎng)絡重構(gòu)情況下呈現(xiàn)短時間的環(huán)網(wǎng)結(jié)構(gòu).這兩類拓撲結(jié)構(gòu)都可以用有向圖表示.有向圖的存儲方法主要有鄰接矩陣法、鄰接表法、十字鏈表法.鄰接矩陣法占用空間大；鄰接表法具有修改簡單、存儲方便的特點，調(diào)整較少的支路和節(jié)點信息即可實現(xiàn)配電網(wǎng)網(wǎng)絡重構(gòu)，但難以判斷是否有環(huán)網(wǎng)生成；十字鏈表是鄰接表法和逆鄰接表法結(jié)合起來得到的一種鏈表，不僅繼承了鄰接表法的優(yōu)點，而且能根據(jù)逆鄰接表更易判斷出構(gòu)成環(huán)網(wǎng)結(jié)構(gòu)的負荷節(jié)點，便于破環(huán)后潮流計算的求解.為此采用十字鏈表對復雜配電網(wǎng)物理拓撲結(jié)構(gòu)進行存儲.

考慮到配電系統(tǒng)的結(jié)構(gòu)特點和潮流計算的需求，對十字鏈表相應改進：配電網(wǎng)中的每個節(jié)點、每條支路在十字鏈表中都有一個頂節(jié)點、邊節(jié)點類型的存儲節(jié)點相對應.頂節(jié)點由四部分組成：“verlink”為指針，指向下一頂為數(shù)據(jù)域點；info為數(shù)據(jù)域，用于對vi、負荷節(jié)點標識flag、負荷節(jié)點阻抗值Zi、閉環(huán)負荷節(jié)點破環(huán)后增加的vr等數(shù)據(jù)信息的存儲；“firstin”為指針，指向第一個以該點為終點的支路；firstout為指針，指向第一個以該點為出發(fā)點的支路.邊節(jié)點由五部分組成：bpoint指針和epoint指針指示支路的出發(fā)點和終點，info存儲支路阻抗，elink為指針指向終點相同的下一條支路，blink為指針，指向出發(fā)點相同的下一條支路.

在構(gòu)建十字鏈表前，需在配電網(wǎng)原始接線圖上進行一次節(jié)點編號且編號方案任意（保證電源節(jié)點編號為零，其他節(jié)點無遺漏即可），此后再也無須額外的節(jié)點優(yōu)化編號.在根據(jù)節(jié)點數(shù)目構(gòu)建十字鏈表后，還需按照節(jié)點編號由小到大的順序依次填充所設置的十字鏈表.對于頂節(jié)點，若為負荷節(jié)點，置負荷節(jié)點標識flag＝vi（負荷節(jié)點編號）、vr＝0，同時記錄負荷節(jié)點阻抗值Zi；若為內(nèi)節(jié)點，上述參數(shù)統(tǒng)一設置為零.對于各邊節(jié)點，記錄以l為起點、k為終點的支路阻抗Zl，k，完成十字鏈表的初始化.

當配電網(wǎng)絡物理結(jié)構(gòu)改變重構(gòu)網(wǎng)絡時，只需在記錄其物理拓撲結(jié)構(gòu)的十字鏈表中添加或刪除頂節(jié)點、邊節(jié)點，并改變相關指針指向即可.現(xiàn)以7節(jié)點系統(tǒng)為例說明在十字鏈表上的網(wǎng)絡重構(gòu)過程.圖3a為配電網(wǎng)絡原始接線方式，對所有節(jié)點編號，編號順序任意；圖3b為其對應的十字鏈表.當接入2～5支路構(gòu)成環(huán)網(wǎng)時，需新建一邊節(jié)點（2，5）；修改v5頂點的firstin指針指向（2，5）邊節(jié)點；v2頂節(jié)點所連接的（2，3）邊節(jié)點的blink指針指向（2，5）邊節(jié)點，同時，（2，5）邊節(jié)點的elink指針指向（4，5）邊節(jié)點.

2.1.2 關聯(lián)矩陣

生成節(jié)點阻抗矩陣的關鍵是構(gòu)建關聯(lián)矩陣，使方程組中節(jié)點阻抗矩陣中的元素Zi，j與各輸入支路阻抗、節(jié)點負荷阻抗建立起正確的對應組合關系.

從圖1看出，若有n個節(jié)點，除去根節(jié)點（表示配電網(wǎng)供電電源）后，則有n－1個節(jié)點、n－1條支路，其支路編號可記為支路電流流入支路末端節(jié)點的節(jié)點編號.

根據(jù)配電網(wǎng)各節(jié)點之間的電流流向關系，參照文獻［12］對關聯(lián)矩陣重新定義使之更易理解和應用：關聯(lián)矩陣A為（n－1）×（n－1）階矩陣，i，j對應于除去根節(jié)點以外的節(jié)點編號，當存在流入節(jié)點j的電流且流經(jīng)節(jié)點（支路）i時，Ai，j＝1，否則 Ai，j＝0.表1即為圖3a所對應的關聯(lián)矩陣.

2.1.3 阻抗矩陣的建立與修改

假定每一負荷節(jié)點處為一個兩端網(wǎng)絡，其阻抗為R＋jΩ，則

其中：有功功率P＝｜U｜｜I｜cosφ；無功功率Q＝｜U｜｜I｜sinφ；｜U｜，｜I｜，φ分別是網(wǎng)絡兩端口處U，I的模和兩者的相角差.若已知vi的P，Q，U，則

同時，假定U＝1.0＋0j，則vi的初始阻抗Zi＝R＋j Ω＝cos2φ／P＋j sin2φ／Q.

由十字鏈表、關聯(lián)矩陣可構(gòu)建包含除供電電源之外的所有配電網(wǎng)負荷節(jié)點的初始阻抗矩陣Z，用于第一次算法求解.Z為一復系數(shù)對稱矩陣，由對角線元素Zi，i和非對角線元素Zi，j構(gòu)成.

當j是i的虛擬節(jié)點編號vr時，

Zi，vr可由式（9）計算.Ak，i是關聯(lián)矩陣的k 行i列元素，Zl，k是十字鏈表中起點為l、終點為k的支路阻抗.

以上是結(jié)合十字鏈表與關聯(lián)矩陣提出的一種阻抗矩陣生成算法.無論配電網(wǎng)是否由環(huán)網(wǎng)構(gòu)成，該算法都能完成阻抗矩陣元素的自動化計算.該算法主要包括判斷是否由環(huán)網(wǎng)構(gòu)成和阻抗矩陣元素計算兩個步驟：① 判斷是否由環(huán)網(wǎng)構(gòu)成：依次遍歷十字鏈表中頂節(jié)點，根據(jù)負荷節(jié)點標識flag≠0訪問其中的負荷節(jié)點.若每一負荷節(jié)點的firstin指針所指支路的elink指針均為null，則可判定該配電網(wǎng)無環(huán)網(wǎng)構(gòu)成；若某一負荷節(jié)點vi的firstin指針所指支路的elink指針不為空，則可判定在該負荷節(jié)點處形成閉環(huán).則將其info數(shù)據(jù)域中的閉環(huán)負荷節(jié)點破環(huán)后增加的虛擬節(jié)點編號vr＝h＋＋（h＝N－1，N為配電網(wǎng)節(jié)點總數(shù)），并在關聯(lián)矩陣中添加關于vr的行和列，并調(diào)整關聯(lián)矩陣中的元素值.② 阻抗矩陣元素計算：依次遍歷十字鏈表中的負荷節(jié)點，若無環(huán)網(wǎng)構(gòu)成，根據(jù)式（7），（9）依次計算阻抗矩陣中關于每一負荷節(jié)點的回路方程系數(shù)，直至所有負荷節(jié)點被訪問完畢即可生成阻抗矩陣；若在vi處有環(huán)網(wǎng)生成，則根據(jù)式（7），（10），（9），（8），（10），（9）順序，先后計算關于vi，vr的回路方程系數(shù)，其他節(jié)點的回路方程系數(shù)仍由式（7），（9）完成.隨著負荷節(jié)點遍歷結(jié)束，初始阻抗矩陣也就生成建立.

2.2 基于GPU的回路阻抗方程組求解

長期以來，GPU一直局限于處理圖形渲染的計算任務，通過在硬件上增加并行處理單元和存儲控制單元的方式，具備了同代CPU無可比擬的強大并行運算能力和高內(nèi)存讀寫帶寬.尤其是2007年，NVIDA公司的統(tǒng)一計算架構(gòu)（compute unified device architecture，CUDA）的推出，使GPU具有更好的可編程性，可廣泛用于圖形渲染以外領域的計算.與傳統(tǒng)的并行計算技術相比，基于GPU的并行計算具有低成本、高性價比的顯著優(yōu)勢.

從式（1）看出，基于回路阻抗法的配電網(wǎng)的潮流計算可以歸結(jié)為求解一組復系數(shù)線性方程組，其矩陣表達式為B＝Z｜VS（B 是n×（n＋1）矩陣，n是負荷節(jié)點的數(shù)目），是一種高密度數(shù)據(jù)運算，適用于GPU快速求解.

2.2.1 CUDA編程模型

在CUDA編程架構(gòu)下（如圖4所示），一個應用程序分為主機（host）端和設備（device）端兩部分：主機端是指在CPU上可執(zhí)行的部分，設備端是在顯卡上執(zhí)行的部分，其中，運行在GPU上的并行計算函數(shù)稱為內(nèi)核（kernel）函數(shù).在程序運行過程中，主機端程序準備好數(shù)據(jù)并復制到顯存中；設備端程序執(zhí)行主機端設置的一個個內(nèi)核函數(shù)，每個內(nèi)核函數(shù)都將按照線程網(wǎng)格的概念在GPU上運行（每個線程網(wǎng)格可包含多個線程塊，每個線程塊可包含多個線程，同一線程塊內(nèi)的線程可通過共享存儲器（shared memory）通信、同步）；爾后，主機端程序從顯卡內(nèi)存中取回運算結(jié)果［16］.

圖4 基于CUDA的編程模型結(jié)構(gòu)Fig.4 CUDA programming model

2.2.2 基于CUDA的復系數(shù)線性方程組求解

求解線性方程組的方法有兩大類：直接法和迭代法.迭代法速度快，但有誤差；直接法準確、可靠，具體有LU分解法和高斯消元法.從本質(zhì)上講，LU分解法是高斯消元法的一種變相表達形式，實質(zhì)上是將矩陣通過初等行變換變?yōu)?個上三角矩陣，變換矩陣是其單位下三角矩陣的過程.LU分解算法多，根據(jù)數(shù)據(jù)的訪問模式不同，可劃分為left－looking和right－looking兩類算法.left－looking算法是：每次迭代運算先計算矩陣中的1列各元素值，后對該列元素調(diào)整.求解過程不僅與當前列的左邊矩陣中的部分元素相關，而且各元素之間也存在計算依存關系，難以并行處理.right－looking算法是：每次迭代運算先計算矩陣中的1列，再用該列中的元素去計算該列右邊的子矩陣的各元素值.其求解計算不相關，可并行處理［17］.

為此，筆者根據(jù)right－looking算法可并行處理的特性，結(jié)合CUDA編程模型以及復系數(shù)矩陣的結(jié)構(gòu)特點，提出了基于right－looking LU分解法的并行高斯消去法（如圖5所示），將大量計算集中于GPU處理，充分發(fā)揮其并行運算的能力，將回代過程減為一次，在保證求解精度不變的情況下使求解速度大為提高.

圖5 基于right－looking LU分解法的并行高斯消去算法Fig.5 Parallel Gaussvan elimination algorithm based on rightlooking LU factorization

在LU分解法求解過程中，選取絕對值大的元素作為主元是確保矩陣算法穩(wěn)定求解應用最廣泛的技術，其中列主元法是其最常用的方法［18］.這一處理過程在一般LU分解算法中是一個必不可少的環(huán)節(jié).但在配電網(wǎng)潮流計算的直接求解法所形成的阻抗矩陣中，Zi，i＞Zj，i，在完成第（i－1）步迭代運算過后，第i列元素中的在該列的第i行之后的所有元素中絕對值始終最大.故這一矩陣置換選取列主元的步驟在本配電網(wǎng)潮流計算算法中可省略，節(jié)省了大量的運算時間.

在算法流程設計中，GPU上主要進行主列元上各元素和子矩陣中各元素的并行計算，可按照CUDA編程規(guī)則分別編寫對應的核函數(shù)column＿element（）和submatrix＿element（）以實現(xiàn)，每個核函數(shù)由多個線程并行完成.在主列元的各元素求解計算時，安排每個線程計算1個元素值；而在子矩陣各元素的求解計算時，讓每個線程來計算子矩陣中1列元素值.這是因為如果用每個線程來計算1個元素值，則所需的線程數(shù)將超越1個塊所允許使用的最大線程數(shù).而在CUDA中有一硬件特性：同1塊中的線程均可在同1個流多處理器中同步執(zhí)行相同的指令序列，且可共用共享存儲器，明顯提高運算速度.

為表述算法方便，在算法流程圖5中，用Bi，j來表示復系數(shù)矩陣中的某一復數(shù)元素，但在實際編程過程中，用 Ri，j，Ii，j對應其實部和虛部（即用與矩陣B同等階數(shù)的矩陣R和I來描述復數(shù)矩陣B）.因為每個塊的共享存儲器大小為16～64kbit，通常情況下難以將R，I整體讀入共享存儲器，但在每次循環(huán)計算時，可將主元列和主元列上第1個元素的同行元素從全局存儲器上讀入共享存儲器，而R和I存于全局存儲器，以此充分發(fā)揮共享存儲器的訪問速度遠比全局存儲器快得多的優(yōu)勢，從而減小通信延遲、提高數(shù)據(jù)訪問效率.

隨著2個內(nèi)核函數(shù)分別調(diào)用N－1次，整個LU分解過程結(jié)束，但還需由CPU調(diào)用函數(shù)back＿substitution（ ）回代計算，才能完成對復系數(shù)線性方程組的求解.

2.3 支路電流、內(nèi)節(jié)點電壓計算

經(jīng)第一次GPU求解可得到各負荷節(jié)點電流，如果負荷節(jié)點功率不滿足給定的精度要求，則需按式（6）對自阻抗Zi，i中的Zi調(diào)整.隨之，式（1）和（4）也需對Zi，i重新賦值，然后迭代求解下一次的GPU，直至滿足精度要求.然后，利用所得負荷節(jié)點電流計算配電網(wǎng)各支路電流、內(nèi)節(jié)點電壓，分別以Ii′，Uj′標識，i′，j′為內(nèi)節(jié)點編號.根據(jù)基爾霍夫定律和所定義的關聯(lián)矩陣與十字鏈表，可分別推導出它們的計算式

當j′為負荷節(jié)點時，可由式（1）迭代計算得到Ij′；Ii′，Uj′中的各值可按任意順序求解.

3 仿真信息可視化

信息可視化是指利用計算機技術，對抽象數(shù)據(jù)進行交互與可視化，以提高人的認知功能［19］.其相應的可視化技術有二維表格、三維圖、幾何變換技術、層次技術等幾大類［20］，但缺乏時空分析能力.因此，將GIS和虛擬現(xiàn)實技術相結(jié)合，用于進行潮流計算仿真信息可視化，使仿真抽象結(jié)果在三維虛擬空間相關實體設備上得以具體、直觀地體現(xiàn)；從三維角度全面反映仿真配電網(wǎng)的運行狀況，從而有助于調(diào)度人員的分析決策和培訓人員的學習理解.

3.1 基于GIS技術的配電網(wǎng)三維仿真系統(tǒng)

實現(xiàn)三維GIS系統(tǒng)開發(fā)的軟件工具有IMAGINE VirtualGIS， ArcGis 3DAnalyst，GeoMedia Terrain，PAMAP GIS Topographer，Java與Vrml，Creator與Vega Prime，前四個主要用于三維地理環(huán)境的分析與瀏覽，后兩類軟件組合可用于基于GIS系統(tǒng)的視景仿真.在最后一類軟件組合中，Creator軟件具有強大的多邊形三維建模能力，大面積地形精確生成功能，擁有針對實時應用優(yōu)化的openflight模型數(shù)據(jù)格式；Vega Prime支持Creator模型文件，提供眾多附加仿真功能模塊，迅速創(chuàng)建、渲染各種實時交互的三維仿真環(huán)境.為此，選取上述兩款軟件進行3DGIS動態(tài)視景仿真開發(fā).

實現(xiàn)過程見圖6.首先，使用Photoshop對谷歌地球衛(wèi)星圖片裁剪，然后在creator terrain pro模塊中完成地形建模并以＊.mft格式存儲；建筑物建模過程較為復雜：先要搜集各建筑物圖紙，在autocad軟件中提取出需要的建筑物現(xiàn)狀圖，后以.dxf文件格式導入creator軟件中，用face工具進行地物匹配、geometry工具設置建筑物高度，并用裁剪后的實物照片映射紋理；對于配電網(wǎng)電力設備、花草、樹木等其他要素建模，可根據(jù)所采集的實物照片在Creator軟件中直接進行；這兩個所建模型均以＊.flt格式存儲.此后，根據(jù)實際建筑物和其他要素在三維地理空間中的位置關系，在Vega Prime軟件lynx prime場景編輯器中，將上述＊.flt和＊.mft模型合成，設置生成＊.acf文件格式存儲.最后，在.net vc環(huán)境中調(diào)用Vega Prime API接口函數(shù)讀?。?acf文件，并根據(jù)潮流計算結(jié)果編程設置，實現(xiàn)視景驅(qū)動和仿真信息可視化.

圖6 3DGIS潮流計算仿真信息可視化Fig.6 3DGIS power flow calculation simulation information visualization process

3.2 潮流計算結(jié)果可視化顯示

在實際電力生產(chǎn)作業(yè)中，運行人員通過讀取配電網(wǎng)各變電站一次設備的儀表指示器來檢測配電網(wǎng)的運行情況.如今，三維虛擬現(xiàn)實技術提供了這樣的技術條件，可使使用人員在能夠反映真實場景的虛擬環(huán)境中模擬觀測，使實習者與管理人員將潮流計算中各節(jié)點的計算結(jié)果與監(jiān)測其變化的實體設備儀表一一對應起來，從而深刻理解電網(wǎng)運行原理，提高業(yè)務能力.其中，儀表指針讀數(shù)和控制該儀表的信號數(shù)據(jù)之間的參數(shù)轉(zhuǎn)換公式為

配電網(wǎng)虛擬設備儀表的指針讀數(shù)可視化顯示通過自由度（degrees of freedom，DOF）計算機動畫技術實現(xiàn).該技術主要包括DOF節(jié)點設置和DOF節(jié)點驅(qū)動兩個步驟.DOF節(jié)點設置在Creator中完成，DOF節(jié)點動畫特效在Vega Prime中觸發(fā).在該項技術中，設置成DOF的節(jié)點可控制它的所有子節(jié)點在設置的自由度范圍內(nèi)移動或者旋轉(zhuǎn)，運動范圍由局部坐標系來指定.真實場景中當設備儀表讀數(shù)變化時，其物理表現(xiàn)形式是儀表指針繞著指針軸旋轉(zhuǎn).因此，在Creator中首先建立一個以指針軸與指針鏈接處為局部坐標系原點中心的DOF節(jié)點，然后儀表指針作為它的子節(jié)點，同時設置滿刻度量程旋轉(zhuǎn)的角度.在潮流計算結(jié)束、用戶進入3DGIS配電網(wǎng)虛擬環(huán)境時，由在Vega Prime視景驅(qū)動軟件中的編程語句rotate（儀表指針讀數(shù)）觸發(fā)旋轉(zhuǎn)效果.

從潮流計算的結(jié)果來看，配電網(wǎng)能量分布細化到了最底層的變壓器一級.而現(xiàn)實世界中，此類設備并不配置測量儀表，所以，在3DGIS虛擬環(huán)境中就無法直觀表征其數(shù)值變化.此外，為方便用戶精確了解設備運行狀態(tài)，此類虛擬設備旁邊均設文字說明.

4 實驗結(jié)果與討論

本算法在一臺Intel Core2Quad Q95502.8 GHZ CPU，3GB內(nèi)存，Nvidia GeForce GTX480顯卡的PC機上實現(xiàn)了配電網(wǎng)多種節(jié)點類型設置的潮流計算與仿真信息可視化，使用Visual 2005作為軟件開發(fā)平臺，CUDA SDK 3.2作為GPU程序的主要開發(fā)工具，Vega Prime作為3DGIS仿真可視化工具.

4.1 潮流計算算例

基于文獻［9］和本文所提方法，通過編程，分別實現(xiàn)了基于CPU和基于GPU并行計算的配電網(wǎng)潮流計算程序，并采用文獻［9，11，21］提供的6，20（包括環(huán)網(wǎng)），45，33負荷節(jié)點數(shù)的配網(wǎng)系統(tǒng)原始數(shù)據(jù)以驗證本算法的高效性.在達到10－5級的計算精度條件下，兩類求解算法單次收斂所需的迭代次數(shù)和100次運算花費時間對比如表2所示.

表2 CPU串行與GPU并行求解算法比較Tab.2 Comparison between CPU serial algorithm and GPU parallel algorithm

表2數(shù)據(jù)顯示，本方法對無環(huán)、弱環(huán)配電網(wǎng)潮流計算都能快捷、精確求解.尤其是采用GPU后，在達到CPU串行求解同等精度條件下，加速比隨著節(jié)點數(shù)目的增大而增大：45負荷節(jié)點的潮流計算運算速度已提高到2倍，當負荷節(jié)點數(shù)目大于45時，還可獲得2倍以上的加速比；但當矩陣規(guī)模小于20時，GPU方法實現(xiàn)時間大于CPU.這是由于計算量小、受顯卡數(shù)據(jù)總線帶寬所限.顯存與內(nèi)存數(shù)據(jù)通信所花時間較長，在整個程序運行完成總時間中占據(jù)了較大比重，無法發(fā)揮GPU高效并行計算的能力.

4.2 仿真信息可視化顯示

以上海洋山港110kV供電系統(tǒng)中的23kV配電網(wǎng)作為仿真系統(tǒng)物理設備模型，線路阻抗和節(jié)點負荷功率取自文獻［9］的表5a，b.

表3即是本方法與Goswami所提方法在文獻［9］中所得節(jié)點電壓的求解結(jié)果比較.兩者的計算數(shù)值相近，說明本方法的正確性；在同樣迭代求解3次的情況下，本方法收斂精度誤差是0.9758×10－5p.u（標幺值），遠遠小于文獻［9］方法的3.6716×10－5p.u.說明本方法的計算精度較高.圖7a所示為面積為8.14km2洋山港配電系統(tǒng)3DGIS仿真整體繪制效果圖；圖7b為配電系統(tǒng)的地下電纜三維仿真繪制效果圖；圖7c為23kV配電網(wǎng)中監(jiān)測第17母線的電力開關柜，電流儀表指針指示為34A，并附有該設備所處位置的經(jīng)度、緯度地理信息和潮流計算所得的該母線的電流幅值和角度.圖7d為23kV配電網(wǎng)中第7負荷點對應的變壓器.此設備沒有數(shù)值顯示儀，本系統(tǒng)在設備旁以文字說明其位置地理信息和潮流計算所得的電壓幅值、角度等仿真計算信息.4個圖均為固定視點抓圖，但在程序運行中可通過自主漫游方式或設定路線方式對潮流計算仿真信息進行可視化展示.本方法所繪制圖像具有高度的真實感，人機交互方式靈活，能正確、生動、形象地反映仿真配電網(wǎng)設備三維空間位置關系及其運行狀態(tài).

表3 潮流計算結(jié)果Tab.3 Results of the power flow calcalation

5 結(jié)語

在回路阻抗法的基礎上，針對配電系統(tǒng)，提出了一種快速潮流計算與仿真信息可視化方法.其創(chuàng)新性主要在于：十字鏈表用于配電網(wǎng)的網(wǎng)絡拓撲結(jié)構(gòu)存儲，便于網(wǎng)絡重構(gòu)與構(gòu)成弱環(huán)網(wǎng)的負荷節(jié)點的判定，與關聯(lián)矩陣相結(jié)合，可完成對無環(huán)或有環(huán)配電網(wǎng)阻抗矩陣的自動生成，并可不受節(jié)點編號順序限制計算各內(nèi)節(jié)點電壓、支路電流；基于CUDA編程模型提出了基于right－looking LU分解法的并行高斯消去算法，利用GPU對潮流計算快速求解；并將計算結(jié)果在3DGIS系統(tǒng)中進行仿真可視化.算例仿真實驗證明，該方法對節(jié)點編號無特殊要求，求解計算快速、精確、穩(wěn)定、交互性能好，能使調(diào)度人員、培訓人員對各種節(jié)點類型參數(shù)設置的無環(huán)或有環(huán)配電網(wǎng)進行潮流計算和三維可視化仿真，從中直觀、形象地認識到地理和設備之間的對應關聯(lián)信息，理順電力設備之間的運行邏輯關系，分析掌握配電網(wǎng)的現(xiàn)行運行狀態(tài)，切實提高管理水平和學習效率，可廣泛應用于配電網(wǎng)規(guī)劃、分析、人員培訓等領域.

今后，將基于該方法對配電網(wǎng)短路、三相潮流等特殊性質(zhì)的潮流計算進行仿真模擬與可視化，同時結(jié)合高壓輸電網(wǎng)的潮流計算方法，求解大規(guī)模電力系統(tǒng)的潮流計算，充分發(fā)揮該算法并行計算的高效性.

［1］Godwin N.Application of modern load flow techniques to electric power systems［M］.Saarbrucken：LAP LAMBERT Academic Publishing，2010.

［2］Teng J.A modified Gauss－Seidel algorithm of three－phase power flow analysis in distribution networks［J］.International Journal of Electrical Power and Energy Systems，2002，24（2）：97.

［3］Jegatheesan R，Nor N M，Romlie M F.Newton－Raphson power flow solution employing systematically constructed Jacobian matrix［C］／／2008 IEEE2nd International Power and Energy Conference.Piscataway：IEEE，2008：180－185.

［4］Kulworawanichpong T.Simplified Newton Raphson power－flow solution method［J］.International Journal of Electrical Power and Energy Systems，2010，32（6）：551.

［5］Yang H，Wen F，Wang L，et al.Newton－Downhill algorithm for distribution power flow analysis［C］／／2008 IEEE 2nd International Power and Energy Conference.Piscataway：IEEE，2008：16L28－1632.

［6］Aravindhababu P，Ashokkumar R.A fast decoupled power flow for distribution systems［J］.Electric Power Components and Systems，2008，36（9）：932.

［7］Chang G W，Chu S Y，Wang H L.An improved backward／forward sweep load flow algorithm for radial distribution systems［J］.IEEE Transactions on Power Systems，2007，22（2）：882.

［8］Eminoglu U，Hocaoglu M H.Distribution systems forward／backward sweep－based power flow algorithms：a review and comparison study ［J］. Electric Power Components and Systems，2009，37（1）：91.

［9］Goswami S K，Basu S K.Direct solution of distribution systems［J］.IEE Proceedings C：Generation，Transmission and Distribution，1991，138（1）：78.

［10］王守相，阮同軍，劉玉田.配電網(wǎng)潮流計算的回路阻抗法［J］.電力系統(tǒng)及其自動化學報.1998，10（1）：12.WANG Shouxiang，RUAN Tongjun，LIU Yutian.Distribution system power flow based on loop－impedance equations［J］.Proceedings of the International Conference on Energy Management and Power Delivery，1998，10（1）：12.

［11］王丹，常寶立.一種用于配網(wǎng)潮流計算的節(jié)點編號新方法［J］.電力系統(tǒng)及其自動化學報.2003，15（1）：22.WANG Dan，CHANG Baoli.Distribution system power flow based on loop－impedance equations［J］.Proceedings of the International Conference on Energy Management and Power Delivery，2003，15（1）：22.

［12］劉耀年，豈小梅，李國鵬，等.基于回路阻抗法的配電網(wǎng)潮流計算［J］.繼電器，2004，32（8）：8.LIU Yaonian，QI Xiaomei，LI Guopeng，et al.Load flow algorithm for distribution network based on loop－impedance method［J］.Relay，2004，32（8）：8.

［13］汪宇霆，張焰，張益波.基于改進回路電流法的配電網(wǎng)潮流通用算法［J］.電力系統(tǒng)保護與控制.2010，38（20）：57.WANG Yuting，ZHANG Yan，ZHANG Yibo.A general load flow method for distribution systems based on improved loop－current method［J］.Power System Protection and Control，2010，38（20）：57.

［14］Sun Y， Overbye T J. Visualizations for power system contingency analysis data［J］.IEEE Transactions on Power Systems，2004，19（4）：1859.

［15］Milano F.Three－dimensional visualization and animation for power systems analysis［J］.Electric Power Systems Research，2009，79（12）：1638.

［16］NVIDIA.NVIDIA CUDA programming guide 3.2［EB／OL］（2010－10－22）.http：／／developer.download.nvidia.com.

［17］Dongarra J J，Hammarling S，Walker D W.Key concepts for parallel out－of－core LU factorization ［J］. Computers and Mathematics with Applications，1998，35（7）：13.

［18］Tomov S，Dongarra J，Baboulin M.Towards dense linear algebra for hybrid GPU accelerated manycore systems［J］.Parallel Computing，2010，36（5－6）：232.

［19］Wakita A，Matsumoto F.Information visualization with Web 3D spatial visualization of human activity area and its condition［J］.Computer Graphics（ACM），2003，37（3）：29.

［20］Ltifi H，Ayed M B，Alimi A M，et al.Survey of information visualization techniques for exploitation in KDD［C］／／7th IEEE／ACS International Conference on Computer Systems and Applications.Piscataway：IEEE，2009：218－225.

［21］Baran M E，Wu F F.Network reconfiguration in distribution systems for loss reduction and load balancing［J］.IEEE Transactions on Power Delivery，1989，4（2）：1401.