帶有精英保留機制的混合差分化學(xué)反應(yīng)算法

魏民，楊明磊，錢鋒

（華東理工大學(xué)化工過程先進控制與優(yōu)化技術(shù)教育部重點實驗室，上海200237）

引 言

智能優(yōu)化算法因其良好的全局搜索能力和通用性，一經(jīng)提出就受到廣泛關(guān)注。目前已經(jīng)有遺傳算法(GA)[1]、模擬退火(SA)[2]、粒子群算法(PSO)[3]和蟻群算法(ACO)[4]等一些經(jīng)典理論。近年來，對于經(jīng)典算法的改進和新型智能算法的開發(fā)一直是國內(nèi)外研究的熱點，化學(xué)反應(yīng)算法就是一種新型的智能優(yōu)化算法。

化學(xué)反應(yīng)優(yōu)化(chemical reaction optimization,CRO)算法是由Lam等[5]于2010年提出的一種基于種群的智能優(yōu)化算法，它以化學(xué)反應(yīng)過程中能量守恒為理論基礎(chǔ)，以尋找最低的能量消耗為目的，模擬反應(yīng)中物質(zhì)間的能量傳遞，最終在實際應(yīng)用中取得了良好的效果。

近些年關(guān)于化學(xué)反應(yīng)算法的研究引起了廣泛的關(guān)注，基本的CRO[5]用于求解組合優(yōu)化問題，Lam等[6]針對這一問題對化學(xué)反應(yīng)算法的收斂性進行了驗證。經(jīng)過改進的 CRO算法被廣泛應(yīng)用于求解各種實際問題。Sun等[7]對該算法進行了改進，使其與Lin-Kernighan局部搜索混合，加強其局部搜索能力，并把改進算法用于求解TSP問題。Lam等[8]提出了實數(shù)編碼的化學(xué)反應(yīng)算法(CRCRO)，使其適合求解連續(xù)函數(shù)的優(yōu)化問題。Pan等[9]將化學(xué)反應(yīng)算法用于優(yōu)化網(wǎng)絡(luò)節(jié)點,并取得了很好的效果。在智能計算方面，CRO還可用于進行神經(jīng)網(wǎng)絡(luò)訓(xùn)練[10]和模糊學(xué)習(xí)[11]。對于能源和環(huán)境等實際問題，CRO還可用來解決電力系統(tǒng)的最優(yōu)潮流計算[12]和空氣檢測器的傳感器分布[13]問題等。

雖然 CRO具有良好的全局搜索能力，但是這也導(dǎo)致了其在求解連續(xù)函數(shù)的優(yōu)化問題時收斂速度不高，求解精度下降[7]。因此，本文提出混合差分的化學(xué)反應(yīng)算法(HDECRO)，利用差分算法收斂速度快的特點來彌補 CRO在這方面的缺陷。對于在計算過程中，一些優(yōu)秀分子可能被反應(yīng)消耗掉的現(xiàn)象，本文還在分子群中加入了精英保留策略，以保證算法始終向著最優(yōu)的方向進化。

1 化學(xué)反應(yīng)算法

在化學(xué)反應(yīng)中，反應(yīng)過程往往是沿著最小的能量消耗路徑來進行的?；瘜W(xué)反應(yīng)算法就是根據(jù)這一自然現(xiàn)象與最優(yōu)化問題中尋求極值點的共同特點開發(fā)出來的?；瘜W(xué)反應(yīng)算法模仿化學(xué)反應(yīng)中分子所發(fā)生變化的情況，其目的在于捕捉到反應(yīng)過程中能量變化最小的那條反應(yīng)路徑。該算法包括兩個關(guān)鍵因素：分子和基本反應(yīng)。

1.1 分子

分子是執(zhí)行算法尋優(yōu)操作的個體，每個分子包括3個重要的組成部分：① 分子結(jié)構(gòu)；②潛在能量(potential energy, PE)；③動力能量(kinetic energy,KE)。各個組成部分的含義如下。

① 分子結(jié)構(gòu)：分子結(jié)構(gòu)用來表示每個分子所特有的原子組成和結(jié)構(gòu)，用ω表示。分子結(jié)構(gòu)的確定取決于目標(biāo)函數(shù)可行解的維度，即問題的操作變量有n個，那么相應(yīng)分子中所具有的原子就有n個；

