永磁同步風(fēng)力發(fā)電機(jī)隨機(jī)分岔現(xiàn)象的全局分析?

楊黎暉 葛揚(yáng) 馬西奎

(西安交通大學(xué)電氣工程學(xué)院,電力設(shè)備電氣絕緣國(guó)家重點(diǎn)實(shí)驗(yàn)室,西安 710049)

永磁同步風(fēng)力發(fā)電機(jī)隨機(jī)分岔現(xiàn)象的全局分析?

楊黎暉?葛揚(yáng) 馬西奎

(西安交通大學(xué)電氣工程學(xué)院,電力設(shè)備電氣絕緣國(guó)家重點(diǎn)實(shí)驗(yàn)室,西安 710049)

永磁同步風(fēng)力發(fā)電機(jī)在運(yùn)行過程中不可避免地會(huì)受到風(fēng)能的隨機(jī)干擾,本文建立了在輸入機(jī)械轉(zhuǎn)矩存在隨機(jī)干擾情況下永磁同步風(fēng)力發(fā)電機(jī)的數(shù)學(xué)模型,采用胞映射方法分析了隨機(jī)干擾強(qiáng)度變化時(shí)系統(tǒng)全局結(jié)構(gòu)的演化行為,并通過數(shù)值模擬對(duì)理論分析進(jìn)行驗(yàn)證.研究結(jié)果表明,隨著隨機(jī)干擾強(qiáng)度的增大,系統(tǒng)中會(huì)出現(xiàn)隨機(jī)內(nèi)部激變和隨機(jī)邊界激變,即由于隨機(jī)吸引子與其吸引域內(nèi)的隨機(jī)鞍發(fā)生碰撞而產(chǎn)生的隨機(jī)分岔現(xiàn)象和由于隨機(jī)吸引子與其吸引域邊界發(fā)生碰撞而產(chǎn)生的隨機(jī)分岔現(xiàn)象.研究結(jié)果揭示了隨機(jī)干擾對(duì)永磁同步風(fēng)力發(fā)電機(jī)運(yùn)行性能影響的機(jī)理,為永磁同步風(fēng)力發(fā)電系統(tǒng)的運(yùn)行和設(shè)計(jì)提供了理論依據(jù).

永磁同步風(fēng)力發(fā)電機(jī),隨機(jī)分岔,非線性,胞映射

1 引 言

近年來(lái)風(fēng)力發(fā)電快速發(fā)展,永磁同步發(fā)電機(jī)作為主流機(jī)型在風(fēng)力發(fā)電領(lǐng)域中被廣泛應(yīng)用,因此對(duì)永磁同步發(fā)電機(jī)穩(wěn)定性和可靠性的研究受到了人們的廣泛關(guān)注.永磁同步電機(jī)是一個(gè)強(qiáng)非線性系統(tǒng),其中必然會(huì)存在豐富的非線性動(dòng)力學(xué)行為.近幾年來(lái),國(guó)內(nèi)外許多文獻(xiàn)應(yīng)用非線性動(dòng)力學(xué)理論和方法對(duì)永磁同步電機(jī)進(jìn)行了深入研究［1?9］.文獻(xiàn)［1,2］發(fā)現(xiàn)在某些參數(shù)及工作條件下永磁同步電機(jī)中會(huì)出現(xiàn)Hopf分岔和混沌現(xiàn)象,從而導(dǎo)致電機(jī)振蕩失穩(wěn).文獻(xiàn)［3］發(fā)現(xiàn)倍周期分岔是永磁同步電機(jī)進(jìn)入混沌的途徑之一.文獻(xiàn)［4］給出了預(yù)測(cè)永磁同步電機(jī)發(fā)生混沌的參數(shù)空間,為系統(tǒng)的參數(shù)和控制器設(shè)計(jì)提供了理論參考.文獻(xiàn)［5—9］發(fā)現(xiàn)用于風(fēng)力發(fā)電的直驅(qū)永磁同步風(fēng)力發(fā)電機(jī)中也會(huì)出現(xiàn)混沌現(xiàn)象,從而導(dǎo)致系統(tǒng)發(fā)生“無(wú)序”振蕩.上述研究工作揭示了永磁同步電機(jī)失穩(wěn)的動(dòng)力學(xué)機(jī)理,為其設(shè)計(jì)和優(yōu)化提供了理論基礎(chǔ).但是,這些研究工作都是在確定性系統(tǒng)理論的范圍內(nèi)進(jìn)行的,而實(shí)際上永磁同步電機(jī)應(yīng)用于風(fēng)力發(fā)電時(shí),在運(yùn)行過程中不可避免地會(huì)受到風(fēng)能的隨機(jī)干擾.隨機(jī)干擾可能引起發(fā)電機(jī)強(qiáng)烈的非線性隨機(jī)振蕩、失穩(wěn)甚至崩潰.隨機(jī)干擾對(duì)永磁同步風(fēng)力發(fā)電機(jī)動(dòng)力學(xué)行為的影響將直接關(guān)系到對(duì)其穩(wěn)定工作范圍的正確評(píng)估,是一個(gè)不能忽略的實(shí)際問題.因此,研究永磁同步風(fēng)力發(fā)電機(jī)在隨機(jī)干擾下的動(dòng)力學(xué)行為,特別是全局動(dòng)力學(xué)性態(tài)是使其安全可靠運(yùn)行的重要保證.然而,目前尚未見到關(guān)于永磁同步風(fēng)力發(fā)電機(jī)隨機(jī)動(dòng)力學(xué)行為的研究報(bào)道.

為了對(duì)隨機(jī)動(dòng)力系統(tǒng)全局的拓?fù)湫再|(zhì)進(jìn)行分析,就需要對(duì)系統(tǒng)進(jìn)行相空間結(jié)構(gòu)分析.為了克服直接數(shù)值模擬法對(duì)系統(tǒng)全局分析時(shí)效率低、存在有限步截?cái)鄦栴}的不足,Hsu［10］在20世紀(jì)80年代提出了胞映射方法,該方法一步轉(zhuǎn)移矩陣決定了相空間的最終性態(tài),無(wú)時(shí)間截?cái)?是一種對(duì)隨機(jī)系統(tǒng)進(jìn)行全局分析的有效工具［11］.近年來(lái),在最早提出的簡(jiǎn)單胞映射方法的基礎(chǔ)上發(fā)展形成了多種改進(jìn)方法［12?15］,胞映射方法的理論已逐漸成熟.目前,已有文獻(xiàn)采用胞映射方法對(duì)非線性系統(tǒng)的隨機(jī)動(dòng)力學(xué)行為進(jìn)行全局分析,但大多是針對(duì)一些經(jīng)典模型［16?18］,將其應(yīng)用于分析實(shí)際系統(tǒng)的研究工作并不多見.

