动力系统建模

1、动力系统建模唐 云（清华大学数学科学系）引言1。动力系统的基本概念与方法：从蝴蝶泉到蝴蝶效应2。机械与电力系统中的数学模型3。生命科学的数学建模4。分形模型目 录引言 动力系统是关于系统演化规律的数学学科。数学建模是连接动力系统理论与应用的桥梁。 特点：（1）广泛的应用前景；（2）深刻的数学理论基础；（3）计算技术。 动力系统建模作用：（1）教学：数学建模；（2）科学研究：与国民经济及科学前沿关系密切。从蝴蝶泉到蝴蝶效应1.1 蝴蝶的生态问题 云南大理蝴蝶泉，是影片五朵金花里阿鹏和金花对歌的地方，也有名的游览胜地。泉内蝴蝶种类繁多，每年农历4月15日白族的“蝴蝶会”前后，蝴蝶大的大如巴掌，小的

2、小如蜜蜂，成串悬挂于泉边的合欢树上，盛况空前。明代徐霞客在其“游记”中称：“真蝶万千，连须钩足，自树巅倒悬而下及于泉面”。郭沫若在1961年游蝴蝶泉时也曾留下“首尾联接数公尺，自树下垂疑花序”的诗句。 令人惋惜的是，近十数年，人们已经很难看到美丽的蝴蝶盛会，有时，虽有蝴蝶聚集，但数量已少。据当地父老传言：蝴蝶泉边，原有一蓬枝叶茂密、开白花、发清香的茨蓬，花枝缠在横斜泉面的树干上，蝴蝶沿着这些下垂的花枝连成串。如今，茨蓬已除，泉面树干叶枯，加上周围自然环境受到破坏，田野大量使用农药，误伤不少蝴蝶，那连须钩足悬于泉面的奇观，久已不见。1.2 数学模型: Logistic映射).1 ( Logest

3、ic).(1)()()() 1( Verhust . .)( 1 (0).)( ).() 1( Malthus , 1).()() 1( )()() 1()()( 12nnnnxxxnPabnaPnPbnaPnPnnPaPanPnaPnPkankPnPnPnPnPnPnnP模型经尺度变换，得的修正于是有实际这显然不合无限增大随时，当易见模型则得命成正比，即与若增量开始时的总数。比如一个月时间段表示蝴蝶在第设f(x)= ax(1- x)， x 在0,1内变化xn+1= f(xn) 从0,1内点x0出发，由Logistic映射的迭代形成xn= f n(x0), n = 0,1,2,序列xn称为x0

4、的轨道。由轨道迭代的极限集可组成分岔图种群数的模型简化：相应的迭代为了一个序列，即Logistic映射Logistic映射分岔图关于吸引子 吸引子是动力系统相空间中的一类特殊集合，它能把周围的轨道都“吸引”（收敛）过来。 吸引子分类： （1）平衡点 （2）周期轨 （3）拟周期轨 （4）混沌吸引子 当0a 1时，由于 当1a3时，任何(0,1)中初始值的轨道趋于 两个不动点x1*, x2* ,一个稳定(吸引),另一个xn 0 物种逐渐灭亡x*=1-1/a其中x*是方程f(x)=x的解，为映射f 的不动点（周期1点)例：a =1.5时 xn 1/3.不稳定，轨道xn趋向稳定点。数值迭代：倍周期分岔

