變循環(huán)發(fā)動(dòng)機(jī)數(shù)學(xué)模型求解及有效性研究

吳正佳 孟榮華 余 剛 何海洋 張 屹

（三峽大學(xué) 機(jī)械與動(dòng)力學(xué)院，湖北 宜昌 443002）

變循環(huán)發(fā)動(dòng)機(jī)具有大涵道比性能優(yōu)勢(shì)，可以在不同航段分別實(shí)現(xiàn)高推力和低油耗，因此受到了各航空強(qiáng)國(guó)的廣泛重視，是當(dāng)今航空發(fā)動(dòng)機(jī)的重要研究方向［1］.目前，變循環(huán)發(fā)動(dòng)機(jī)的性能研究主要采用實(shí)驗(yàn)法和計(jì)算機(jī)仿真法.實(shí)驗(yàn)法需要研制復(fù)雜的設(shè)備、投入巨額的資金，因此不可能經(jīng)常采用.而隨著計(jì)算機(jī)仿真技術(shù)的發(fā)展及發(fā)動(dòng)機(jī)數(shù)學(xué)模型研究的不斷深入，計(jì)算機(jī)仿真法可在一定程度上彌補(bǔ)實(shí)驗(yàn)法的不足.計(jì)算機(jī)仿真法的核心問(wèn)題是變循環(huán)發(fā)動(dòng)機(jī)數(shù)學(xué)模型的建立與求解.目前數(shù)學(xué)模型的建立多采用部件法［2］；數(shù)學(xué)模型的求解通常采用牛頓-拉夫遜法［3］、Broyden法［4］、遺傳算法［5-7］等.

國(guó)內(nèi)關(guān)于變循環(huán)發(fā)動(dòng)機(jī)的研究工作主要集中在變循環(huán)發(fā)動(dòng)機(jī)數(shù)學(xué)模型的建立［2，8］，而對(duì)變循環(huán)發(fā)動(dòng)機(jī)數(shù)學(xué)模型的求解算法研究相對(duì)較少.現(xiàn)有的通用求解算法往往有其自身的局限性，這些算法在變循環(huán)發(fā)動(dòng)機(jī)模型求解中的有效性及適用性未見(jiàn)討論.本文采用部件法建立了變循環(huán)發(fā)動(dòng)機(jī)的多維非線性隱式方程組模型，選用了遺傳算法對(duì)變循環(huán)發(fā)動(dòng)機(jī)數(shù)學(xué)模型進(jìn)行求解，提出了算法的有效性評(píng)價(jià)指標(biāo)：初值敏感性、計(jì)算效率、收斂性和穩(wěn)定性，對(duì)遺傳算法進(jìn)行評(píng)價(jià)，為變循環(huán)發(fā)動(dòng)機(jī)模型求解算法的選擇與設(shè)計(jì)提供依據(jù).

1 變循環(huán)發(fā)動(dòng)機(jī)數(shù)學(xué)模型

變循環(huán)發(fā)動(dòng)機(jī)的工作原理如圖1所示.當(dāng)變循環(huán)發(fā)動(dòng)機(jī)處于穩(wěn)定工作狀態(tài)時(shí)，其各部件狀態(tài)互相制約，需滿足共同工作條件：低壓軸功率平衡、高壓軸功率平衡、高壓渦輪進(jìn)口截面流量平衡、低壓渦輪進(jìn)口截面流量平衡、后混合器靜壓平衡、尾噴管面積平衡、風(fēng)扇出口流量平衡［6］.

圖1 變循環(huán)發(fā)動(dòng)機(jī)工作原理圖

由圖1可知，變循環(huán)發(fā)動(dòng)機(jī)為雙轉(zhuǎn)子發(fā)動(dòng)機(jī)，風(fēng)扇與低壓渦輪相連，CDFS、高壓壓氣機(jī)與高壓渦輪相連.

按氣流流過(guò)的順序，從進(jìn)氣道入口至尾噴管出口截面對(duì)變循環(huán)發(fā)動(dòng)機(jī)各部件模型模擬，完成所有部件氣動(dòng)熱力計(jì)算，推導(dǎo)出7個(gè)共同工作方程，根據(jù)上述方程即可建立變循環(huán)發(fā)動(dòng)機(jī)的共同工作方程組模型：

當(dāng)發(fā)動(dòng)機(jī)在某一設(shè)定狀態(tài)下穩(wěn)定工作，會(huì)存在一組確定的獨(dú)立變量滿足隱式方程組.但由于氣動(dòng)熱力計(jì)算、部件特性、數(shù)值計(jì)算等存在誤差，很難得到恰好完全滿足7個(gè)平衡方程的獨(dú)立變量，即共同工作方程組等式很難成立.因此，需要采用數(shù)值算法，使共同工作方程等式右邊接近0，得到一定精度范圍的數(shù)值解.本文將7個(gè)方程轉(zhuǎn)化為殘差方程組.由于7個(gè)方程的數(shù)量級(jí)不同，為降低殘差數(shù)量級(jí)對(duì)方程求解過(guò)程的影響，需要把7個(gè)方程殘差進(jìn)行歸一化.