本文考慮風(fēng)能的隨機(jī)性對(duì)永磁同步風(fēng)力發(fā)電機(jī)全局動(dòng)力學(xué)特性的影響,在建立永磁同步風(fēng)力發(fā)電機(jī)非線性隨機(jī)模型的基礎(chǔ)上,將胞映射方法應(yīng)用于該系統(tǒng)非線性隨機(jī)動(dòng)力學(xué)行為的研究,分析隨機(jī)干擾作用下系統(tǒng)全局結(jié)構(gòu)性態(tài)的演化,給出系統(tǒng)中隨機(jī)分岔的發(fā)生和發(fā)展過程,揭示隨機(jī)干擾對(duì)永磁同步風(fēng)力發(fā)電機(jī)運(yùn)行性能影響的動(dòng)力學(xué)機(jī)理.

2 隨機(jī)輸入轉(zhuǎn)矩作用下永磁同步風(fēng)力發(fā)電機(jī)的數(shù)學(xué)模型

直驅(qū)式永磁同步風(fēng)力發(fā)電機(jī)組的基本結(jié)構(gòu)如圖1所示,機(jī)組由風(fēng)力機(jī)、永磁同步電機(jī)、電力電子變換器及其控制器組成.風(fēng)力機(jī)在風(fēng)能的作用下旋轉(zhuǎn)并帶動(dòng)永磁同步發(fā)電機(jī)轉(zhuǎn)動(dòng),當(dāng)達(dá)到切入風(fēng)速后永磁同步發(fā)電機(jī)將機(jī)械能轉(zhuǎn)化為電能并通過電力電子變換器將符合電網(wǎng)頻率和電壓要求的電能饋入電網(wǎng).為了簡(jiǎn)化分析,忽略電力電子變換器及其控制器的動(dòng)態(tài)過程,得到在轉(zhuǎn)子磁場(chǎng)定向坐標(biāo)系(d-q坐標(biāo)系)下永磁同步風(fēng)力發(fā)電機(jī)的動(dòng)態(tài)方程為

其中,isd,isq和Usd,Usq分別為定子電流和定子電壓的d軸和q軸分量;Rs為定子電阻;Lsd和Lsq分別為d軸和q軸電感,ψr為轉(zhuǎn)子磁鏈,ωr為發(fā)電機(jī)轉(zhuǎn)子轉(zhuǎn)速,Tm為機(jī)械轉(zhuǎn)矩,J為發(fā)電機(jī)的轉(zhuǎn)動(dòng)慣量,β為黏性阻尼系數(shù),np為轉(zhuǎn)子極對(duì)數(shù).

假定發(fā)電機(jī)氣隙均勻,則Lsd=Lsq,通過仿射和時(shí)間變換［1］,可得(1)式的無(wú)量綱狀態(tài)方程為

風(fēng)能的隨機(jī)波動(dòng)會(huì)引起永磁同步風(fēng)力發(fā)電機(jī)的機(jī)械轉(zhuǎn)矩呈隨機(jī)變化.考慮風(fēng)能的隨機(jī)性,永磁同步風(fēng)力發(fā)電機(jī)在隨機(jī)輸入轉(zhuǎn)矩作用下的數(shù)學(xué)模型為

式中,ξ(t)為標(biāo)準(zhǔn)高斯過程,k為隨機(jī)強(qiáng)度.因?yàn)樵谳^短的時(shí)間內(nèi),轉(zhuǎn)矩是圍繞某一個(gè)均值波動(dòng)的,所以這種波動(dòng)在一般情況下可以近似為高斯過程.

3 隨機(jī)輸入轉(zhuǎn)矩作用下永磁同步風(fēng)力發(fā)電機(jī)的隨機(jī)分岔

3.1 隨機(jī)分岔的定義

隨機(jī)干擾導(dǎo)致的系統(tǒng)定性性態(tài)的變化稱為隨機(jī)分岔.對(duì)于隨機(jī)分岔的研究尚處于起步階段,其理論體系尚不完善,缺少嚴(yán)格的一般性定理與準(zhǔn)則.目前,描述隨機(jī)分岔的方式主要有三類:系統(tǒng)穩(wěn)態(tài)概率密度形狀的突變(P-分岔)［19］,系統(tǒng)最大Lyapunov指數(shù)符號(hào)的變化(D-分岔)［19］和系統(tǒng)隨機(jī)吸引子(或隨機(jī)鞍)形態(tài)的突變［16］.其中,第三類定義可以表征在隨機(jī)干擾作用下隨機(jī)系統(tǒng)的全局拓?fù)湫再|(zhì)變化,并且能夠描述具有多吸引子的系統(tǒng)在隨機(jī)激勵(lì)作用下動(dòng)力學(xué)特性的變化,而前兩類隨機(jī)分岔的定義對(duì)此幾乎無(wú)法討論［11］.因此,針對(duì)永磁同步風(fēng)力發(fā)電機(jī)的系統(tǒng)特征和進(jìn)行全局分析的要求,本文基于系統(tǒng)隨機(jī)吸引子(或隨機(jī)鞍)形態(tài)(包括尺寸、周期)的突變這一定義進(jìn)行隨機(jī)分岔分析.

3.2 胞映射方法