② 潛在能量(PE)：PE表示當(dāng)前分子結(jié)構(gòu)ω所具有的目標(biāo)函數(shù)值，即PEω=f(ω)；

③ 動力能量(KE)：由于算法的評價機制可以歸納為 PEω+KE≥PEω，因此，KE表示當(dāng)前ω具有的跳出局部最優(yōu)，開發(fā)新的搜索區(qū)間的能力。

1.2 基本反應(yīng)

在化學(xué)反應(yīng)的過程中，一系列的碰撞會在分子之間發(fā)生，其碰撞不止在分子之間進行，也會在分子與容器壁間進行。因此，化學(xué)反應(yīng)算法模擬了 4種基本反應(yīng)過程，包括：碰壁反應(yīng)、分解反應(yīng)、分子碰撞和化合反應(yīng)。

① 碰壁反應(yīng)：單個分子碰撞容器并被彈回的過程。反應(yīng)導(dǎo)致分子結(jié)構(gòu)ω的輕微變化，屬于非劇烈反應(yīng)。如果當(dāng)前分子結(jié)構(gòu)為ω，則反應(yīng)產(chǎn)生的新分子結(jié)構(gòu)ω′=Neighbor(ω)一定在其附近。反應(yīng)發(fā)生的條件為

可以得到

其中，(PEω+KEω?PEω′)×(1?q)表示在碰壁反應(yīng)中分子所消耗掉的能量，q∈[KELossRate,1]，這部分能量儲存在buffer中

② 分解反應(yīng)：單個分子碰撞容器被彈回，并分解成為兩個或更多分子的過程。由于產(chǎn)生全新分子，反應(yīng)過程必然伴隨大量能量轉(zhuǎn)移，產(chǎn)生的新分子勢必會有與反應(yīng)前截然不同的分子結(jié)構(gòu)。如果當(dāng)前分子結(jié)構(gòu)為ω，新產(chǎn)生的分子結(jié)構(gòu)為ω′1和ω′2，那么反應(yīng)條件必須滿足

不難看出，想要滿足式(4)的反應(yīng)條件是非常困難的，因此過程中也允許buffer協(xié)助反應(yīng)的進行，即

buffer也隨之更新

③ 分子碰撞：兩個分子互相碰撞，然后各自被彈開的過程。反應(yīng)的劇烈程度和碰壁反應(yīng)相似，反應(yīng)結(jié)果只對各自的分子結(jié)構(gòu)有輕微的影響。假使原始分子結(jié)構(gòu)為ω1和ω2，通過反應(yīng)可以在各自鄰域產(chǎn)生新的分子結(jié)構(gòu)ω′1和ω′2。則反應(yīng)的條件為

如果條件公式不成立，則分子保持原有的屬性不變。

④ 化合反應(yīng)：多個分子發(fā)生碰撞并一起組成一個新的分子的過程。如果原始分子的結(jié)構(gòu)為ω1和ω2，兩者合成新的分子結(jié)構(gòu)ω′，由于化合反應(yīng)十分劇烈，所以ω′與反應(yīng)物分子結(jié)構(gòu)有很大不同。該反應(yīng)發(fā)生的條件為

可以得到

如果反應(yīng)條件不成立，則分子保持原來的屬性不變。由于反應(yīng)產(chǎn)物所帶有的 KEω′要比原分子的KEω1和KEω2大得多，所以化合反應(yīng)所得到的新分子ω′擁有更強的跳出局部最優(yōu)的能力。

2 混合差分化學(xué)反應(yīng)算法(HDECRO)

化學(xué)反應(yīng)算法被證明具有很強的全局搜索能力，但是由于算法中操作較多，導(dǎo)致其收斂速度較慢。而差分算法是一種經(jīng)典的隨機搜索算法，其結(jié)構(gòu)簡單，收斂速度快[10]。

2.1 混合差分變異策略

本文采用帶有三角變異的差分算法（TDE），該算法由Fan等[14]提出，相比于標(biāo)準DE算法，TDE具有更快的收斂速度和更高的求解精度，并已經(jīng)成功應(yīng)用到一些實際問題[15-16]，改進算法中TDE變異策略公式如下

其中

以概率CR進行交叉變異，最終得到新的分子結(jié)構(gòu)

2.2 精英保留策略