其中，E1，E2，E3，…，E7為歸一化殘差.

2 遺傳算法求解變循環(huán)發(fā)動(dòng)機(jī)模型

變循環(huán)發(fā)動(dòng)機(jī)計(jì)算機(jī)仿真技術(shù)的核心問(wèn)題最終歸結(jié)為非線性方程組求解.遺傳算法（genetic algorithm，GA）是一種基于自然群體遺傳演化機(jī)制的全局優(yōu)化算法.該算法具有無(wú)需梯度信息、不要求顯式的函數(shù)形式、對(duì)初猜值無(wú)要求的特點(diǎn)，尤其適合于大規(guī)模、高度非線性及無(wú)解析表達(dá)式的目標(biāo)函數(shù)優(yōu)化問(wèn)題，但是運(yùn)行效率不佳.一般通用求解算法往往有其自身的局限性，在變循環(huán)發(fā)動(dòng)機(jī)模型求解中的有效性及適用性未見(jiàn)討論.本文主要采用GA對(duì)上述模型進(jìn)行求解，將求解結(jié)果與牛頓-拉夫遜求解進(jìn)行對(duì)比，并進(jìn)一步討論GA的有效性.

2.1 遺傳算法目標(biāo)函數(shù)設(shè)計(jì)

本文以2013年數(shù)學(xué)建模A題的問(wèn)題二為例，變循環(huán)發(fā)動(dòng)機(jī)的工作狀態(tài)由式（2）中的獨(dú)立變量控制，各獨(dú)立變量及飛行環(huán)境條件（導(dǎo)葉角度αL＝αCDFS＝αH＝αCH＝0°，低壓轉(zhuǎn)速nL＝0.85）下的約束見(jiàn)表1.

表1 獨(dú)立變量及其約束

由式（2）及表1可得變循環(huán)發(fā)動(dòng)機(jī)多維非線性隱式方程組模型為

約束條件

式中，λ為速度系數(shù)，λ＞0，π（λ）、q（λ）為氣動(dòng)函數(shù).

采用GA求解非線性方程組的關(guān)鍵是把方程組求解問(wèn)題轉(zhuǎn)化為如下函數(shù)的優(yōu)化問(wèn)題：

式中，x＝（x1，x2，…，x7），E 取極小值時(shí)的x 即為數(shù)學(xué)模型的解.

2.2 遺傳算法求解發(fā)動(dòng)機(jī)數(shù)學(xué)模型的過(guò)程

Step1：初始種群的生成.

首先進(jìn)行編碼，二進(jìn)制編碼方式利于交叉和變異操作，且可以滿足精度要求，故采用二進(jìn)制編碼；確定7個(gè)自變量x的搜索范圍，0＜x1＜1，0＜x2＜1，0＜x3＜1，800＜x4＜1800，0＜x5＜1，0＜x6＜1，0＜x7＜1；為了保證平衡方程的解滿足一定的精度要求，經(jīng)過(guò)多次調(diào)試與仿真計(jì)算，最終確定編碼長(zhǎng)度為10；根據(jù)編碼規(guī)則隨機(jī)產(chǎn)生N個(gè)初始編碼串，每個(gè)編碼串就是一個(gè)個(gè)體，N個(gè)個(gè)體在一起構(gòu)成初始種群.GA以這N個(gè)編碼串作為起點(diǎn)開(kāi)始迭代.設(shè)置種群大小N＝20，進(jìn)化代數(shù)T＝50.

Step2：適應(yīng)度值評(píng)價(jià).

GA是一種搜索優(yōu)化算法，采用GA求解非線性方程組的關(guān)鍵是如何把方程組求解問(wèn)題轉(zhuǎn)化為函數(shù)優(yōu)化問(wèn)題，即構(gòu)造可以表明個(gè)體或解的優(yōu)劣性的適應(yīng)度函數(shù)，本文中適應(yīng)度函數(shù)與目標(biāo)函數(shù)相同，將λi均取值為1／7.根據(jù)適應(yīng)度函數(shù)，計(jì)算群體P（t）中個(gè)體的適應(yīng)度.適應(yīng)度函數(shù)為

其中，x＝（x1，x2，…，x7），顯然，若數(shù)學(xué)模型有解，則適應(yīng)度函數(shù)（5）的最小值為0，如果求得的適應(yīng)度函數(shù)的值越接近于0，那么對(duì)應(yīng)的解就越精確.

