求解黑箱全局優(yōu)化問題的一種改進(jìn)的隨機(jī)算法

周 喆

(重慶師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，重慶 401331)

求解黑箱全局優(yōu)化問題的一種改進(jìn)的隨機(jī)算法

周 喆

(重慶師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，重慶 401331)

提出了一個(gè)改進(jìn)的隨機(jī)算法求解黑箱全局優(yōu)化問題，算法中使用的徑向基函數(shù)插值是昂貴的目標(biāo)函數(shù)的近似.算法在可靠的信賴域內(nèi)通過簡單而有效的方法選擇下一個(gè)迭代點(diǎn).在早期的迭代過程中,算法傾向于在響應(yīng)面模型的全局最優(yōu)值附近改善響應(yīng)面模型的局部逼近精度.它可能暫時(shí)陷入局部最優(yōu),但它有能力在之后的迭代中探索可行域中的其他區(qū)域并改善響應(yīng)面的全局逼近精度.數(shù)值實(shí)驗(yàn)結(jié)果說明了該算法的有效性.

全局優(yōu)化; 昂貴的目標(biāo)函數(shù); 響應(yīng)面模型; 徑向基函數(shù)插值

本文主要考慮如下形式的全局優(yōu)化問題：

minf(x)
s.t.a≤x≤b,x∈Rd

其中:f(x)∶Rd→R是一個(gè)昂貴的黑箱函數(shù)，也就是說它的函數(shù)估值代價(jià)非常大，a,b∈Rd.一般來說目標(biāo)函數(shù)f(x)的導(dǎo)數(shù)信息是不可用的.令D=[a,b]表示上述問題的可行域.全局優(yōu)化算法的任務(wù)是在可行域D中找到一個(gè)全局最優(yōu)點(diǎn)x*使得對(duì)所有的x∈D都有f(x)≤f(x*).在實(shí)際中通過全局優(yōu)化算法找到問題的精確的全局最優(yōu)點(diǎn)是很困難的，所以一般只要求找到全局最優(yōu)點(diǎn)x*的一個(gè)逼近點(diǎn).

響應(yīng)面方法作為一種流行的方法經(jīng)常被用來求解此類優(yōu)化問題.一般，響應(yīng)面方法通過建立響應(yīng)面模型來代替原有的模型來參與優(yōu)化過程，響應(yīng)面模型是真實(shí)的目標(biāo)函數(shù)的比較廉價(jià)的近似.早期的響應(yīng)面方法[1-2]通常采用低階的多項(xiàng)式回歸以及多變量插值并結(jié)合信賴域方法[3-4].后來，基于Kriging[5-6]、徑向基函數(shù)(RBF)[7-9]的全局優(yōu)化方法也逐漸被研究.此外，一些通過在一組隨機(jī)生成的候點(diǎn)中選擇下一個(gè)迭代點(diǎn)的隨機(jī)響應(yīng)面(SRS)方法[10-12]也被提出.對(duì)于求解黑箱全局優(yōu)化問題，這些隨機(jī)響應(yīng)面(SRS)方法和其他流行的算法相比較顯得更加有效.

本文基于徑向基函數(shù)插值模型提出了一種改進(jìn)的隨機(jī)響應(yīng)面方法，該方法的框架分為兩個(gè)階段：第一階段，通過檢測連續(xù)的兩次響應(yīng)面模型的全局最優(yōu)值之間的歐式距離來確定是否選擇當(dāng)前的響應(yīng)面模型的全局最優(yōu)點(diǎn)作為下一個(gè)迭代點(diǎn).如果沒有執(zhí)行這一步驟，則將下一個(gè)迭代點(diǎn)的探索區(qū)域限制在當(dāng)前的響應(yīng)面模型的全局最優(yōu)點(diǎn)的附近，并在該區(qū)域內(nèi)從一簇隨機(jī)生成且服從正態(tài)分布的候選點(diǎn)中選擇下一個(gè)迭代點(diǎn)；第二階段，當(dāng)連續(xù)多次迭代中得到的迭代點(diǎn)所對(duì)應(yīng)的目標(biāo)函數(shù)值沒有產(chǎn)生實(shí)質(zhì)性的改進(jìn)的次數(shù)超過自定義的閾值時(shí)，則將探索區(qū)域擴(kuò)充到整個(gè)可行域并在一簇均勻分布的隨機(jī)候選點(diǎn)中選擇下一個(gè)迭代點(diǎn).另外，在這兩個(gè)階段中將采用一些基于極大極小距離的策略來保證選擇的下一個(gè)迭代點(diǎn)不會(huì)太靠近已有的取樣點(diǎn).

