趙 璇，常思江*, 張哲瑋，趙林林

(1. 南京理工大學(xué) 能源與動(dòng)力工程學(xué)院，南京 210094；2. 中國(guó)兵器工業(yè)試驗(yàn)測(cè)試研究院，陜西 華陰 714204； 3.遼沈工業(yè)集團(tuán)有限公司，沈陽(yáng) 110045)

0 引 言

滑翔制導(dǎo)炮彈是在常規(guī)炮彈的基礎(chǔ)上增加制導(dǎo)控制系統(tǒng)而形成的遠(yuǎn)程精確打擊武器，憑借效費(fèi)比高、 反應(yīng)快且使用靈活等優(yōu)點(diǎn)逐漸受到各國(guó)重視。制導(dǎo)炮彈由火炮以一定初速發(fā)射后，尾翼張開保持穩(wěn)定飛行，隨后啟動(dòng)火箭發(fā)動(dòng)機(jī)開始助推，發(fā)動(dòng)機(jī)熄火后炮彈則像普通尾翼彈一樣繼續(xù)在升弧段上爬升。當(dāng)炮彈到達(dá)彈道頂點(diǎn)附近時(shí)，鴨舵張開，同時(shí)制導(dǎo)控制系統(tǒng)開始工作，在降弧段通過(guò)偏轉(zhuǎn)鴨舵使全彈產(chǎn)生一個(gè)正攻角，由此增大升力向前滑翔，實(shí)現(xiàn)增程。由于制導(dǎo)炮彈成本低、 體積小，帶來(lái)了一系列外形設(shè)計(jì)和彈道設(shè)計(jì)的問題。若不能合理地開展參數(shù)設(shè)計(jì)，將直接影響制導(dǎo)炮彈的滑翔和機(jī)動(dòng)能力，甚至難以完成期望的作戰(zhàn)任務(wù)。

氣動(dòng)外形設(shè)計(jì)與彈道優(yōu)化是影響制導(dǎo)炮彈總體性能的關(guān)鍵技術(shù)，國(guó)內(nèi)外已開展了相關(guān)研究。Vasile等通過(guò)DATCOM軟件在給定工況下形成氣動(dòng)數(shù)據(jù)庫(kù)，并利用粒子群算法優(yōu)化了增程炮彈的舵面尺寸； Yang等利用自適應(yīng)遺傳算法以射程最大為目標(biāo)函數(shù)，對(duì)鴨式布局制導(dǎo)火箭彈的鴨舵和尾翼進(jìn)行了外形優(yōu)化； Fowler等在2015年提出了一種基于試驗(yàn)設(shè)計(jì)和Kriging模型的制導(dǎo)炮彈外形設(shè)計(jì)優(yōu)化方法，并通過(guò)遺傳算法進(jìn)行多目標(biāo)尋優(yōu)得到一組Pareto最優(yōu)設(shè)計(jì)集，盡管該文獻(xiàn)考慮了6自由度彈道方程的耦合，但并未涉及方案彈道的優(yōu)化設(shè)計(jì)； Vasile等又在2020年提出了一種基于空氣動(dòng)力學(xué)、 飛行動(dòng)力學(xué)和制導(dǎo)控制等多學(xué)科的設(shè)計(jì)方法，用于評(píng)估遠(yuǎn)程制導(dǎo)炮彈的性能，但是其彈道設(shè)計(jì)是在外形確定后再進(jìn)行的，氣動(dòng)和彈道沒有緊密耦合。國(guó)內(nèi)方面，雷娟棉等設(shè)計(jì)了增程制導(dǎo)炮彈的氣動(dòng)外形，并通過(guò)風(fēng)洞實(shí)驗(yàn)驗(yàn)證了該外形的氣動(dòng)性能，但未采用優(yōu)化方法； 史金光等建立了鴨式制導(dǎo)炮彈的氣動(dòng)外形優(yōu)化數(shù)學(xué)模型，并通過(guò)一定的優(yōu)化設(shè)計(jì)方法得到較佳的外形參數(shù)； 孫世巖等結(jié)合參數(shù)化建模方法、 非線性面元法及遺傳算法，對(duì)制導(dǎo)炮彈的氣動(dòng)外形進(jìn)行了優(yōu)化設(shè)計(jì)； 趙璇等提出了一種基于多可信度代理模型技術(shù)的氣動(dòng)外形快速優(yōu)化方法，對(duì)制導(dǎo)炮彈的鴨舵外形進(jìn)行了優(yōu)化，有效提升了炮彈升阻比。但上述研究均沒有涉及彈道學(xué)科。單繼祥等將氣動(dòng)和彈道耦合建立了代理模型，并利用遺傳算法以最大落地速度為目標(biāo)對(duì)激波針外形進(jìn)行了優(yōu)化，但不涉及彈道的動(dòng)態(tài)優(yōu)化。僅就彈道優(yōu)化方面而言，史金光、 易文俊等分別利用序列二次規(guī)劃法和最大升阻比法對(duì)滑翔增程彈的最遠(yuǎn)滑翔彈道進(jìn)行了優(yōu)化，并對(duì)比分析了兩種方法的各自優(yōu)缺點(diǎn)，但并未研究對(duì)固定目標(biāo)方案的彈道設(shè)計(jì)方法； 陳琦等以飛行時(shí)間為性能指標(biāo)對(duì)固定目標(biāo)的滑翔彈道進(jìn)行了優(yōu)化設(shè)計(jì)，考慮了動(dòng)力學(xué)滯后并引入了虛擬控制量作為優(yōu)化變量，利用Gauss偽譜法求解了彈道優(yōu)化模型。

