序列近似優化方法及其在火箭外形快速設計中的應用*

彭　科，胡　凡，張為華，武澤平

(國防科技大學 航天科學與工程學院, 湖南 長沙　410073)

?

序列近似優化方法及其在火箭外形快速設計中的應用*

彭科，胡凡，張為華，武澤平

(國防科技大學 航天科學與工程學院, 湖南 長沙410073)

摘要：針對序列近似優化方法在代理模型構造與采樣策略等方面的不足，基于采樣點局部密度，引入與局部密度成反比的樣本點影響體積概念，通過優化總影響體積確定徑向基函數最優核寬度，滿足序列近似優化過程不同規模、非均勻樣本條件下的函數近似需要；建立潛在可行域最大距離加點準則，并與潛在最優加點準則結合平衡算法的探索能力與開發能力；建立三步收斂判定準則；構建改進序列近似優化的算法流程。對于Golinski減速器的優化設計問題，算法在目標函數調用42次后便搜索到全局最優解，體現了其良好的全局尋優能力與搜索效率。以“天航二號”火箭為例，建立其外形優化問題數學模型，所提優化方法在調用原始計算模型165次之后便搜索到全局最優解，大大提高了設計效率，同時飛行試驗也表明設計結果滿足要求。

關鍵詞：序列近似優化；代理模型；采樣策略；局部密度；收斂判定準則；“天航二號”火箭；外形優化

序列近似優化(Sequential Approximate Optimization，SAO)方法基于少量初始采樣點構造初始代理模型，采取一定的加點策略更新采樣點，逐步提高代理模型對最優解的近似精度，根據收斂準則終止算法并輸出最優解[1-4]。相比智能優化算法，序列近似優化方法可在保證全局最優的前提下大幅降低原模型調用次數，顯著提高優化效率[1]。代理模型構造與采樣策略是序列近似優化方法的兩項關鍵技術，國內外學者對此進行了大量研究[3-7]。

徑向基插值模型在精度和魯棒性方面皆較為可靠[8]，是序列近似優化方法中被廣泛使用的代理模型，其核寬度的合理確定對近似精度有決定性影響。文獻[6]解決了采樣點較為均勻條件下的徑向基函數核寬度確定問題；文獻[4]提出了一種采樣點個數趨于無窮情況下的非均勻核寬度確定方法；文獻[1]提出了基于采樣點局部密度的核寬度確定方法，但算法中總影響體積這一關鍵參數依賴經驗確定，難以實現最優的近似精度。序列近似優化過程中樣本點數量由少逐漸增多，樣本分布較不均勻，建立不同規模、非均勻樣本的徑向基函數(Radial Basis Function, RBF)核寬度確定方法，對提高序列近似優化過程代理模型精度、提高優化效率具有重要意義。

合理的加點策略是平衡算法探索能力與開發能力以保證全局最優前提下提高算法收斂速率的關鍵因素，目前常用的再采樣策略主要有：EI準則[9]、潛在最優加點準則、最大曲率準則、最小密度加點準則等，現有準則各有特點，但皆難以同時保證優化結果的全局性與收斂速率。

火箭外形設計問題是總體設計過程的重要問題，不斷提高外形設計效率是設計人員的不懈追求。本研究基于樣本點局部密度，引入與局部密度成反比的樣本點影響體積概念，通過優化總影響體積確定徑向基函數最優核寬度，并驗證近似效果；建立潛在可行域最大距離加點準則并與潛在最優加點準則結合，平衡序列近似優化過程的收斂速度與全局搜索精度；建立收斂判定準則；構建改進序列近似優化的算法流程并通過算例分析算法性能；以“天航二號”(TianHang-Ⅱ, TH-Ⅱ)火箭為例，建立其外形優化問題數學模型，采用序列近似優化方法求解并分析結果。

1徑向基函數最優核寬度確定方法

基本徑向基函數的數學模型為：

(1)

取基函數為Gauss函數，即：

(2)

式中，σi為基函數的核寬度。

基本徑向基函數模型是一種非線性代理模型，對線性模型預測誤差較大，故將線性項引入徑向基代理模型，可以增加其線性泛化能力。其數學模型為：

(3)

式中，n為設計空間維數，λj(j=0,1,2,…,n)為線性回歸系數。

式(3)中前兩項為樣本點的線性回歸模型：

(4)

(5)

此時將N個訓練樣本代入式(5)，求解線性方程組即可求得基函數的權系數wi(i=1,2,…,N)。

下面選取線性模型

y=x,-3

(6)

與非線性模型

(7)

對基本徑向基函數與含線性項徑向基函數的近似能力進行比較，結果如圖1、圖2所示，后者對線性函數的近似能力顯著高于前者，對非線性函數的近似能力也略優于前者，體現了其很好的線性和非線性泛化能力。

