混沌區(qū)間多目標粒子群優(yōu)化算法及其應用

李俊,李濟順*,,HAL Gurgenci,李倫,GUAN Zhiqiang,楊芳

(1.河南科技大學 機電工程學院,河南洛陽 471003;2.河南科技大學 河南省機械設計及傳動系統(tǒng)重點實驗室,河南洛陽 471003;3.University of Queensland,School of Mechanical and Mining Engineering，Brisbane QLD4072,Australia)

超臨界二氧化碳(sCO2)循環(huán)介質與傳統(tǒng)的蒸汽介質相比具有高溫下腐蝕性小,所需壓縮功少,熱效率高等優(yōu)點。因此sCO2布雷頓循環(huán)被認為是下一代應用于包括核能、化石燃料和集中太陽能發(fā)電的技術[1-3]。為了解決高溫下發(fā)電機轉子系統(tǒng)散熱問題,引入了sCO2氣體冷卻裝置。冷卻裝置同時具有輔助支承轉子葉輪懸臂和散熱的作用。為了避免冷卻裝置的剛度為轉子的動力學特性帶來負面影響,同時提升轉子系統(tǒng)動態(tài)性能,本文基于轉子動力學對軸系的相關參數(shù)進行全面的優(yōu)化設計。

軸系的優(yōu)化模型是工程中典型的多目標優(yōu)化問題,各個目標函數(shù)可能是相互沖突的,通常只能得到全局非劣最優(yōu)解或次優(yōu)解,這為傳統(tǒng)數(shù)值優(yōu)化算法的應用帶來了挑戰(zhàn)。為了解決多目標優(yōu)化問題,包括多尺度量子諧振算法[4]、多目標遺傳法[5]、仿生粒子群算法[6]等智能優(yōu)化算法相繼被用在實際工程上。李超等[7]基于試驗應用多目標遺傳算法對某航空發(fā)動機軸系的初始結構布局進行了優(yōu)化設計,并論證了優(yōu)化的結構布局可以提升轉子的動力學特性。張微等[8]綜合轉子動力學特性與阻尼器相關參數(shù),應用胞映射方法對擠壓油膜阻尼器的彈簧剛度和油膜間隙進行多目標優(yōu)化設計,并得到了Pareto最優(yōu)解集。為了充分考慮風力發(fā)電機葉片氣動特性與其結構性能之間的耦合效應,Wang等[9]提出了一種多目標形狀拓撲耦合優(yōu)化方法,同時優(yōu)化了轉子葉片的外形和內部結構布局。Ramanujam等[10]以直升機旋翼為基線構型,通過多變量粒子群優(yōu)化技術在轉子動力學分析中的應用,以實現(xiàn)變速變幾何轉子系統(tǒng)輸出功率最小化。Xu等[11]將RBF神經(jīng)網(wǎng)絡、徑向基函數(shù)和進化算法相結合,提出了一種考慮支承剛度不確定性的轉子系統(tǒng)動態(tài)響應優(yōu)化設計方法。

Lu等[12]利用粒子群優(yōu)化算法和支持向量機,在任意載荷下確定空間鋼構架的測量響應和估計響應之間的映射關系。Mijbas等[13]采用混沌粒子群(CPSO)算法優(yōu)化后的阻尼控制器對電力系統(tǒng)穩(wěn)定器(PSS)和PID進行增強,從而抑制單機無窮大系統(tǒng)的低頻振蕩,取得了較好的控制效果。Norsahperi等[14]研究了基于分數(shù)階神經(jīng)基(NNFOPID)和粒子群(PSOFOPID)的雙旋翼氣動系統(tǒng)的變槳距控制方法,并實驗驗證了兩種控制器的魯棒性和低功耗特性。Nowak等[15]采用響應面-遺傳(RSM-GA)算法結合有限元技術對汽輪機的啟動曲線進行了優(yōu)化,使汽輪機整個啟動過程所用的時間不到優(yōu)化前的一半。Bernardini等[16]將二元遺傳優(yōu)化算法用于識別直升機前飛旋翼的葉片形狀和結構特性等設計參數(shù)。

