王 軍，趙 肅

（中航工業(yè)沈陽(yáng)發(fā)動(dòng)機(jī)設(shè)計(jì)研究所，沈陽(yáng)110015）

0 引言

目前，航空發(fā)動(dòng)機(jī)高精度穩(wěn)態(tài)數(shù)學(xué)模型均按變比熱計(jì)算方法[1]建立，具有高度非線性、基于部件特性的特點(diǎn)。模型一般采用迭代方法求解發(fā)動(dòng)機(jī)共同工作方程，常用的迭代方法有Newton-Raphson（N-R）法、Broyden秩1法[2]、N+1點(diǎn)殘量法[3]和最速下降法，上述方法具有嚴(yán)格的數(shù)理局部收斂性，且收斂性依賴初值的選取，雖通過(guò)采用阻尼因子、松弛因子、初值的有限域優(yōu)化探索和迭代步長(zhǎng)的線性探索與回溯[4]等改進(jìn)方法可擴(kuò)大其收斂范圍，但一般改善程度都較?。昏b于發(fā)動(dòng)機(jī)工作包線寬廣及在特殊工作條件（如幾何面積突變等）下，發(fā)動(dòng)機(jī)各部件共同工作時(shí)，因初值偏離最優(yōu)解較大導(dǎo)致模型出現(xiàn)不收斂的情況。

為了克服上述算法在收斂性方面的不足，粒子群算法（ParticleSwarmOptimization，PSO）等進(jìn)化算法被引入發(fā)動(dòng)機(jī)部件模型等非線性方程的求解[5-6]中，取得了較好效果。PSO具有全局收斂的能力，在進(jìn)化初期收斂速度快，運(yùn)算簡(jiǎn)單，易于實(shí)現(xiàn)，但其計(jì)算效率和收斂精度偏低，同時(shí)存在最優(yōu)解不穩(wěn)定的問(wèn)題。結(jié)合常規(guī)迭代算法計(jì)算效率高和PSO全局收斂的優(yōu)點(diǎn)的混合算法改善了發(fā)動(dòng)機(jī)部件模型求解過(guò)程中的收斂性。

本文采用混合粒子群算法解決在改變發(fā)動(dòng)機(jī)渦輪導(dǎo)向器面積對(duì)性能影響計(jì)算中迭代不收斂的問(wèn)題，以滿足穩(wěn)態(tài)性能仿真的需要。

1 粒子群算法及其改進(jìn)

1.1 基本的粒子群算法

基本的粒子群算法[7]以模擬鳥(niǎo)的群體智能為特征，以求解優(yōu)化問(wèn)題為背景。每只鳥(niǎo)被稱為1個(gè)粒子，每個(gè)粒子用其幾何位置和速度向量表示，參考各自的既定方向、個(gè)體所經(jīng)歷的最優(yōu)方向和整個(gè)鳥(niǎo)群所經(jīng)歷的最優(yōu)方向來(lái)確定飛行。假設(shè)在1個(gè)D 維的目標(biāo)探索空間中有n個(gè)粒子，其中第i個(gè)粒子表示為1個(gè)D 維的向量，xi=（xi1，xi2，…，xiD）（i=1，2，…，n），表示第i個(gè)粒子在此探索空間中的位置，vi=（vi1，vi2，…，viD），表示第i個(gè)粒子的速度。假設(shè)第i個(gè)粒子迄今為止探索到的最優(yōu)位置為pi=（pi1，pi2，…，piD），整個(gè)粒子群迄今為止探索到的最優(yōu)位置為pg=（pg1，pg2，…，pgD），采用下列公式對(duì)粒子群進(jìn)行速度和位置更新

式中：i=1，2，…，n，d=1，2，…，D；c1，c2為非負(fù)常數(shù)的學(xué)習(xí)因子；r1，r2為 [0,1]間的隨機(jī)數(shù)；vid=[vmin，vmax]，vmin/vmax為更新速度的最大/最小邊界；ω 為慣性權(quán)重。

1.2 改進(jìn)的粒子群算法

在基本粒子群算法的基礎(chǔ)上，為提高粒子全局探索能力，發(fā)展了帶“被動(dòng)聚集壓力”因子的PSO、自適應(yīng)慣性權(quán)重的PSO[8]、混合探索粒子群算法（MSPSO）[9]及加速收斂的PSO（ACPSO）[10]等。

帶“被動(dòng)聚集壓力”因子的PSO