胞映射的基本思想是把動(dòng)力系統(tǒng)的狀態(tài)空間離散化為大量狀態(tài)胞的集合,系統(tǒng)中狀態(tài)的轉(zhuǎn)移對(duì)應(yīng)胞與胞的轉(zhuǎn)移,通過對(duì)胞之間轉(zhuǎn)移關(guān)系的研究來(lái)描述原系統(tǒng)的動(dòng)力學(xué)行為［11］.在眾多版本的胞映射方法中,廣義胞映射圖論方法因具有顯著優(yōu)點(diǎn)而被認(rèn)為是一種用于研究隨機(jī)動(dòng)力系統(tǒng)的有效工具.特別是該方法可以得到動(dòng)力系統(tǒng)的不穩(wěn)定解集(鞍),這是其他方法不易得到的.為了更全面地反映永磁同步風(fēng)力發(fā)電機(jī)的隨機(jī)動(dòng)力學(xué)全局特性,本文采用廣義胞映射圖論方法作為基本胞映射方法.由于本文研究的永磁同步風(fēng)力發(fā)電機(jī)的模型為3階,因此為了降低需要計(jì)算的胞空間的維數(shù)、減小計(jì)算量,在進(jìn)行胞映射數(shù)值計(jì)算之前,先采用Poincarè型胞映射方法［20］將系統(tǒng)轉(zhuǎn)化為Poincarè截面上的點(diǎn)映射系統(tǒng),再求出系統(tǒng)在Poincarè截面上的胞映射關(guān)系.在分析多吸引子共存的情況時(shí),由于系統(tǒng)的全局特性較復(fù)雜,為了既能準(zhǔn)確反映系統(tǒng)的復(fù)雜空間結(jié)構(gòu)又能提高計(jì)算效率,本文采用基于復(fù)合胞化空間的迭代方法［21］,根據(jù)狀態(tài)空間的不同采用不同尺度的胞化方案,以得到更精細(xì)的胞映射動(dòng)力系統(tǒng).本文胞映射算法實(shí)現(xiàn)的主要過程為:

1)選取合適的Poincarè截面Σ, 將(3)式轉(zhuǎn)化為Poincarè截面Σ上的點(diǎn)映射系統(tǒng)x(n+1)=g(x(n)),其中

2)在Poincarè截面Σ上選取研究范圍R,R為二維平面,在R內(nèi)分別沿坐標(biāo)軸x1和x2劃分為間距分別為h1和h2的N維胞狀態(tài)空間,記為Γ={z1,z2,···,zN},其中zi為胞狀態(tài)空間中的胞元素;

3)在每個(gè)胞內(nèi)取M個(gè)采樣點(diǎn),對(duì)每個(gè)點(diǎn)記其坐標(biāo)為(x1j,x2j),以每個(gè)點(diǎn)為初值,經(jīng)過映射g后的像點(diǎn)為根據(jù)每個(gè)胞內(nèi)采樣點(diǎn)和其像點(diǎn)的關(guān)系,求出每個(gè)胞的像胞及一步轉(zhuǎn)移概率pij.pij構(gòu)成轉(zhuǎn)移概率矩陣P,用pj(n)表示系統(tǒng)在第n步映射中位于胞j的概率,p(n)是以pj(n)為元素的N維矢量,則得到廣義胞映射關(guān)系p(n+1)=P·p(n);

4)將每個(gè)胞對(duì)應(yīng)于有向圖的一個(gè)頂點(diǎn),若兩個(gè)胞之間存在一步可達(dá)關(guān)系,則在相應(yīng)的頂點(diǎn)之間建立一條有向邊,從而根據(jù)P建立有向圖D=(V,E),其中V表示頂點(diǎn)集合,E表示有向邊的集合;

5)對(duì)形成的有向圖D=(V,E)進(jìn)行分析得到強(qiáng)聯(lián)通子圖,根據(jù)強(qiáng)聯(lián)通子圖的特征得到系統(tǒng)的吸引子、吸引域、鞍和吸引域邊界等信息;

6)若所分析的全局特性較復(fù)雜,則適當(dāng)選取需要精細(xì)分割的胞狀態(tài)集合Γ1,對(duì)Γ1進(jìn)行細(xì)分為n1×n2的復(fù)合胞化空間.根據(jù)上述方法生成有向圖D′=(V′,E′),并對(duì)其進(jìn)行強(qiáng)聯(lián)通分析以得到相應(yīng)的結(jié)果;

7)如有必要,返回6),否則停止.

在分析多吸引子共存的系統(tǒng)的全局特性時(shí),判斷每個(gè)吸引子的穩(wěn)定裕度是十分必要的.本文采用廣義胞映射中保護(hù)層厚度這個(gè)概念來(lái)定量刻畫吸引子被其吸引域保護(hù)的程度,即可反映吸引子的穩(wěn)定裕度.通過比較保護(hù)層厚度這個(gè)值的大小來(lái)判斷在多個(gè)吸引子共存的系統(tǒng)中隨著隨機(jī)強(qiáng)度的增加,哪個(gè)吸引子將會(huì)首先消失.第i個(gè)永久族(吸引子)的保護(hù)層厚度的定義為［22］

它表示第i個(gè)永久族集合A(i)與多駐處胞集合(吸引域邊界)B(i,j)之間的距離.

3.3 隨機(jī)內(nèi)部激變

圖2 (網(wǎng)刊彩色)參數(shù)選取sd=?75,sq=?15,m=?50時(shí),隨著隨機(jī)強(qiáng)度k的變化永磁同步風(fēng)力發(fā)電機(jī)全局性態(tài)的變化情況(紅色區(qū)域?yàn)槲覣(1),灰色區(qū)域?yàn)锳(1)的吸引域B(1),藍(lán)色區(qū)域?yàn)榘癝(1)) (a)k=0;(b)k=0.11;(c)k=0.12;(d)k=0.13;(e)k=0.14;(f)k=0.15Fig.2.(color online)The changes in global character of the permanent magnet synchronous generator with the variations of the stochastic intensity k,when the parameters are set assd=?75,sq=?15,m=?50(the attractor A(1)is marked as red area,its attraction basin B(1)is marked as gray area,the saddle S(1)is marked as blue area):(a)k=0;(b)k=0.11;(c)k=0.12;(d)k=0.13;(e)k=0.14;(f)k=0.15.

