馬洋 榮偉 江長(zhǎng)虹 張青斌 霍霖

（1國(guó)防科學(xué)技術(shù)大學(xué)航天科學(xué)與工程學(xué)院，長(zhǎng)沙 410073）

（2北京空間機(jī)電研究所，北京 100094）

1 引言

采用球冠倒錐外形的返回艙在再入過(guò)程中將經(jīng)歷十分嚴(yán)酷而復(fù)雜的氣動(dòng)力熱環(huán)境。球冠表面受氣流強(qiáng)烈壓縮，承受著高溫、高壓氣流，氣流在球冠和倒錐聯(lián)結(jié)的肩部附近發(fā)生劇烈膨脹，在后體錐面上形成高速低壓的流動(dòng)。精心設(shè)計(jì)返回艙外形可以有效地減少嚴(yán)酷力熱環(huán)境對(duì)其再入過(guò)程產(chǎn)生的不利影響，也是保證成功完成再入任務(wù)的關(guān)鍵點(diǎn)之一。返回艙在大氣層內(nèi)利用空氣阻力減阻，要獲得盡可能大的阻力以盡快減小飛行速度可以通過(guò)外形設(shè)計(jì)達(dá)到目的，但阻力的增大會(huì)受到諸多因素的制約。首先由于阻力與升阻比存在著天然的“矛盾關(guān)系”，而升阻比又是控制再入彈道、改進(jìn)落點(diǎn)精度最關(guān)鍵的因素，因而在增大阻力時(shí)必須考慮到升阻比的要求。另外在高速再入過(guò)程中產(chǎn)生的氣動(dòng)加熱也會(huì)對(duì)返回艙外形設(shè)計(jì)提出特殊要求。

氣動(dòng)特性的精確預(yù)示是進(jìn)行返回艙外形優(yōu)化設(shè)計(jì)的基礎(chǔ)，陳河梧采用試驗(yàn)手段研究了 Ma=5~8，攻角在–27°~2°條件下，類神舟飛船返回艙的高超聲速氣動(dòng)特性，認(rèn)為返回艙球冠主導(dǎo)小攻角的氣動(dòng)力變化，大攻角時(shí)倒錐的氣動(dòng)力貢獻(xiàn)加大，升力和俯仰力矩出現(xiàn)比較明顯的非線性增量[1]。紀(jì)楚群和呂俊明采用計(jì)算流體力學(xué)(Computational Fluid Dynam ics, CFD)數(shù)值仿真手段對(duì)鐘形返回艙和火星再入返回艙進(jìn)行數(shù)值模擬，對(duì)比分析了基于Euler方程和Navier-Stokes(NS)方程的CFD模擬結(jié)果，并分析了化學(xué)非平衡效應(yīng)的影響，研究表明，在高超聲速情況下，由于物體迎風(fēng)面的壓力遠(yuǎn)遠(yuǎn)大于背風(fēng)面壓力，背風(fēng)面壓力計(jì)算誤差對(duì)整體氣動(dòng)力特性影響很小，同時(shí)氣流在物面的分離情況不嚴(yán)重，因此用 Euler方程可足夠準(zhǔn)確地預(yù)測(cè)返回艙高超聲速情況的流場(chǎng)及氣動(dòng)特性[2-3]。化學(xué)非平衡效應(yīng)對(duì)返回艙氣動(dòng)特性影響較小，但它會(huì)對(duì)激波形狀和駐點(diǎn)位置產(chǎn)生影響。返回艙外形優(yōu)化設(shè)計(jì)研究方面，李治宇采用Euler方法在ISIGHT集成軟件環(huán)境下對(duì)返回艙外形展開(kāi)優(yōu)化設(shè)計(jì)，得到了升阻比和駐點(diǎn)熱流的優(yōu)化前緣[4]。John E Theisinger對(duì)火星再入返回艙外形進(jìn)行了詳細(xì)的氣動(dòng)力熱分析，通過(guò)對(duì)外形進(jìn)行優(yōu)化分析，對(duì)阻力、配平狀態(tài)下的穩(wěn)定性、表面熱流、駐點(diǎn)熱流、容積率、容積、升阻比等各種相互制約的性能進(jìn)行協(xié)調(diào)[5-6]。

