基于全局靈敏度分析的改進微分進化算法

何小龍, 白俊強, 李宇飛

(西北工業大學 航空學院, 陜西 西安　710072)

?

基于全局靈敏度分析的改進微分進化算法

何小龍, 白俊強, 李宇飛

(西北工業大學 航空學院, 陜西 西安710072)

摘要：為了提高微分進化算法在問題維度較高、計算量有限的情況下的尋優能力，提出了基于全局靈敏度分析方法的改進微分進化算法。使用Morris-One-at-a-Time (MOAT)全局靈敏度方法對典型函數進行靈敏度分析，與另一種全局靈敏度分析方法Sobol方法進行了對比，證明MOAT方法計算效率和精度較高。基于MOAT方法分別構造了考慮靈敏度信息的改進交叉算子和變異算子，從而得到2種改進微分進化算法GSADE1和GSADE2。使用5個顯式50維測試函數對2種改進算法進行了測試，并與基準微分進化算法進行了對比，發現2種改進算法在收斂速度和魯棒性方面都有所提高。

關鍵詞：算法; 全局優化; 微分進化; 計算機仿真; 靈敏度分析

微分進化算法(differential evolution, DE)與進化規劃、遺傳規劃、遺傳算法類似，都使用了優勝劣汰、適者生存的思想，是近年來發展的一種性能突出的進化算法(evolutionary algorithm)。DE由Rainer Stom 和 Kenneth Price于1995年提出[1]，實踐發現其收斂速度快、可調參數少、魯棒性好、算法簡單，在計算機、控制工程、化工、電力等眾多學科中得到了廣泛應用。國內外有許多學者進行過研究，蘇海軍、劉波、Swagatam Das等人[2-4]分別對DE算法研究進行了綜述。

但是DE還存在一些問題，一方面是收斂性問題，DE也存在著早熟和局部搜索能力弱的缺點[5]，而且在問題維度較高時更為嚴重，需要極高的計算量來保證優化效果；另一方面是如何針對特定問題設置合理的參數來保證良好的算法性能。為了繼續提高DE的性能，有許多學者致力于對DE進行改進，一種思路是對已有操作算子進行改進，如Fan等人提出的三角法變異算子[6]，另一種思路是提出新的操作算子，如Wang等提出的加速和遷移操作[7]研究。

本文針對以上兩方面問題，提出將全局靈敏度分析方法(global sensitivity analysis,GSA)與DE相結合的改進思路。首先，計算靈敏度的過程就是從問題中提取更多信息的過程，用靈敏度信息來和優化算法結合，有助于算法捕捉問題的特征，從而提高算法性能。其次，通過構造新的算子公式，使用靈敏度信息對算法參數進行控制，一定程度上解決了算法參數設置的問題。

1基本微分進化算法

1.1變異算子

(1)

在種群進化過程中,如果已經靠近最優解,那么個體之間的差分值會逐漸減小,那么這種擾動就會自動變小。

放縮因子F較小會使得算法變異能力不足、對優化空間的探索能力不足,從而導致容易陷入局部最優;較大的F提高了算法搜索到全局最優的概率,但是當F較大時,算法的收斂速度會明顯降低,這是因為當差分向量乘以F后會使變異個體距離父代較遠,對父代個體的繼承性較差,收斂速度會明顯下降。根據一些學者針對測試函數的使用經驗,放縮因子的選取范圍為0.5～0.9。

Price和Storn對基本DE的變異公式進行微調,提出了幾種不同的變異操作算子,具體細節的不同可用符號DE/X/Y/Z來區分。Price曾提出DE/BEST/2/BIN是一種比較值得研究的算子,本文采用這種算子作為基礎進行改進。

1.2交叉算子

交叉運算的目的是通過變異個體和第n代個體的隨機重組提高種群的多樣性。標準微分進化算法的交叉公式如下

(2)

1.3選擇算子

(3)

2MOAT全局靈敏度分析方法

2.1靈敏度分析方法

靈敏度分析方法的主要研究內容是:定性或定量地分析輸入變量(或參數)的變化如何影響輸出值,簡單說來就是分析輸入參數對輸出值影響的程度大小[8]。靈敏度分析的對象可以是任意一個系統或模型,只需要輸入和輸出數據即可,本質上是通過系統的數學性質對系統的內在特性進行研究。靈敏度分析方法有很多種,各種方法出于不同的目的而使用了不同的構造思想和數學工具,導致這些方法各具特色。對特定的問題而言,可根據其適用條件和使用要求選取最合適的方法。