5、,01nnaxx 当1+61/2a3.5440903506时, 从任意的点x0出 也称为周期2点,对应轨道称周期2轨道.(原来周期点失稳)发的轨道将逐渐沿着四个数值振动，它们满足),(),( ,6132*4*3xfxxfxxxxan它们满足振动，绕着两个数时，当称为周期4点，对应轨道称周期4轨道(原有周期点 若a再增大，周期4点又会失稳，而产生新的稳定又失稳)周期8点，这个周期不断加倍的过程将重复无限次，会依次出现周期16点，周期32点，. ，这种过程称为倍周期分岔.相应的分岔值c1=3, c2=1+61/2构成一个单调增加的数列ck.其极限值为c*=3.569945557391。),(),(

6、),(),(324xfxxfxxfxxfx 当c*a4时，Logistic映射进入混沌区域.反映出 遍历性：点 x0的轨道不趋向任何稳定的周期性,即不同初始值，即使它们离得非常近，它们的的是：轨道, 它的轨道在(0,1)(或其中某些区间)内的任何一个子区间(a,b)内都会出现无数次. 敏感性：轨道表现出对初始条件的强烈敏感轨道也终将以某种方式分离.混沌的特点 Feigenbau常数(ck-ck-1)/(ck+1-ck)在k趋于无穷时，趋于常数 q =4.6692016这常数的意义在于普适性，例如周期3窗口，还有 存在周期窗口：混沌区域内某些地方仍有倍周期分岔，例如a3.835附近其他映射 任取

7、(0,1)中的点x0,可以通过作图来取得迭代 在以xn为横坐标、xn+1为纵坐标的第一象限作抛物线弧：xn+1a xn(1- xn)的数值序列xn,从而也通过图象直观地看出由x0出发的轨道的变化. 这作图的过程颇象蜘蛛织网,故称为蛛网迭代. 图像方法：蛛网迭代11xnxn1x0 x1x1x2 1a3 从(0,1)中任何初值出发的轨道趋向不动点 (周期1点) 3a61/2+1 从任何初值出发的轨道趋向周期2点61/2+1a 3.54409035从任何初值出发的轨道趋向周期4点 a=3.58轨道进入浑沌状态 a= 4 轨道的浑沌性表现充分 蛛网迭代的优点是轨道非常直观形象.缺点是当周期数较大时不易

8、看清轨道变化细节 密度分布图：密度分布图： 从密度：从一个初始点 x0出发，由迭代所 产生的序列xn (n一般很大)在区间 0,1上的概率分布密度. 将具体算法：将0,1区间分成m个长度为h=1/m的小区间，序列xnnN=0 落在各个小区间ih,(i+1)h的个数为ki,则该序列落在各小区间的概率（即密度)为pi=ki/N i=0,1,2,m 密度图：横轴为区间 0,1, 纵轴为概率 p.每个小区间上的细柱线的高度等于该区间上密度 a=3.2 (m=100 N=10000 x0= 0.1)(这是周期2情况) a=3.45(这是周期4情况) a=3.55(周期8的情况) 以上密度图显示在 0ac

9、*的情况下,xn只有极少数落在周期点以外的小区间,而最终以几乎相等的概率落在周期点所在的小区间。 a=3.6(进入混沌区) (最混沌状态) a= 4LogisticLogistic模型的混沌自相似（分形）模型的混沌自相似（分形）图一图一 图二图二 图三：图二局部放大图图三：图二局部放大图图四：图三局部放大图图四：图三局部放大图 1961年冬天年冬天E. Lorenz 进行关于天气预报的计进行关于天气预报的计算。他考虑下面加热的流体由热传导进入对流，算。他考虑下面加热的流体由热传导进入对流，然后产生湍流的过程，然后产生湍流的过程， 对对Rayleigh-Bernard方方程进行约化，得到下面的程

10、进行约化，得到下面的Lorenz方程。方程。D. Gulick, Encounters with Chaos, Mc-Graw Hill, Inc., New York, 1992.()xyxyxyxzzzxy 和和 为正数为正数. , , 和和 与流体的物理性质相关与流体的物理性质相关. Lorenz 取取 , , 和和 ., xyz102883Lorenz 吸引子吸引子对初值条件的敏感依赖性对初值条件的敏感依赖性 10,000 个几乎相同的个几乎相同的初值条件初值条件 从这从这10,000 个初值出个初值出发的每条轨道经过同样发的每条轨道经过同样时间后其终点很不一致时间后其终点很不一致.

11、这些点都在同一个吸引这些点都在同一个吸引子范围内子范围内.S. H. Strogatz, Nonlinear Dynamics and Chaos with Applications to Physics, Biology, Chemistry and Engineering, Addison-Wesley Publishing Company, New York, 1994.所谓“蝴蝶效应”是指： 初始值很小的改变会引起绝然不同的结果.巴西的一只蝴蝶扇动翅膀会引起明年在得克萨斯的大风暴吗？E. Lorenz 数学的伟大使命在于从混沌中发现秩序。 倍尔在那个混沌的体制中，结构上的微小差异几乎都

