基于粒子群和人工蜂群混合算法的氣動優化設計

楊梅花，夏露，張欣，米百剛

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

基于粒子群和人工蜂群混合算法的氣動優化設計

楊梅花，夏露，張欣，米百剛

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

現代啟發式智能算法存在全局與局部搜索能力的平衡問題，針對此問題，采用雙種群進化策略和信息交流機制，提出一種基于粒子群算法和人工蜂群算法相結合的新型混合優化算法——MABCPSO，并分別進行函數測試和翼型的氣動優化設計驗證。結果表明：MABCPSO新型混合優化算法具有更好的尋優能力，相比粒子群算法和人工蜂群算法，該算法能以更少的進化代數分別提高1.7%和2.2%的減阻效果。

粒子群算法；人工蜂群算法；混合算法；氣動優化設計

0 引 言

隨著計算機技術和計算流體力學的快速發展，采用數值優化方法進行飛行器的氣動設計，可以大幅提高優化設計的效率并節約設計成本。然而，傳統的基于目標函數梯度信息的算法都屬于局部優化方法，難以滿足實際工程對全局最優解的需要。現代啟發式智能算法雖然具有比傳統方法更好的全局性，但其求得的全局最優與局部最優之間卻多存在嚴重的沖突[1]。因此，探索發展具有良好的全局/局部搜索平衡能力的優化算法十分必要。

目前，被廣泛應用的、具有代表性的啟發式智能算法主要包括粒子群算法(Particle Swarm Optimization，簡稱PSO)和人工蜂群算法(Artificial Bee Colony，簡稱ABC)。粒子群算法是J.Kennedy和R.Eberhart[2]于1995年提出的進化算法，該算法作為一種基于群體智能的新型演化計算技術，具有簡單易實現、設置參數少、優化能力強等優點；但在進化后期，粒子群算法側重于開采能力，其局部搜索能力強，卻存在種群多樣性消失及早熟停滯等現象，從而陷入局部最優的缺陷。為了克服粒子群算法的上述缺點，國內外又先后提出了大量改進的粒子群優化算法，例如，具有繁殖子群能力的粒子群優化算法[3]、基于細菌覓食性的改進粒子群優化算法[4]、基于社會模型改進的粒子群優化算法[5]等。人工蜂群算法是D.Karaboga等[6-8]于2005年提出的一種模擬蜜蜂蜂群智能搜索行為的新型全局尋優算法，具有操作簡單、設置參數少、魯棒性高、收斂速度快、收斂精度高等優點；但其犧牲了算法的開采能力，側重于提高開發探索能力，局部搜索能力弱但全局搜索能力較強[9]。針對人工蜂群算法的缺點，國內外也提出了一些改進算法，例如，基于局部搜索的人工蜂群算法[10]、基于粒子群和人工蜂群的混合算法[11]、基于混沌局部搜索算子的人工蜂群算法[12]等。總體來看，對粒子群算法和人工蜂群算法的改進研究在一定程度上改善了原算法的固有缺陷，在全局/局部搜索能力之間做了一定的平衡。但是仍然存在優化后期全局/局部搜索能力相互沖突，容易陷入局部最優等問題。

針對現代啟發式智能算法在全局/局部搜索能力之間存在的如何相互平衡問題，本文提出一種基于PSO算法和ABC算法的新型混合優化算法——MABCPSO，并通過函數測試和RAE2822跨音速翼型減阻優化設計對MABCPSO算法進行分析驗證。

1 算法原理

1.1 粒子群算法

粒子群優化算法通過模擬鳥群的飛行覓食行為，將解空間中的設計變量看作一群沒有體積和質量的“鳥”(也稱為“粒子”)，將所求問題的最優解看作鳥群尋找的“食物”，所有粒子通過不斷改變飛行方向和距離，并以個體極值pbest和群體極值gbest來引導自己調整飛行狀態，更新速度和位置，從而趨向問題的最優解。PSO算法的流程如圖1所示。

粒子可根據式(1)～式(2)來更新：

(1)

(2)

式中：x為粒子當前的位置；v為粒子的速度；w為慣性權重因子，它是影響算法性能的重要參數，其大小決定了粒子對當前速度繼承的多少，w的取值范圍通常為[0.4，0.9]，其在進化過程中采用線性遞減策略；c1、c2為學習因子，用來調節粒子飛向pbest和gbest的步長；r1、r2為產生于(0，1)的隨機數；下標i和j為第i個粒子和第j維；上標t為進化代數。