在化學(xué)反應(yīng)算法中，因為有KE的存在，導(dǎo)致原本結(jié)構(gòu)優(yōu)良的分子也可能在反應(yīng)過程中生成新的分子。這種現(xiàn)象一方面保持了種群的多樣性，使算法具有更好的全局搜索能力；但是另一方面也減緩了算法的收斂速度，在有限的計算次數(shù)時降低了計算的精確性。因此在改進算法中，引入了精英保留策略。

精英保留策略的具體實現(xiàn)方法是，在計算之初對所有分子所帶有的 PE進行升序排列，選擇最小的前10個分子存入精英數(shù)據(jù)庫，在這一代計算結(jié)束后更新精英分子群的數(shù)據(jù)庫，并把所得到的精英分子并入種群Pop中進行下一次迭代。

在這里需要注意的是，由于反應(yīng)前后能量守恒的原則，加入精英種群后需要在buffer中減去精英分子群增加的能量

2.3 改進算法概述

CRO算法總共由4種反應(yīng)組成：碰壁反應(yīng)、分解反應(yīng)、分子碰撞和化合反應(yīng)。這幾個反應(yīng)的根本不同點在于能量轉(zhuǎn)移的規(guī)模不同，這就使不同反應(yīng)所產(chǎn)生的新分子具有不同的搜索能力。

① 碰壁反應(yīng)能量轉(zhuǎn)移最少，新產(chǎn)生的分子和原分子的差異不大，可以利用這一特點進行鄰域搜索，提高碰壁反應(yīng)的概率能夠加強整個算法的局部搜索能力，提高搜索精度。在這里，利用TDE變異策略來獲得碰壁反應(yīng)中的新分子

② 分子碰撞也屬于弱反應(yīng)，用于鄰域搜索。因此在這里采用正態(tài)分布擾動來作為獲得新分子的變異策略，以此保持局部搜索中種群的多樣性，避免結(jié)果出現(xiàn)早熟現(xiàn)象

③ 分解反應(yīng)和化合反應(yīng)都屬于劇烈反應(yīng)，其中帶有大量的能量轉(zhuǎn)移，能夠產(chǎn)生截然不同的新分子，利用這一特點來保持整個算法的全局搜索能力。在這兩個反應(yīng)中，使用隨機的變異策略

其中，rand是函數(shù)約束范圍內(nèi)的隨機數(shù)。

可以看出 KE′decomposition數(shù)值較小，而 KE′synthesis數(shù)值較大，這也意味著兩種反應(yīng)所產(chǎn)生的新分子具有的穩(wěn)定性不同，即在接下來的計算中保存下來的能力不同，這也是同時使用這兩種反應(yīng)而不是有所取舍的原因[1]。

在這里需要注意的是，分解反應(yīng)和化合反應(yīng)都會導(dǎo)致整個分子種群的總數(shù)的變化，所以在實際應(yīng)用中應(yīng)該盡量控制這兩種反應(yīng)發(fā)生的概率。

算法的整個流程如圖1所示，在算法的初始化階段主要設(shè)定MoleColl（決定單分子反應(yīng)和多分子反應(yīng)進行的比例）、decomposition_criterion(單分子中分解反應(yīng)的比例)和 synthsis_criterion（多分子中化合反應(yīng)的比例）。由于分解和化合操作的存在，會造成整個種群數(shù)量的波動，因此在程序迭代過程中進行規(guī)定，如果種群數(shù)量大于1.5倍初始種群數(shù)量，則 decomposition_criterion=0（關(guān)閉分解反應(yīng)）；如果種群數(shù)量小于 0.8倍初始種群數(shù)量，則synthsis_criterion=0（關(guān)閉化合反應(yīng)）。

圖1 改進算法流程圖Fig. 1 Flow chart of modified algorithm

3 仿真測試與結(jié)果分析

3.1 測試函數(shù)與實驗環(huán)境

本文采用的測試函數(shù)如表1所示，其中D表示求解問題的維度，S表示函數(shù)的約束范圍，fmin表示函數(shù)的全局最小值。8個測試函數(shù)分別具有不同的特點，主要可以分為單峰測試函數(shù)、多峰測試函數(shù)和復(fù)合多峰函數(shù)3類。

表1 測試函數(shù)Table 1 Benchmark functions from CEC 2005