所謂的氣動(dòng)-彈道綜合優(yōu)化，是通過(guò)某種優(yōu)化方法將氣動(dòng)、 彈道等多個(gè)學(xué)科分析模型進(jìn)行綜合考量，將靜態(tài)參數(shù)優(yōu)化與動(dòng)態(tài)過(guò)程優(yōu)化相結(jié)合，以獲得具有更優(yōu)整體性能的一種技術(shù)或方法。由上述文獻(xiàn)可知，無(wú)論是制導(dǎo)炮彈的外形優(yōu)化(靜態(tài)參數(shù)優(yōu)化)，還是彈道優(yōu)化(動(dòng)態(tài)過(guò)程優(yōu)化)，大多是分別獨(dú)立研究，而關(guān)于制導(dǎo)炮彈氣動(dòng)-彈道綜合優(yōu)化雖近期有所關(guān)注，但研究相對(duì)較少。張海瑞等為了研究某升力體構(gòu)型飛行器在大空域、 寬速域范圍內(nèi)的氣動(dòng)適應(yīng)性問題，提出一種基于氣動(dòng)-彈道一體化模型的優(yōu)化設(shè)計(jì)方法，并在多約束條件下實(shí)現(xiàn)了飛行器氣動(dòng)外形全局優(yōu)化。由文獻(xiàn)[14]可知，氣動(dòng)-彈道綜合優(yōu)化技術(shù)能夠考慮到大空域、 寬速域飛行剖面特征對(duì)氣動(dòng)外形適應(yīng)性的影響，有助于全面提升飛行器氣動(dòng)外形的設(shè)計(jì)精度和水平，因而具有很好的應(yīng)用前景。

對(duì)于制導(dǎo)炮彈，受限于發(fā)射平臺(tái)(如火炮等)，其舵面小、 控制能力有限，但是其炮口初速可達(dá)馬赫數(shù)2.5以上，滑翔速度一般為亞、 跨聲速，飛行高度可達(dá)20 km左右。發(fā)射環(huán)境惡劣、 飛行空域大使得炮彈動(dòng)力學(xué)特性變化較大，且存在諸多不確定擾動(dòng)，這就對(duì)制導(dǎo)炮彈的氣動(dòng)外形和控制參數(shù)設(shè)計(jì)提出了較高的要求。根據(jù)上述文獻(xiàn)，制導(dǎo)炮彈的氣動(dòng)外形優(yōu)化一般是對(duì)其幾何外形參數(shù)的靜態(tài)優(yōu)化，而彈道優(yōu)化則是對(duì)滑翔過(guò)程中舵偏角、 攻角等控制變量的變化規(guī)律進(jìn)行動(dòng)態(tài)優(yōu)化。通常，控制參數(shù)設(shè)計(jì)是建立在炮彈氣動(dòng)參數(shù)的基礎(chǔ)之上，二者具有較強(qiáng)的耦合關(guān)系，若引入氣動(dòng)-彈道綜合優(yōu)化技術(shù)，則更有利于滑翔制導(dǎo)炮彈的方案彈道設(shè)計(jì)及總體性能參數(shù)優(yōu)化。因此，本文以某大口徑滑翔制導(dǎo)炮彈為對(duì)象，綜合利用氣動(dòng)力工程算法、 自適應(yīng)hp-Radau偽譜法，建立以射程和飛行時(shí)間等性能指標(biāo)為目標(biāo)函數(shù)、 舵翼外形參數(shù)和攻角為設(shè)計(jì)變量的氣動(dòng)-彈道綜合優(yōu)化模型，將炮彈外形參數(shù)優(yōu)化與彈道動(dòng)態(tài)優(yōu)化過(guò)程深度融合，并通過(guò)采用基于Kriging模型和組合加點(diǎn)法則的代理優(yōu)化算法，高效獲取該氣動(dòng)-彈道綜合優(yōu)化問題的解。

1 氣動(dòng)-彈道綜合優(yōu)化模型

滑翔制導(dǎo)炮彈的外形參數(shù)決定了炮彈的氣動(dòng)特性，并最終影響彈道特性。為了體現(xiàn)這些外形參數(shù)對(duì)彈道的影響，需要研究氣動(dòng)與彈道學(xué)科之間的耦合關(guān)系，建立各學(xué)科的分析計(jì)算模型。參照文獻(xiàn)[15]給出滑翔制導(dǎo)炮彈的氣動(dòng)、 彈道學(xué)科與炮彈外形耦合關(guān)系的設(shè)計(jì)結(jié)構(gòu)矩陣(Design Structure Matrix, DSM)，如圖1所示，對(duì)角線元素表示各個(gè)學(xué)科，右上方表示數(shù)據(jù)的前饋傳遞，左下方表示數(shù)據(jù)的反饋傳遞。

圖1 滑翔制導(dǎo)炮彈設(shè)計(jì)結(jié)構(gòu)矩陣

本文的計(jì)算模型包括氣動(dòng)分析模型和彈道分析模型。氣動(dòng)分析模型是根據(jù)炮彈的外形參數(shù)、 實(shí)時(shí)飛行狀態(tài)(如馬赫數(shù)、 高度、 攻角)等計(jì)算氣動(dòng)系數(shù)(如零升阻力系數(shù)、 升力系數(shù)導(dǎo)數(shù)等)，為彈道分析提供數(shù)據(jù)支撐。彈道分析模型則需要外形、 氣動(dòng)等學(xué)科計(jì)算結(jié)果與彈道方程相結(jié)合并優(yōu)化設(shè)計(jì)出最優(yōu)控制規(guī)律，完成炮彈全過(guò)程飛行仿真。