靈敏度分析方法大體上可分為局部靈敏度分析方法、全局靈敏度分析方法和基于人工神經網絡的方法這幾大類[9]。局部靈敏度分析方法把導數信息、梯度作為靈敏度信息,這種信息是局限于某個基準點附近的信息,包括手工求導、有限差分、自動微分、復變量方法等幾種。全局靈敏度分析方法是在整個定義域內,研究設計變量對系統輸出的影響,又分為篩選法、蒙特卡羅方法、基于方差的方法等幾類。篩選法通常用來處理含有大量輸入變量的模型,包括OAT方法、COTTER方法、IFFD方法等,計算量相對比較小,MOAT[10]是其中較優秀的一種。蒙特卡羅方法包括散點圖法、相關系數法、回歸分析法,大多是使用統計學的方法,對于非線性、非單調的問題適用性較差。方差法包括重要性估計法、傅里葉幅值靈敏度測試法(FAST)等,是近年來研究發展最為迅速、應用最為廣泛的一類方法。其中比較特殊的是Sobol方法[11],具有蒙特卡羅方法和方差法兩方面的特點。關于這些方法的性能對比和分析,在若干文獻[12]里面都可以找到其他學者的研究結果。靈敏度分析已廣泛應用于大氣化學、氣體排放、魚群動力學、綜合指標分析、金融投資、放射性廢物處理、地質信息系統、固體物理學等等方向。

近年來全局靈敏度分析方法的發展和應用遠遠超過了局部靈敏度分析方法,而基于神經網絡的方法較新但研究較少。其中MOAT、Sobol等幾種方法的研究最為廣泛,MOAT方法的特點是適用于較高維度問題,計算量遠遠低于Sobol,但精度稍低于Sobol方法。

2.2MOAT靈敏度分析方法

Morris于1991年提出了一種One-at-a-time(OAT)方法,以下簡稱為MOAT方法,在GSA中是一種相對高效的方法。OAT方法的最主要特點是每次分析一個分量對輸出的影響,如果需要分析所有分量的影響,那么需要多次計算。MOAT實際上是對各種局部靈敏度分析方法的拓展,主要是在可行域內隨機取值,從而在一定程度上消除了對初始點的依賴,最后通過取平均值體現了全局性。MOAT方法可以用于識別具有多個輸入變量的模型中的少數關鍵變量,一些學者把有這種作用的方法稱為篩選法。

MOAT首先定義了Elementary Effect(EE),例如對于設計變量的第i個分量

EEi(x)=

EE的定義式類似于偏導數,但是由于Δ的取值并非接近于無窮小,因此EE不同于偏導數。MOAT方法中,求解EE時的基準點x可在整個定義域內隨機選取,因此在一定程度上擺脫了對于x的依賴。假設設計變量x=[x1,x2,…xk],共有k個分量。如果計算所有分量的影響,那么至少需要取k+1個點x1,x2…xk+1,最簡單情況就是xi+1-xi=Δ·ei,ei=[0,0…1…0]是第i個分量為1的單位向量,Δ是MOAT方法選取的一個參數。

(4)

那么計算k+1個函數值,兩兩相減然后除以Δ就可求得k個EE。Δ是影響EE的一個重要參數。MOAT把設計變量的定義域離散化來求得Δ。例如分量xi的初始定義域是實數域,而現在使用MOAT只從{0,1/(p-1),2/(p-1),…,1}這p個值中隨機選取某一個作為xi。這樣一來,整個k維的定義域空間Ω被劃分成k×p個網格點,x只能從這些網格點中隨機選取。Δ為1/(p-1)的倍數。而且需要保證,當x為區域Ω中的任意值,點x+Δ仍在Ω內。MOAT的總計算量為r(k+1),每k+1個點組成一個路徑(path),重復r次然后取平均即得到靈敏度。

MOAT方法的具體計算步驟如下:

實際使用中需要考慮基準點在定義域內的隨機選取。隨機化處理的步驟如下:

1) 令D*是k維對角矩陣,對角元素是1或者-1,概率各是0.5。設Jm,k是m×k的全1矩陣,得到(1/2)[(2B-Jm,k)D*+Jm,k]。

2) 設x*是隨機選擇的基準值,每一個分量都是從{0,1/(p-1),2/(p-1),…,1-Δ}中隨機等概率地選擇一個值。

