張 強(qiáng)， 黃宵寧

(南京工程學(xué)院電力工程學(xué)院， 南京 211167)

應(yīng)用數(shù)值積分法計(jì)算電力系統(tǒng)混沌閾值

張 強(qiáng)， 黃宵寧

(南京工程學(xué)院電力工程學(xué)院， 南京 211167)

Melnikov函數(shù)是分析同(異)宿軌道出現(xiàn)混沌的最有效方法，用該函數(shù)的數(shù)值積分法計(jì)算單機(jī)無(wú)窮大電力系統(tǒng)在周期性負(fù)荷擾動(dòng)下的混沌閾值。通過(guò)相應(yīng)無(wú)擾系統(tǒng)的Hamilton函數(shù)求得時(shí)間與功角的關(guān)系式，使Melnikov函數(shù)由對(duì)時(shí)間的積分變成對(duì)功角的積分形式，再用復(fù)化Simpson公式求得閾值。該方法避免了求解無(wú)擾系統(tǒng)的同宿軌道參數(shù)解析式，并且，無(wú)需將系統(tǒng)輸入的機(jī)械功率設(shè)定為小量。雖然不能得到Melnikov函數(shù)的解析式，但閾值曲線顯示出了可能產(chǎn)生的混沌區(qū)域。

電力系統(tǒng)； 混沌閾值； 同宿軌道； 數(shù)值積分法

電力系統(tǒng)是典型的非線性動(dòng)力系統(tǒng)，存在著豐富的動(dòng)力學(xué)行為，如分岔、混沌和分形等[1，2]，它們對(duì)電力系統(tǒng)安全運(yùn)行有著重要的影響。在這些現(xiàn)象中，混沌是研究最多的之一，它是電力系統(tǒng)低頻振蕩研究對(duì)象的自然延伸，人們發(fā)現(xiàn)電力系統(tǒng)可以通過(guò)倍周期分岔、準(zhǔn)周期分岔、環(huán)面分岔和周期3等途徑進(jìn)入混沌，采用的分析方法主要有Melnikov函數(shù)、Poincaré映射、Lyapunov指數(shù)和頻譜分析等。

Melnikov函數(shù)是研究混沌現(xiàn)象的解析方法，它給出了一類非線性動(dòng)力系統(tǒng)Smale馬蹄變化意義下出現(xiàn)混沌現(xiàn)象的判據(jù)[3]。解析的Melnikov函數(shù)是研究Hamilton系統(tǒng)在弱周期策動(dòng)力激勵(lì)下混沌運(yùn)動(dòng)的最實(shí)用和簡(jiǎn)便的方法，對(duì)于自治可積系統(tǒng)周期擾動(dòng)的混沌性質(zhì)的研究是少有的精確方法之一[4]。

單機(jī)無(wú)窮大電力系統(tǒng)在周期性負(fù)荷擾動(dòng)下的數(shù)學(xué)模型是研究混沌現(xiàn)象的主要模型之一[5，6]，該模型看似簡(jiǎn)單，但它的軌道結(jié)構(gòu)卻較為復(fù)雜。在應(yīng)用Melnikov函數(shù)分析該模型時(shí)，將其視為弱周期擾動(dòng)的Hamilton系統(tǒng)，由于非線性項(xiàng)中含有常數(shù)項(xiàng)，無(wú)法求得對(duì)應(yīng)無(wú)擾系統(tǒng)的同宿軌道參數(shù)表達(dá)式，進(jìn)而無(wú)法求得Melnikov函數(shù)的表達(dá)式，不得不將系統(tǒng)輸入的機(jī)械功率設(shè)定為小量，使得無(wú)擾系統(tǒng)的軌道變?yōu)楫愃捃壍?，影響了分析的?zhǔn)確性，所以，存在著較大的局限性[7]。

Melnikov函數(shù)的數(shù)值積分法可以較好地解決上述問(wèn)題，同時(shí)也可以簡(jiǎn)化計(jì)算過(guò)程。該方法的核心是把對(duì)時(shí)間的積分變換成對(duì)狀態(tài)變量的積分，再用數(shù)值積分法得到可能出現(xiàn)的混沌閾值[8～11]。本文應(yīng)用Melnikov函數(shù)的數(shù)值積分法得到在周期性負(fù)荷擾動(dòng)下電力系統(tǒng)的混沌閾值，避免了求解同宿軌道的參數(shù)表達(dá)式和Melnikov積分的解析式，無(wú)需將系統(tǒng)輸入的機(jī)械功率設(shè)定為小量，從而，較好地反映了電力系統(tǒng)的實(shí)際情況。雖然這種數(shù)值積分法不能得到Melnikov函數(shù)的解析表達(dá)式，但由閾值曲線能清楚地顯示出可能產(chǎn)生混沌的區(qū)域。