當(dāng)k=0時(shí),即輸入機(jī)械轉(zhuǎn)矩不存在隨機(jī)擾動(dòng)時(shí)永磁同步風(fēng)力發(fā)電機(jī)的全局性態(tài)如圖2(a)所示.此時(shí),系統(tǒng)有1個(gè)吸引子A(1)和1個(gè)鞍S(1),且鞍存在于吸引子的吸引域B(1)內(nèi)部.考慮隨機(jī)強(qiáng)度k的變化對(duì)永磁同步風(fēng)力發(fā)電機(jī)全局特性產(chǎn)生的影響.隨著隨機(jī)強(qiáng)度k的增大,隨機(jī)吸引子和隨機(jī)鞍變大.當(dāng)k=0.11時(shí),系統(tǒng)的全局特性如圖2(b)所示,隨機(jī)吸引子逐漸伸向隨機(jī)鞍,其指端即將接觸到隨機(jī)鞍的右上半支.進(jìn)一步增大k,當(dāng)k=0.12時(shí),隨機(jī)吸引子與隨機(jī)鞍發(fā)生碰撞,導(dǎo)致隨機(jī)鞍的右上半支消失,成為隨機(jī)吸引子的一部分,吸引子的形態(tài)發(fā)生突變,系統(tǒng)發(fā)生了隨機(jī)分岔,如圖2(c)所示.繼續(xù)增大隨機(jī)強(qiáng)度,隨機(jī)吸引子再次逐漸靠近隨機(jī)鞍的右下半支和左半支,如圖2(d)和圖2(e)所示.當(dāng)隨機(jī)強(qiáng)度k由0.14增至0.15時(shí),隨機(jī)吸引子與隨機(jī)鞍再次發(fā)生碰撞,隨機(jī)鞍突然消失并成為隨機(jī)吸引子的一部分,在關(guān)心的區(qū)域內(nèi)僅有一個(gè)大的隨機(jī)吸引子,如圖2(f)所示.此時(shí)鞍變?yōu)槲?吸引子的幾何空間體積突然增大,形態(tài)發(fā)生了突變,系統(tǒng)再次發(fā)生隨機(jī)分岔.表1給出了隨著隨機(jī)強(qiáng)度k的變化,系統(tǒng)的吸引子和鞍的胞的變化情況.這種隨機(jī)分岔實(shí)際上是一種由于隨機(jī)吸引子與其吸引域內(nèi)的隨機(jī)鞍發(fā)生碰撞而產(chǎn)生的隨機(jī)激變現(xiàn)象,即隨機(jī)內(nèi)部激變［24］.

由上述永磁同步風(fēng)力發(fā)電機(jī)全局性態(tài)的變化可以看出,當(dāng)輸入轉(zhuǎn)矩存在較小強(qiáng)度的隨機(jī)擾動(dòng)時(shí),永磁同步風(fēng)力發(fā)電機(jī)會(huì)圍繞一個(gè)穩(wěn)定運(yùn)行點(diǎn)做小幅的隨機(jī)振蕩,隨著隨機(jī)強(qiáng)度的增大,隨機(jī)吸引子的胞數(shù)增多,以致隨機(jī)振蕩的幅值也會(huì)增大.當(dāng)隨機(jī)吸引子與隨機(jī)鞍首次碰撞而發(fā)生隨機(jī)分岔時(shí),隨機(jī)吸引子胞數(shù)的突增導(dǎo)致永磁同步風(fēng)力發(fā)電機(jī)會(huì)突然發(fā)生更為明顯的隨機(jī)振蕩.當(dāng)隨機(jī)吸引子與隨機(jī)鞍再次碰撞,隨機(jī)鞍突然消失,隨機(jī)吸引子形態(tài)突變,系統(tǒng)再次發(fā)生隨機(jī)分岔時(shí),永磁同步風(fēng)力發(fā)電機(jī)將突然劇烈、無(wú)序地振蕩,甚至發(fā)電機(jī)的轉(zhuǎn)速r會(huì)出現(xiàn)負(fù)值,系統(tǒng)失去穩(wěn)定性.可見,隨機(jī)擾動(dòng)強(qiáng)度增大而引起的隨機(jī)內(nèi)部激變是導(dǎo)致永磁同步風(fēng)力發(fā)電機(jī)失穩(wěn)的動(dòng)力學(xué)本質(zhì).

表1 隨著隨機(jī)強(qiáng)度k的變化吸引子和鞍的胞數(shù)目的變化情況Table 1.The changes in the number of the attractor and saddle with the variations of the stochastic intensity k.

3.4 隨機(jī)邊界激變

當(dāng)k=0時(shí),即不存在隨機(jī)擾動(dòng)時(shí)永磁同步風(fēng)力發(fā)電機(jī)的全局性態(tài)如圖3(a)所示.此時(shí),系統(tǒng)有2個(gè)共存的吸引子A(1)和A(2),2個(gè)吸引子對(duì)應(yīng)的2個(gè)吸引域B(1)和B(2)的交界處存在多駐處胞,即吸引域邊界B(1,2).根據(jù)(4)式可計(jì)算出吸引子A(1)的保護(hù)層厚度為d1=18.35,吸引子A(2)的保護(hù)層厚度為d2=310.82.由于d1＜d2,因此可預(yù)判隨著隨機(jī)強(qiáng)度的增大,吸引子A(1)將會(huì)首先消失.

考慮隨機(jī)強(qiáng)度k的變化對(duì)永磁同步風(fēng)力發(fā)電機(jī)全局特性產(chǎn)生的影響.隨著隨機(jī)強(qiáng)度k的增大,2個(gè)隨機(jī)吸引子和吸引域邊界都逐漸增大,如圖3(b)所示.當(dāng)k=0.0033時(shí),系統(tǒng)的全局特性如圖3(c)所示,隨機(jī)吸引子A(1)逐漸伸向吸引域邊界,其指端即將接觸到吸引域邊界.進(jìn)一步增大k,當(dāng)k=0.0034時(shí),隨機(jī)吸引子A(1)與吸引域邊界B(1,2)發(fā)生碰撞,該吸引子連同其吸引域及吸引域邊界消失,而在消失的吸引子的位置上出現(xiàn)了一個(gè)隨機(jī)鞍S(1),此時(shí)在研究的區(qū)域內(nèi)系統(tǒng)僅存在一個(gè)隨機(jī)吸引子A(2).由于吸引子的數(shù)量和性質(zhì)都發(fā)生了突變,系統(tǒng)發(fā)生隨機(jī)分岔,如圖3(d)所示.這里驗(yàn)證了前面得到的結(jié)論,即當(dāng)隨機(jī)強(qiáng)度增大到一定程度時(shí),隨機(jī)吸引子A(1)首先消失.表2給出了隨著隨機(jī)強(qiáng)度的變化,系統(tǒng)的吸引子和鞍的胞數(shù)的變化情況.上述隨機(jī)分岔是一種由于隨機(jī)吸引子與其吸引域邊界發(fā)生碰撞而產(chǎn)生的隨機(jī)激變現(xiàn)象,即隨機(jī)邊界激變［24］.