1.1 氣動(dòng)分析模型

目前獲得氣動(dòng)參數(shù)的方法包括工程算法、 數(shù)值模擬、 風(fēng)洞試驗(yàn)和靶道自由飛行試驗(yàn)等。在無(wú)法開展風(fēng)洞試驗(yàn)或自由飛行試驗(yàn)的初步設(shè)計(jì)階段，氣動(dòng)工程算法和計(jì)算流體力學(xué)數(shù)值模擬是重要的氣動(dòng)參數(shù)預(yù)測(cè)手段。一般而言，數(shù)值模擬的計(jì)算精度高于工程算法，但數(shù)值計(jì)算耗時(shí)很長(zhǎng)，難以實(shí)現(xiàn)快速優(yōu)化，而工程算法一般僅給出整體和部件的氣動(dòng)力和力矩，計(jì)算速度快且易與其他計(jì)算程序連接進(jìn)行一體化計(jì)算。為兼顧精度和快速性，本文采用文獻(xiàn)[17]中的氣動(dòng)力工程算法來(lái)分析滑翔制導(dǎo)炮彈的氣動(dòng)特性。表1給出了該算法對(duì)不同部件所選用的計(jì)算方法，從而計(jì)算出所需的氣動(dòng)參數(shù)。

表1 氣動(dòng)力工程計(jì)算方法

1.2 彈道分析模型

對(duì)于滑翔制導(dǎo)炮彈，為了實(shí)現(xiàn)遠(yuǎn)程精確打擊，往往需要在發(fā)射前預(yù)先設(shè)計(jì)好一條方案彈道，飛行時(shí)控制機(jī)構(gòu)按照既定的方案彈道操縱舵面調(diào)整飛行姿態(tài)。由于制導(dǎo)炮彈在滑翔段為無(wú)動(dòng)力飛行，機(jī)動(dòng)能力有限，那么方案彈道設(shè)計(jì)的好壞將直接影響制導(dǎo)炮彈的增程效果及作戰(zhàn)性能，因此，滑翔制導(dǎo)炮彈的方案彈道設(shè)計(jì)也是總體設(shè)計(jì)中的一個(gè)重要環(huán)節(jié)?；瓒蔚膹椀涝O(shè)計(jì)實(shí)際上可以看作是一個(gè)單階段的軌跡優(yōu)化問題。本文以最大射程、 飛行時(shí)間短/能量消耗少為優(yōu)化目標(biāo)，在狀態(tài)約束及過(guò)程約束的條件下利用自適應(yīng)hp-Radau偽譜法實(shí)現(xiàn)軌跡優(yōu)化。

1.2.1 炮彈運(yùn)動(dòng)方程

在初步設(shè)計(jì)階段，為便于研究彈道特性，采用質(zhì)點(diǎn)運(yùn)動(dòng)方程作為彈道計(jì)算模型，并作如下假設(shè)： (1)制導(dǎo)炮彈在飛行過(guò)程中是瞬時(shí)平衡的； (2)氣象條件是標(biāo)準(zhǔn)無(wú)風(fēng)雨的，且略去外界干擾的影響； (3)制導(dǎo)炮彈在飛行過(guò)程中無(wú)側(cè)向運(yùn)動(dòng)?；柚茖?dǎo)炮彈在鉛垂平面內(nèi)的質(zhì)心運(yùn)動(dòng)方程組如下：

(1)

式中：為時(shí)間；為射程；為高度；為速度；為彈道傾角；為炮彈質(zhì)量；()為攻角控制指令；和分別為阻力和升力，有

(2)

1.2.2 自適應(yīng)hp-Radau偽譜法

Radau偽譜法的原理是將每個(gè)網(wǎng)格子區(qū)間[-1,]上的狀態(tài)變量和控制變量在一系列LGR(Le-gendre-Gauss-Radau)點(diǎn)上離散化，并以離散點(diǎn)為節(jié)點(diǎn)構(gòu)造Lagrange插值多項(xiàng)式來(lái)逼近狀態(tài)變量和控制變量，把連續(xù)時(shí)間最優(yōu)控制問題轉(zhuǎn)化為具有一系列代數(shù)約束的非線性規(guī)劃問題，再通過(guò)某種非線性規(guī)劃求解器得到數(shù)值最優(yōu)解。

本文采用文獻(xiàn)[20]中的自適應(yīng)hp-Radau偽譜法，將求解最優(yōu)控制問題的h法(局部配點(diǎn)法)和p法(全局偽譜法)結(jié)合起來(lái)，綜合考慮解的精度和計(jì)算效率，以微分代數(shù)約束在特定采樣點(diǎn)上的殘差作為解的近似精度評(píng)估準(zhǔn)則，對(duì)于需要改進(jìn)求解精度的區(qū)間，根據(jù)軌跡曲率通過(guò)細(xì)化網(wǎng)格區(qū)間或增加區(qū)間內(nèi)插值多項(xiàng)式的維數(shù)兩種方式提高求解精度。在整個(gè)迭代求解過(guò)程中，子區(qū)間的個(gè)數(shù)、 區(qū)間寬度和區(qū)間內(nèi)插值多項(xiàng)式的維數(shù)都是變化的。

