动力系统建模课件

动力系统建模

唐云（清华大学数学科学系）

ytang@引言1。动力系统的基本概念与方法：从蝴蝶泉到蝴蝶效应2。机械与电力系统中的数学模型3。生命科学的数学建模4。分形模型目录引言动力系统是关于系统演化规律的数学学科。数学建模是连接动力系统理论与应用的桥梁。特点：（1）广泛的应用前景；（2）深刻的数学理论基础；（3）计算技术。动力系统建模作用：（1）教学：数学建模；（2）科学研究：与国民经济及科学前沿关系密切。从蝴蝶泉到蝴蝶效应

1.1蝴蝶的生态问题

云南大理蝴蝶泉，是影片《五朵金花》里阿鹏和金花对歌的地方，也有名的游览胜地。泉内蝴蝶种类繁多，每年农历4月15日白族的“蝴蝶会”前后，蝴蝶大的大如巴掌，小的小如蜜蜂，成串悬挂于泉边的合欢树上，盛况空前。明代徐霞客在其“游记”中称：“真蝶万千，连须钩足，自树巅倒悬而下及于泉面”。郭沫若在1961年游蝴蝶泉时也曾留下“首尾联接数公尺，自树下垂疑花序”的诗句。

令人惋惜的是，近十数年，人们已经很难看到美丽的蝴蝶盛会，有时，虽有蝴蝶聚集，但数量已少。据当地父老传言：蝴蝶泉边，原有一蓬枝叶茂密、开白花、发清香的茨蓬，花枝缠在横斜泉面的树干上，蝴蝶沿着这些下垂的花枝连成串。如今，茨蓬已除，泉面树干叶枯，加上周围自然环境受到破坏，田野大量使用农药，误伤不少蝴蝶，那连须钩足悬于泉面的奇观，久已不见。1.2数学模型:Logistic映射f(x)=ax(1-x)，x在[0,1]内变化xn+1=f(xn)

从[0,1]内点x0出发，由Logistic映射的迭代形成xn=fn(x0),n=0,1,2,…序列{xn}称为x0的轨道。由轨道迭代的极限集可组成分岔图种群数的模型简化：相应的迭代为了一个序列，即Logistic映射Logistic映射分岔图关于吸引子

吸引子是动力系统相空间中的一类特殊集合，它能把周围的轨道都“吸引”（收敛）过来。吸引子分类：（1）平衡点（2）周期轨（3）拟周期轨（4）混沌吸引子■当0<a<1时，由于■当1<a<3时，任何(0,1)中初始值的轨道趋于

两个不动点x1*,

x2*,一个稳定(吸引),另一个xn

→0物种逐渐灭亡x*=1-1/a其中x*是方程f(x)=x的解，为映射f

的不动点（周期1点)例：a=1.5时

xn→1/3.不稳定，轨道{xn}趋向稳定点。数值迭代：倍周期分岔

■当1+61/2<a<3.5440903506…时,从任意的点x0出

也称为周期2点,对应轨道称周期2轨道.(原来周期点失稳)发的轨道将逐渐沿着四个数值振动，它们满足称为周期4点，对应轨道称周期4轨道(原有周期点

若a再增大，周期4点又会失稳，而产生新的稳定又失稳)周期8点，这个周期不断加倍的过程将重复无限次，会依次出现周期16点，周期32点，….，这种过程称为倍周期分岔.相应的分岔值c1=3,c2=1+61/2…构成一个单调增加的数列{ck}.其极限值为c*=3.569945557391…。

当c*<a<4时，Logistic映射进入混沌区域.反映出■遍历性：点

x0的轨道不趋向任何稳定的周期性,即不同初始值，即使它们离得非常近，它们的的是：轨道,它的轨道在(0,1)(或其中某些区间)内的任何一个子区间(a,b)内都会出现无数次.■敏感性：轨道表现出对初始条件的强烈敏感轨道也终将以某种方式分离.混沌的特点■

Feigenbau常数(ck-ck-1)/(ck+1-ck)在k趋于无穷时，趋于常数

q=4.6692016这常数的意义在于普适性，例如周期3窗口，还有■存在周期窗口：混沌区域内某些地方仍有倍周期分岔，例如a＝3.835附近其他映射

任取(0,1)中的点x0,可以通过作图来取得迭代

在以xn为横坐标、xn+1为纵坐标的第一象限作抛物线弧：xn+1＝a

xn(1-xn)的数值序列{xn},从而也通过图象直观地看出由x0出发的轨道的变化.这作图的过程颇象蜘蛛织网,故称为蛛网迭代.

图像方法：蛛网迭代11xnxn＋1x0x1x1x2■

1<a<3

从(0,1)中任何初值出发的轨道趋向不动点(周期1点)

■3<a<61/2+1