① 單峰測試函數(shù)：F1(x)～F4(x)，用來測試算法的收斂速度和求解精度；

② 多峰測試函數(shù)：F5(x)、F6(x)，用來測試算法全局搜索能力；

③ 復(fù)合多峰函數(shù)：F7(x)、F8(x)，這是一類結(jié)構(gòu)復(fù)雜又難于求解的測試函數(shù)。其中F7(x)的最優(yōu)點在邊界范圍內(nèi)，而F8(x)的全局最優(yōu)點在邊界上，圖形分別如圖2、圖3所示，這兩個函數(shù)用來進一步測試算法求解復(fù)雜多峰問題的全局搜索能力。

為了直觀比較改進算法的各項性能，本文選取了標(biāo)準差分算法（DE）[14,17]、改進差分算法(jDE、TDE)[18]、改進粒子群算法(LDW-PSO、PSOTVAC)[19-20]，以及實數(shù)編碼的化學(xué)反應(yīng)算法(RCCRO)[7]進行對比實驗。在計算過程中，每個算法均選用相同的進化次數(shù)和種群規(guī)模。

整個算法的實驗環(huán)境為 AMD Quad-Core 2.20 GHz, 4.00GB RAM。算法使用 Matlab編程，在MATLAB 7.11.0 平臺上運行。實驗過程中，對于每個測試函數(shù)，優(yōu)化算法均獨立求解30次，最終得到實驗結(jié)果。

圖2 測試函數(shù)F7(x)Fig.2 Benchmark function F7(x)

圖3 測試函數(shù)F8(x)Fig. 3 Benchmark function F8(x)

3.2 改進算法參數(shù)設(shè)定

圖4 F1(x)各參數(shù)調(diào)節(jié)值Fig. 4 Parameter value for F1(x)

改進算法中，需要預(yù)先設(shè)定的參數(shù)分別為InitialKE，MoleColl，KELossRate，decom_crite&syn_crite。由于分解反應(yīng)概率（decom_crite）和化合反應(yīng)概率（syn_crite）存在著相互制約的關(guān)系，前者使分子總量增多，后者使分子總量減少，因此這里兩個參數(shù)一同設(shè)定。

由于上述參數(shù)的選取對于改進算法的有效性起到關(guān)鍵作用，因此有必要對各項參數(shù)的功能進行介紹。

① InitialKE：決定分子群總能量的大小，即算法跳出局部最優(yōu)的能力，取值在區(qū)間[0,10000]內(nèi)，取值越小算法全局搜索能力越弱，取值過大會導(dǎo)致改進算法不收斂。

② MoleColl：判定分子進行單分子反應(yīng)或多分子反應(yīng)，在區(qū)間[0,1]內(nèi)取值，MoleColl的值越大代表算法進行多分子反應(yīng)的概率升高，即避免早熟的性能越好，但是由于算法主要依靠單分子反應(yīng)來進行收斂，因此應(yīng)該給予執(zhí)行單分子反應(yīng)更高的概率，建議取值在0.1～0.3之間。

③ KELossRate：KE存入buffer的比率，取值區(qū)間為[0,1]。KELossRate取值越大，分子群總KE的值減少越快，有助于改進算法的快速收斂，但是同時會導(dǎo)致算法全局搜索能力變差，一般取值 0.2左右。

④ decom_crite&syn_crite：判定算法執(zhí)行化合及分解反應(yīng)的概率，取值在區(qū)間[0,1]內(nèi)，為了保證算法的收斂性，該參數(shù)不宜過大，一般取0.05～0.2之間。

為了更準確地設(shè)定操作參數(shù)，本文分別選取單峰和多峰測試函數(shù)各1個來比較不同參數(shù)下改進算法的尋優(yōu)性能。測試中，被測算法的種群數(shù)量N=200，迭代次數(shù)G=1000。

圖4和圖 5分別為改進算法中不同參數(shù)對于F1(x)和F5(x)的測試結(jié)果。其中表現(xiàn)最好的一組數(shù)據(jù)所對應(yīng)的參數(shù)在圖中已經(jīng)用圓圈標(biāo)出。根據(jù)對結(jié)果的分析，可以確定算法的設(shè)定參數(shù)如表2所示。

表2 HDECRO參數(shù)設(shè)定Table 2 Parameter setting of HDECRO

3.3 實驗結(jié)果與分析

