潘世賢，李志剛，鄭杰輝，季天瑤，陳思思

(華南理工大學(xué) 電力學(xué)院，廣東 廣州 510641)

如何以經(jīng)濟(jì)、安全、可靠的方式調(diào)度現(xiàn)有的發(fā)電機(jī)組以滿足用戶的負(fù)荷需求，一直是電力系統(tǒng)最關(guān)注的問題之一。這類優(yōu)化問題被稱為經(jīng)濟(jì)調(diào)度問題，即在滿足負(fù)荷需求、機(jī)組出力和運(yùn)行限制等一系列等式和不等式約束條件下，確定電力系統(tǒng)不同發(fā)電機(jī)組的出力計劃，最小化機(jī)組的總?cè)剂铣杀?。根?jù)是否考慮爬坡約束等動態(tài)約束，可以分為靜態(tài)經(jīng)濟(jì)調(diào)度和動態(tài)經(jīng)濟(jì)調(diào)度[1]。用戶的負(fù)荷需求波動很大，而發(fā)電機(jī)組的出力受到機(jī)組自身及其他條件限制，難以完全實(shí)時跟隨負(fù)荷波動，需要在一段時期內(nèi)進(jìn)行優(yōu)化，因此動態(tài)經(jīng)濟(jì)調(diào)度更符合實(shí)際。提高經(jīng)濟(jì)調(diào)度問題求解算法的準(zhǔn)確性、高效性、魯棒性，可以為電力系統(tǒng)高效運(yùn)行創(chuàng)造巨大的經(jīng)濟(jì)效益。

在求解經(jīng)濟(jì)調(diào)度問題的算法中：數(shù)學(xué)規(guī)劃算法計算速度較快，有嚴(yán)格的理論基礎(chǔ)，但是其性能對初始猜測很敏感，容易收斂到局部最優(yōu)，出現(xiàn)發(fā)散時難以確定解的存在性[2]；啟發(fā)式算法(如遺傳算法[3]、差分演化[4]、模擬退火法[5]等)在一定程度上可以跳出局部最優(yōu)，但計算時間長，算法針對性比較強(qiáng)，具有大量特定參數(shù)，魯棒性差。

近年來，非迭代的全純函數(shù)嵌入式方法(holomorphic embedding method，HEM)憑借良好的性能，受到不少學(xué)者的關(guān)注，已經(jīng)被應(yīng)用到電力系統(tǒng)潮流分析[6]、交直流潮流分析[8]、電壓穩(wěn)定性分析[10]、系統(tǒng)薄弱點(diǎn)識別[11]、連續(xù)潮流計算[12]和電氣聯(lián)合網(wǎng)絡(luò)的動態(tài)仿真[13]等方面。HEM本質(zhì)是一種遞歸的等式方程組求解器，可以用于求解經(jīng)濟(jì)調(diào)度問題的Karush-Kuhn-Tucker(KKT)方程，尋找最優(yōu)的調(diào)度方案。在KKT方程嵌入復(fù)變量，構(gòu)造全純的嵌入式方程；代入展開的變量級數(shù)，推導(dǎo)出系數(shù)間的線性遞歸方程，并計算級數(shù)的各階系數(shù)；為了加快收斂速度，可以利用帕德近似擴(kuò)大冪級數(shù)的收斂半徑，使原問題盡快落在收斂域內(nèi)，得到滿足要求的近似解。當(dāng)收斂精度高或者系統(tǒng)規(guī)模比較大的時候，HEM將消耗較多時間來進(jìn)行帕德近似計算[6]，而高階帕德近似計算可能導(dǎo)致近似解停滯，無法進(jìn)一步收斂到最優(yōu)[14]。