1 數(shù)學(xué)模型

經(jīng)過(guò)時(shí)間變換后在周期性負(fù)荷擾動(dòng)下單機(jī)無(wú)窮大電力系統(tǒng)的數(shù)學(xué)模型為

(1)

式中：δ為功角；ω為角速度增量；Pm為輸入的機(jī)械功率；D為阻尼系數(shù)；Pd、Ω分別為擾動(dòng)功率幅值和擾動(dòng)頻率，這些物理量都以標(biāo)幺值表示。

(2)

系統(tǒng)(2)是Hamilton系統(tǒng)，其Hamilton函數(shù)為

(3)

取Pm=0.7時(shí)系統(tǒng)(2)的相平面如圖1所示，圖中(δc，0)為中心點(diǎn)，(δs，0)為鞍點(diǎn)，通過(guò)鞍點(diǎn)的曲線為同宿軌道為

{q(t)|t∈R}={δ(t)，ω±(t)}

(4)

(δp，0)是軌道的轉(zhuǎn)折點(diǎn)，該點(diǎn)可通過(guò)式(3)計(jì)算得到，當(dāng)t0=0時(shí)，δ0=δp。當(dāng)Pm=0時(shí)，無(wú)擾系統(tǒng)的軌道為異宿軌道，本文研究Pm>0的情況。

圖1 無(wú)擾系統(tǒng)(2)相平面(Pm=0.7)

2 Melnikov函數(shù)

無(wú)擾系統(tǒng)(2)的同宿軌道的Melnikov函數(shù)為

(5)

其中

由于ω±(t)取決于Pm，所以廣義積分A±、B±也與Pm有關(guān)。

(6)

則可能出現(xiàn)Smale馬蹄變換意義下的混沌。

3 Melnikov函數(shù)的數(shù)值積分法

如果要得到式(6)的閾值，首先必須求出A±(Pm)、B±(Pm，Ω)的解析式，進(jìn)而求出M±(t0)的解析式。但是，在多數(shù)情況下，A±(Pm)、B±(Pm，Ω)的解析式難以求得，這是因?yàn)闊o(wú)法求得同宿軌道的解析式，即便能夠得到，也未必能求出A±(Pm)、B±(Pm，Ω)的解析式，因此，Melnikov函數(shù)的數(shù)值積分法應(yīng)運(yùn)而生。本文采用文獻(xiàn)[8]的方法，其在文獻(xiàn)[10，11]中得到了較好的應(yīng)用，分別用于同宿和異宿軌道的分析。

由式(2)的對(duì)稱性可知，δ(t)、ω(t)分別是時(shí)間t的偶函數(shù)和奇函數(shù)。對(duì)于式(3)，當(dāng)t>0時(shí)，在同宿軌道上有

(7)

分離上式中的變量，并對(duì)等式兩邊積分

(8)

將式(8)代入A±(Pm)、B±(Pm，Ω)，得

(9)

(10)

應(yīng)用復(fù)化Simpson公式對(duì)式(9)、式(10)進(jìn)行數(shù)值積分，得到A±(Pm)、B±(Pm，Ω)。在對(duì)B±(Pm，Ω)數(shù)值積分時(shí)，由于積分上下限δp、δs是被積函數(shù)的奇異點(diǎn)，則要對(duì)它們進(jìn)行適當(dāng)修正，即上下限為δp+Δδ、δs-Δδ，Δδ為足夠小的正數(shù)，本文采用Δδ=1×10-5[8，10]。圖2是Pd/D-Ω的關(guān)系曲線，即混沌閾值曲線，曲線上方為混沌可能存在的區(qū)域。雖然這種數(shù)值積分法不能得到Melnikov函數(shù)的解析表達(dá)式，但由閾值曲線圖卻能清楚地顯示出可能產(chǎn)生混沌的區(qū)域。只要給定Ω，就可以得到唯一的可能使系統(tǒng)出現(xiàn)混沌的閾值。

圖2 系統(tǒng)(1)的混沌閾值曲線

4 數(shù)值仿真

對(duì)于系統(tǒng)(1)取Pm=0.7、D=0.02、Ω=0.4，代入式(6)得|A±/B±|=0.717 2，則當(dāng)Pd>D|A±/B±|=0.014 3時(shí)，系統(tǒng)可能出現(xiàn)混沌。仿真時(shí)初始狀態(tài)取(arcsinPm，0)，Pd取不同數(shù)值時(shí)得到的略去非穩(wěn)態(tài)后的功角曲線δ(t)和相平面δ-ω如圖3～圖5所示。

(a) 功角曲線δ(t)

(b) 相平面δ-ω