（1）F1(x)～F4(x)單峰測試函數(shù)的實驗結(jié)果如表 3所示，所有結(jié)果均由算法設(shè)定G=5000，并且獨立運行30次得到。表中，min表示30次計算中最小解，mean表示所有結(jié)果的平均值，stdDev表示所有結(jié)果的均方差，rank表示總結(jié)了這幾項指標(biāo)后對算法的排名。

圖5 F5(x)各參數(shù)調(diào)節(jié)Fig. 5 Parameter value for F5(x)

由表3可以看出，HDECRO在求解高維的單峰問題時，具有良好的計算精度，在一些問題的求解上，效果甚至優(yōu)于經(jīng)典的智能算法。

圖6是F1(x)～F4(x)測試函數(shù)優(yōu)化值隨進化代數(shù)的變化曲線，從圖中可以看出，由于KE的存在，HDECRO的進化曲線具有輕微振蕩的特點。

從對于單峰測試函數(shù)的仿真實驗中可以看出，未改進算法 RCCRO的收斂速度在所有算法中最慢，均排在最后一位，而混合了 TDE的改進算法HDECRO則很好地改善了收斂速度的缺陷，其求解精度與TDE相仿，均排在所有算法的前列。

表3 F1(x)～F4(x)單峰測試函數(shù)結(jié)果Table 3 Results of unimodal benchmark functions F1(x)—F4(x)

圖6 求解F1(x)～F4(x)函數(shù)收斂速度曲線Fig. 6 Convergence speed for benchmark functions F1(x)—F4(x)

（2）F5(x)和F6(x)多峰測試函數(shù)的測試結(jié)果如表4和圖7所示，數(shù)據(jù)由算法設(shè)定G=5000，并且獨立運行30次得到。從表中可以看出，HDECRO在處理高維的多峰問題時依然具有很高的求解精度，并能夠找到全局最優(yōu)解。

表4 F5(x),F6(x)多峰測試函數(shù)運行結(jié)果Table 4 Results of multi-modal benchmark functions F5(x), F6(x)

圖7 求解F5(x)，F(xiàn)6(x)函數(shù)收斂速度曲線Fig. 7 Convergence speed for benchmark functions F5(x), F6(x)

表5 F7(x),F8(x)復(fù)雜多峰函數(shù)運行結(jié)果Table 5 Results of composition multi-modal benchmark functions F7(x), F8(x)

需要注意的是，一些在單峰問題上表現(xiàn)突出的優(yōu)化算法在求解多峰問題上效果不佳，如TDE。但是結(jié)合了TDE的改進算法HDECRO依然具有較好的求解效果。這也驗證了改進算法在混合 TDE和 RCCRO時，能夠達到取長補短效果的必要性和有效性。

（3）F7(x)和F8(x)復(fù)雜多峰測試函數(shù)的實驗結(jié)果如表 5所示，結(jié)果均由算法設(shè)定G=1000，并且獨立運行30次得到。圖8和圖9分別表示對于F7(x)和F8(x)，每一次計算中優(yōu)化算法所能夠求得的最終優(yōu)化值。

可以很直觀地看出，經(jīng)典智能算法在處理復(fù)雜多峰問題時很容易陷入局部最優(yōu)，尤其是求解如F8(x)這樣全局最優(yōu)點在邊界上的函數(shù)。在處理帶有明顯欺騙性的復(fù)雜多峰問題時，傳統(tǒng)算法幾乎無法找到全局最優(yōu)，而 HDECRO得到目標(biāo)函數(shù)的全局最優(yōu)解的概率卻很高。從圖形中還可以看到，LDW-PSO同樣具有良好的全局收斂性，但是結(jié)合單峰問題的求解結(jié)果看，其求解精度不佳，總體來說效果不如本文提出的HDECRO算法。

圖8 F7(x)計算次數(shù)與最優(yōu)解曲線Fig. 8 Function value of each simulation for F7(x)

圖9 F8(x)計算次數(shù)與最優(yōu)解曲線Fig. 9 Function value of each simulation for F8(x)

4 結(jié) 論