1.2 人工蜂群算法

人工蜂群算法包含三個組成部分：采蜜蜂、觀察蜂和偵察蜂，采蜜蜂和觀察蜂執行開采過程，偵察蜂執行探索過程。每個蜜源的位置代表一個可能解，每個蜜源的花蜜量代表相應解的質量或適應度。首先，隨機產生初始群體P，即SN個初始解，SN為采蜜蜂的數量或蜜源數量(采蜜蜂的數量=蜜源數量)。每個解xi(i=1,2,3,...,SN)是一個D維向量(D為搜索空間的維數)。初始化之后，以初始解為基礎，采蜜蜂、觀察蜂和偵察蜂開始進行循環搜索。采蜜蜂記住以前的最優解，在采蜜源的鄰域內搜索產生新解，搜索公式為

(3)

式中：k為不同于i的蜜源；j為隨機選擇的下標；φij為隨機數，其范圍是[-1,1]，控制著xij鄰域內蜜源位置的產生。

采蜜蜂采用貪婪準則，當搜索解優于記憶最優解時，替換記憶解；反之，則保持記憶解不變。采蜜蜂完成搜索過程后，將記憶中的蜜源信息與觀察蜂共享。觀察蜂按照與花蜜量相關的概率選擇一個蜜源位置，并像采蜜蜂一樣對記憶中的位置做一定的改變，并檢查新候選位置的花蜜量，若新位置優于記憶中的位置，則用新位置替換原來的蜜源位置，否則保留原位置。一個觀察蜂選擇某個蜜源的概率為

(4)

式中：fiti為解xi的適應度值。

假如蜜源位置xi經過“limit”次采蜜蜂和觀察蜂的循環搜索后，無法被改進，那么該位置將被放棄，此時采蜜蜂轉變為偵察蜂，并隨機搜索一個蜜源替換原蜜源。偵察蜂搜索新蜜源的操作公式為

(5)

1.3MABCPSO混合優化算法

通過分析上述兩種算法可以看出：PSO算法在進化后期，側重于開采能力，局部搜索能力強，易陷入局部最優；而ABC算法犧牲了開采能力，側重于提高開發探索能力，局部搜索能力弱但全局搜索能力較強。

本文結合PSO算法和ABC算法的優勢，提出一種新的混合優化算法——MABCPSO。MABCPSO將種群隨機地等分為兩組，一組種群中的個體按照PSO算法進化尋優，另一組種群中的個體按照ABC算法進行尋優。為了保證該算法的探索開發能力，同時提高開采能力、加快尋優速度，在ABC算法中提出新的搜索方程：

(6)

式中：φij為隨機數，其范圍是[0，1.5]；pij為個體極值；pg為群體極值。

新的搜索方程由于有群體極值pg的引導，在保證探索開發能力的同時也提高了其開采能力。MABCPSO混合優化算法在全局搜索能力和局部搜索能力之間建立一個相對較好的平衡。

MABCPSO算法在進化過程，利用一種信息交流機制來交流搜尋信息，使信息能夠在兩個種群中傳遞，有助于個體避免錯誤的信息判斷而陷入局部最優。

兩個種群的信息交流機制如下：

(1) 式(1)中的gbest用兩個種群中的最優個體來替換；(2) 式(6)中的xkj從采用PSO算法進化的種群中隨機選取；(3) 式(6)中的pij、pg用兩個種群中的最優個體來替換。

MABCPSO算法的具體步驟如下：

(1) 設置群體規模N，最大迭代次數maxgen，限定循環次數limit，學習因子c1、c2，慣性權重的初始值w1和終值w2；

(2) 將群體等分為兩個種群，種群1中的個體按照PSO算法尋優進化，種群2中的個體按照改進的ABC算法尋優進化；

(3) 設置初始迭代步數iter=0；

(4) 種群1和種群2分別按照式(1)、式(2)和式(6)產生新解，并計算其適應度的值，利用信息交流機制交換信息；

(5) 更新種群1中的pbest和gbest，如果種群2中新解的適應度值優于xi，則用其替換xi，否則不變。