3) 設P*是一個k維隨機擾動矩陣,每一列包含一個1,其他所有元素都是0,而且任意2列的1不在同一行。

其中D*、x*、P*是各自獨立地隨機生成。最后就得到了

(5)

2.3MOAT方法與Sobol方法對比分析

本小節對MOAT與Sobol 2種方法進行對比,使用這2種方法對測試函數進行靈敏度計算并對比二者在計算效率、計算精度方面的問題。由于MOAT方法和Sobol方法的原理不同,這2種方法使用了2種不同的數學方式來度量靈敏度信息。但是Campolongo發現,MOAT方法定義的Elementary Effects可以作為對Sobol方法定義的Sti(total sensitivity indice)一個較好的近似,理想情況下二者會逼近于同樣一個期望,因此具有一定的可比性[13]。

此處給出的是Sobol gstar function的測試結果對比,這個函數是Sobol和很多學者在靈敏度研究時經常使用的測試函數中的一個。gstar函數中的系數采用a=[0,0.1,0.2,0.3,0.4,0.8,1,2, 3,4]。

由于MOAT和Sobol方法都依賴于樣本進行計算,因此樣本數量的影響十分重要。在具體操作過程中,應逐漸提高樣本數量直至計算收斂,結果才是可信的。MOAT設置參數為:level=4,重復10次,總計110個樣本。Sobol方法使用拉丁超立方取樣,樣本總數1 200個。計算結果如圖1～圖2所示,橫坐標是設計變量的序號,縱坐標分別是MOAT和Sobol定義的靈敏度值。

圖1　gstar函數使用MOAT計算靈敏度結果

圖2　gstar函數使用Sobol計算靈敏度結果

從這2種方法的計算結果可以看出,MOAT的Elementary Effects和Sobol的Total Sensitivity Indice有近似相同的趨勢,但是具體數值有所差異。g函數的結果中對x2和x4的靈敏度計算差別較大。gstar函數的結果MOAT不如Sobol的結果精確。從參數設置來看,x1至x10的權重系數依次遞減,那么靈敏度應該也是依次遞減的,Sobol正確體現出了這種趨勢。MOAT由于樣本生成方法的影響,在每個變量的區間中取定若干點,而不是完全隨機取值,這樣在靈敏度計算時,設計變量的每一維的取值比較固定。可能就是這種原因導致了MOAT的靈敏度精度比Sobol稍差一些。另外,圖中顯示的圓圈表示求解靈敏度方差。

測試算例表明:MOAT方法的理論基礎弱于定量分析,但優點是需要樣本數量少。Sobol方法理論較嚴密,但計算量需求極大。考慮到本研究的應用背景是針對計算量受到限制的問題,Sobol方法實用性不強,而且對于改進DE而言,對靈敏度的精度要求并非很高,因此選擇MOAT方法來展開后續研究。

3改進微分進化算法

針對某些實際工程問題,需要在給定計算量約束的情況下,盡量獲取更好的優化結果。這也就限制了種群的規模和迭代次數,對優化算法的性能提出了更苛刻的要求。

本文研究提出了2種改進思路,分別對微分進化算法的交叉和變異公式進行修改,主要考慮靈敏度信息是對設計變量的重要程度的一種度量,靈敏度值越大,表明對輸出的影響越大,那么在算法迭代過程中,可以考慮將靈敏度信息耦合到變異算子或者交叉算子中去。本質上,這是通過靈敏度信息來控制交叉或者變異算子,其思路與自適應DE的一些方法類似。

第一種改進方式是修改交叉算子,改進后的DE以下簡稱為GSADE1。根據靈敏度計算的結果,可以據此修改每個分量的交叉概率。具體操作方式為:將CR擴展為矢量,CR=[CR1,…CRD],其維度和設計變量的維度相同,CR的每個分量的大小使用靈敏度信息來控制。設靈敏度矢量S=[S1,…,Si,…SD],其中靈敏度的最小值和最大值分別是Smin=min(S)、Smax=max(S),那么可以得到CR的計算公式如下

(6)

式中，α和β相加等于1,例如α=0.4、β=0.6,另外假設靈敏度的值是大于零的,那么CRi最終是在[0.6,1]之間變化,即式中的Si首先被歸一化,然后被縮放到指定的區間內。最后,微分進化算法的交叉公式如下,基本形式保持不變。考慮到MOAT方法計算的精度不是很高的問題,通過設置參數能保證交叉算子的值位于一定區間內,即使靈敏度精度不足,也能保證一定量的交叉概率,減弱了因靈敏度計算精度不足的影響。改進變異算子的思路也是同理。