針對上述問題，本文提出一種用于求解電力系統(tǒng)經(jīng)濟(jì)調(diào)度的重啟動全純函數(shù)嵌入式算法(restarted holomorphic embedding method，RHE)，在HEM中引入重啟動機(jī)制和啟發(fā)式規(guī)則，使算法僅涉及低階帕德近似計算，并對不等式約束進(jìn)行篩選，減小計算規(guī)模，同時確定合適的初始值以充分發(fā)揮HEM的性能。最后，結(jié)合不同規(guī)模的靜態(tài)和動態(tài)經(jīng)濟(jì)調(diào)度問題，對所提出的RHE進(jìn)行性能測試。

1 經(jīng)濟(jì)調(diào)度模型

假設(shè)求解的電力系統(tǒng)包含NB個節(jié)點(diǎn)、NG個發(fā)電機(jī)組和NF條輸電線路。

1.1 目標(biāo)函數(shù)

經(jīng)濟(jì)調(diào)度以最小化調(diào)度期間的發(fā)電機(jī)組總?cè)剂铣杀綟為目標(biāo)[15]，即

(1)

式中：t為調(diào)度時段；T為調(diào)度時段個數(shù)；Q、C、k分別為發(fā)電成本的二次項、一次項和常數(shù)項系數(shù)；Pt為發(fā)電機(jī)組在時段t的有功出力(下標(biāo)t表示時段t，下同)。

1.2 約束條件

經(jīng)濟(jì)調(diào)度問題的約束條件可以表示為[16]：

(2)

式中：Bbus、BF分別為節(jié)點(diǎn)導(dǎo)納矩陣、支路導(dǎo)納矩陣；列向量Θt、Dt分別為節(jié)點(diǎn)電壓相角和有功負(fù)荷；CG為發(fā)電機(jī)關(guān)聯(lián)矩陣，當(dāng)發(fā)電機(jī)g連接在節(jié)點(diǎn)h，其元素(g,h)為1，否則為0；Fmax為線路允許傳輸功率上限；Pmin、Pmax分別為發(fā)電機(jī)組出力的下限、上限；Rt為發(fā)電機(jī)組的最大爬坡功率。式(2a)為基于直流潮流的功率平衡方程，式(2b)為線路的傳輸約束，式(2c)、(2d)分別為發(fā)電機(jī)組的出力約束和爬坡速率約束。

相比于僅考慮系統(tǒng)輸出總有功功率與總負(fù)荷平衡關(guān)系的經(jīng)濟(jì)調(diào)度模型[18]，本文采用基于直流潮流的功率平衡方程，考慮了電網(wǎng)架構(gòu)對潮流分布的影響，模型更準(zhǔn)確。

1.3 模型求解

通過引入非負(fù)的松弛變量Z將不等式約束轉(zhuǎn)換為等式約束，得到經(jīng)濟(jì)調(diào)度問題形式為：

(3)

式中：λt、μt分別為對應(yīng)約束條件的等式、不等式拉格朗日乘子；下標(biāo)F、G、R分別表示線路傳輸約束、發(fā)電機(jī)容量約束以及爬坡約束，上標(biāo)+、-分別表示約束的上限、下限，下同。

對于時段t，令Ht為等式約束、Gt為不等式約束，則經(jīng)濟(jì)調(diào)度問題的拉格朗日函數(shù)

(4)

問題對應(yīng)的KKT方程為：

(5)

2 HEM

對于KKT方程，可以利用HEM求解，主要包含4個步驟。

2.1 步驟1：構(gòu)造嵌入式方程

在式(5)中嵌入復(fù)變量s，構(gòu)造全純的嵌入式方程為

L(X(s))=L(X(0))+s(RHS-L(X(0))).

(6)

式中X(s)為關(guān)于s的函數(shù)。s=0時，式(6)恒成立，X(0)為可以任意選取的變量初值。

X(s)可以展開成麥克勞林級數(shù)[19]：

(7)

式中X[k]為級數(shù)的第k階系數(shù)?？梢酝茖?dǎo)出，變量級數(shù)的常數(shù)項等于變量初值，即X[0]=X(0)。

2.2 步驟2：求解冪級數(shù)系數(shù)

將式(7)代入到式(6)中，可以推導(dǎo)出線性的遞歸方程為