从任何初值出发的轨道趋向周期2点■61/2+1<a<

3.54409035从任何初值出发的轨道趋向周期4点■

a=3.58轨道进入浑沌状态■

a=4轨道的浑沌性表现充分

蛛网迭代的优点是轨道非常直观形象.缺点是当周期数较大时不易看清轨道变化细节

密度分布图：■从密度：从一个初始点

x0出发，由迭代所产生的序列{xn}(n一般很大)在区间[0,1]上的概率分布密度.■将具体算法：将[0,1]区间分成m个长度为h=1/m的小区间，序列{xn}nN=0落在各个小区间[ih,(i+1)h]的个数为ki,则该序列落在各小区间的概率（即密度)为pi=ki/N

i=0,1,2,…,m■密度图：横轴为区间[0,1],纵轴为概率

p.每个小区间上的细柱线的高度等于该区间上密度■

a=3.2(m=100N=10000x0=0.1)(这是周期2情况)■

a=3.45(这是周期4情况)■

a=3.55(周期8的情况)

以上密度图显示在

0<a<c*的情况下,{xn}只有极少数落在周期点以外的小区间,而最终以几乎相等的概率落在周期点所在的小区间。■

a=3.6(进入混沌区)

(最混沌状态)■

a=4Logistic模型的混沌自相似（分形）图一

图二图三：图二局部放大图图四：图三局部放大图1.3关于“蝴蝶效应”1961年冬天E.Lorenz进行关于天气预报的计算。他考虑下面加热的流体由热传导进入对流，然后产生湍流的过程，对Rayleigh-Bernard方程进行约化，得到下面的Lorenz方程。D.Gulick,EncounterswithChaos,Mc-GrawHill,Inc.,NewYork,1992.Lorenz方程

和

为正数.

,,和

与流体的物理性质相关.Lorenz取

,,和

.Lorenz吸引子对初值条件的敏感依赖性10,000个几乎相同的初值条件从这10,000个初值出发的每条轨道经过同样时间后其终点很不一致.这些点都在同一个吸引子范围内.S.H.Strogatz,NonlinearDynamicsandChaoswithApplicationstoPhysics,Biology,ChemistryandEngineering,Addison-WesleyPublishingCompany,NewYork,1994.所谓“蝴蝶效应”是指：

初始值很小的改变会引起绝然不同的结果.

巴西的一只蝴蝶扇动翅膀会引起

明年在得克萨斯的大风暴吗？E.Lorenz

数学的伟大使命在于从混沌中发现秩序。—倍尔在那个混沌的体制中，结构上的微小差异几乎都会造成行为方式上的巨大变化，可控制的行为似乎已被排除。—斯图尔特.考夫曼关于混沌一个动力系统是混沌的，如果它满足：(1)有一个由周期轨道组成的稠密集合；(2)轨道敏感地依赖于其初值条件；(3)为拓扑传递的.混沌的度量性质：(1)正Lyapunov指数(2)拓扑熵与测度熵(3)分形结构混沌出现在各个领域的一种现象:数学、物理、由此引起的复杂而有趣的现象“侏罗纪公园”中的恐龙重现生物、金融、经济、管理等等：宇宙的起源

龙卷风的产生、厄尔尼诺现象东南亚金融危机爆发可以从某些简单的离散的数学模型开始,讨论2.机械和电力系统的数学模型

2.1动力学模型Newton力学体系是第一个，也是最基本的动力系统数学模型。建模过程：（1）数据积累：第谷(TychoBrahe,1546-1601)（2）经验公式：开普勒(J.Kepler,1571-1630)的行星“三大定律”。（3）数学模型：牛顿(I.Newton,1642-1727)的“万有引力”。（4）验证：如哈雷彗星和海王星的发现。非线性振动Duffing方程Duffing方程位移x位移x时间Tdx/dt位移x时间Tdx/dt位移x

耗散系统相体积的演化

解释初值敏感和奇怪吸引子要用到相体积的伸展与折叠，定量描述这一特征的量是李雅普诺夫指数，这要证明三维（以上）相体积

而中要有正有负。讲混沌时专门解释或用初等的例子说明不严格，并且要花费一定学时。但只要在前面讲正则方程和刘维定理时稍做改动就可以自然引出耗散系统相体积演化公式，概念的引入严格、定量，学时反而减少。通常讲刘维定理，对相空间保守体系其证明用到过去教材都是代入保守体系的正则方程。实际上，对一力学体系，如果除保守力外还含非保守力，则正则方程应写为其中代表非保守力，代入前式后得出，积分马上得到容易证明，对保守力对非保守力（如）对高维耗散系统，自然导出指数形式形式其中（i=1～f）可正可负总和为正。极限环与分岔考虑一个数学例子：

用二维极坐标表示，作代换得到当时有方程有三个根：

