楊 薈,周建棟,李 想

(北京化工大學(xué)經(jīng)濟(jì)管理學(xué)院,北京100029)

1 引 言

隨著我國(guó)城市化進(jìn)程的不斷推進(jìn),人們對(duì)省時(shí)高效,安全便捷的交通出行工具的需求不斷增加.地鐵作為現(xiàn)代化城市公共交通系統(tǒng)的重要組成部分,相比于其他交通形式(如公交車,出租車等),具有準(zhǔn)時(shí)性高,載客量大,發(fā)車頻率和集約化程度高的特性,逐漸成為人們出行方式的首選.然而,地鐵運(yùn)營(yíng)在滿足高運(yùn)營(yíng)效率和安全性要求的同時(shí),供電需求和耗電量也在大幅增加.據(jù)資料顯示,2012年北京地鐵系統(tǒng)總耗電量超過(guò)10億kW·h,其中列車牽引能耗占總耗電量的70%左右[1].據(jù)《中國(guó)城市軌道交通年度報(bào)告》課題組和中國(guó)土木工程學(xué)會(huì)城市軌道交通技術(shù)工作委員會(huì)初步統(tǒng)計(jì),截至2015年底我國(guó)內(nèi)地有25個(gè)城市建成并運(yùn)營(yíng)110條城市軌道交通線路,運(yùn)營(yíng)總里程達(dá)3 293 km.隨著城市地鐵網(wǎng)絡(luò)的快速擴(kuò)展,高能源消耗帶來(lái)的高運(yùn)行成本問(wèn)題將越來(lái)越嚴(yán)重.無(wú)論是從滿足可持續(xù)發(fā)展要求還是減輕地鐵運(yùn)營(yíng)成本的角度出發(fā),地鐵節(jié)能優(yōu)化都勢(shì)在必行.再生制動(dòng)是將列車制動(dòng)過(guò)程中的動(dòng)能轉(zhuǎn)換成電能并提供給處于同一供電站的其他列車加速的一種能量回收技術(shù).在確保地鐵運(yùn)營(yíng)安全性和服務(wù)水平的前提下,通過(guò)充分利用再生制動(dòng)技術(shù)對(duì)地鐵時(shí)刻表進(jìn)行優(yōu)化,從而提高再生能量的利用率,成為實(shí)現(xiàn)地鐵節(jié)能運(yùn)行的主要途徑之一.也成為時(shí)刻表優(yōu)化問(wèn)題的另一目的.

目前,國(guó)內(nèi)外關(guān)于列車再生能量利用的研究,主要包括運(yùn)行控制策略優(yōu)化和時(shí)刻表優(yōu)化兩個(gè)方面[2].在列車運(yùn)行控制策略優(yōu)化研究中,Gordon等[3]和Chen等[4]從再生能量存儲(chǔ)的角度出發(fā),研究了基于地鐵能耗最小化的列車運(yùn)行控制方法.Li等[5]提出了地鐵系統(tǒng)的動(dòng)態(tài)調(diào)度和控制策略,并構(gòu)建凸優(yōu)化模型以提高地鐵節(jié)能性能.

在時(shí)刻表優(yōu)化研究中,周劍斌等[6]基于廣州地鐵1號(hào)線的真實(shí)運(yùn)營(yíng)數(shù)據(jù)提出一種新的時(shí)刻表編制方法,通過(guò)實(shí)行該時(shí)刻表可使每列列車日平均用電節(jié)約0.95 W·h.顏邦杰等[7]兼顧運(yùn)營(yíng)成本和服務(wù)水平,通過(guò)調(diào)整行車預(yù)留時(shí)間得到地鐵節(jié)能時(shí)刻表,并在文章最后給出了地鐵節(jié)能運(yùn)行的相關(guān)建議.Ramos等[8]通過(guò)建立混合整數(shù)規(guī)劃模型對(duì)同一供電區(qū)間內(nèi)列車的加速和制動(dòng)時(shí)間進(jìn)行協(xié)調(diào),從而達(dá)到提高再生能量利用率的目的.Pe?a-Alcaraz等[9]以最大化加速列車和制動(dòng)列車間的重疊時(shí)間為目標(biāo)建立數(shù)學(xué)規(guī)劃模型,實(shí)驗(yàn)結(jié)果表明該最優(yōu)時(shí)刻表不僅能夠有效地節(jié)約能耗,而且不會(huì)對(duì)乘客服務(wù)質(zhì)量和投資成本造成消極影響.Yang等[10]考慮更多的現(xiàn)實(shí)因素(如:真實(shí)速度曲線,供電系統(tǒng)結(jié)構(gòu)等)并在此基礎(chǔ)上提出一種新的節(jié)能調(diào)度方法.該方法通過(guò)調(diào)整列車的停站時(shí)間可使能耗年均降低6.97%.考慮均勻的發(fā)車間隔,步兵等[11]從節(jié)能運(yùn)行的角度出發(fā),將單車節(jié)能運(yùn)行和多車協(xié)調(diào)運(yùn)行進(jìn)行統(tǒng)一規(guī)劃,構(gòu)建以最小化總能耗為目標(biāo)的時(shí)刻表優(yōu)化模型.Yang等[12]以最大化重疊時(shí)間為目標(biāo)建立整數(shù)規(guī)劃模型,通過(guò)優(yōu)化停站時(shí)間和固定的發(fā)車間隔可使早高峰和晚高峰時(shí)段的重疊時(shí)間分別提高22.06%和15.19%.Li等[13]和Yang等[14]均對(duì)速度曲線進(jìn)行了優(yōu)化,并以降低能耗為目標(biāo)建立了時(shí)刻表優(yōu)化模型.Albrecht等[15,16]肯定了Li等[13]提出的一體化節(jié)能優(yōu)化方法,并指出調(diào)整繁忙線路的停站時(shí)間可顯著提高再生能量利用率.Xu等[17]首先對(duì)站間的乘客流量變化進(jìn)行分析并對(duì)乘客等待時(shí)間進(jìn)行定義,進(jìn)而建立了以最小化乘客等待時(shí)間和最小化能量消耗為目標(biāo)的雙目標(biāo)時(shí)刻表優(yōu)化模型.此外,Li等[18]將站間運(yùn)行時(shí)間作為隨機(jī)變量,基于傳統(tǒng)動(dòng)力學(xué)方程所描述的速度曲線建立了以最大化再生能量利用值的期望值為目標(biāo)的隨機(jī)時(shí)刻表優(yōu)化模型.