Step3：選擇操作.

采用競(jìng)標(biāo)賽法對(duì)產(chǎn)生的個(gè)體進(jìn)行選擇操作.隨機(jī)地從種群中挑選3個(gè)個(gè)體，然后將最好的個(gè)體選擇做父本個(gè)體，重復(fù)進(jìn)行這個(gè)過(guò)程，完成個(gè)體的選擇.在交叉和變異操作完成后，采用精英選擇法選擇20個(gè)最優(yōu)個(gè)體，作為下一代的父本.

Step4：交叉操作.

采用二進(jìn)制編碼中最常見(jiàn)的單點(diǎn)交叉方式，交叉點(diǎn)的范圍為［1，19］，在該點(diǎn)為分界相互交換變量.按照交叉概率Pc，選擇種群中個(gè)體進(jìn)行交叉，本文中交叉概率為0.8.

Step5：變異操作.

變異算子以較小的概率對(duì)個(gè)體的某個(gè)或者某些位進(jìn)行一些特殊的改變，進(jìn)而生成新的個(gè)體，本文的變異操作采用單點(diǎn)變異方法，變異概率Pm為0.1.群體P（t）經(jīng)過(guò)選擇、交叉、變異運(yùn)算后得到下一代群體P（t＋1）.

Step6：終止條件判斷.

若t≤T，則t←t＋1，轉(zhuǎn)到步驟（2）；若t＞T，則以進(jìn)化過(guò)程中所得到的具有最大適應(yīng)度的個(gè)體作為最優(yōu)解輸出，終止運(yùn)算.

2.3 遺傳算法求解發(fā)動(dòng)機(jī)模型的流程圖

圖2 遺傳算法流程圖

3 仿真與結(jié)果分析

3.1 仿真結(jié)果

將表1的約束條件代入GA的Matlab［9-10］程序，運(yùn)行GA程序求解變循環(huán)發(fā)動(dòng)機(jī)數(shù)學(xué)模型，運(yùn)行時(shí)間為3598.22s，得到相對(duì)誤差0.4454時(shí)目標(biāo)函數(shù)取得最小值，求得各自變量的值見(jiàn)表2.

表2 遺傳算法求解多維非線性隱式方程組結(jié)果

3.2 算法有效性評(píng)價(jià)

本文選用了GA對(duì)變循環(huán)發(fā)動(dòng)機(jī)數(shù)學(xué)模型進(jìn)行了求解，為評(píng)價(jià)該算法的有效性，特提出算法有效性評(píng)價(jià)指標(biāo)：初值敏感性、計(jì)算效率、收斂性和穩(wěn)定性.通過(guò)上述運(yùn)行結(jié)果分析對(duì)遺傳算法進(jìn)行有效性評(píng)價(jià).

初值敏感性評(píng)價(jià)：GA求解僅需要設(shè)定搜索范圍，不需選定初值，多次運(yùn)行均可獲得較滿意的結(jié)果；而牛頓-拉夫遜本身對(duì)初始值的選擇依賴性較大，始值設(shè)計(jì)不合理時(shí)，容易形成局部收斂，產(chǎn)生局部最優(yōu)解，甚至無(wú)解.GA運(yùn)行與各種對(duì)于牛頓-拉夫遜的研究對(duì)比表明，GA在對(duì)初值的選擇上要求較低.

計(jì)算效率評(píng)價(jià)：程序在同樣配置的計(jì)算機(jī)上運(yùn)行，配置為：CPU 2.9GHz，RAM 4.0G.如果選用牛頓-拉夫遜運(yùn)行效率較快，一般運(yùn)行時(shí)間是90～100s，而GA相對(duì)運(yùn)行較慢，一般需要3600s左右，兩種算法效率相差比較明顯.主要是因?yàn)镚A在開(kāi)始進(jìn)化的時(shí)候，種群的適應(yīng)度函數(shù)曲線收斂速度很快，隨著遺傳代數(shù)的不斷增加，種群朝著最優(yōu)個(gè)體適應(yīng)度收斂的速度會(huì)迅速下降，導(dǎo)致在一定遺傳代數(shù)之后，種群的適應(yīng)度變化很緩慢，需要很長(zhǎng)時(shí)間才能進(jìn)化到最優(yōu)個(gè)體，計(jì)算效率在后期很不理想.