，，，点是稳定定态（焦点）

（硬激励）稳定定态（焦点）；

不稳定的极限环；

稳定的极限环。黑洞（软激励）

是稳定的极限环不稳定的焦点极限环分岔叉式分岔：定性举例，旋转单摆

Hopf分岔：点到极限环的突变

稳定定态；

稳定极限环；

终态稳定点与初始条件有关（亚临界Hopf分岔）。

倍周期分岔:周期成倍突变

阻尼单摆的强迫振动方程引入无因次量后化为

:时，稳定的定态，，不稳定点，

摆长为l，小球质量为m的单摆，相对与平衡的下垂位置的角位移为θ，重力加速度为g，则其运动方程为：

(1)或

(2)

等式右边是周期性的驱动力，其角频率为。把方程式无量纲化，用去除每一项，将无量纲的时间叫做t，即得

(3)式中都是无量纲化的。

阻尼单摆的强迫振荡旋转数其中Henon-Heiles星体势模型等势图粒子真实轨迹粒子的庞加莱截面激励转子2.3关于电力系统的数学模型电力工业是国民经济的支柱产业，并且直接影响到人民的生活。我国电力工业发展迅速：194919902001装机容量（亿千瓦）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时许，俄亥俄州北部34.5万伏超高压突然发生故障。一小时内波及包括纽约和多伦多在内的美加东中部大停电，5000万人陷入黑暗之中。至16日10时基本恢复正常。估计经济损失每天达300亿美元。这两张卫星照片分别显示了美国和加拿大部分地区当地时间8月13日晚9时21分（左），及14日9时03分的夜间光亮强度，从中可以看出停电前后这一地区的夜间照明情况的差异。美国东部时间１４日下午，美国东北部和加拿大部分地区发生大面积停电，波及美加两国的许多城市，给当地交通、通信和居民的生活造成严重影响。纽约市目前已有８５％的地区恢复了电力供应。电力系统的建模负荷发电机组电力网法拉第电磁感应定律克希荷夫电流定律（节点）和电压定律（回路）电力系统的数学模型含n+m+1条母线且第m+1条母线为无穷大的电力系统方程为：微分代数方程DAE(Differential-AlgebraicEquation)在电力系统中为动态状态变量，一般是发电机电压和转角；为瞬时变量，一般是母线电压及其他潮流变量；参数通常是系统参数，元件参数及负荷和电压设定值等操作参数。关于DAE的稳定性与分岔考虑DAE单机无穷大系统SMIB分岔图参数空间奇异诱导分岔SIB

(singularityinducedbifurcation)障碍点(impassepoint)(0,0)xy电力系统的稳定性电力系统的稳定性是指系统在受到扰动后恢复到原始稳态，或达到新的稳态运行的能力。它主要研究电压稳定性，前面所提到的一些重大事故也都是由电压崩溃引起的。可以把电力系统的稳定性分成静态和动态两大类。它们所研究的对象和方法各有不同，如下表。电力系统的稳定性问题静态稳定性电力系统的静态稳定性分析归于求解一组非线性代数方程的潮流分析法。在这里重要的是计算其稳定区域，即可行域（或静态安全域）。如Venkatasubramanian等［１９９６］所指出的，可行域的边界通常是由上述鞍结分岔(SNB)，Hopf分岔(HB)和奇异诱导分岔(SIB)组成。动态稳定性对电力系统的动态稳定性研究可按动态过程所经历的时间长短而引起电压失稳分成三类：（１）零秒～（约）10秒，为暂态电压稳定。（２）１～５分钟（多为２～３分钟），为中期电压稳定。（３）几分钟～几十分钟，为长期电压稳定。其中（３）常可归结静态稳定性来研究，而（１）是当前电力系统稳定性研究的基本课题。混合系统与暂态稳定性混合系统(Hybridsystem)为研究电力系统的暂态稳定性，可将DAE分时间段定义，即看成一个混合系统，其一般形式如下。电力系统暂态稳定性分析方法(1)

Lyapunov函数法

(2)

TEF（暂态能量函数法）（Hsiao-DongChiang）

关于BCU法（boundaryofstabilityregionbasedcontrollingunstableequilibriumpoint）

(3)

EAC（等面积准则）和EEAC（扩展的EAC）（薛禹胜）

关于同步稳定性与CCEBC（“互补簇簇标能量壁垒准则”）

ISD（孤立稳定域）与NARI（由ISD导诱导的邻域吸引子）

为论证EEAC法的合理性，可以把电力系统的暂态过程近似地看作一个分时间段的简单Hamilton系统，