式中：c3為被動(dòng)聚集壓力因子；Prd為粒子群中隨機(jī)選擇的1個(gè)粒子；r3為[0,1]間的隨機(jī)數(shù)。

自適應(yīng)慣性權(quán)重的PSO

式中:ωi為慣性權(quán)重，根據(jù)適應(yīng)度函數(shù)值或者迭代次數(shù)自動(dòng)調(diào)整。

混合探索粒子群算法（MSPSO）

式中：α 為[0,1]之間的1個(gè)常數(shù)；k 為迭代次數(shù)；Pld是第k 代種群中粒子最好位置。

加速收斂的PSO算法（ACPSO）

式中:Θ 為三角函數(shù)算子，一般取Θ=sin；α 為角度值，一般取α∈[0，π/8]；β 為大于零的常數(shù)，一般取β=3。

目前對(duì)PSO的改進(jìn)主要集中在算法參數(shù)和粒子更新結(jié)構(gòu)的調(diào)整上，目的是使粒子跳出局部最優(yōu)，使其全局和局部探索能力達(dá)到最佳平衡，提高算法的性能。但從航空發(fā)動(dòng)機(jī)部件模型求解實(shí)例來(lái)看，因部件模型高度非線性化，導(dǎo)致模型的收斂速度和精度均低于傳統(tǒng)N-R等算法的收斂速度和精度，為進(jìn)一步提高模型的收斂效率，在上述研究基礎(chǔ)上發(fā)展了粒子群混合算法。

1.3 粒子群混合算法

為綜合傳統(tǒng)算法和粒子群算法的優(yōu)點(diǎn)，提出了PSO-NR（粒子群-牛頓拉夫遜）和PSO-N+1（粒子群-N+1點(diǎn)殘量法）[11]等混合算法。

以PSO-NR混合算法為例說(shuō)明混合算法的工作原理：在該算法中，N-R法仍為求解發(fā)動(dòng)機(jī)部件非線性模型的主體算法，在性能計(jì)算時(shí)對(duì)在設(shè)定循環(huán)迭代次數(shù)內(nèi)不收斂的工作點(diǎn)采用PSO算法，變量初值采用N-R法最后1次循環(huán)獲得的數(shù)值。當(dāng)循環(huán)迭代次數(shù)達(dá)到設(shè)定值后（根據(jù)計(jì)算精度設(shè)置合適的迭代次數(shù)），再次采用N-R法進(jìn)行迭代計(jì)算，達(dá)到后期快速收斂的目的。如果計(jì)算仍不收斂，考慮到PSO獲得的最優(yōu)解不穩(wěn)定，再次采用混合算法，使用次數(shù)一般不大于10次，以免陷入死循環(huán)。

2 應(yīng)用粒子群混合算法求解發(fā)動(dòng)機(jī)部件模型

2.1 發(fā)動(dòng)機(jī)部件模型的建立

以某型軍用混合排氣渦扇發(fā)動(dòng)機(jī)為研究對(duì)象，按照輸入的控制規(guī)律和變量初值及部件間遵循的流量、壓力和功率平衡原則建立發(fā)動(dòng)機(jī)共同工作方程[12]，并將其轉(zhuǎn)換為誤差方程（殘差方程）

PSO是1種優(yōu)化算法，采用PSO求解非線性方程組，需要將方程組的求解問(wèn)題轉(zhuǎn)化為函數(shù)的優(yōu)化問(wèn)題。應(yīng)用無(wú)約束優(yōu)化方法求解非線性方程組（式（6））時(shí)，通常將其轉(zhuǎn)化為非線性最小二乘問(wèn)題：

2.2 粒子群混合算法參數(shù)的設(shè)置

2.2.1 學(xué)習(xí)因子

對(duì)于學(xué)習(xí)因子c1和c2，關(guān)系到個(gè)體最優(yōu)與全局最優(yōu)對(duì)粒子的影響程度。數(shù)學(xué)研究顯示，c1+c2>4且c1>c2時(shí)收斂性較好。學(xué)習(xí)因子對(duì)收斂性影響對(duì)比如圖1所示，采用文獻(xiàn)[10]中3組取值求解部件模型的收斂情況。

圖1 學(xué)習(xí)因子對(duì)收斂性影響對(duì)比

2.2.2 慣性權(quán)重