1.3 綜合優(yōu)化模型

在以上兩個(gè)學(xué)科分析模型的基礎(chǔ)上，建立氣動(dòng)-彈道綜合優(yōu)化模型。

1.3.1 目標(biāo)函數(shù)

任務(wù)需求不同則優(yōu)化目標(biāo)不同。一般而言，為了考察滑翔制導(dǎo)炮彈的滑翔能力，需要對(duì)其最大射程進(jìn)行優(yōu)化研究。另外，當(dāng)攻擊固定目標(biāo)時(shí)，為了提高炮彈的快速反應(yīng)能力，往往要求飛行時(shí)間盡可能的短； 同時(shí)，為了減小被攔截的概率，提高末端攻擊的機(jī)動(dòng)性，也需要制導(dǎo)炮彈在滑翔過(guò)程中盡量節(jié)省能量?；柽^(guò)程中消耗的能量，采用控制量平方對(duì)時(shí)間的積分來(lái)表示。考慮到制導(dǎo)炮彈的控制能力有限，同時(shí)避免彈體失速，設(shè)定為攻角的最大允許值，并在目標(biāo)函數(shù)中進(jìn)行歸一化處理。因此，目標(biāo)函數(shù)為

(3)

式中：，，為權(quán)重系數(shù)，視不同需求分配不同的權(quán)重系數(shù)。

1.3.2 設(shè)計(jì)變量

以某滑翔制導(dǎo)炮彈的鴨舵和尾翼為設(shè)計(jì)對(duì)象，依據(jù)舵面參數(shù)確定的原則，確定后掠梯形舵面的外形參數(shù)，主要包括舵面的半展長(zhǎng)/2、 根弦長(zhǎng)、 舵面前緣至彈頂?shù)木嚯x和前緣后掠角； 矩形尾翼的外形參數(shù)包括尾翼的半展長(zhǎng)/2、 根弦長(zhǎng)以及尾翼前緣至彈頂?shù)木嚯x。以上7個(gè)主要幾何參數(shù)就是外形優(yōu)化的靜態(tài)設(shè)計(jì)變量，如圖2所示。

圖2 某滑翔制導(dǎo)炮彈的外形設(shè)計(jì)變量

與外形參數(shù)不同，滑翔過(guò)程中的控制參數(shù)需要隨時(shí)間連續(xù)變化。本文選取攻角作為控制變量，總的設(shè)計(jì)變量包含7個(gè)靜態(tài)變量和1個(gè)動(dòng)態(tài)變量:

=[2,,,,2,,,]

(4)

1.3.3 約束條件

對(duì)于滑翔制導(dǎo)炮彈，所參加調(diào)優(yōu)的設(shè)計(jì)方案除了滿足設(shè)計(jì)變量的范圍約束，還應(yīng)保證幾何約束、 穩(wěn)定性約束、 控制約束以及滑翔終端約束等9個(gè)條件:

(5)

綜上，滑翔制導(dǎo)炮彈的氣動(dòng)-彈道綜合優(yōu)化數(shù)學(xué)模型可表示為

(6)

式中:()表示式(1)中制導(dǎo)炮彈質(zhì)點(diǎn)彈道模型的6個(gè)方程。

2 代理優(yōu)化算法

由于耦合了彈道分析模型，增大了優(yōu)化迭代的計(jì)算量，增加了計(jì)算時(shí)間。對(duì)此，引入代理模型技術(shù)來(lái)代替優(yōu)化過(guò)程中耗時(shí)的精確數(shù)值分析過(guò)程。所謂代理模型，是指在分析和優(yōu)化過(guò)程中可替代那些復(fù)雜且費(fèi)時(shí)的數(shù)值分析的近似數(shù)學(xué)模型。目前比較常用的代理模型有響應(yīng)面模型、 徑向基函數(shù)模型、 人工神經(jīng)網(wǎng)絡(luò)模型和Kriging模型等。近年來(lái)，代理模型不再是簡(jiǎn)單替代，而是發(fā)展為一種基于歷史數(shù)據(jù)來(lái)驅(qū)動(dòng)新樣本點(diǎn)加入，并逼近局部或全局最優(yōu)解的優(yōu)化機(jī)制，即代理優(yōu)化(Surrogate-Based Optimization，SBO)算法。

本文采用基于Kriging模型的代理優(yōu)化算法作為氣動(dòng)-彈道綜合優(yōu)化模型的優(yōu)化方法，與文獻(xiàn)[14]方法不同的是，在更新代理模型的時(shí)候，依據(jù)兩個(gè)不同的加點(diǎn)準(zhǔn)則能夠一次性增加兩個(gè)樣本點(diǎn)，同時(shí)兼顧了全局搜索和局部挖掘，收斂效果更好，其優(yōu)化流程如圖3所示。首先，采用一定的試驗(yàn)設(shè)計(jì)方法確定一組樣本點(diǎn)，并通過(guò)氣動(dòng)、 彈道分析模型分別計(jì)算其響應(yīng)值構(gòu)成初始樣本集。然后，建立設(shè)計(jì)變量和目標(biāo)函數(shù)之間的Kriging模型來(lái)計(jì)算任意點(diǎn)處的改善期望和預(yù)估響應(yīng)值。其次，利用組合加點(diǎn)準(zhǔn)則篩選出改善期望最大及目標(biāo)函數(shù)預(yù)估值最小的兩個(gè)樣本點(diǎn)。最后，判斷模型預(yù)估最小值處的誤差是否滿足收斂精度，若不滿足則對(duì)新樣本點(diǎn)進(jìn)行氣動(dòng)-彈道綜合計(jì)算，從而構(gòu)成新樣本集并更新代理模型； 若滿足精度要求，則輸出最優(yōu)方案。