AX[n]=B[n-1].

(8)

式中：X[n]為變量冪級數(shù)的第n階系數(shù)；其余項僅涉及前n-1階系數(shù)，移到等式右側(cè)構(gòu)成右端項B[n-1]；A為系數(shù)矩陣，僅包含常數(shù)項，不會隨著n改變，其表達(dá)式如式(9)所示，其中μt、Zt取初值，D(a)表示將向量a對角化，EG為單位矩陣。當(dāng)n=1時，B[0]=RHS-L(X(0))；當(dāng)n≥2時，右端項B[n-1]如式(10)所示。

(9)

(10)

2.3 步驟3：計算帕德近似

為了將冪級數(shù)的收斂域擴(kuò)展到最大，解析延拓一般選取對角或近似對角的帕德近似，即

(11)

2.4 步驟4：判斷收斂條件

采用常規(guī)的方程組最大不匹配度Mi來判斷收斂，對于不等式約束只需判斷是否滿足條件，收斂判據(jù)為

max(Gt),‖μtZt‖∞)≤ε.

(12)

(13)

如果近似解滿足收斂條件，將其作為原問題的解輸出；如果不滿足收斂條件，則返回步驟2計算更高階的系數(shù)，重新計算新的帕德近似值。

3 RHE

HEM的運(yùn)行時間主要消耗在帕德近似計算，當(dāng)收斂精度高或者計算大規(guī)模系統(tǒng)時，帕德近似需要計算到非常高階，不僅耗時長，還可能在不滿足收斂精度的近似解停滯不動[14]。動態(tài)經(jīng)濟(jì)調(diào)度問題需要優(yōu)化整個調(diào)度時段的機(jī)組出力，涉及的計算量很大，直接利用HEM進(jìn)行求解時，求解速度可能很慢。考慮到在求解過程中，有部分不等式約束在多數(shù)情況下不起作用，可以先不考慮它們，放開相應(yīng)的限制，求解一個可行范圍更大但計算規(guī)模較小的新問題。通過降低問題的求解規(guī)模，提高算法效率。

本文在HEM中引入重啟動機(jī)制和啟發(fā)式規(guī)則，提出了用于求解動態(tài)經(jīng)濟(jì)調(diào)度問題的啟發(fā)式全純函數(shù)嵌入算法。RHE設(shè)定冪級數(shù)系數(shù)的最高次數(shù)nmax，將計算過程分為初始過程和重啟動過程。初始過程中，當(dāng)系數(shù)計算到最高次數(shù)算法仍未收斂時，停止計算更高階系數(shù)，進(jìn)入重啟動過程。重啟動過程中啟發(fā)式規(guī)則利用上一過程得到的近似解，篩除優(yōu)化問題的不等式約束，簡化當(dāng)前過程中需要求解的問題，并為構(gòu)建的新問題選取合適的初始值。

3.1 初始過程

(14)

將近似解代入式(12)判斷是否收斂，不滿足時算法進(jìn)入重啟動過程。

3.2 重啟動過程

當(dāng)I為空集時，重啟動過程暫時無需考慮額外的不等式約束，變量初始值選取為上一過程的近似值。

將考慮的不等式約束和對應(yīng)的互補(bǔ)條件添加到方程組式(14)中，構(gòu)成當(dāng)前過程求解的方程：

(15)

確定好方程和初值后，式(15)可依據(jù)第2章的HEM進(jìn)行求解。構(gòu)造與式(6)同形式的嵌入式方程，則推導(dǎo)出遞歸方程AactX[n]=Bact，其中系數(shù)矩陣

(16)

(17)

右端項Bact由式(10)的對應(yīng)行組成，

(18)

在重啟動過程中，障礙因子不斷更新，逐步逼近0，即

(19)