隨著迭代次數(shù)的增加，最優(yōu)解的探索范圍將逐漸縮小，對(duì)于慣性權(quán)重，變慣性權(quán)重的收斂效果要比定慣性權(quán)重的好。1種方法采用遞減函數(shù)[10]來(lái)保證算法不會(huì)因粒子運(yùn)動(dòng)慣性過(guò)大而造成收斂緩慢，另1種方法是慣性權(quán)重隨著粒子適應(yīng)度的變化而變化，適應(yīng)度值增大慣性權(quán)重也增大，反之隨其減小而減小。

前者ω 函數(shù)定義為

式中：ωmax、ωmin分別為慣性權(quán)重的上、下限；T 為迭代總次數(shù)；n 為當(dāng)前迭代次數(shù)；x 為函數(shù)的凸凹形態(tài)。

后者ω 函數(shù)定義為

式中：Fitk為某個(gè)粒子第k 次迭代時(shí)的適應(yīng)度值。

對(duì)于所研究的部件模型，2種變慣性權(quán)重方法的收斂性對(duì)比如圖2所示。從圖中可見(jiàn)，以適應(yīng)度值為自變量的慣性權(quán)重的收斂速度更快、收斂精度更高。

圖2 慣性權(quán)重對(duì)收斂性影響對(duì)比

2.2.3 局部改進(jìn)的PSO

利用第2.1節(jié)建立的發(fā)動(dòng)機(jī)部件模型，測(cè)試基本PSO、帶“被動(dòng)聚集壓力”因子的PSO、MSPSO和ACPSO等粒子群算法的收斂性，結(jié)果如圖3所示。從圖中可見(jiàn)，3 種改進(jìn)的PSO收斂速度較基本PSO的快，收斂精度差別不大，本文采用混合探索粒子群算法（MSPSO）。

圖3 局部改進(jìn)方法對(duì)收斂性影響對(duì)比

2.2.4 迭代誤差限

對(duì)于粒子群混合算法，需設(shè)置2個(gè)迭代誤差限，即PSO和N-R 算法的迭代誤差限。在一般情況下，PSO的誤差限要大于N-R 算法的，主要因?yàn)镻SO后段收斂緩慢，較小的誤差限會(huì)導(dǎo)致迭代次數(shù)增加，計(jì)算效率下降，PSO的誤差限可取收斂曲線的拐點(diǎn)。N-R 算法的計(jì)算精度較高，可以取目標(biāo)誤差限作為其誤差限。本文PSO的誤差限取0.03，N-R 算法的誤差限取0.003。

3 變導(dǎo)向器面積計(jì)算實(shí)例與試驗(yàn)驗(yàn)證

一般來(lái)說(shuō)，實(shí)際發(fā)動(dòng)機(jī)很難完全滿足設(shè)計(jì)要求，這就需要發(fā)動(dòng)機(jī)在調(diào)試階段為滿足性能匹配和優(yōu)化的要求，具備一定的調(diào)整能力。主要體現(xiàn)在風(fēng)扇、壓氣機(jī)可調(diào)角度的優(yōu)化，以及渦輪導(dǎo)向器排氣面積、噴管喉道面積的微調(diào)上。壓氣機(jī)可調(diào)葉片角度及噴管喉道面積一般為發(fā)動(dòng)機(jī)控制參數(shù)，其調(diào)節(jié)規(guī)律易于實(shí)現(xiàn)。渦輪導(dǎo)向器在固定涵道比的發(fā)動(dòng)機(jī)上是不可調(diào)的，為了滿足性能優(yōu)化需求，需要生產(chǎn)不同組別的渦輪導(dǎo)向器供試驗(yàn)用，費(fèi)用高、周期長(zhǎng)。利用數(shù)值仿真可在發(fā)動(dòng)機(jī)生產(chǎn)之前確定生產(chǎn)組別的數(shù)量和大小，指導(dǎo)調(diào)試方向，減少試制和試驗(yàn)費(fèi)用。

現(xiàn)有的發(fā)動(dòng)機(jī)穩(wěn)態(tài)數(shù)學(xué)模型大多采用經(jīng)典的N-R 算法，在進(jìn)行渦輪導(dǎo)向器面積變化對(duì)發(fā)動(dòng)機(jī)性能影響計(jì)算時(shí)，因特性圖或折合流量差別較大，導(dǎo)致誤差突變，出現(xiàn)迭代不收斂的現(xiàn)象，采用PSO-NR 算法能夠很好地解決。

3.1 高壓渦輪導(dǎo)向器面積變化

