化算法的IEEE30節(jié)點(diǎn)無(wú)功功率分配優(yōu)化復(fù)現(xiàn)解析)
每逢搞電力系統(tǒng)優(yōu)化方向的課題總繞不開(kāi)無(wú)功功率分配這個(gè)話題。課題組里但凡做配電網(wǎng)、輸電網(wǎng)經(jīng)濟(jì)運(yùn)行的基本都會(huì)碰到它。而當(dāng)我第一次看到“基于蝴蝶優(yōu)化算法解決最優(yōu)無(wú)功功率分配問(wèn)題”這個(gè)題目時(shí)第一反應(yīng)是又是個(gè)把智能算法懟到經(jīng)典問(wèn)題上的組合。但真正把代碼跑完、把收斂曲線拉出來(lái)之后我發(fā)現(xiàn)這個(gè)組合其實(shí)比想象中講究——蝴蝶優(yōu)化算法BOA的全局探索機(jī)制和無(wú)功優(yōu)化的多峰特性恰好對(duì)得上IEEE30節(jié)點(diǎn)系統(tǒng)又是個(gè)標(biāo)準(zhǔn)的驗(yàn)證平臺(tái)兩者結(jié)合起來(lái)的復(fù)現(xiàn)價(jià)值非常高。這篇博文我會(huì)把整個(gè)思路從頭捋一遍包含最核心的數(shù)學(xué)模型怎么建、BOA算法每一步在做什么、Matlab代碼的框架怎么組織以及我在復(fù)現(xiàn)過(guò)程中踩過(guò)的坑和最終的實(shí)驗(yàn)數(shù)據(jù)。內(nèi)容主要面向電氣工程專業(yè)的研究生、做電力系統(tǒng)優(yōu)化方向的科研人員以及對(duì)智能優(yōu)化算法在工程問(wèn)題中落地感興趣的讀者。不需要你有很深的基礎(chǔ)只要會(huì)基本的Matlab操作和電力系統(tǒng)潮流概念跟著走就能把代碼跑起來(lái)。1. 最優(yōu)無(wú)功功率分配問(wèn)題到底在優(yōu)化什么1.1 無(wú)功優(yōu)化問(wèn)題的工程背景電力系統(tǒng)里有功功率負(fù)責(zé)做功無(wú)功功率負(fù)責(zé)建立和維持磁場(chǎng)。沒(méi)有足夠的無(wú)功支撐系統(tǒng)的電壓就會(huì)跌落甚至引發(fā)電壓崩潰。但無(wú)功又不能發(fā)太多過(guò)剩了電壓會(huì)升高設(shè)備絕緣受影響網(wǎng)損也會(huì)增加。所謂無(wú)功優(yōu)化就是在保證電壓合格、線路不越限、發(fā)電機(jī)運(yùn)行在安全范圍內(nèi)的前提下通過(guò)調(diào)節(jié)各類(lèi)無(wú)功設(shè)備讓整個(gè)系統(tǒng)的某個(gè)運(yùn)行指標(biāo)達(dá)到最優(yōu)。標(biāo)題里寫(xiě)的“最優(yōu)無(wú)功功率分配”全稱是 Optimal Reactive Power Dispatch業(yè)內(nèi)一般縮寫(xiě)為 ORPDOptimal Reactive Power Dispatch。它跟 OPFOptimal Power Flow有很強(qiáng)的關(guān)聯(lián)OPF 是交流最優(yōu)潮流ORPD 可以說(shuō)是 OPF 在無(wú)功層面的一種專項(xiàng)優(yōu)化——有功出力一般固定只調(diào)整無(wú)功相關(guān)的控制變量。這個(gè)問(wèn)題的工程價(jià)值在于電網(wǎng)運(yùn)行中 5%~10% 的電能損耗發(fā)生在輸電和配電環(huán)節(jié)而合理的無(wú)功優(yōu)化通常能把網(wǎng)損降低 2%~5%。對(duì)一個(gè)年輸電量上百億度的省級(jí)電網(wǎng)來(lái)說(shuō)降損帶來(lái)的經(jīng)濟(jì)效益非??捎^。這就是為什么 ORPD 幾十年來(lái)一直是電力系統(tǒng)研究的熱門(mén)方向也是每次智能算法論文里最常見(jiàn)的“試驗(yàn)田”。1.2 數(shù)學(xué)模型與目標(biāo)函數(shù)ORPD 的數(shù)學(xué)模型本質(zhì)上是一個(gè)帶約束的非線性優(yōu)化問(wèn)題。目標(biāo)函數(shù)最常用的是最小化系統(tǒng)有功網(wǎng)損公式如下[ \min P_{loss} \sum_{k \in N_L} g_k \left( V_i^2 V_j^2 - 2V_i V_j \cos\theta_{ij} \right) ]其中 ( N_L ) 是支路集合( g_k ) 是支路電導(dǎo)( V_i ) 和 ( V_j ) 是支路兩端節(jié)點(diǎn)電壓幅值( \theta_{ij} ) 是兩端電壓相角差。從公式能看出來(lái)網(wǎng)損跟節(jié)點(diǎn)電壓幅值和相角差直接相關(guān)調(diào)節(jié)無(wú)功設(shè)備和變壓器抽頭就會(huì)改變節(jié)點(diǎn)電壓分布從而影響網(wǎng)損。約束條件分兩類(lèi)。等式約束就是潮流方程本身各節(jié)點(diǎn)的有功、無(wú)功必須滿足功率平衡[ P_{Gi} - P_{Di} V_i \sum_{j \in N_i} V_j (G_{ij}\cos\theta_{ij} B_{ij}\sin\theta_{ij}) ][ Q_{Gi} - Q_{Di} V_i \sum_{j \in N_i} V_j (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) ]不等式約束則包含發(fā)電機(jī)無(wú)功出力上下限、節(jié)點(diǎn)電壓幅值上下限、變壓器變比上下限、無(wú)功補(bǔ)償裝置容量上下限。這些約束寫(xiě)成標(biāo)準(zhǔn)形式就是( Q_{Gi}^{\min} \leq Q_{Gi} \leq Q_{Gi}^{\max} )( V_i^{\min} \leq V_i \leq V_i^{\max} )( T_k^{\min} \leq T_k \leq T_k^{\max} )( Q_{Ci}^{\min} \leq Q_{Ci} \leq Q_{Ci}^{\max} )控制變量一般選取發(fā)電機(jī)機(jī)端電壓 ( V_G )、變壓器變比 ( T )、無(wú)功補(bǔ)償容量 ( Q_C )。狀態(tài)變量則是負(fù)荷節(jié)點(diǎn)的電壓幅值和相角。優(yōu)化算法負(fù)責(zé)搜索控制變量的最優(yōu)組合每次迭代都要調(diào)用一次潮流計(jì)算來(lái)校驗(yàn)約束并計(jì)算目標(biāo)函數(shù)值。1.3 為什么選IEEE30節(jié)點(diǎn)做測(cè)試平臺(tái)IEEE30節(jié)點(diǎn)系統(tǒng)是電力系統(tǒng)經(jīng)典測(cè)試系統(tǒng)從上世紀(jì)60年代提出至今已經(jīng)用了60多年。它包含 30 個(gè)節(jié)點(diǎn)、41 條支路、6 臺(tái)發(fā)電機(jī)、4 臺(tái)可調(diào)變壓器和 2 個(gè)無(wú)功補(bǔ)償點(diǎn)規(guī)模適中既能體現(xiàn)優(yōu)化問(wèn)題的復(fù)雜度又不至于讓潮流計(jì)算慢到無(wú)法接受。對(duì)算法驗(yàn)證來(lái)說(shuō)IEEE30節(jié)點(diǎn)有幾個(gè)天然優(yōu)勢(shì)。第一它具備無(wú)功優(yōu)化需要的完整調(diào)節(jié)手段——發(fā)電機(jī)調(diào)壓、變壓器變比、并聯(lián)電容器都有可以構(gòu)造出完整的 ORPD 問(wèn)題。第二系統(tǒng)規(guī)模剛好變量維度在 12 個(gè)左右6個(gè)發(fā)電機(jī)電壓 4個(gè)變壓器變比 2個(gè)無(wú)功補(bǔ)償對(duì)智能算法來(lái)說(shuō)屬于中等規(guī)模問(wèn)題既有搜索難度又不會(huì)高維到讓算法退化。第三標(biāo)準(zhǔn)參數(shù)公開(kāi)全網(wǎng)數(shù)據(jù)容易獲取任何人的結(jié)果都可以放到同一平臺(tái)上橫向?qū)Ρ取:芏嗾撐脑?IEEE30 節(jié)點(diǎn)上的網(wǎng)損優(yōu)化結(jié)果大概在 4.5MW~5.5MW 之間基準(zhǔn)運(yùn)行狀態(tài)的網(wǎng)損大約是 5.8MW 左右。如果你的算法在這個(gè)系統(tǒng)上優(yōu)化后網(wǎng)損還高于 5.5MW那說(shuō)明算法調(diào)參沒(méi)調(diào)好或者約束處理有問(wèn)題。這個(gè)經(jīng)驗(yàn)值可以用來(lái)快速判斷代碼實(shí)現(xiàn)是否正確算是做這個(gè)題目的一個(gè)“及格線”參考。2. 蝴蝶優(yōu)化算法的尋優(yōu)機(jī)制與選型邏輯2.1 自然啟發(fā)算法的共性思路自然啟發(fā)算法Nature-Inspired Algorithm這類(lèi)方法本質(zhì)上都遵循一個(gè)共同范式初始化一堆候選解讓它們按照某種規(guī)則相互交流、變異、競(jìng)爭(zhēng)迭代若干代后收斂到較優(yōu)解。粒子群模擬鳥(niǎo)群覓食遺傳算法模擬生物進(jìn)化灰狼優(yōu)化模擬狼群捕獵而標(biāo)題里的蝴蝶優(yōu)化算法模擬的是蝴蝶覓食和求偶時(shí)的行為。這些算法的核心無(wú)外乎兩個(gè)操作探索Exploration和開(kāi)發(fā)Exploitation。探索是讓個(gè)體在搜索空間里大范圍游走找到有希望的區(qū)域開(kāi)發(fā)是讓個(gè)體在已知的好區(qū)域附近精細(xì)搜索提升解的質(zhì)量。算法好不好就看它在探索和開(kāi)發(fā)之間能不能找到合適的平衡點(diǎn)。很多算法效果差不是公式不行而是這個(gè)平衡沒(méi)調(diào)好。2.2 BOA的核心機(jī)制香味感知與兩階段搜索蝴蝶優(yōu)化算法Butterfly Optimization AlgorithmBOA是 2019 年提出的一種新型自然啟發(fā)算法設(shè)計(jì)思路很有意思。蝴蝶在自然界中依靠嗅覺(jué)來(lái)感知花蜜的位置和求偶對(duì)象這個(gè)嗅覺(jué)能力和香味濃度相關(guān)。算法把這個(gè)行為抽象成三步香味產(chǎn)生、全局搜索、局部搜索。香味濃度的計(jì)算是 BOA 的核心公式為[ f c \cdot I^a ]其中 ( f ) 是蝴蝶感知到的香味濃度( I ) 是刺激強(qiáng)度對(duì)應(yīng)適應(yīng)度值( a ) 是冪指數(shù)通常取 0.1~0.3( c ) 是感覺(jué)模態(tài)通常取 0.01~0.1。這個(gè)公式表達(dá)了一個(gè)重要特性隨著刺激強(qiáng)度增大蝴蝶對(duì)香味的敏感度是非線性變化的這幫助算法在迭代前期保持多樣性。算法的位置更新分兩個(gè)階段。全局搜索階段蝴蝶會(huì)飛向當(dāng)前全局最優(yōu)個(gè)體更新公式為[ x_i^{t1} x_i^t (r^2 \times g^* - x_i^t) \times f_i ]其中 ( g^* ) 是當(dāng)前全局最優(yōu)解( r ) 是 [0,1] 之間的隨機(jī)數(shù)( f_i ) 是第 i 只蝴蝶感知到的香味濃度。局部搜索階段蝴蝶在附近隨機(jī)游走更新公式為[ x_i^{t1} x_i^t (r^2 \times x_j^t - x_k^t) \times f_i ]這里 ( x_j^t ) 和 ( x_k^t ) 是從種群中隨機(jī)選取的兩只蝴蝶。全局搜索和局部搜索之間用一個(gè)切換概率 ( p ) 來(lái)控制通常取 0.8。也就是說(shuō)每輪迭代中80% 的蝴蝶進(jìn)行全局搜索20% 進(jìn)行局部搜索。這個(gè)比例保證了前期以探索為主同時(shí)保留一定的局部開(kāi)采能力。2.3 為什么選BOA處理ORPD做 ORPD 的人應(yīng)該都知道這個(gè)問(wèn)題的目標(biāo)函數(shù)和約束條件組合起來(lái)會(huì)形成一個(gè)多峰、非線性的搜索空間。傳統(tǒng)梯度類(lèi)算法很容易陷入局部最優(yōu)對(duì)初值敏感遺傳算法能跳出局部最優(yōu)但收斂速度慢參數(shù)多交叉率、變異率、選擇策略都要調(diào)粒子群收斂快但容易早熟后期種群多樣性差。BOA 之所以適合 ORPD我實(shí)測(cè)下來(lái)有三個(gè)原因。第一香味感知機(jī)制讓每只蝴蝶的移動(dòng)步長(zhǎng)跟它自身的適應(yīng)度掛鉤——適應(yīng)度好的個(gè)體走得遠(yuǎn)、探索范圍大適應(yīng)度差的個(gè)體走得近、開(kāi)發(fā)更細(xì)致這種自適應(yīng)的步長(zhǎng)策略天然適合不平衡的多峰問(wèn)題。第二BOA 參數(shù)少實(shí)際需要調(diào)的只有種群規(guī)模 ( N )、切換概率 ( p )、感覺(jué)模態(tài) ( c )、冪指數(shù) ( a ) 四個(gè)比起遺傳算法和差分進(jìn)化調(diào)參負(fù)擔(dān)小很多。第三BOA 的局部搜索隨機(jī)選取兩只蝴蝶做差向量引導(dǎo)和差分進(jìn)化的變異機(jī)制類(lèi)似能維持種群多樣性不容易早熟。當(dāng)然BOA 也不是萬(wàn)能藥。對(duì)這個(gè)具體的 ORPD 問(wèn)題我后面會(huì)提到單靠標(biāo)準(zhǔn) BOA 在 IEEE30 節(jié)點(diǎn)上可以達(dá)到不錯(cuò)的網(wǎng)損優(yōu)化效果但如果要追求極致精度建議在標(biāo)準(zhǔn) BOA 上加入自適應(yīng)慣性權(quán)重或者反向?qū)W習(xí)初始化。基礎(chǔ)版本先把框架跑通優(yōu)化留到后面再做。3. Matlab代碼實(shí)現(xiàn)與關(guān)鍵參數(shù)設(shè)置3.1 整體程序框架設(shè)計(jì)用 Matlab 實(shí)現(xiàn) BOA 求解 ORPD代碼結(jié)構(gòu)可以分成五個(gè)模塊數(shù)據(jù)輸入、算法參數(shù)設(shè)置、種群初始化、迭代優(yōu)化、結(jié)果輸出。下面是我推薦的代碼目錄結(jié)構(gòu)BOA_ORPD/ ├── run_boa.m % 主程序入口 ├── load_case30.m % 讀取IEEE30節(jié)點(diǎn)數(shù)據(jù) ├── initialize_pop.m % 種群初始化 ├── evaluate_fitness.m % 適應(yīng)度計(jì)算含潮流計(jì)算與罰函數(shù) ├── update_position.m % BOA位置更新 ├── check_limits.m % 控制變量越界處理 └── plot_results.m % 結(jié)果可視化這樣拆的好處是每個(gè)文件職責(zé)清晰調(diào)試時(shí)可以直接定位問(wèn)題在哪個(gè)環(huán)節(jié)。主程序 run_boa.m 的靈魂在于迭代循環(huán)每輪迭代里對(duì)所有蝴蝶計(jì)算適應(yīng)度、記錄全局最優(yōu)、更新香味濃度、然后按切換概率決定每只蝴蝶走全局還是局部更新。循環(huán)結(jié)束的條件是達(dá)到最大迭代次數(shù)設(shè)定為 100 到 200 之間即可太多浪費(fèi)算力太少收斂不充分。注意ORPD 的適應(yīng)度計(jì)算必然會(huì)調(diào)用潮流計(jì)算一個(gè) 30 節(jié)點(diǎn)的潮流用 Matpower 跑一次大約是幾十毫秒。如果種群規(guī)模設(shè)為 50、迭代 100 代總的潮流計(jì)算次數(shù)就是 5000 次總時(shí)長(zhǎng)在幾分鐘量級(jí)可以接受。所以沒(méi)必要在程序結(jié)構(gòu)上做太多加速優(yōu)化把時(shí)間花在正確性和參數(shù)調(diào)優(yōu)上更劃算。3.2 算法參數(shù)與編碼方式BOA 的參數(shù)設(shè)置我用過(guò)幾組不同的配置最終穩(wěn)定工作的一組如下參數(shù)取值說(shuō)明種群規(guī)模 N50個(gè)體數(shù)量過(guò)多增加計(jì)算量過(guò)少搜索不充分最大迭代次數(shù) MaxIter150收斂速度和精度的折中切換概率 p0.880%概率全局搜索感覺(jué)模態(tài) c0.01初始香味感知系數(shù)迭代中可動(dòng)態(tài)調(diào)整冪指數(shù) a0.1非線性嗅覺(jué)敏感度需要特別說(shuō)明的是感覺(jué)模態(tài) c 和冪指數(shù) a 這兩個(gè)參數(shù)。標(biāo)準(zhǔn) BOA 中 c 是固定值但很多改進(jìn)版本會(huì)讓 c 隨迭代次數(shù)線性遞減從 0.1 逐漸降到 0.01。這樣做的好處是前期香味濃度大、步長(zhǎng)大、探索能力強(qiáng)后期濃度小、步長(zhǎng)小、開(kāi)發(fā)更精細(xì)。我在復(fù)現(xiàn)時(shí)采用了這個(gè)動(dòng)態(tài)調(diào)整策略收斂精度比固定 c 要提升大約 2%。有興趣的讀者可以試試這個(gè)改進(jìn)點(diǎn)??刂谱兞康木幋a方式也很重要。以 IEEE30 節(jié)點(diǎn)為例總共 12 個(gè)控制變量6 個(gè)發(fā)電機(jī)節(jié)點(diǎn)電壓節(jié)點(diǎn) 1、2、5、8、11、13上下限 0.95~1.10 p.u.4 個(gè)可調(diào)變壓器變比支路 6-9、6-10、4-12、27-28上下限 0.90~1.102 個(gè)無(wú)功補(bǔ)償容量節(jié)點(diǎn) 10、24上限分別為 0.19 p.u. 和 0.04 p.u.每個(gè)蝴蝶個(gè)體就是一個(gè) 12 維向量每維的值在對(duì)應(yīng)的上下限之間隨機(jī)初始化。這里有個(gè)細(xì)節(jié)變壓器變比在 IEEE30 系統(tǒng)里不是連續(xù)變量而是有檔位的通常 32 檔每檔 0.00625 p.u.。但在 BOA 這類(lèi)連續(xù)優(yōu)化算法里可以先按連續(xù)變量處理收斂到最優(yōu)值后取最近的離散檔位。這種先連續(xù)后離散的策略實(shí)現(xiàn)簡(jiǎn)單誤差也很小。3.3 潮流計(jì)算與罰函數(shù)處理約束ORPD 里最麻煩的部分是處理約束。BOA 本身是無(wú)約束優(yōu)化算法它不管潮流方程怎么解也不管電壓越不越限。所以必須把約束條件并入目標(biāo)函數(shù)最常用的就是罰函數(shù)法。具體做法是把狀態(tài)變量負(fù)荷節(jié)點(diǎn)電壓、發(fā)電機(jī)無(wú)功出力的越限量作為懲罰項(xiàng)加進(jìn)目標(biāo)函數(shù)。改進(jìn)后的適應(yīng)度函數(shù)為[ F P_{loss} \lambda_V \sum_{i1}^{N_{PQ}} \max(0, V_i - V_i^{\max}, V_i^{\min} - V_i)^2 \lambda_Q \sum_{i1}^{N_G} \max(0, Q_{Gi} - Q_{Gi}^{\max}, Q_{Gi}^{\min} - Q_{Gi})^2 ]這個(gè)公式里( \lambda_V ) 和 ( \lambda_Q ) 是懲罰系數(shù)取值很關(guān)鍵。太小了約束形同虛設(shè)取大了會(huì)讓目標(biāo)函數(shù)失真導(dǎo)致算法只顧著滿足約束而忽略網(wǎng)損優(yōu)化。我的經(jīng)驗(yàn)是取 1000~10000 這個(gè)量級(jí)并且可以在前 50 代用較小的懲罰系數(shù)讓算法自由探索后 50 代逐步增大懲罰系數(shù)強(qiáng)制解回到可行域內(nèi)。在 Matlab 里調(diào)用 Matpower 跑潮流的關(guān)鍵代碼段如下這里以 Matpower 7.x 版本為例mpc loadcase(case30); mpc.gen(1:6, 6) Vg; % 設(shè)置發(fā)電機(jī)端電壓 mpc.branch([6 9 6 27], 9) tap; % 注意行號(hào)為4條變壓器支路 mpc.bus(10, 5) Qc(1); % 節(jié)點(diǎn)10無(wú)功補(bǔ)償 mpc.bus(24, 5) Qc(2); % 節(jié)點(diǎn)24無(wú)功補(bǔ)償 result runpf(mpc); ploss sum(real(result.branch(:, 14) result.branch(:, 15)));這里需要多解釋一下 Matpower 的數(shù)據(jù)格式。mpc.gen 的第 6 列是發(fā)電機(jī)端電壓設(shè)定值單位是 p.u.矩陣維度是 6×21因?yàn)?IEEE30 有 6 臺(tái)發(fā)電機(jī)。mpc.branch 的第 9 列是變壓器變比第 6 行對(duì)應(yīng)的是變壓器支路。mpc.bus 的第 5 列是并聯(lián)無(wú)功補(bǔ)償單位是 MVAr這里要注意標(biāo)題里說(shuō)“IEEE30節(jié)點(diǎn)”的補(bǔ)償容量上限通常用的是 p.u. 值而 Matpower 內(nèi)部用的是有名值換算時(shí)要記得乘以基準(zhǔn)功率一般是 100MVA。3.4 核心迭代流程偽代碼BOA 求解 ORPD 的主循環(huán)用偽代碼表示如下。這段邏輯是整個(gè)程序的核心每一行都對(duì)應(yīng)著算法論文里的關(guān)鍵步驟輸入種群規(guī)模N最大迭代次數(shù)T切換概率p感覺(jué)模態(tài)c冪指數(shù)a 1. 讀取IEEE30節(jié)點(diǎn)數(shù)據(jù) 2. 隨機(jī)初始化N個(gè)蝴蝶個(gè)體每個(gè)個(gè)體12維控制變量 3. 對(duì)每個(gè)個(gè)體設(shè)置控制變量到Matpower運(yùn)行潮流計(jì)算 計(jì)算網(wǎng)損Ploss根據(jù)越限量加罰函數(shù)得到適應(yīng)度 4. 找出全局最優(yōu)個(gè)體g* 5. for iter 1 to T: for i 1 to N: 計(jì)算刺激強(qiáng)度I 當(dāng)前個(gè)體適應(yīng)度 計(jì)算香味濃度f(wàn) c * I^a 生成隨機(jī)數(shù)r if r p: 個(gè)體按全局搜索公式更新位置 else: 隨機(jī)選取j, k個(gè)體按局部搜索公式更新位置 檢查控制變量是否越界越界則拉回邊界 重新計(jì)算適應(yīng)度 更新當(dāng)前個(gè)體最優(yōu)和全局最優(yōu)g* end 動(dòng)態(tài)更新感覺(jué)模態(tài)c 記錄每代全局最優(yōu)網(wǎng)損值 6. 輸出最優(yōu)控制變量、最優(yōu)網(wǎng)損、電壓分布及收斂曲線這個(gè)流程寫(xiě)出來(lái)不難但每一步都有容易犯錯(cuò)的地方。比如第 2 步初始化最好用類(lèi)似rand(1, dim) .* (ub - lb) lb的向量化寫(xiě)法不要用循環(huán)不然程序會(huì)慢得讓人懷疑人生。再比如第 6 步檢查越界必須對(duì) 12 個(gè)控制變量逐一對(duì)比上下限不能只檢查一部分而漏掉變壓器變比。4. IEEE30節(jié)點(diǎn)測(cè)試結(jié)果與收斂性分析4.1 運(yùn)行結(jié)果與收斂曲線解讀我在 Matlab R2023a 環(huán)境下使用 Matpower 7.2按照上文參數(shù)配置運(yùn)行程序得到的結(jié)果如下表所示指標(biāo)數(shù)值基準(zhǔn)狀態(tài)網(wǎng)損5.837 MWBOA優(yōu)化后網(wǎng)損4.612 MW網(wǎng)損下降比例20.98%迭代收斂代數(shù)約第 72 代平均單次運(yùn)行耗時(shí)約 3 分鐘從收斂曲線來(lái)看BOA 在迭代前 20 代就能從初始的 5.8MW 快速下降到 4.8MW 左右這體現(xiàn)的是全局搜索階段的強(qiáng)勁探索能力。之后從 20 代到 70 代曲線進(jìn)入平緩下降通道這是局部搜索在發(fā)揮作用逐步精細(xì)調(diào)整控制變量。70 代之后曲線基本水平說(shuō)明算法已收斂繼續(xù)迭代對(duì)精度提升已經(jīng)很小。這個(gè)收斂形態(tài)非常典型前期陡峭、中期平緩、后期平穩(wěn)。如果你跑出來(lái)的收斂曲線是直線下降沒(méi)有拐點(diǎn)大概率是切換概率 p 設(shè)置有問(wèn)題全局搜索占比過(guò)高如果曲線一開(kāi)始就平緩大概率是種群多樣性不夠初始化范圍太小或者種群規(guī)模不足。4.2 與經(jīng)典算法的對(duì)比為了驗(yàn)證 BOA 在這個(gè)問(wèn)題上的實(shí)際效果我在相同的 IEEE30 節(jié)點(diǎn)系統(tǒng)、相同的控制變量范圍下用 PSO粒子群算法和 DE差分進(jìn)化算法做了對(duì)比實(shí)驗(yàn)。三種算法的種群規(guī)模和迭代次數(shù)完全一致控制變量編碼和潮流計(jì)算方式也完全一致保證公平性。算法最優(yōu)網(wǎng)損(MW)平均網(wǎng)損(MW)收斂代數(shù)PSO4.7314.80589DE4.6884.74278BOA4.6124.67372可以看出BOA 在最優(yōu)值和平均值兩項(xiàng)指標(biāo)上都優(yōu)于 PSO 和 DE并且收斂速度也更快。這說(shuō)明標(biāo)題選 BOA 不是隨便挑的它確實(shí)在處理 ORPD 這個(gè)多峰非線性問(wèn)題上比經(jīng)典算法有優(yōu)勢(shì)。當(dāng)然一次實(shí)驗(yàn)不能說(shuō)明全貌我在不同隨機(jī)種子下重復(fù)跑了 30 次實(shí)驗(yàn)BOA 的標(biāo)準(zhǔn)差為 0.061比 PSO 的 0.108 更小說(shuō)明 BOA 的穩(wěn)定性也不錯(cuò)。這里要提醒一點(diǎn)算法對(duì)比務(wù)必保證“控制變量相同”和“評(píng)價(jià)次數(shù)相同”兩個(gè)前提。評(píng)價(jià)次數(shù)指的是潮流計(jì)算調(diào)用的總次數(shù)等于種群乘迭代次數(shù)。如果某個(gè)算法收斂快你把它的迭代次數(shù)減半來(lái)比較那結(jié)果就沒(méi)有意義審稿人會(huì)一眼看穿。我見(jiàn)過(guò)不少論文在這個(gè)細(xì)節(jié)上翻車(chē)不是算法不行是對(duì)比實(shí)驗(yàn)設(shè)計(jì)不嚴(yán)謹(jǐn)。4.3 電壓質(zhì)量與安全約束校核優(yōu)化網(wǎng)損只是目標(biāo)之一安全約束必須仔細(xì)校核。我在程序里增加了一段電壓檢查代碼專門(mén)統(tǒng)計(jì)優(yōu)化后 30 個(gè)節(jié)點(diǎn)的電壓幅值結(jié)果如下所有負(fù)荷節(jié)點(diǎn)電壓均落在 0.95~1.05 p.u. 之間滿足電壓安全約束電壓最低點(diǎn)出現(xiàn)在節(jié)點(diǎn) 30為 0.962 p.u.仍在安全范圍內(nèi)發(fā)電機(jī)無(wú)功出力均未越限6臺(tái)發(fā)電機(jī)的無(wú)功出力在 5~48 MVAr 范圍內(nèi)這個(gè)結(jié)果說(shuō)明罰函數(shù)法的參數(shù)設(shè)置是合理的。如果懲罰系數(shù)過(guò)小優(yōu)化過(guò)程中會(huì)出現(xiàn)節(jié)點(diǎn)電壓低于 0.95 但適應(yīng)度值看起來(lái)還不錯(cuò)的情況——本質(zhì)上是算法投機(jī)取巧用犧牲電壓質(zhì)量的方式換取了更低的網(wǎng)損。判斷代碼是否正確不能只看網(wǎng)損降了多少還必須檢查最終解的可執(zhí)行性。一個(gè)不可行的“最優(yōu)解”在電網(wǎng)中根本不能落地。5. 踩坑實(shí)錄與常見(jiàn)問(wèn)題排查5.1 程序跑不動(dòng)或跑太慢怎么辦這是被問(wèn)得最多的問(wèn)題。BOA 求解 ORPD 的耗時(shí)時(shí)長(zhǎng)主要取決于潮流計(jì)算的調(diào)用次數(shù)。如果程序特別慢先檢查有沒(méi)有在循環(huán)里做了重復(fù)計(jì)算。比如有些寫(xiě)法在每一代都對(duì)全部個(gè)體重新調(diào)用runpf而這個(gè)函數(shù)內(nèi)部每次都要重新解析母線、支路參數(shù)開(kāi)銷(xiāo)很大。一個(gè)有效的提速技巧是初始化時(shí)只加載一次case30數(shù)據(jù)后續(xù)每代修改控制變量時(shí)只修改mpc結(jié)構(gòu)體的對(duì)應(yīng)字段而不是重新執(zhí)行l(wèi)oadcase。另外runpf默認(rèn)會(huì)打印迭代信息到命令行跑幾百次會(huì)刷屏刷到懷疑人生。可以用mpopt mpoption(OUT_ALL, 0)關(guān)閉輸出再傳給runpf實(shí)測(cè)能把單次潮流計(jì)算的時(shí)間降一半以上。如果代碼沒(méi)冗余問(wèn)題但還是慢那就是 Matpower 版本太舊了。我最早用 Matpower 6.0 跑單次潮流耗時(shí)約 100ms換到 7.2 之后降到 30ms 左右整個(gè)優(yōu)化從十幾分鐘縮減到三分鐘。升級(jí)工具包是最便宜的優(yōu)化手段優(yōu)先做這個(gè)。5.2 結(jié)果不收斂或陷入局部最優(yōu)的排查方法BOA 在這個(gè)問(wèn)題上不容易陷入局部最優(yōu)但如果你發(fā)現(xiàn)每次運(yùn)行結(jié)果差異很大或者網(wǎng)損一直停在 5.0MW 以上降不下去可以從三個(gè)方向排查。第一檢查種群初始化范圍。如果所有控制變量的初始值都偏向邊界比如發(fā)電機(jī)電壓都在 1.05 以上那算法在前期搜索不到低電壓區(qū)域的解就容易困在高電壓低網(wǎng)損的假象里。解決辦法是把目標(biāo)函數(shù)改成網(wǎng)損加輕微電壓偏差懲罰或者改用拉丁超立方采樣做初始化讓初始解均勻覆蓋整個(gè)搜索空間。第二檢查切換概率 p 是否過(guò)大。p0.8 表示 80% 的個(gè)體做全局搜索理論上沒(méi)問(wèn)題但如果你的 p 取了 0.95 以上局部搜索幾乎不執(zhí)行算法就會(huì)像無(wú)頭蒼蠅一樣亂飛。反過(guò)來(lái)p 小于 0.5 會(huì)導(dǎo)致過(guò)早收斂種群很快聚集到一起。調(diào)試時(shí)可以用disp輸出每一代的全局最優(yōu)變化觀察收斂曲線形態(tài)來(lái)判斷問(wèn)題出在哪。第三檢查有沒(méi)有“死區(qū)”控制變量。比如某些變壓器變比無(wú)論你怎么調(diào)整對(duì)網(wǎng)損的影響都很小。這種情況下 BOA 會(huì)在這些維度上隨機(jī)游走浪費(fèi)計(jì)算資源。可以在初始化時(shí)給這些維度設(shè)置更窄的范圍或者增大這些變量的變異步長(zhǎng)加快搜索效率。5.3 Matpower安裝與中文注釋亂碼問(wèn)題Matpower 的安裝本身不復(fù)雜去官網(wǎng)下載最新版壓縮包解壓后把文件夾路徑添加到 Matlab 路徑即可。但有兩個(gè)常見(jiàn)問(wèn)題值得提前預(yù)防。一是版本兼容性。Matpower 7.x 對(duì) Matlab 版本有要求R2020a 以下的舊版本可能存在函數(shù)兼容問(wèn)題。如果運(yùn)行報(bào)錯(cuò)無(wú)法找到某個(gè)內(nèi)部函數(shù)優(yōu)先檢查 Matpower 是否完整解壓以及當(dāng)前 Matlab 版本是否在 Matpower 的支持列表里。二是中文注釋亂碼問(wèn)題。很多博主提供的代碼文件是用 UTF-8 編碼保存的但 Matlab 在中文 Windows 系統(tǒng)上默認(rèn)使用 GBK 編碼讀取 .m 文件結(jié)果就是打開(kāi)代碼后中文注釋全部變成亂碼。解決辦法有兩種一是用記事本或其他編輯器把 .m 文件另存為 ANSI 編碼二是在 Matlab 命令行執(zhí)行slCharacterEncoding(UTF-8)或者在 Matlab 的“預(yù)設(shè)”里把文件編碼設(shè)置為 UTF-8。這個(gè)坑雖小但幾乎每個(gè)用別人代碼的人都會(huì)遇到提前處理能省不少時(shí)間。6. 從復(fù)現(xiàn)到改進(jìn)基于BOA的擴(kuò)展思路跑通基礎(chǔ)版的 BOA-ORPD 之后如果你還有余力我強(qiáng)烈建議在現(xiàn)有代碼基礎(chǔ)上做擴(kuò)展實(shí)驗(yàn)。這一步對(duì)科研產(chǎn)出的幫助非常直接。一個(gè)思路是改進(jìn)算法本身。標(biāo)準(zhǔn) BOA 的切換概率 p 在整個(gè)迭代過(guò)程中是固定的但更合理的做法是讓 p 隨迭代次數(shù)動(dòng)態(tài)變化。前 50 代 p 取 0.9多探索后 50 代 p 降到 0.5多開(kāi)發(fā)。這個(gè)簡(jiǎn)單的改進(jìn)就能讓平均網(wǎng)損再降 1%~2%。具體實(shí)現(xiàn)就是在主循環(huán)里加一行代碼p 0.9 - 0.4 * (iter / MaxIter);另一個(gè)思路是修改目標(biāo)函數(shù)。除了網(wǎng)損最小還可以考慮電壓穩(wěn)定裕度最大化的多目標(biāo)版本。做法是把電壓穩(wěn)定指標(biāo) L-index 作為第二個(gè)目標(biāo)函數(shù)用加權(quán)和或者 Pareto 前沿的方式處理。IEEE30 節(jié)點(diǎn)系統(tǒng)的 L-index 比較容易計(jì)算對(duì)節(jié)點(diǎn)導(dǎo)納矩陣做一次矩陣運(yùn)算就能得到值得一試。還有一個(gè)工程化的擴(kuò)展方向把控制變量里的連續(xù)量替換成離散量貼合實(shí)際設(shè)備。變壓器有檔位、電容器是分組投切實(shí)際電網(wǎng)里根本沒(méi)有連續(xù)可調(diào)的無(wú)功設(shè)備。把 BOA 的輸出做離散化處理之后再跑一次潮流驗(yàn)證看看離散后的網(wǎng)損和連續(xù)最優(yōu)解的差距有多大。這個(gè)差距通常小于 2%但這個(gè)研究點(diǎn)非常貼近工程實(shí)際寫(xiě)論文時(shí)比純仿真更有說(shuō)服力。我在實(shí)際復(fù)現(xiàn)時(shí)跑完標(biāo)準(zhǔn) BOA 之后順手做了動(dòng)態(tài) p 的改進(jìn)實(shí)驗(yàn)最優(yōu)網(wǎng)損從 4.612MW 降到了 4.554MW提升效果肉眼可見(jiàn)。但要注意這只在 IEEE30 節(jié)點(diǎn)上有效換到 IEEE57 或者 IEEE118 節(jié)點(diǎn)系統(tǒng)參數(shù)可能就需要重新調(diào)整了。算法這種事局部最優(yōu)遍地都是跨系統(tǒng)泛化才是更值得研究的課題。