以上的研究成果表明,群體智能技術在工程優(yōu)化、參數(shù)識別和控制等領域得到廣泛的應用。對于高速軸系優(yōu)化工作所面臨的主要問題包括兩個方面。首先,轉子的工作轉速可能包含或跨越多個臨界轉速,并且要求不平衡響應限制在安全和可以被接受的水平。其次,冷卻裝置的剛度是非線性的而且隨著轉子的轉速改變而變化,因此冷卻裝置的剛度是一個非線性的二維矩陣。針對這兩個問題本文提出了混沌區(qū)間多目標粒子群優(yōu)化(CIMPSO)算法,并對工作轉速為5.0×104r/min,功率為300 kW超臨界二氧化碳微型透平機的軸系相關參數(shù)進行了優(yōu)化設計。

1 微型透平機優(yōu)化問題

1.1 微型透平機主軸系統(tǒng)基本結構

在圖1所示的透平機轉子系統(tǒng)結構中,S1、S2和S3分別是葉輪、冷卻裝置和密封機構軸向尺寸,S1=2D,S2=S3=D;冷卻裝置的介質采用超臨界二氧化碳氣體。

圖1 微型透平機轉子示意圖

相距為L的主要支承的復合剛度可以采用串聯(lián)剛度阻尼系統(tǒng)的表達式[17]來計算，即

(1)

對于采用柔性軸承座的可傾瓦軸承結構,Koil和Coil分別為軸承的剛度矩陣和阻尼矩陣,關于Koil和Coil的計算方法在參考文獻[18]中有詳細的介紹,這里不再贅述。Khub和Chub分別為軸承座的剛度和阻尼矩陣。轉子的渦動頻率表示為復數(shù),即

A=iω

(2)

若轉子的直徑為D,定義軸承間距為

L=bD

(3)

式中b為待優(yōu)化系數(shù)。作為輔助支承其剛度是轉子轉速和轉子位移的非線性函數(shù),可以表示為

(4)

為了建立有限元模型和傳遞矩陣模型,葉輪被等效為直徑d=100 mm,厚度h=10 mm的剛性圓盤,其質量為0.71 kg。此時圓盤和葉輪具有相同的直徑轉動慣量和極轉動慣量。若葉輪的平衡等級按照汽輪機轉子標準取G2.5級[17],轉子的工作轉速(即工頻)為nwork,則其殘余質量最大偏心距e=75/(πnwork)。

1.2 轉子系統(tǒng)優(yōu)化目標函數(shù)

(5)

式中:輸出參數(shù)ei和fi為分別表示復振型和截面的內力(彎矩和剪力);u為盤軸單元對應的傳遞矩陣,關于傳遞矩陣的詳細介紹可以參考文獻[17]。轉子系統(tǒng)設計是一個多目標優(yōu)化問題,可以作為適應度評價函數(shù)有

1) 殘余質量的不平衡響應最大的峰值

(6)

2) 距離工作轉速最近的臨界轉速與工作轉速(工頻)之間的差值

(7)

3) 工頻點處不平衡響應幅值變化率

(8)

在實際數(shù)值計算中可以寫成差分形式

(9)

式中:An+1和An-1分別表示與工頻點相鄰的轉速點處的不平衡響應的幅值,同樣可以由復振型ei求得。Δn為掃頻間隔。則其約束優(yōu)化數(shù)學模型可以寫為

f(x)=f(A(x),Δn(x),ADC(x))

(10)

由于多個目標函數(shù)可能相互沖突,而且目標函數(shù)值并不在同一數(shù)量級,因此種群的適應度值的計算可以采用各函數(shù)的曼哈頓距離加權求和(MMD)的方式。上述A(x)的作用是抑制轉子的共振幅值,而Δn和ADC(x)的作用是使工頻遠離共振頻率且確保轉子運行在穩(wěn)態(tài)渦動區(qū)。