在調(diào)整某型發(fā)動(dòng)機(jī)高壓渦輪導(dǎo)向器面積的計(jì)算時(shí)，可采用小偏差流量不變或更換部件特性的方法。前者主要考慮在慢車轉(zhuǎn)速以上，高壓渦輪導(dǎo)向器處于臨界狀態(tài)，導(dǎo)向器面積的變化可以認(rèn)為只是流過(guò)渦輪的折合流量的變化，且忽略渦輪效率變化的影響，在一定條件下能夠滿足精度需要；后者的計(jì)算精度較高，但需要部件提供精確的計(jì)算或試驗(yàn)修正特性。

在上述計(jì)算條件下，采用PSO-NR 算法計(jì)算高壓渦輪導(dǎo)向器相對(duì)于基準(zhǔn)值偏小3.5%對(duì)轉(zhuǎn)差、低壓渦輪出口排氣溫度T6和推力F 的影響，其計(jì)算和試驗(yàn)結(jié)果的對(duì)比如圖4～6所示。

圖4 高壓渦輪導(dǎo)向器面積變化對(duì)轉(zhuǎn)差的影響

圖5 高壓渦輪導(dǎo)向器面積變化對(duì)排氣溫度的影響

圖6 高壓渦輪導(dǎo)向器面積變化對(duì)推力的影響

從圖中可見(jiàn)，當(dāng)高壓渦輪導(dǎo)向器面積減小3.5%時(shí)，發(fā)動(dòng)機(jī)轉(zhuǎn)差增大0.6%～1.0%，排氣溫度降低10～15K，推力減小1.2%～1.9%。

3.2 低壓渦輪導(dǎo)向器面積變化

在進(jìn)行調(diào)整某型發(fā)動(dòng)機(jī)低壓渦輪導(dǎo)向器面積的計(jì)算時(shí)，考慮到低壓渦輪導(dǎo)向器僅在高轉(zhuǎn)速范圍內(nèi)才能處于臨界狀態(tài)，采用更換部件特性的方法計(jì)算全轉(zhuǎn)速特性。計(jì)算低壓渦輪導(dǎo)向器相對(duì)基準(zhǔn)值增大4%對(duì)轉(zhuǎn)差、低壓渦輪出口排氣溫度T6和推力F 的影響，其計(jì)算和試驗(yàn)結(jié)果的對(duì)比如圖7～9所示。

從圖中可見(jiàn)，當(dāng)?shù)蛪簻u輪導(dǎo)向器面積增加4%時(shí)，發(fā)動(dòng)機(jī)轉(zhuǎn)差增大0.5%～0.7%，排氣溫度在低轉(zhuǎn)速下略有降低，在高轉(zhuǎn)速下基本一致，推力在全轉(zhuǎn)速范圍內(nèi)一致。

圖7 低壓渦輪導(dǎo)向器面積變化對(duì)轉(zhuǎn)差的影響

圖8 低壓渦輪導(dǎo)向器面積變化對(duì)排氣溫度的影響

圖9 低壓渦輪導(dǎo)向器面積變化對(duì)推力的影響

4 總結(jié)

為滿足快速收斂，提高收斂性和收斂效率的要求，設(shè)置適合的學(xué)習(xí)因子、慣性權(quán)重、迭代誤差限及選擇合適的粒子群改進(jìn)算法，可改善粒子群算法的前段收斂速度和后段的收斂精度。

采用PSO-NR求解本文建立的發(fā)動(dòng)機(jī)共同工作方程組，可有效解決高、低壓渦輪導(dǎo)向器面積改變對(duì)性能影響計(jì)算不收斂的問(wèn)題；采用換特性計(jì)算方法的計(jì)算結(jié)果與試驗(yàn)結(jié)果基本一致。改變渦輪導(dǎo)向器面積對(duì)發(fā)動(dòng)機(jī)性能影響的仿真計(jì)算可為發(fā)動(dòng)機(jī)調(diào)試提供技術(shù)支持。

[1]童凱生.航空渦輪發(fā)動(dòng)機(jī)性能變比熱計(jì)算方法[M].北京：航空工業(yè)出版社，1991：9-18.TONG Kaisheng.Calculation method of performance with variable specific heat for aviation turbine engine[M].Beijing：Aviation Industry Press，1991：9-18.(in Chinese)