由其原理及優(yōu)化流程可以看出，初始樣本的選擇、 代理模型及其訓(xùn)練、 優(yōu)化加點(diǎn)準(zhǔn)則，是代理優(yōu)化算法的三大要素。下面從這三個(gè)方面入手逐一介紹，并對(duì)代理模型的精度進(jìn)行驗(yàn)證。

2.1 試驗(yàn)設(shè)計(jì)

代理優(yōu)化算法的第一步就是設(shè)計(jì)初始樣本點(diǎn)并建立初始代理模型。不同的試驗(yàn)設(shè)計(jì)方法選取的樣本點(diǎn)分布不同，往往導(dǎo)致代理模型的近似精度不同，從而影響優(yōu)化效率。針對(duì)全局優(yōu)化問題，通常采用現(xiàn)代試驗(yàn)設(shè)計(jì)方法選取一組樣本點(diǎn)。一般而言，基于Kriging模型的代理優(yōu)化算法對(duì)初始樣本點(diǎn)數(shù)目并不敏感，初始樣本點(diǎn)數(shù)與后期增加的樣本點(diǎn)數(shù)之比在1∶2以下。本文采用拉丁超立方設(shè)計(jì)在外形參數(shù)的7維設(shè)計(jì)空間中進(jìn)行抽樣。

圖3 氣動(dòng)-彈道綜合優(yōu)化流程圖

2.2 Kriging代理模型

Kriging模型是一種源于地質(zhì)統(tǒng)計(jì)學(xué)的高斯隨機(jī)過(guò)程模型，與其他代理模型顯著不同的是，其不僅能給出對(duì)未知函數(shù)的預(yù)估值，還能給出預(yù)估值的誤差估計(jì)，因此，在各學(xué)科的優(yōu)化設(shè)計(jì)中逐漸受到重視。

Kriging模型假設(shè)未知函數(shù)的響應(yīng)值與設(shè)計(jì)變量之間的真實(shí)關(guān)系為

=()+()

(7)

式中：()為一個(gè)已知的關(guān)于的回歸模型，提供了設(shè)計(jì)空間中的全局趨勢(shì)模型，一般為多項(xiàng)式模型；()為一個(gè)靜態(tài)隨機(jī)過(guò)程，其均值為0，方差為。在設(shè)計(jì)空間的不同位置處，這些隨機(jī)變量具有某種相關(guān)性，該相關(guān)性可用協(xié)方差表示為

[(),()]=(,)

(8)

式中：(,)表示任意兩點(diǎn)的相關(guān)函數(shù)，當(dāng)兩點(diǎn)距離為0時(shí)等于1，當(dāng)距離無(wú)窮大時(shí)等于0，即距離越大相關(guān)性越小。

(9)

(10)

此外，Kriging模型還能給出預(yù)估值的均方差估計(jì):

(1-)/]

(11)

該均方差可以用于指導(dǎo)加入新樣本點(diǎn)來(lái)提高代理模型精度或逼近優(yōu)化問題的最優(yōu)解。

2.3 組合加點(diǎn)準(zhǔn)則

通過(guò)初始樣本點(diǎn)建立代理模型后，即可通過(guò)某種法則加入新的樣本點(diǎn)循環(huán)更新代理模型，直至優(yōu)化問題收斂。這個(gè)過(guò)程中采用傳統(tǒng)優(yōu)化算法確定由加點(diǎn)法則產(chǎn)生的新樣本點(diǎn)的過(guò)程，即為“子優(yōu)化”。主優(yōu)化每一次循環(huán)加點(diǎn)都需要進(jìn)行一次完整的子優(yōu)化迭代收斂過(guò)程。在子優(yōu)化中只需調(diào)用代理模型預(yù)估而無(wú)需精確數(shù)值分析，可以快速計(jì)算。目前，基于Kriging模型的代理優(yōu)化算法已經(jīng)發(fā)展了多種加點(diǎn)準(zhǔn)則，其中最大改善期望(Expected Improvement，EI)準(zhǔn)則因?yàn)閮?yōu)異的全局搜索特性而得到廣泛應(yīng)用。

(12)

當(dāng)標(biāo)準(zhǔn)差()>0時(shí)，改善期望計(jì)算如下:

(13)

式中:和分別為標(biāo)準(zhǔn)正態(tài)累計(jì)分布函數(shù)和標(biāo)準(zhǔn)正態(tài)分布概率密度函數(shù)。當(dāng)標(biāo)準(zhǔn)差()=0時(shí)，[()]=0。

通過(guò)求解最大化EI值的子優(yōu)化問題：

max[()]

s.t.≤≤

(14)

從而可得到新的樣本點(diǎn)。

理論上，EI方法是一種全局優(yōu)化算法，但在優(yōu)化后期收斂較慢。因此，為彌補(bǔ)其缺陷，參考了文獻(xiàn)[24]引入最小化代理模型預(yù)測(cè)(Minimize Surrogate Prediction，MSP)準(zhǔn)則，在增加改善期望最大處樣本點(diǎn)的基礎(chǔ)上，再搜索當(dāng)前代理模型中預(yù)估值最小的樣本點(diǎn)一并加入樣本集中，即

(15)