通過(guò)仿真發(fā)現(xiàn)，隨著Pd增大系統(tǒng)出現(xiàn)倍周期分岔，軌道經(jīng)歷了周期1→周期2(Pd=0.06)→周期4(Pd=0.294)的變化，當(dāng)Pd=0.304時(shí)系統(tǒng)失去穩(wěn)定，功角δ(t)趨于無(wú)窮，混沌僅存在于Pd一個(gè)很窄的數(shù)值范圍內(nèi)。倍周期分岔是通向混沌的途徑之一，即由穩(wěn)定的不動(dòng)點(diǎn)→周期1→周期2→周期4→周期8→……無(wú)限倍周期→混沌狀態(tài)[12]。此外，從純數(shù)學(xué)的角度來(lái)說(shuō)δ(t)可以趨于無(wú)窮，但在物理上是不可實(shí)現(xiàn)的，實(shí)際出現(xiàn)的就是混沌運(yùn)動(dòng)[13]。

需要說(shuō)明的是，應(yīng)用Melnikov函數(shù)判定存在Smale馬蹄變換意義下混沌的參數(shù)值與數(shù)值仿真出現(xiàn)混沌的參數(shù)值存在一定誤差，仿真出現(xiàn)混沌的Pd/D值比理論計(jì)算的要大，Pd的量級(jí)比D的量級(jí)要高。實(shí)際上，Melnikov函數(shù)確定的混沌閾值僅是一階近似解，過(guò)大地估計(jì)了混沌的區(qū)域，主要原因是[3]：產(chǎn)生真實(shí)混沌運(yùn)動(dòng)僅有Poincaré映射的第一層次的周期解是不夠的，還需要更高層次的周期解，所有這些層次周期解在構(gòu)成反映真實(shí)混沌運(yùn)動(dòng)的奇異吸引子時(shí)都有不可忽略的作用；Melnikov函數(shù)的自身缺陷造成了Smale馬蹄意義下的混沌并非真正物理意義下的混沌，因而與數(shù)值仿真下的混沌存在差異。但是，無(wú)容置疑，Melnikov函數(shù)已從定性上很好地解釋了從周期軌道到混沌態(tài)的機(jī)理[4]。

(a) 功角曲線δ(t)

(b) 相平面δ-ω

(a) 功角曲線δ(t)

(b) 相平面δ-ω

5 結(jié)語(yǔ)

電力系統(tǒng)混沌現(xiàn)象是電力系統(tǒng)低頻振蕩研究對(duì)象的自然延伸，本文應(yīng)用Melnikov函數(shù)的數(shù)值積分法計(jì)算在周期性負(fù)荷擾動(dòng)下電力系統(tǒng)的混沌閾值，避免了求解同宿軌道的參數(shù)表達(dá)式和Melnikov函數(shù)的解析式，無(wú)需再將輸入的機(jī)械功率設(shè)定為小量，較之以往的方法有了改進(jìn)。雖然這種數(shù)值積分法不能得到Melnikov函數(shù)的解析表達(dá)式，但由閾值曲線能清楚地顯示出可能產(chǎn)生混沌的區(qū)域，只要給定擾動(dòng)頻率，就可以得到唯一的可能使系統(tǒng)出現(xiàn)混沌的閾值。該方法可以擴(kuò)展到對(duì)在準(zhǔn)周期擾動(dòng)下的電力系統(tǒng)混沌研究[14]。

由于無(wú)法得到所研究模型的同宿軌道參數(shù)表達(dá)式，所以難以對(duì)應(yīng)用Melnikov函數(shù)的解析法和數(shù)值積分法進(jìn)行比較。然而，文獻(xiàn)[11]將二者應(yīng)用于同一異宿軌道(該異宿軌道可用解析式表示)的分析和比較，得到的閾值曲線基本一致。因此，只要數(shù)值積分的方法選擇得當(dāng)，Melnikov函數(shù)的數(shù)值積分法即為一種簡(jiǎn)單和有效的方法，可以大大地減少解析推導(dǎo)的工作量。

[1] 王寶華，楊成梧，張強(qiáng)(Wang Baohua，Yang Chengwu，Zhang Qiang).電力系統(tǒng)分岔與混沌研究綜述(Summary of bifurcation and chaos research in electric power system)[J].電工技術(shù)學(xué)報(bào)(Transactions of China Electrotechnical Society)，2005，20(7)：1-10.

[2] 張強(qiáng)，王寶華，楊成梧(Zhang Qiang，Wang Baohua，Yang Chengwu).電力系統(tǒng)安全盆的分形侵蝕及其控制(Fractal erosion of safe basins in power system and its control)[J].電網(wǎng)技術(shù)(Power System Technology)，2005，29(24)：63-67.

[3] 劉曾榮.混沌的微擾判據(jù)[M].上海：上?？萍冀逃霭嫔?，1994.