(a) 線性函數(a) Linear function

(b) 非線性函數(b) Nonlinear function圖1　基本徑向基函數對測試函數的近似Fig.1　Elementary RBF for testing function approximation

式(2)中基函數的核寬度σi對代理模型預測精度有決定性影響，文獻[1]提出了基于樣本點局部密度的核寬度σi計算方法，局部密度定義如式(8)所示。

(a) 線性函數(a) Linear function

(b) 非線性函數(b) Nonlinear function圖2　含線性項徑向基函數對測試函數的近似Fig.2　RBF containing linear terms for testing function approximation

(8)

(9)

(10)

在樣本點分布較密的區域應使每個樣本點的影響適當減小，即核寬度σi應相對較小，而在樣本點分布較為稀疏區域應取相對較大的核寬度σi。 算法步驟如下[1]：

1)根據式(8)計算每個樣本點的密度ρ(xi)；

2)按Vi/Vj=ρj/ρi計算每個樣本點的影響體積之比；

Vt為所有樣本點的影響體積總和，其取值直接決定各樣本點對應核寬度，從而對徑向基代理模型精度產生影響，文獻[1]中Vt取值根據經驗給出。圖3以函數f(x)=xsin(1.5x)為例，給出了10個隨機采樣點情況下不同Vt值對應代理模型與原函數的對比情況；代理模型均方根誤差(Root Mean Square Error, RMSE)隨Vt變化規律如圖4所示。由圖可知，當Vt過小或過大皆不能實現較好的近似效果，確定最優Vt即可得到徑向基函數最優核寬度，對提高代理模型精度具有關鍵意義。

圖3　不同Vt值對應代理模型與原函數的對比Fig.3　Comparison between surrogate models applying different Vtvalues and original model

圖4　不同Vt值對應代理模型均方根誤差Fig.4　RMSE of different Vt

最優Vt值與樣本在設計空間的分布有關，對樣本數量不敏感[1]。將樣本點均分為分區情況基本一致的兩份分別作為近似建模樣本與驗證樣本，以驗證樣本點代理模型預測值均方根誤差最小為目標，采用黃金分割法[10]迭代搜索最優Vt。 近似建模樣本各點局部密度ρ(xi)確定后，黃金分割法搜索最優Vt的單步計算過程如下：

1)根據Vt值計算近似建模樣本點的核寬度ri；

2)求解式(3)中λj與wi得代理模型s(x)；

3)根據代理模型s(x)計算驗證樣本點預測值；

4)驗證樣本點計算預測值均方根誤差。

采用以上方法得圖3對應樣本點最優Vt值為1.17，對應核寬度如圖5所示，代理模型與原函數的對比如圖6所示，二者均方根誤差為0.25，吻合較好，表明給出的徑向基函數最優核寬度確定方法有效。

圖5　最優樣本點核寬度Fig.5　Optimized kernel width of each sampling points

圖6　最優Vt值對應代理模型與原函數的對比Fig.6　Comparison between surrogate model applying optimized Vt and original model

2改進序列近似優化方法

2.1序列近似優化方法基本流程

文獻[1]建立了如圖7所示的序列近似優化方法基本流程，主要步驟為：①將設計空間線性映射到n維單位立方體內，采用拉丁超立方方法得到初始采樣點，并計算采樣點目標函數與約束值；②基于徑向基函數構建目標函數與各約束值的代理模型；③基于代理模型，求解樣本點更新的優化問題；④收斂判斷。

圖7　序列近似優化方法基本流程Fig.7　Flow chart of sequential approximate optimization

2.2采樣點更新準則

一般優化問題可表述為如式(11)所示形式：

(11)

式中，i=0表示無不等式約束，j=0表示無等式約束。

最大距離加點準則是序列近似優化過程常用的加點準則，可最大程度保證算法的探索能力。

令求解域中點x到所有樣本點的最小距離dmin(x)為：

(12)

(13)

為平衡算法的探索能力與開發能力，兼顧序列近似優化過程的收斂速度與全局搜索精度，將潛在最優加點準則與潛在可行域最大距離加點準則同時使用。潛在最優加點將式(14)的最優解作為新的樣本點加入到模型中。

(14)

采樣點更新涉及式(13)、式(14)定義的兩個優化問題，皆采用粒子群優化算法求解。

2.3收斂判定準則

按以下3個步驟判斷算法是否收斂：

(15)

式中，i=N-Sxl+1,N-Sxl+2,…,N-1。

(16)

以上三步皆滿足時則近似序列優化收斂。

改進后的序列近似優化方法流程圖如圖8所示。

圖8　改進序列近似優化算法流程圖Fig.8　Flow chart of enhanced sequential approximate optimization

2.4算例分析