采用遺傳算法求解式(14)～(15)的子優(yōu)化問題，再對(duì)最優(yōu)解進(jìn)行精確數(shù)值分析，從而實(shí)現(xiàn)在一步迭代中添加兩個(gè)新樣本，直至整個(gè)優(yōu)化過(guò)程收斂。研究表明，該組合加點(diǎn)方法對(duì)于氣動(dòng)優(yōu)化設(shè)計(jì)而言非常實(shí)用。

2.4 模型精度驗(yàn)證

為了檢驗(yàn)基于EI+MSP組合加點(diǎn)準(zhǔn)則的代理優(yōu)化算法的性能，利用某二維函數(shù)對(duì)該算法進(jìn)行測(cè)試； 此外，為了驗(yàn)證Kriging模型對(duì)氣動(dòng)-彈道綜合計(jì)算問題的精度，也進(jìn)行了相應(yīng)測(cè)試以充分檢驗(yàn)其擬合能力。

2.4.1 函數(shù)測(cè)試

以改進(jìn)的雙變量Branin函數(shù)為例:

(16)

式中：∈[-5, 10];∈[0, 15]。

圖4展示了將其變量縮放至[0，1]區(qū)間后的三維全貌。該函數(shù)在右下角有兩個(gè)局部最優(yōu)解，在左上角有一個(gè)全局最優(yōu)解。初始樣本點(diǎn)為5個(gè)，分別采用EI+MSP組合加點(diǎn)和單獨(dú)EI加點(diǎn)兩種代理優(yōu)化算法來(lái)求解該函數(shù)的最小值，其收斂過(guò)程如圖5所示(為方便比較，對(duì)目標(biāo)函數(shù)值取對(duì)數(shù))。

圖4 Branin函數(shù)三維圖

圖5 收斂歷程比較

可以看出，利用EI+MSP組合加點(diǎn)準(zhǔn)則只需12次循環(huán)就可收斂至全局最優(yōu)解附近，而EI準(zhǔn)則需要13次循環(huán)才可以將最小值降低到2以下。在優(yōu)化后期，由于EI準(zhǔn)則具有較強(qiáng)的全局搜索能力，使得其很難在最優(yōu)解附近進(jìn)行充分挖掘，加點(diǎn)過(guò)程基本變?yōu)殡S機(jī)搜索。因此，在EI準(zhǔn)則的基礎(chǔ)上增加MSP準(zhǔn)則可以很好地彌補(bǔ)EI加點(diǎn)難以收斂的局限性，使得全局搜索和局部挖掘得以平衡，達(dá)到更好的效果。

2.4.2 氣動(dòng)彈道計(jì)算測(cè)試

為了檢驗(yàn)Kriging模型對(duì)不同外形方案下彈道計(jì)算問題的擬合能力，分別采用100個(gè)、 60個(gè)、 30個(gè)樣本點(diǎn)，構(gòu)建滑翔制導(dǎo)炮彈外形參數(shù)和最大射程之間的Kriging模型，并利用這三個(gè)模型對(duì)另50個(gè)樣本點(diǎn)(與構(gòu)建模型所用的樣本點(diǎn)不重復(fù))進(jìn)行預(yù)測(cè)。樣本點(diǎn)均采用拉丁超立方設(shè)計(jì)抽樣獲得。同時(shí)，對(duì)這50個(gè)預(yù)測(cè)樣本點(diǎn)采用第1節(jié)的氣動(dòng)-彈道分析模型進(jìn)行精確數(shù)值計(jì)算，其計(jì)算結(jié)果作為理論計(jì)算值與代理模型的預(yù)測(cè)值進(jìn)行比較。

理論計(jì)算值與Kriging模型預(yù)測(cè)值的對(duì)比如圖6所示，實(shí)線表示計(jì)算值與模型預(yù)測(cè)值完全相等的線，點(diǎn)距離實(shí)線越近表示計(jì)算值與預(yù)測(cè)值的偏差越小。由圖6可以看出，除個(gè)別樣本點(diǎn)預(yù)測(cè)值與計(jì)算值偏離較大以外，大部分樣本點(diǎn)的Kriging模型預(yù)測(cè)結(jié)果與精確計(jì)算結(jié)果基本一致。當(dāng)構(gòu)建Kriging模型的初始樣本點(diǎn)個(gè)數(shù)為100時(shí)，模型預(yù)測(cè)結(jié)果的平均誤差為0.885%，當(dāng)初始樣本點(diǎn)分別為60個(gè)和30個(gè)時(shí)，平均誤差分別為1.492%和1.850%，均小于2%。經(jīng)過(guò)分析比較發(fā)現(xiàn)，若個(gè)別誤差較大的點(diǎn)距離最優(yōu)解較遠(yuǎn)，則對(duì)最終的優(yōu)化結(jié)果影響較小，且增加相應(yīng)的樣本點(diǎn)數(shù)目可以有效提高Kriging模型的局部預(yù)測(cè)精度。因此，Kriging模型在一定的精度允許范圍內(nèi)能夠代替復(fù)雜的氣動(dòng)-彈道計(jì)算程序。

圖6 預(yù)測(cè)值與計(jì)算值對(duì)比

3 算例分析

依據(jù)不同的任務(wù)需求，第1.3節(jié)所建立的滑翔制導(dǎo)炮彈氣動(dòng)-彈道綜合優(yōu)化數(shù)學(xué)模型也有著不同的目標(biāo)函數(shù)和約束條件。本文分別以射程最大、 飛行時(shí)間短/控制能量消耗小為目標(biāo)開展氣動(dòng)-彈道綜合優(yōu)化設(shè)計(jì)。