1.3 約束空間

為輸出部件預留一定的空間的同時保證轉子的剛度,因此參數(shù)b的約束范圍是

3≤b≤6

(11)

軸承座的相關參數(shù)約束范圍取

(12)

(13)

式中:Kupper為工程實際應用中冷卻裝置剛度的上限,由潤滑介質壓強,冷卻裝置的結構參數(shù)所決定。對于軸系來說,每一次掃頻分析,輸入系統(tǒng)的冷卻裝置剛度都是一個與轉速和轉子位移有關二維矩陣。由于氣膜剛度的非線性,在數(shù)值分析的過程中冷卻裝置的剛度是未知的。這為常規(guī)的遺傳算法、模擬退火算法和粒子群算法等在軸系中優(yōu)化的應用帶來了挑戰(zhàn)。針對該問題,本文提出基于粒子群的混沌區(qū)間多目標優(yōu)化算法(CIMPSO),并對一種超臨界二氧化碳太陽能微型透平機的轉子系統(tǒng)進行了優(yōu)化設計。在上述冷卻裝置剛度約束范圍內取特征數(shù)Kge作為粒子向量中的一個元素,則待優(yōu)化的粒子向量可被描述為

(14)

2 轉子系統(tǒng)優(yōu)化方法

2.1 標準粒子群優(yōu)化算法

基本粒子群優(yōu)化(Particle swarm optimization,PSO)計算方法是一種基于自我認知和社會心理學的計算技術。在尋優(yōu)空間中的每一個潛在可行解都被抽象為個體粒子(Particle)向量xi。粒子群算法首先在約束空間中隨機初始種群,每個粒子由速度和位移進行定義。在搜索空間中個體粒子在局部尋優(yōu)的同時實現(xiàn)信息共享機制。粒子正確的行為得到加強,同時被其他粒子所模仿,最終收斂到最優(yōu)解或者非劣最優(yōu)目標域(Pareto front)。

標準粒子群優(yōu)化算法迭代更新粒子[6]的表達式為:

vi(t+1)=wvvi(t)+c1r1[pi(t)-xi(t)]+

c2r2[pg(t)-xi(t)]

(15)

xi(t+1)=wxxi(t)+vi(t+1)

(16)

設初始化粒子的個數(shù)為N,i=1,2…,N,表示粒子的編號。每個粒子由目標函數(shù)f(x)決定其適應度。在迭代過程中各個粒子記憶、追蹤當前個體最優(yōu)粒子pi(t)和全局最優(yōu)粒子pg(t)。wv和wx為慣性權重,進一步?jīng)Q定了算法的局部和全局尋優(yōu)能力。r1、r2取0到1區(qū)間的隨機數(shù),表示粒子的認知能力。c1、c2為學習因子,表示學習能力。

2.2 混沌區(qū)間多目標粒子群模型

非線性粒子指得是每個粒子中的若干元素是一個二維非線性變化的元素。對于本文軸系來說冷卻裝置的剛度是一個根據(jù)轉軸轉速和位移變化的參數(shù)。粒子中非線性元素的取值受到工程應用的約束。設約束的上限和下限分別為：

(17)

式中:α和β是基于粒子元素特征數(shù)的拓展系數(shù)。為了實現(xiàn)含有非線性元素的粒子群尋優(yōu),這里引入混沌變量。由典型的混沌系統(tǒng)Logistic方程[6]，可以得到具有隨機性的混沌變量?；煦缱兞吭谝欢ǚ秶鷥润w現(xiàn)規(guī)律性、遍歷性和隨機性。

zn+1=μzn(1-zn)n=1,2,…,N

(18)

式中：μ為控制參量；zn∈[0,1]。每一迭代步對粒子xi更新后,將冷卻裝置剛度特征數(shù)Kge映射到Logistic方程,進行混沌變量初始化