Golinski減速器優化設計問題是美國國家航空航天局提供的一個多學科設計優化測試問題[11]，其目標函數為重量最小，共有7個設計變量：齒輪面寬度x1、牙膜x2、小齒輪牙齒數目x3、軸1長度x4、軸2長度x5以及軸1與軸2的直徑x6和x7。其數學模型表述如式(17)所示。

(17)

式中：

各設計變量取值方位如表1所示。

表1　Golinski減速器設計變量取值范圍

注：圓點為可行解；三角形為不可行解。圖9　Golinski減速器優化問題目標函數與誤差收斂曲線Fig.9　Convergence curve of objective function and relative error for Golinski retarder optimization

所用方法在迭代至第4步，目標函數調用15次之后便搜索到全局最優解附近；迭代至17步，目標函數調用42次后搜索到全局最優解：x=[3.5,0.7,17,7.3,7.715 3,3.350 2,5.286 7]T；代理模型函數值與真實模型函數值的相對誤差僅為1.02×10-9。文獻[12]提出的并行模擬退火單純形算法是目前搜索效率相對很高的智能優化算法，該方法搜索到Golinski減速器優化問題最優解附近約需調用1000次原始模型。本文方法能在保證全局搜索精度的前提下顯著減少模型調用次數，可大幅提高火箭外形設計優化問題求解效率。

3“天航二號”火箭外形快速設計優化

3.1優化問題

根據總體要求與任務特點，“天航二號”火箭氣動布局為“—·×”布局，總長為2000 mm，頭部長度為450 mm，發動機長度取為320 mm，彈翼、空氣舵皆為矩形，空氣舵緊貼儀器艙后端面安裝，如圖10所示。“天航二號”滿載質量為60 kg，對應質心系數0.56，耗盡質量為56 kg，對應質心系數為0.53。

圖10　“天航二號”火箭外形布局Fig.10　 TH-Ⅱ rocket configuration

1)設計變量為：翼距頭部距離、翼弦長、翼展長、舵弦長、舵展長，取值范圍依次為[1000, 1250]、[150, 190]、[350, 500]、[100, 150]、[150, 200](單位mm)。

綜上，“天航二號”火箭外形優化問題數學模型為：

(18)

式中，g4(x)為幾何約束。

3.2優化結果與分析

優化過程氣動特性計算采用工程方法[13]完成，取初始樣本點為100個，迭代至第33步，原始計算模型調用165次之后便搜索到全局最優解，所用機時僅為190.6 s(AMD Athlon 四核處理器，主頻2.8 GHz)，大大提高了設計效率。優化過程目標函數收斂曲線如圖11所示，相對誤差收斂曲線如圖12所示，外形優化結果如圖13所示。需要說明的是，式(18)中約束條件與目標函數計算分三步，依次為：先g1(x)和g2(x)，然后g3(x)和g4(x)，最后f(x)，圖11中三角形對應樣本點因其g2(x)值不滿足約束，無須計算目標函數f(x)值，f(x)為人為給定值。

注：圓點為可行解；三角形與方行為不可行解。圖11　火箭外形優化目標函數收斂曲線Fig.11　Convergence curve of objective function for rocket shape optimization

圖12　火箭外形優化相對誤差收斂曲線Fig.12　Convergence curve of relative error for rocket shape optimization

圖13　“天航二號”火箭外形優化結果Fig.13　TH-II rocket optimized shape

圖14為采用本軟件系統設計優化氣動外形得到的“天航二號”火箭飛行試驗概況，飛行試驗取得圓滿成功，表明本文設計優化得到的“天航二號”試驗火箭靜穩定性、操縱性、配平升阻力特性等性能指標符合設計要求，滿足工程需要。

4結論

研究改進了序列近似優化方法，并將其應用于火箭快速設計，主要研究工作和結論如下：

(a) 總裝后外形(a) Shape after assembling

(b) 點火出架(b) Firing　　(c) 開傘回收(c) Recovery using parachute圖14　“天航二號”火箭飛行試驗概況Fig.14　TH-Ⅱ rocket flight test

1)建立了含線性項的徑向基函數，提高了其對線性函數的近似能力；基于樣本點的局部密度，采用黃金分割法得到了最優總影響體積數值，進而計算最優核寬度，仿真結果表明了該方法的有效性。

2)建立了潛在可行域最大距離加點準則并與潛在最優加點準則結合，兼顧了序列近似優化算法的探索能力與開發能力；提出了序列近似優化方法三步收斂判定準則；構建了改進序列近似優化的算法流程；通過Golinski減速器優化設計問題分析了算法性能，目標函數調用42次后便搜索到全局最優解，體現了良好的全局尋優能力與搜索效率。