[4] 李月，楊寶俊.混沌振子系統(tǒng)(L-Y)與檢測(cè)[M].北京：科學(xué)出版社，2007.

[5] 張強(qiáng)(Zhang Qiang).電力系統(tǒng)非線性振蕩研究(Study on nonlinear oscillations of electric power system )[J].電力自動(dòng)化設(shè)備(Electric Power Automation Equipment)，2002，22(5)：17-19.

[6] 朱志宇，蔡立勇，劉維亭(Zhu Zhiyu，Cai Liyong， Liu Weiting). 基于Melnikov方法的電力系統(tǒng)混沌振蕩參數(shù)計(jì)算(Computation of chaotic oscillation parameter in electrical power system based on Melnikov method)[J]. 電力系統(tǒng)及其自動(dòng)化學(xué)報(bào)(Proceedings of the CSU-EPSA)，2008，20(3)：41-45.

[7] 柳明，吳捷(Liu Ming，Wu Jie).微擾電力系統(tǒng)中的次諧及混沌軌道(Chaos and sub-harmonic orbit in power system under perturbation)[J].電力系統(tǒng)自動(dòng)化(Automation of Electric Power Systems)，2002，26(15)：9-14.

[8] Yagasaki K. Chaos in a pendulum with feedback control[J]．Nonlinear Dynamics，1994，6(2)：125-142.

[9] Xu Jianxin，Yan Rui，Zhang Weinian. An algorithm for Melnikov functions and application to a chaotic rotor[J].SIAM Journal on Scientific Computing，2005，26(5)：1525-1546.

[10]李亞峻，李月(Li Yajun，Li Yue).用Melnikov函數(shù)的數(shù)值積分法估計(jì)混沌閾值(Estimate of chaotic threshold by numerical integral method of Melnikov function)[J].系統(tǒng)仿真學(xué)報(bào)(Journal of System Simulation)，2004，16(12)：2692-2695.

[11]Zhuang D，Yu F，Lin Y. Chaotic threshold analysis of nonlinear vehicle suspension by using a numerical integral method[J].International Journal of Automotive Technology，2007，8(1)：33-38.

[12](德)海因茨·奧托·佩特根，哈特穆特·于爾根斯，迪特馬爾·紹柏.混沌與分形——科學(xué)的新疆界[M].2版.田逢春譯.北京：國(guó)防工業(yè)出版社，2008.

[13]裴欽元，李驪(Pei Qinyuan，Li Li).含二次非線性項(xiàng)的受迫振子在主共振曲線上表現(xiàn)的渾沌特性(Chaotic behavior of forced oscillator containing a square nonlinear term on principal resonance curves)[J].應(yīng)用數(shù)學(xué)和力學(xué)(Applied Mathematics and Mechanics)，1995，16(3)：217-223.

[14]張強(qiáng)，王寶華，楊成梧(Zhang Qiang，Wang Baohua，Yang Chengwu).基于二階平均法和Melnikov法準(zhǔn)周期負(fù)荷擾動(dòng)電力系統(tǒng)混沌振蕩分析(Analysis of chaotic oscillation in power system under quasi-periodical load disturbance by using second order averaging method and Melnikov method)[J].電工技術(shù)學(xué)報(bào)(Transactions of China Electrotechnical Society)，2006，21(6)：115-121.

ComputtionofChaoticThreshholdinPowerSystembyNumericalIntegralMethod

ZHANG Qiang， HUANG Xiao-ning

(Electric Power Engineering School， Nanjing Institute of Technology， Nanjing 211167， China)

Melnikov function is the most effective method for homoclinic (heteroclinic) orbits，and a numerical method of the function is applied to compute the chaotic threshold in a single machine infinite bus system disturbed by a periodical load. The relationship between time and power angle is formulated based on the Hamilton function of the undisturbed system and an integral of time variable is changed into one of power angle in the function. Finally, the threshold is reached by the compound Simpson rule. The presented method can avoid the parametric formulation of homoclinic orbit of an undisturbed system, and doesn't need set the mechanical input of the system to be a small quantity. Although the analytical Melnikov function cannot be obtained, a possible chaotic area is shown by a chaotic threshold curve.

power system； chaotic threshold； homoclinic orbit； numerical integral method

2010-08-18；

2010-09-08

江蘇省高校自然科學(xué)基金資助項(xiàng)目(09KJB470002)

TM 712

A

1003-8930(2011)05-0035-04

張 強(qiáng)(1959-)，男，教授，研究方向?yàn)殡娏ο到y(tǒng)運(yùn)行與控制。Email:zhqnj0511@126.com 黃宵寧(1972-)，男，副教授，研究方向?yàn)橛?jì)算機(jī)視覺(jué)及虛擬仿真在電力系統(tǒng)中應(yīng)用。Email:hun_njit@163.com