(19)

式中t表示當前迭代步。將z1(t)代入Logistic方程進行迭代產(chǎn)生混沌變量數(shù)列zi+1(t)(i=1,2,…,N)。把混沌變量數(shù)列zi+1(t)逆映射到求解空間中得到

(20)

(21)

pg(t)=comfit[pi(t)]

(22)

式中comfit(·)指的是當前迭代步綜合目標函數(shù)最優(yōu)的粒子。每個粒子的子目標函數(shù)的歸一化為：

(23)

(24)

通過加權的方式得到每個個體粒子的綜合目標函數(shù)值為

(25)

(26)

(27)

而當前迭代步中被替換后的粒子中剛度特征數(shù)元素任然采用Kge(t)。最后按照標準粒子群算法式(15)和式(16)迭代得到新粒子

(28)

上文所描述的算法稱為混沌區(qū)間多目標粒子群優(yōu)化算法(Chaos interval multi-objective particle swarm optimization algorithm,CIMPSO)。與標準混沌粒子群優(yōu)化算法(Chaos particle swarm optimization,ChPSO)相比,本文提出的CIMPSO可以對多目標和含有二維非線性元素的粒子群模型進行空間尋優(yōu)計算。當個體最優(yōu)粒子和全局最優(yōu)粒子不在隨迭代步進化,此時粒子群的各個粒子獨立更新,最終所有粒子趨于全局最優(yōu)解或者非劣最優(yōu)解。而粒子群中所有粒子都包含有二維非線性元素特征數(shù)的混沌映射,這就保證了在映射區(qū)間范圍內二維非線性元素的區(qū)間得到有效的空間尋優(yōu)計算。

2.3 CIMPSO算法歸納和收斂性分析

CIMPSO的主要優(yōu)勢在于解決含有非線性、二維變化的元素的多目標優(yōu)化問題。算法主要分為以下幾個步驟。

步驟1 初始化相關參數(shù),包括學習因子c1、c2,粒子群規(guī)模N,尋優(yōu)次數(shù)T,速度慣性權重wv和位置遺忘慣性權重wx,拓展系數(shù)α和β,控制參量μ,加權系數(shù)φ1,φ2,…,φN。在約束范圍內隨機初始化粒子群中的所有粒子位置。

步驟2 粒子位置記為xi(t)。確定粒子群非線性元素的特征數(shù)xi(t,j)(j表示粒子向量中非線性元素的位置編號,t表示迭代步數(shù))。將特征數(shù)映射迭代產(chǎn)生混沌變量序列

(29)

(30)

用MMD算法計算當前層級中粒子的每個目標函數(shù)的Manhattan距離Mi。由加權的方式得到每個粒子的綜合適應度

fi(x)=φ1M1+φ2M2+…+φkMk

(31)

則全局最優(yōu)解即為綜合目標函數(shù)的極值所對應的粒子。

pg(t)=comfit[Pi(t)]

(32)

步驟4 生成混沌變量序列

(33)

Zitzler提出的ZDT問題是常用的兩目標測試函數(shù)[19]。本文選取ZDT1和ZDT2作為被測對象:

f1(x)=x1

(34)

f2(x)=g(x)[1-(f1(x)/g(x))m]

(35)

xi∈[0,1],i=1,2,…n;n=30

(36)

當m=1/2時上式為ZDT1(凸函數(shù)),當m=2時為ZDT2(凹函數(shù))。它們的理想Pareto前沿端已經(jīng)被給出,如圖2中的點劃線和虛線所示。元素x1同時影響兩個目標函數(shù),因此被選為特征數(shù)元素。圖中的藍色短線段為數(shù)值實驗所得到的解區(qū)間。對應的初始化粒子個數(shù)為50,迭代步數(shù)為200步,尋優(yōu)計算次數(shù)為200次。從圖中與理想Pareto前沿對比可以明顯看出實驗解以一定的誤差范圍收斂于非劣最優(yōu)解(PFtrue)。其中ZDT1和ZDT2的測試解(EXPERIMENT)的平均誤差分別為e1=16.1%和e2=11.69%。