傳統(tǒng)研究大多基于動(dòng)力學(xué)方程描述列車的站間運(yùn)動(dòng)軌跡.然而列車在實(shí)際運(yùn)行過(guò)程中受線路狀況、機(jī)車條件、電力供應(yīng)及人為操作等因素的影響,動(dòng)力學(xué)方程無(wú)法準(zhǔn)確地描述其在站間的真實(shí)行為.因此,為了更加貼合地鐵運(yùn)行的實(shí)際情形,本文運(yùn)用真實(shí)運(yùn)動(dòng)軌跡數(shù)據(jù)替代動(dòng)力學(xué)方程描述的速度曲線,提出基于數(shù)據(jù)驅(qū)動(dòng)構(gòu)建模型的方法.另外,由于列車的發(fā)車間隔與再生能量利用之間存在潛在關(guān)系,合理選擇發(fā)車間隔也成為地鐵時(shí)刻表節(jié)能優(yōu)化問(wèn)題的主要優(yōu)化途徑之一.因此,本文還研究了非均勻發(fā)車間隔的時(shí)刻表對(duì)再生能量利用的影響問(wèn)題.采用隨機(jī)變量的樂(lè)觀值測(cè)度對(duì)不同延遲情形下再生能量利用值進(jìn)行數(shù)學(xué)描述,以確保含有隨機(jī)參數(shù)的約束條件以一定的置信水平成立,從而建立了基于數(shù)據(jù)驅(qū)動(dòng)的地鐵節(jié)能時(shí)刻表隨機(jī)機(jī)會(huì)約束規(guī)劃模型.考慮到所建模型含有較多的決策變量,且目標(biāo)函數(shù)和約束條件的構(gòu)成較復(fù)雜,運(yùn)用傳統(tǒng)優(yōu)化方法進(jìn)行求解將耗費(fèi)大量計(jì)算時(shí)間,難以滿足地鐵時(shí)刻表優(yōu)化的敏捷性要求,因此本文采用“優(yōu)化過(guò)程設(shè)計(jì)科學(xué),全局尋優(yōu)能力強(qiáng),能夠適應(yīng)大規(guī)模多峰優(yōu)化問(wèn)題且計(jì)算耗費(fèi)時(shí)間合理”的猴群智能算法對(duì)模型進(jìn)行了求解.結(jié)果表明基于本模型的地鐵運(yùn)營(yíng)時(shí)刻表可使早高峰和晚高峰時(shí)段的再生能量利用量樂(lè)觀值顯著提高.與北京地鐵亦莊線現(xiàn)行時(shí)刻表相比,通過(guò)優(yōu)化非均勻的發(fā)車間隔可至少提高再生能量利用量8.44%.

2 隨機(jī)機(jī)會(huì)約束規(guī)劃模型的建立

2.1 基礎(chǔ)知識(shí)

為了求解含有隨機(jī)變量的數(shù)學(xué)規(guī)劃問(wèn)題,Charnes等[19]在1959年首次提出了機(jī)會(huì)約束規(guī)劃的概念,其核心思想是確保含有隨機(jī)參數(shù)的約束條件以一定的置信水平成立.本節(jié)重點(diǎn)介紹由Liu[20]給出的maximax型機(jī)會(huì)約束規(guī)劃模型.

先介紹隨機(jī)變量的樂(lè)觀值的定義.

定義1若ξ是一個(gè)隨機(jī)變量,且α∈(0,1],則ξ的α?樂(lè)觀值定義為

其中Pr為概率測(cè)度.

從定義1可知隨機(jī)變量ξ至少以概率α大于或等于ξsup(α),也即α?樂(lè)觀值是滿足條件的的最大值.

設(shè)x是一個(gè)決策向量,ξ是一個(gè)隨機(jī)向量,f(x,ξ),gj(x,ξ),j=1,2,...,p分別是含有隨機(jī)參數(shù)的目標(biāo)函數(shù)與約束函數(shù).若決策者希望最大化目標(biāo)函數(shù)的α?樂(lè)觀值,可構(gòu)建如下maximax型機(jī)會(huì)約束規(guī)劃模型

其中是目標(biāo)函數(shù)f(x,ξ)的α?樂(lè)觀值.

隨機(jī)機(jī)會(huì)約束規(guī)劃模型的求解方法主要分為兩類,第一類是分析目標(biāo)函數(shù)和約束函數(shù)的確定表達(dá)式,進(jìn)而將其轉(zhuǎn)化為確定的等價(jià)模型進(jìn)行求解.在大多數(shù)情況下,由于隨機(jī)機(jī)會(huì)約束規(guī)劃模型中的函數(shù)具有不確定性,其近似表達(dá)式往往很難得到,因此這一類方法不具有廣泛的應(yīng)用性.第二類方法是將神經(jīng)網(wǎng)絡(luò),Monte Carlo模擬和啟發(fā)式算法相結(jié)合對(duì)原模型進(jìn)行求解,該方法也稱之為混合智能算法[21].典型的啟發(fā)式算法包括遺傳算法[21,22],粒子群算法[23?25],猴群算法[26]等.

2.2 模型建立

當(dāng)前地鐵供電大多采用單邊供電的模式.在該模式下,電能在不同供電區(qū)間內(nèi)無(wú)法傳遞,并且上行和下行的供電是相互獨(dú)立的.如圖1所示,列車3再生的電能不能用于列車4的加速,列車1和列車2之間的能量無(wú)法互相傳遞.相反,如果某列車處于牽引階段時(shí),恰好有另一輛同行列車位于同一供電區(qū)間且處于制動(dòng)階段,那么重疊時(shí)間內(nèi)制動(dòng)列車再生的電能就可以供給牽引列車用于加速.如圖2所示為同一供電區(qū)間內(nèi)兩輛連續(xù)列車在站間運(yùn)行的速度曲線,圖中,列車i正在加速離開j站,同一變電區(qū)間內(nèi)的另一列車i+1正在制動(dòng)進(jìn)入j站,此時(shí)陰影部分代表的就是再生電能被利用的重疊時(shí)間部分.為了提高再生能量的利用率,本文的目標(biāo)就是最大化牽引列車和制動(dòng)列車的重疊時(shí)間,使得更多的再生能量得以被牽引列車?yán)?從而降低列車運(yùn)行的整體能耗.文中所用參數(shù)及符號(hào)見表1.