本文提出了一種帶有精英保留機制的混合差分化學(xué)反應(yīng)優(yōu)化算法(HDECRO)，它利用了實數(shù)編碼的化學(xué)反應(yīng)優(yōu)化(RCCRO)算法優(yōu)良的全局搜索能力，并在此基礎(chǔ)上將碰壁反應(yīng)的變異策略改進為帶有三角變異的差分進化策略(TDE)，提高算法的求解精度。在每次迭代之初，還利用精英分子群合并的方法來保持進化過程中分子群的優(yōu)良品質(zhì)。通過對一系列具有不同特點的測試函數(shù)的仿真實驗可以證明，改進算法在相同的計算量的情況下具有良好的計算精確度，并且在處理復(fù)雜多峰問題時，依然具有良好的全局搜索能力。

[1]Holland J H. Adaptation in Natural and Artificial Systems[M].University Michigan Press, 1975

[2]Metropolis N, Rosenbluth A, Rosenbluth M, Teller A, Teller E.Equation of state calculations by fast computing machines [J].Chem.Phys., 1953, 21: 1087-1092

[3]Kennedy J, Eberhart R. Particle swarm optimization//Proceedings of the 4th IEEE International Conference on Neural Networks [C].Piscataway: IEEE Service Center, 1995:1942-1948

[4]Dorigo M，Gambardella L M. Ant colony system：a cooperative learning approach to the traveling salesman problem [J].IEEE Transactions on Evolutionary Computation, 1997(1): 53-66

[5]Lam A Y S, Li V O K. Chemical-reaction-inspired metaheuristic for optimization [J].IEEE Transactions on Evolutionary Computation,2010, 14(3): 381-399

[6]Lam A Y S, Li V O K, Xu J. On the convergence of chemical reaction optimization for combinatorial optimization [J].IEEE Transactions on Evolutionary Computation, 2013, 17(5):605-620

[7]Sun Jian, Wang Yuting. Hybrid algorithm based on chemical reaction optimization and Lin-Kernighan local search for the traveling salesman problem//IEEE 2011 Seventh International Conference on Natural Computation[C]. 2011

[8]Lam A Y S, Li V O K, Yu J J Q. Real-coded chemical reaction optimization [J].IEEE Transactions on Evolutionary Computation,2012, 16(3): 339-353

[9]Pan B, Lam A Y S, Li V O K. Network coding optimization based on chemical reaction optimization//IEEE Global Commun. Conf. [C].2011

[10]Yu J, Lam A Y S. Evolutionary artificial neural network based on chemical reaction optimization//IEEE Congress on Evolutionary Computation [C]. 2011

[11]Lam A Y S, Li V O K, Wei Z. Chemical reaction optimization for the fuzzy rule learning problem//WCCI 2012 IEEE World Congress on Computational Intelligence[C]. Australia, 2012

[12]Sun Y, Lam A Y S, Xu J. Chemical reaction optimization for the optimal power flow problem//WCCI 2012 IEEE World Congress on Computational Intelligence[C]. Australia, 2012

[13]Yu J, Li V O K, Lam A Y S. Sensor deployment for air pollution monitoring using public transportation system//WCCI 2012 IEEE World Congress on Computational Intelligence[C]. Australia, 2012

[14]Fan H Y, Lampinen J. A trigonometric mutation operation to differential evolution [J].Journal of Global Optimization, 2003, 27:105-129

[15]Angira R, Santosh A. Optimal control of chemical engineering processes using trigonometric differential evolution algorithm//Proceedings of International Symposium & 58th Annual Session of IIChE in association with International Partners[C].2005

[16]Angira R, Santosh A. Optimization of dynamic systems: a trigonometric differential evolution approach [J].Computers and Chemical Engineering, 2007, 31:1055-1063

[17]Price Kenneth V. Differential evolution: a fast and simple numerical optimizer//1996 Biennial Conference of the North American[C]. 1996

[18]Brest J, Greiner S, Boskovic B. Self-adapting control parameters in differential evolution: a comparative study on numerical benchmark problems [J].IEEE Transactions on Evolutionary Computation, 2006,10(6): 646-657

[19]Shi Y, Eberhart R C. Empirical study of particle swarm optimization//Proceedings of the 1999 Congress on Evolutionary Computation[C].1999

[20]Krishna Teerth Chaturvedi, Manjaree Pandit, Laxmi Srivastava.Particle swarm optimization with time varying acceleration coef fi cients for non-convex economic power dispatch [J].Electrical Power and Energy Systems, 2009, 31:249-257