(6) 計算xi的適應度值，并根據式(4)計算概率pi；觀察蜂根據pi選擇蜜源，并按照式(6)產生新解，如果新解的適應度值優于xi，則用其替換xi，否則不變；

(7) 經過limit次循環后，判斷是否有需要丟棄的解，若有，則根據式(5)產生滿足約束條件的新解；

(8) 更新迭代步數iter=iter+1，并記錄當前整個群體中的最佳個體，如果滿足精度要求或進化已達到最大迭代次數，則終止算法，否則轉至步驟(4)。

2 算法函數測試

為了測試并分析算法的尋優能力，分別對PSO算法、ABC算法和MABCPSO算法進行函數測試。選取的測試函數，為函數測試中常用的基準測試函數，基準測試函數共有12個[13-15]，本文選取三個典型的測試函數。

(1)Sphere函數：

全局最優值：minf1(x)=f1(0,0,…,0)

該函數為非線性對稱單峰函數，大多數算法都能輕松達到優化效果，主要用于測試算法的尋優精度。

(2)Rastrigin函數

(-5.12≤xi≤5.12)

全局最優值：min f2(x)=f2(0,0,…,0)

該函數具有大量正弦拐點并且排列得很深的局部最優點，優化算法很容易在通往全局最優點的路徑上陷入一個局部最優點。

(3)Ackley函數

(-32≤xi≤32)

全局最優值：min f3(x)=f3(0,0,…,0)

該函數為多峰函數，存在大量局部最優點，優化算法很容易陷入局部最優。

測試中，三個函數的變量取15維，三種算法的種群規模均為60，慣性權重因子w采用線性遞減策略隨進化代數由0.9線性遞減到0.4[16]，c1=c2=2，limit=100，分別用三種算法進行100次優化，每次優化算法各自進化400代，統計這100次優化的平均目標函數值和搜索到的最小目標函數值。函數測試結果如表1所示。

表1 函數測試結果

從表1可以看出：在調用目標函數次數相同的情況下，采用MABCPSO算法所求得的平均值和最小值，其精度都更高，收斂效果更好，表明該新型混合算法具有更好的尋優能力。

為了探索MABCPSO新型混合算法在高維問題中的尋優能力，將三個測試函數的變量取40維，三種算法的種群規模均為90，慣性權重因子w采用線性遞減策略隨進化代數由0.9線性遞減到0.4，c1=c2=2，limit=100，分別用三種算法進行100次優化，每次優化算法各自進化400代，統計這100次優化的平均目標函數值和搜索到的最小目標函數值。函數測試結果如表2所示。

表2 高維函數測試結果

從表2可以看出：在高維尋優問題中，三種算法的尋優精度都不同程度地降低，但是在調用目標函數次數相同的情況下，采用MABCPSO算法所求得的平均值和最小值，其精度仍然是三種算法中最高的，收斂效果也最好，表明該新型混合算法對于高維問題仍然有較好的尋優能力。

3 翼型優化設計

將MABCPSO新型混合算法應用于實際二維翼型氣動優化設計中，進一步分析和驗證其優化效果。

優化設計模型包括三個基本要素：目標函數、約束和設計變量。對于翼型的氣動優化設計而言，一般確定了設計點后，可選取表示升阻特性或力矩特性的氣動系數作為目標函數，對翼型的厚度和面積等進行限制作為約束，將翼型的幾何形狀作為設計變量。

3.1 翼型的幾何表示

翼型的幾何表示方法會直接影響到翼型優化設計的結果，采用CST參數化方法[17]描述翼型的幾何形狀。用12個設計變量進行翼型的幾何表示，上下翼面各有6個設計變量進行控制。用CST方程直接表示翼型的幾何坐標，CST方程的表達式為

上表面

yu=C(x)·Su(x)+x·yTEu

下表面

yl=C(x)·Sl(x)+x·yTEl

式中：下標u、l分別為翼型的上下表面；yTEu、yTEl分別為上下表面后緣的y坐標。

類函數C(x)定義為

C(x)=xN1·(1-x)N2

型函數S(x)定義為

對于一般的翼型，類函數C(x)中的N1和N2分別取0.5和1.0；Si(x)為Bernstein多項式函數；Aui和Ali為待定系數，只要確定了Aui和Ali，即確定了整個翼型。