為了防止問題無解時算法一直發(fā)生重啟，需設(shè)定最大重啟動次數(shù)Nmax。當(dāng)重啟動次數(shù)達(dá)到最大依舊未計算到收斂解，輸出結(jié)果不收斂。RHE無需事先確定初始值，初始過程的線性方程組可以直接求解，重啟動過程的初值由啟發(fā)式規(guī)則確定。

重啟動機(jī)制使該算法只涉及到低階帕德近似計算，避免了高階帕德近似計算的高耗時和停滯解問題；啟發(fā)式規(guī)則能夠根據(jù)上一過程的近似解，篩選可能起作用的約束加入方程組進(jìn)行求解，同時確定合適的初始點(diǎn)，降低了解方程的規(guī)模，提高了HEM的計算效率。

RHE的計算流程如圖1所示。

圖1 RHE的計算流程

4 算例測試

為了驗(yàn)證RHE的性能，結(jié)合不同規(guī)模的靜態(tài)和動態(tài)經(jīng)濟(jì)調(diào)度問題進(jìn)行測試。所有測試都在MATLAB平臺上建模實(shí)現(xiàn)，并在硬件配置為Intel 3.0 GHz四核處理器和8 GB RAM的計算機(jī)上進(jìn)行仿真。為了驗(yàn)證本文提出方法的準(zhǔn)確性，以MATPOWER包的內(nèi)點(diǎn)法(MATPOWER interior point solver，MIPS)作為基準(zhǔn)[21]，同時結(jié)合基于傳統(tǒng)經(jīng)濟(jì)調(diào)度的HEM(HE method bases on the primitive economic dispatch，EDHE)進(jìn)行分析[22]。設(shè)收斂精度ε=10-5，冪級數(shù)系數(shù)最高次數(shù)nmax=15，重啟動最高次數(shù)Nmax=5。

4.1 靜態(tài)經(jīng)濟(jì)調(diào)度

為了研究RHE的有效性和重啟動機(jī)制對算法的影響，利用RHE求解不考慮爬坡約束的靜態(tài)經(jīng)濟(jì)調(diào)度問題。為了便于與EHDE對比，使用與文獻(xiàn)[22]相同的算例進(jìn)行測試，系統(tǒng)的信息匯總見表1。每新增3個冪級數(shù)系數(shù)進(jìn)行1次帕德近似，視為1步遞歸計算。

表1 靜態(tài)經(jīng)濟(jì)調(diào)度的測試算例

圖2、圖3分別是RHE和EDHE的不匹配度變化情況。2種方法都能在10步內(nèi)達(dá)到收斂精度，RHE的步數(shù)明顯少于EDHE，說明重啟動機(jī)制能有效提高收斂速度。

圖3 EDHE的不匹配度

比較RHE和MIPS、EDHE在發(fā)電成本和有功出力上的差異，結(jié)果見表2、表3。發(fā)電成本的最大誤差在4.1×10-5美元以內(nèi)，有功出力誤差不超過5.6×10-8MW，這個計算誤差可以忽略不計，說明3種方法都能準(zhǔn)確收斂到同一最優(yōu)解。

表2 RHE與MIPS的計算誤差

表3 RHE與EDHE的計算誤差

表4統(tǒng)計了RHE、EDHE、MIPS求解不同算例的執(zhí)行時間。RHE的耗時基本少于EDHE和MIPS。求解57節(jié)點(diǎn)系統(tǒng)時，RHE耗時略高于EDHE，但相較于MIPS，兩者的執(zhí)行時間都較短。相比于其他4個算例，RHE計算IEEE 118節(jié)點(diǎn)系統(tǒng)的時間明顯更長，這是因?yàn)樵谇蠼膺^程中多執(zhí)行了一個重啟動過程，但耗時依舊少于EDHE。

通過與MIPS、EDHE對比，驗(yàn)證了RHE的計算準(zhǔn)確性和良好的收斂特性。引入重啟動機(jī)制，明顯提高了收斂速度，減少了算法的執(zhí)行時間。

4.2 動態(tài)經(jīng)濟(jì)調(diào)度

