第27卷第3期 2O10年O6月 工 程 数 学 学 报 v01.27 No.3 June 2010 CHINESE JOURNAL OF ENGINEERING MATHEMATICS 文章编 ̄-:1005—3085(2010)03—0479—08 非线性微分动力系统稳定域计算的波形松弛方法术 蔺小林 (西安交通大学理学院,西安710049) 摘要:非线性微分动力系统稳定域计算是在许多领域具有实际应用的问题。本文对非线性微分动力系统 稳定域的计算方法进行了总结,通过对稳定域边界流形的全面分析,提出了用波形松弛方法计算 稳定域边界流形的思想,给出了计算稳定域边界流形的波形松弛算法。第一步,找出微分动力系 统的所有平衡点,确定渐近稳定平衡点稳定域边界上的平衡点;第二步,用微分动力系统的反方 向系统确定原系统渐近稳定平衡点稳定域边界上平衡点的稳定流形;第三步,渐近稳定平衡点稳 定域的边界是由边界上平衡点和该平衡点稳定流形的并集构成。最后用例子来说明波形松弛方法 在计算稳定域边界流形的有效性。 关键词:非线性微分动力系统:稳定域;波形松弛方法 分类号:AMS(20001 65L07 中图分类号:O175.14 文献标识码:A 1 引言 非线性微分动力系统平衡点的稳定域f或称为吸引域)在诸如电力系统,经济系统,生 物系统等方面都有非常重要的应用。从应用的观点来看希望得到最大的稳定域,这样在 电力系统的设计【 】及相关参数的选取中可以节省由于不必要的过保护设计而产生的浪费。 同样在经济系统、化学反应器、生物系统f2】等方面可以节省大量的人力物力,实现最大 的经济效益。确定稳定域的方法有Lyapunov函数方法、拓扑动力系统方法和其它数值方 法等,此方面比较早的综述文章有Genesio,Tartaglia,Vicino[3J和Varaiya,Wu,Chen[4】, 他it"]总结了二十世纪八十年代以前稳定域的研究工作,给出了计算稳定域的方法和所 取得的成果,同时基于轨线的拓扑结构提出了“轨线逆转”方法。Lyapunov函数方法就 是通过对所研究的微分系统构造相应的Lyapunov函数,应用稳定性方面的结论得到系统 渐近稳定平衡点的稳定域,这种稳定域具有较大的保守性,通常是稳定域的子域,如 通过对Lyapunov函数进行逐次迭代来得到稳定域的一部分以及通过构造最大Lyapunov函 数方法[ ]来得到稳定域等。最近,Chesi[6,7]提出了用一系列连续Lyapunov函数而不是一 个Lyapunov函数来估计多项式微分系统的稳定域。拓扑动力系统方法主要是基于微分拓 扑和动力系统的基本理论,在相空间对轨线的全局结构进行分析,进而给出稳定域的计算 方法,Chiang,Hirsch和Wu[8】,Chiang和Thorp[9j,Lee和Chiang[10】以及Zaborszky,Huang, Zheng和Leung[11]等从拓扑动力学的观点出发对稳定域进行了深入研究,取得了显著的结果。 计算稳定域还有其他方法,如能量函数方法、数值方法等。 波形松弛方法是微分方程数值解法,它将整个微分系统分成若干个子系统,各自独立进行 仿真计算,子系统之间的耦合关系通过子系统中的一些波形信号的相互传递来实现,在系统之 收稿日期:2008—03—07.作者简介:蔺小林(1961年6月生),男,教授.研究方向:微分动力系统稳定域计算 基金项目:国家自然科学基金(10771168). 工 程 数 学 学 报 第27卷 间联系相对较弱的情况下,比较容易得到收敛结果。 本文对微分动力系统的稳定域进行了综合分析,得到微分动力系统一个渐近稳定平衡点稳 定域的边界是由闭轨线或由不稳定平衡点以及该不稳定平衡点的稳定或不稳定流形构成。基于 这些结果给出了寻找稳定域边界流形的波形松弛数值算法,并用例子说明波形松弛方法在计算 稳定域边界流形的有效性。 2 动力系统基本概念与波形松驰算法 对于非线性自治微分动力系统 警: )出 、 , :一 (1) (2) ), 其中X∈M c R ,f:M—M,我们假设向量场.厂是连续函数并且对其变量具有一阶连续偏 导数,以保证系统(1)和(2)解的存在唯一性,在相同时间变化情况下,系统(1)和(2)的轨线 走向方向相反,因此称系统(2)为系统(1)的反方向系统。 方程 f(x)=0 (3) 的一个解称为系统(1)的一个平衡点,我们用集合E={ :.厂 )=0}来表示系统(1)的平衡点 集合,E中既有孤立点也可能含有连续可微闭轨线。 在局部坐标系下,如果在平衡点X 处的雅可比矩阵 .,的特征值没有零实部,则称平 衡点 为系统(1)的双曲平衡点。对一个双曲平衡点X ,可以把系统(1)在点 处的切空 间 (M)唯一地分解成两个子空间 和E“的直和,即E 0 E“,每个子空间在线性算 子 ,下是不变的。对于E 来说, ,的特征值具有负实部;对于E 来说,厶 ,的特征值 具有正实部。记Es的维数为n ,记E“的维数为礼 子空间Es和Eu可表示如下 稳定的子空间 E =span{v ,V ,…,vna), 不稳定的子空间E“=span{w ,W ,.一,W u) 这里V ,u。,…,V 是佗 个(正交化)特征向量,它们的特征值具有负实部, 是n 个(正交化)特征向量,它们的特征值具有正实部,显然n +n =n。 对系统(1),称平衡点x 为n 型。一个零型的平衡点称为汇,一个n型的平衡点称为源, 其它型的平衡点称为鞍点,显然汇是渐近稳定的平衡点,源和鞍点是不稳定的平衡点。 设圣是一个双曲平衡点,它的稳定流形 。(圣)和不稳定流形 “( )定义如下 ( )={ E M:Ot(x)一圣,当t一+。。), (4) “(岔)=( E M:西t( )一圣,当t一一。。). 类似地,一个双曲型闭轨 的稳定流形 (7)和不稳定流形 “ )定义如下 (7)={ E M:Ot(x)一7,当t一+。o), “(7)={ ∈M:西t( )一,y,当t一一。。}. (5) (6) (7) 第3期 蔺小林:非线性微分动力系统稳定域计算的波形松弛方法 481 由于流圣£(.)临界元素的稳定流形对应于 一c(.)临界元素的不稳定流形,通过这种性质可 以把稳定流形的性质转换为不稳定流形的性质。对于一个集Sc ,如果从 中出发的 微分动力系统(1)的每一条轨线对任何时刻t始终属于S,则称 为(1)的不变集。显然集 合 s(.)和 “(.)都是(1)的不变集。 轨线的无穷远性质可以用它的 一极限集和 一极限集来研究。如果存在 中的时问序 列{£i,,当t 一+。。时,使得 Y=.1iratt—+o。 . tt( ), 则称Y属于 的 一极限集(记为 ( ))。类似地,Q一极限集Q( )由tt一一。。来定义。可以证明 这些极限集是集合M的闭不变子集。例如,一个平衡点是它自己的u一极限集,轨线的u一极 限集在它的稳定流形内,轨线的a一极限集属于它的不稳定流形。一个闭轨 既是 上任何点 的 极限集,也是 上任何点的OL一极限集。 横截性的概念是动力系统的基本概念,如果4,B是 中的单射浸入流形,且满足在每一 点 ∈AnB处, 和B的切空间张成z处M的切空间,即 (A)+ (B)= (M), ∈A n B, 或 , 完全不相交,则称 和B满足横截性条件。 个双曲平衡点 的重要特性之一就是它的稳定和不稳定流形在 处横截相交。横截相交 ~的重要性是由于在向量场的扰动下它保持不变。 对系统(1)和系统(2)作如下假设: A11稳定边界上的所有平衡点是双曲型; A21在稳定边界上所有平衡点的稳定和不稳定流形满足横截相交条件; A31稳定边界上每条轨线当t一。o时趋于平衡点之一。 考虑微分方程的初值问题 ), (8) I x(to)=X0,t∈[to, , 其中,:[to,T]× 一R , ∈ ,假设 满足解的存在唯一性条件。 微分动力系统(8)一般的波形松弛迭代格式为 1)( )){ (91 I ( + )(to)=X0,t∈[to, ], 此处函数F满足 F:[to,T]×Ⅱ ×Ⅱ 一Ⅱ “, F(t, , ):f(t, ). 用波形松弛方法来计算微分动力系统初值问题已有许多较好的结论,Jiang[12,13】和Gander, Ruehli[141等用波形松弛方法研究了微分动力系统初值问题,给出了迭代格式收敛的充分性条 件,用波形松弛方法、单调波形松弛方法和Krylov子空间波形松弛方法等对微分动力系统初 值问题以及周期边值问题进行了研究,得到了迭代序列收敛及单调收敛的充分性条件,这些充 分性条件保证本文用波形松弛方法计算微分动力系统一个渐近稳定平衡点稳定域边界流形算法 的收敛性。 482 工 程 数 学 学 报 第27卷 3 确定稳定域边界流形的波形松弛算法 根据渐近稳定平衡点稳定域边界流形的结构特征,在假设条件A1)到A3)成立情况下,有 如下计算稳定域边界的波形松弛算法。 算法确定渐近稳定平衡点X 的稳定域边界OA(x )。 第1步:找出系统(1)的所有平衡点,确定渐近稳定平衡点X。稳定域边界上的平衡点; 第2步:用系统(2)确定系统(1)渐近稳定平衡点 。稳定域边界上这些平衡点的稳定流形: 第3步:渐近稳定平衡点X 稳定域的边界是由第2步确定的平衡点和该平衡点稳定流形的 并集构成。 算法中第1步包括找出f(z)=0的所有解。第2步可以用波形松弛方法来完成,过程如 下: 首先,确定系统(1)渐近稳定平衡点X 稳定域边界上的平衡点。 i)找出系统(1)的所有平衡点,零型平衡点记为X ,否则记为岔 ,并求出系统(1)在平衡 点圣t处的雅可比矩阵; ii)找出雅可比矩阵具有正实部特征值且具有单位长度的不稳定特征向量,记为Yi; iii)找出这些具有单位长度的不稳定特征向量 和平衡点岔 的E一球边界的每一个交截(交 截点是 +eyi和 i一£ ); iv)从每个交截点开始对系统(2)用波形松弛方法计算积分轨线到某一步,如果轨线 保留在这个£.球内,则进行下一步;否则,用 E代替E,相应的交截点为幺士ol ̄yi,此 处0<ol<1,重复这一步; v)从这些交截点开始用波形松弛方法计算积分轨线: vi)重复iji)到v)步。如果任何轨线趋于X ,则平衡点在稳定域边界上。 其次,确定系统(1)渐近稳定平衡点 稳定域边界上平衡点的稳定流形。 对平面系统,稳定域边界上平衡点的类型是一个鞍点或两个源。在这种情型下,1.型平衡 点的稳定流形是一维,它可以用波形松弛方法确定如下: a)找出在平衡点圣处雅可比矩阵单位长度的稳定特征向量Y; b)找出这个单位长度的稳定特征向量Y与平衡点岔的 一球边界的交截(交截点为仝+ 和仝一cy); C)从这些交截点开始对系统(2)用波形松弛方法计算向量场的积分轨线到某一步。如果轨 线保留在这个£一球内,则进行下一步;相反,我们用OlE代替£,相应的交截点为 士oily,此 处0<Ol<1,重复这一步; d)从这些交截点开始对系统(2)用波形松弛方法计算积分轨线; e)所得轨线是系统(2)平衡点的不稳定流形。 对高维系统,计算过程类似于上述过程,计算结果能提供渐近稳定平衡点稳定域边界的稳 定流形轨线集。 注如果系统(1)只有一个渐近稳定平衡点,则在该平衡点较远处取几个点,从这些点处开 始用波形松弛方法计算积分轨线,可以得到渐近稳定平衡点稳定域边界。 4 数值算例 例1考虑如下系统【8】 第3期 蔺小林:非线性微分动力系统稳定域计算的波形松弛方法 dXl = ’ dXl ’ (10) (11) 容易得到系统(10)的两个平衡点(0,0)和(1,2),其中(0,0)是渐近稳定平衡点,而(1,2)是1-型 平衡点,其部分不稳定流形收敛到平衡点(0,0)。对系统(11)在点(1,2)处,由于 Of 0,2)=(2 )j(1j ( 特征值A1= 所对应的单位特征向量为Q1=(一1/3,v ̄/3)T,特征值A2=一 所对应的单 位特征向量为OL2:(1/3, /3) ,取e=0.1,则在点 ( (0 , )=(0.966667,2.047133)和(面(0 , ’)=(1.033333,1.952867) 处用波形松弛方法对系统(ii)进行求解,其雅可比迭代格式为 f 一 卅 ,, l : 一 ,七 1)2,…. 画出轨线图如图1所示,其中虚线为用波形松弛方法计算出的系统(10)渐近稳定平衡 点(0,0)的稳定域边界流形,所得稳定域与系统实际稳定域相同。 一 ———— \\、\'}——————一 —/。 图1: 系统(10)渐近稳定平衡点(0,0)的稳定域 例2考虑如下系统f3,sl I, 一卅 I dx2:0 一2m 一 --O.1x13, (12) 工 程 数 学 学 报 第27卷 (1x1 _dt 二1一 2, l+2.。 2+ +0.1x3・ (13) dx2=一 系统有三个平衡点:(0,0)是系统(12)的局部渐近稳定平衡点,(一2.55,。2.55)是系统(12)的 1一型平衡点,(一7.45,一7.45)是系统(12)的局部渐近稳定平衡点。对于系统(13)在点(-2・55,-2.55) 处,由于 、 , \ of l 一 =. ,一 .(一。. +2 1 +。.3 2-.。1)j 一。. 。,一 . (一3. 925-21), 特征值 1:3.370628所对应的特征向量为 1:(1,一2.370628) ,特征值 2=一0・370628所对 应的特征向量为a2=(1,0.629372)T,取E=0.1,则在点 ), :(-2.450000,一2.823063)和(童(0 , )=(-2.650000,一2.312937) 处用波形松弛方法对系统(13)进行求解,其雅可比迭代格式为 f : i ,, l dx(k+1)-0.1x ̄州)+2m +(z +0.1( , 1j2,…・ 画出轨线图如图2所示,其中虚线为用波形松弛方法计算出的系统(12)渐近稳定平衡点 的(0,0)和(一7.45, 7.45)的稳定域边界流形,所得稳定域与系统实际稳定域相同。 y 图2: 系统02)渐近稳定平衡点(0,0)和(一7.45,一7.45)的稳定域 例3考虑Van der pol振子方程[。 d 1 —d—t —一 。, (14) dx2———d t : 1+ (+A【 一J、 i一1)X2,一 , (15) dx2-面Xl--) ̄(z _1) 第3期 蔺小林:非线性微分动力系统稳定域计算的波形松弛方法 系统仅有一个平衡点(0,0),取参数 =1,在(0,0)较远处取一点(x7 , ),比如(-1.8,0.3), 从此点开始对系统(15)用波形松弛方法求解,其雅可比迭代格式为 f dx ̄k +1): , 【 dx(k+1):一 ; )+ +¨一( i ) , :。,1,2,…. 图3: =l时系统(14)渐近稳定平衡点(0,0)的稳定域 5 小结 我们用波形松弛方法计算了微分动力系统渐近稳定平衡点稳定域的边界流形,得到了渐近 稳定平衡点的稳定域,由于已经得到微分代数系统多指标【15,16,18]以及泛函微分方程【l7】波形松 弛算法的收敛性结论,我们可以尝试用波形松弛方法计算微分代数系统(含多指标)以及泛函微 分方程的稳定域或者收敛域问题。 参考文献: [1】Praprost K L,Loparo K A.A stability theory for constrained dynamic systems with applications to electric power systems[J】.IEEE Trans on Automatic Control,1996,41(11):1605—1617 【2】Ratledge C,Kristiansen B.Basic Biotechnology[M].New York:Cambridge,2001 【3】Genesio R,Tartaglia M,Vicino A.On the estimation of symptaotic stability regions:state of the art and new proposals[J】.IEEE Transactions on Automatic Control,1985,ac一3(8):747—755 【4】Varaiya P,Wu F F,Chen R L.Direct method for transient stability analysis of power systems:recent results[J].Proceedings of the IEEE,1985,73:1703-1715 【5】Vannelli A,Vidyasagar M.Maximal Lyapunov functions and domains of attrctiaon for autonomous non— linear systems[J].Automatica,1985,21:69—80 【6】Chesi G.Estimating the domain of attraction for uncertain polynomial systems[J].Automatica,2004,40: 1981—1986 【7】Chesi G.Estimating the domain of attraction via union of continuous families of Lyapunov estimates[J]. Systems Control Letters,2007,56:326 333 工 程 数 学 学 报 第27卷 【8】Chiang H D,Hirsch M,Wu F F.Stability regiOIlS of nonlinear autonomous dynamical systems[J].IEEE Transactions on Automatic Control,1988,33:16—27 【9】Chiang H D,Thorp J S.Stability regions of nonlinear dynamical systems:a constructive methodology[J]. IEEE Transactions on Automatic Control,1989,ac一34:1229—1241 【10】Lee J,Chiang H D.Theory of stability region for a class of no hyperbolic dynamical systems and its appli— cation to constraint satisfctaion problems[J1.IEEE Transctaions on Circuits and Systems I:Fundamental Theory and Applications,2002,49(2):196-209 [11】Zaborszky J,Huang G,Zheng B,et a1.On the phase portrait of a clssa of large nonlinear dynamic systems such as the power system[J].IEEE Transactions on Automatic Control,1988,33:4-15 [12】Jiang Y L,Chen R M M,Wang O.Waveform relaxation of nonlinear second—order diferential equations[J]. IEEE ransTactions on Circuits and Systems I:Fundamental Theory and Applications,2001,48(11):134 1347 【13]Jiang Y L.Periodic waveorfm relxataion solutions of nonlinear dynamic equations[J].IEEE Applied Math- ematics and Computation,2003,135(2—3):219-226 f141 Gander M,Ruehli A E.Optimized waveform relaxation methods orf RC type circuits[J1.IEEE ransTactions on Circuits and Systems I:Fundamental Theory and Applications,2004,51(4):755—768 【15】李宝成,蒋耀林.关于动力学系统伪谱与扰动系统谱之间关系的研究[J].工程数学学报,2004,21(6):936—940 Li B C,Jiang Y L.A research on relationships for pseud0spectra of dynamical systems and spectra of its perturbed systems[JJ.Chinese Journal of Engineering Mathematics,2004,21(6):936—940 『161孙卫,樊晓光,李立.关于微分.代数系统Runge-Kutta迭代方法的研究[J].工程数学学报,2005,22(6):1070. 1075 Sun W,Fan X G,Li L.Research on iterative Runge-Kutta methods for DAEs[J].Chinese Journal of Engineering Mathematics,2005,22(6):1070—1075 [17]王宇莹,赵炜,王峰.一类泛函微分代数系统的波形松弛方法[J].工程数学学报,2006,23(1):71—79 Wang Y Y,Zhao W,Wang F.The waveform relaxation method for functional diferential・algebraic equa- tions[J】.Chinese Journal of Engineering Mathematics,2006,23(1):71—79 [18】孙卫,邹建华,王彤.非线性指标.3微分.代数系统波形松弛算法的实现【J].工程数学学报,2007,24(3):539—542 Sun W,Zou J H,Wang T.Implementing waveform relaxation method for nonliear differential—algebraic systems of index一3[J】.Chinese Journal of Engineering Mathematics,2007,24(3):539—542 The Waveform Relaxation Method for Estimating the Asymptotic Stability Region of Nonlinear Differential Dynamic Systems LIN Xiao-lin (School of Science,Xi’an Jiaotong University,Xi’an 710049) Abstract:Estimating the asymptotic stability region of nonlinear differential dynamic systems is a problem applied in many fields.The paper outlines the computing methods for estimating the asymp- totic stability region of nonlinear differential dynamic systems.By analyzing the boundary manifolds of the asymptotic stability region,we present the idea of calculating the boundary manifolds of asymp ̄ totic stability region with the waveform relaxation method and offer the algorithms.Firstly,find out all equilibrium points and determine which equilibrium point states in the boundary of asymptotic stability region of the stable equilibrium point.Secondly,use the negative system of nonlinear differ— ential dynamic systems to fix stable manifolds of equilibrium point setting in the boundary.Thirdly, the boundary of the asymptotic stability region of stable equilibrium point is composed of stable equi— librium points and their stable manifolds.Numerical examples further illustrate the correctness of the theoretical results in this paper. Keywords:nonlinear differential dynamic systems;region of asymptotic stability;waveform relaxation method Received:07 Mar 2008. Accepted:07 May 2008. Foundation item:The National Natural Science Foundation of China(10771 168)