粒子群混合算法在變導向器面積計算中的應用

王 軍，趙 肅

（中航工業(yè)沈陽發(fā)動機設計研究所，沈陽110015）

0 引言

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

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

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

1 粒子群算法及其改進

1.1 基本的粒子群算法

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

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

1.2 改進的粒子群算法

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

帶“被動聚集壓力”因子的PSO

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

自適應慣性權重的PSO

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

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

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

加速收斂的PSO算法（ACPSO）

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

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

1.3 粒子群混合算法

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

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

2 應用粒子群混合算法求解發(fā)動機部件模型

2.1 發(fā)動機部件模型的建立

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

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

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

2.2.1 學習因子

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

圖1 學習因子對收斂性影響對比

2.2.2 慣性權重

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

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

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

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

式中：Fitk為某個粒子第k 次迭代時的適應度值。

對于所研究的部件模型，2種變慣性權重方法的收斂性對比如圖2所示。從圖中可見，以適應度值為自變量的慣性權重的收斂速度更快、收斂精度更高。

圖2 慣性權重對收斂性影響對比

2.2.3 局部改進的PSO

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

圖3 局部改進方法對收斂性影響對比

2.2.4 迭代誤差限

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

3 變導向器面積計算實例與試驗驗證

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

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

3.1 高壓渦輪導向器面積變化

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

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

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

圖5 高壓渦輪導向器面積變化對排氣溫度的影響

圖6 高壓渦輪導向器面積變化對推力的影響

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

3.2 低壓渦輪導向器面積變化

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

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

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

圖8 低壓渦輪導向器面積變化對排氣溫度的影響

圖9 低壓渦輪導向器面積變化對推力的影響

4 總結(jié)

為滿足快速收斂，提高收斂性和收斂效率的要求，設置適合的學習因子、慣性權重、迭代誤差限及選擇合適的粒子群改進算法，可改善粒子群算法的前段收斂速度和后段的收斂精度。

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

[1]童凱生.航空渦輪發(fā)動機性能變比熱計算方法[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ā)動機動態(tài)數(shù)學模型非線性方程組解法研究[J].長沙航空職業(yè)技術學院學報，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]李家瑞，孫健國，張紹基.航空發(fā)動機總體性能數(shù)學模型的1種收斂算法[J].航空發(fā)動機，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),王如根 等.航空燃氣渦輪發(fā)動機數(shù)值仿真[M].北京：國防工業(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]陳長憶，葉永春.基于粒子群算法的非線性方程組求解[J].計算機應用與軟件，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ā)動機數(shù)學模型解法[J].推進技術，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]劉志雄，梁華.粒子群算法中隨機數(shù)參數(shù)的設置與試驗分析[J].控制理論與應用，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]張強，李本威，馬力.粒子群算法在航空發(fā)動機部件模型求解中的應用[J].系統(tǒng)仿真學報，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].控制理論與應用,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]任子暉，王堅.加速收斂的粒子群優(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ā)動機數(shù)學模型解法[J].燃氣渦輪試驗與研究，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ā)動機原理[M].西安：西北工業(yè)大學出版社，2005：252-255.LIAN Xiaochun，WU Hu.Theory of aeroengine[M].Xi’an：Northwest Polytechnical University Press，2005：252-255.(in Chinese)