本文氣動(dòng)特性分析采用三維 NS方程，不考慮真實(shí)氣體的化學(xué)非平衡效應(yīng)，引入代理模型以減少計(jì)算開(kāi)銷，采用Fay-Riddell公式計(jì)算駐點(diǎn)熱流，并利用成熟的改進(jìn)非劣分類遺傳算法(Nondom inated Sorting Genetic A lgorithm-II, NSGA-II)[7]，對(duì)參數(shù)化設(shè)計(jì)的返回艙外形進(jìn)行多目標(biāo)優(yōu)化設(shè)計(jì)，該優(yōu)化設(shè)計(jì)方法在返回艙方案設(shè)計(jì)階段具有一定的應(yīng)用價(jià)值。

2 研究方法

在氣動(dòng)外形優(yōu)化設(shè)計(jì)中，設(shè)計(jì)結(jié)果的精度和計(jì)算效率往往相互沖突。結(jié)果的精度主要由氣動(dòng)特性計(jì)算精度和尋優(yōu)算法精度決定，本文選擇經(jīng)典遺傳算法尋優(yōu)，尋優(yōu)算法精度和效率不在本文的研究之列。氣動(dòng)特性計(jì)算方法主要有工程估算方法、CFD數(shù)值仿真以及兩種方法相結(jié)合的混合方法等。工程估算效率最高，其精度也最差，較少直接用于優(yōu)化設(shè)計(jì)。CFD方法計(jì)算精度很高，許多情況下可以替代試驗(yàn)研究，但在飛行器氣動(dòng)外形優(yōu)化設(shè)計(jì)中，需要反復(fù)修改氣動(dòng)外形，多次進(jìn)行氣動(dòng)特性計(jì)算，采用CFD方法的計(jì)算成本即使是對(duì)最先進(jìn)的計(jì)算資源來(lái)說(shuō)仍顯得力不從心?；旌戏椒ㄓ捎诰哂泄こ坦浪愀咝Ш虲FD仿真高精度的潛質(zhì)，在優(yōu)化設(shè)計(jì)中正受到一些研究人員的重視[8-9]，但其計(jì)算精度受CFD結(jié)果和所采取的擬合算法的限制，其應(yīng)用還需要進(jìn)一步的研究。

本文引入代理模型技術(shù)以減小CFD計(jì)算的開(kāi)銷，并采用改進(jìn)的EI（Expected Improvement）函數(shù)加點(diǎn)策略[10]，充分發(fā)揮代理模型的高效性，并以高效的加點(diǎn)策略保證結(jié)果的精度，很好地解決了氣動(dòng)外形優(yōu)化設(shè)計(jì)中精度與效率的問(wèn)題，具體優(yōu)化設(shè)計(jì)方法如圖1所示。首先采用試驗(yàn)設(shè)計(jì)方法產(chǎn)生構(gòu)建代理模型的樣本點(diǎn)，然后采用CFD方法計(jì)算樣本點(diǎn)的氣動(dòng)特性，并基于計(jì)算結(jié)果構(gòu)建代理模型，再判斷代理模型計(jì)算結(jié)果與CFD計(jì)算結(jié)果的誤差是否滿足精度要求，如果不滿足則按加點(diǎn)策略加點(diǎn)重新進(jìn)行CFD計(jì)算并構(gòu)建代理模型，直至滿足收斂要求，最后用構(gòu)建好的代理模型代替CFD仿真進(jìn)行外形優(yōu)化設(shè)計(jì)。

圖1 優(yōu)化設(shè)計(jì)流程Fig.1 Optimization design process

2.1 氣動(dòng)特性分析方法及驗(yàn)證

氣動(dòng)性能計(jì)算采用成熟的Fluent軟件。選取可實(shí)現(xiàn)的k-ε兩方程湍流模型，并配合使用非平衡壁面函數(shù)，這樣更適合于模擬返回艙附近擾流流場(chǎng)中出現(xiàn)的大分離和大漩渦特性，對(duì)流項(xiàng)離散采用二階AUSM格式，返回艙表面滿足無(wú)滑移邊界條件，進(jìn)口取來(lái)流參數(shù)，出口數(shù)值邊界條件采用外推方式獲得。