1 預(yù)備知識(shí)

由于模型構(gòu)造的簡易性以及在全局逼近上的可靠性，徑向基函數(shù)插值模型將被用來作為算法中的響應(yīng)面模型.假設(shè)給定一個(gè)點(diǎn)集An={x1,x2,…,xn}?D以及相應(yīng)的目標(biāo)函數(shù)值f(x1),f(x2),…,f(xn).徑向基函數(shù)插值模型被定義成如下形式[7]：

(1)

表1 徑向基函數(shù)的常見形式

未知的參數(shù)λ,c能夠通過求解如下線性方程組:

(2)

來得到，其中:

Φij=φ(‖x-xi‖),i,j=1,2,…,n,

F=(f(x1),f(x2),…,f(xn))T.

2 算法框架

接下來提出了一個(gè)基于響應(yīng)面模型的改進(jìn)的隨機(jī)算法ISARS，該算法提供了非常有效的迭代搜索策略并能夠保證局部搜索和全局搜索之間的平衡.在下面的符號(hào)中，n表示已有的取樣點(diǎn)的數(shù)量，An={x1,x2,…,xn}表示已有的n個(gè)取樣點(diǎn)的集合，離散的數(shù)據(jù)點(diǎn)集Bn={(x1,f(x1)),(x2,f(x2)),…,(xn,f(xn))}被用來構(gòu)造響應(yīng)面模型sn(x).

ISARS算法：

輸入：

1)一個(gè)確定的定義在閉超矩形D=[a,b]?Rd上的目標(biāo)函數(shù)f.

2)一個(gè)特別的響應(yīng)面模型，比如帶有尾部多項(xiàng)式的Cubic徑向基函數(shù)插值模型.

3)一個(gè)采用對(duì)稱拉丁超立方設(shè)計(jì)SLHD[15]生成的初始取樣點(diǎn)集An0={x1,x2,…,xn0}?D.

6)允許進(jìn)行的最多的函數(shù)估值次數(shù)Nmax.

7)每一次迭代中隨機(jī)生成的試驗(yàn)點(diǎn)的數(shù)量k.(默認(rèn)值：k=min(500 d,5 000).)

8)連續(xù)失敗的迭代次數(shù)的閾值Tfail.(默認(rèn)值：Tfail=min(5d+1,20).)

輸出：算法停止迭代后計(jì)算出的最優(yōu)解.

Step3.當(dāng)不滿足迭代終止條件(比如：當(dāng)nlt;Nmax)時(shí)，執(zhí)行

Step3.3 (選擇下一個(gè)迭代點(diǎn)) 如果Cfaillt;Tfail，那么

Step3.3(b)(i) (在整個(gè)可行域中隨機(jī)生成試驗(yàn)點(diǎn)) 隨機(jī)生成一組試驗(yàn)點(diǎn)集Ωn={yn,1,yn,2,…,yn,k}.對(duì)于j=1,2,…,k，yn,j=a+(b-a)?qj，其中qj∈Rd是一個(gè)隨機(jī)生成的各分量在(0,1)上服從均勻分布的向量，?表示兩個(gè)向量對(duì)應(yīng)的分量相乘的運(yùn)算.

結(jié)束

Step3.4 (進(jìn)行耗時(shí)的函數(shù)計(jì)算) 計(jì)算迭代點(diǎn)xn+1處的目標(biāo)函數(shù)f(xn+1).

結(jié)束

3 算法分析

接下來分析算法中比較關(guān)鍵的一些步驟.在ISARS算法中，在擬合或者更新響應(yīng)面模型(Step3.1)之前通常需要檢測是否存在極端的目標(biāo)函數(shù)值，因?yàn)闃O大或者極小的目標(biāo)函數(shù)值會(huì)導(dǎo)致響面模型的波動(dòng)劇烈而影響逼近精度.為了避免這一情況，當(dāng)maxf(xi)-minf(xi)gt;2·103,xi∈An時(shí)，采用文獻(xiàn)[9]中的對(duì)數(shù)轉(zhuǎn)換函數(shù)對(duì)所有的目標(biāo)函數(shù)值進(jìn)行穩(wěn)定化處理.φ(f)的函數(shù)表達(dá)式如下：

