应用数值积分法计算电力系统混沌阈值

张 强， 黄宵宁

(南京工程学院电力工程学院， 南京 211167)

应用数值积分法计算电力系统混沌阈值

张 强， 黄宵宁

(南京工程学院电力工程学院， 南京 211167)

Melnikov函数是分析同(异)宿轨道出现混沌的最有效方法，用该函数的数值积分法计算单机无穷大电力系统在周期性负荷扰动下的混沌阈值。通过相应无扰系统的Hamilton函数求得时间与功角的关系式，使Melnikov函数由对时间的积分变成对功角的积分形式，再用复化Simpson公式求得阈值。该方法避免了求解无扰系统的同宿轨道参数解析式，并且，无需将系统输入的机械功率设定为小量。虽然不能得到Melnikov函数的解析式，但阈值曲线显示出了可能产生的混沌区域。

电力系统； 混沌阈值； 同宿轨道； 数值积分法

电力系统是典型的非线性动力系统，存在着丰富的动力学行为，如分岔、混沌和分形等[1，2]，它们对电力系统安全运行有着重要的影响。在这些现象中，混沌是研究最多的之一，它是电力系统低频振荡研究对象的自然延伸，人们发现电力系统可以通过倍周期分岔、准周期分岔、环面分岔和周期3等途径进入混沌，采用的分析方法主要有Melnikov函数、Poincaré映射、Lyapunov指数和频谱分析等。

Melnikov函数是研究混沌现象的解析方法，它给出了一类非线性动力系统Smale马蹄变化意义下出现混沌现象的判据[3]。解析的Melnikov函数是研究Hamilton系统在弱周期策动力激励下混沌运动的最实用和简便的方法，对于自治可积系统周期扰动的混沌性质的研究是少有的精确方法之一[4]。

单机无穷大电力系统在周期性负荷扰动下的数学模型是研究混沌现象的主要模型之一[5，6]，该模型看似简单，但它的轨道结构却较为复杂。在应用Melnikov函数分析该模型时，将其视为弱周期扰动的Hamilton系统，由于非线性项中含有常数项，无法求得对应无扰系统的同宿轨道参数表达式，进而无法求得Melnikov函数的表达式，不得不将系统输入的机械功率设定为小量，使得无扰系统的轨道变为异宿轨道，影响了分析的准确性，所以，存在着较大的局限性[7]。

Melnikov函数的数值积分法可以较好地解决上述问题，同时也可以简化计算过程。该方法的核心是把对时间的积分变换成对状态变量的积分，再用数值积分法得到可能出现的混沌阈值[8～11]。本文应用Melnikov函数的数值积分法得到在周期性负荷扰动下电力系统的混沌阈值，避免了求解同宿轨道的参数表达式和Melnikov积分的解析式，无需将系统输入的机械功率设定为小量，从而，较好地反映了电力系统的实际情况。虽然这种数值积分法不能得到Melnikov函数的解析表达式，但由阈值曲线能清楚地显示出可能产生混沌的区域。

1 数学模型

经过时间变换后在周期性负荷扰动下单机无穷大电力系统的数学模型为

(1)

式中：δ为功角；ω为角速度增量；Pm为输入的机械功率；D为阻尼系数；Pd、Ω分别为扰动功率幅值和扰动频率，这些物理量都以标幺值表示。

(2)

系统(2)是Hamilton系统，其Hamilton函数为

(3)

取Pm=0.7时系统(2)的相平面如图1所示，图中(δc，0)为中心点，(δs，0)为鞍点，通过鞍点的曲线为同宿轨道为

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

(4)

(δp，0)是轨道的转折点，该点可通过式(3)计算得到，当t0=0时，δ0=δp。当Pm=0时，无扰系统的轨道为异宿轨道，本文研究Pm>0的情况。

图1 无扰系统(2)相平面(Pm=0.7)

2 Melnikov函数

无扰系统(2)的同宿轨道的Melnikov函数为

(5)

其中