采用結(jié)構(gòu)網(wǎng)格離散流場(chǎng)，由于流場(chǎng)的對(duì)稱性，只對(duì)一半流場(chǎng)進(jìn)行網(wǎng)格劃分，在靠近返回艙表面對(duì)網(wǎng)格進(jìn)行加密，返回艙和用于數(shù)值精度驗(yàn)證的計(jì)算網(wǎng)格如圖2和圖3所示。

圖2 返回艙計(jì)算網(wǎng)格Fig.2 Computing grid of reentry capsul

圖3 數(shù)值驗(yàn)證算例計(jì)算網(wǎng)格Fig.3 Computing grid of CFD validation

為驗(yàn)證本文CFD計(jì)算模型的準(zhǔn)確性，對(duì)文獻(xiàn)[11]提供的火星探測(cè)實(shí)驗(yàn)飛行器外形進(jìn)行數(shù)值模擬，并與文獻(xiàn)[3]和[11]計(jì)算結(jié)果進(jìn)行對(duì)比（如表1所示，其中α、CL、CD、L/D和CM分別表示攻角，升、阻力系數(shù)，升阻比和俯仰力矩系數(shù)），可以看出，本文 CFD計(jì)算得到的氣動(dòng)特性與文獻(xiàn)結(jié)果的誤差約為2~9%，數(shù)值仿真模型具有較高的可信度。

表1 火星探測(cè)實(shí)驗(yàn)飛行器氣動(dòng)特性計(jì)算驗(yàn)證Tab.1 Aerodynam ic computing accuracy verification of Mars Science Laboratory

2.2 返回艙駐點(diǎn)熱流計(jì)算

由于熱流數(shù)值仿真需要比氣動(dòng)力數(shù)值模擬更多的計(jì)算資源，而工程經(jīng)驗(yàn)公式能在趨勢(shì)上給出氣動(dòng)熱與其它性能的大致協(xié)調(diào)關(guān)系，因此在方案設(shè)計(jì)階段采用經(jīng)驗(yàn)公式計(jì)算駐點(diǎn)熱流不失為一種折衷的選擇。返回艙駐點(diǎn)熱流采用簡(jiǎn)化的Fay-Riddell公式[12]：

2.3 代理模型構(gòu)建

代理模型技術(shù)主要包括多項(xiàng)式響應(yīng)面(Polynom ial Response Surface)、徑向基函數(shù)(Radial Basis Function)、神經(jīng)網(wǎng)絡(luò)(artificial neural networks)和Kriging方法等，前兩種方法在處理非線性強(qiáng)和高維數(shù)問(wèn)題時(shí)表現(xiàn)不好，神經(jīng)網(wǎng)絡(luò)方法雖然較適合處理非線性強(qiáng)和高維數(shù)問(wèn)題，但其預(yù)測(cè)精度和魯棒性在一定程度上不及Kriging方法[13]。

本文采用Kriging方法構(gòu)建替代CFD分析的代理模型[14]，為了增強(qiáng)代理模型的尋優(yōu)精度和效率，借鑒經(jīng)典的EI函數(shù)加點(diǎn)方法[10]并加以改進(jìn)，提出在現(xiàn)有樣本點(diǎn)基礎(chǔ)上，同時(shí)加入EI值最大的點(diǎn)和Kriging模型預(yù)測(cè)方差最大的點(diǎn)來(lái)改善模型預(yù)測(cè)精度，EI值和 Kriging模型預(yù)測(cè)方差的具體計(jì)算方法參考文獻(xiàn)[10,14]，這樣的加點(diǎn)尋優(yōu)方法一方面具有EI函數(shù)方法兼顧局部和全局尋優(yōu)的特點(diǎn)，另一方面在模型預(yù)測(cè)方差最大處加點(diǎn)更注重模型的全局精度，并且多點(diǎn)加點(diǎn)方式也適合于并行計(jì)算，有利于提高優(yōu)化設(shè)計(jì)的效率。

2.4 尋優(yōu)算法