圖2 理想Pareto前沿(PFtrue)和測試解集(EXPERIMENT)

3 轉子系統(tǒng)模型尋優(yōu)計算

如上文所述超臨界二氧化碳透平機轉子系統(tǒng)待優(yōu)化參數(shù)包括主要支承系統(tǒng)軸承間距L=bD,軸承座質量mhub、剛度Khub和阻尼Chub,冷卻冷卻裝置的剛度Kg。其中Kg是轉子轉速和位移的函數(shù),其特征數(shù)取Kge。各元素約束條件為

(37)

表1 粒子群模型參數(shù)

表1中參數(shù)的取值范圍均在0～1之間,其中拓展系數(shù)、控制參量以及曼哈頓權值均由決策者根據(jù)具體的優(yōu)化問題的側重點自行設定。學習因子c1和c2分別影響算法的局部搜索和全局搜索能力,為了避免尋優(yōu)過程陷入局部極值和提高收斂的精度,c1和c2均取中間偏大值0.6。慣性權重wv和wx是用來進一步平衡局部尋優(yōu)和全局尋優(yōu)之間的矛盾,在算法完成后通過測試調整慣性權重的大小。

按照設計經(jīng)驗,軸承的間距取軸徑的3～5倍,這里取3;軸承座的參振質量取轉子質量(2.49 kg)的一半,即1.2 kg;軸承座的剛度取可傾瓦軸承剛度的二分之一,阻尼取可傾瓦軸承的阻尼系數(shù);冷卻裝置的剛度任取約束空間的中間值。則經(jīng)驗粒子參數(shù)為

(38)

圖3所示的點劃線表示的是在該參數(shù)下轉子系統(tǒng)的不平衡響應。

圖3 撓性軸承座可傾瓦軸承轉子系統(tǒng)優(yōu)化設計前后轉子不平衡響應

由圖3可知優(yōu)化前相對峰值為13.51-2.0=11.51。在約束空間內隨機初始化粒子,連續(xù)尋優(yōu)計算25次,分別將25組非劣最優(yōu)解帶入轉子系統(tǒng)所產(chǎn)生的不平衡響應如圖中實線所示。與點劃線相比,實線表現(xiàn)出的最小峰值衰減率為69.6%,且轉子的響應曲線更加平緩,體現(xiàn)了CIMPSO算法的魯棒性。工作轉速(nwork=5.0×104r/min)與距離最近的臨界轉速的差值約占工作轉速的30%,保證了轉子工作在平滑的穩(wěn)態(tài)渦動區(qū)。

柔性轉子系統(tǒng)的臨界轉速和不平衡響應與軸承的特性密切相關,為了驗證算法可以適用于不同的轉子系統(tǒng),現(xiàn)對主動磁懸浮(AMB)軸系進行優(yōu)化計算。AMB軸系的待優(yōu)化參數(shù)為軸承間距bD和冷卻裝置的剛度kcb,與撓性軸承座可傾瓦軸系相比可以優(yōu)化的元素較少,且在所研究的頻率范圍內,AMB軸系出現(xiàn)了懸臂主導的多個模態(tài)頻率。圖4為 AMB支撐的微型透平機優(yōu)化前后葉輪的頻響曲線圖。磁懸浮軸承的剛度和阻尼分別為kAMB=8.2×105N/m、CAMB=100 Nm/s。圖中實線和虛線分別為多次尋優(yōu)計算得到的Pareto前端和優(yōu)化前的經(jīng)驗粒子xe所確定的諧響應曲線。與虛線相比,實線的最大幅值有了明顯的衰減。另一方面,隨著模態(tài)頻率向低頻方向移動,部分優(yōu)化解使得工頻范圍內懸臂的模態(tài)頻率削減為單個。

