劉毅,陳吉寧,杜鵬飛
(清華大學環境科學與工程系環境模擬與污染控制國家重點聯合實驗室,北京 100084) 摘要:模型參數優化是通過極小化目標函數使得模型輸出和實際觀測數據之間達到最佳的擬合程度. 由于環境模型本身的復雜性,常規優化算法難以達到參數空間上的全局最優. 近年來,隨著計算機運算效率的快速提高,直接優化方法得到了進一步開發與廣泛應用. 本文比較了CRS、SCE UA、SA 和Annealing2Simplex 等4 種算法應用于環境模型參數優化的結果和計算效率。
關鍵詞:參數優化;環境模型;CRS 算法;SCE UA 算法;Simulated2Annealing 算法;Annealing2Simplex 算法
中圖分類號:X11 文獻標識碼:A 文章編號:025023301 (2002) 0220620001 Comparison of Parameter Optimization Algorithms for Environmental Model
Yi Liu ,J ining Chen ,Pengfei Du ( Environmental Simulation and Pollution Control State Key Joint Laboratory ,Department of Environmental Science and Engineering, Tsinghua University , Beijing 100084 ) Abstract:Parameters identification is achieved through the minimization of objective function based on model outputs and the observed data. Because of ever increasing complexity of environmental2models, there are significant difficulty for conventional optimal methods to present a global optimization. On the contrast , however , direct optimization algorithms are widely developed in recent years due to increasing computer efficiency and show promising applications. Four direct optimal algorithms , i. e. CRS algorithm , SCE UA algorithm , SA algorithm and Annealing2Simplex algorithm , were thus selected in this paper to compare their performances via case studies.
Keywords:parameter optimization ; environmental model ; CRS algorithm; SCE UA algorithm; Simulated2Annealing algorithm; Annealing2Simplex algorithm 人們對環境系統的深入研究是建立在環境模型的廣泛應用基礎上的. 為了更加精確地刻畫環境系統的行為,環境模型在近10 年里表現出了強烈的復雜化趨勢;不同空間尺度、不同時間過程模型的耦合,進一步加劇了這一過程. 環境模型的復雜性導致了模型結構和參數可識別性問題的提出,并成為當今環境建模理論研究的熱點[1 ,2 ] . 其中在不確定性的框架下,模型參數的優化是研究的一個重要方面.
解決優化問題的難度主要取決于模型參數的空間維數和模型本身的非線性特征. 一般來說,參數越多、非線性越強,優化時間和精度就越差,同時也越不能夠保證優化算法是否收斂到整體最優. 傳統經驗表明,求解優化問題的困難主要體現為[4 ,5 ] : ①全局搜索可能收斂到多個不同的吸引域; ②每一個吸引域可能包含一個或多個局部最小值; ③目標函數在n 維參數空間上不連續; ④參數及相互間存在高度靈敏性和顯著非線性干擾; ⑤在最優解的附近,目標函數往往不具有凸性。
優化算法可以分為直接算法和間接算法2大類. 間接算法(如牛頓法以及各種以牛頓法為基礎的改進算法) 的局限性主要在于要求目標函數在相關值域上必須是可微的;而直接算法僅涉及目標函數值的計算,不需要計算目標函數的導數. 因此盡管后者的計算效率相對較低,但在環境問題的實際應用中,它可以有效和簡潔地解決由于模型復雜性所衍生的不可微函數的優化問題. 同時,與求解問題所耗費的時間相比,通常更為關注解的結果能夠在多大程度上描述系統的行為[1 , 2 ] . 所有這些使得直接算法在近年來得到了迅速發展和廣泛的應用.
本文的目的在于分析和比較幾種近年來漸為接受的參數直接優化算法的計算效率、計算精度及其算法穩定性. 由于直接算法本質上的隨機性,以及對于不同優化問題所表現出的算法特性上的差異,因此本文僅以經典測試函數和環境水文模型為例,對幾種算法進行詳細的比較研究. 1 參數優化算法 模型參數的優化即是尋求一組參數,使模型的輸出與實際觀測數據之間按給定目標函數的度量方式達到最佳擬合,即: Min f ( ys , yobv ,θ) = f ( ys , yobv ,θ3 ) θ ∈ S (1)
式中, f ( y , θ) 為目標函數; ys 表示模型輸出變量, yobv為系統觀測值;θ3表示參數可行域S 上的最優參數向量. 最早提出的直接算法均是簡單隨機方法,這些算法除了計算效率低外,主要缺陷在于要求目標函數在鄰域上必須是不相關的.
控制隨機搜索算法(CRS) [10 , 11 ] 有效地克服了簡單隨機算法的主要缺點,在計算過程中保存指定數目的參數樣本及其對應的目標函數值,引入幾何學中“重心”的概念,即考慮了新點產生的隨機性,又在一定程度上保證了搜索的整體性. 復合形混合演化算法( SCE UA) [4 , 5 ]所采用的競爭演化和復合形混合的概念繼承了CRS 算法中全局搜索和復合形演化的思想,該方法將生物自然演化過程引入到數值計算中,模擬了生物進化的過程,提高了計算效率和全局搜索整體最優的能力. 模擬退火算法(SA) [7 , 12 ]則假設優化問題的解及其目標函數分別與固體物質的微觀狀態及其能量所對應,采用蒙特卡洛(Monte Carlo) 隨機方法模擬固體穩定“退火”的過程,并假設優化過程中遞減目標函數值的控制參數t 與“退火”過程中的溫度T 所對應,對于控制參數t 的每一個值,算法持
續進行“產生新解2判斷2接受/ 舍棄”的迭代過程,應用該算法的關鍵在于確定合理的冷卻進度表. 退火單純形算法(AS) [8 , 9 ]綜合了下山單純形方法和模擬退火法2 種優化算法,充分利用單純形的形變信息,以一定溫度T 所對應的概率接受準則來指導單純形的映射、壓縮和擴展等過程,提高了計算效率和算法穩定性.
運用直接算法進行參數優化首先要正確選取算法的內部控制參數. 一般根據具體問題的特點,采取數值實驗的方法不斷測試參數“性能”,從而確定既保證算法收斂到滿意的極值,又不至于使計算過于耗費時間而導致算法失效. 選擇控制參數本身實際就是一個“優化”過程. 2 測試函數的優化 首先采用一個廣泛使用的、具有多個局部極小值的經典函數,對上述幾種優化算法的全局搜索能力進行較為全面的測試,其函數形式為:  盡管測試函數僅有x 、y 2 個變量,但其等值線卻具有復雜的空間結構,即多個局部極小、值域空間上的不連續和全局最小值周圍的香蕉狀狹長低谷(圖1) . 現已知參數全局最小值為(5 ,5) ,所對應的目標函數值為10.0。 
如前所述,為了保證直接算法的優化精度和避免算法失效,在算法的運行中需要在如下３個優化終止準則中進行選擇和組合: ①前m組最佳估計參數收斂到可接受的參數空間,即最優參數的不確定性小于ε1 ; ②最佳目標函數值的優化改進速率小于給定值ε2 ; ③搜索過程中參數采樣總數不超過給定的計算次數n.
表1 給出了4 種直接算法分別基于終止準則1 和2 ,以及不同初值假設條件下的測試函數參數最優值和計算效率. 為了盡可能消除算法隨機性對于算法比較的影響,表格中數據除了計算次數的偏差外,均為相同條件下5 次運算結果的算術平均值.
由于AS 算法中的下山單純形方法具有很強的方向性,因而在兩個終止準則約束下,算法均可迅速達到全局最優化,且結果具有較高精度. 盡管SCE UA 算法引入競爭演化思想后提高了樣本空間的搜索效率,但是不同種群從不同方向上向全局最優點逼近,以及同一種群內部仍然采用的單純形控制隨機搜索方法,使得該算法在計算效率上遠不及AS 算法. 然而SCEUA 算法在各種初值條件下均可得到十分精確的全局優化結果,可見該算法的穩定性和可靠性比AS 算法具有更大的優勢. 與前2 種算法相比,CRS 算法和SA 算法在計算效率和精度方面均表現不佳. 特別是當參數初值發生變化時,必須相應地調整控制參數才能得到較好的計算結果,這說明2 種算法的適應能力較差. 表1 4 種直接算法的參數優化結果1)
Table 1 Parameter optimization results of four direct optimal algorithms | 優化算法 | 編號 | 假設初值 | 參數優化結果 | 目標函數值 | 計算次數 | | 準則1 | 準則2 | | x | y | x | y | x | y | 準則1 | 準則2 | 準則1 | 準則2 | 偏差/ % | | CRS | 1 | 0.40 | 0.60 | 0.5000 | 0.5001 | 0.5000 | 0.5001 | 10.0001 | 10.0001 | 565 | 584 | -3.25 | | 2 | 0.10 | 0.20 | 0.5000 | 0.4999 | 0.5000 | 0.5000 | 10.0001 | 10.0000 | 806 | 841 | -4.16 | | 3 | 0.20 | 0.80 | 0.5001 | 0.5001 | 0.5000 | 0.5000 | 10.0002 | 10.0000 | 585 | 618 | -5.34 | | 4 | 0.70 | 0.70 | 0.4998 | 0.5001 | 0.5000 | 0.4999 | 10.0002 | 10.0000 | 541 | 571 | -5.25 | | 5 | 0.80 | 0.85 | 0.5000 | 0.5000 | 0.5000 | 0.5000 | 10.0000 | 10.0000 | 595 | 583 | 2.06 | | SCE UA | 1 | 0.40 | 0.60 | 0.5000 | 0.5000 | 0.5000 | 0.5000 | 10.0000 | 10.0000 | 479 | 356 | 34.55 | | 2 | 0.10 | 0.20 | 0.5000 | 0.5000 | 0.5000 | 0.5000 | 10.0000 | 10.0000 | 494 | 350 | 41.14 | | 3 | 0.20 | 0.80 | 0.5000 | 0.5000 | 0.5000 | 0.5000 | 10.0000 | 10.0000 | 501 | 359 | 39.55 | | 4 | 0.70 | 0.70 | 0.5000 | 0.5000 | 0.5000 | 0.5000 | 10.0000 | 10.0000 | 431 | 305 | 41.31 | | 5 | 0.80 | 0.85 | 0.5000 | 0.5000 | 0.5000 | 0.5000 | 10.0000 | 10.0000 | 480 | 335 | 43.28 | | SA | 1 | 0.40 | 0.60 | 0.5001 | 0.5002 | 0.5000 | 0.5000 | 10.0013 | 10.0000 | 230 | 349 | -34.10 | | 2 | 0.10 | 0.20 | 0.5041 | 0.4973 | 0.5000 | 0.5000 | 10.1300 | 10.0000 | 408 | 504 | -19.05 | | 3 | 0.20 | 0.80 | 0.4999 | 0.5002 | 0.5000 | 0.4999 | 10.0002 | 10.0000 | 468 | 579 | -19.17 | | 4 | 0.70 | 0.70 | 0.5002 | 0.4999 | 0.5000 | 0.4999 | 10.0005 | 10.0001 | 468 | 582 | -19.59 | | 5 | 0.80 | 0.85 | 0.5001 | 0.5000 | 0.5001 | 0.5000 | 10.0000 | 10.0000 | 490 | 544 | -9.93 | | AS | 1 | 0.40 | 0.60 | 0.4999 | 0.5001 | 0.4999 | 0.5001 | 10.0000 | 10.0000 | 73 | 63 | 15.87 | | 2 | 0.10 | 0.20 | 0.5000 | 0.5001 | 0.5000 | 0.5001 | 10.0001 | 10.0001 | 150 | 242 | -38.02 | | 3 | 0.20 | 0.80 | 0.5000 | 0.5000 | 0.5000 | 0.5000 | 10.0000 | 10.0000 | 108 | 135 | -20.00 | | 4 | 0.70 | 0.70 | 0.5001 | 0.4999 | 0.5001 | 0.4999 | 10.0000 | 10.0000 | 99 | 191 | -48.17 | | 5 | 0.80 | 0.85 | 0.4999 | 0.5001 | 0.5001 | 0.4999 | 10.0000 | 10.0001 | 157 | 157 | 0.00 |
1) 計算次數偏差的計算公式為( n準則1 - n準則2) / n準則2
表2 目標函數計算次數的統計特性
Table 2 Statistical characteristics of the calculating times of objective function
計算次數的統計量 | CRS | SCE- UA | SA | AS | 準則1的標準方差 | 106.89 | 27.36 | 106.65 | 35.46 | 的標準方差 | 114.05 | 22. 15 | 96.24 | 66.52 | 準則1和準則2偏差均值/% | -3.19 | 39.97 | –20.37 | –18.06 |
可以采用各種初始條件下對應同一終止準則計算次數的標準方差來表征算法的隨機性.盡管基于控制性隨機搜索思想構造的直接優化算法不可能完全消除計算結果的隨機性,但是算法結構的精心設計將在一定程度上減小隨機性影響,同時也就意味著算法穩定性的增大. 表2 中結果顯示SCE UA 算法的計算次數在各種初始條件下的標準方差遠小于其他3 種直接算法,再次表明了算法所具有的良好穩定性. 這一結果的產生是由于SCE UA 算法結構中引入的復合形概念及其競爭混合算法大大削弱了計算結果的隨機性.
對應準則1 和2 的計算次數偏差,表征了算法對于選用不同終止準則的敏感性. 根據表2 結果可以得到4 種算法的敏感性排序為(從大到小) : SCE UA > SA > AS > CRS ,特別是CRS 算法的計算次數偏差的最大值和統計平均值分別為- 5.34 %和- 3.19 %,遠小于其他3種算法. 更進一步發現,算法結構復雜性與算法對于終止準則選用的敏感性之間高度相關,表現為本例中算法敏感性排序與復雜性順序完全一致. 其原因在于,直接算法在結構復雜性所形成的內部約束與目標函數的外部約束的共同作用下趨向全局最優,其結果必然是一個不斷協調和折衷的過程. 算法越復雜,內部約束就越強,其計算結果就越發表現出與終止準則之間的高度相關性. 因此簡單算法對于外部終止準則控制性約束的適應能力較之復雜算法更為靈活. 3 環境模型參數優化:水文模型實例 隨著環境模型的不斷開發和廣泛應用,環境模型的種類和數量日益豐富,模型本身所表現出的結構特征也日趨復雜. 因此研究和比較直接算法的參數優化特性,既無可能也無必要將其應用于每一個環境模型的參數優化評估中. 現實而合理的方法是選擇具有典型結構復雜性和非線性特征的環境模型進行實例研究.由于篇幅所限,本研究僅以箱式水文模型為例,對4 種參數直接優化算法進行比較研究.
研究中采用2 個箱式的模型結構模擬徑流在流域中的運動情況(例如,Cooper 等人[3 ]在其研究中給出的水文箱式模型介紹),該模型共有11 個參數,表3 給出了本文中模型參數的取值范圍及其物理意義.
表3 水文模型參數取值及其物理意義
Table 3 The parameter value and its physical meaning of
| 模型參數 | 最大取值 | 最小取值 | 說明 | | k1 | 0.05 | 0.3 | 不同箱體出口處的徑流流速/ h-1 | | k2 | 0.008 | 0.1 | | k3 | 0.01 | 0.1 | | k4 | 0.005 | 0.1 | | k5 | 0.0005 | 0.002 | | h1 | 5.0 | 20.0 | 不同箱體出口處的高度/ mm | | h2 | 1.0 | 10.0 | | h3 | 700.0 | 1500.0 | | h4 | 500 | 1000.0 | | H1 | 0.1 | 20.0 | 上層和下層水箱高度/ mm | | H2 | 5.0 | 200.0 |
根據不同最優化準則構造的目標函數將對參數優化結果會產生一定影響[2 , 6 ] . 本研究以最小二乘法為優化準則構造的目標函數如下: 
其中, f obv為觀測流量, f sim表示模擬輸出流量; N為模擬點數. 本研究取N = 168 ,模擬暴雨期7 天的小時觀測流量. 終止準則采用準則1 和3.
圖2 是水文模型參數優化后的流量模擬曲線,可以直觀的看到4 條優化曲線均較好的擬合了觀測流量. 表4 給出了4 種直接算法分別在相同條件下運算5 次,以不確定性表示的最優參數值和最小二乘的目標函數值, 以及以CPU 時間表示的平均計算效率.
目標函數值和最佳估計參數的統計特征值共同表征了優化算法的穩定性和可靠性. 根據表4 中的計算結果, SA 算法和CRS 算法的目標函數值和最佳估計參數均具有較大不確定性,盡管2 種算法分別得到了相對較小的目標函數值( ZSA = 117774 , ZCRS = 118648) ,但是目標函數的均值和標準方差2 個統計量卻均較大. 相反,AS算法和SCE UA算法則表現出更為優秀的計算特性,即算法既能夠保證一定的計算精度,又具有可靠的重復性. 與SCE UA算法相比,AS 算法的目標函數值和最佳估計參數的不確定性更為縮小,進一步說明AS 算法全局參數尋優過程所具備的計算穩定性和可靠性. 
表4 4 種算法計算結果和計算效率
Table 4 Calculating results and its efficiency of four algorithms | 項 目 | CRS 算法 | SCE UA 算法 | SA 算法 | AS 算法 | | 模型參數 | k1/h-1 | 0.1890～0.2272 | 0.2114～0.2166 | 0.1920～0.2199 | 0.2133～0.2136 | | k2/h-1 | 0.0396～0.0453 | 0.0416～0.0433 | 0.0354～0.0454 | 0.0433～0.0434 | | k3/h-1 | 0.0106～0.0227 | 0.0101～0.0116 | 0.0100～0.0215 | 0.0100～0.0100 | | k4/h-1 | 0.0100～0.0789 | 0.0292～0.0639 | 0.0105～0.0679 | 0.0051～0.0988 | | k5/h-1 | 0.0010～0.0019 | 0.0010～0.0013 | 0.0009～0.0020 | 0.0008～0.0010 | | h1/mm | 18.303～19.121 | 17.700～19.970 | 18.109～19.993 | 18.870～19.945 | | h2/mm | 1.6035～7.9821 | 1.5460～4.3554 | 1.5281～9.3075 | 2.8560～4.1662 | | h3/mm | 930.99～1404.2 | 927.73～1307.4 | 701.88～1383.2 | 856.38～1447.9 | | h4/mm | 603.94～948.72 | 545.10～855.93 | 505.29～977.14 | 634.18～993.52 | | H1/mm | 2.2234～7.2466 | 6.0477～8.5449 | 1.3114～8.0765 | 7.3351～8.4248 | | H2/mm | 96.759～192.79 | 120.84～171.94 | 156.24～199.49 | 159.50～199.89 | | 目標函數值 | 函數值范圍 | 1.8648～1.9841 | 1.9377～1.9453 | 1.7774～2.1551 | 1.9330～1.9343 | | 均值 | 1.9537 | 1.9405 | 1.9469 | 1.9335 | | 標準方差 | 0.0504 | 0.0029 | 0.1347 | 0.0005 | | 平均CPU 時間(5000 次) | 2.7085e + 004 | 2.8219e + 004 | 2.5927e + 004 | 1.9965e + 004 |
對應于一定退火溫度下的接受準則,AS 算法中的下山單純形法在目標函數下降的方向上搜索全局最優,不僅有效避免了算法陷入局部極小值,而且大大提高了計算效率. 表4 的計算結果表明AS 算法的計算效率遠高于其他3 種直接算法. 由于SCE UA 的算法設計是在CRS算法基礎上引入生態學競爭演化模式,盡管全局優化的可靠性和穩定性有較大提高,但卻是以算法復雜性的增加和計算效率的下降為代價.
此外,研究中發現SA 算法引入的物理結晶過程使得算法本身控制參數的選擇具有較大困難. 對于不同優化問題,初始溫度、退火率和馬爾可夫鏈長度所組成的冷卻進度表有很大差異. 使用者必須通過數值實驗精確調整冷卻進度表使之能夠適應優化問題,這就降低了SA 算法的適應性并增加了實際應用過程中的困難.
實例研究表明,環境模型的參數優化過程遠較測試函數復雜,4 種直接算法的最佳估計參數和目標函數無一收斂到單一數值,并且采用統計特征量描述優化結果的不確定性也具有較大差異. 直接算法優化結果的不確定性分析不僅為進一步比較和理解算法特性提供了有效途徑,更為重要的是其揭示了優化算法不能完 全解決復雜模型參數的識別問題,因此基于環境模型結構和參數的不確定性研究成為必然. 4 結論 本文通過2 個模型實例的參數優化計算,對4 種直接算法的基本優化特性進行了分析比較,研究表明AS 算法和SCE UA 算法具有較高的可靠性和穩定性,AS 算法同時還具有較高的計算效率.
理論上,對于任何在整個參數可行域上進行隨機搜索的優化算法,當計算次數足夠大并適當選取了控制參數時,算法都能夠收斂到唯一的全局最優值. 然而實踐遠較之復雜:不僅需要在計算精度和計算效率之間做出艱難的選擇,更必須面對模型復雜性增加所導致的參數可識別問題. 由于算法隨機性的影響,特別是模型參數具有的不確定性,4 種直接算法得到的最佳參數估計值與目標函數無法較好吻合. 由此可見,基于參數優化算法框架下的模型參數仍然不可識別,即優化算法不能解決模型參數的不確定性問題. 參考文獻:
1 Beck M B. Water quality modeling : a review of the analysis of uncertainty. Wat . Res. Res. , 1987 , 23 ( 8) : 1393 ～1442.
2 Chen J , Wheater H S. Identification and uncertainty analysis of soil water retention models using lysimeter data. Wat .Res. Res. , 1999 , 35 (8) : 2401～2414.
3 Cooper V A et al. Evaluation of global optimization methods for conceptual rainfall2runoff model calibra tion. Wat . Sci.Tech. , 1997 , 36 (5) 53～60.
4 Duan Q Y et al. Shuffled complex evolution approach for ef2 fective and efficient global minimization. J . Optimization Theory and Applications , 1993 , 76 (3) : 501～521.
5 Duan Q Y et al. Optimal use of the SCE UA global opti2 mization method for calibrating watershed models. J . Hy－drol. , 1994 , 158 : 265～284.
6 Freedman V L et al. Parameter identifiability for catchment－scale erosion modelling : a comparison of optimization algo－rithms. J . Hydrol. , 1998 , 207 : 83～97.
7 康立山等. 非數值并行算法(第1 冊) ———模擬退火算法.北京:科學出版社,1994. 32～128.
8 Pan L H , Wu L S. A hybrid global optimization method for inverse estimation of hydraulic parameters : Annealing－Simplex method. WR972418 , Department of Soil & Environ－mental Science , University of California , Riverside. 1997. 87～196.
9 Press W H et al. Numerical Recipes in Fortran. Second Edi2 tion. New York : Cambridge Univ. Press , 1992. 6～97.
10 Price W L. A controlled random search procedure for global optimization. The Comput . J . , 1977 , 20 : 367～370.
11 Price WL. Global optimization by controlled random search.J . of Optimization Theory and Applications , 1983 , 40 : 333～348.
12 姚姚. 蒙特卡洛非線性反演方法及應用. 北京:冶金工業出版社, 1997. 43～212.
基金項目:高等學校優秀青年教師教學科研獎勵計劃資助項目 作者簡介:劉毅(1975～) ,男,博士研究生,主要從事環境系統分析方向的研究工作.
收稿日期:2001201210 ;修訂日期:2001204227 |