為了進(jìn)一步驗(yàn)證RHE的性能，結(jié)合改進(jìn)IEEE 14節(jié)點(diǎn)和300節(jié)點(diǎn)系統(tǒng)的動態(tài)經(jīng)濟(jì)調(diào)度問題進(jìn)行測試，考慮調(diào)度時長為24個時段(1 h為1個時段)的調(diào)度問題。發(fā)電機(jī)成本參數(shù)和負(fù)荷需求見附錄A。

附錄A

4.2.1 IEEE 14節(jié)點(diǎn)系統(tǒng)

在IEEE 14節(jié)點(diǎn)系統(tǒng)中，為了測試RHE求解動態(tài)經(jīng)濟(jì)調(diào)度問題的性能，設(shè)計3種爬坡約束限制，分別為發(fā)電機(jī)組當(dāng)前容量的25%(算例1)、18%(算例2)、10%(算例3)。

利用RHE和MIPS分別求解以上3個算例，仿真結(jié)果表明2種方法的最大不匹配度都能達(dá)到收斂精度要求。表5給出RHE與MIPS的計算誤差，發(fā)電成本的誤差在3.1×10-6美元以內(nèi)，有功出力的誤差不超過2.0×10-6MW，說明2種方法成功收斂到同一個最優(yōu)解。

表5 RHE與MIPS的計算誤差(算例1、2、3)

各個發(fā)電機(jī)組在每個調(diào)度時段的有功出力如圖4所示。算例1中，發(fā)電機(jī)組的爬坡速率限制足夠?qū)捤?，沒有動態(tài)約束起作用，能夠根據(jù)負(fù)荷需求以最經(jīng)濟(jì)的方式進(jìn)行功率分配。算例2中，機(jī)組4的爬坡約束在5個時段起作用，機(jī)組5的爬坡約束只在1個時段起作用。算例3中，5個機(jī)組的爬坡約束均有起作用。表6列出了機(jī)組爬坡約束起作用的具體時段。

表6 發(fā)電機(jī)組起作用的爬坡約束

隨著發(fā)電機(jī)組爬坡約束限制不斷收緊，相鄰時段內(nèi)功率可變化的范圍縮窄。為了保持系統(tǒng)功率平衡，功率分配占比高的機(jī)組的出力曲線波動較大，峰谷更加明顯，如機(jī)組2；功率分配占比低的機(jī)組的出力曲線相對平緩，如機(jī)組4。動態(tài)約束起作用后，系統(tǒng)功率分配的經(jīng)濟(jì)性下降，發(fā)電成本有所上升，從算例1到算例3的成本逐漸升高。

接下來分析2種方法的收斂速度。圖5給出了RHE、MIPS算法的最大不匹配度隨遞歸步數(shù)、迭代步數(shù)的變化情況。隨著動態(tài)約束起作用，2種方法所需的計算步數(shù)都會增多，其中RHE的變化更明顯：MIPS需要8～10步才能收斂；RHE的步數(shù)從不到4步增加到7步?？傮w表現(xiàn)上，RHE具有更好的收斂速度。

圖5 IEEE 14節(jié)點(diǎn)系統(tǒng)最大不匹配度

RHE和MIPS計算3個算例的執(zhí)行時間見表7?？紤]到帕德近似計算的高度可并行性，可以從總執(zhí)行時間中剔除其計算時間，得到調(diào)整時間，這個時長更具有代表性[6]。無論忽略與否，RHE的耗時都明顯少于MIPS；此外，RHE求解算例3的耗時明顯增加，是因?yàn)槠鹱饔玫膭討B(tài)約束增多，需要更多重啟動過程才能收斂。

表7 IEEE 14節(jié)點(diǎn)系統(tǒng)的執(zhí)行時間

4.2.2 IEEE 300節(jié)點(diǎn)系統(tǒng)

在IEEE 300節(jié)點(diǎn)系統(tǒng)中，同樣設(shè)計3個不同爬坡約束限制的算例，測試RHE求解大規(guī)模系統(tǒng)的性能。速率限制分別設(shè)定為發(fā)電機(jī)組當(dāng)前容量的25%(算例4)、20%(算例5)、10%(算例6)。

