林金榕，張 逸，李傳棟，邵振國，李寬宏

（1. 福州大學(xué) 電氣工程與自動化學(xué)院，福建 福州 350108；2. 國網(wǎng)福建省電力有限公司電力科學(xué)研究院，福建 福州 350007；3. 國網(wǎng)福建省電力有限公司福州供電公司，福建 福州 350009）

0 引言

2020 年我國提出2030 年碳達(dá)峰和2060 年碳中和的“雙碳”目標(biāo)。在該背景下，清潔能源的消納以及分布式電源（DG）的高比例接入，都對配電網(wǎng)的規(guī)劃和運(yùn)行提出了更高的要求，由此導(dǎo)致的配電網(wǎng)諧波、三相不平衡等電能質(zhì)量問題的影響越來越大，這些電能質(zhì)量問題對配電網(wǎng)的安全運(yùn)行造成了很大的危害［1-4］。為了評估電能質(zhì)量問題的影響并為電能質(zhì)量治理提供依據(jù)，三相諧波潮流計算將是一種重要的分析手段。

現(xiàn)有的配電網(wǎng)潮流仿真分析工具［5-6］在處理環(huán)網(wǎng)、高滲透率DG 接入引起的多重電能質(zhì)量問題時建模繁瑣、效率較低，而且難以滿足在線實時應(yīng)用需求。文獻(xiàn)［7-9］基于前推回代法對三相諧波潮流進(jìn)行研究，但由于前推回代法固有的局限，在處理各類DG 和局部存在的環(huán)網(wǎng)時需要對該方法進(jìn)行改進(jìn)［10-12］，這不可避免地增加了算法復(fù)雜度，使計算次數(shù)增多，收斂性下降。文獻(xiàn)［13］基于PQ分解法對三相諧波潮流進(jìn)行研究，該方法簡化了牛頓-拉夫遜法的實現(xiàn)過程，但PQ分解法的應(yīng)用前提是電抗遠(yuǎn)大于電阻，而且PQ分解法在算法收斂過程中不修正搜索方向，收斂性差，這導(dǎo)致其在配電網(wǎng)中難以應(yīng)用。綜合現(xiàn)有的研究情況可知，配電網(wǎng)三相諧波潮流的研究方法主要有2 種：一是基于前推回代法的潮流算法，其實現(xiàn)簡單，但由于算法固有的梯度簡化問題，在配電網(wǎng)DG 高比例接入和電壓控制策略日益復(fù)雜的情況下，其優(yōu)勢日益削弱，缺陷日益凸現(xiàn)；二是基于牛頓-拉夫遜法的潮流算法，其優(yōu)勢在于搜索方向用的是完整的梯度，收斂性好，同時算法可直接擴(kuò)展到各類潮流控制策略，其缺點是三相諧波導(dǎo)納矩陣和雅可比矩陣構(gòu)造十分復(fù)雜，難以構(gòu)建規(guī)范的實現(xiàn)形式，模型擴(kuò)展和應(yīng)用難，因此相關(guān)研究成果極少。

為了適應(yīng)未來配電網(wǎng)發(fā)展中DG 高比例接入和潮流控制方式日益復(fù)雜的情況，本文提出一種通用的規(guī)范型式下的三相諧波潮流算法。首先，采用矩陣索引技術(shù)生成牛頓-拉夫遜法相關(guān)的矩陣；其次，利用矩陣范型生成三相諧波導(dǎo)納矩陣；然后，利用生成的矩陣計算三相諧波潮流；最后，分別利用IEEE 14 和IEEE 123 節(jié)點系統(tǒng)對本文算法的準(zhǔn)確性、有效性、效率以及處理環(huán)網(wǎng)的能力進(jìn)行驗證分析。

1 三相諧波潮流基礎(chǔ)模型

參考文獻(xiàn)［14-15］中電氣元件對應(yīng)的諧波等值阻抗，限于篇幅，本文只介紹后文仿真中所使用的模型。

1）線路模型。

當(dāng)計算的諧波次數(shù)不高并且線路長度不大于300 km時，線路單相諧波阻抗近似為：

式中：h為諧波次數(shù)；r、x、z、b分別為電阻、電抗、阻抗、電納，下標(biāo)L表示線路，下標(biāo)h表示第h次諧波。

2）變壓器模型。

式中：下標(biāo)G表示發(fā)電機(jī)。

2 矩陣范型

2.1 矩陣范型算法原理