由上述永磁同步風(fēng)力發(fā)電機(jī)全局性態(tài)的變化可以看出,當(dāng)輸入轉(zhuǎn)矩中的隨機(jī)擾動(dòng)強(qiáng)度較小時(shí),永磁同步風(fēng)力發(fā)電機(jī)會(huì)圍繞一個(gè)穩(wěn)定運(yùn)行點(diǎn)做小幅的隨機(jī)振蕩,隨著隨機(jī)強(qiáng)度的增大,隨機(jī)振蕩的幅值逐漸增大,且對(duì)該穩(wěn)定工作點(diǎn)起到保護(hù)作用的吸引域邊界也在增大,導(dǎo)致隨機(jī)振蕩的范圍逐漸靠近吸引域邊界,系統(tǒng)的穩(wěn)定裕度變小.當(dāng)此隨機(jī)吸引子與其吸引域邊界碰撞而發(fā)生隨機(jī)分岔后,原來(lái)的穩(wěn)定解集(吸引子A(1))消失,取而代之的是不穩(wěn)定解集(鞍S(1)),在研究的區(qū)域內(nèi)僅存在一個(gè)隨機(jī)吸引子A(2).由于僅存的隨機(jī)吸引子所對(duì)應(yīng)的發(fā)電機(jī)轉(zhuǎn)速為負(fù),所以對(duì)于工程實(shí)際中的永磁同步風(fēng)力發(fā)電機(jī)來(lái)說(shuō),發(fā)電機(jī)將發(fā)生運(yùn)行的突然崩潰.可見,隨機(jī)擾動(dòng)強(qiáng)度增大而引起的隨機(jī)外部激變是導(dǎo)致永磁同步風(fēng)力發(fā)電機(jī)運(yùn)行崩潰的動(dòng)力學(xué)本質(zhì).

圖3 (網(wǎng)刊彩色)參數(shù)選取sd=10,sq=?50,m=?55時(shí),隨著隨機(jī)強(qiáng)度k的變化永磁同步風(fēng)力發(fā)電機(jī)全局形態(tài)的變化情況(紅色區(qū)域?yàn)槲覣(1)和A(2),白色區(qū)域?yàn)锳(1)的吸引域B(1),灰色區(qū)域?yàn)锳(2)的吸引域B(2),綠色區(qū)域?yàn)槲蜻吔鏐(1,2),藍(lán)色區(qū)域?yàn)榘癝(1)) (a)k=0;(b)k=0.0025;(c)k=0.0033;(d)k=0.0034Fig.3.(color online)The changes in global character of the permanent magnet synchronous generator for wind turbine with the variations of the stochastic intensity k,when the parameters aresd=10,sq=?50,m=?55(the attractor A(1)and A(2)is marked as red area,the attraction basin B(1)of A(1)is marked as white area,the attraction basin B(2)of A(2)is marked as gray area,the boundary of the attraction basin B(1,2)is marked as green area,the saddle S(1)is marked as blue area):(a)k=0;(b)k=0.0025;(c)k=0.0033;(d)k=0.0034.

表2 隨著隨機(jī)強(qiáng)度k的變化,吸引子、鞍和吸引域邊界的胞數(shù)目的變化情況Table 2.The changes in the number of the attractor,saddle and boundary of the attraction basin with the variations of the stochastic intensity k.

4 數(shù)值模擬

為了更直觀地描述第3節(jié)中所分析的永磁同步風(fēng)力發(fā)電機(jī)的隨機(jī)分岔現(xiàn)象,并驗(yàn)證采用胞映射方法得到的分析結(jié)果的正確性,本節(jié)采用Euler-Maruyama(EM)方法對(duì)隨機(jī)輸入轉(zhuǎn)矩作用下永磁同步風(fēng)力發(fā)電機(jī)的隨機(jī)響應(yīng)進(jìn)行數(shù)值模擬.從統(tǒng)計(jì)學(xué)角度來(lái)看,用EM方法求解隨機(jī)微分方程組時(shí),每一條維納過程路徑上的解過程可以看成是總體中的一個(gè)抽樣樣本.為了闡釋數(shù)值模擬的結(jié)果,本節(jié)中分別給出仿真結(jié)束時(shí)刻發(fā)電機(jī)轉(zhuǎn)速的分布直方圖［25］(仿真結(jié)束時(shí)刻轉(zhuǎn)速的所有隨機(jī)路徑在不同轉(zhuǎn)速區(qū)間的個(gè)數(shù))以描述轉(zhuǎn)速的分布情況;轉(zhuǎn)速的樣本均值［25］曲線(所有隨機(jī)路徑上解過程的平均值)、第一四分位數(shù)和第三四分位數(shù)［25］曲線(每個(gè)時(shí)刻所有隨機(jī)路徑的數(shù)據(jù)排序后處于25%和75%位置上的值所構(gòu)成的曲線)以描述隨機(jī)響應(yīng)的波動(dòng)情況;某一隨機(jī)路徑的時(shí)域波形圖.

4.1 隨機(jī)內(nèi)部激變

永磁同步風(fēng)力發(fā)電機(jī)參數(shù)γ和σ、定子電壓sd和sq以及機(jī)械轉(zhuǎn)矩m的取值與3.3節(jié)相同,仿真時(shí)間50 s,取2000個(gè)隨機(jī)路徑,初值取在吸引子A(1)附近,為當(dāng)隨機(jī)強(qiáng)度k=0.1,系統(tǒng)未發(fā)生隨機(jī)分岔時(shí),發(fā)電機(jī)轉(zhuǎn)速的隨機(jī)響應(yīng)數(shù)值模擬結(jié)果如圖4所示.可以看出,發(fā)電機(jī)的轉(zhuǎn)速分布于13.9—14.7之間,風(fēng)力發(fā)電機(jī)圍繞r=14.29的平衡點(diǎn)以較小幅度平緩地隨機(jī)波動(dòng),處于穩(wěn)定運(yùn)行狀態(tài).

當(dāng)隨機(jī)強(qiáng)度k=0.5,系統(tǒng)發(fā)生隨機(jī)內(nèi)部激變后,發(fā)電機(jī)轉(zhuǎn)速的隨機(jī)響應(yīng)情況如圖5所示.可以看出,與圖4相比發(fā)電機(jī)轉(zhuǎn)速的隨機(jī)振蕩幅值明顯增大,分布于?20—30之間,在劇烈振蕩過程中轉(zhuǎn)速出現(xiàn)負(fù)值,實(shí)際系統(tǒng)無(wú)法正常運(yùn)行,失去穩(wěn)定性.與圖2對(duì)比可以看出,發(fā)生隨機(jī)分岔后轉(zhuǎn)速的分布與圖2(f)中吸引子所覆蓋的轉(zhuǎn)速范圍相同.圖5所示的失穩(wěn)現(xiàn)象正是由于隨機(jī)吸引子與其吸引域內(nèi)的隨機(jī)鞍發(fā)生碰撞導(dǎo)致隨機(jī)鞍消失、隨機(jī)吸引子的幾何空間體積突然增大而引起的.