以某155 mm滑翔制導(dǎo)炮彈為例，炮彈質(zhì)量40.4 kg，總長(zhǎng)995 mm，頭部長(zhǎng)405 mm，其鴨舵和尾翼外形參數(shù)的變化范圍如表2所示，構(gòu)成了一個(gè)7維設(shè)計(jì)空間。

自適應(yīng)hp-Radau偽譜法優(yōu)化彈道時(shí)，假設(shè)滑翔啟控點(diǎn)為彈道頂點(diǎn)，高度為20 km，彈道傾角為0°，速度為300 m/s。

優(yōu)化是否收斂可根據(jù)代理模型在最優(yōu)解附近的近似精度來(lái)判斷：

(17)

表2 設(shè)計(jì)變量的取值范圍

3.1 單目標(biāo)優(yōu)化問題

單目標(biāo)優(yōu)化可以看作是多目標(biāo)優(yōu)化中的特殊情況。當(dāng)以射程最大為優(yōu)化目標(biāo)時(shí)，取=1，=0，=0，即目標(biāo)函數(shù)為

=-()

(18)

約束條件如下：

(19)

首先，采用拉丁超立方設(shè)計(jì)在表1所示的設(shè)計(jì)空間中抽取樣本數(shù)20個(gè)； 然后，對(duì)每一個(gè)外形樣本采用第1節(jié)所述氣動(dòng)分析模型來(lái)計(jì)算其氣動(dòng)特性，再通過(guò)彈道分析模型進(jìn)行方案彈道的設(shè)計(jì)。經(jīng)過(guò)對(duì)20個(gè)樣本點(diǎn)的學(xué)習(xí)構(gòu)建出制導(dǎo)炮彈外形參數(shù)和最大射程之間的Kriging模型，并利用組合加點(diǎn)準(zhǔn)則不斷增加樣本點(diǎn)，提高模型精度并搜尋最優(yōu)解。

通過(guò)求解氣動(dòng)-彈道綜合優(yōu)化模型，收斂過(guò)程如圖7所示。在優(yōu)化過(guò)程中，每次循環(huán)給Kriging模型增加改善期望最大處和模型預(yù)測(cè)最小處兩個(gè)樣本點(diǎn)。當(dāng)模型預(yù)測(cè)最小值的相對(duì)誤差小于1.0×10時(shí)，輸出為最優(yōu)解。綜合優(yōu)化方案(Opt-1方案)的外形如圖8所示，表3為綜合優(yōu)化結(jié)果。

圖7 綜合優(yōu)化收斂過(guò)程(Opt-1方案)

圖8 綜合優(yōu)化方案外形(Opt-1方案)

表3 綜合優(yōu)化結(jié)果(Opt-1方案)

作為對(duì)比，對(duì)基準(zhǔn)外形也以射程最大為目標(biāo)函數(shù)進(jìn)行彈道優(yōu)化。由表3可以看出，和基準(zhǔn)方案相比，綜合優(yōu)化方案Opt-1的鴨舵和尾翼面積均有所增大，且兩者至彈頂?shù)木嚯x都縮短了，從而可以提高升力并減小靜穩(wěn)定度。與此同時(shí)，Opt-1方案的鴨舵和尾翼的展弦比更大，可以有效提升炮彈的升阻比，符合“理想狀況下，升阻比越大，射程越遠(yuǎn)”的規(guī)律。自適應(yīng)hp-Radau偽譜法對(duì)Opt-1方案的彈道優(yōu)化結(jié)果如圖9所示。

圖9(a)為制導(dǎo)炮彈的飛行軌跡圖，制導(dǎo)炮彈經(jīng)過(guò)彈道最高點(diǎn)后，近乎以直線滑翔。其中，Opt-1方案的最大射程可以達(dá)到60.9 km，基準(zhǔn)方案的最大射程為44.2 km，通過(guò)綜合優(yōu)化使得滑翔距離增加了37.8%。圖9(b)為控制變量攻角的變化曲線。剛開始滑翔時(shí)攻角逐漸增加，這是因?yàn)榕趶椝幐叨容^高，空氣密度小，大攻角下可以產(chǎn)生足夠的升力。而Opt-1方案的攻角幅值小于基準(zhǔn)方案，攻角曲線變化更為平緩?？紤]到需用的攻角越小，炮彈的控制裕量就越大，Opt-1方案有著更好的控制效果。在彈道末段，攻角逐漸增加至最大值，這是為了延長(zhǎng)落地時(shí)間，提高滑翔距離。此外，在控制量變化較為劇烈的區(qū)域，節(jié)點(diǎn)比較密集，而在控制量變化平緩的區(qū)域，節(jié)點(diǎn)較稀疏。這主要是hp自適應(yīng)配點(diǎn)的作用，通過(guò)采樣點(diǎn)曲率來(lái)判斷是否增加配點(diǎn)數(shù)或細(xì)化區(qū)間。這樣做可以兼顧求解精度和效率，提高Radau偽譜法的收斂速度。圖9(c)為滑翔過(guò)程中彈道傾角的變化曲線。相較于基準(zhǔn)方案，Opt-1方案的彈道傾角幅值較小，且變化均勻，這也可以從圖9(a)中直觀地看出，Opt-1方案的滑翔軌跡更為平直。圖9(d)為速度變化曲線。在開始滑翔后，飛行速度逐漸減小，但由于重力分量的作用，速度又略微回升，最終緩慢下降，且Opt-1方案的落速小于基準(zhǔn)方案。這表明射程的提高在一定程度上是以犧牲制導(dǎo)炮彈的動(dòng)能為代價(jià)的。