RHE和MIPS的計算誤差見表8。雖然算例6的發(fā)電成本誤差數(shù)量級在10-4美元，但相對于數(shù)量級在107美元的總成本，這樣的誤差可以忽略不計。2種方法的有功出力誤差也非常小，它們成功收斂到同一個最優(yōu)解。

表8 RHE與MIPS的計算誤差(算例4、5、6)

表9列出了機(jī)組爬坡約束起作用的時段。25%爬坡速率限制下，測試系統(tǒng)中有13個機(jī)組的爬坡約束起作用，在無限制下速率可能超過25%，影響到發(fā)電機(jī)的運(yùn)行壽命，實(shí)際運(yùn)行中有必要對爬坡速率進(jìn)行限制。隨著爬坡約束的收緊，約束起作用的機(jī)組明顯增多，算例6中只有13個機(jī)組沒有受限；同一機(jī)組約束起作用的時間也逐步延長。

表9 發(fā)電機(jī)組起作用的爬坡約束

2種方法的最大不匹配度如圖6所示。MIPS達(dá)到收斂的計算步數(shù)主要分布在8～10步或16～20步；RHE的步數(shù)分布較分散，主要在4步或者16步左右，需要的步數(shù)一般少于MIPS，表現(xiàn)出更好的收斂速度。

圖6 IEEE 300節(jié)點(diǎn)系統(tǒng)最大不匹配度

2種方法的執(zhí)行時間見表10。RHE的總執(zhí)行時間多于MIPS，但是其忽略帕德近似計算后的調(diào)整時間略少于MIPS，兩者相差不大，這也說明帕德近似計算占據(jù)了HEM的大部執(zhí)行時間(尤其系統(tǒng)規(guī)模較大時)。RHE采用的是常規(guī)的帕德近似方法，已有文獻(xiàn)提出更高效的帕德近似方法[23]，計算時長可以壓縮到原來的1/6以下。考慮并行運(yùn)算時，仍然可以采用調(diào)整時間作為近似耗時。

表10 IEEE 300節(jié)點(diǎn)系統(tǒng)的執(zhí)行時間

5 結(jié)束語

本文提出了一種求解電力系統(tǒng)經(jīng)濟(jì)調(diào)度問題的RHE。在常規(guī)全純嵌入算法中引入了重啟動機(jī)制和啟發(fā)式規(guī)則，該算法只涉及到低階帕德近似計算，減少了計算高階帕德近似的時間，同時避開了近似解可能停滯的情況；在每個重啟動過程中，篩選可能起作用的約束加入方程組進(jìn)行求解，無需計算完整的KKT方程，降低了方程的求解規(guī)模，同時可以為HEM求解方程提供合適的初始點(diǎn)。

經(jīng)過測試，驗(yàn)證了該算法能高效、準(zhǔn)確地求解靜態(tài)和動態(tài)經(jīng)濟(jì)調(diào)度問題。求解靜態(tài)經(jīng)濟(jì)調(diào)度問題時，RHE的表現(xiàn)優(yōu)于EDHE和MIPS。在求解動態(tài)經(jīng)濟(jì)調(diào)度問題方面，與MIPS相比：計算小規(guī)模系統(tǒng)時，RHE在收斂速度和執(zhí)行時間上具有明顯的優(yōu)勢；計算大規(guī)模系統(tǒng)時，RHE需要消耗很多時間去計算帕德近似，不考慮這部分可并行運(yùn)算的時間，RHE在收斂速度和執(zhí)行時間上仍有優(yōu)勢。

未來的工作可以著眼于：①使用更準(zhǔn)確的交流潮流模型，考慮合約出力約束等更多動態(tài)約束；②提高算法的計算性能，開發(fā)更高效的帕德近似計算方法，實(shí)現(xiàn)并行運(yùn)算。