圖1 供電系統(tǒng)結(jié)構(gòu)示意圖Fig.1 The illustration on power supply system

圖2 同一供電區(qū)間內(nèi)兩輛連續(xù)列車間的速度–時(shí)間分布圖Fig.2 Speed-time profile of two successive trains within a substation

表1 符號(hào)系統(tǒng)Table 1 Symbol system

通常將列車的運(yùn)行軌跡分為加速,巡航,惰行和制動(dòng)四個(gè)工況[13].由于人為操縱因素和外界因素的影響,在列車行駛過(guò)程中加速,巡航和惰行工況是頻繁交替出現(xiàn)的.然而在對(duì)再生能量利用值進(jìn)行計(jì)算時(shí),再生能量的收集僅來(lái)源于制動(dòng)工況,再生能量的利用主要來(lái)源于加速工況,巡航和惰行工況極少涉及.因此,將復(fù)雜的工況交替過(guò)程按四個(gè)工況進(jìn)行逐一劃分是沒有必要的,并且還會(huì)增加運(yùn)算的復(fù)雜度和計(jì)算時(shí)間,不符合地鐵時(shí)刻表優(yōu)化的敏捷性要求.因此,本文將只包含加速,巡航和惰行工況交替出現(xiàn)的運(yùn)行階段稱為牽引階段,將制動(dòng)工況出現(xiàn)的運(yùn)行階段稱為制動(dòng)階段.并以制動(dòng)階段結(jié)束時(shí)刻為結(jié)點(diǎn)對(duì)站間的運(yùn)動(dòng)軌跡進(jìn)行劃分,使每個(gè)站間由若干個(gè)連續(xù)的區(qū)間構(gòu)成,并且每個(gè)區(qū)間由一個(gè)牽引階段和一個(gè)制動(dòng)階段構(gòu)成,以圖2為例,列車i在站間j的運(yùn)動(dòng)軌跡劃分為區(qū)間1和區(qū)間2,其中,斜紋代表列車處于加速工況,黑色代表列車處于巡航或惰行工況,斜紋和黑色部分合起來(lái)稱為牽引階段,波點(diǎn)代表列車處于制動(dòng)工況,也即制動(dòng)階段,并且每個(gè)區(qū)間都由一個(gè)牽引階段和一個(gè)制動(dòng)階段組成.

2.3 數(shù)據(jù)驅(qū)動(dòng)

列車在實(shí)際運(yùn)行過(guò)程中往往受線路狀況,機(jī)車條件及人為因素等影響,動(dòng)力學(xué)方程不能準(zhǔn)確地描述其運(yùn)動(dòng)軌跡的這一變化特點(diǎn).以北京地鐵亦莊線榮昌東路–同濟(jì)南路站間為例,圖3(a)和圖3(b)分別展示了由傳統(tǒng)動(dòng)力學(xué)方程描述的速度曲線

圖3 傳統(tǒng)動(dòng)力學(xué)方程所描述速度曲線和真實(shí)速度曲線對(duì)比圖Fig.3 Comparison chart of speed profile described by traditional motion equation and real-life data

和真實(shí)的速度曲線.因此,為了使得所建模型可靠性高并且具有現(xiàn)實(shí)意義,本文基于數(shù)據(jù)驅(qū)動(dòng)構(gòu)建地鐵時(shí)刻表優(yōu)化模型,用真實(shí)的速度數(shù)據(jù)取代動(dòng)力學(xué)方程描述的速度曲線.

2.4 非均勻發(fā)車間隔

地鐵運(yùn)營(yíng)公司通常為列車分配均勻的發(fā)車間隔.實(shí)際上,列車的發(fā)車間隔與再生能量利用之間存在潛在關(guān)系.因此,優(yōu)化非均勻的發(fā)車間隔也是提高再生能量利用率的途徑之一.并且,隨著優(yōu)化能力與運(yùn)算速度的提升,實(shí)現(xiàn)列車發(fā)車間隔的實(shí)時(shí)調(diào)整也成為了必然趨勢(shì).綜合以上因素,為了盡可能地提高再生能量利用率,本文考慮非均勻的發(fā)車間隔.

2.5 延遲分析

在地鐵的實(shí)際運(yùn)營(yíng)過(guò)程中,地鐵時(shí)刻表常常會(huì)受到乘客客流量增大等因素的影響而發(fā)生延遲.為了盡可能避免延遲現(xiàn)象的出現(xiàn),地鐵運(yùn)營(yíng)公司通常會(huì)對(duì)站間運(yùn)行時(shí)間進(jìn)行調(diào)節(jié),從而縮小與既定時(shí)刻表的差距.本文考慮三種不同的延遲情景,分別為無(wú)延遲,輕微延遲和嚴(yán)重延遲,并用pqj(q=1,2,3)分別表示它們?cè)诟髡鹃g的發(fā)生概率.在每種情景下,地鐵運(yùn)營(yíng)公司將分別制定相應(yīng)的站間運(yùn)行策略和列車需要遵循的速度曲線.由于延遲情況的出現(xiàn)是隨機(jī)的,從而導(dǎo)致不同策略的選擇是隨機(jī)的.因此,分別給出不同情景下第j站站間的第n個(gè)區(qū)間的牽引時(shí)間和制動(dòng)時(shí)間的表達(dá)式為

2.6 目標(biāo)函數(shù)分析

在地鐵系統(tǒng)中,每輛列車均由動(dòng)車和拖車車廂編制而成,本文給出列車運(yùn)行所受阻力的公式[27]如下

定義決策變量 X=(h1,h2,...,hI,d1,d2,...,dI?1),隨機(jī)向量 ξ =(ξ1,ξ2,...ξJ?1),其中 hi為每輛列車的發(fā)車間隔,dj為每站的停站時(shí)間,ξj是隨機(jī)變量,i=1,2,...,I;j=1,2,...,J?1.當(dāng)列車離開j站,即n=1時(shí),定義列車在各站的出發(fā)時(shí)間為