12、会造成行为方式上的巨大变化，可控制的行为似乎已被排除。 斯图尔特.考夫曼关于混沌 一个动力系统是混沌的，如果它满足： (1) 有一个由周期轨道组成的稠密集合； (2) 轨道敏感地依赖于其初值条件； (3) 为拓扑传递的. 混沌的度量性质： (1)正Lyapunov指数 (2)拓扑熵与测度熵 (3)分形结构混沌出现在各个领域的一种现象:数学、物理、 由此引起的复杂而有趣的现象 “侏罗纪公园”中的恐龙重现生物、金融、经济、管理等等： 宇宙的起源 龙卷风的产生、厄尔尼诺现象 东南亚金融危机爆发 可以从某些简单的离散的数学模型开始, 讨论2. 机械和电力系统的数学模型2.1 动力学模型 Newton力

13、学体系是第一个，也是最基本的动力系统数学模型。建模过程： （1）数据积累：第谷(Tycho Brahe, 1546-1601) （2）经验公式：开普勒(J. Kepler, 1571-1630)的行星“三大定律”。 （3）数学模型：牛顿(I. Newton, 1642-1727)的“万有引力”。 （4）验证：如哈雷彗星和海王星的发现。非线性振动均为向量形式。和多自由度系统：上述单自由度系统： ).,(fxtxxfx Duffing 方程tfxxxxcos23 得，现取 ,31. 01.20.5 ftxxxx2 . 1cos31. 03 . 03 Duffing 方程方程位移位移x位移位移x时间

14、时间Tdx/dt0001. 0) 0(,0001. 2) 0(.0000. 0) 0(,0000. 2) 0( xxxx蓝线红线位移位移x时间时间Tdx/dt位移位移x000001. 0) 0(,000001. 2) 0(.000000. 0) 0(,000000. 2) 0( xxxx蓝线红线 耗散系统相体积的演化 解释初值敏感和奇怪吸引子要用到相体积的伸展与折叠，定量描述这一特征的量是李雅普诺夫指数，这要证明三维（以上）相体积 而 中要有正有负。讲混沌时专门解释或用初等的例子说明不严格，并且要花费一定学时。但只要在前面讲正则方程和刘维定理时稍做改动就可以自然引出耗散系统相体积演化公式，概念

15、的引入严格、定量，学时反而减少。 通常讲刘维定理，对相空间保守体系ttteeeVxxxV3210321321,00dtVddtd或其证明用到过去教材都是代入保守体系的正则方程。实际上，对一力学体系，如果除保守力外还含非保守力，则正则方程应写为其中 代表非保守力，代入前式后得出，积分马上得到)(1ppqqdtddtVdfQqHppHqVpQdtVdf1dtpQfeVV10Q容易证明，对保守力 对非保守力（如 ）对高维耗散系统，自然导出指数形式 形式其中 （i=1f）可正可负总和为正。01fpQkxxxF)(teVV20teVV0f321i极限环与分岔考虑一个数学例子： 2222122221221

16、222221122211211)()(2)()(2xxxxxxcxxxxxxxxxxcxx用二维极坐标表示，作代换.sin,cos21rxrx2242rrcrr得到 当 时有 0r 0)2(42rrcr方程有三个根： ， 00rcr1121cr1122-1.5-1-0.500.51-1-0.500.511.5, 1c00r，点是稳定定态（焦点） 01c（硬激励）00 rr稳定定态（焦点）； 2rr 1rr 不稳定的极限环； 稳定的极限环。1rr 黑洞 0c（软激励） 1rr 0r 1rr 0r 1rr 是稳定的极限环 00 rr不稳定的焦点极限环极限环分岔分岔叉式分岔叉式分岔：定性举例，旋转单

17、摆 Hopf Hopf 分岔分岔：点到极限环的突变 1c稳定定态； 稳定极限环； 0c01c终态稳定点与初始条件有关（亚临界Hopf 分岔）。 倍周期分岔倍周期分岔: : 周期成倍突变 阻尼单摆的强迫振动方程引入无因次量后化为 :tfcossin2 00sin , 0时， 0稳定的定态， ，不稳定点， 0 摆长为摆长为l ，小球质量为，小球质量为m的单摆，相对与平衡的下垂位置的角位的单摆，相对与平衡的下垂位置的角位移为移为，重力加速度为，重力加速度为g ，则其运动方程为：，则其运动方程为： (1)或或 (2) 等式右边是周期性的驱动力，其角频率为。把方程式无量纲化，等式右边是周期性的驱动力，其

18、角频率为。把方程式无量纲化，用去除每一项，将无量纲的时间叫做用去除每一项，将无量纲的时间叫做t ，即得，即得 (3)式中式中都是无量纲化的。都是无量纲化的。 阻尼单摆的强迫振荡阻尼单摆的强迫振荡tmlFmDcossin20 tfcossin2 tFmglmlDcossin)(sintmglI )(2mod)sin(1nKnnCK旋转数mnnmmK01lim21),(其中)(KCHenon-Heiles 星体势模型星体势模型322231)(21yyxyxV32222231)(21)(21yyxyxyxH等势图等势图粒子真实轨迹粒子真实轨迹粒子的庞加莱截面粒子的庞加莱截面)(cos22nTtIJH

19、).(sinnTtJIJ激励转子激励转子2.3 关于电力系统的数学模型 电力工业是国民经济的支柱产业，并且直接影响到 人民的生活。 我国电力工业发展迅速：19491990 2001装机容量（亿千瓦）0.0181.353.36年发电量（万亿千瓦时） 0.00430.6181.45世界排名2542重大停电事故时间地点造成损失1978年12月法国电网电压崩溃停电4-7小时，直接损失2亿美元1982年12月加拿大魁北克停电8.5小时1987年7月日本东京停电3小时1996年7月美国西部200万用户停电3小时1996年8月美国西部750万用户停电3-6小时北美东部大停电 2003年8月14日下午3时许，

20、俄亥俄州北部34.5万伏超高压突然发生故障。一小时内波及包括纽约和多伦多在内的美加东中部大停电，5000万人陷入黑暗之中。至16日10时基本恢复正常。 估计经济损失每天达300亿美元。 这两张卫星照片分别显示了美国和加拿大部分地区当地时间8月13日晚9时21分（左），及14日9时03分的夜间光亮强度，从中可以看出停电前后这一地区的夜间照明情况的差异。美国东部时间日下午，美国东北部和加拿大部分地区发生大面积停电，波及美加两国的许多城市，给当地交通、通信和居民的生活造成严重影响。纽约市目前已有的地区恢复了电力供应。电力系统的建模负荷发电机组电力网法拉第电磁感应定律克希荷夫电流定律（节点）和电压定律

21、（回路）电力系统的数学模型含n+m+1条母线且第m+1条母线为无穷大的电力系统方程为：微分代数方程DAE(Differential-Algebraic Equation),(0),(pyxgpyxfx nqmnRRf:mqmnRRg:mRy在电力系统中 为动态状态变量，一般是发电机电压和转角； 为瞬时变量，一般是母线电压及其他潮流变量；参数 通常是系统参数，元件参数及负荷和电压设定值等操作参数。nRxqRp关于DAE的稳定性与分岔考虑DAE),(0),(pyxgpyxfx .,),()1(,),(.)1()(.)(,()(),(0),(det|),(0),(|),(:),0000000称为奇异

22、面或不存在附近可能解不唯一在则但若可得系统的唯一解代入将使附近存在唯一在时，由隐函数定理则当记光滑。其中（pppppypppmnqmnSyxSyxxyyxyxxyyxxSyxpyxgyxSpyxgyxRRgf单机无穷大系统SMIB分岔图参数空间奇异诱导分岔SIB(singularity induced bifurcation)000.00,2SIByxyxx处出现在障碍点(impasse point)障碍点，可得后改为类似，将障碍点称为前，时不能再往前，则当的解考虑例)(11)(0011, 1)(, 1)(1)0(, 1)0(0, 12backwardxxforwardttttyttxyxyx

23、x(0,0)xy电力系统的稳定性电力系统的稳定性是指系统在受到扰动后恢复到原始稳态，或达到新的稳态运行的能力。它主要研究电压稳定性，前面所提到的一些重大事故也都是由电压崩溃引起的。可以把电力系统的稳定性分成静态和动态两大类。它们所研究的对象和方法各有不同，如下表。电力系统的稳定性问题静态稳定性电力系统的静态稳定性分析归于求解一组非线性代数方程的潮流分析法。在这里重要的是计算其稳定区域，即可行域（或静态安全域）。如Venkatasubramanian等所指出的，可行域的边界通常是由上述鞍结分岔(SNB)，Hopf分岔(HB)和奇异诱导分岔(SIB)组成。动态稳定性 对电力系统的动态稳定性研究可按

24、动态过程所经历的时间长短而引起电压失稳分成三类： （）零秒（约）10秒，为暂态电压稳定。 （）分钟（多为分钟），为中期电压稳定。 （）几分钟几十分钟，为长期电压稳定。 其中（）常可归结静态稳定性来研究，而（）是当前电力系统稳定性研究的基本课题。混合系统与暂态稳定性混合系统(Hybrid system)为研究电力系统的暂态稳定性，可将DAE分时间段定义，即看成一个混合系统，其一般形式如下。电力系统暂态稳定性分析方法(1) Lyapunov函数法(2) TEF（暂态能量函数法）（Hsiao-Dong Chiang） 关于BCU法（boundary of stability region based

25、 controlling unstable equilibrium point）(3) EAC（等面积准则）和EEAC（扩展的EAC） （薛禹胜） 关于同步稳定性与CCEBC（“互补簇簇标能量壁垒准则”） ISD（孤立稳定域）与NARI（由ISD导诱导的邻域吸引子） 为论证EEAC法的合理性，可以把电力系统的暂态过程近似地看作一个分时间段的简单Hamilton系统， 其中势能函数 V为 nRVM , .),(,0),(, 0),()(postinpretVtVtVV1111*)cos()(niiinjijiijPBV2.4 关于大型发电机组轴承的转子动力学研究 事故：如大同、秦岭等地20万千瓦

26、发电机组的七次（含一次在国外）重大事故。 轴系特点： （1）多自由度； （2）高速、强震动； （3）非线性，如分岔与混沌。 数学模型研究。3.生命科学中的数学建模 3.1 生态模型 设有n个种群（密度为） 相互作用，则一般形式为（Kolmogorov） （1） 其中，记 ，则对 表示 即 （如食饵）对 起促进作用，而 表示 （如捕食者）对 起阻碍作用。, 1,nixinixgxxfxiiii, 2 , 1),()(jiijxgg/0,ijgjiijxxjxixjxix0ijg 通常 。显然， 都是不变的超平面。记 我们感兴趣的是正平衡解 的稳定性（吸引性）。0iig, 2 , 1, 0nixi

27、, 2 , 1, 0;nixRxRinnnRx*捕食者食饵模型 11221122121111111222222212(1),(),(),c xyc xmyxxrxKa xymya xymye c xyyda xymye c mxyyda xymy 考虑两个捕食者和一个食饵的三维情形，这里两个捕食者是对称的，它们之间形成竞争的关系，且三物种严格依赖于比率： 分析表明，该模型通常没有孤立的严格正平衡解，即至少有一个物种会灭绝。但在一定的条件下会出现由平衡解组成的一条平衡曲线。该曲线一边稳定，另一边不稳定，在中间出现一类“无参数分岔”现象，使得系统从总体上是稳定的。捕食系统的整体稳定性3.2 病毒感

28、染模型 SARS病毒概述及其生命周期 病毒在细胞内各生化反应的动态模拟 病毒在细胞间传播的动态模拟 进一步的工作SARS病毒概述 一种新型冠状病毒,属单链正义RNA病毒 特点 增长迅速,易变异 主要结构 基因组RNA 结构蛋白SARS病毒的生命周期病毒在细胞内各生化反应的动态模拟符号:均为大于零的常数。的参数病毒粒子；模型中涉及未知蛋白；蛋白；蛋白；蛋白；蛋白；聚合酶；负链；正链987654321xxSxMxExNxRNAxRNAxRNAx假设RNA聚合酶的合成先于其他生化反应,即设细胞内RNA聚合酶为常数 病毒生命周期的数学模型.)(,8 ,7,6,5,4 ,)(99184991122323

29、121221118432321211xdxxckdtdxixdxxcxkdtdxxdxcxxkxvdtdxxdxxcxcxxkxvdtdxiiiiiiiiiiii系统由平衡点(0,0,0,0,0,0,0,0),且当时该平衡点稳定.0212121kkvvdd该系统有零平衡点和一个正平衡解,)(1,8,7,6,5,4,)(849911122122eiieiieeiieieieeexxcdxixcdxkxxkdxvx8321111221121.0iieiieiieeidxcdxkcxkdkxvvx）（满足：其中数值结果 正链RNA N-蛋白病毒在细胞间传播的动态模拟 考虑易感染细胞, 已感染细胞,自

30、由病毒粒子和免疫细胞间的关系 数学模型.)() )(;,),(;,1849的病毒粒子数是每个被感染细胞释放为衰减率为免疫参数是病毒感染细胞的参数细胞的生产速率为易感染是单调增加的且关于描述免疫细胞的变化胞的数量细自由的病毒粒子和免疫已感染细胞分别表示易感染细胞其中txtxckbdpkrVVmfmVISiiiii,),(,4312111mdpVmVmfmVdpVmVSkbIVIdVSkISdVSkrS . 0),(.,0), 0 , 0 ,(. 0,),(.)(,11111112321*221849VImVISrbkrddddRmSakVkVaVmfxxckbeeiii满足出其他平衡点另外可进一

31、步由系统解该平衡点稳定时当则系统有平衡点取为常数即令间的传播度远远大于病毒在细胞设细胞内病毒的增长速数值结果 易感染细胞数量 已感染细胞数量 病毒粒子数量 免疫细胞数量S,I,V,m四种粒子与参数 的依赖关系 S与 的依赖关系 1kp和1kp和I与 的依赖关系 1kp和3.3 关于生物信息学和系统生物学 生物信息 系统生物学 模型：随机微分方程等4.分形模型 分形是简单空间(如欧氏空间)中具有某种精细结构的复杂集合，其特点为： (1) 具有自相似结构； (2) 不同于传统的几何图形, 不是某些简单方程的解集, 但常可通过对较简单的变换作迭代来产生. (3) 需用分维数来度量, 其维数通常大于相

32、应的拓扑维； (4) 具有混沌性质.法国的Mandelbrot.B 开创了分形几何1967年的论文：“英国海岸线的长度不确定” （fractal geometry）的研究（1）具有无限嵌套层次的精细结构对自然几何形态的数学研究海岸线的长度随测量尺度变化（2）在不同尺度下具有某种相似特性科赫雪花科赫雪花 维数d=log4/log3=1.26186Koch 雪花曲线设E0为单位直线段三等分后，中间一段用与其组成等边三角形的另两边代替，得到E1对E1的4条线段的每一条重复以上做法，得到E2以此方法重复，可得En当n趋于无穷，得到的极限曲线就是Koch 曲线用Mathematica 画koch曲线re

33、dokochptlist_List := Blocktmp = , i, pnum = Lengthptlist, Fori = 1, i Sqrt3/6自相似性精细结构：复杂性不随尺度减小而消失处处不光滑，每一点是尖点长度：En的长度(4/3)n趋于无穷本身定义方式简单Koch 曲线的特点Koch曲线在有限区域却长度无限，它具有分维数。单参数的函数曲线是一维的吗？设是平面上边长为1/2的正三角形，构造 fnf1f2f3以此方式得到 fn ，在0,1一致收敛到极限函数 f的象将为整个三角形分维数将单位边长的线段，正方形，立方体分成边长为1/2的同样几何物体，得到21，22，23个小线段，正方形，立方体注意指数给出了几何物体的维数若将几何物体的长度（线度）缩小为1/r，定义分形维数得到N个相似小几何物体，那么维数d满足N=rdd=logN/log rKoch曲线的维数？约1.2618Cantor集从单位区间0,1出发，三分去中段，得E1，E1两个区间三分去中得E2 ，极限集合为Cantor集 这是一个完备的、完全不连通、具有自相似的精细结构的集合，其长度为0。康托尔三分集合维数d=l