[2]熊純.航空發(fā)動(dòng)機(jī)動(dòng)態(tài)數(shù)學(xué)模型非線性方程組解法研究[J].長(zhǎng)沙航空職業(yè)技術(shù)學(xué)院學(xué)報(bào)，2002，2（3）：39-42.XIONG Chun.Study on the solution of nonlinear equation group of dynamic mathematical model of aeroengine[J].Changsha Aeronautical Vocational and Technical Collegle Journal,2002，2（3）：39-42.(in Chinese)

[3]李家瑞，孫健國(guó)，張紹基.航空發(fā)動(dòng)機(jī)總體性能數(shù)學(xué)模型的1種收斂算法[J].航空發(fā)動(dòng)機(jī)，2005，31（4）：48-50.LI Jiarui,SUN Jianguo,ZHANG Shaoji.A convergence algorithm of aeroengine performance mathematical model[J].Aeroengine，2005，31（4）：48-50.(in Chinese)

[4]駱廣琦,桑增產(chǎn),王如根 等.航空燃?xì)鉁u輪發(fā)動(dòng)機(jī)數(shù)值仿真[M].北京：國(guó)防工業(yè)出版社，2007：6-8.LUO Guangqi，SANG Zengchan，WANG Rugen，et al.Numerical methods for aviation gas turbine engine simulation[M].Beijing：National Defense Industry Press，2007：6-8.(in Chinese)

[5]陳長(zhǎng)憶，葉永春.基于粒子群算法的非線性方程組求解[J].計(jì)算機(jī)應(yīng)用與軟件，2006，23（5）：137-139.CHEN Changyi，YE Yongchun.Solving nonlinear systems of equations based on particle swarm optimization[J].Computer Application and Software,2006，23（5）：137-139.(in Chinese)

[6]蘇三買，陳永琴.基于混合遺傳算法的航空發(fā)動(dòng)機(jī)數(shù)學(xué)模型解法[J].推進(jìn)技術(shù)，2007，28（6）：661-664.SU Sanmai，CHEN Yongqin.Hybrid genetic algoritham in solving aeroengine nonlinear mathematical model[J].Journal of Propulsion Technology，2007，28（6）：661-664.(in Chinese)

[7]劉志雄，梁華.粒子群算法中隨機(jī)數(shù)參數(shù)的設(shè)置與試驗(yàn)分析[J].控制理論與應(yīng)用，2010，27（11）：1489-1496.LIU Zhixiong，LIANG Hua.Parameter setting and experimental analysis of the random number in particle swarm optimization algorithm [J].Control Theory&Applications，2010，27（11）：1489-1496.(in Chinese)

[8]張強(qiáng)，李本威，馬力.粒子群算法在航空發(fā)動(dòng)機(jī)部件模型求解中的應(yīng)用[J].系統(tǒng)仿真學(xué)報(bào)，2009，21（12）：3584-3587.ZHANG Qiang,LI Benwei,MA Li.Application of particle swarm optimization in solution to aeroengine component-level model[J].Journal of System Simulation (S1004-731X)，2009，21（12）：3584-3587.(in Chinese)

[9]連志剛，焦斌.一種混合探索的粒子群算法[J].控制理論與應(yīng)用,2010，27（10）：1404-1410.LIAN Zhigang，JIAO Bin.A particle swarm algorithm based on hybrid exploration[J].Control Theory and Application，2010，27(10)：1404-1410.(in Chinese)

[10]任子暉，王堅(jiān).加速收斂的粒子群優(yōu)化算法[J].控制與決策，2011，26（2）：201-206.REN Zihui,WANG Jian.Accelerate convergence particle swarm optimization[J].Control and Decision Making，2011，26（2）：201-206.(in Chinese)

[11]駱廣琦，劉波，宋頔源.基于混合粒子算法的航空發(fā)動(dòng)機(jī)數(shù)學(xué)模型解法[J].燃?xì)鉁u輪試驗(yàn)與研究，2011，24（2）：5-8.LUO Guangqi，LIU Bo，SONG Diyuan.Hybrid particle swarm optimization in solving aeroengine nonlinear mathematical model[J].Gas Turbine Experiment and Research,2011，24（2）：5-8.(in Chinese)

[12]廉筱純，吳虎.航空發(fā)動(dòng)機(jī)原理[M].西安：西北工業(yè)大學(xué)出版社，2005：252-255.LIAN Xiaochun，WU Hu.Theory of aeroengine[M].Xi’an：Northwest Polytechnical University Press，2005：252-255.(in Chinese)