假設(shè)第1輛列車從0時(shí)刻出發(fā),那么第i輛列車經(jīng)過(guò)了i?1個(gè)發(fā)車間隔,并且經(jīng)過(guò)了前j?1個(gè)站的運(yùn)行時(shí)間和前j個(gè)站的停站時(shí)間后,從第j站出發(fā).因此,當(dāng)n=1時(shí),式(5)表示列車i在第j站的出發(fā)時(shí)間;當(dāng)n=2,3,...,Nj時(shí),Tijn(X,ξ)表示列車i在第j站站間的第n個(gè)區(qū)間開始牽引的時(shí)間.不難發(fā)現(xiàn),列車i在第j站站間的第n個(gè)區(qū)間開始牽引的時(shí)間與列車i在第j站站間的第n?1個(gè)區(qū)間結(jié)束制動(dòng)的時(shí)間相同.進(jìn)而,可以給出列車在區(qū)間n開始制動(dòng),結(jié)束制動(dòng),以及在下一區(qū)間開始牽引的時(shí)間為

已知列車i在各站間的速度曲線,根據(jù)動(dòng)能定理,可以得到t時(shí)刻列車i在第j站站間的第n個(gè)區(qū)間產(chǎn)生的動(dòng)能為

阻力消耗的能量為

因此,對(duì)于每一個(gè)時(shí)間單元[t,t+1],列車i加速需要的能量可表示為

根據(jù)式(7)可知,只有當(dāng) Tijn(X,ξ,t)≤t≤Bijn(X,ξ,t)時(shí),即列車處于牽引階段時(shí),才會(huì)有加速工況出現(xiàn),并需要為加速過(guò)程提供能量.其他時(shí)間均處于制動(dòng)階段,加速需要的能量為0.

牽引階段不僅包括加速工況,還有巡航和惰行工況,其中巡航和惰行工況動(dòng)能為非正值.為了只求得列車加速時(shí)需要的能量,只需得到Fijn(X,ξ,t)中非負(fù)的值即可.因此,將式(7)修正為

同樣的,對(duì)于每一個(gè)時(shí)間單元[t,t+1],列車i制動(dòng)可以回收的能量可表示為

在式(12)中,只有當(dāng)列車處于 Bijn(X,ξ,t)≤t≤Aijn(X,ξ,t)時(shí),即制動(dòng)階段時(shí),才會(huì)產(chǎn)生回收能量,其他時(shí)間列車沒有再生能量回收.

若考慮共有I輛列車,每輛列車在同一區(qū)間具有相同的速度曲線且在同一車站具有相同的停站時(shí)間,則整條線路上全部列車的再生能量利用量可表示為

其中

2.7 隨機(jī)機(jī)會(huì)約束規(guī)劃模型

本文以最大化再生能量利用量為目標(biāo)函數(shù),即f(X,ξ).通過(guò)最大化目標(biāo)函數(shù)f(X,ξ)的α?樂(lè)觀值,達(dá)到提高再生能量利用率的目的.本文建立的隨機(jī)機(jī)會(huì)約束規(guī)劃模型如下.

模型(14)以最大化再生能量利用量的α?樂(lè)觀值為目標(biāo).其中α和β是事先給定的置信水平,是再生能量利用量函數(shù)的α?樂(lè)觀值

第二個(gè)約束保證總的旅行時(shí)間以β的置信水平滿足時(shí)間窗約束;第三個(gè)和第四個(gè)約束確保停站時(shí)間和每輛車的發(fā)車間隔符合時(shí)間窗的約束,以更好地滿足乘客的需求;第五個(gè)約束是為了保證最優(yōu)解為整數(shù).

本文所研究的時(shí)刻表優(yōu)化方法的創(chuàng)新點(diǎn)在于:1)發(fā)車間隔是非均勻的,通過(guò)調(diào)整列車的發(fā)車間隔,使得再生能量的利用率提高;2)建立基于數(shù)據(jù)驅(qū)動(dòng)的模型,以真實(shí)數(shù)據(jù)取代傳統(tǒng)動(dòng)力學(xué)方程描述的速度曲線,使得模型的建立更加準(zhǔn)確和可靠;3)從樂(lè)觀值的角度出發(fā),以最大化再生能量利用量樂(lè)觀值為目標(biāo),建立隨機(jī)機(jī)會(huì)約束規(guī)劃模型.

3 猴群算法

本文建立的基于數(shù)據(jù)驅(qū)動(dòng)的地鐵節(jié)能時(shí)刻表隨機(jī)機(jī)會(huì)約束規(guī)劃模型具有決策變量多,目標(biāo)函數(shù)和約束條件復(fù)雜的特點(diǎn),因此使用傳統(tǒng)優(yōu)化方法進(jìn)行求解將花費(fèi)大量計(jì)算時(shí)間.為了求解所建隨機(jī)機(jī)會(huì)約束規(guī)劃模型,需要使用一種優(yōu)化過(guò)程設(shè)計(jì)科學(xué)有效,不會(huì)陷入局部最優(yōu),能夠適應(yīng)所建模型結(jié)構(gòu)且計(jì)算時(shí)間合理的智能算法.猴群算法(Monkey Algorithm,MA)是由Zhao等[28]于2008年提出的一種在求解多維,大規(guī)模,多峰優(yōu)化問(wèn)題時(shí)體現(xiàn)出明顯優(yōu)勢(shì)的新型群智能優(yōu)化算法.猴群算法通過(guò)模擬自然界中猴群爬山的過(guò)程,設(shè)計(jì)了初始化,爬過(guò)程,望過(guò)程,跳過(guò)程和終止條件等步驟.其中初始化利用隨機(jī)方法產(chǎn)生每個(gè)猴子的初始位置,爬過(guò)程利用目標(biāo)函數(shù)的偽梯度來(lái)逐步改善優(yōu)化模型的目標(biāo)函數(shù)值,望過(guò)程是在猴群達(dá)到局部最優(yōu)后尋找相鄰區(qū)域內(nèi)更優(yōu)解的過(guò)程,而跳過(guò)程則迫使猴群由當(dāng)前區(qū)域轉(zhuǎn)移到新的區(qū)域,以避免算法局部最優(yōu).猴群算法優(yōu)越的全局優(yōu)化性能已被廣泛應(yīng)用于求解各種復(fù)雜優(yōu)化問(wèn)題,如:輸電線擴(kuò)展規(guī)劃問(wèn)題[29],入侵檢測(cè)問(wèn)題[30],傳感器優(yōu)化問(wèn)題等[31,32].本文在猴群算法望過(guò)程之后設(shè)計(jì)信息共享機(jī)制來(lái)避免猴群中信息不對(duì)稱造成的局部最優(yōu)解優(yōu)劣差異,從而更好的找尋全局最優(yōu)解.猴群算法的流程圖如圖4所示.