遺傳算法模擬生物在自然環(huán)境中的遺傳和進(jìn)化過(guò)程，是一種自適應(yīng)、全局優(yōu)化、概率搜索算法，它具有良好的魯棒性、全局性和并行性，在當(dāng)今工程技術(shù)領(lǐng)域的多目標(biāo)優(yōu)化問(wèn)題中應(yīng)用非常廣泛。在眾多的多目標(biāo)遺傳算法中，NSGA采用簡(jiǎn)潔明晰的非優(yōu)超排序機(jī)制，使算法具有逼近Pareto最優(yōu)解的能力，采用排擠機(jī)制保證得到的Pareto最優(yōu)解具有良好的分布，Deb在NSGA基礎(chǔ)上改進(jìn)得到NSGA-II[15]，引入精英保留策略對(duì)產(chǎn)生的非劣解進(jìn)行非劣快速分類和密度估計(jì)操作，比較其擁擠度，確定是否可接受為Pareto最優(yōu)解。NSGA2算法提出了新的基于分級(jí)的快速非劣排序算法，大大降低了計(jì)算復(fù)雜度；同時(shí)為使Pareto最優(yōu)解散布范圍較大，盡可能均勻遍布，引入了擁擠距離的概念，并采用擁擠比較算子代替需要計(jì)算共享參數(shù)的適應(yīng)度共享方法；最后引入精英保留機(jī)制，不但擴(kuò)大了采樣空間，而且有利于保持優(yōu)良個(gè)體，迅速提高了種群的整體水平。由于NSGA-II算法的上述優(yōu)點(diǎn)和可靠性，本文以其作為優(yōu)化算法。

3 多目標(biāo)優(yōu)化與分析

3.1 返回艙外形描述

返回艙外形如圖4所示，可以采用最大直徑dm、最大高度H、端框直徑dD和高度HD、大底半徑RN、前體半徑RF、大底倒圓半徑 Rc和后體傾角θ來(lái)表示。

3.2 問(wèn)題描述

在返回艙再入大氣的過(guò)程中，需要依靠大氣阻力使其減速才能完成再入返回的任務(wù)，同時(shí)返回艙也需要盡可能大的升阻比，這樣能實(shí)現(xiàn)橫向和縱向的機(jī)動(dòng)，得到更寬的再入走廊，從而改善落點(diǎn)的精度。但是阻力和升阻比從來(lái)就存在著“矛盾”，要想得到滿足再入返回要求的返回艙外形必須同時(shí)考慮這兩個(gè)目標(biāo)，進(jìn)行多目標(biāo)優(yōu)化。另外在返回艙高速再入過(guò)程中，嚴(yán)重的氣動(dòng)加熱是必須考慮的問(wèn)題，同時(shí)為滿足載人以及設(shè)備存放的要求，返回艙內(nèi)的空間也是設(shè)計(jì)需要關(guān)注的地方，本文將返回艙表面最大熱流（駐點(diǎn)熱流）和容積率作為約束條件加入到優(yōu)化設(shè)計(jì)過(guò)程中，以體現(xiàn)多種約束條件下的設(shè)計(jì)。

在圖4所示的8個(gè)外形參數(shù)中，最大高度H和最大直徑dm受返回艙總體規(guī)模的限制變化較小，端框直徑dD和高度HD由于處在背風(fēng)區(qū)對(duì)返回艙的氣動(dòng)力和氣動(dòng)熱特性影響較小，對(duì)容積特性的影響也較小，它們都不適宜作為本文優(yōu)化問(wèn)題的設(shè)計(jì)變量。后體傾角θ雖然可變范圍較大且對(duì)返回艙性能有較大影響，但在大底半徑RN、前體半徑RF和大底倒圓半徑Rc確定的前提下，θ是確定的，即θ、RN、RF和Rc只有是3個(gè)參數(shù)是獨(dú)立的。本文選取RN、RF和Rc作為設(shè)計(jì)變量，它們的取值范圍如表2。

圖4 返回艙外形參數(shù)Fig.4 Configuration parameters of reentry capsule

表2 設(shè)計(jì)變量及其變化范圍Tab.2 Design variables and their variation range

在某一典型的再入條件下[1]，返回艙外形優(yōu)化設(shè)計(jì)問(wèn)題可以這樣描述：

3.3 優(yōu)化結(jié)果