(7)

第二種改進方式是修改變異算子,改進后的DE以下簡稱為GSADE2。本文使用的變異算子形式如下,從中可以看出主要影響因素就是個體的選取和放縮因子F的影響。對F的修改方式與CR類似,將F擴展為矢量,F=[F1,…FD],其維度和設計變量的維度相同,F的每個分量的大小使用靈敏度信息來控制。設靈敏度矢量S=[S1,…,Si,…SD],其中靈敏度的最小值和最大值分別是Smin=min(S)、Smax=max(S),那么可以得到F的計算公式如下

(8)

例如λ=0.3、ω=0.7,另外優于靈敏度的值是大于零的,那么Fi最終是在[0.7,1]之間變化,即式中的Si被變換到指定的區間內。最終得到改進后的變異公式如下

(9)

4函數測試及分析

本節使用5個解析函數作為測試問題,采用基本DE和改進算法分別尋優,以驗證改進算法的性能。使用2005年IEEE進化計算大會的標準測試函數[14]中的中的ShiftedSchwefel′s、ShiftedRotatedHighConditionedEllipticFunction、ShiftedRosenbrock′sFunction、ShiftedRotatedRastrigin′sFunction、ShiftedRotatedExpandedScaffer′s,仿照該文獻中的編號依次稱為函數2、3、6、10、14,設計變量范圍和函數中的參數均使用文獻中所給出的值。其中,函數2是單峰值函數,函數3是高條件數的單峰值函數,函數6是多峰值函數,函數10和函數14具有旋轉性和極多局部峰值。這幾個函數包含了單峰值、多峰值、帶旋轉等多種特性,因此具有一定的代表性。這幾個函數的維度可以人為指定,考慮到本文面向的工程問題的維度在50維左右,因此針對50維的情況來進行測試。

各函數均使用MOAT方法計算靈敏度信息,樣本總數510個,如果收斂較快則減少樣本數目。

初始DE參數設置使用DE/BEST/2/BIN變異算子,放縮因子F=0.5,交叉因子CR=0.9。GSADE1參數設置是α=0.1、β=0.9、放縮因子F=0.5。GSADE2參數設置是λ=0.2、ω=0.5、交叉因子CR=0.9。除了以上區別,改進算法的其他參數都與初始DE相同。種群大小統一設定為50和100,限定進化50代。

對每個函數重復測試20次,將收斂歷史求和取平均,得到平均的收斂歷史。收斂的最終結果的20次平均值如表1所示。

表1　收斂平均值對比

圖3～圖7顯示了3個較復雜測試函數的最后20步平均收斂曲線。

種群數100的優化結果顯示:2種改進算法完全優于初始DE。考慮到計算靈敏度信息使用的樣本數量,GSADE1和GSADE2的收斂速度也比基本DE更快。在函數2、函數6、函數14中改進算法優勢明顯。而在函數3、函數10算例中,改進算法在優化前期優勢不明顯,但是在算法后期優于DE,最終GSADE1和GSADE2的收斂結果都優于DE。

在種群數50的情況下,GSADE1只在函數2、函數6、函數14的測試中優于DE。主要原因可能是在種群50嚴重不足地情況下,由于靈敏度對交叉概率的影響,導致在函數3、函數10的測試中后期開發能力弱化。雖然在早期搜索中收斂較快,但是在算法后期開發能力不足。GSADE2優于GSADE1,只有函數3的測試中表現差于DE。這個算例表明交叉算子和變異算子對開發能力的影響有所不同。

圖3　函數6收斂結果　　　　　　圖4　函數10收斂結果　　　　　　　　　　圖5　函數14收斂結果

5結論

本文提出2種基于全局靈敏度分析的改進DE算法GSADE1和GSADE2并進行了算例驗證,得出以下結論:

1)MOAT方法計算效率高,精度稍低于Sobol方法,將之與優化算法結合對于改善優化搜索結果、提高優化效率具有積極作用。

2) 2種改進DE算法在種群數等于100的情況下完全優于初始DE,而在種群數50的情況下GSADE2優于GSADE1,但是仍有所不足,需要進一步針對極小種群的情況進行研究。

參考文獻：

[1]StornR,PriceK.DifferentialEvolution-aSimpleandEfficientHeuristicforGlobalOptimizationoverContinuousSpaces[J].JournalofGlobalOptimization, 1997, 11(4): 341-359