收斂性與穩(wěn)定性評(píng)價(jià)：GA具有全局收斂性，對(duì)初始點(diǎn)沒(méi)有任何要求與限制，而且只要有解存在，總能以概率1求得方程組的解；不要求組成方程的函數(shù)連續(xù)、可導(dǎo)，可求解常規(guī)算法難以求解的超越方程組.在求解變循環(huán)發(fā)動(dòng)機(jī)數(shù)學(xué)模型時(shí)，多次運(yùn)行GA均能得到結(jié)果，且每次運(yùn)行所得結(jié)果較接近，該方法收斂性與穩(wěn)定性良好.

3.3 求解算法改進(jìn)展望

根據(jù)GA和牛頓-拉夫遜兩種算法的優(yōu)缺點(diǎn)，可以將兩種方法相結(jié)合，避免GA計(jì)算效率低和牛頓-拉夫遜容易局部收斂的不足，高效快速地求出最優(yōu)解.算法的具體步驟為：

1）采用GA產(chǎn)生初始種群，經(jīng)過(guò)選擇、交叉、變異等操作，不斷產(chǎn)生新種群，開(kāi)始進(jìn)化時(shí)，種群向高適應(yīng)度方向快速收斂，當(dāng)進(jìn)化到G代數(shù)后，停止進(jìn)化，按照精英選擇法則選擇適應(yīng)度的N個(gè)個(gè)體.

2）將N個(gè)個(gè)體作為初始值，按照牛頓-拉夫遜的迭代法則，求出N個(gè)有效解，最后從中選擇最優(yōu)解.

改進(jìn)方法的有效性說(shuō)明：開(kāi)始采用GA求解N個(gè)個(gè)體，不僅對(duì)初始值的選取沒(méi)有任何要求，而且收斂速度也很快，求出的N個(gè)個(gè)體采用精英選擇法則選取的，適應(yīng)度高，更接近于最優(yōu)解.以這些在最優(yōu)解附近的個(gè)體作為牛頓-拉夫遜算法的初始值，避免了由于初始值設(shè)計(jì)不合理而陷入局部最優(yōu)解，所以很容易得到全局最優(yōu)解，同時(shí)收斂速度很快.

4 結(jié) 論

本文采用GA求解變循環(huán)發(fā)動(dòng)機(jī)的數(shù)學(xué)模型，將方程組求解問(wèn)題轉(zhuǎn)化為函數(shù)優(yōu)化問(wèn)題，設(shè)計(jì)了算法的求解過(guò)程，獲取了算法的運(yùn)行結(jié)果，提出了算法的有效性評(píng)價(jià)指標(biāo)：初值敏感性、計(jì)算效率、收斂性與穩(wěn)定性.經(jīng)過(guò)上述分析可知GA在計(jì)算效率上與牛頓-拉夫遜相比有很大差距，但是初值敏感性較低，收斂性較好，GA的優(yōu)勢(shì)比較明顯.求解算法改進(jìn)展望較好地結(jié)合了遺傳算法和牛頓-拉夫遜的優(yōu)點(diǎn)，為后期的研究提出了思路.

［1］Willis E A，Welliver A D.Variable Cycle Engines for Super－Sonic Cruising Aircraft［R］.Palo Alto，California：AIAA／SAE 12th Joint Propulsion Conference，1981.

［2］茍學(xué)中，周文祥，黃金泉.變循環(huán)發(fā)動(dòng)機(jī)部件級(jí)建模技術(shù)［J］.航空動(dòng)力學(xué)報(bào)，2013，28（1）：104－111.

［3］戴耀松.固體火箭-沖壓發(fā)動(dòng)機(jī)的研究進(jìn)展［J］.推進(jìn)技術(shù)，1987，8（5）.

［4］Broyden C G.Quasi－Newton Methods and Their Application to Function Minimization［J］.Mathematics of Computation，1967，21：368－381.

［5］蘇三買(mǎi).遺傳算法及其在航空發(fā)動(dòng)機(jī)非線性數(shù)學(xué)模型中的應(yīng)用研究［D］.西安：西北工業(yè)大學(xué)，2002：61－65.

［6］周 明，孫樹(shù)棟.遺傳算法原理及應(yīng)用［M］.北京：國(guó)防工業(yè)出版社，1999.

［7］蘇三買(mǎi).遺傳算法在航空發(fā)動(dòng)機(jī)非線性數(shù)學(xué)模型中的應(yīng)用［J］.推進(jìn)技術(shù)，2004，25（3）：237－240.

［8］茍學(xué)中.變循環(huán)發(fā)動(dòng)機(jī)建模及控制規(guī)律研究［D］.南京：南京航空航天大學(xué)，2011：23－27.

［9］張志涌.精通 Matlab6.5［M］.北京：北京航空航天大學(xué)出版社，2003.

［10］陳秋蓮.基于遺傳算法工具箱的優(yōu)化計(jì)算實(shí)現(xiàn)［J］.現(xiàn)代電子技術(shù)，2007（2）：124－129.