圖4 猴群算法流程圖Fig.4 The flow chart of Monkey Algorithm

3.1 解的表示及初始化

設(shè)猴群規(guī)模為正整數(shù)M,第m只猴子的當(dāng)前位置用向量Xm=(xm1,xm2,...,xmn)表示,m=1,2,...,M.每只猴子的位置實(shí)際上對(duì)應(yīng)所建立隨機(jī)機(jī)會(huì)約束規(guī)劃模型的一個(gè)決策向量,也即代表模型的一個(gè)可行解.Xm由兩部分決策變量組成:列車的發(fā)車間隔時(shí)間和列車在車站j的停站時(shí)間,也即Xm=(h1,h2,...,hI;d1,d2,...,dJ?1)且n=I+J?1.

從可行域中隨機(jī)產(chǎn)生M個(gè)可行解(猴子)作為猴群算法的初始種群,初始化過(guò)程的算法設(shè)計(jì)如下:

步驟1對(duì)第m 只猴子,隨機(jī)產(chǎn)生解向量其中且i=1,2,...,I,產(chǎn)生隨機(jī)列車發(fā)車間隔;且j=1,2,...,J?1,產(chǎn)生隨機(jī)列車停留時(shí)間.

步驟2若Xtemp滿足模型的約束條件,則令Xm=Xtemp為第m只猴子(可行解);若不滿足,則重復(fù)步驟1直至得到第m只猴子,也即第m個(gè)可行解;

步驟3重復(fù)步驟1至步驟2,直至產(chǎn)生M只猴子.

3.2 爬過(guò)程

猴群算法的爬過(guò)程是一個(gè)通過(guò)迭代不斷改變猴子所處位置以改善模型目標(biāo)函數(shù)值的過(guò)程.在本文中,爬過(guò)程可通過(guò)連續(xù)迭代來(lái)優(yōu)化所建隨機(jī)機(jī)會(huì)約束規(guī)劃模型的可行解,使其不斷逼近最大再生能量利用量的取值.為避免過(guò)多地使用目標(biāo)函數(shù)的梯度信息,與采用梯度法的牛頓法不同,猴群算法采用偽梯度法[33,34]的思想,僅計(jì)算當(dāng)前相鄰位置的目標(biāo)函數(shù)值,以逐步移動(dòng)的方式設(shè)計(jì)算法的爬過(guò)程.設(shè)第m只猴子的當(dāng)前位置為Xm=(xm1,xm2,...,xmn),m=1,2,...,M.爬過(guò)程設(shè)計(jì)如下:

步驟 1對(duì)第 m 只猴子,隨機(jī)產(chǎn)生向量 ΔXm=(Δxm1,Δxm2,...,Δxmn)T,其中 Δxmk是從區(qū)間[?φ,φ]中隨機(jī)產(chǎn)生的整數(shù),φ為爬過(guò)程的步長(zhǎng),m=1,2,...,M,k=1,2,...,n;

步驟2若則

步驟3重復(fù)步驟1至步驟2,直到達(dá)到設(shè)定的爬次數(shù)Nc.

猴群算法中定義的參數(shù)φ稱為爬步長(zhǎng),φ越小,求解出的結(jié)果越精確,但往往導(dǎo)致較長(zhǎng)計(jì)算時(shí)間.對(duì)于規(guī)模較大的求解模型,φ值可適當(dāng)設(shè)置得較大或采取變動(dòng)的爬步長(zhǎng)值.

3.3 望過(guò)程

在爬過(guò)程結(jié)束之后,猴群中的每只猴子都到達(dá)了各自所處位置附近的最高峰,也即求解的隨機(jī)機(jī)會(huì)約束規(guī)劃模型的目標(biāo)函數(shù)達(dá)到了局部最優(yōu)解,但未必是全局最優(yōu).猴群算法的望過(guò)程是使每只猴子在局部最優(yōu)的位置上,在一定視野范圍ρ內(nèi)向周圍眺望,觀察在其鄰近區(qū)域內(nèi)是否有比當(dāng)前所處的山峰更高的山峰,也即尋找是否存在比目標(biāo)函數(shù)的局部最優(yōu)解更優(yōu)的更優(yōu)解.若猴子發(fā)現(xiàn)更高的山峰,就會(huì)從當(dāng)前的山峰跳向更高峰所在的區(qū)域,繼續(xù)通過(guò)爬過(guò)程找到新區(qū)域內(nèi)的最高峰.對(duì)于第m只猴子,m=1,2,...,M,其望過(guò)程如下:

步驟1設(shè)猴子眺望的視野為ρ;

步驟2設(shè)隨機(jī)產(chǎn)生整數(shù)1,2,...,n;

步驟3若滿足模型約束條件,則使

步驟4重復(fù)步驟2至步驟3,直到達(dá)到設(shè)定的望次數(shù)Nw.

猴群算法望過(guò)程中定義的視野范圍參數(shù)ρ代表猴子所能眺望到的最遠(yuǎn)距離.ρ越小,求解出的結(jié)果越精確,但往往增加運(yùn)算時(shí)間.對(duì)于規(guī)模較大的求解模型,ρ值可適當(dāng)設(shè)置得較大或采取變動(dòng)的視野范圍.

3.4 信息共享機(jī)制

猴群算法在經(jīng)過(guò)爬-望過(guò)程后,猴群中的每只猴子各占山頭(達(dá)到局部最優(yōu)解),然而每只猴子位置的不同造成局部最優(yōu)解也存在優(yōu)劣差別.通過(guò)信息共享機(jī)制,以便讓猴群之間的信息進(jìn)行交流,從而使處于較高位置的猴子與處于較低位置的猴子進(jìn)行合作,增強(qiáng)了算法的尋優(yōu)能力和解的多樣性.

設(shè)在當(dāng)次迭代的所有猴子中的最高位置為第m只猴子的當(dāng)前位置Xm=(xm1,xm2,...,xmn),m=1,2,...,M,設(shè)計(jì)猴群算法的信息共享機(jī)制如下:

步驟1隨機(jī)產(chǎn)生實(shí)數(shù)γ∈[0,1];

步驟2記1,2,...,n;

步驟3若滿足模型的所有約束條件,則使

步驟4對(duì)更新的X1,X2,...,XM重復(fù)爬過(guò)程.

3.5 跳過(guò)程