多目標(biāo)優(yōu)化結(jié)果如圖5所示。圖中空心圓圈為優(yōu)化得到的滿足約束條件的前緣點(diǎn)，它們分布均勻，滿足約束條件，并體現(xiàn)了兩個(gè)優(yōu)化目標(biāo)間的沖突和妥協(xié)，設(shè)計(jì)者可以根據(jù)返回艙再入飛行的實(shí)際需求在前緣點(diǎn)上選擇相應(yīng)的氣動(dòng)外形；實(shí)心圓點(diǎn)為同時(shí)滿足 CD和 L/D最大的理想狀態(tài)（圖中標(biāo)注“Ideal”），由于優(yōu)化目標(biāo)相互沖突，該理想狀態(tài)是達(dá)不到的；“＋”標(biāo)記點(diǎn)表示用于構(gòu)建代理模型的樣本點(diǎn)；前緣上首尾標(biāo)注“Max CD”和“Max L/D”的點(diǎn)為滿足約束前提下，分別以CD和L/D為目標(biāo)的單目標(biāo)優(yōu)化結(jié)果；圖中標(biāo)注“TPF”的點(diǎn)對(duì)應(yīng)于某一典型的優(yōu)化前緣(Typical Pareto Front，TPF)。圖5還給出了樣本點(diǎn)和前緣點(diǎn)上典型狀態(tài)對(duì)應(yīng)的返回艙外形。

圖5 多目標(biāo)優(yōu)化結(jié)果Fig.5 Multi-objective optim ization results

3.4 優(yōu)化設(shè)計(jì)精度與效率分析

在構(gòu)建代理模型時(shí)初始樣本點(diǎn)為30個(gè)，經(jīng)過(guò)3次加點(diǎn)迭代和模型修正后，Kriging模型最大預(yù)測(cè)方差滿足收斂要求，然后隨機(jī)產(chǎn)生4個(gè)測(cè)試點(diǎn)對(duì)建立的代理模型進(jìn)行驗(yàn)證，升阻比、阻力系數(shù)和容積率的最大誤差分別為6.52%、3.58%和1.24%。另外在優(yōu)化前緣上選取3個(gè)標(biāo)注代號(hào)為“Max CD”、“Max L/D”和“TPF”的外形，采用CFD方法計(jì)算對(duì)應(yīng)的氣動(dòng)特性，與采用代理模型預(yù)測(cè)的結(jié)果比較，表3為阻力系數(shù)和升阻比的比較結(jié)果，可見(jiàn)，CFD結(jié)果與代理模型響應(yīng)差別最大不超過(guò)7%。由此可見(jiàn)基于本文介紹的方法構(gòu)建的代理模型具有較好的精度，可以替代真實(shí)物理模型進(jìn)行優(yōu)化設(shè)計(jì)。

表3 典型前緣點(diǎn)特性驗(yàn)證Tab.3 Aerodynam ic character verification of typical optim ization front

本文在優(yōu)化過(guò)程中，總共調(diào)用代理模型20 051次，這意味著如果氣動(dòng)分析直接采用CFD方法的話，總共需要進(jìn)行20 051次CFD數(shù)值計(jì)算，而采用代理模型替代數(shù)值仿真后，代理模型的計(jì)算時(shí)間與CFD計(jì)算時(shí)間相比幾乎可以忽略不計(jì)，整個(gè)優(yōu)化過(guò)程的計(jì)算時(shí)間主要體現(xiàn)在構(gòu)建代理模型時(shí)的40次CFD計(jì)算時(shí)間上，可見(jiàn)引入代理模型后，大約可以減少20 000次CFD分析時(shí)間，這樣高效的分析方法十分適合返回艙總體優(yōu)化設(shè)計(jì)。

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

通過(guò)對(duì)返回艙進(jìn)行氣動(dòng)外形優(yōu)化設(shè)計(jì)，得到了如下結(jié)論：

1）引入代理模型能極大提高氣動(dòng)特性的計(jì)算效率，同時(shí)采用改進(jìn)的EI加點(diǎn)策略能將樣本點(diǎn)的數(shù)目大大減少而不降低代理模型的精度，優(yōu)化計(jì)算結(jié)果表明，代理模型計(jì)算結(jié)果與CFD計(jì)算結(jié)果誤差可以控制在7%以內(nèi)，基于代理模型技術(shù)的優(yōu)化過(guò)程可以節(jié)省絕大部分（95%以上）的計(jì)算開(kāi)銷；

2）以返回艙升阻比和阻力系數(shù)為目標(biāo)，以駐點(diǎn)熱流和容積率為約束條件進(jìn)行多目標(biāo)多約束優(yōu)化，得到了優(yōu)化前緣，為返回艙總體設(shè)計(jì)提供了有益的借鑒。

References)