圖4 參數(shù)選取sd=?75,sq=?15,m=?50,系統(tǒng)穩(wěn)定時(shí)永磁同步風(fēng)力發(fā)電機(jī)轉(zhuǎn)速的隨機(jī)響應(yīng)(a)t=50 s時(shí)在2000條路徑上的分布直方圖;(b)在2000條路徑上解過程的樣本均值(Mean)曲線、第一四分位數(shù)(Q1)和第三四分位數(shù)(Q3)曲線;(c)某一隨機(jī)路徑的時(shí)域波形Fig.4.The stochastic responses of the rotor speed of the permanent magnet synchronous generator for wind turbine when the system is stable,with the parameterssd=?75,sq=?15,m=?50:(a)The histogram at t=50 s over 2000 trajectories;(b)the mean and the fi rst-third quartile quartiles over 2000 trajectories;(c)the waveform of one sample trajectory.

4.2 隨機(jī)邊界激變

發(fā)電機(jī)參數(shù)γ和σ、定子電壓sd和sq以及機(jī)械轉(zhuǎn)矩m的取值與3.4節(jié)相同,仿真時(shí)間為50 s,取10000個(gè)隨機(jī)路徑,初值設(shè)在吸引子A(1)的吸引域B(1)內(nèi),為當(dāng)隨機(jī)強(qiáng)度k=0.0025,系統(tǒng)未發(fā)生隨機(jī)分岔時(shí),發(fā)電機(jī)轉(zhuǎn)速的隨機(jī)響應(yīng)情況如圖6所示.可以看出,由于初值設(shè)在吸引子A(1)的吸引域內(nèi),因此發(fā)電機(jī)圍繞r=8.07的平衡點(diǎn)以較小幅度隨機(jī)振蕩,轉(zhuǎn)速的振蕩幅值在7.85—8.3之間,永磁同步風(fēng)力發(fā)電機(jī)處于穩(wěn)定運(yùn)行狀態(tài).

圖5 參數(shù)選取sd=?75,sq=?15,m=?50,發(fā)生隨機(jī)分岔后永磁同步風(fēng)力發(fā)電機(jī)轉(zhuǎn)速的隨機(jī)響應(yīng) (a)t=50 s時(shí)在2000條路徑上的分布直方圖;(b)在2000條路徑上解過程的樣本均值(Mean)曲線、第一四分位數(shù)(Q1)和第三四分位數(shù)(Q3)曲線;(c)某一隨機(jī)路徑的時(shí)域波形Fig.5.The stochastic responses of the rotor speed of the permanent magnet synchronous generator for wind turbine after the stochastic bifurcation happens,with the parameterssd=?75,sq=?15,m=?50:(a)The histogram at t=50 s over 2000 trajectories;(b)the mean and the fi rst-third quartile quartiles over 2000 trajectories;(c)the waveform of one sample trajectory.

圖6 參數(shù)選取sd=10,sq=?50,m=?55,系統(tǒng)穩(wěn)定時(shí)永磁同步風(fēng)力發(fā)電機(jī)轉(zhuǎn)速的隨機(jī)響應(yīng) (a)t=50 s時(shí)在10000條路徑上的分布直方圖;(b)在10000條路徑上解過程的樣本均值(Mean)曲線、第一四分位數(shù)(Q1)和第三四分位數(shù)(Q3)曲線;(c)某一隨機(jī)路徑的時(shí)域波形Fig.6.The stochastic responses of the rotor speed of the permanent magnet synchronous generator for wind turbine when the system is stable,with the parameterssd=10,sq=?50,m=?55:(a)The histogram at t=50 s over 10000 trajectories;(b)the mean and the fi rst-third quartile quartiles over 10000 trajectories;(c)the waveform of one sample trajectory.

當(dāng)隨機(jī)強(qiáng)度k=0.015,系統(tǒng)發(fā)生隨機(jī)邊界激變后,發(fā)電機(jī)轉(zhuǎn)速的隨機(jī)響應(yīng)情況如圖7所示.可見,當(dāng)t=50 s時(shí)發(fā)電機(jī)轉(zhuǎn)速r在10000條路徑上的分布直方圖出現(xiàn)了雙峰,轉(zhuǎn)速的樣本均值曲線不斷下降.具體到某一隨機(jī)路徑,發(fā)電機(jī)在運(yùn)行一段時(shí)間后,轉(zhuǎn)速突然降為負(fù)值,實(shí)際系統(tǒng)失穩(wěn)崩潰.與圖3對(duì)比可以看出,轉(zhuǎn)速分布直方圖中兩個(gè)雙峰的中心r=8.07和r=?1.44,分別與圖3(a)中吸引子A(1)和A(2)所對(duì)應(yīng)的轉(zhuǎn)速符合.當(dāng)隨機(jī)強(qiáng)度增大而導(dǎo)致隨機(jī)吸引子A(1)與其吸引域邊界碰撞而發(fā)生隨機(jī)分岔后,原來(lái)的穩(wěn)定解集(吸引子A(1))消失并變?yōu)椴环€(wěn)定解集(鞍S(1)),在研究的區(qū)域內(nèi)僅存在一個(gè)隨機(jī)吸引子A(2).由于數(shù)值模擬的初值設(shè)在原吸引子A(1)的吸引域內(nèi),即發(fā)生隨機(jī)分岔后的鞍S(1)區(qū)域內(nèi),因此在運(yùn)行一段時(shí)間后,發(fā)電機(jī)的運(yùn)行狀態(tài)被吸引子A(2)所吸引,以A(2)為平衡點(diǎn)做隨機(jī)振蕩.由于吸引子A(2)所對(duì)應(yīng)的轉(zhuǎn)速為負(fù),因此永磁同步風(fēng)力發(fā)電機(jī)的運(yùn)行崩潰.