猴群算法跳過(guò)程的主要目的是使猴子從當(dāng)前的區(qū)域轉(zhuǎn)向新的搜索區(qū)域,防止算法陷入局部極值導(dǎo)致搜索停滯.以當(dāng)前猴群的重心所在位置為跳躍支點(diǎn),每只猴子沿著當(dāng)前所處位置指向支點(diǎn)的方向跳向新的搜索區(qū)域.設(shè)猴子跳躍的最大范圍為[c,d],跳過(guò)程設(shè)計(jì)如下:

步驟1隨機(jī)產(chǎn)生實(shí)數(shù)η∈[c,d];

步驟2計(jì)算=1,2,...,n,記g=(g1,g2,...,gn)T為跳躍支點(diǎn);

步驟3計(jì)算xmk+Round(η|gk? xmk|),m=1,2,...,M,k=1,2,...,n;

步驟 4若滿足所求解模型的所有約束條件,則使并重復(fù)爬過(guò)程.

3.6 終止條件

在經(jīng)過(guò)以上爬過(guò)程,望過(guò)程,信息共享,跳過(guò)程步驟后,猴群算法完成第一次迭代,猴群中所有猴子的位置得到更新.重復(fù)迭代循環(huán)一定次數(shù)Imax直到最優(yōu)解保持不變后,算法終止并輸出最優(yōu)解.

4 數(shù)值算例

將北京地鐵亦莊線(見圖5)的數(shù)據(jù)應(yīng)用于模型(13),并給出數(shù)值算例證明模型(13)的有效性和意義.

圖5 北京地鐵亦莊線Fig.5 Beijing Yizhuang subway line

北京地鐵亦莊線將北京經(jīng)濟(jì)技術(shù)開發(fā)區(qū)與北京主城區(qū)緊緊連接在一起,為高端產(chǎn)業(yè)的發(fā)展提供了堅(jiān)實(shí)的交通保障,也使開發(fā)區(qū)的企業(yè)員工和社區(qū)居民有了更為便捷的交通方式.亦莊開發(fā)區(qū)19萬(wàn)余名企業(yè)員工,近4萬(wàn)社區(qū)居民可以近距離搭乘地鐵直達(dá)市區(qū).亦莊線共有14個(gè)車站和6個(gè)變電站,其中亦莊文化園,榮昌東街和同濟(jì)南路站客流量較大.因此,針對(duì)以上三個(gè)車站設(shè)置了無(wú)延遲,輕微延遲和嚴(yán)重延遲三種情景,并設(shè)置了相應(yīng)的可能性,即p1j=0.7,p2j=0.2和p3j=0.1.亦莊線當(dāng)前的運(yùn)行時(shí)刻表如表2所示.根據(jù)亦莊線的真實(shí)速度曲線,表3給出了各站間的牽引和制動(dòng)時(shí)間.

表2 亦莊線當(dāng)前時(shí)刻表Table 2 The current timetable for Beijing Yizhuang subway line

模型中的其他參數(shù)設(shè)置見表4.

表4 模型中的參數(shù)設(shè)置Table 4 The parameter setting of model

在運(yùn)行猴群算法之前,首先對(duì)猴群算法進(jìn)行收斂性實(shí)驗(yàn).通過(guò)程序迭代,得到猴群算法的相關(guān)參數(shù)設(shè)置如表5所示,算法性能如圖6所示.

表5 猴群算法中的參數(shù)Table 5 The parameters setting of MA

圖6 猴群算法算法迭代性能分析Fig.6 The analysis on iterative performance of MA

然后,分別針對(duì)置信水平α和β進(jìn)行實(shí)驗(yàn),研究其對(duì)最優(yōu)解的影響.固定α=90%,求解當(dāng)β取95%,90%,80%,70%和60%時(shí)的最優(yōu)解.同樣的,再固定β=90%,求解當(dāng)α取95%,90%,80%,70%和60%時(shí)的最優(yōu)解,結(jié)果如表6所示.

從表6中可以看出,隨著α和β的不斷減小,不同延遲情景下再生能量利用量樂(lè)觀值會(huì)不斷增大.因此,本文一方面為了提高再生能量利用量的最壞估計(jì),選取α=95%;另一方面,為了使得總運(yùn)行時(shí)間以更高的置信水平符合時(shí)間窗的約束,提高時(shí)刻表的魯棒性,選取β=95%.但是在其他領(lǐng)域或模型的應(yīng)用中,可以根據(jù)策略和實(shí)際情況的不同選取不同的置信水平.

表6 置信水平α和β對(duì)最優(yōu)解的影響研究Table 6 Influence of confidence level on the optimal solution

根據(jù)乘客流量和需求等特點(diǎn),北京地鐵亦莊線主要出現(xiàn)早高峰,晚高峰和平峰三種現(xiàn)象.早高峰和晚高峰時(shí)期乘客需求較大,列車發(fā)車間隔會(huì)縮短.同時(shí),發(fā)車頻率上升也使得這段時(shí)間內(nèi)的能耗較大.因此,如何降低早高峰和晚高峰的能耗也是地鐵節(jié)能優(yōu)化的重要環(huán)節(jié).本文分別針對(duì)早高峰和晚高峰時(shí)段,設(shè)置發(fā)車間隔的變動(dòng)范圍為[330,390]和[290,350],并設(shè)計(jì)最優(yōu)時(shí)刻表,對(duì)此時(shí)段內(nèi)的再生能量利用量進(jìn)行優(yōu)化.通過(guò)運(yùn)行猴群算法,計(jì)算出早高峰時(shí)段的再生能量利用量樂(lè)觀值為168.33 kW·h.晚高峰時(shí)期的再生能量利用量樂(lè)觀值為171.94 kW·h,其最優(yōu)時(shí)刻表分別如表7和表8所示.

表7 早高峰時(shí)段最優(yōu)時(shí)刻表Table 7 The optimal timetable of morning peak hours

表8 晚高峰時(shí)段最優(yōu)時(shí)刻表Table 8 The optimal timetable of evening peak hours

從表中可以看出,相比當(dāng)前時(shí)刻表,為了增大牽引階段和制動(dòng)階段的重疊時(shí)間,某些車站會(huì)做出調(diào)整使得站間運(yùn)行時(shí)間增大,比如:小紅門站,舊宮站,亦莊火車站站等;同時(shí)也會(huì)有一些車站的停站時(shí)間會(huì)減小,比如:亦莊橋站,榮京東街站,同濟(jì)南路站等.此外,列車的發(fā)車間隔也會(huì)做出相應(yīng)調(diào)整,以早高峰時(shí)段為例,第3,4,7等輛列車的發(fā)車間隔增大,第1,2,5等輛列車的發(fā)車間隔減小;以晚高峰時(shí)段為例,第1,2,4等輛列車的發(fā)車間隔增大,第10輛列車發(fā)車間隔減小.

