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

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

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

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

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

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

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

圖1 變循環(huán)發(fā)動機工作原理圖

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

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

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

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

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

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

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

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

表1 獨立變量及其約束

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

約束條件

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

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

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

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

Step1：初始種群的生成.

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

Step2：適應(yīng)度值評價.

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

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

Step3：選擇操作.

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

Step4：交叉操作.

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

Step5：變異操作.

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

Step6：終止條件判斷.

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

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

圖2 遺傳算法流程圖

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

3.1 仿真結(jié)果

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

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

3.2 算法有效性評價

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

［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ā)動機部件級建模技術(shù)［J］.航空動力學(xué)報，2013，28（1）：104－111.

［3］戴耀松.固體火箭-沖壓發(fā)動機的研究進(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］蘇三買.遺傳算法及其在航空發(fā)動機非線性數(shù)學(xué)模型中的應(yīng)用研究［D］.西安：西北工業(yè)大學(xué)，2002：61－65.

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

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

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

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

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