由于ω±(t)取决于Pm，所以广义积分A±、B±也与Pm有关。

(6)

则可能出现Smale马蹄变换意义下的混沌。

3 Melnikov函数的数值积分法

如果要得到式(6)的阈值，首先必须求出A±(Pm)、B±(Pm，Ω)的解析式，进而求出M±(t0)的解析式。但是，在多数情况下，A±(Pm)、B±(Pm，Ω)的解析式难以求得，这是因为无法求得同宿轨道的解析式，即便能够得到，也未必能求出A±(Pm)、B±(Pm，Ω)的解析式，因此，Melnikov函数的数值积分法应运而生。本文采用文献[8]的方法，其在文献[10，11]中得到了较好的应用，分别用于同宿和异宿轨道的分析。

由式(2)的对称性可知，δ(t)、ω(t)分别是时间t的偶函数和奇函数。对于式(3)，当t>0时，在同宿轨道上有

(7)

分离上式中的变量，并对等式两边积分

(8)

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

(9)

(10)

应用复化Simpson公式对式(9)、式(10)进行数值积分，得到A±(Pm)、B±(Pm，Ω)。在对B±(Pm，Ω)数值积分时，由于积分上下限δp、δs是被积函数的奇异点，则要对它们进行适当修正，即上下限为δp+Δδ、δs-Δδ，Δδ为足够小的正数，本文采用Δδ=1×10-5[8，10]。图2是Pd/D-Ω的关系曲线，即混沌阈值曲线，曲线上方为混沌可能存在的区域。虽然这种数值积分法不能得到Melnikov函数的解析表达式，但由阈值曲线图却能清楚地显示出可能产生混沌的区域。只要给定Ω，就可以得到唯一的可能使系统出现混沌的阈值。

图2 系统(1)的混沌阈值曲线

4 数值仿真

对于系统(1)取Pm=0.7、D=0.02、Ω=0.4，代入式(6)得|A±/B±|=0.717 2，则当Pd>D|A±/B±|=0.014 3时，系统可能出现混沌。仿真时初始状态取(arcsinPm，0)，Pd取不同数值时得到的略去非稳态后的功角曲线δ(t)和相平面δ-ω如图3～图5所示。

(a) 功角曲线δ(t)

(b) 相平面δ-ω

通过仿真发现，随着Pd增大系统出现倍周期分岔，轨道经历了周期1→周期2(Pd=0.06)→周期4(Pd=0.294)的变化，当Pd=0.304时系统失去稳定，功角δ(t)趋于无穷，混沌仅存在于Pd一个很窄的数值范围内。倍周期分岔是通向混沌的途径之一，即由稳定的不动点→周期1→周期2→周期4→周期8→……无限倍周期→混沌状态[12]。此外，从纯数学的角度来说δ(t)可以趋于无穷，但在物理上是不可实现的，实际出现的就是混沌运动[13]。

需要说明的是，应用Melnikov函数判定存在Smale马蹄变换意义下混沌的参数值与数值仿真出现混沌的参数值存在一定误差，仿真出现混沌的Pd/D值比理论计算的要大，Pd的量级比D的量级要高。实际上，Melnikov函数确定的混沌阈值仅是一阶近似解，过大地估计了混沌的区域，主要原因是[3]：产生真实混沌运动仅有Poincaré映射的第一层次的周期解是不够的，还需要更高层次的周期解，所有这些层次周期解在构成反映真实混沌运动的奇异吸引子时都有不可忽略的作用；Melnikov函数的自身缺陷造成了Smale马蹄意义下的混沌并非真正物理意义下的混沌，因而与数值仿真下的混沌存在差异。但是，无容置疑，Melnikov函数已从定性上很好地解释了从周期轨道到混沌态的机理[4]。

(a) 功角曲线δ(t)

(b) 相平面δ-ω

(a) 功角曲线δ(t)

(b) 相平面δ-ω

5 结语