圖7 參數(shù)選取sd=10,sq=?50,m=?55,發(fā)生隨機(jī)分岔后永磁同步風(fēng)力發(fā)電機(jī)轉(zhuǎn)速的隨機(jī)響應(yīng) (a)t=50 s時(shí)在10000條路徑上的分布直方圖;(b)在10000條路徑上解過程的樣本均值(Mean)曲線、第一四分位數(shù)(Q1)和第三四分位數(shù)(Q3)曲線;(c)某一隨機(jī)樣本的時(shí)域波形Fig.7.The stochastic responses of the rotor speed of the permanent magnet synchronous generator for wind turbine after the stochastic bifurcation happens,with the parameterssd=?75,sq=?15,m=?50:(a)The histogram at t=50 s over 10000 trajectories;(b)the mean and the fi rst-third quartile quartiles over 10000 trajectories;(c)the waveform of one sample trajectory.

4.3 結(jié)果討論

需要說(shuō)明的是,一方面,與確定性系統(tǒng)的分岔行為不同,由于受到隨機(jī)因素的作用,隨機(jī)系統(tǒng)即使?jié)M足一定的分岔?xiàng)l件也并非一定會(huì)發(fā)生隨機(jī)分岔,而是具有發(fā)生分岔的可能.因此,第3節(jié)中采用胞映射法分析得到的永磁同步風(fēng)力發(fā)電機(jī)發(fā)生分岔時(shí)隨機(jī)強(qiáng)度的臨界值僅能表示當(dāng)隨機(jī)強(qiáng)度大于臨界值時(shí)系統(tǒng)發(fā)生隨機(jī)分岔的概率較大.其意義在于可指導(dǎo)系統(tǒng)的運(yùn)行和設(shè)計(jì)從而降低發(fā)生隨機(jī)分岔的概率,但無(wú)法通過此臨界值對(duì)隨機(jī)系統(tǒng)的分岔做出確定性的預(yù)測(cè).

另一方面,由于對(duì)隨機(jī)響應(yīng)的直接數(shù)值模擬需要計(jì)算大量樣本軌線的長(zhǎng)期行為才能得到更精確的結(jié)果,而受到計(jì)算機(jī)計(jì)算能力的限制,總是對(duì)有限的樣本進(jìn)行有限步截?cái)?因此直接數(shù)值模擬得到的結(jié)果必然存在誤差.

綜上,本節(jié)對(duì)隨機(jī)響應(yīng)的數(shù)值模擬驗(yàn)證了永磁同步風(fēng)力發(fā)電機(jī)中隨機(jī)內(nèi)部激變和隨機(jī)邊界激變的存在性.然而,由于隨機(jī)系統(tǒng)的復(fù)雜性和直接數(shù)值模擬法的不足,無(wú)法對(duì)發(fā)生隨機(jī)分岔時(shí)的隨機(jī)強(qiáng)度臨界值進(jìn)行準(zhǔn)確的驗(yàn)證.

5 結(jié) 論

本文針對(duì)風(fēng)能的隨機(jī)性,建立了在輸入機(jī)械轉(zhuǎn)矩存在隨機(jī)干擾情況下永磁同步風(fēng)力發(fā)電機(jī)的非線性隨機(jī)模型,采用胞映射方法分析了隨機(jī)干擾強(qiáng)度變化時(shí)系統(tǒng)全局結(jié)構(gòu)的變化過程及出現(xiàn)的隨機(jī)分岔現(xiàn)象,并通過對(duì)系統(tǒng)的隨機(jī)響應(yīng)進(jìn)行數(shù)值模擬驗(yàn)證了理論分析的正確性.研究結(jié)果表明,隨著隨機(jī)干擾強(qiáng)度的增大,系統(tǒng)中出現(xiàn)了兩種隨機(jī)分岔現(xiàn)象:

1)當(dāng)隨機(jī)強(qiáng)度增至臨界值時(shí),隨機(jī)吸引子與其吸引域內(nèi)的隨機(jī)鞍發(fā)生碰撞,導(dǎo)致隨機(jī)鞍突然消失、隨機(jī)吸引子的幾何空間體積突增、形態(tài)發(fā)生突變,在此變化過程中,永磁同步風(fēng)力發(fā)電機(jī)發(fā)生強(qiáng)烈的隨機(jī)振蕩并最終失去穩(wěn)定性;

2)當(dāng)隨機(jī)強(qiáng)度增至臨界值時(shí),一個(gè)隨機(jī)吸引子與其吸引域邊界發(fā)生碰撞,該吸引子連同其吸引域及吸引域邊界消失,在原吸引子的位置上出現(xiàn)了一個(gè)隨機(jī)鞍,系統(tǒng)吸引子的數(shù)量發(fā)生了突變,在此變化過程中,永磁同步風(fēng)力發(fā)電機(jī)的隨機(jī)振蕩幅值不斷增大,最終由于穩(wěn)定解集的消失,系統(tǒng)失穩(wěn)崩潰.

總之,小強(qiáng)度的隨機(jī)干擾就會(huì)導(dǎo)致系統(tǒng)穩(wěn)定性結(jié)構(gòu)的完全破壞,使原本處于穩(wěn)定狀態(tài)的系統(tǒng)發(fā)生強(qiáng)烈的無(wú)序振蕩和運(yùn)行崩潰.因此,本文研究隨機(jī)干擾對(duì)系統(tǒng)動(dòng)力學(xué)行為演化的影響對(duì)于系統(tǒng)的安全可靠運(yùn)行具有重要的指導(dǎo)意義.

［1］Zhang B,Li Z,Mao Z Y,Pang M X 2001Proc.CSEE21 13(in Chinese)［張波,李忠,毛宗源,龐敏熙2001中國(guó)電機(jī)工程學(xué)報(bào)21 13］

［2］Xue W,Guo Y L,Chen Z Q 2009Acta Phys.Sin.58 8146(in Chinese)［薛薇,郭彥嶺,陳增強(qiáng) 2009物理學(xué)報(bào)58 8146］

［3］Wei D Q,Luo X S,Fang J Q,Wang B H 2006Acta Phys.Sin.55 54(in Chinese)［韋篤取,羅曉曙,方錦清,汪秉宏2006物理學(xué)報(bào)55 54］

［4］Rasoolzadeh A,Tavazoei M S 2012Phys.Lett.A377 73

［5］Yang G L,Li H G 2009Acta Phys.Sin.58 7552(in Chinese)［楊國(guó)良,李惠光 2009物理學(xué)報(bào) 58 7552］

