吳正佳 孟榮華 余 剛 何海洋 張 屹
(三峽大學(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ù).
變循環(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為歸一化殘差.
變循環(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的有效性.
本文以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é)模型的解.
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 遺傳算法流程圖
將表1的約束條件代入GA的Matlab[9-10]程序,運行GA程序求解變循環(huán)發(fā)動機數(shù)學(xué)模型,運行時間為3598.22s,得到相對誤差0.4454時目標(biāo)函數(shù)取得最小值,求得各自變量的值見表2.
表2 遺傳算法求解多維非線性隱式方程組結(jié)果
本文選用了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)定性良好.
根據(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)解,同時收斂速度很快.
本文采用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.
三峽大學(xué)學(xué)報(自然科學(xué)版)2014年4期