由于三相諧波潮流模型復(fù)雜，相關(guān)的工程應(yīng)用一直難以推進(jìn)。本文基于矩陣范型框架，提出一種將三相諧波潮流計算中的矩陣元素生成、修改和索引等過程轉(zhuǎn)化為矩陣運(yùn)算的算法實現(xiàn)框架。該框架組合矩陣內(nèi)積、關(guān)聯(lián)矩陣和矩陣索引等技術(shù)，并通過標(biāo)準(zhǔn)的矩陣運(yùn)算過程實現(xiàn)電網(wǎng)潮流算法和計算機(jī)底層數(shù)據(jù)結(jié)構(gòu)的解耦。本質(zhì)上，矩陣范型是一種散列映射算法，其實現(xiàn)了從單相電網(wǎng)的設(shè)備序參數(shù)列表到三相電網(wǎng)導(dǎo)納矩陣元素的規(guī)范映射。矩陣范型示意圖如圖1所示。圖中：a1—a3為向量的元素；l為矩陣非零元素行坐標(biāo)；q為矩陣非零元素列坐標(biāo)；g為矩陣非零元素數(shù)值；f為元件首端節(jié)點；t為元件末端節(jié)點；下標(biāo)f f表示首端節(jié)點對首端節(jié)點；下標(biāo)tt表示末端節(jié)點對末端節(jié)點；下標(biāo)f t表示首端節(jié)點對末端節(jié)點；下標(biāo)t f表示末端節(jié)點對首端節(jié)點；?表示克羅內(nèi)克積。此外，由于工程化的需要，本文三相諧波潮流算法中考慮三序參數(shù)和三相參數(shù)的相互轉(zhuǎn)換、線路和變壓器參數(shù)到三相諧波導(dǎo)納矩陣的轉(zhuǎn)換以及矩陣的重排和隨機(jī)訪問。若采用傳統(tǒng)的面向元素級別的稀疏矩陣方法來實現(xiàn)以上過程，則編譯和調(diào)試難度大，且難以實現(xiàn)模型的便捷擴(kuò)展和算法的功能增加。

圖1 矩陣范型示意圖Fig.1 Schematic diagram of matrix paradigm

2.2 矩陣范型過程的實現(xiàn)

為降低配電網(wǎng)三相建模的復(fù)雜性，減少計算量，本文將序網(wǎng)參數(shù)矩陣作為輸入，并通過相應(yīng)的矩陣范型變換得到本文所需的三相諧波導(dǎo)納矩陣。范型變換主要分為如下3步。

1）為了便于數(shù)據(jù)的輸入，通過轉(zhuǎn)換矩陣Tr整合同一元件同次諧波的三序參數(shù)，并進(jìn)行相序轉(zhuǎn)換操作，得到三相諧波基礎(chǔ)導(dǎo)納矩陣YHbase，其表示元件首端到末端節(jié)點之間的三相支路導(dǎo)納矩陣，本文稱為三相諧波基礎(chǔ)導(dǎo)納矩陣。具體過程如下：

式中：下標(biāo)H 表示矩陣包含各次諧波分量；下標(biāo)base表示基礎(chǔ)分量參數(shù)；下標(biāo)120 表示三序；Ip為p維單位矩陣，Ep為元素全為1的p維列向量，p=n(m-1)，n為元件總數(shù)，m為諧波的最大次數(shù)；1∶1∶3p表示以1 為起點，1 為步長，3p為終點的列向量，其他同理；Tabc120、T120abc為相序轉(zhuǎn)換矩陣；函數(shù)gsparse(L，M，N)表示以向量L的元素為行坐標(biāo)、向量M的元素為列坐標(biāo)、向量N的元素為數(shù)值，生成稀疏矩陣。

2）在三相諧波基礎(chǔ)導(dǎo)納矩陣的基礎(chǔ)上繼續(xù)推導(dǎo)出三相諧波子導(dǎo)納矩陣，包括元件首、末端節(jié)點的自導(dǎo)納矩陣和互導(dǎo)納矩陣。

式中：K、U分別為變壓器原邊、副邊分接頭變比矩陣。

3）通過構(gòu)造索引矩陣的方式，將三相諧波子導(dǎo)納矩陣整合為三相諧波導(dǎo)納矩陣YH。具體過程如下：

圖2 給出了利用本文范型生成三相導(dǎo)納矩陣的過程。圖中：fi、ti(i=1，2，…，n) 分別為第i個元件首、末端節(jié)點編號。

圖2 三相諧波導(dǎo)納矩陣生成圖Fig.2 Generation diagram of three-phase harmonic admittance matrix

3 三相諧波潮流過程的實現(xiàn)

3.1 三相諧波導(dǎo)納矩陣

在具備上述矩陣范型的基礎(chǔ)上，生成各種模型的三相諧波導(dǎo)納矩陣。