3)以“天航二號”火箭為例，建立了其外形優化問題數學模型，采用改進的序列近似優化方法迭代至第33步，原始計算模型調用165次之后便搜索到全局最優解，所用機時僅為190.6 s，大大提高了設計效率。“天航二號”火箭飛行試驗取得圓滿成功，表明設計優化得到的“天航二號”試驗火箭各項氣動性能指標符合設計要求，滿足工程需要。

參考文獻(References)

[1]Wang D H, Wu Z P, Fei Y, et al. Structural design employing a sequential approximation optimization approach[J]. Computers and Structures, 2014, 134:75-87.

[2]Kitayama S, Arakawa M, Yamazaki K. Sequential approximate optimization for discrete design variable problems using radial basis function network[J]. Applied Mathematics and Computation, 2012, 219(8): 4143-4156.

[3]Nakayama H, Arakawa M, Sasaki R. Simulation-based optimization using computational intelligence[J]. Optimization and Engineering, 2002, 3(2): 201-214.

[4]Kitayama S, Arakawa M, Yamazaki K. Sequential approximate optimization using radial basis function network for engineering optimization[J]. Optimization and Engineering, 2011, 12(4): 535-557.

[5]Deng Y M, Lam I C, Tor S B, et al. A CAD-CAE integrated injection molding design system [J].Engineering with Computers, 2002, 18: 80-92.

[6]Kitayama S, Yamazaki K. Simple estimate of the width in Gaussian kernel with adaptive scaling technique[J]. Applied Soft Computing, 2011, 11(8): 4726-4737.

[7]Luo C T, Zhang S L, Wang C, et al. A metamodel-assisted evolutionary algorithm for expensive optimization[J]. Journal of Computational and Applied Mathematics, 2011, 236(5): 759-764.

[8]Jin R, Chen W, Simpson T W. Comparative studies of metamodeling techniques under multiple modeling criteria[J]. Journal of Structural and Multidisciplinary Optimization, 2001, 23(1): 1-13.

[9]Jones D R. A taxonomy of global optimization methods based on response surfaces[J]. Journal of Global Optimization, 2001, 21: 345-383.

[10]謝政， 李建平， 陳摯. 非線性最優化理論與方法[M]. 北京: 高等教育出版社, 2010.

XIE Zheng, LI Jianping, CHEN Zhi. Nonlinear optimization theory and methods[M]. Beijing: Higher Education Press, 2010. (in Chinese)

[11]Padula S L, Alexandrov N, Green L L. MDO test suite at NASA langley research center[C]//Proceedings of 6th AIAA/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, Bellevue, 1996.

[12]Luo Y Z,Tang G J. Parallel simulated annealing using simplex method [J].AIAA Journal, 2006, 44(12): 3143-3146.

[13]嚴恒元. 飛行器氣動特性分析與工程計算[M]. 西安: 西北工業大學出版社, 1990.

YAN Hengyuan. Aircraft aerodynamic characteristics analysis and engineering calculation[M]. Xi′an: Northwestern Polytechnic University Press, 1990. (in Chinese)

Sequential approximate optimization method and its application in rapid design of rocket shape

PENGKe,HUFan,ZHANGWeihua,WUZeping

(College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China)

Abstract：Sequential approximate optimization method has shortcomings in several respects, such as surrogate model establishing and infill strategy at present. Basing on local density of sampling points, the influence volume concept which is inversely proportional to local density was introduced and then the optimal kernel width of radial basis function was obtained by means of total influence volume optimization, thus, the function approximation needs in sequential approximate optimization process under the conditions of different scales and heterogeneous samples were satisfied. Potential feasible region infill strategy was proposed and potential optimal strategy was applied together, both exploration and exploitation capacity of the algorithm were satisfied. Three-step convergence criterion was set up. The algorithm flow process of sequential approximate optimization was constructed. For Golinski reducer optimization problem,the global optimal solution was solved after calculating original model 42 times, which embodied the good global optimization capacity and searching efficiency of the algorithm. Shape optimization mathematical model was established for TH-II rocket, global optimal shape was gained after 165 times of original model calling using the proposed method. The design efficiency was increased greatly and TH-II rocket aerodynamic shape was proved reliable by flight testing.

Key words：sequential approximate optimization; surrogate model; infill strategy; local density; convergence criterion; TH-II rocket; shape optimization

中圖分類號：V421.1；TP391.7

文獻標志碼:A

文章編號：1001-2486(2016)01-129-08

作者簡介：彭科(1989—)，男，四川鄰水人，博士研究生，E-mail：pengke_pk@163.com；張為華(通信作者)，男，教授，博士，博士生導師，E-mail：zhangweihua@nudt.edu.cn

基金項目：國家自然科學基金資助項目(51105368)

*收稿日期：2015-02-03

doi:10.11887/j.cn.201601021

http://journal.nudt.edu.cn