电力系统混沌现象是电力系统低频振荡研究对象的自然延伸，本文应用Melnikov函数的数值积分法计算在周期性负荷扰动下电力系统的混沌阈值，避免了求解同宿轨道的参数表达式和Melnikov函数的解析式，无需再将输入的机械功率设定为小量，较之以往的方法有了改进。虽然这种数值积分法不能得到Melnikov函数的解析表达式，但由阈值曲线能清楚地显示出可能产生混沌的区域，只要给定扰动频率，就可以得到唯一的可能使系统出现混沌的阈值。该方法可以扩展到对在准周期扰动下的电力系统混沌研究[14]。

由于无法得到所研究模型的同宿轨道参数表达式，所以难以对应用Melnikov函数的解析法和数值积分法进行比较。然而，文献[11]将二者应用于同一异宿轨道(该异宿轨道可用解析式表示)的分析和比较，得到的阈值曲线基本一致。因此，只要数值积分的方法选择得当，Melnikov函数的数值积分法即为一种简单和有效的方法，可以大大地减少解析推导的工作量。

[1] 王宝华，杨成梧，张强(Wang Baohua，Yang Chengwu，Zhang Qiang).电力系统分岔与混沌研究综述(Summary of bifurcation and chaos research in electric power system)[J].电工技术学报(Transactions of China Electrotechnical Society)，2005，20(7)：1-10.

[2] 张强，王宝华，杨成梧(Zhang Qiang，Wang Baohua，Yang Chengwu).电力系统安全盆的分形侵蚀及其控制(Fractal erosion of safe basins in power system and its control)[J].电网技术(Power System Technology)，2005，29(24)：63-67.

[3] 刘曾荣.混沌的微扰判据[M].上海：上海科技教育出版社，1994.

[4] 李月，杨宝俊.混沌振子系统(L-Y)与检测[M].北京：科学出版社，2007.

[5] 张强(Zhang Qiang).电力系统非线性振荡研究(Study on nonlinear oscillations of electric power system )[J].电力自动化设备(Electric Power Automation Equipment)，2002，22(5)：17-19.

[6] 朱志宇，蔡立勇，刘维亭(Zhu Zhiyu，Cai Liyong， Liu Weiting). 基于Melnikov方法的电力系统混沌振荡参数计算(Computation of chaotic oscillation parameter in electrical power system based on Melnikov method)[J]. 电力系统及其自动化学报(Proceedings of the CSU-EPSA)，2008，20(3)：41-45.

[7] 柳明，吴捷(Liu Ming，Wu Jie).微扰电力系统中的次谐及混沌轨道(Chaos and sub-harmonic orbit in power system under perturbation)[J].电力系统自动化(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函数的数值积分法估计混沌阈值(Estimate of chaotic threshold by numerical integral method of Melnikov function)[J].系统仿真学报(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](德)海因茨·奥托·佩特根，哈特穆特·于尔根斯，迪特马尔·绍柏.混沌与分形——科学的新疆界[M].2版.田逢春译.北京：国防工业出版社，2008.

[13]裴钦元，李骊(Pei Qinyuan，Li Li).含二次非线性项的受迫振子在主共振曲线上表现的浑沌特性(Chaotic behavior of forced oscillator containing a square nonlinear term on principal resonance curves)[J].应用数学和力学(Applied Mathematics and Mechanics)，1995，16(3)：217-223.

[14]张强，王宝华，杨成梧(Zhang Qiang，Wang Baohua，Yang Chengwu).基于二阶平均法和Melnikov法准周期负荷扰动电力系统混沌振荡分析(Analysis of chaotic oscillation in power system under quasi-periodical load disturbance by using second order averaging method and Melnikov method)[J].电工技术学报(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

江苏省高校自然科学基金资助项目(09KJB470002)

TM 712

A

1003-8930(2011)05-0035-04

张 强(1959-)，男，教授，研究方向为电力系统运行与控制。Email:zhqnj0511@126.com 黄宵宁(1972-)，男，副教授，研究方向为计算机视觉及虚拟仿真在电力系统中应用。Email:hun_njit@163.com