在Step3.3(a)(v)及Step3.3(b)(iii)中，從試驗(yàn)點(diǎn)集中選擇下一個(gè)迭代點(diǎn)的過程采用文獻(xiàn)[11]中的LMSRS方法和文獻(xiàn)[12]中的DYCORS方法，但是ISARS算法預(yù)先對(duì)試驗(yàn)點(diǎn)進(jìn)行了一次篩選排除了部分試驗(yàn)點(diǎn)(見Step3.3(a)(iv)以及Step3.3(b)(ii)).比如在二維平面上，以已有的取樣點(diǎn)為中心的圓形之內(nèi)的取樣點(diǎn)都被排除掉.在進(jìn)行全局搜索的情況下，當(dāng)ρn比較小時(shí)，下一個(gè)迭代點(diǎn)更傾向于在已有的取樣點(diǎn)附近來選擇以提高局部逼近精度；而當(dāng)ρn比較大時(shí)，選取的下一個(gè)迭代點(diǎn)更多地被限制在未完全探索的區(qū)域來提高全局逼近精度.

在主要的迭代過程中，ISARS算法通過閾值Tfail分為兩個(gè)階段.第一階段，試驗(yàn)點(diǎn)在當(dāng)前響應(yīng)面模型的全局最優(yōu)點(diǎn)附近隨機(jī)生成并服從正態(tài)分布；第二階段，試驗(yàn)點(diǎn)在整個(gè)可行域中隨機(jī)生成并服從均勻分布.

此外，當(dāng)算法陷入局部最優(yōu)時(shí)，重啟動(dòng)策略作為一種有效地改進(jìn)全局優(yōu)化方法的策略很容易應(yīng)用到ISARS算法中，在算法過程中僅需要設(shè)置一個(gè)重啟動(dòng)標(biāo)志Rflag并初始化Rflag=0，當(dāng)計(jì)數(shù)器Cfailgt;Tfail時(shí)，重啟動(dòng)標(biāo)志重置為Rflag=1，這時(shí)ISARS算法將拋棄之前所有的數(shù)據(jù)信息重新開始執(zhí)行算法的整個(gè)流程.

4 數(shù)值實(shí)驗(yàn)

本節(jié)將使用ISARS算法來求解維度在2到6之間的Dixon和Szeg?[16]測試問題并將它與文獻(xiàn)[11]中的LMSRS方法和文獻(xiàn)[12]中的DYCORS方法的數(shù)值結(jié)果進(jìn)行比較.這些問題的一些信息包括名稱、維數(shù)、定義域、局部極小點(diǎn)的個(gè)數(shù)、全局最優(yōu)點(diǎn)的個(gè)數(shù)以及全局最優(yōu)值在表2中給出.其中，Branin測試問題的所有的局部極小點(diǎn)都是全局最優(yōu)點(diǎn)，而對(duì)于其他的測試問題都只有一個(gè)全局最優(yōu)點(diǎn).

在進(jìn)行比較的三種算法中，采用的響應(yīng)面模型均是帶有尾部多項(xiàng)式的Cubic徑向基函數(shù)插值模型，并將采用了RBF插值模型的ISARS，LMSRS，DYCORS算法分別記作ISARBF，LMSRBF，DYCORBF算法.此外，結(jié)合了重啟動(dòng)策略的三種算法分別記為ISARBF-R，LMSRBF-R，DYCORBF-R算法.

表2 測試問題及其屬性

對(duì)于同一個(gè)問題，每一種算法都運(yùn)行30次.為了保證算法的公平性，對(duì)于求解同一個(gè)問題，不同的算法中初始取樣點(diǎn)的位置保證是相同的.另外，LMSRBF算法和DYCORS算法中的參數(shù)都保持和原文獻(xiàn)[11-12]中的一樣.在ISARS算法中，設(shè)置RBF準(zhǔn)則的權(quán)重模式為γ=〈0.02，0.25，0.5，0.95〉.

定義相對(duì)誤差E為：

其中fglobal≠0是函數(shù)f的全局最優(yōu)值.所有測試問題中，二維的問題允許進(jìn)行的函數(shù)估值次數(shù)為200次，其他的問題允許進(jìn)行的函數(shù)估值次數(shù)為500次.

