复杂网络上动力系统同步的研究进展[1]

1复杂网络上动力系统同步的研究进展赵明，汪秉宏,蒋品群，周涛中国科学技术大学近代物理系，合肥230026摘要：本文简要介绍复杂网络的基本概念并详细总结了近年来复杂网络上动力学系统的同步的研究进展，主要内容有复杂网络同步的稳定性分析，复杂网络上动力学系统同步的特点，网络的几何特征量对同步稳定性的影响，以及提高网络同步能力的方法等。最后文章提出了这一领域的几个有待解决的问题及可能的发展方向。关键词:复杂网络；同步；小世界网络；无标度网络中图分类号:O261文献标识码:A0引言自然界中存在的大量复杂系统都可以通过形形色色的网络加以描述[1-6]。一个典型的网络是由许多节点与连接两个节点之间的一些边组成的，其中节点用来代表真实系统中不同的个体，而边则用来表示个体间的关系，往往是两个节点之间具有某种特定的关系则连一条边，反之则不连边，有边相连的两个节点被看作是相邻的。例如，神经系统可以看作大量神经细胞通过神经纤维相互连接形成的网络[7]；计算机网络可以看作是自主工作的计算机通过通信介质如光缆、双绞线、同轴电缆等相互连接形成的网络[8]；类似的还有电力网络[7]、社会关系网络[7,9-10]、交通网络[11]等等。对网络最早进行研究的是数学家，其基本理论是图论。经典的图论总是倾向于用某种规则的拓扑结构模拟真实网络[12-13]，到了二十世纪五十年代末期，Erdös和Rényi建立了随机网络的基本模型[14-16]，后来的近半个世纪，随机图一直是科学家研究真实网络最有力的武器[17]。直到最近几年，由于计算机数据处理和计算能力的飞速发展，科学家们发现大量的真实网络既不是规则网络，也不是随机网络，而是具有与前两者皆不同的统计特征的网络。这样的一些网络被科学家们叫做复杂网络，其诸多统计特征中最重要的是小世界(smallworld)效应[7,18]和无标度(scalefree)特性[19-20]。在网络中，两点间的距离被定义为连接两点的最短路径所包含的边的数目，把所有节点对的距离求平均，就得到了网络的平均距离。另外一个叫做簇系数的参数，专门用来衡量节点集聚成团的情况。对于某个节点，它的簇系数被定义为它所有相邻节点之间连边的数目占可能的最大连边数目的比例。类似的，网络的簇系数C则是所有节点簇系数的平均值。研究表明，规则网络具有大的簇系数和大的平均距离，随机网络则具有小的簇系数和小的平均距离。1998年，Watts和Strogatz通过以某个很小的概率改变规则网络中边的连接方式构造出了一种介于规则网络和随机网络之间的网络（WS型小世界网络），它同时具有大的簇系数和小的平均距离，因此既不能当作规则网络处理，也不能被看作是随机网络[7]。后来物理学家把大的簇系数和小的平均距离两个统计特征合在一起称为小世界效应，具有这种效应的网络就是小世界网络(smallworldnetwork)。大量的实证研究表明，真实网络几乎都具有小世界效应[1-11]，有的甚至具有所谓的超小世界效应[21]，同时科学家还发现大量真实网络的节点度服从幂率分布[1-6,8,10,21-23]。节点度是指一个节点拥有相邻节点的数目，节点度服从幂律分布就是说具有某个特定度的节点数目与这个特定的度之间的关系可以用一个幂函数近似地表示。幂函数曲线是一条下降相对缓慢的曲线，这使得度很大的节点可以在真实网络中存在。对于随机网络和规则网络，度分布区间非常狭窄，几乎找不到偏离节点度均值较大的点，故其节点度均值可以被看作其节点度的一个特征标度。由于幂函数具有标度不变性，因此我们把节点度服从幂律分布的网络叫做无标度网络(scalefreenetwork)，并称这种节点度的幂律分布为网络的无标度特性。对于物理学家而言，研究复杂网络的终极目标之一是理解网络拓扑结构对物理过程的影响[2]。网络上的同步现象是一项重要的研究课题。人们已观测到的同步现象包括夏日夜晚青蛙的齐鸣、萤火虫的同步发光，心肌细胞和大脑神经网络的同步[24-26]，剧场中观众鼓掌频率的逐渐同步[27]，等等。在以前的研究中，人们忽略了网络的拓扑性质，在研究同步问题时，自然地选择了最容易模拟和分析的规则网络或随机网络，而没有仔细思考和研究这种选择是不是应该的，不同的选择会不会对物理过程产生不可忽略的影响。当然，如果理论研究和实验结果都说明复杂网络的同步情况与规则或随机网络别无二致，那么我们至少暂时还可以心安理得地使用以前的结论。但是，不幸的是，复杂网络上的同步能力与规则或随机网络相比确实存在明显的不同。类似的情况还出现在其他的物理过程中，在接下来的几节中我们将详细地介绍网络拓扑性质对网络同步能力的影响。1网络同步的基本概念如果在网络的每个节点上加上一个动力学系统，这个动力学系统既可以是极限环也可以是混沌的；而让有边相连的两个节点的动力学系统之间存在相互的耦合作用，就形成了一个动力学网络。严格地说，设网络有N个节点，第i个节点在n时刻的m维状态变量是，单个节点在不考虑耦合作用的时候所满足的状态方程是。是每个节点状态变量的函数，用于对其它节点进行耦合。这样，在存在耦合作用的情况下，第i个节点所满足的状态方程是)(nix))(()1(nniixFx=+mmRR→:H))(())(()1(nGnnjjijiixHxFx∑+=+σ，(1)对于连续系统)()(jjijiiGxHxFx∑+=σ，(2)其中σ是耦合强度，表示耦合矩阵G的矩阵元，定义如下：ijG⎪⎩⎪⎨⎧−=01iijkG，(3)otherwisejjiiΛ∈=其中是节点i的度，是与节点i相邻的节点的集合。耦合矩阵G包含了网络结构的全部信息。例如，最近邻耦合网络、星型网络以及完全网络的耦合矩阵分别为：ikiΛ⎟⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎜⎝⎛−−−−=21010210012110121#####G，，⎟⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎜⎝⎛−−−+−=10010101001111112#####NG⎟⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎜⎝⎛+−+−+−+−=11111111111111113NNNN#####G。在耦合的作用下，经过一段时间的演化，使得，网络就进入了同步状态。sxxx====N21当然并不是所有的网络在任意耦合强度或耦合方式下都能实现同步。在下面的几节中，我们将仔细讨论网络同步的稳定性问题，网络实现同步的条件和网络拓扑结构对同步的影响等等。23jz2复杂网络同步的稳定性分析Pecora和Carroll研究了线性耦合网络同步的稳定性问题，给出了主稳定函数判据[28]。他们首先假设：(1)所有的耦合振子都是完全相同的，(2)从每个振子提取的用于耦合其它振子的函数也是完全相同的，(3)同步流形是不变流形，(4)节点的耦合方式使在同步流形附近可以线性化。假设(1)和(3)是为了保证相空间中同步超平面的存在，假设(2)是为了使动力学系统和网络结构的稳定性图象更加清晰具体，假设(4)是为了更好地应用线性近似这一研究耦合系统最常用的方法。在此基础上，Pecora和Carroll逐步完成了当网络上的耦合振子系统的同步混沌态存在短波分岔时的同步稳定性分析，提出用主稳定性函数方法确定动力学网络同步的稳定性[28-30]。下面我们介绍主稳定性函数的定义及其使用方法。首先对动力学网络的同步流形进行线性稳定性分析，很多人已经做了这方面的工作[28-38]，我们主要采用Chen等人[38]的叙述过程。当存在耦合作用时第i个节点所满足的状态方程是(2)，在同步状态s附近对其进行线性化，得到()()iiijjDGDσ=+∑zFszHs，(4)其中是节点i在同步流形上的变分，izs)(⋅FD和)(⋅HD分别是函数和的阶雅可比(Jacobian)矩阵。利用阶矩阵重写(4)式，得FHmm×Nm×),,,(21Nzzzz=()()TDDσ=+⋅⋅zFszHszG。(5)根据约丹规范型(Jordancanonicalforms)理论，上式的稳定性是由G的特征值γ决定的，设其对应的特征向量为e，并且令，将e右乘(5)式，得到=⋅uze(()()DDσγ=+uFsHs)⋅u。(6)这样原来要讨论的维空间的稳定性问题被简化到Nm×mm×维空间，并且通常情况下Nm。以上是利用连续系统进行讨论，离散系统可给出类似(6)式的结论。我们可以利用(6)式计算单个系统的李雅普诺夫指数(Lyapunovexponents)，设这些指数分别为mλλλλ≥≥≥=2max1。注意到由于∑==NjijG10，0=γ总是G的一个特征值，相应的特征向量是，它对应同步流形模式。其它T)111(1−N个特征向量所张成的子空间横截于同步流形，如果所有这些横截李雅普诺夫指数都小于零，系统稳定。设βασγi+=(采用复数是因为耦合矩阵G的特征值中可能有复数，如非对称耦合)并代入(6)中，(()()()DiDαβ=++uFsHs)⋅u，(7)计算最大李雅普诺夫指数maxλ随α和β的变化关系，这就是Pecora和Carroll定义的主稳定性函数(masterstabilityfunction)。图1(a)(摘自文献[28])给出了耦合振子系统的同步混沌态存在短波分岔时它们的关系图，0maxλ定义了稳定区域。下面举个例子来说明这个过程。以Rössler系统[39])()(cxzbdtdzayxdtdyzydtdx−+=+=+−=(8)作为节点上的动力学系统，令a＝b＝0.2，c＝7.0，使系统处于混沌区。通过变量x进行耦合，即⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛=000000001H。(9)图1给出了通过变量x进行耦合的Rössler振子的主稳定函数，该图是关于实轴对称的，所以只画出了实轴上方部分。可以看到半圆形的曲线包围着稳定区域。对于0αβ==，因为各个节点之间没有耦合，而Rössler系统又是混沌的，所以max0λ。对于0β=，当(0)α减小时，maxλ先是由正变为负，而当α进一步减小时，maxλ又由负变为正。这就表明，过弱和过强的耦合都会使同步状态不稳定。如果在同步区域内固定α的值而增大虚部β的绝对值，maxλ也会由负变为正。这也表明，耦合的不对称性增强也会使同步状态不稳定。以上三种使网络失去同步的方式分别叫做长波分岔(long-wavelengthbifurcationorLWB)，短波分岔(short-wavelengthbifurcationorSWB)和中波分岔(intermediate-wavelengthbifurcationorIWB)。图1Rössler振子通过变量x进行耦合的主稳定函数[参考文献28]Barahona和Pecora[37]深入研究了对称耦合的情况。他们是用拉普拉斯算子(Laplacian)来表示网络的结构的，与前面的耦合矩阵相差一个负号，其特征值是实数且011ma0Nxγγγγ−=≤≤≤=。图2(摘自文献[37])清楚地表达了maxλ随α的变化规律，那么只要满足12(,),1,2,,1kkNσγαα∈=−，整个网络就能够实现同步。因此，线性稳定同步状态的存在可以用一个几何条件来表示：⎩⎨⎧2max11ασγασγ，(10)即βααγγ≡121max，(11)β的大小是由单个节点的状态方程F和变量的耦合方式H决定的，对于大多数混沌振子其值介于5到100之间[37]。这样，在研究动力学网络的同步稳定性问题时，可以将节点的动力学性质与网络的耦合矩阵的特征值谱分开来讨论：首先确定β的取值范围，然后计算网络的拉普拉斯矩阵的特征值比max1/γγ，根据(11)式就能确定在这样的网络上对于选定的动力学系统，是否能够实现动力学网络的同步。在分析过程中会发现这样的现象：能够使一个动力学系统同步的网络不一定能够使另一个系统实现同步，[40]中给出了这样的例子。汪小帆和陈关荣[36]研究了耦合振子是连续系统的复杂网络的同步稳定性问题，得到下面的定理(定理中一些符号的使用与上文有所不同，尊重原文我们采用所摘文献中使用的符号，但同时用文字表明其物理意义，下同)。设网络的状态方程是()iiijjca=+Γ∑xfxxj，(12)4其中nnR×Γ∈是