其中势能函数V为2.4关于大型发电机组轴承的转子动力学研究事故：如大同、秦岭等地20万千瓦发电机组的七次（含一次在国外）重大事故。轴系特点：（1）多自由度；（2）高速、强震动；（3）非线性，如分岔与混沌。数学模型研究。3.生命科学中的数学建模3.1生态模型设有n个种群（密度为）相互作用，则一般形式为（Kolmogorov）（1）其中，记，则对表示即（如食饵）对起促进作用，而表示（如捕食者）对起阻碍作用。通常。显然，都是不变的超平面。记我们感兴趣的是正平衡解的稳定性（吸引性）。捕食者－食饵模型：

考虑两个捕食者和一个食饵的三维情形，这里两个捕食者是对称的，它们之间形成竞争的关系，且三物种严格依赖于比率：

分析表明，该模型通常没有孤立的严格正平衡解，即至少有一个物种会灭绝。但在一定的条件下会出现由平衡解组成的一条平衡曲线。该曲线一边稳定，另一边不稳定，在中间出现一类“无参数分岔”现象，使得系统从总体上是稳定的。捕食系统的整体稳定性3.2病毒感染模型SARS病毒概述及其生命周期病毒在细胞内各生化反应的动态模拟病毒在细胞间传播的动态模拟进一步的工作SARS病毒概述一种新型冠状病毒,属单链正义RNA病毒特点

增长迅速,易变异主要结构

基因组RNA

结构蛋白SARS病毒的生命周期病毒在细胞内各生化反应的动态模拟符号:假设RNA聚合酶的合成先于其他生化反应,即设细胞内RNA聚合酶为常数病毒生命周期的数学模型系统由平衡点(0,0,0,0,0,0,0,0),且当时该平衡点稳定.该系统有零平衡点和一个正平衡解数值结果正链RNAN-蛋白病毒在细胞间传播的动态模拟考虑易感染细胞,已感染细胞,自由病毒粒子和免疫细胞间的关系数学模型

数值结果易感染细胞数量已感染细胞数量病毒粒子数量免疫细胞数量S,I,V,m四种粒子与参数的依赖关系S与的依赖关系I与的依赖关系3.3关于生物信息学和系统生物学生物信息系统生物学模型：随机微分方程等4.分形模型

分形是简单空间(如欧氏空间)中具有某种精细结构的复杂集合，其特点为：

(1)具有自相似结构；

(2)不同于传统的几何图形,不是某些简单方程的解集,但常可通过对较简单的变换作迭代来产生.(3)需用分维数来度量,其维数通常大于相应的拓扑维；

(4)具有混沌性质.法国的Mandelbrot.B开创了分形几何1967年的论文：“英国海岸线的长度不确定”

（fractalgeometry）的研究（1）具有无限嵌套层次的精细结构对自然几何形态的数学研究海岸线的长度随测量尺度变化（2）在不同尺度下具有某种相似特性科赫雪花

维数d=log4/log3=1.26186Koch雪花曲线设E0为单位直线段三等分后，中间一段用与其组成等边三角形的另两边代替，得到E1对E1的4条线段的每一条重复以上做法，得到E2以此方法重复，可得En当n趋于无穷，得到的极限曲线就是Koch曲线用Mathematica画koch曲线redokoch[ptlist_List]:=Block[{tmp={},i,pnum=Length[ptlist]},For[i=1,i<pnum,i=i+1,tmp=Join[tmp,{ptlist[[i]],ptlist[[i]]*2/3+ptlist[[i+1]]/3,(ptlist[[i]]+ptlist[[i+1]])/2+{ptlist[[i]][[2]]-ptlist[[i+1]][[2]],ptlist[[i+1]][[1]]-ptlist[[i]][[1]]}*Sqrt[3]/6,ptlist[[i]]/3+ptlist[[i+1]]*2/3,ptlist[[i+1]]}]];tmp]Inko01={{0,0},{1,0}};Show[Graphics[Line[Nest[redokoch,Inko01,4]],AspectRatio->Sqrt[3]/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/logrKoch曲线的维数？约1.2618Cantor集从单位区间[0,1]出发，三分去中段，得E1，E1两个区间三分去中得E2

，极限集合为Cantor集

这是一个完备的、完全不连通、具有自相似的精细结构的集合，其长度为0。康托尔三分集合维数d=log2/log3=0.6309Sierpinski集合三角形四等分去中间小三角形所得极限图形维数＝？Weierstrass

函数W(x)=(s－2)ksin(kx),>1,1<s<2

数学分析中的著名例子：处处连续，但无处可微lambda=2;nmax=20;s=1.2;Plot[Sum[lambda^((s-2)k)Sin[(lambda^k)x],{k,1,nmax}],{x,-1,1}]使用Mathematica给s以不同的值的函数，自仿射S=1.2S=1.5S=1.99S=1.7复变函数的迭代Julia集：固定考虑Zk+1=Zk2+给定复数初值Z0,

，得到无穷复数序列{Zk}J