3.2 翼型優化模型及約束處理

以RAE2822跨音速翼型為研究對象，考慮其在11km高空的巡航能力，對其進行單點減阻優化設計。設計狀態為：Ma=0.73，迎角α=2.3°，Re=5.51×106[18-20]。

優化目標為：最小化翼型阻力系數CD；約束條件為：翼型的面積、最大厚度不減小，力矩系數不減小，升力系數CL基本保持不變。氣動特性的分析計算采用基于多重網格的可壓縮N-S方程數值解法。一般的，處理約束條件的方法有：懲罰函數法、解碼器法、修補法、算法修正法和可行解搜索法等，目前使用最廣泛的是懲罰函數法，即罰函數法[21]，故本文采用罰函數法來處理約束問題。

3.3 優化算法

分別采用PSO算法、ABC算法和MABCPSO算法進行翼型優化。三種算法的種群規模均為60且初始種群相同，進化代數為50代，即三種算法經過相同次數的流場計算，慣性權重因子w隨進化代數由0.9線性遞減到0.4，c1=c2=2，limit=100。

3.4 優化結果及分析

三種算法計算得到的阻力系數的優化時間歷程如圖2所示，可以看出：在相同的初始種群下，僅經過一代進化，即從第一代起，MABCPSO算法的尋優結果就明顯優于PSO算法和ABC算法；在之后的進化過程中，MABCPSO算法在達到給定的阻力系數時，所需進化的代數明顯少于PSO算法和ABC算法；MABCPSO算法收斂到的最優值也明顯優于PSO算法和ABC算法。表明MABCPSO算法收斂速度更快、效率更高、尋優能力更強。

優化后翼型的幾何形狀與初始翼型的對比如圖3所示，可以看出：優化后的翼型，其上表面頂點向后移動，PSO算法優化后的翼型最大厚度位置較初始翼型向后，MABCPSO算法和ABC算法優化后的翼型最大厚度位置與初始翼型基本相同。

優化前后翼型的表面壓力分布如圖4所示，可以看出：優化后，前加載和后加載增大，翼型上表面壓力過渡更為平滑，在很大程度上減弱了激波，其中MABCPSO算法得到的翼型上表面壓力分布過渡得最為平滑，顯著地消除了激波，表現出更好的氣動性能。

三種算法的優化效果對比如表3所示，可以看出：在相同的計算量下，采用PSO算法進行優化，阻力系數降低了13.9%；采用ABC算法進行優化，阻力系數降低了13.4%；而采用本文提出的MABCPSO算法，阻力系數降低了15.6%。

表3 優化效果對比

綜上所述，與初始翼型相比，三種算法均能達到減阻目的，且均滿足規定的約束要求；對于像氣動設計這樣高復雜性和強非線性的問題，MABCPSO算法將粒子群算法與改進的人工蜂群算法有機地結合起來，共享最優信息，保持種群的多樣性，使得該算法能夠有效地跳出局部最優，加快收斂速度，在全局/局部搜索能力之間展現出很好的平衡，具有更好的尋優性能。

4 結 論

(1) 本文提出的基于PSO算法和ABC算法相結合的新型混合優化算法——MABCPSO，在全局和局部搜索能力之間建立了一個良好的平衡，該算法在進化過程中，利用改進的ABC算法的開發能力和波動性，保持了種群的多樣性，避免種群陷入局部最優；同時，利用PSO算法的開采能力，重點搜索具有較高適應度值的個體，加快了收斂速度。

(2)MABCPSO新型混合優化算法能夠很好地應用于二維翼型的跨聲速氣動優化設計，相比傳統優化方法，該算法能夠更加快速地得到阻力更小的新翼型。

(3)MABCPSO新型混合優化算法盡管尋優性能良好，但是對于像氣動設計這樣的強非線性復雜問題，還需進一步明確其內在機理，尋求改進，建立更為高效的優化模型。

[1] 李丁, 夏露. 改進的粒子群優化算法在氣動設計中的應用[J]. 航空學報, 2012, 33(10): 1809-1816.LiDing,XiaLu.Applicationofimprovedparticleswarmoptimizationalgorithmtoaerodynamicdesign[J].ActaAeronauticaetAstronauticaSinica, 2012, 33(10): 1809-1816.(inChinese)