圖4 主動磁懸浮軸承轉子系統(tǒng)優(yōu)化設計前后轉子不平衡響應

任取數(shù)值計算所得的非劣優(yōu)化解,即

x=[3.8 1.6 1.01×1085.05×1042.43×105]

(39)

特征數(shù)x(5)所對應的冷卻裝置的剛度變化范圍是Kg=[Kge(1-α),Kge(1+β)]。為了直觀分析算法的收斂過程,令全局適應度為

f=1/A(x)+Δn(x)

(40)

該組優(yōu)化解的收斂過程如圖5所示。由圖可知,全局適應度在算例迭代110步后穩(wěn)定。

圖5 CIMPSO算法的收斂性

圖6中的實線和虛線分別表示目標函數(shù)A(x)/e和Δn(x)隨著迭代計算的收斂過程。由曲線的變化趨勢可知,目標函數(shù)A(x)和Δn(x)是相互沖突的,虛線位置向下波動的同時實線必然相對向上波動。整個優(yōu)化過程伴隨著A(x)和Δn(x)之間的博弈。最終全局收斂時,虛線所代表的目標函數(shù)Δn(x)妥協(xié),適量增加以換取峰值目標函數(shù)A(x)的減小。再一次證明了CIMPSO算法對有沖突的多目標轉子優(yōu)化問題可以進行有效的尋優(yōu)計算。

圖6 子目標函數(shù)變化趨勢

4 仿真驗證

傳遞矩陣法的計算精度取決于轉子被分成盤軸單元的數(shù)量。盤軸單元數(shù)愈多計算精度越高,但是累積的數(shù)值誤差越大。本節(jié)以轉子諧響應為目標建立其有限元模型,對優(yōu)化結果進行對比驗證。如圖7所示,在有限元三維模型中葉輪等效為相同質量、轉動慣量的質量點。柔性軸承座有效參振質量為mhub,在圖中顯示為圓環(huán)形輪轂。輪轂的內側和外側分別是虛擬的可傾瓦油膜軸承和剛度阻尼系統(tǒng)。

圖7 轉子系統(tǒng)三維有限元模型

分別將xopt和xcon帶入有限元模型中,諧響應分析結果如圖8所示,分別對應點劃線和虛線。優(yōu)化后相對峰值(10.04-1.65=8.39)同優(yōu)化前(2.94-1.6=1.34)相比,衰減84.02%。另一方面,目標函數(shù)Δn(x)與優(yōu)化前相比增加了12%,進一步提高了轉子系統(tǒng)設計制造的容錯能力,同時驗證了采用CIMPSO算法得出的優(yōu)化解的有效性。圖中實線為算法中采用傳遞矩陣法得到的轉子不平衡響應。與有限元法相比,最大幅值誤差為20.1%,穩(wěn)態(tài)渦動半徑和臨界轉速基本保持一致,這表明采用CIMPSO算法得到的優(yōu)化解對于轉子系統(tǒng)的設計、控制轉子的振動具有一定的指導意義。

圖8 轉子系統(tǒng)諧響應分析結果

5 結論

針對超臨界二氧化碳微型透平機的轉子系統(tǒng)進行了多目標優(yōu)化設計?；诨煦鐓^(qū)間粒子群優(yōu)化算法,提出了可以考慮設計變量為二維非線性元素的優(yōu)化設計模型,并歸納總結了混沌區(qū)間粒子群優(yōu)化算法(CIMPSO)。使用該算法對轉子系統(tǒng)進行了優(yōu)化設計,得到了全局非劣最優(yōu)或者次優(yōu)解集,驗證了算法的穩(wěn)定性。最后使用商業(yè)軟件建立有限元模型進行仿真驗證,結果表明,經(jīng)過CIMPSO算法優(yōu)化后的參數(shù)可以顯著改善轉子的動力學性能。