另一方面,當(dāng)前北京地鐵亦莊線的時(shí)刻表在早高峰和晚高峰時(shí)段的發(fā)車間隔分別固定在350 s和300 s.為了分析非均勻發(fā)車間隔的時(shí)刻表對(duì)于提高再生制動(dòng)能量利用量的意義.令北京地鐵亦莊線在各站的停站時(shí)間與當(dāng)前時(shí)刻表保持一致,運(yùn)用猴群算法優(yōu)化列車非均勻的發(fā)車間隔,使得時(shí)刻表的再生能量利用量樂(lè)觀值最大化.通過(guò)計(jì)算,得到最優(yōu)時(shí)刻表在早高峰和晚高峰時(shí)段的再生能量利用量樂(lè)觀值分別為160.39 kW·h和158.59 kW·h,對(duì)應(yīng)的發(fā)車間隔見表9和表10.

表9 早高峰時(shí)段最優(yōu)發(fā)車間隔Table 9 The optimal headway time of morning peak hours

表10 晚高峰時(shí)段最優(yōu)發(fā)車間隔Table 10 The optimal headway time of evening peak hours

此外,通過(guò)計(jì)算,得到均勻發(fā)車間隔的北京地鐵亦莊線當(dāng)前時(shí)刻表在早高峰和晚高峰時(shí)段的再生能量利用量樂(lè)觀值分別為147.90 kW·h和40.47kW·h.因此,本文的最優(yōu)時(shí)刻表相比于現(xiàn)行時(shí)刻表的再生能量利用量可以提高至少8.44%,其對(duì)比結(jié)果如圖7所示.算例分析的結(jié)果表明優(yōu)化非均勻的發(fā)車間隔有利于提高地鐵系統(tǒng)的再生能量利用率.

圖7 均勻發(fā)車間隔與非均勻發(fā)車間隔再生能量利用量樂(lè)觀值對(duì)比圖Fig.7 Comparison chart of optimistic value between uniform and non-uniform headway time

5 結(jié)束語(yǔ)

本文以最大化再生能量利用量為目標(biāo),基于地鐵真實(shí)的運(yùn)動(dòng)軌跡,構(gòu)建數(shù)據(jù)驅(qū)動(dòng)的時(shí)刻表隨機(jī)機(jī)會(huì)約束規(guī)劃模型,為地鐵運(yùn)營(yíng)管理提供了新的思路和方法.該模型通過(guò)調(diào)整停站時(shí)間和非均勻的發(fā)車間隔,提高再生能量的利用率,從而實(shí)現(xiàn)地鐵節(jié)能運(yùn)行.此外,本文還設(shè)計(jì)了猴群算法進(jìn)行求解并基于北京地鐵亦莊線的真實(shí)數(shù)據(jù)開展了模擬仿真.結(jié)果表明基于本模型的地鐵運(yùn)營(yíng)早高峰和晚高峰時(shí)段的再生能量利用量樂(lè)觀值分別為168.33 kW·h和171.94 kW·h.與北京地鐵亦莊線現(xiàn)行時(shí)刻表相比,通過(guò)優(yōu)化非均勻的發(fā)車間隔可至少提高再生能量利用量8.44%.

本文將機(jī)會(huì)約束規(guī)劃思想與運(yùn)營(yíng)調(diào)度問(wèn)題相結(jié)合,對(duì)列車時(shí)刻表進(jìn)行優(yōu)化,不僅適用于地鐵系統(tǒng),還可擴(kuò)展應(yīng)用于高鐵,城際列車的時(shí)刻表優(yōu)化問(wèn)題中.另外,對(duì)于地鐵系統(tǒng)而言,在進(jìn)行節(jié)能優(yōu)化的同時(shí)還要考慮乘客需求,時(shí)刻表魯棒性等多種因素.因此,如何協(xié)調(diào)這些因素之間的關(guān)系將成為下一步研究的主要方向.

[1]Ogasa M.Energy saving and environmental measures in railway technologies:Example with hybrid electric railway vehicles.IEEJ Transactions on Electrical and Electronic Engineering,2008,3(1):15–20.

[2]Yang X,Li X,Ning B,et al.A survey on energy-efficient train operation for urban rail transit.IEEE Transactions on Intelligent Transportation Systems,2016,17(1):2–13.

[3]Gordon S P,Lehrer D G.Coordinated train control and energy management control strategies//Proceedings of the 1998 ASME/IEEE Joint.IEEE,1998:165–176.

[4]Chen J F,Lin R L,Liu Y C.Optimization of an MRT train schedule:Reducing maximum traction power by using genetic algorithms.IEEE Transactions on Power Systems,2005,20(3):1366–1372.

[5]Li X,Lo H K.Energy minimization in dynamic train scheduling and control for metro rail operations.Transportation,Part B:Methodological,2014,70:269–284.

[6]周劍斌,蘇 浚,何泳斌.地鐵列車運(yùn)行再生能利用的研究.城市軌道交通研究,2004,7(4):33–35.Zhou J B,Su J,He Y B.Utilization of train’s regenerative energy in metro system.Urban Mass Transit,2004,7(4):33–35.(in Chinese)

[7]顏邦杰,張辰秋,林志銘.捷運(yùn)列車排點(diǎn)與節(jié)能.交通信息與安全,2011,29(1):139–144.Yan B J,Zhang C Q,Lin Z M.The MRT train scheduling and energy saving.Journal of Transport Information and Safety,2011,29(1):139–144.(in Chinese)

[8]Ramos A,Pena M T,Fernández A,et al.Mathematical programming approach to underground timetabling problem for maximizing time synchronization//Proceedings of the International Conference on Industrial Engineering and Industrial Management,IEEE,Madrid,Spain,2007:88–95.

[9]Pe?a-Alcaraz M,Fernández A,Cucala A P,et al.Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy//Proceedings of the Institution of Mechanical Engineers,Part F:Journal of Rail and Rapid Transit,2012,226(4):397–408.