[2]KennedyJ,EberhartR.Particleswarmoptimization[C]∥IEEEInternationalConferenceonNeuralNetworks. 1995(4): 1942-1948.

[3]LvbjergM,GroupEP,RasmussenTK,etal.Hybridparticleswarmoptimiserwithbreedingandsubpopulations[EB/OL]. (2001)[2017-02-11].https:∥www.researchgate.net/publication/2490119_Hybrid_Particle_Swarm_Optimiser_with_Breeding_and_Subpopulations.

[4] 夏露, 李丁. 基于細菌覓食性的改進粒子群優化算法翼型設計[J]. 空氣動力學學報, 2012, 30(4): 533-538.XiaLu,LiDing.Animprovedparticleswarmoptimizationbasedonbacterialforagingtoairfoildesign[J].ActaAerodynamicaSinica, 2012, 30(4): 533-538.(inChinese)

[5] 夏露, 孫騰騰, 李丁. 基于社會模型改進粒子群算法的翼型優化設計[J]. 航空計算技術, 2013, 43(5): 1-6.XiaLu,SunTengteng,LiDing.Airfoilaerodynamicoptimizationdesignusingimprovedparticleswarmoptimizationbasedonsocialmodel[J].AeronauticalComputingTechnique, 2013, 43(5): 1-6.(inChinese)

[6]KarabogaD.Anideabasedonhoneybeeswarmfornumericaloptimization[R].TechnicalReport-TR06, 2005.

[7]KarabogaD,BasturkB.Apowerfulandefficientalgorithmfornumericalfunctionoptimization:artificialbeecolony(ABC)algorithm[J].JournalofGlobalOptimization, 2007, 39(3): 459-471.

[8]KarabogaD,BasturkB.Artificialbeecolony(ABC)optimizationalgorithmforsolvingconstrainedoptimizationproblems[J].FoundationsofFuzzyLogic&SoftComputing, 2007, 11(3): 789-798.

[9] 余勝威.MATLAB優化算法案例分析與應用[J]. 北京: 清華大學出版社, 2014.YuShengwei.Matlaboptimizationalgorithmcaseanalysisandapplication[M].Beijing:TsinghuaUniversityPress, 2014.(inChinese)

[10] 劉三陽, 張平, 朱明敏. 基于局部搜索的人工蜂群算法[J]. 控制與決策, 2014(1): 123-128.LiuSanyang,ZhangPing,ZhuMingmin.Artificialbeecolonyalgorithmbasedonlocalsearch[J].Control&Decision, 2014(1): 123-128.(inChinese)

[11] 王志剛. 基于粒子群和人工蜂群算法的混合優化算法[J]. 科學技術與工程, 2012, 12(20): 4921-4925.WangZhigang.Hybridoptimizationalgorithmbasedonparticleswarmoptimizationandartificialbeecolonyalgorithm[J].ScienceTechnology&Engineering, 2012, 12(20): 4921-4925.(inChinese)

[12] 王翔, 李志勇, 許國藝, 等. 基于混沌局部搜索算子的人工蜂群算法[J]. 計算機應用, 2012, 32(4): 1033-1036.WangXiang,LiZhiyong,XuGuoyi,etal.Artificialbeecolonyalgorithmbasedonchaoslocalsearchoperator[J].JournalofComputerApplications, 2012, 32(4): 1033-1036.(inChinese)

[13] 欒麗君, 譚立靜, 牛奔. 一種基于粒子群優化算法和差分進化算法的新型混合全局優化算法[J]. 信息與控制, 2007, 36(6): 708-714.LuanLijun,TanLijing,NiuBen.Anovelhybridglobaloptimizationalgorithmbasedonparticleswarmoptimizationanddifferentialevolution[J].Information&Control, 2007, 36(6): 708-714.(inChinese)

[14] 高鷹, 謝勝利. 基于模擬退火的粒子群優化算法[J]. 計算機工程與應用, 2004, 40(1): 47-50.GaoYing,XieShengli.Particleswarmoptimizationalgorithmsbasedonsimulatedannealing[J].ComputerEngineering&Applications, 2004, 40(1): 47-50.(inChinese)