圖9 彈道優(yōu)化結(jié)果(Opt-1方案)

3.2 多目標(biāo)優(yōu)化問題

當(dāng)攻擊固定目標(biāo)時(shí)，為了提高末端機(jī)動(dòng)性，縮短飛行時(shí)間的同時(shí)往往需要兼顧能量消耗，在滑翔過(guò)程中應(yīng)盡量節(jié)省能量。因此，目標(biāo)函數(shù)取為時(shí)間和能量的函數(shù)：

(20)

式中：=0;=06;=04。

約束條件如下：

(21)

采用相同的優(yōu)化流程，收斂歷程如圖10所示。觀察圖7和圖10可以發(fā)現(xiàn)，在第四次循環(huán)時(shí)，相對(duì)誤差都存在跳躍式增大，分析認(rèn)為這是由于組合加點(diǎn)準(zhǔn)則的全局搜索特性導(dǎo)致的，在第四次循環(huán)加點(diǎn)時(shí)從某個(gè)局部最優(yōu)區(qū)域搜索進(jìn)入了另一個(gè)區(qū)域。圖11為多目標(biāo)綜合優(yōu)化方案(Opt-2方案)對(duì)應(yīng)的外形圖。對(duì)基準(zhǔn)外形也進(jìn)行同樣多目標(biāo)的彈道優(yōu)化作為對(duì)比，優(yōu)化結(jié)果如表4所示。

圖10 綜合優(yōu)化收斂過(guò)程(Opt-2方案)

圖11 綜合優(yōu)化方案外形(Opt-2方案)

表4 綜合優(yōu)化結(jié)果(Opt-2方案)

不難看出，和基準(zhǔn)方案相比，Opt-2方案的鴨舵前緣后掠角增大，導(dǎo)致鴨舵的展弦比增大為5.42，鴨舵的根梢比提高至2.12。Opt-2方案尾翼的展弦比降低為2.74，其尾翼面積接近面積上限。Opt-2方案在馬赫數(shù)0.5～0.9下的靜穩(wěn)定度約為4%～6%，符合設(shè)計(jì)要求。其飛行時(shí)間為112.33 s，略小于基準(zhǔn)方案； 同時(shí)，其能量消耗降低了約46%。Opt-2方案的彈道優(yōu)化結(jié)果如圖12所示。

從圖12(a)可以看出，Opt-2方案的彈道曲線比基準(zhǔn)方案更接近一條直線。圖12(b)為攻角隨時(shí)間變化的規(guī)律，在整個(gè)滑翔過(guò)程中，Opt-2方案的攻角始終小于基準(zhǔn)方案，攻角最大幅值也低于基準(zhǔn)方案且更加平滑，這意味著優(yōu)化后的制導(dǎo)炮彈減少了能量消耗，有著更大的控制裕量。而末端攻角的快速衰減是由于對(duì)末端彈道傾角進(jìn)行約束所致，如圖12(c)所示，最終彈道傾角降低至-40°，滿足約束條件。觀察圖12(d)發(fā)現(xiàn)，由于滑翔初始攻角較小，受重力的作用，飛行速度快速上升，隨著攻角增大，速度增加放緩并在45 s左右開始下降。Opt-2方案的速度在降低的過(guò)程中始終大于基準(zhǔn)方案，最終落地速度為267 m/s，因此綜合優(yōu)化后有著更好的攻擊效果。

圖12 彈道優(yōu)化結(jié)果(Opt-2方案)

4 結(jié) 論

針對(duì)滑翔制導(dǎo)炮彈的氣動(dòng)、 彈道設(shè)計(jì)問題，本文研究了一種基于Kriging模型的氣動(dòng)-彈道綜合優(yōu)化方法，綜合考慮了外形參數(shù)的靜態(tài)優(yōu)化和攻角控制規(guī)律的動(dòng)態(tài)優(yōu)化。以某155 mm滑翔制導(dǎo)炮彈的鴨舵和尾翼為優(yōu)化設(shè)計(jì)對(duì)象，采用工程算法進(jìn)行氣動(dòng)分析以及自適應(yīng)hp-Radau偽譜法進(jìn)行彈道分析，并通過(guò)基于Kriging模型的代理優(yōu)化算法求解最優(yōu)外形和最優(yōu)彈道。分別以射程最大、 飛行時(shí)間短/能量消耗小為目標(biāo)函數(shù)，通過(guò)仿真驗(yàn)證了該氣動(dòng)-彈道綜合優(yōu)化方法的可行性。

結(jié)果表明： (1)本文所構(gòu)建的Kriging模型能很好地代替氣動(dòng)-彈道計(jì)算過(guò)程，對(duì)不同方案最大射程的預(yù)測(cè)精度較高?；贙riging模型和EI+MSP組合加點(diǎn)準(zhǔn)則的代理優(yōu)化算法收斂效果較好，可滿足約束條件。(2)和基準(zhǔn)方案相比，經(jīng)氣動(dòng)-彈道綜合優(yōu)化后的Opt-1方案使射程提高了37.8%，Opt-2方案在攻擊固定目標(biāo)時(shí)能有效減少飛行時(shí)間并使能量消耗降低了46%。本文研究可為滑翔制導(dǎo)炮彈的設(shè)計(jì)與研制提供一定的理論參考。