1）三相線路諧波導(dǎo)納矩陣。

首先，構(gòu)造范型式（8）中的輸入矩陣：

令范型中的n取線路數(shù)，由范型式（8）得到三相線路諧波基礎(chǔ)導(dǎo)納矩陣YLHbase，然后由YLHbase推導(dǎo)出三相線路諧波子導(dǎo)納矩陣YLHsub，再根據(jù)范型式（13）即可獲得三相線路諧波導(dǎo)納矩陣YLH。

2）三相變壓器諧波導(dǎo)納矩陣。

由于本文所采用的三相變壓器模型無需進(jìn)行相序轉(zhuǎn)換，因此范型式（8）中的相序轉(zhuǎn)換矩陣可省去。范型式（8）中的輸入矩陣為：

令范型中的n取變壓器數(shù)，由范型式（8）得到三相變壓器諧波基礎(chǔ)導(dǎo)納矩陣YTHbase，然后由YTHbase推導(dǎo)出三相變壓器諧波子導(dǎo)納矩陣YTHsub，再根據(jù)范型式（13）即可獲得三相變壓器諧波導(dǎo)納矩陣YTH。

3）三相負(fù)荷／發(fā)電機(jī)諧波導(dǎo)納矩陣。

由于本文不考慮負(fù)荷與發(fā)電機(jī)之間的相互影響，其對應(yīng)的三相諧波導(dǎo)納矩陣是對角矩陣，可以直接得到三相負(fù)荷諧波導(dǎo)納矩陣YSH和三相發(fā)電機(jī)諧波導(dǎo)納矩陣YGH，無需利用范型。

3.2 牛頓-拉夫遜三相基波潮流

牛頓-拉夫遜潮流計算的基礎(chǔ)方程形式為：

式中：U為電壓，上標(biāo)set 表示設(shè)定值，下標(biāo)x表示實部，下標(biāo)y表示虛部；P為有功功率；Q為無功功率；Δ表示差分。

3.3 諧波潮流算法

本文諧波源采用的是諧波電流含有率模型，諧波源相角、幅值與基波相角、幅值的關(guān)系參考文獻(xiàn)［7］。諧波部分的算法使用交替迭代法，可將基波和諧波潮流分開交替求解，即求解基波潮流時以上一次迭代的諧波潮流結(jié)果為邊界條件，求解諧波潮流時以上一次迭代的基波潮流結(jié)果為邊界條件，并基于功率守恒原理實現(xiàn)基波和諧波潮流的交替迭代。具體步驟如圖3所示。

圖3 三相諧波潮流算法流程圖Fig.3 Flowchart of three-phase harmonic power flow algorithm

式中：IH為諧波三相電流列向量。

4 算例分析

4.1 有效性驗證

本節(jié)以IEEE 14 節(jié)點系統(tǒng)為例，在節(jié)點7、14 接入如表1所示的非線性負(fù)荷，即諧波源H1和H2，接線如附錄A 圖A1 所示，對應(yīng)的電流頻譜如表2 所示。仿真軟件為MATLAB 9.5，硬件為CPU inter（R）Core（TM）i7-8750H、8 GB內(nèi)存。

表1 非線性負(fù)荷節(jié)點信息Table 1 Bus information of nonlinear loads

表2 2種諧波源的電流頻譜Table 2 Current spectra of two harmonic sources

為了驗證本文三相諧波潮流計算結(jié)果的準(zhǔn)確性和有效性，分別用本文算法和OpenDSS 進(jìn)行對比計算。潮流收斂精度取為10-10，將本文矩陣范型中的線路元件數(shù)n設(shè)置為13，最大諧波次數(shù)m設(shè)置為13，得到線路的三相諧波導(dǎo)納矩陣，再結(jié)合諧波源電流頻譜即可計算三相諧波潮流。限于篇幅，算法范型實現(xiàn)的中間過程如附錄B 所示，各節(jié)點三相電壓的總諧波畸變率（THD）如圖4 所示，以節(jié)點13 為例，2種算法的三相諧波電壓如圖5所示。

圖4 本文算法與OpenDSS三相電壓的THD計算結(jié)果對比Fig.4 Comparison of THD calculation results for threephase voltage between proposed algorithm and OpenDSS

圖5 節(jié)點13諧波電壓結(jié)果對比Fig.5 Comparison of harmonic voltage for Bus 13

由圖4和圖5可以看出，本文算法計算得到的諧波電壓和電壓THD 的偏差（電壓THD 的平均相對偏差為2.83%，諧波電壓平均相對偏差為4.07%）均在合理范圍內(nèi)，滿足工程需求，驗證了本文算法的準(zhǔn)確性與有效性。算法計算結(jié)果的細(xì)微差異主要源于本文在元件模型、計算精度上與OpenDSS存在差異。