[10]Yang X,Chen A,Li X,et al.An energy-efficient scheduling approach to improve the utilization of regenerative energy for metro systems.Transportation Research,Part C:Emerging Technologies,2015,57:13–29.

[11]步 兵,丁 奕,李辰嶺,等.列車控制與行車調(diào)度一體化節(jié)能方法的研究.鐵道學(xué)報(bào),2013,35(12):64–71.Bu B,Ding Y,Li C L,et al.Research on integration of train control and train scheduling.Journal of the China Railway Society,2013,35(12):64–71.(in Chinese)

[12]Yang X,Li X,Gao Z,et al.A cooperative scheduling model for timetable optimization in subway systems.IEEE Transactions on Intelligent Transportation Systems,2013,14(1):438–447.

[13]Li X,Lo H K.An energy-efficient scheduling and speed control approach for metro rail operations.Transportation Research,Part B:Methodological,2014,64:73–89.

[14]Yang X,Li X,Ning B,et al.An optimisation method for train scheduling with minimum energy consumption and travel time in metro rail systems.Transportmetrica B:Transport Dynamics,2015,3(2):79–98.

[15]Albrecht A,Howlett P,Pudney P,et al.The key principles of optimal train control:Part 1:Formulation of the model,strategies of optimal type,evolutionary lines,location of optimal switching points.Transportation Research,Part B:Methodological,2016,94:482–508.

[16]Albrecht A,Howlett P,Pudney P,et al.The key principles of optimal train control:Part 2:Existence of an optimal strategy,the local energy minimization principle,uniqueness,computational techniques.Transportation Research,Part B:Methodological,2016,94:509–538.

[17]Xu X, Li K, Li X. A multi-objective subway timetable optimization approach with minimum passenger time and energy consumption.Journal of Advanced Transportation.Journal of Advanced Transportation,2016,50(1):69–95.

[18]Li X,Yang X.A stochastic timetable optimization model in subway systems.International Journal of Uncertainty,Fuzziness and Knowledge-Based Systems,2013,21(S):1–15.

[19]Charnes A,Cooper W W.Chance-constrained programming.Management Science,1959,6(1):73–79.

[20]Liu B.Uncertain Programming.Hoboken,New Jersey:John Wiley&Sons,Inc.,1999.

[21]劉寶碇,趙瑞清,王 剛.不確定規(guī)劃及應(yīng)用.北京:清華大學(xué)出版社,2006.Liu B D,Zhao R Q,Wang G.Uncertain Programming and Application.Beijing:Tsinghua University Press,2006.(in Chinese)

[22]周愛蓮,李旭宏,毛海軍.一類企業(yè)物流中心動(dòng)態(tài)選址模型研究.系統(tǒng)工程學(xué)報(bào),2011,26(3):360–366.Zhou A L,Li X H,Mao H J.Research on a multiple-period dynamic location model of enterprise logistics centers.Journal of Systems Engineering,2011,26(3):360–366.(in Chinese)

[23]肖 寧.求解隨機(jī)機(jī)會(huì)約束規(guī)劃的混合智能算法.計(jì)算機(jī)工程與應(yīng)用,2010,46(22):43–46.Xiao N.Stochastic chance-constrained programming problems with hybrid intelligent algorithm.Computer Engineering and Applications,2010,46(22):43–46.(in Chinese)

[24]石 莉,楊善林,馬 英,等.一種新的組合權(quán)重集結(jié)方法及合理性評(píng)價(jià)研究.系統(tǒng)工程學(xué)報(bào),2012,27(4):481–491.Shi L,Yang S L,Ma Y,et al.A novel method of combination weighting for multiple attribute decision making.Journal of Systems Engineering,2012,27(4):481–491.(in Chinese)

[25]王洪峰,王 娜,汪定偉,等.一種求解多峰優(yōu)化問(wèn)題的改進(jìn)Species粒子群算法.系統(tǒng)工程學(xué)報(bào),2012,27(6):854–864.Wang H F,Wang N,Wang D W,et al.Improved species-based particle swarm optimizer for multi-modal optimization problems.Journal of Systems Engineering,2012,27(6):854–864.(in Chinese)

[26]趙瑞清,郝士鵬.一類新的模糊約束滿足問(wèn)題的建模與求解.系統(tǒng)工程學(xué)報(bào),2010,25(3):415–420.Zhao R Q,Hao S P.Modeling and solving a novel class of fuzzy constraint satisfaction problems.Journal of Systems Engineering,2010,25(3):415–420.(in Chinese)

[27]Xiao X.The Record of Vehicle Condition for Beijing Metro Yizhuang line.Technical Report.Beijing Mass Transit Railway Operation Corporation Limited,2014.

[28]Zhao R,Tang W.Monkey algorithm for global numerical optimization.Journal of Uncertain Systems,2008,2(3):165–176.

[29]王靖然,余貽鑫,曾 沅.離散猴群算法及其在輸電網(wǎng)擴(kuò)展規(guī)劃中的應(yīng)用.天津大學(xué)學(xué)報(bào),2010,43(9):798–803.Wang J R,Yu Y X,Zeng Y.Discrete monkey algorithm and its application in transmission network expansion planning.Journal of Tianjin University,2010,43(9):798–803.(in Chinese)

[30]張佳佳,張亞平,孫濟(jì)洲.基于猴群算法的入侵檢測(cè)技術(shù).計(jì)算機(jī)工程,2011,37(14):131–133.Zhang J J,Zhang Y P,Sun J Z.Intrusion detection technology based on monkey algorithm.Computer Engineering,2011,37(14):131–133.(in Chinese)

[31]Yi T H,Li H N,Zhang X D.A modified monkey algorithm for optimal sensor placement in structural health monitoring.Smart Materials and Structures,2012,21(10):105033.

[32]伊廷華,張旭東,李宏男.基于改進(jìn)猴群算法的傳感器優(yōu)化布置方法研究.計(jì)算力學(xué)學(xué)報(bào),2013,30(2):218–223.Yi T H,Zhang X D,Li H N.Optimal sensor placement based on monkey algorithm.Chinese Journal of Computational Mechanics,2013,30(2):218–223.(in Chinese)

[33]Kiefer J,Wolfowitz J.Stochastic estimation of the maximum of a regression function.The Annals of Mathematical Statistics,1952,23(3):462–466.

[34]Spall J C.An overview of the simultaneous perturbation method for efficient optimization.Johns Hopkins APL Technical Digest,1998,19(4):482–492.