[2]蘇海軍,楊煜普,王宇嘉. 微分進化算法的研究綜述[J]. 系統工程與電子技術, 2008, 30(9): 1793-1797

SuHaijun,YangYupu,WangYujia.ResearchonDifferentialEvolutionAlgorithm:aSurvey[J].SystemsEngineeringandElectronic, 2008, 30(9): 1793-1797 (inChinese)

[3]劉波,王凌,金以慧. 差分進化算法研究進展[J]. 控制與決策, 2007, 22(7): 721-729

LiuBo,WangLing,JinYihui.AdvancesinDifferentialEvolution[J].ControlandDecision, 2007, 22(7): 721-729 (inChinese)

[4]DasS,SuganthanPN.DifferentialEvolution:aSurveyoftheState-of-the-Art[J].IEEETransonEvolutionaryComputation, 2011, 15(1): 4-31

[5]趙光權. 基于貪婪策略的微分進化算法及其應用研究[D]. 哈爾濱：哈爾濱工業大學, 2007

ZhaoGuangquan.DifferentialEvolutionAlgorithmwithGreedyStrategyandItsApplications[D].Harbin,HarbinInstituteofTechnology, 2007 (inChinese)

[6]FanHY,LampinenJ.ATrigonometricMutationOperationtoDifferentialEvolution[J].JournalofGlobalOptimization, 2003, 27(1): 105-129

[7]WangFS,JingCH,TsaoGT.Fuzzy-Decision-MakingProblemsofFuelEthanolProductionUsingaGeneticallyEngineeredYeast[J].Industrial&EngineeringChemistryResearch, 1998, 37(8): 3434-3443

[8]SaltelliA,RattoM,AndresT,etal.GlobalSensitivityAnalysis[M].ThePrimer,JohnWileyandSons, 2008

[9]朱亞濤. 基于敏度分析的飛行器氣動/隱身綜合優化設計策略研究[D]. 上海：上海交通大學, 2010

ZhuYatao.TheStrategyInvestigationonAircraftAerodynamic——StealthCompositeOptimalDesignBasedonSensitivityAnalysis[D].Shanghai,ShanghaiJiaotongUniversity, 2010 (inChinese)

[10]MorrisMD.FactorialSamplingPlansforPreliminaryComputationalExperiments[J].Technometrics, 1991, 33: 161-174

[11]Sobol′IM.SensitivityEstimatesforNonlinearMathematicalModels[J].MatModel, 1990, 2(1): 112-118

[12]GanY,DuanQ,GongW,etal.AComprehensiveEvaluationofVariousSensitivityAnalysisMethods:ACaseStudywithaHydrologicalModel[J].EnvironmentalModelling&Software, 2014, 51: 269-285

[13]CampolongoF,CariboniJ,SaltelliA.AnEffectiveScreeningDesignforSensitivityAnalysisofLargeModels[J].EnvironmentalModelling&Software, 2007, 22(10): 1509-1518

[14]SuganthanPN,HansenN,LiangJJ,etal.ProblemDefinitionsandEvaluationCriteriafortheCEC2005SpecialSessiononReal-ParameterOptimization[J].NanyangTechnologicalUniversity, 2005(1): 341-357

收稿日期：2015-10-22

基金項目：國家“973”計劃(2014CB744804)資助

作者簡介：何小龍(1989—)，西北工業大學博士研究生，主要從事飛行器氣動外形優化設計及優化方法研究。

中圖分類號：TP301.6

文獻標志碼：A

文章編號:1000-2758(2016)03-0411-07

A Global Sensitivity Analysis Enhanced Differential Evolution Algorithm

He Xiaolong, Bai Junqiang, Li Yufei

(College of Aeronautics, Northwestern Polytechnical University, Xi′an 710072, China)

Abstract:An improved differential evolution (DE) algorithm based on global sensitivity analysis is proposed to enhance performance in high dimension problems with limited computation resources. Morris-One-at-a-Time(MOAT) method was firstly tested on a typical function and compared with Sobol sensitivity method, showing high efficiency and acceptable result. Then MOAT method is used to calculate sensitivity for each dimension of input vector, and new crossover and mutation operators are proposed to incorporate sensitivity for two improved algorithm GSADE1 and GSADE2. Five 50-dimension functions were used for test, showing both two new algorithms are better than the original DE.

Keywords:algorithm; global optimization; differential evolution; computer simulation; sensitivity analysis