[1]陳河梧. 飛船返回艙高超聲速氣動(dòng)特性的風(fēng)洞實(shí)驗(yàn)分析 [J]. 航天器工程, 2008, 17(5): 77-81.CHEN Hewu. Wind-tunnel Test Analysis on Hypersonic Aerodynam ic Characteristics of Returnable Module [J]. Spacecraft Engineering, 2008, 17(5): 77-81. (in Chinese)

[2]紀(jì)楚群, 周偉江. 鐘形返回艙空氣動(dòng)力特性數(shù)值模擬 [J]. 航天返回與遙感, 2001, 22(1): 33-37.JI Chuqun, ZHOU Weijiang. A Numerical Simulation of Aerodynamic Characters for Bell Capsule [J]. Spacecraft Recovery and Remote Sensing, 2001, 22(1): 33-37. (in Chinese)

[3]呂俊明, 程曉麗, 王強(qiáng). 火星科學(xué)實(shí)驗(yàn)室氣動(dòng)特性數(shù)值分析 [J]. 力學(xué)與實(shí)踐, 2013, 35(1): 31-35.LV Junm ing, CHENG Xiaoli, WANG Qiang. Numerical Aerodynam ic Analysis of Mars Science Laboratory [J]. Mechanics in Enginering, 2013, 35(1): 31-35. (in Chinese)

[4]李治宇, 楊彥廣, 袁先旭. 基于 Euler方程的返回艙氣動(dòng)外形優(yōu)化設(shè)計(jì)方法研究 [J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2012, 30(5):653-657.LI Zhiyu, YANG Yanguang, YUAN Xianxu. The Study of the Reentry Capsule Shape Optimization Method Based on the Solving of the Euler Equations. [J]. Acta Aerodynam ic Sinica, 2012, 30(5): 653-657. (in Chinese)

[5]John E T, Robert D B. Multi-Objective Hypersonic Entry Aeroshell Shape Optim ization [J]. Journal of Spacecraft and Rockets, 2009, 46(5): 957-966.

[6]John E T, Robert D B. Aerothermodynam ic Shape Optimization of Hypersonic Entry Aeroshells [C]. 13th AIAA/ISSMO Multidisciplinary Analysis Optimization Conference, Fort Worth, Texas, 13-15 September, 2010.

[7]Deb K. Multi-Objective Optimization Using Evolutionary Algorithm [M]. Chichester: John W: ley & Sons, 2002.

[8]Slawomir K, Leifur L. Transonic Airfoil Shape Optim ization Using Variable-Resolution Models and Pressure Distribution Alignment [C]. 29th AIAA Applied Aerodynam ics Conference, Honolulu, Hawaii, 2011.

[9]Li H, Lin P, Yang Z. Modeling and Design Optimization of a Common Aero Vehicle w ith Parameterized Configuration [C].AIAA Modeling and Simulation Technologies Conference, M inneapolis, M innesota, 2012.

[10]Donald R J, Matthias S, William J W. Efficient Global Optim ization of Expensive Black-Box Functions [J]. Journal of Global Optimization, 1998, 13: 455-492.

[11]Artem A D, Edquist K, Shoenenberger M. Influence of the Angle of Attack on the Aerothermodynam ic Environment of the Mars Science Laboratory [J]. AIAA paper, 2006, 2006-3889.

[12]車競(jìng), 唐碩, 何開(kāi)鋒. 類乘波體飛行器氣動(dòng)加熱的工程計(jì)算方法 [J]. 彈道學(xué)報(bào), 2006, 18(4): 93-96.CHE Jing, TANG Shuo, HE Kaifeng. Aerothemo Engineering Method for Quasi-Waverider Vehicle [J]. Journal of Ballistics,2006, 18(4): 93-96. (in Chinese)

[13]Ricardo M P, André R D C, Curran C. Comparison of Surrogate Models in a Multidisciplinary Optimization Framework for Wing Design [J]. AIAA Journal, 2010, 48(5): 995-1006.

[14]Lophaven S N, Nielsen H B, Sondergaard J. DACE a Matlab Kriging Toolbox [R]. Technical University of Denmark, 2002.

[15]Deb K, Pratap A, Agrawal S. A Fast and Elitist Multi-objective Genetic Algorithm: NSGA-II [R]. IEEE Transactions on Evolutionary Computation, 2000, 6(2): 182-197.