表3中顯示了各個(gè)算法對(duì)于求解同一問題在30次運(yùn)行計(jì)算之后得到的最優(yōu)解關(guān)于一些指標(biāo)的統(tǒng)計(jì)結(jié)果.這些指標(biāo)包括最優(yōu)解中最好的值fglobal，最差的值fWorst，中值fMedian，平均值fmean，標(biāo)準(zhǔn)誤差RMSE以及計(jì)算終止后找到的最優(yōu)解與問題的最優(yōu)解之間的相對(duì)誤差為E≤1%下的成功率δSucc.

表3 算法運(yùn)行30次之后得到的最優(yōu)目標(biāo)函數(shù)值關(guān)于各指標(biāo)的統(tǒng)計(jì)結(jié)果

續(xù)表3

測試問題指標(biāo)ISARBFISARBF-RLMSRBFLMSRBF-RDYCORSDYCORS-RShekel7fBest10．400010．401610．399010．399010．402710．3958fWorst5．116410．33621．83703．72201．83735．0880fMedian10．357710．38573．722510．38653．245010．3845fmean9．489110．37955．23859．59144．65629．6817RMSE1．95280．01693．20922．04652．98101．7908δSucc24/3030/308/3026/306/3026/30Shekel10fBest10．530810．531610．525210．532510．536310．5338fWorst5．164110．46601．67643．83322．421710．5035fMedian10．489310．51784．501410．52283．835410．5236fmean9．954210．51306．113710．11905．758210．5232RMSE1．59610．01913．69171．51143．48700．0075δSucc26/3030/3012/3028/3010/3030/30Hartman6fBest3．32213．32233．32233．32233．32243．3223fWorst3．10273．31963．17043．32163．20313．2015fMedian3．32003．32163．32233．32223．32233．3223fmean3．28883．32153．27943．32223．28663．3182RMSE0．06510．00060．06090．00010．05460．0217δSucc24/3030/3020/3030/3021/3029/30

分析表3中的數(shù)據(jù)可知，對(duì)于Branin測試問題，幾種算法均找到了響應(yīng)面模型的全局最優(yōu)點(diǎn)，其中ISARBF-R算法以及DYCORS算法得到的最優(yōu)值的平均值均為0.397 9并且標(biāo)準(zhǔn)誤差均為0.000 1.對(duì)于Goldstein-Price測試問題，ISARBF-R算法的在所有的算法中得到的最優(yōu)值平均值最低，并且在沒有加入重啟動(dòng)策略的情況下，ISARBF算法在每一次運(yùn)行中都找到了相對(duì)誤差低于1%的全局最優(yōu)值的逼近解.對(duì)于Hartman3測試問題，DYCORBF-R算法表現(xiàn)得最好.

對(duì)于Shekel類的測試問題，在沒有加入重啟動(dòng)策略的情況下，ISARBF算法的表現(xiàn)要遠(yuǎn)優(yōu)于LMSRBF算法和DYCORBF算法，因?yàn)槠渌麅煞N算法在30次運(yùn)行計(jì)算中找到相對(duì)誤差低于1%的全局最優(yōu)值的逼近解的運(yùn)行次數(shù)比ISARBF算法少得多.另外，對(duì)于Shekel5和Shekel7測試問題，ISARBF-R算法在平均值，標(biāo)準(zhǔn)誤差以及成功率這些指標(biāo)上要優(yōu)于其他的算法.對(duì)于Shekel10測試問題，DYCORBF-R算法找到的全局最優(yōu)值更加精確并且在30次運(yùn)行計(jì)算中都找到了相對(duì)誤差低于1%的全局最優(yōu)值的逼近解.對(duì)于Hartman6測試問題，雖然LMSRBF-R算法求得的最優(yōu)解的標(biāo)準(zhǔn)誤差是最低的，但是比較最差的最優(yōu)解以及平均值可知，ISARBF-R算法的表現(xiàn)更好.

綜合以上的分析可知，在沒有加入重啟動(dòng)策略的情況下三種的算法中，ISARBF算法的表現(xiàn)要遠(yuǎn)遠(yuǎn)優(yōu)于LMSRBF算法以及DYCORBF算法.雖然加入重啟動(dòng)策略后各個(gè)算法都相應(yīng)的得到了改進(jìn)，但是改進(jìn)之后的算法在整體表現(xiàn)上來看還是ISARBF-R算法表現(xiàn)的更好.