［6］Zheng G,Zou J X,Xu H B,Qin G 2011Acta Phys.Sin.60 060506(in Chinese)［鄭剛,鄒見效,徐紅兵,秦鋼2011物理學(xué)報(bào)60 060506］

［7］Wang L,Li Y H,Lu G L,Zhu X H 2011Electric Power Automation Equipment31 45(in Chinese)［王磊, 李穎暉,逯國(guó)亮,朱喜華2011電力自動(dòng)化設(shè)備31 45］

［8］Ren L N,Liu F C,Jiao X H,Li J Y 2012Acta Phys.Sin.61 060506(in Chinese)［任麗娜,劉福才,焦曉紅,李俊義2012物理學(xué)報(bào)61 060506］

［9］Messadi M,Mellit A,Kemih K,Ghanes M 2015Chin.Phys.B24 010502

［10］Hsu C S 1980J.Appl.Mech.47 931

［11］Xu W 2013Numerical Analysis Methods for Stochastic Dynamical System(Beijing:Science Press)p43(in Chinese)［徐偉2013非線性隨機(jī)動(dòng)力學(xué)的若干數(shù)值方法及應(yīng)用(北京:科學(xué)出版社)第43頁(yè)］

［12］Hsu C S 1981J.Appl.Mech.48 634

［13］Tongue B H,Gu K 1988J.Sound Vib.125 169

［14］Jiang J,Xu J X 1994Phys.Lett.A188 137

［15］Hsu C S 1995Int.J.Bifurcat.Chaos5 1085

［16］Xu W,He Q,Rong H W,Fang T 2003Acta Phys.Sin.52 1365(in Chinese)［徐偉,賀群,戎海武,方同 2003物理學(xué)報(bào)52 1365］

［17］Xu W,Yue X L 2010Sci.China:Technol.Sci.53 664

［18］Li Z G,Jiang J,Hong L 2015Int.J.Bifurcat.Chaos25 1550109

［19］Arnold L 1998Random Dynamical Systems(Berlin,Heidelberg,New York:Springer)pp34–35

［20］Levitas J,Weller T,Singer J 1994J.Sound Vib.176 641

［21］He Q,Xu W,Li S,Xiao Y Z 2008Acta Phys.Sin.57 4021(in Chinese)［賀群,徐偉,李爽,肖玉柱 2008物理學(xué)報(bào)57 4021］

［22］Gong P L,Xu J X 1998Appl.Math.Mech.19 1087(in Chinese)［龔璞林,徐健學(xué)1998應(yīng)用數(shù)學(xué)和力學(xué)19 1087］

［23］Zhu W Q,Yu J S 1987J.Sound Vib.117 421

［24］Hong L,Xu J X 2002Acta Phys.Sin.51 2694(in Chinese)［洪靈,徐健學(xué) 2002物理學(xué)報(bào) 51 2694］

［25］Guan Y H 2009Statistics(Beijing:Higher Education Press)pp66–83(in Chinese)［管于華 2009統(tǒng)計(jì)學(xué) (北京:高等教育出版社)第66—83頁(yè)］

Global analysis of stochastic bifurcation in permanent magnet synchronous generator for wind turbine system?

Yang Li-Hui?Ge Yang Ma Xi-Kui
(State Key Laboratory of Electrical Insulation and Power Equipment,School of Electrical Engineering,Xian Jiaotong University,Xi’an 710049,China)

The permanent magnet synchronous generator(PMSG)for wind turbine system operating under inevitable stochastic disturbance from wind power is a nonlinear stochastic dynamical system.With the random interaction and nonlinearity,the intense nonlinear stochastic oscillation is likely to happen in such a system,causing the system to be unstable or even collapse.However,the PMSG is usually considered as a deterministic system when analyzing its nonlinear dynamic behaviors in the past researches.Such a simpli fi cation can lead to wrong predictions for the system stability and reliability.This paper aims to discuss the e ff ect of the stochastic disturbance on the nonlinear dynamic behavior of the PMSG.Based on the derived PMSG model considering the stochastic disturbance from the input mechanical torque,the evolution of the system global structure with the stochastic intensity is investigated using the generalized cell mapping digraph method.Meanwhile,the occurrence process and development process of the stochastic bifurcation are illustrated.Based on this global analysis,the intrinsic mechanism for the e ff ect of the stochastic disturbance on the operating performances of the PMSG is revealed.Finally,the numerical simulations based on the Euler-Maruyama algorithm are carried out to validate the results of the theoretical analysis.The results present that as the intensity of the stochastic disturbance increases,two kinds of stochastic bifurcations can be observed in the PMSG system according to the de fi nition of a sudden change in characteristic of the stochastic attractor.One is the stochastic interior crisis that occurs when a stochastic attractor collides with a stochastic saddle in its attraction basin interior,leading to the abrupt increase of the attractor and the disappearance of the saddle.This kind of bifurcation results in the intense stochastic oscillation and instability of the PMSG system.Another stochastic bifurcation is the stochastic boundary crisis which occurs when a stochastic attractor collides with the boundary of its attraction basin and results in the disappearance of the attractor.This sudden change of the number of stochastic attractors induces the stable solution set to vanish and thus the PMSG system to collapse.In a word,even the stochastic disturbance with small intensity may lead to the complete destruction of the stable structure of the PMSG,inducing the system to su ff er a strong disordered oscillation or the operation to collapse.The results of this paper can provide signi fi cant theoretic reference for both practically operating and designing the PMSG for wind turbine systems.

permanent magnet synchronous generator,stochastic bifurcation,nonlinear,cell mapping

5 March 2017;revised manuscript

2 July 2017)

(2017年3月5日收到;2017年7月2日收到修改稿)

10.7498/aps.66.190501

?國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):51207122)和陜西省自然科學(xué)基礎(chǔ)研究計(jì)劃(批準(zhǔn)號(hào):2016JQ5034)資助的課題.

?通信作者.E-mail:lihui.yang@mail.xjtu.edu.cn

?2017中國(guó)物理學(xué)會(huì)Chinese Physical Society

PACS:05.45.–a,88.50.G–

10.7498/aps.66.190501

*Project supported by National Natural Science Foundation of China(Grant No.51207122)and Natural Science Basic Research Plan in Shaanxi Province of China(Grant No.2016JQ5034).

?Corresponding author.E-mail:lihui.yang@mail.xjtu.edu.cn