唐建林,余 濤,張孝順,李卓環(huán),陳俊斌
(1.華南理工大學電力學院,廣東廣州 510640;2.汕頭大學工學院,廣東汕頭 515063)
隨著化石燃料的消耗以及日趨嚴重的全球變暖形勢,可再生能源和高效能源實用技術正引起廣泛關注[1].為了解決當前能源短缺、環(huán)境惡化問題,世界眾多國家都開始調(diào)整能源布局,需找新的能源發(fā)展方式、提高能源利用率.
在轉(zhuǎn)變能源發(fā)展方式、調(diào)整能源結構上,人們開始更多的關注以風力發(fā)電和光伏發(fā)電為代表的可再生能源發(fā)電方式[2].這種分布式發(fā)電方式具有投資小、占地少、節(jié)能環(huán)保等優(yōu)點[3].由于風力發(fā)電和光伏發(fā)電時輸出功率的隨機性和波動性,傳統(tǒng)的經(jīng)濟調(diào)度方法已經(jīng)不再適用[4].以往的多數(shù)研究中沒有考慮風力發(fā)電、光伏發(fā)電產(chǎn)生的成本,在分析計算中均將風電、光電等效成只發(fā)電卻不產(chǎn)生成本的發(fā)電機[5].這種處理方式并沒有考慮風電、光電的隨機性、波動性對電力系統(tǒng)裝機容量提出的額外要求,這勢必會對經(jīng)濟調(diào)度分析結果的客觀性產(chǎn)生影響.
電力系統(tǒng)和熱力系統(tǒng)具有天然的互補性,對二者進行統(tǒng)一調(diào)度能夠在更廣的時空范圍內(nèi)提高資源優(yōu)化配置的能力和能量的利用效率.我國大型風電場、光伏電場集中的北方地區(qū),現(xiàn)有常規(guī)機組中50%~70%的發(fā)電機組為熱電聯(lián)產(chǎn)[6].高滲透率風電、光電與高比例熱電聯(lián)產(chǎn)機組占主導的電源結構,給電力系統(tǒng)調(diào)度和控制帶來了巨大挑戰(zhàn).
為了實現(xiàn)風力機組、光伏機組和熱電聯(lián)產(chǎn)機組的綜合利用,對熱電經(jīng)濟調(diào)度(combined heat and power economic dispatch,CHPED)進行優(yōu)化,在滿足系統(tǒng)所有不等式和等式約束的同時,實現(xiàn)了最優(yōu)的發(fā)電調(diào)度和供熱調(diào)度,使供電和供熱的總成本降至最低.處理該問題的方法主要有兩大類,第一類是以二次規(guī)劃法[7]、內(nèi)點法[8]、牛頓法[9]為代表的經(jīng)典優(yōu)化算法.它們在求解經(jīng)濟調(diào)度問題依賴具體的數(shù)學模型且如果系統(tǒng)是非線性的,均在不連續(xù)函數(shù)和約束時,容易得到局部最優(yōu)解.雖然應用Cplex等成熟數(shù)學求解器可以實現(xiàn)對經(jīng)濟調(diào)度的快速求解,但是其模型應用設計比較繁瑣.第二類是啟發(fā)式人工智能算法.遺傳算法(genetic algorithm,GA)[10]、粒子群算法(particle swarm optimization,PSO)[11]、差分進化算法(differential evolution algorithm,DE)[12]、灰狼優(yōu)化算法(grey wolf optimizer,GWO)[13]等智能算法對具體的數(shù)學模型要求程度低,應用更簡便,在電力系統(tǒng)各類優(yōu)化問題上已有較為成熟的應用.然而大多數(shù)算法都有其各自的弱點,主要特征如下:
1)種群的初始化完全隨機,無法提高解的多樣性,從而無法避免減少過早收斂問題.
2)求解優(yōu)化問題效率低下,大部分智能算法采用單一的群體搜索一個潛在的全局最優(yōu)解,很容易陷入低質(zhì)量的局部最優(yōu)解.考慮到發(fā)電機的閥點效應和熱電聯(lián)產(chǎn)系統(tǒng)的非凸不等式約束,CHPED問題是一個高度非線性、非光滑、非凸的問題,需要強大的優(yōu)化方法來求解,以避免陷入局部最優(yōu)解.本文提出一種新的基于混沌映射理論的優(yōu)化算法.該算法因其內(nèi)部具有多層次的搜索器結構,被稱為“多搜索器優(yōu)化算法”.根據(jù)3種不同類型的測試函數(shù)和27機組熱電系統(tǒng)仿真得到的結果,表明多搜索器優(yōu)化算法在非線性、非光滑、非凸的優(yōu)化問題上是有用的,并且效果顯著.
這一部分主要介紹風速和太陽光照的模擬情況與功率輸出的關系.
2.1.1 風速和功率輸出的關系
風速的概率分布參數(shù)是最能體現(xiàn)當?shù)仫L資源統(tǒng)計特性的指標,同時也是風力發(fā)電研究中的重要參數(shù).有很多種模型可用來近似擬合風速概率分布,主要包括威布爾分布、瑞利分布、對數(shù)正態(tài)分布和正態(tài)分布等多種函數(shù).由于兩參數(shù)威布爾模型可調(diào)整參數(shù)來適應風速的周期變化,所以在風速的概率建模分中應用非常廣泛.服從威布爾分布的概率密度函數(shù)定義如下[14]:
相應的累積分布函數(shù)為
式中:V為風速隨機變量;v為風速;k為風速概率分布函數(shù)的形狀參數(shù),它決定了風速概率分布密度曲線的總體形狀;c為尺度參數(shù),反應風場的平均風速,起放大或縮小曲線的作用,但不影響分布密度曲線的形狀.
假設風力發(fā)電機的輸出特性是一個簡化的分段線性函數(shù),對于給定的風速v,其輸出功率定義如下[15]:
式中:pwt為單臺風力機組的輸出功率;為額定功率;vin為風力機組的切入風速;vr為風力機組的額定風速;vout為風力機組的切出風速.
通過如下的線性變換,將風速的概率分布轉(zhuǎn)變?yōu)轱L電的概率分布[16]:
式中:T為一般情況的變化函數(shù)形式;PWT為風電隨機變量;V為風速隨機變量.
對于威布爾函數(shù)來說,風電輸出功率隨機變量的離散部分如下:
當風機的輸出功率處于連續(xù)變化的范圍內(nèi)時,其威布爾的概率密度函數(shù)如下式:
式中φ=(vr?vin)/vin.
2.1.2 太陽光照和功率輸出的關系
根據(jù)統(tǒng)計數(shù)據(jù),太陽的光照強度在一定時間段內(nèi)可以近視地認為服從Beta分布,它的概率密度函數(shù)如下所示[17]:
式中:rmax為最大光照強度;r為實際光照強度;μ為數(shù)學期望;σ為標準差.
光伏系統(tǒng)的輸出功率與光照強度近視成線性關系,不考慮功率損耗,且忽略環(huán)境溫度的影響,則可以得到下式[18]:
式中:A為光伏電池陣列總面積;η為光伏組件總的光電轉(zhuǎn)化效率.從而,可以推導出光伏電池的概率密度函數(shù)如下:
式中:ppv為光伏系統(tǒng)的輸出功率;為光伏系統(tǒng)的最大輸出系統(tǒng).
2.2.1 目標函數(shù)
整個系統(tǒng)的架構如圖1所示,供給側(cè)主要由燃氣發(fā)電廠、熱電聯(lián)產(chǎn)裝置、供熱廠、風電場和光伏電廠所構成,它們產(chǎn)生的電和熱分別通過線路和管道進行傳輸?shù)竭_需求側(cè).
圖1 含可再生能源的系統(tǒng)框架Fig.1 System framework with renewable energy
其目標是使整個系統(tǒng)的總生產(chǎn)成本最小.
2.2.2 僅供電成本
通常情況下,燃料成本可以用二次函數(shù)來表示,但是在實際操作中,大型發(fā)電機組通常含有多氣門汽輪機,當蒸汽進氣閥開始啟動時,由于拉絲效應,燃油損耗將會急劇增加從而導致燃油成本的增加,會在機組耗量曲線上疊加一個脈動效果,產(chǎn)生了閥點效應.為了更加精確地表示發(fā)電機組的發(fā)電費用,實際上需要在原來二次函數(shù)的基礎上,再疊加一個正弦函數(shù),這將增加非光滑和非凸的特性,其耗量特性曲線方程如圖2所示[19].
圖2 閥點效應Fig.2 Valve-point loading effect
2.2.3 熱電聯(lián)產(chǎn)成本
熱電聯(lián)產(chǎn)機組不僅能實現(xiàn)熱與電負荷的供應,還能減少污染氣體排放,具有良好的社會效益和經(jīng)濟效益,其成本函數(shù)如下[20]:
2.2.4 僅供熱成本
產(chǎn)熱機組的成本函數(shù)如下式[21]:
式中:ak,bk,ck為成本系數(shù);為產(chǎn)熱量;Nk為僅供熱機組的數(shù)量.
2.2.5 風力發(fā)電成本
通常情況下,風力發(fā)電廠由私人所有,電網(wǎng)運營商需要根據(jù)電力購買協(xié)議從風力發(fā)電商購買電力,需要支付一定的費用[22].其費用包含3個部分.第1部分是購買風力的直接成本.
式中dwt,l為風電的直接成本系數(shù).
第2部分是低估了風電出力而導致的懲罰成本.由于風力發(fā)電的波動性很大,風電預測通常會有一些誤差.當計劃的風力發(fā)電量低于可利用的風力發(fā)電量時,多余的風能就會被浪費.
式中:Kue,wt,l分別為風電的低估系數(shù);Pwt,rate,l為第l臺風力機組的額定功率;Pwt,l為第l臺風力機組的計劃輸出功率.
第3部分是高估風電出力而產(chǎn)生的輔助服務成本.如果按照計劃的風力發(fā)電高于實際可用的風力發(fā)電,那么就應該通過其他地方獲得備用電源來滿足負荷需求.
式中Koe,wt,l為風電的高估成本.
總的風力發(fā)電成本如下:
2.2.6 光伏發(fā)電成本
光伏發(fā)電的成本函數(shù)類似于風力發(fā)電,也是由3部分構成[22].
式中:dpv,m為光電的直接成本系數(shù);Kue,pv,m為光電的低估系數(shù);Koe,pv,m為光電的高估系數(shù).
總的光伏發(fā)電成本如下:
2.3.1 等式約束
燃油機組、熱電聯(lián)產(chǎn)機組、風力機組和光伏發(fā)電機組的電力平衡約束由等式約束來表示,熱電聯(lián)產(chǎn)機組和產(chǎn)熱機組的熱平衡約束也由等式約束來表示.
式中:Pd為電負荷需求功率;Hd為熱負荷需求功率.
2.3.2 不等式約束
1)純發(fā)電機組的運行約束為式(28)、熱電聯(lián)產(chǎn)機組中關于電和熱的運行約束為式(32)和式(33)、風力發(fā)電的運行約束為式(29)、光伏發(fā)電的運行約束為式(30)、純發(fā)熱機組的運行約束為式(31):
2)對于熱電聯(lián)產(chǎn)機組,圖3表示了熱電聯(lián)產(chǎn)機組的熱電耦合關系,是由定點ABCDEF連線所圍成的封閉區(qū)域,沿著區(qū)域的邊界線段BC,機組發(fā)熱量遞增,出力遞減,而沿著線段CD機組發(fā)熱量則遞減[23].
圖3 熱電聯(lián)產(chǎn)裝置運行區(qū)域Fig.3 Feasible operation region of the CHP units
本節(jié)主要介紹多搜索器優(yōu)化算法.大多數(shù)啟發(fā)式算法的靈感來自大自然,人們從大自然的運行規(guī)律中找到了許多解決實際問題的方法.但是多跟蹤器優(yōu)化算法并不是受自然的啟發(fā),而是根據(jù)其他算法的優(yōu)點和缺點設計出來的,并且整合了混沌理論.
在許多實際的系統(tǒng)中,混沌是一種無處不在的非線性現(xiàn)象.混沌運動可以根據(jù)自身的規(guī)律性和遍歷性,達到一定規(guī)模的每一個狀態(tài),這比簡單的隨機搜索算法要好.混沌搜索具有隨機性、遍歷性和規(guī)律性.因此,在許多優(yōu)化算法中,通常利用混沌理論引入初始化解或替換算法中的隨機參數(shù)以增強遍歷性并加速全局最優(yōu)收斂.由于混沌的非重復特性,基于混沌映射的算法比原來的算法更具有更全面的搜索能力.該特性是算法從局部最優(yōu)中逃逸的關鍵,避免了過早收斂[24].
多搜索器優(yōu)化算法在初始化的過程中使用了Logistic映射公式來產(chǎn)生具有混沌狀態(tài)的混沌變量,混沌迭代公式如下式所示:
式中:0 ≤Z0≤1,Zt個變量;κ為控制變量,當選κ=4時,系統(tǒng)處于完全混沌狀態(tài),其混沌空間為[0,1].
隨機產(chǎn)生一個n維且各分量值均在0~1之間的混沌矢量Zl=(zl1,zl2,···,zln)作為第1個種群.以Z1為初始值按照式(34)進行N ?1迭代搜索計算,就得到N個矢量Z1,Z2,···,ZN.將產(chǎn)生的混沌矢量Zi(i=1,2,···,N)的各分量按式(35)映射到式(28)–(31)的約束范圍,就可以得到各機組的發(fā)電功率和發(fā)熱量的初始化狀態(tài),完成種群的初始化過程:
式中xmax和xmin為優(yōu)化變量的上、下限.
在搜索區(qū)域內(nèi)隨機分布一系列搜索器,被稱為“全局搜索器”,并命名為GS.圍繞在每個全局搜索器附近,并在指定半徑內(nèi)以此搜索器為中心,隨機分布了另外一些搜索器,被稱為“局部搜索器”,命名為LS,這些局部搜索器服從全局搜索器.設立局部搜索器的目的是為了在每個全局搜索器附近能夠感知周圍的動態(tài)變化過程然后進行局部搜索,當全局搜索器搜索完后,局部搜索器便展開后續(xù)搜索.這一特性能夠提高收斂速度,避免局部最優(yōu).
全局搜索器的搜索半徑的計算如式(36),它取決于跟全局最優(yōu)點的距離和每個全局搜索器的等級:
式中:rmin和rmax分別為預先設定的最小半徑和最大半徑;kGS定義為全局搜索器的等級,當全局搜索器優(yōu)化后的目標值最小時,kGS為0,以此類推,當全局搜索器得到的目標值最大時,kGS為最大;nGS為全局搜索器的個數(shù);Gbest定義為全局最優(yōu)解;dist(GS,Gbest)定義了全局搜索器跟全局最優(yōu)點之間的距離.
當局部搜索器的搜索結束后,為了全局搜索器能夠找到最優(yōu)點,全局搜索器開始進行隨機游走.隨機游走規(guī)則見圖4,已經(jīng)成功應用于很多算法里面[25].從一個點隨機游走到另一個點,需要考慮兩個變量:距離和方向角,其中距離用d來表示,方向角用Θ來表示.
其中w為π/4.
圖4 全局搜索器游走規(guī)則Fig.4 Global searcher walking rules
當所有的全局搜索器按照隨機游走規(guī)則移動就可以得到新的位置,所得到的解跟目前的全局最優(yōu)解Gbest進行比較,如果效果更好,全局最優(yōu)解將會被替換.
本文采用罰函數(shù)法把系統(tǒng)的約束條件引入到目標函數(shù)中,通過罰函數(shù)法處理不等式約束,并保留等式約束的可行解,以避免算法出現(xiàn)“越限”現(xiàn)象.將等式約束優(yōu)化轉(zhuǎn)換為無等式約束問題,采用一個懲罰系數(shù)M構建的適應度函數(shù)為
顯然,當滿足等式約束條件時,此時罰函數(shù)為0,適應度函數(shù)的值即為成本函數(shù).
本文所提算法的求解流程圖如圖5所示.
圖5 多搜索器優(yōu)化算法流程圖Fig.5 Flow chart of multi searcher algorithm
為了測試所提算法在電力系統(tǒng)經(jīng)濟調(diào)度模型中的求解性能,本文引入10種不依賴于具體數(shù)學模型的智能算法與多搜索器優(yōu)化算法進行求解效率的對比,分別是生物地理學優(yōu)化算法(biogeography-based optimization,BBO)[26]、文化算法(cultural algorithms,CA)[27]、螢火蟲算法(firefly algorithm,FA)[28]、遺傳算法(genetic algorithm,GA)[10]、灰狼優(yōu)化算法(grey wolf optimizer,GWO)[13]、飛蛾撲火算法(moth-flame optimization,MFO)[29]、粒子群優(yōu)化算法(particle swarm optimization,PSO)[11]、模擬退火算法(simulated annealing,SA)[30]、教與學優(yōu)化算法(teach-learn based optimization algorithm,TLBO)[31]和改進遺傳算法(proposed genetic algorithms,PGA)[32].
根據(jù)上述的仿真模型以及大量的測試結果,本文的所有的算法參數(shù)設置如表1 所示.仿真計算均在CPU為英特爾Core i7–6700、主頻3.4 GHz、內(nèi)存16 GB的計算機運行.
表1 算法主要參數(shù)設置Table 1 Parameter setting of algorithms
為了驗證算法的有效性,下面將選擇3種不同類型的函數(shù)進行數(shù)值實驗.迭代步數(shù)均設為100.
Sphere 函數(shù)是較為簡單的單峰函數(shù),可以考察算法的收斂速度:
全局最優(yōu)點為xi=0,f(x)=0.單峰函數(shù)中的適應度收斂曲線對比結果如圖6所示,可以看出MFO無法收斂到全局最優(yōu),MSO的收斂速度最快,PGA的收斂速度次之.
圖6 11個算法在單峰測試函數(shù)中的適應度函數(shù)收斂曲線對比Fig.6 Comparison of convergence curves of eleven algorithms in single peak test function
Step函數(shù)是個不連續(xù)的階梯函數(shù),可以用來檢驗算法的有效性:
全局最優(yōu)點xi=0,f(x)=0.不連續(xù)階梯函數(shù)的適應度收斂曲線對比結果如圖7所示,可以看出MSO和FA的收斂速度最快.
圖7 11種算法在不連續(xù)階梯測試函數(shù)中的適應度函數(shù)收斂曲線對比Fig.7 Contrast of convergence curves of fitness functions of eleven algorithms in discontinuous step test functions
Rastrigin 函數(shù)是多峰函數(shù),可以用來檢驗算法跳出局部最優(yōu)能力:
全局最優(yōu)點xi=0,f(x)=0.不連續(xù)階梯函數(shù)的適應度收斂曲線對比結果如圖8所示,可以看出MSO的全局尋優(yōu)能力體現(xiàn)較為明顯,其他對比算法都不同程度地陷入局部最優(yōu)點.
圖8 11種算法在多峰測試函數(shù)中的適應度函數(shù)收斂曲線對比Fig.8 Contrast of convergence curves of fitness functions of eleven algorithms in multimodal test functions
本系統(tǒng)含有27臺機組包含13臺燃油機組并且都考慮了閥點效應、4臺熱電聯(lián)產(chǎn)機組、4臺只供熱機組、3臺風電機組、3臺光伏機組.對于燃油機組、熱電聯(lián)產(chǎn)機組、供熱機組來說,其數(shù)據(jù)分別如表2–4,熱電聯(lián)產(chǎn)機組的可行運行區(qū)域如圖9所示[33].對于風電機組來說,風力發(fā)電機組的容量分別為130 MW,94 MW和94 MW;形狀參數(shù)k為2,尺度參數(shù)c為15;切入風速為15 m/s,額定風速為5 m/s,切出風速為45 m/s[23];對于光伏機組來說,光伏電站的光伏列陣總面積為80000 m2,光電轉(zhuǎn)換效率為14%,最大光照強度為700 W/m2;形狀參數(shù)為0.95,尺度參數(shù)為0.95[23].
風力機組直接成本系數(shù)為120$/MWh、低估成本系數(shù)為15$/MWh、高估成本系數(shù)為20$/MWh.光伏機組的直接成本系數(shù)為200$/MWh、低估成本系數(shù)為15$/MWh、高估成本系數(shù)為20$/MWh[23].總的電力需求、熱需求分別為3837.136 MW,615.372 MW[23].所有的網(wǎng)絡損耗均被忽略掉.為兼顧算法最佳性能、計算速度和公平性,迭代步數(shù)設置為300.
表2 僅發(fā)電系統(tǒng)參數(shù)Table 2 Parameters of thermal units
表3 熱電聯(lián)產(chǎn)系統(tǒng)參數(shù)Table 3 Parameters of CHP units
表4 僅產(chǎn)熱系統(tǒng)參數(shù)Table 4 Parameters of heat-only units
圖9 熱電聯(lián)產(chǎn)裝置運行區(qū)域Fig.9 Feasible operation region of the CHP unit
圖10代表了各種算法分別運行10次的盒須圖,從盒須圖對比中可以發(fā)現(xiàn)MSO收斂穩(wěn)定性不如GA和GWO,但是其變化范圍仍然是在可以接受的,且優(yōu)化效果要好于GA和GWO.
為了驗證MSO的速度,圖11給出各種算法分別運行10次得到的平均運算時間對比條形圖.從對比圖中可以看出,MSO算法的計算時間用時最少,比SA算法快了接近4.96倍,能夠很好地滿足時間上要求.
圖12、表5和表6給出了不同算法適應度函數(shù)的收斂情況和優(yōu)化結果對比,其中:Power代表機組類型為僅發(fā)電的燃氣機組;CHP代表熱電聯(lián)產(chǎn)機組;Heat代表僅發(fā)熱的熱機組;WT 代表風力發(fā)電機組;PV代表光伏發(fā)電機組.從最后優(yōu)化的結果可以看出,MSO得到的適應度函數(shù)值最小,從而可以驗證MSO算法能夠跳出局部最優(yōu)的缺點,具有很好的全局尋優(yōu)能力.并且由圖可得MSO從114代就開始收斂,極大提高了優(yōu)化效率.
圖10 11種算法分別運行10次盒須圖對比Fig.10 Box and whisker plots of total operating costs obtained by eleven algorithms in 10 runs
圖11 11種算法運行10次平均計算時間Fig.11 Average execution time obtained by eleven algorithms in 10 runs
本文以降低系統(tǒng)運行成本為目標建立起風–光–熱–電聯(lián)合調(diào)度的優(yōu)化模型,并提出了一種多搜索器優(yōu)化算法用來解決含可再生能源協(xié)同優(yōu)化調(diào)度問題,本文的主要貢獻如下:
1)引入混沌理論可有效保持初始種群的不重復性和搜索器的多樣性,改善傳統(tǒng)啟發(fā)式算法易陷入局部最優(yōu)的缺陷.
2)雙層搜索器的設置能夠有效地解決風光熱電協(xié)同調(diào)度問題,既能夠快速收斂到最優(yōu)點,又能避免陷入局部最優(yōu),提高優(yōu)化效果.
3)提出一種多搜索器算法,通過標準測試函數(shù)和算例的驗證,表明該算法能夠獲得更好的經(jīng)濟調(diào)度結果,該算法對優(yōu)化模型的依賴程度較低,能夠有效解決含發(fā)電機閥點效應的非凸優(yōu)化問題.
筆者下一步的研究工作將把遷移學習合理地整合到多搜索器優(yōu)化算法中,并且考慮網(wǎng)絡的實際損耗和風光的出力不確定性,完善風光熱電模型和應用到更加復雜的場景中.
表6 11種算法優(yōu)化結果對比(二)Table 6 Comparative results of optimal solutions obtained by different algorithms(part 2)