[15] 周丹, 葛洪偉, 張歡慶, 等. 基于健康度的人工蜂群粒子群算法[J]. 計算機工程與應用, 2016, 52(18): 62-67.ZhouDan,GeHongwei,ZhangHuanqing,etal.Particlehealthdegreebasedartificialbeecolonyparticleswarmoptimization[J].ComputerEngineering&Applications, 2016, 52(18): 62-67.(inChinese)

[16] 胡建秀, 曾建潮. 微粒群算法中慣性權重的調整策略[J]. 計算機工程, 2007, 33(11): 193-195.HuJianxiu,ZengJianchao.Selectiononinertiaweightofparticleswarmoptimization[J].ComputerEngineering, 2007, 33(11): 193-195.(inChinese)

[17] 卜月鵬, 宋文萍, 韓忠華, 等. 基于CST參數化方法的翼型氣動優化設計[J]. 西北工業大學學報, 2013, 31(5): 829-836.BuYuepeng,SongWenping,HanZhonghua,etal.AerodynamicoptimizationdesignofairfoilbasedonCSTparameterizationmethod[J].JournalofNorthwesternPolytechnicalUniversity, 2013, 31(5): 829-836.(inChinese)

[18] 王偉, 白俊強, 張揚, 等. 基于人工魚群算法及Hicks-Henne型函數的翼型優化設計研究[J]. 科學技術與工程, 2011, 11(24): 5870-5874.WangWei,BaiJunqiang,ZhangYang,etal.AerodynamicoptimizationdesignofairfoilbasedonAFSAandHicks-Henneshapefunction[J].ScienceTechnology&Engineering, 2011, 11(24): 5870-5874.(inChinese)

[19] 李丁, 夏露. 一種混合粒子群優化算法在翼型設計中的應用[J]. 航空計算技術, 2010, 40(6): 66-71.LiDing,XiaLu.Applicationofahybridparticleswarmoptimizationtoairfoildesign[J].AeronauticalComputingTechnique, 2010, 40(6): 66-71.(inChinese)

[20] 劉南, 白俊強, 邱亞松, 等. 基于本征正交分解的氣動外形設計空間重構方法研究[J]. 西北工業大學學報, 2015, 33(2): 171-177.LiuNan,BaiJunqiang,QiuYasong,etal.Investigatingaerodynamicshapedesignspacereconstructionusingproperorthogonaldecomposition(POD)[J].JournalofNorthwesternPolytechnicalUniversity, 2015, 33(2): 171-177.(inChinese)

[21] 張利彪. 基于粒子群優化算法的研究[D]. 長春: 吉林大學, 2004.ZhangLibiao.Researchbasedonparticleswarmoptimization[D].Changchun:JilinUniversity, 2004.(inChinese)

(編輯：馬文靜)

Aerodynamic Optimization Design Based on Hybrid Optimization Algorithm of Particle Swarm Optimization and Artificial Bee Colony Algorithm

Yang Meihua, Xia Lu, Zhang Xin, Mi Baigang

(School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China)

The modern intelligent algorithm is of the balance problem of global and local search ability. Because of this, a new hybrid optimization algorithm based on particle swarm optimization and artificial bee colony algorithm(it is called hybrid method of ABC and PSO-MABCPSO) is proposed, which is based on the two-group evolutionary strategy and information sharing mechanism, and the function tests and airfoil aerodynamic optimization design validation are carried out respectively. The results show that MABCPSO new hybrid optimization algorithm has better searching ability, which can improve the drag reduction effect by 1.7% and 2.2% with less evolutionary algebra than particle swarm optimization(PSO) and artificial bee colony algorithm(ABC) algorithm respectively.

particle swarm optimization(PSO); artificial bee colony algorithm(ABC); hybrid algorithm; aerodynamic optimization design

2017-02-11；

2017-04-12

楊梅花,541620411@qq.com

1674-8190(2017)02-182-08

V211.4

A

10.16615/j.cnki.1674-8190.2017.02.010

楊梅花(1992-)，女，碩士研究生。主要研究方向：飛行器設計。

夏 露(1977-)，女，副教授，碩導。主要研究方向：飛行器總體設計、飛行器氣動外形隱身一體化設計。

張 欣(1993-)，男，碩士研究生。主要研究方向：飛行器設計。

米百剛(1989-)，男，博士研究生。主要研究方向：飛行器氣動設計、數值模擬。