4.2 算法對比

本節(jié)在IEEE 123 節(jié)點系統(tǒng)中對本文算法與傳統(tǒng)前推回代法進(jìn)行比較，潮流收斂精度取為10-10。為了驗證本文算法的計算性能，設(shè)置3 種不同運(yùn)行方式進(jìn)行對比，不同方式所對應(yīng)的負(fù)荷功率情況如表3所示，2種算法結(jié)果對比如圖6所示。

表3 IEEE 123節(jié)點系統(tǒng)不同運(yùn)行方式的負(fù)荷功率情況Table 3 Load power situation of IEEE 123-bus system under different operation modes

由圖6 可知，隨著負(fù)荷增大，本文算法的迭代次數(shù)和運(yùn)行時間幾乎維持不變，而傳統(tǒng)前推回代法的迭代次數(shù)和運(yùn)行時間均有不同程度的增長。這是由于牛頓-拉夫遜法采用的是完整的梯度信息，而前推回代法采用的是近似的梯度，搜索效率相對低下。此外，在重負(fù)荷情況下，由于潮流解接近靜態(tài)安全極限，梯度的搜索效率必然會下降。由此可見，本文算法具有計算速度快、效率高和收斂性好的優(yōu)點。

圖6 本文算法和傳統(tǒng)前推回代法收斂情況對比Fig.6 Convergence situation comparison between proposed algorithm and traditional forward and backward algorithm

在上述IEEE 123 節(jié)點系統(tǒng)的基礎(chǔ)上，改變系統(tǒng)連接方式，在保持基礎(chǔ)連接不變的情況下，連接系統(tǒng)中線路11-27、17-96、36-51、54-94、39-66、119-120、114-122、75-85，構(gòu)成多個環(huán)網(wǎng)，如附錄A圖A2所示。由于傳統(tǒng)前推回代法不具備計算環(huán)網(wǎng)的能力，僅給出本文算法的收斂情況，如表4所示。

表4 具有環(huán)網(wǎng)情況下的收斂情況Table 4 Convergence situation in case of ring network

由表4 可知，在存在環(huán)網(wǎng)的情況下，本文算法的收斂時間與迭代次數(shù)與正常情況下相比幾乎不變，說明本文算法有較好的處理環(huán)網(wǎng)的能力。穩(wěn)定的迭代次數(shù)和較短的運(yùn)行時間得益于本文矩陣范型較高的運(yùn)行效率和牛頓-拉夫遜法自身的特性。

由以上算例可知，由于本文算法是基于牛頓-拉夫遜法和散列映射原理構(gòu)建三相諧波潮流的網(wǎng)絡(luò)模型，因此算法在環(huán)網(wǎng)和放射型網(wǎng)絡(luò)下的處理效率是相同的。此外，由于本文算法的潮流計算時間為ms級，因此，即使在配電網(wǎng)拓?fù)浣Y(jié)構(gòu)的頻繁切換過程中，本文算法的收斂性和速度仍可以得到保證，在工程實際中具有廣泛應(yīng)用的潛力。

5 結(jié)論

本文提出基于牛頓-拉夫遜法并按照矩陣范型運(yùn)算的三相諧波潮流新算法，該算法保留了牛頓-拉夫遜法二次收斂的穩(wěn)定性，適用于所有電網(wǎng)結(jié)構(gòu)和運(yùn)行方式，得到如下結(jié)論。

1）基于各序網(wǎng)絡(luò)獨立的原理以及矩陣范型的思想進(jìn)行算法設(shè)計，實現(xiàn)了三相諧波潮流仿真的快速開發(fā)和高效計算。將各序網(wǎng)絡(luò)用矩陣范型的方式進(jìn)行描述，并充分運(yùn)用多重索引矩陣等技術(shù)，極大地降低了相應(yīng)算法的復(fù)雜性。

2）基于面向矩陣運(yùn)算的算法，新設(shè)備模型的擴(kuò)展可以通過矩陣以及向量疊加的方式實現(xiàn)，該方式靈活方便，且可以充分運(yùn)用成熟的矩陣計算程序?qū)崿F(xiàn)程序的速度優(yōu)化。

3）算法兼容樹形網(wǎng)絡(luò)和環(huán)形網(wǎng)絡(luò)，能快速適應(yīng)目前日益復(fù)雜的配電網(wǎng)仿真。

筆者后續(xù)將嘗試采用并行計算方式進(jìn)一步提高本文算法的效率。

附錄見本刊網(wǎng)絡(luò)版（http：//www.epae.cn）。