[1]KHURIAI,CORNELLJA.Responsesurfaces[M].MDekker,1987.

[2]MYERSRH,MONTGOMERYDC.Responsesurfacemethodology:processandproductinoptimizationusingdesignedexperiments[J].Technometrics,1996,59(3):185-186.

[3]CONNAR,SCHEINBERGK,TOINTPL.Recentprogressinunconstrainednonlinearoptimizationwithoutderivatives[J].MathematicalProgramming,1997,79(1):397-414.

[4]POWELLMJD.Uobyqa:unconstrainedoptimizationbyquadraticapproximation[J].MathematicalProgramming,2002,92(3):555-582.

[5]HUANGD,ALLENTT,NOTZWI,etal.Globaloptimizationofstochasticblack-boxsystemsviasequentialkrigingmeta-models[J].JournalofGlobalOptimization,2006,34(3):441-466.

[6]JONESDR,SCHONLAUM,WELCHWJ.Efficientglobaloptimizationofexpensiveblack-boxfunctions[J].JournalofGlobalOptimization,1998,13(4):455-492.

[7]GUTMANNHM.Aradialbasisfunctionmethodforglobaloptimization[J].JournalofGlobalOptimization,2001,19(3):201-227.

[8] REGIS R G,SHOEMAKER C A.Constrained global optimization of expensive black box functions using radial basis functions[J].Journal of Global Optimization,2005,31(1):153-171.

[9] REGIS R G,SHOEMAKER C A.A quasi-multistart framework for global optimization of expensive functions using response surface models[J].Journal of Global Optimization,2013,56(4):1719-1753.

[10] REGIS R G.Stochastic radial basis function algorithms for large-scale optimization involving expensive black-box objective and constraint functions[J].Computers and Operations Research,2011,38(5):837-853.

[11] REGIS R G,SHOEMAKER C A.A stochastic radial basis function method for the global optimization of expensive functions[J].European Journal of Operational Research,2007,19(4):497-509.

[12] REGIS R G,SHOEMAKER C A.Combining radial basis function surrogates and dynamic coordinate search in high-dimensional expensive black-box optimization[J].Engineering Optimization,2013,45(5):1-27.

[13] BJRKMAN M,HOLMSTRM K.Global optimization of costly nonconvex functions using radial basis functions[J].Optimization and Engineering,2000,1(4):373-397.

[14] POWELL M J D.The theory of radial basis function approximation in 1990[C]//Light,W.(ed.),Advances in Numerical Analysis,Volume 2:Wavelets,Subdivision Algorithms and Radial Basis Functions,pp.105-210.Oxford University Press,1992.

[15] YE K Q,LI W,SUDJIANTO A.Algorithmic construction of optimal symmetric Latin hypercube designs[J].Journal of Statistical Planning and Inference,2000,90(1):145-159.

[16] DIXON L C W, G P.Towards global optimization 2[M].North-Holland Pub Co,1978.

責(zé)任編輯：時(shí)凌

AnImprovedStochasticAlgorithmforSolvingBlack-boxGlobalOptimizationProblems

ZHOU Zhe

(School of Mathematics Sciences，Chongqing Normal University，Chongqing 401331，China)

This paper proposes an improved stochastic algorithm for solving black-box global optimization by using radial basis function (RBF) interpolation which is an approximation of the expensive objective function.The algorithm selects the next iteration point by a simple and effective method in a reliable trust region.In the process of early iterations,algorithm tends to improve the local approximation accuracy around the global minimizer of response surface (RS) model.It may be temporarily trap in a local optimum,but it is capable of exploring other global region in latter iterations and improves the global approximation accuracy of response surface model.The numerical experiment results show the effectiveness of the algorithm.

global optimization;expensive objective functions;response surface model;radial basis function interpolation

2017-05-11.

重慶師范大學(xué)研究生科研創(chuàng)新項(xiàng)目(YKC16001).

周喆(1991-)，男，碩士，主要從事最優(yōu)化計(jì)算方法的研究.

1008-8423(2017)04-0403-06

10.13501/j.cnki.42-1569/n.2017.12.011

O241.6

A