一门崭新的交叉科学_网络科学_下篇_

文章编号:100020542(2007)0420361288

物　理　学　进　展

PROGRESSINPHYSICS

Vol.27,No.4Dec.,2007

一门崭新的交叉科学:网络科学(下篇)

方锦清1,汪小帆2,郑志刚3,李　翔2,狄增如3　毕　桥1

(国家自然科学基金重点项目联合网络研究组1中国原子能科学研究院,北京102413;

2上海交通大学,上海200240;3北京师范大学,北京100875)

摘　要:在上篇综述里,我们评述了网络科学的发展简史,基本概念和网络分类,以及国内外迄

今提出的复杂网络的主要理论模型及其拓扑特性。着重总结了我们“一院二校”开展的国家自然科学基金重点项目的研究进展,上篇的研究成果主要涉及加权网络模型、和谐统一的混合模型和量子信息网络模型及其纳米网络,简单说,我们探索了从宏观网络到微观网络的不同层次的若干特点和结果。在下篇综述里,我们将系统概述5个重要课题的研究进展,从第7章到第11章,课题内容包括:复杂网络的动力学完全同步与网络的拓扑结构之间的关系;网络拓扑结构的对称性破缺导致动力学部分同步;给出判断给定网络部分同步斑图稳定性的方法和一般判据,讨论了由网络的多种拓扑对称导致的部分同步斑图的竞争和选择问题以及李指数谱的简并性质。进一步,分别以具有小世界和无标度拓扑结构的束流输运网络为典型的“混沌复杂网络”,提出了实现束晕-混沌的同步和控制的若干方法,特别是实现分区网络的多目标的同步控制问题。同时,基于Vicsek模型和Boid模型,研究了生物体集群动态网络的拓扑结构和动力学方程,提出复杂多智能体网络的协调控制策略。另外,从复杂网络的不同拓扑结构对合作涌现和资源分配的作用角度,评述了三个主要的博弈模型—少数者博弈、囚徒困境和雪堆博弈,研究了复杂网络的群聚系数和网络异质性等结构特性对个体的博弈行为的重要影响。最后第12章,我们简介其它课题的进展,讨论了网络科学面临的挑战,并展望了其应用与发展前景。

关键词:复杂网络的完全同步和部分同步;具有小世界和无标度的束流传输网;多目标控制与同

步;多智能体网络;网络上博弈问题;挑战性课题与展望

中国分类号:O231.3;N93　　　文献标识码:A

7　复杂网络同步的若干研究进展

同步现象广泛存在于自然界和社会生活,所以在自然科学和社会科学中同步问题的理论

收稿日期:2007206215

基金项目:国家自然科学基金重点项目(批准号:70431002)

362物理学进展　　27卷

及其应用被广泛研究。如,二个靠近的钟摆出现同步,大群萤火虫能够同步闪光和熄灭[192];观众掌声从杂乱到有节奏地同步鼓掌[193,194],以及耦合激光同步,等等。同步现象的研究历史由来已久,特别是,在20世纪90年代发现混沌同步后,由于具有应用潜力,国内外掀起了研究混沌同步的热潮。现在,研究复杂网络的同步是混沌同步研究的深化和拓广,一直是普遍重视的一个重要的研究课题,已有不少有关同步的综述和著作[21,22,195]。近期《Nature》杂志上

[196,197]

的文章指出,纳米耦合振子之间发生了同步行为,这有可能用于研制新的无线通信元件。但是,同步现象有时产生有害作用,如成千上万的人同时过桥会引起了桥体的同步振动具有危险性;Internet上路由器周期性发送路由信息而导致网络通信拥塞。假定一个网络中所有个体的状态都是周期变化的,这种同步现象都可用数学语言来描述,其中每个个体是一个动力学系统,而诸多的动力学个体之间存在着多种特定的耦合关系,产生复杂的斑图结构及多样的动力学行为。迄今同步大致分为完全同步和部分同步(广义同步)两大类。本章侧重于讨论复杂网络的完全同步问题,下章将探讨耦合混沌振子网络中的部分同步现象。近年来,由于复杂网络的两个重要发现—小世界特性[1]和无标度特性[3],引起人们关注大规模实际网络的拓扑结构特性对网络同步能力的影响,主要考虑的是无权无向的网络。本章重点介绍我们自己在该领域的研究进展,特别是在加权网络上的同步[5,198,199]和如何提高网络的同步化性能方面的成果[200,201]。7.1　同步的基本概念为

[3,51,66,198～202]

考虑一个由N个相同的节点构成的连续时间动态网络,其中第i个节点的状态方程

:

・

xi=f(xi)+c

j=1

6

N

aijH(xj)(7.1)

其中,xi=(xi1,xi2,…xin)∈Rn为节点i的状态变量;常数c>0为网络的耦合强度;

):Rn→Rn为各个节点状态变量之间的内部耦合函数,也称为各节点的输出函数,这里H(・

假设每个节点的输出函数是完全相同的;耦合矩阵A=(aij)∈RN×N表示网络的拓扑结构。

当耦合矩阵A描述了一个一般的拓扑结构时,具体定义为:若节点i和节点j(i≠j)之间有连接,则aij>0;否则aij=0(i≠j)。对角线上的元素为

aii=-j=1,j≠i

()()()

6

N

aij　　i=1,2,…,N(7.2)

假设网络是连通的,那么耦合矩阵A有(i)且仅有一个重数为1的零特征值,而A其余的特征值的实部均为负实数。

如果当t→∞时有

(7.3)x1(t)=x2(t)=…=xN(t)=s(t)

就称动态网络(7.1)达到完全(渐近)同步。由于耗散耦合条件,同步状态s(t)∈Rn必为单个

・

孤立节点的解,即满足s(t)=f(s(t))。这里s(t)可以是孤立节点的平衡点、周期轨道、甚至是混沌轨道,这样就有多种同步。7.1.1　同步判据

Pecora和Carroll研究了系统(7.1)线性耦合网络同步的稳定性问题,给出了主稳定函数

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

363

判据[66]。Motter等人在他们工作的基础上扩展网络的同步判据,研究网络的Lapalcae矩阵不能对角化时的同步[203,204]。

将系统(7.1)重写为

・

xi=f(xi)-c

j=1

6

N

lijH(xj)(7.4)

其中,L=(lij)=(-aij)是网络的Lapalcian矩阵。注意,这里的矩阵L不一定是对称的,因为在这里的网络并不局限于无向网络。

稳定性分析可以通过考虑等式(7.4)的同步解的变分而得到:

ξi=DF(s)ξi-c6LijDH(s)ξj

j=1

・

N

(7.5)

等式(7.5)同样可以写成矩阵形式

ξ=DF(s)ξ-cDH(s)ξL・

T

(7.6)

T

ξξ其中ξ=(在Pecora2Carroll的分析中,1,…,N)是一个N×N维的矩阵,L代表L的转置。

假设Lapalcian矩阵L是可以对角化的,从而对等式(7.5)进行对角化。这里并不假设L是可

以对角化的。相反,利用L的Jordan标准形转换。对于任意的N×N维的矩阵L,存在一个可逆阵P,可以将矩阵L转换为Jordan标准形

0

PLP=J=

-1

B1

…

B1(7.7)

其中Bi是块状的形式为

-λ

Bi=

1

-λ

…

…1

(7.8)

-λT

(P-1)和-λ是L的特征值(可能为复数)。在等式(7.6)中令η=ξ,可以得到

(7.9)

η=DF(s)η-cDH(s)ηJ

・

T

作为一个所有ξi的线性组合,η的每一列通常代表一种对整个网络的扰动,但不是对任何特定

的振子。因此,当且仅当在等式(7.9)中的所有的这些列都收敛于0,同步解才是稳定的。7.1.2　Lapalcian矩阵L可对角化

在讲述一般情况以前,先考虑一下L可以对角化的情况。在这种情况下,矩阵J是一个

λ对角阵,在对角线上的特征值依次为λ1,…N。于是对η的每一列,这个等式变为相互独立,并且取得形式为

・

y=[DF(s)-αDH(s)]y(7.10)

λi。其中D代表Jacobian运算符,当y代表η的第i列α=c把α看作一个可以调整的复数参数,

等式(7.10)被称作一个主稳定等式,它的稳定截面作为α的函数决定系统同步解的线性稳定

364物理学进展　　27卷

)通常用来检查稳定性,并被称作主稳定函性。对于y=0的解的最大Lyapunov指数Λ(αλλ数[201]。L的特征值可以写为0=λ-Re1<-Re2≤…≤N,则同步解线性稳定的条件为

Λ(cλi)<0,　　i=2,3,…,N(7.11)

下面考虑Laplacian矩阵的特征值均为实数的情况,并且可以排列为0=-λ-λ1≤2≤…≤-λ用SRΑR代表主稳定函数为负值的区域。根据动力学函数f,输出函数H和同步状态N。

。

λ类型I网络:SR=(-∞,α对于这种类型的网络,如果d1),其中-∞<α1<0。2<α1或

[205,206]

者d>|α,则同步状态是稳定的。这就意味着类型I网络的同步化能力可以用耦1/λ2|

合矩阵的第二大特征值来刻画。较小的λ2(绝对值大)导致较强的同步化能力。α类型II网络:SR=(α2,3),其中-∞<

αλN>对于这种类型的网络,如果d2<α3<0。

s的不同,存在着三类网络

[66,201]αλ,则2和d2<α3,也就是λN/λ2<α2/α3同步状态是稳定的。对于典型的振子有α2/α3>1。这就意味着类型II网络的同步化能力可以用耦合矩阵的特征值之比λ较小N/λ2来刻画。的λN/λ2导致较强的同步化能力。[204]

类型III网络:SR=空集。对于任意的耦合强度和外耦合矩阵,这种类型的网络是不可能达到同步的。

图7.1是这三类网络的主稳定函数的示意图。当内部耦合矩阵为对角阵H=diag(1,0,

图7.1　三种网络的主稳定函数示意图[11]0)时,一个典型的类型I的网络为:当节点的

[205]

动力学为Chua电路时,α;一个典型的类型II的网络为:当节点的动力学为1=-10

[196]Rossler振子时α,2/α。3≈37.85

7.1.3　Lapalcian矩阵L不可对角化

现在考虑L不一定可以对角化的情况[203,204]。Jordan标准形的每一块对应于在η中的列的组合,并且满足等式(7.9)的一个子集。例如,如果Bi是k×k维的,并且相应的η的列用η1,ηη2,…,k表示,则这些等式的形式为

ηηDH(s)]1=[DF(s)-α1

ηηDH(s)]2=[DF(s)-α2-cDH(s)η1

・・

・

(7.12)(7.13)

ηη(7.14)DH(s)]k=[DF(s)-αk-cDH(s)ηk-1

ηηλ相关的广义特征值空间中的扰动模式。这里η等式(7.12)与主稳定1,2…k代表了与特征值

λ)<0,η等式(7.9)有完全相似的形式,因此当t→∞且仅当Λ(c对于等式1指数收敛到0。

(7.13)稳定的条件显然牵扯更多,但是可以按照下面的方式确定。λ)<0和DH(s)是假设Λ(c

λ)<0保证第一项和第二项有界的,则等式(7.13)的第二项也是指数小的。则同样的条件Λ(c

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

365

λ)<0,的稳定效果,这就导致当t→∞,η重复运用上述的讨论表明:如果Λ(c2指数收敛到0。ηηλ)<0是对应于每个块Bi的等式的线性稳定性条件。这就说明Λ(c3,…,k也必然收敛到0。

当所有的Jordan块都被考虑,我们可以看到在一般的非对角化的条件下,同步解的稳定性条件同样也是由等式(7.11)给出。

尽管对于可以对角化的和非对角化的两种情况稳定性条件相同,但值得注意的是有一个本质区别。如果矩阵化L是可以对角的,于是每个扰动模式是与其他模式相解耦的,则η的每一列同时指数收敛和独立于其他列。另一方面,如果是L是不可以对角化的,一些扰动模式可能受到长期振荡的影响,因为它们与具有同一特征值的其他模式相耦合。在等式(7.12)2(7.14)中,η2在开始收敛前必须等待η1变得充分小,η3必须等待η2,依次类推,则ηk在开始收敛前必须等待很长时间。Jordan块的尺寸k越大,我们所预期的这个过渡时间越长。7.2　网络结构对网络同步化性能的影响

关于描述复杂网络基本结构的主要参数以及它们对网络同步化能力的影响,可以参考有关复杂网络同步研究的综述文章[21,22,195]。与具有规则的最近邻的耦合拓扑的网络相比,已经发现具有小世界拓扑结构的网络的同步化能力大大的增强[202,207]。这种同步化能力的提高最初人们认为是由于增加了少量的长程边而导致的网络节点之间的平均距离降低而获得的。Nishikawa等人[208]发现具有均匀度分布的无权网络比度分布不均匀的无权网络更加容易同步,即使前者的最短距离比后者的最短距离大。Hong等人[209]的研究表明网络的最大的中介数比越小,网络越容易达到同步。近几年来,加权网络的同步化也被广泛研究。例如,Motter等人[210]发现适当的加权无标度网络可以提高网络的同步化能力,Chavez等人[211]发现在加权无标度网络中在网络的路径的基础上的加权过程提高了同步化能力。7.2.1　群聚系数对同步性能的影响

群聚系数是用来表征网络中的一个节点的邻居同样也相互为邻居的趋势,它是另外一个影响网络同步化能力的重要参数[1]。McGraw等人[212]研究了群聚系数对无权网络的相位同步的影响。他们发现群聚系数可以提高连通度高的节点的同步化能力,即使它抑制了无标度网络的全局同步化能力。由于已经发现对一个网络进行加权可以提高网络的同步化能力,那很自然的要问这样一个问题:群聚系数是否抑制加权网络的同步化能力?

根据边交换的方法,Kim在不改变网络的度分布的条件下,提出了一个通过重连边机制来改变网络的群聚系数的模型[213]。这个算法的具体操作:随机选择两条边,分别连接节点A和B,以及节点C和D。这两条边相互交换端点,即原来的两条边A2B和C2D被修改为两条新边A2D和B2C。很明显,这个过程并不改变每个节点的度,只是改变了网络中的三角形的个数。边的交换只有在新的网络具有一个高(或者低)的群聚系数时才被接受。

我们研究了无权网络和加权网络的WS小世界网络的同步。具体做法是,首先生成一个WS小世界网络,然后通过交换边的方法增加网络的群聚系数。从图7.2和图7.3可以看出,随着网络群聚系数的增加,根据网络的同步化判据,网络的第二大特征值λ类网络的同步能2增加,第Ⅰ力可增强,但是由于特征值之比R=λ类网络的同步化性能将降低。N/λ2增加,因此,第Ⅱ

366物理学进展　　27卷

图7.2　在不同的重连概率下无权WS小世界网络的同步化能力。插图为在不同的群聚系数下,无

权WS小世界网络的同步化能力。数据是50次独立试验的平均结果,网络的规模为N=

[198]

500。(a)λ2vs.群聚系数C;(b)Rvs.C

图7.3　在不同的重连概率下,β=1时的加权WS小世界网络的同步化能力。插图为在不同的群聚系数下,β=1时的加权WS小世界网络的同步化能力。其他的条件与图7.2相同。

[3]

(a)λ2vs.群聚系数C;(b)Rvs.C

同样,我们对BA无标度网络和LC局域世界网络也作了同样的研究,发现不管对无权网

络和加权网络,群聚系数越大,度分布具有指数分布和幂律分布的网络的同步化性能越差。

吴翔等人研究了第二类网络在不同群聚系数情况下的同步化能力,发现网络的群聚系数越大,网络的同步化能力越差[213]。这一点与我们得到结论是相符合的。7.2.3　网络的异质性对网络的影响

我们知道,一个无权网络越异质,即网络的度分布越不均匀,网络越难以达到同步[207]。然而,Motter等人的研究表明,一个加权网络变得越异质,加权网络的同步化能力却没发生改变[209]。值得注意的是,在Motter等人研究的加权网络在网络生成过程中并没有改变边的权值。但是,我们通常所遇到网络边的权值会在网络的生长过程中会发生改变,如交通流网络中,随着新边的加入,网络中的原有边的权重会作相应的调整。

最近,Barrat,Barthelemy和Vespignani(BBV)提出了一个加权无标度网络模型[56]。它同样有两个生长步骤:(1)优先增长:一个节点i与新结点连接的概率为

∏

n→i

=si/

∑s

j

j

,其

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

367

中si=

的权值被初始赋为1。一个新边的加入只会触动节点i的邻居j∈v(i)的权值的局部变化:wij

γ2

→wij+Δwij,其中Δwij=δwij/si,δ为实数。BBV网络的度分布p(k)～k,其中p(k)为节δ+3)/(2δ+1)为幂率指数[56]。点的度为k的概率,γ=(4图7.4中显示了不同的δ时,BBV网络的度分布和幂率指数,网络的规模N=10000和m=m0=3。我们不难看出:随着δ的增加,网络的幂率指数降低。也就是说,随着δ的增加,网络的异质性增加。

6

j

(2)边的权值的调整:每条新边(n,i)wij,wij是连接节点i和节点j的边的权值。

图7.4　数据是在100次试验的平均结果,网络的规模N=10000和m=m0=3。(a)度分布p(k)vs.

k[198];(b)幂率指数γvs.δ

等式(7.1)的耦合阵子取为lij=wij或者lij=wij=wβij

v(i)是阵子i的邻居的集合。矩阵W=(wij)

N×N

则BBV模型所决定的加权为wij

wii=-j=1,j≠1

的定义如下:如果振子i和振子j存在一条边,

=wji(i≠j);否则,wij=wji=0(i≠j)。我们取

6

j∈v(t)

wij,其中β是个实数和

β

6

N

wij(i=1,2,…,N)。

如果等式(7.1)的右端的耦合因素为lij=wij,则矩阵G是一个对称矩阵,并且它的特征

值为满足0=-λ-λ-λ如果等式1≤2≤…≤N,此时我们说这个加权网络是对称耦合的。

(7.1)的右端的耦合因素为lij=wβij/

们说网络是非对称耦合的。然而,在这种情况下,矩阵G的特征值仍然满足不等式0=-λ-1≤

[209,211]

λ-λ。2≤…≤N

图7.5　研究了对称耦合的加权网络在不同网络规模下的同步能力的变化。从图7.5(a)中可见,对于同一个网络规模,我们可以发现第二大特征值λ由于δ的2随着δ的增加而降低。

增加,网络变得更加异质,这就意味着类型I对称耦合加权网络的同步化能力随着网络变得异质而增强。然而,从图7.5(b)中可以看出,对于同一个网络规模,耦合矩阵的特征值之比R随着δ的增加而迅速增加。因此,网络变得非异质,这就意味着类型II对称耦合加权网络随着网络变得异质而不容易达到同步。

下面我们研究加权网络为非对称耦合时,网络的同步化能力的变化。图7.6描述了类型I非对称耦合的加权网络的同步化能力为β的函数,网络是按照BBV模型生成的,网络规模N=1000,m=m0=3和δ=5。可以观察到β≈1时第二大特征值λ这就意味着当β≈2最小。1,类型Ⅰ非对称耦合的加权网络的同步化能力最强。于是,我们着重研究一下在β=1时,其他

6

j∈v(i)

wij,则矩阵G对所有的实数β是非对称,此时我

β

368物理学进展　　27卷

图7.5　在不同的网络规模N和m=m0=3下,对称耦合加权网络的同步化能力[198]。(a)λ2vs.δ

;(b)R=λN/λ2vs.δ

参数的变化对类型I非对称耦合的加权网络的同步化能力的影响。图7.7(a)描述了类型I非对称耦合非加权网络在不同的网络规模下的同步化能力随着δ的变化情况,参数m=m0=3和β=1。对于同一个网络规模,我们可以发现第二大特征值λ2随着δ的增加而降低(绝对值变大)。因此,类型I非对称耦合加权网络的同步化能力随着网络变得异质而增加。很明显在图7.7(b)中,对于同一个m值,λ因此,随着δ增加,类2随着δ的增加而降低。

型I非对称耦合加权网络的同步化能力提高。

[198]图7.8描述了类型II非对称耦合的加权图7.6　第二大特征值λ2vs.β。数据是100

网络的同步化能力为β的函数,网络是按照次试验的平均值,并且N=1000,m=

m0=3和δ=5

图7.7　类型I非对称耦合加权网络的同步化能力[198]。(a)在不同的网络规模N和m=m0=3时,λ2

vs.δ;(b)在不同的m值和N=1000时,λ所有的数据是50次数据的平均值2vs.δ。

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

369

BBV模型生成的,网络规模N=1000,m=m0=3和δ=5。可以观察到β≈1时,特征值之比R最小。这一点与参

考文献[204,211]的结论一致。接下来,我们着重研究一下在β=1时,其他参数的变化对类型II非对称耦合的加权网络的同步化能力的影响。

图7.9(a)描述了类型II非对称耦合加权网络的同步化能力。对于同一个网络规模,耦合矩阵的特征值之比R随着δ的增加而迅速减小。在图7.9(b)中,对于同一个m

值,R随着δ的增加而降低。因此,随着δ的增加,网络的异

质性减弱,类型II非对称耦合加权网络的同步化能力提高。图7.8　特征值之比[198]

R=λN/λ2vs.

7.3　确定性因素对网络同步的影响

正如上章所述,无权的BA无标度模型以及有权的β。数据是100次试验

的平均值,并且N=1000,m=m0=3andδ=5

图7.9　类型II非对称耦合加权网络的同步化能力[198]。(a)在不同的网络规模N和m=m0=3时,R

=λ;(b)在不同的m值和N=1000时R=λ。所有的数据是50次数据的N/λ2vs.δN/λ2vs.δ

平均值

BBV和TDE模型,都不能完全描述很多实际网络的性质

[20]

.因为真实世界是一个随机因素

和确定因素和谐的统一的世界,随机因素是否起主导作用要视具体网络情况而定。具有复杂相互作用的确定性因素在实际的网络中有时起到决定性的作用,特别是对于非线性动力学和混沌动态复杂网络。

我们用对无权网络混合优先连接机制代替BBV模型中的随机优先连接机制[50]。这个新的模型叫做混合优先BBV模型[199],就是上章提出的和谐统一的混合择优模型(HUHPM),把这个模型的思想方法应用于BBV模型,称为HUHPM2BBV网络。它的生成算法为:

(i)生长:初始的网络是连通的,它有m0个节点和较少的e0条边组成,并且边的权重为1。在每一时步,增加一个新的节点和m(m≤m0)条边相连。这些边分别连接新节点和m个原先存在于网络中的节点。新边的权重为1。

(ii)混合优先连接:一个混合比的定义为

d/r=确定性择优连接(DPA)的时间步数与随机性择优连接(RPA)的时间步数之比,即

370物理学进展　　27卷

d/r=DPA/RPA,混合比d/r取值范围为0到无穷大。

(a)随机优先连接(RPA):对于加权HUHPM2BBV,网络中原有的一个节点i与新节点连

接的概率为

pi(si)=si/

其中si=权值变为:

6

6

j

sj(7.15)

j

wij是节点i的强度,wij是连接节点i和节点j的边的权值。与节点i相连的边的

wij→wij+Δwij

(7.16)(7.17)

其中Δwij依赖于局部动力学,并且是权值wij的一个函数:

Δwij=δwij/si

　　(b)确定性优先连接(DPA):在每次连接后,整个网络按照节点度从大到小进行排序:k1≥k2≥…≥kN。然后,选择m个度最大的不同节点与新节点连接。经过d时步的确定性优先连接(DPA)以及经过上面r时步后,这个过程将生成一个网络规模为N=r+d+m0的网络。依次重复以上算法中的步骤(a)和(b)直到生成我们所需要的网络。我们分析发现了HUHPM2BBV网络同时具有无标度和小世界的性质。随着混合比的增加,网络度分布的幂率指数和网络的群聚系数增加;而网络的平均最短路径减小。

这些新出现的网络的拓扑性质如何影响一个动态网络的动力学性质?图7.10描述了类型I的HUHPM网络的同步化能力。对于无权模型(δ=0),在图7.10(a)中,当混合比d/r<1/1,1/|λ1/1之后,1/|λ2|逐步上升,但是随d/r≥2|迅速降低。因此,随着混合比的增加,类型I的无权HUHPM网络的同步化能力首先慢慢减弱,然后同步化能力增强。

图7.10(a)还表明对于HUHPM2BBV加权模型(δ>0),随着混合比d/r的增加,1/|λ因此,随着混合比的增加,类型I的HUHPM2BBV网2|首先缓慢减弱,然后快速增加。

络随着δ的增加而变得容易同步。

图7.10(b)描述了类型I的HUHPM2BBV网络在不同的混合比下,1/|λ2|随着网络规模N的变化情况。很明显,对于同一个混合比d/r/,1/|λ2|随着网络规模N的增加而增加。当网络规模从200增加到5000时,在d/r=1/1,1/|λ2|从0.3890增加到0.4287;随着d/r继续增加,1/|λ因此,类型I的HUHPM2BBV网络随着网络规模的增加变得2|缓慢增加。

难以同步。我们同样注意到对于相同的网络规模,随着混合比d/r的增加,1/|λ2|首先增加然后下降。因此,随着混合比的增加,类型I的HUHPM2BBV网络的同步化能力先减弱后增强。

图7.11　描述了类型II的HUHPM2BBV网络的同步化能力。从图7.11(a)中可以看出,对于无权模型(δ=0),当混合比从1/199增加到4/1,R从176.9增加到1300.9;但是当混合比从4/1增加到199/1时,R从1300.9减少到1177.2。这就意味着随着混合比的增加,类型II的HUHPM2BBV网络的同步化能力也是先减弱后增强。

从图7.11(a)中还看出,对于加权模型和相同的δ值,随着混合比d/r的增加,R的值先增加后减小。因此,随着混合比的增加,类型II的HUHPM2BBV网络的同步化能力先减弱后增强。我们从图7.11(a)也可以观察到,对于相同的混合比,R随着δ的增加而增加。对于d/r=1/1,当δ从0增加到5,R从1113增加到9094。这意味着类型II的HUHPM2BBV网

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

371

图7.10　类型I的HUHPM2BBV网络的同步化能力[199]。(a)1/|λ2|vs.混合比d/r在不同的δ值,

N=2000和m=m0=3下。(b)1/|λ2|vs.N在不同的混合比d/r,δ=2和m=m0=

3下。数据为50次结果的平均值

图7.11　类型II的HUHPM2BBV网络的同步化能力[199]。(a)R=λN/λ2vs.混合比d/r在不同的δ值,

(b)Rvs.N在不同的混合比d/r,δ=2和m=m0=3下。数N=2000和m=m0=3下。

据为50次结果的平均值

络随着δ的增加而变得难以达到同步。

在图7.11(b)中,对于相同的混合比d/r,特征值之比R随着网络规模的增加而增加。当d/r=1/1,当网络规模N从200增加到5000,R从300增加到8370。这就意味着网络规模的增加,类型II的HUHPM2BBV网络变得难以达到同步。从图7.11(b)还看出,对于相同的网络规模,随着混合比d/r增加,R先增加后减小。因此,随着混合比的增加,类型II的HUHPM2BBV网络的同步化能力先减弱后增强。7.4　如何构建同步化性能高的网络

目前大部分对网络同步化性能提高的改进模型,都是在基于一个已有模型的基础上,通过给网络加入少量的扰动边或者调整网络中边的权值从而提高网络的同步化性能。我们知道,第一类网络的同步化能力与网络耦合矩阵的第二大特征值密切相关,那么,这个第二大特征值对网络的生成过程的影响如何?

372物理学进展　　27卷

7.4.1　同步最优网络模型

我组提出了一个增长的、并且同步化性能有所提高的网络同步最优网络模型,该模型的生成规则如下[200]:

①增长:从一个具有m0个节点的网络开始,每次引入一个新的节点并将其连接到m(≤

m0)个已存在的节点上。

②同步最优连接:新节点与已经存在的节点i相连接时要使得构成的新网络的同步化性能最优,即要使得新的网络耦合矩阵的第二大特征值最小。

在经过t步后,产生一个有N=t+m0个节点、mt条边的同步化性能最优的网络模型。需要注意的是,在同步最优网络生成的过程中,只是在每条新边加入时网络的同步化性能达到最优,并不能保证最终得到的整个网络的同步特性达到全局最优。图7.12　m0=m=5,30个节点的同步最优网络拓扑结构示意图[200]

直接由同步最优网络模型的拓扑图(见7.12)可以很容易的发现,同步最优网络模型拓扑结构类似于多中心网络。在多中心网络模型中,有m个hub节点,其余N-m个节点均与这m个hub节点相连。这使得同步最优网络在恶意攻击下要比BA模型更加容易崩溃。7.4.2　同步优先网络

在同步最优网络的基础上,结合同步优先机制,进一步提出了同步优先网络模型[201]。构造算法描述如下:

①增长:从一个具有m0节点的网络开始,每次引入一个新的节点并将其连接到(m≤m0)个已存在的节点上;

②同步优先连接:新节点与已经存在的节点i相连接的概率与得到的网络的同步化性能有如下关系:

这样经过tµm0时间步后,得到一个具有N=m0+t个节点、mt条边的同步优先网络模型。

同步优先网络与BA无标度网络模型的区别在于,同步优先网络中的优先连接概率与得

∏=λ/6λ

i

22j

(7.18)

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

373

到的网络的同步化性能成正比,而不是与节点的度成正比。这样,在同步优先网络演化过程中,每加入一条新边后,同步化性能高的网络被优先生成得到。BA无标度网络的连接度分布非常不均匀:大部分节点有少量的连边而有小部分节点却具有大量的连边。对于生成的同步最优网络,研究发现该网络模型具有类多中心结构。而同步优先网络模型的拓扑结构比较复杂,它的连接度分布既不服从幂律分布又不服从指数分布。7.4.3　同步最优和同步优先网络的同步化性能

令初始网络的节点个数m0=m;同步最优网络和同步优先网络的耦合矩阵分别为Aso(m,N)和Asp(m,N),λ2so(m,N)和λ2sp(m,N)为Aso(m,N)和Asp(m,N)第二

大特征值。

图7.13给出了m0=m=3,BA无标度网络(带有圆圈的实线)、同步最优网络(带有菱形的虚线)和同步优先网络(带有方形的点　图7.13　当

划线)耦合矩阵的第二大特征值随时间变化曲线。图7.13中所有数据均作了5次平均得到。随着网络中节点个数的增多,三个网络

同步最优m0=m=3时,无标度网络、

网络以及同步优先网络耦合矩阵的第二大

[201]

特征值λ2(m,N)变化趋势图

的第二大特征值都趋于一个负常数。其中同步最优网络模型的同步化性能最强,同步优先网络模型的同步化能力类似于无标度网络模型,只是前者稍稍逊色于后者。当m0=m=5,7时,得到与m0=m=3相似的结论,只不过随着m的增大,网络耦合矩阵的第二大特征值都相应地减小。表7.1给列出了具有500个节点的各网络耦合矩阵第二大特征值。由表中数据可以发现,同步最优网络的同步化性能要优于BA无标度网络和同步优先网络的同步特性[200]。

表7.1　BA、同步最优网络、多中心网络以及同步优先网络模型的耦合矩阵的第二大特征值对比表[201]

网络模型

BA模型

m0=3

m0=5

m0=7

-1.2329-1.9898-3-1.2493

-2.8758-3.9635-5-2.8765

-4.6110-5.8710-7-4.5982

同步最优模型多中心模型同步优先模型

7.4.4　同步最优和同步优先网络同步的鲁棒性和脆弱性

由前面分析已知,同步最优网络拓扑结构类似于多中心网络,当随机去除网络中的部分节

点时,被去除的大多数是度非常小的节点,因此网络的同步化性能基本保持不变(见图7.14中带菱形的点线)。同步优先网络的度分布虽然不严格服从幂律分布,但网络中度的分布也不是均匀的,同样存在大量度小的节点,当我们随机去除同步优先网络中5%的节点后,网络耦合

374物理学进展　　27卷

矩阵的第二大特征值从-1.2493增大为-0.800(图7.14种带方形的点划线),即使再随机去除更多的节点,耦合矩阵第二大特征值的变化幅度也不会太大。图7.14比较了随机去除BA无标度网络、ER随机图、同步最优和同步优先网络中5%的节点后,网络耦合矩阵第二大特征值的变化情况。所以,同步最优和同步优先网络都对抵抗随机故障具有较强的鲁棒性。

同步最优网络中有少数几个度非

同步对随机故障的鲁棒性[200](符号意义同图7.13)常大的hub节点,如果恶意去除这几图7.14　

个(略大于m个)节点,整个网络就会崩溃为一个个孤立的节点。对于同步优先网络,当5%度最大的节点被去除后,第二大特征值从-1.2493变为-0.3802,网络没有出现孤立节点或孤立簇,而BA无标度网络在去除不到1%的最大度节点后,网络中就出现不连通子图,整个网络无法达到同步。如果继续研究最大连通子图的同步化性能会发现,当5%度最大的节点被去除后,同步优先网络耦合矩阵的第二大特征值上升幅值仍小于BA无标度网络中最大连通子图的第二大特征值的上升幅值(见图7.15),因此说明同步优先网络要比BA无标度网络对于恶意攻击更具鲁棒性。也就是说,我们所提出的同步优先网络克服了BA无标度网络的“A2chilles踵(Achilles’Heel)”。我们都知道随机网络的度分布是均质的,去除任意一个节点,都会对网络产生近似的影

响,因此同步优先动态网络具有较强的同步鲁棒性进一步说明了同步优先连接机制要比BA模型的连接度优先连

图7.15　同步对恶意攻击的脆弱性[200]

接带来更均质的连接分布。图7.15中

的插图表示的是去除掉百分之f连接度最大的节点后网络中最大连通簇在整个网络中所占的比例。由图中可以看出当移除同步优先网络中5%度最大的节点后,网络依然保持连通,由于BA无标度网络的度分布要比随机图的异质,因此去除5%度最大节点后BA网络中的最大连通子图要小于随机图中的最大连通子图。7.5　有向网络中的同步

Nishikawa和Motter在文献[203]中指出一个给定的网络满足以下条件:1.网络中可以

嵌入一个有向生成树。2.网络中没有反馈环.3.归一化每个节点的输入信号时,网络耦合矩阵的特征根除有一个为0外,其余的都相同。因此,根据第二类网络的特征值判据得到此时网

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

375

络的特征值之比R=λN/λ2=1,此时第二类网络的同步化性能达到最优。我们说一个网络包含一个生成树是指如果存在一个节点,叫做根节点,并且对网络中每个其他节点j,必须至少存在一条有向路径从根节点指向节点j。这里的网络的同步化准则是按照Lyapunov指数来确定,它是个局部的结果。Lu和Chen研究了有向网络中非线性振子耦合的连续系统的全局同步[214],指出如果网络中包含一个有向生成树,对于充分大的耦合强度,网络可以达到同步。这就意味着如果有一个动力系统直接或者间接地影响其它所有的动力系统,对于足够强的耦合强度,同步是可能的。卢文联和陈天平等人将上述有向网络的连续系统同步判据推广到了由线性耦合映射格子(LCMLs)组成的有向网络离散系统中[215],我们采取原文中的符号描述系统的同步方程:

(t))+εx(t+1)=f(x′

j

j=1

6

m

j

lijf(x(t)),　　i=1,2,…,m(7.19)

其中x(t)=(x1(t),…,x(t))

ii

in

T

∈R是第i个节点的状态变量,t是离散时间,ε是耦合强度,

n

耦合矩阵L的非对角线元素非负,且每行元素和为0。他们指出如果耦合强度满足以下关系[214]-ρi-

1)1)2222

ρρ-ρi-βi(1-i+i-βi(1-22

κκ<ε<22

ββii

(7.20)

βκ的设置参见文献[214],则该所对应的图含有一个有向的生成树,网络可以达到同步。其中ρi,i,下面的例子用来说明LCMLs的同步是通过一个有向的可约图并且Laplacian矩阵:

01

L=

0

-1

00

-1

000

-1

(7.21)

10

012

12

　　耦合映射是一个二维的神经网络,其方程为[214]:

-1-0.5

(7.22)f(x)=kx+Δttanh(125x)+0.8

-0.5-2

其中Δt=0.02,k=1-Δ和t=0.98。根据不等式(7.20),如果(6/7)<ε<(8/7),则耦合系统(7.19)可以达到同步,如图7.13所示。

此外,在不改变网络的节点的度的情况下,我组提出了一个利用重连边方法构建的加权均质网络模型。发现随着均质网络中边的重连概率的增加,网络变得更容易同步。对于相同的边重连概率,随着网络规模的增加,加权均质网络的同步化能力降低,但是无权的NW小世界网络的同步化能力加强。另外,提出了一种新的提高网络同步能力的方法,它适用于网络的节点的介数与度不成正比例关系的情形,研究了此方法对网络其他结构特征的影响。结果表明,该方法可以有效提高此类网络的同步能力,同时增大网络的最短路径,降低网络的聚类性。我们观察到:边的介数的最大值可以影响网络的同步能力,但它不是一个可以独立决定网络同步能力的因素。

376物理学进展　　27卷

图7.16　同步方差与耦合强度的关系[215]。初始值在[-5,5]之间随机选取,同步用方差度

量为:

var=〈(1/(m-1))

6mi=1

2〈・〉代表时间平均‖xj(t)-x(t)‖,,其中2〉7.6　小结

本章首先介绍了完全同步判据等基本概念,然后讨论和总结了我们研究的网络同步的主

要进展,包括网络的结构特性对无权网络和加权网络的影响,以及如何构建同步化能力较强的网络和有向网络中的同步。然而,复杂网络同步方面的研究还有许多问题有待解决:

1.目前虽然已有的大量工作研究了网络的拓扑结构与同步化性能之间的关系,但还缺乏针对不同的网络类型揭示网络同步化能力的拓扑特性的特征量,这是一个很有意义的研究方向,我们将在第10章中进一步讨论。

2.网络的形成与演化机制影响着网络的结构,从而决定了网络的动力学特性,包括同步化能力。今后需要进一步考虑建立更加符合现实网络的演化机制,特别是基于局部信息的演化机制,寻求有效提高或减弱网络同步化性能以及鲁棒性的方法。现有的许多研究局限于网络拓扑固定的情形。关于时变拓扑的研究仍然有待进一步深入。

3.将复杂网络同步化研究与动态网络中的其它相关研究结合起来也许会能有效推动这些相关问题的研究。例如,一致性问题(consensusproblem)、群聚(swarming)、蜂拥(flocking)和网络拥塞等问题都与网络同步密切相关。这类课题综合研究必将有助于将同步化理论应用于实际工程领域。

8　耦合动力学网络的部分同步:现象和稳定性分析

8.1　引言

同步是耦合动力学系统中最简单的合作行为之一,同步现象在物理、化学、生物等领域无所不在[215～223]。耦合混沌振子的同步行为由于与人们对混沌认识的直觉不同而引起了广泛的关注(如众所周知的蝴蝶效应,在混沌轨道上小的扰动即造成了指数分离)。自从Pecora和

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

377

Carroll发现了混沌系统的完全同步[215],人们又相继发现了广义同步[221,224～228]、相同步[221]和

部分同步[222]等多种同步现象并进行了深入的研究。由于实际时空系统的建模需要,耦合混沌系统的同步现象引起了人们广泛的重视。

近年来,复杂网络引起了人们广泛的研究兴趣[1,3,19,20],主要关注各种网络的拓扑结构以及拓扑结构对网络动力学的作用。振子互相耦合形成网络,例如小世界网络、无标度网络等,多种描述拓扑结构的几何量,如平均最短路径、群聚系数、度分布等等,都影响网络的同步能力[202,211,229～237]。然而,由于网络的复杂性,影响耦合混沌振子网络的决定性因素还不清楚,引发了许多争论。

关于复杂网络的全局同步有大量的研究。本章更感兴趣的是在全局同步之前的振子合作行为。完全同步之前的分岔非常重要,对此的研究将为斑图、输运和波等许多合作行为带来新的契机。为了简化研究,对仅有几个非局域连接的规则网络的研究是非常有指导意义的。人们已经发现长程连接对网络的同步能力起着非常重要的作用[238]。为此,我们加入一些长程连接,不同方式的非局域连接将影响系统的合作行为,从而发现了有趣的部分同步现象。

近来,全局耦合混沌系统中的部分同步现象已被广泛研究,但是空间结构并没有被考虑[239]。局域耦合系统中的部分同步现象也已被发现[240]。通过对网络拓扑结构的对称性分析,人们讨论了系统的部分同步解[241]。迄今为止,大部分研究结果都依赖于数值实验。事实上,关于完全同步的理论已非常清楚,同步流形的稳定性取决于最大横向条件李指数。那么提出部分同步的判据是非常迫切的。有了这个判据,可以更好地研究网络拓扑对网络同步的影响。

本章探讨了在耦合混沌振子网络中的部分同步斑图现象,揭示了不同网络中的多种部分同步斑图[242～246]。部分同步斑图的存在依赖于网络拓扑的对称性质,同时部分同步现象的发生又依赖于相应解的稳定性。基于对部分斑图流形和横向流形的李指数谱的比较,本章给出了一种判定给定网络的部分同步斑图稳定性的判据。并从判据的角度研究了网络多种对称性的存在而导致的部分同步斑图现象的竞争和选择,发现对于一些有长程连接的拓扑上对称的环状网络,它的部分同步斑图状态可能是不稳定的。我们给出了表示同步和异步状态的相图,揭示了部分同步斑图系统的李指数谱的表现形式。研究表明:长程连接的增加破坏了环状网络李指数谱的简并性质,这可以很好的解释为网络拓扑对称性的破缺。目前的研究揭示了在达到全局同步之前的耦合动力系统固有的合作分叉行为。我们发现的现象和所作的讨论在揭示时空系统的斑图动力学方面具有理论意义和应用潜力。8.2　部分同步:现象

当若干给定的n维非线性系统

・

x=f(X)　(X=(x1,x2,…,xn))(8.1)

被耦合成一个更大的系统时,我们可以把那些非线性系统看作是一个网络的节点,而那些线性耦合就是这个网络的边,这样耦合网络系统的方程就可以写为

・

Γ󰃁CX　　　　　　　　　　　　　　(8.2)　　　　　　　　　　　　X=F(X)+ε

ε表示连接的耦合强　　X=(x1,x2,…,xn),F=(f1,f2,…,fN),N是网络中的结点数。

度,而Γ:Rm→Rm则是两个系统之间的连接方式。C=M-D,M是网络的邻接矩阵(矩阵元

Mij表示在网络中第i个结点和第j个结点之间的边数),D为一个对角的矩阵,且满足Dii=

378物理学进

C

展　

C

　27卷

6

N

Mij,由此易知j=1

λλ在本章中我们约定矩阵的本征值表示为:λ1≥2≥…≥N

1对于网络中的第i个与第j个振子的轨道,dij=lim‖xi(t)-xj(t)‖dt表示它们

T→∞

6

Nj=1

Cij=0,它的最大本征值λ1=0且λi<0(i=2,3,…,N)。

T

∫的距离。当dij→0表示这两个振子的轨道正收敛到同一个吸引子上,在这个意义上,认为这两个振子达到了同步。

我们以环状的网络上加一条非局域连接为例子来说明网络中两个振子之间的同步过程。取一个有10个节点的环状网络,并在其中的第1节点与第4节点间加入非局域连接(记作(1,4)),

・・・

(y-x),y=γ如8.1(a)所示。以Lorenz振子(x=σx-y-xz,z=xy-βz)为网络中的节点。图8.1　(a)具有一条非局域连接的10个节点的网络图;(b)振子间平均距离与耦合强度的关系

图8.1(b)给出了各对节点(振子)间的平均距离dij随耦合强度ε的变化。这里采用

σ=10,γ=27,β=8/3,而连接方式

0

000

000Γ=1

0

由图8.1可见,随着耦合强度的增加,有些振子之间先达到同步(如图8.1(a)中的d14),有些则后达到同步(如图8.1(a)中的d12),不妨把这种在有的耦合强度下只有部分振子之间发生同步的现象称之为网络中的“部分同步”。图8.2(a)(c)分别采用了Rossler振子和Logistic振子为节点,由图8.2(b)(d)可见,部分同步同样可以被观察到。这初步证明了部分同步这一现象的出现有相当的广泛性和一般性。

从图8.1中还不难发现,部分同步不但发生在那些直接连接在一起的振子之间(如图8.1(b)中的d23和d14),而且也发生在那些没有直接连接的振子之间(如图8.1(a)中的d5,10和

d69),而有的直接连接在一起的振子却不一定会同步(如图8.1(a)中的d56和d67),甚至,各

个振子对之间的部分同步都是在同一耦合强度(记为εPaS)下发生的,即上述几个距离是“同时”变零的。因而,就可以用如图8.1(a)这样的拓扑结构图来表示网络在适当耦合强度(ε>ε其中线代表耦合,其他符号代表各个振子,而相同的符号代表相互同步PaS)下的同步情况。

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

379

的各个振子。

图8.2　(a)具有一条非局域连接的6个节点的网络图;(b)节点动力学为Rossler振子时平均距

离与耦合强度的关系;(c)一条非局域连接的4个节点的网络图;(d)对Logistic映射情况平均距离与耦合强度的关系在图8.1(a)中,整个网络被分成了N/2个小集团,每个集团中有两个振子。而当网络中的非局域连接发生变化的时候,网络同步后的集团数就不一定再是振子数的一半了。图8.3给出了更多的可以发生部分同步的网络(这里我们仍然采用Lorenz振子为网络的节点),它们

(c)的网络有着不同的振子数和连接方式。图8.3(a)为N=10,(1,都与图8.1(a)及8.2(a)、

5)的情况,同步的集团数为6个,而且有两个集团只有一个振子。比较两种情况不难发现,这种同步总是发生在网络中镜像对称的两个振子之间,而在后一种情况下第3和第8个振子正好处在对称轴上,所以它们就只好和自己“同步”了。对于不同振子数和其他的连接方式也得到了同样的结论(如图8.3(b)(c)(d))。

图8.3　几种可以达到部分同步的对称网络

数值实验显示,部分同步总能够在由一条非局域连接的环状网络(记作Ns=1)中实现,而这种网络中又总是有一个镜像的对称性。但是,对于Ns=2的网络以及其他更为一般的网

络,情况就复杂了。例如图8.4(a)的情况就没有发现有部分同步的情况发生。由此可见存在一个镜像的对称性很可能是网络同步的一个必要条件。实际上,在这个网络(动力系统)中,要

380物理学进展　　27卷

发生某种同步就需要每个同步集团中的振子都有相同的动力学(方程的形式),而一个没有对称性的网络是不能做到这一点的。相反的,对于图8.4(b)(c)(d)的情况,它们都有一条镜像的对称轴,部分同步就是可以发生的,且发生的规律与前述的情况是相似的。

图8.4　(a)不满足镜像对称的具有两条非局域连接的网络;(b)平行网络;(c)交叉网络;(d)八字网络但是,上述命题的逆命题是否成立呢?即存在至少一个对称性是不是发生部分同步的充分条件呢?答案是否定的,图8.5(a)(b)的网络各有一个对称轴,但是在这些网络中是观察不到部分同步的。更进一步讲,如果系统中有多个对称性,那么情况又会怎样呢?这时部分同步会发生吗?如果能够,它又会沿着哪一条对称轴呢?例如网络N=8,Ns=2,(1,5),(3,7)(图8.5(c))就有四个镜像的对称性及三个旋转对称性,但是在这种网络中是观察不到部分同步的。由此可见,网络部分同步的发生条件不仅仅包括其对称性,它必然包含更多的关于网络拓扑结构及格点上动力学的信息。

图8.5　几种两条对称非局域连接但无法达到部分同步的网络

8.3　部分同步发生的判据

实际上,从动力学系统的角度来看,部分同步就是系统方程(8.2)的一个解,不妨表为:

x

k

l

=x

2

k

=…=x

kNk

,　　(k=1,2,…,Nc)(8.4)

这里Nc为发生部分同步时的集团数,Nk为第k个同步集团内的振子数,且满足

k=1

6

Nc

Nk=N。

这些等式确定了相空间上的一个mNc维的子流型,同前面讨论完全同步时的情况类似,初值在这个流型上的轨道将始终在其上运行,因而它是个不变子流型。而一定的对称性只能保证这个解或者说这个流型在系统中的存在性,但是它不能保证这个解(部分同步子流型)的稳定性。

A.完全同步:主方程

要弄清楚部分同步的发生条件,先关注一下完全同步的发生条件是有启发性的[229]。对于任意的一个高维的动力学系统(8.2)如果存在一个完全同步态:

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

12

x(t)=x(t)=…=x

N

381

(t)=s(t)(8.5)

　　那么可以在同步态s(t)附近线性化方程,即令线性化xi=s+ui(i=1,2,…,N),可得:

Γ󰃁CU(t)(8.6)U(t)=IN󰃁Df(s)+ε

这里U=(u1,u2,…,uN)而Df为f的雅可比矩阵在s(t)附近的值。注意到方程(8.6)的第一项是一个分块对角化的矩阵,不妨做坐标变换将C对角化。于是方程(8.6)就成为:

・

λkΓvk,(k=1,2,…,N)(8.7)vk=Df+ε

　　至此,方程(8.6)就被化成了N个独立的方程。这时因为λ1=0使得(8.7)中k=1的方程对应于系统在同步流型s(t)上的动力学,而其他的方程(k=2,3,…,N)就表示垂直于该流型的方向上(横向)的动力学。这样就把系统的同步流型和横向系统的动力学分开了,而后者的收敛就表示了系统正在趋于同步态。要判断在网络中完全同步态是否稳定(是否这样的同步现象能够发生),只需找到使方程

Γv　　　　　　　　　　　　　　(8.8)　　　　　　　　　　　　　　　v=Df+α

λk}(k=2,…,N)都落入了这个区域的最大李指数小于零的α平面上的区域。然后如果{ε中,则说明了横向的系统是稳定的,因而完全同步是可以发生在这个网络中的。

B.网络的邻接矩阵

事实上,上面的主方程方法可以用到对部分同步稳定性的分析中。对部分同步的分析,我们可以把网络系统分为两个部分,横向子系统的稳定性决定了部分同步的发生。这里讨论的网络系统满足镜像对称。

如果一个网络的邻接矩阵C在置换变换T下是不变的,即满足

TCT

-1

=C(8.9)

　　那么称这个网络有对称性T。具体的来说,例如图8.3(a)中的网络具有镜像的对称性,T矩阵可表为

T=

F4O6×4

O4×6F6

(8.10)

这里,Fk为k阶反对角矩阵,即满足Fi,K-i+1=1,i=1,2,…,K,且Fij=0对于i+j≠K/2

+1,OK×对称性变换的具体形式因网络的编号方式不同而有所不同,但L为K×L阶零矩阵。

这不改变问题的本质。另一个例子是图8.3(c)的网络,它有两个对称性分别为

01000000101000000001

(8.11)T1=00100,　T2=00100

00

00

00

01

1010

01

00

00

00　　为了使C分块对角化,考虑其对称变换矩阵T的本征矢矩阵S。借助本征矢矩阵S,我

们可以将矩阵T及其逆矩阵(如果存在)对角化,于是有:

(8.12)TS=SΛ或

ST

-1

S

-1

=Λ

-1

(8.13)

这里Λ为对角矩阵。若记

382物理学进

-1

展　　27卷

M=SCS(8.14)

-1

则有:

M=STS

-1

SCS

-1

STS

-1

(8.15)(8.16)

于是

M=ΛMΛ

-1

这样,若设(将M与Λ各分为四块)

M=

acac

-1

bdbd,Λ=

α00β0

(8.17)

则

ΛMΛ

-1

α0=

0β-1α0

-1β=

α

-1βca

βba

-1β(8.18)

将方程(8.18)带入方程(8.16)可得

βba

要使上式成立,只能

b=c=0

(8.20)(8.21)

β=b,ca-1

=c(8.19)

或

α=β

但是Λ作为T的本征值,不一定总能满足方程(8.21)。这意味着M是分块对角化的,即网络邻

接矩阵C的对称变换矩阵T的本征矢矩阵S可以使C分块对角化,或写作

A0-1

(8.22)G=SCS=

0B易知这个变换是一个相似变换,因而这个变换不会改变原来矩阵的本征值,而只是将它们重新

imtv

分配到了两个块A和B中。我们记{λi}(i=1,2,3,…,NA)及{λi}(i=1,2,3,…,NB)分别为A、B两块的本征值,NA、NB为这两个块的维数,且满足NA+NB=N。另一方面,由于两块中

im

会有一个(不妨设为A)有原来矩阵的0本征值,所以λ=0。1

C.部分同步发生的判据

运用上面的变换,

　　　　　　　　　　　　　　　　W=(Γ󰃁S)X　　　　　　　　　　　　　(8.23)考虑到我们只关心部分同步流形附近的动力学,方程(8.2)可以分成两部分:

・

Γ󰃁AWim(8.24)　　　　　　　　　　　　　Wim=FS(Wim)+ε

Γ󰃁BWtv(8.25)Wtv=Fs(W)+ε其中W=(w1,w2,…,wN)=(W

im

,W

tv

),W

im

12

=(w,w,…,w

NA

)表示同步流形上的动

力学,Wtv=(wNA+1,wNA+2,…,wN)表示横向流形上的动力学。

对方程(8.24)在同步流行附近线性化,得到:

(8.26)　　　　　　　　　　　　　w1=w2=…=wNA=w

i

对同步态(8.26)加以微扰,wi(t)=w(t)+δw(t)　(i=1,2,…,NA),我们可以得到

δΓ󰃁AδWim(8.27)Wim=DFs(Wim)+ε

　　那么通过对角化矩阵A,方程(8.27)可以分解成NA个独立的方程:

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

・

383

vk=

λkΓvk,(k=1,2,…,NA)DFs+ε

im(8.28)

　　写成一般的形式

Γv(8.29)　　　　　　　　　　　　　　　　v=DFs+α　　注意到方程(8.27)是一个独立的系统。要讨论它的完全同步只需要像前面方程(8.8)的λim讨论就可以了。也就是说,当增大耦合强度ε的时候,只要{εi}(i=2,3,…,K)都落入α平

im12

λ面的稳定区(ε=w2是最后一个),那么在这个耦合强度下,系统(8.27)就达到了同步态w

・

=…=w

k

。这时,再来看方程

δW

・

tv

=

Γ󰃁BδWDFs(W)+ε

tv

(8.30)

由于方程(8.27)中的w

im

已经相等,所以(8.30)的同步问题也就又和前面完全同步问题一样λtiελtv了,也即要求在增加耦合强度时{εK)都落入到α平面的稳定区(i}(i=1,2,…,N-l)是

tv

最后一个),这里{λ而如前所述,系统(8.30)的收敛意i}(i=1,2,…,N-K)为B的本征值。味着整个系统的部分同步,但是如果需要在系统(8.27)达到同步后再加大耦合强度才能使得系统(8.30)达到同步,那么这个时候整个系统达到的就不是一个部分同步态而是一个完全同步态了。也就是说,要使部分同步在一个网络中发生,那么还需要系统(8.27)达到同步的时候系统(8.30)已经达到了同步,即同步的阈值ε而前面讨论(8.27)与(8.30)的同步问im>εtv。

题时只是需要看是否本征值与耦合强度的积都落入了α平面的稳定区,同时考虑到这些本征值的符合,实际上就需要imtvλ>λ21

(8.31)

至此,我们就可以总结一下如何判断一个给定的网络是否可以到达某个部分同步态:

1.寻找网络的镜像对称性,如果没有则不能达到部分同步,如果有,则适当的编号使得其

邻接矩阵关于其反对角线对称;

2.视情况通过坐标变换将邻接矩阵分块对角化;

3.分别求出两个块的第二大和最大的本征值,然后和(8.31)式比较,如果(8.31)成立,

则部分同步可以发生在这个网络中。

tvλ4.增大系统的耦合强度,当ε1进入α平面的稳定区时的εPaS,即为部分同步的临界耦合强度。

值得强调的是我们的判据具有一般性,在整个讨论中,我们每一要求网络拓扑结构的具体形式,也不局限网络的对称性,也就是说这个判据不仅仅适用于镜像对称的网络。8.4　数值模拟实验结果

下面我们开展数值模拟来验证部分同步现象,有了部分同步发生的判据,部分同步斑图将会更好的得到研究。

A.矩阵分析

欲验证以上的判据,不妨首先以Ns=2的环形网络为例。先看两个具体的例子:图8.4与图8.5。它们都是“交叉网络”,但是数值实验却表明它们中一个可以观察到部分同步,另一个却不能。现试用上述判据说明其区别。对于前一个例子,网络所包含的镜像对称性的对称矩阵就为Tc=F10,满足TcCcTc-1=Cc。而Tc为的本征矢矩阵就为

384物理学进展　　27卷

Sc=

1I52F5

I5-F5(8.32)

即TcSc=Sc(I5󰃁(-I5)),其中󰂽表示两矩阵的直和。用它的本征矢矩阵将Cc分块对角化,即得Sc

-10100-310000-10101-3000

ScCcSc=Ac󰂽Bc=

-1

100

010

-30

-2

21

-3󰂽

000

000

-30

-3-1

-1

02

0

-1

1

-21

这里Ac为行零的矩阵,其最大本征值为0。而Ac的本征值为

im{λk}{0,-0.4592,-1.4875,-2.7719,-5.2814}

Bc的本征值为

tv{λi}={-1,-2,-3,-4,-4}于是我们有

tmtv

λ=-0.4592>λ21=-1

所以(8.31)式成立,因而部分同步是可以发生的。这与前面的结论是相同的。对于后一个例

-1

子,其对称矩阵亦为Tb=F10,满足TbCbTb=Cb,因而就与上例有相同的本征值矩阵。用它分块儿对角化Cb得

-31011-200-1-11-21001-31-11

-1

SbCbSb=Ab󰂽Bb=01-210󰂽01-300

11

00

10

-20

-1-1-1

-100

-20

-4010

　　同样,Ab是行零矩阵,其本征值为

im{λi}={0,-0.8299,-2,-2.6889,-4.4812}

而Bb的本征值为

tv{λi}={-0.6837,-1.4206,-2.8654,-4,-5.0303}这样

tmtv

λ=-0.8299<λ21=-0.6837

　　即(8.31)不成立,因而部分同步不能发生。这也是与前面的数值结果相同的。

上述两个网络都是N=10的交叉网络,但是却有着不同的动力学特性。它们的不同在于有不同的非局域连接,这说明非局域连接对网络动力学产生巨大影响。

B.有两条长程连接的环状网络的部分同步

下面讨论更一般的情况。为此,必须先找到所有Ns=2且有镜像对称性的环状网络。首先,两条非局域连接会有四个端点,而考虑到环状网络的旋转对称性,总可以认为编号为1的振子是其中的某个端点。其次,又由于镜像对称性的要求,不难得到第二个端点与第一个端点之间的距离应该等于第三个端点与第四个端点的距离。由此如果将第二、三个端点的振子编号为j、k(k≥j),那么就得到第四个端点的振子编号n4=k+j-1。最后,同样是由于网络

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

385

图8.6　(a)平行网络标号规则;(b)部分同步区(圆点和方点均在0以上)

的旋转对称性,可知只要讨论

2≤j≤N/2和j≤k≤N/2的情况即可。在确定了四个端点之后在这四个端点上也只有三种连接的方式是符合对称性的要求的,它们是:

1.平行网络(parallelnetworks):两条平行的非局域连接可以被表为(1,n4)及(j,k),如图8.6(a);2.交叉网络(crossnetworks):两条交叉的非局域连接可以被表为(1,k)及(j,n4),如图8.7(a);

3.八字网络(lambdanetwork):两条非局域连接可以被表为(1;j)及(k;n4),如图8.8(a)。

这样只要变化i,j就可以得到所有符合要求的网络了。

图8.7　交叉网络的部分同步区域　　　　　　　　　图8.8　八字网络的部分同步区域

计算中,我们仍然采用Lorenz振子为网络中的振子,其参数和连接方式均与图8.1中所

示的例子相同。考虑到第二节中的理论适用于各种振子以及第二节中所列举的具体例证,有理由认为其他的振子与连接方式将得到类似的结论。

数值实验表明,各种平行网络都是可以发生部分同步的,实际上,平行网络可以看作是Ns

=1情况的一种加强,因而所有平行网络都可以发生部分同步并不难理解。而交叉网络和八字网络则有的可以部分同步有的不可以。这里,对任意网络,记可以使得某个振子与其镜像对

386物理学进展　　27卷

称的像同步的最小耦合强度为εPaS,即部分同步的临界耦合强度(如果一个网络可以发生部分同步,那么各个同步集团会在同一个耦合强度形成,见图8.1(b));同时,另记可以使任意两个

ε振子同步的最小的耦合强度为ε这样就可定义一个物理量Δc,即完全同步的临界耦合强度。

ε>0,则说明部分同步的发生早于完全同步,即当耦合强度增加时部分同步=εc-εPaS,如果Δ

ε=0,则说明部分同步和完全同步是在同一耦合强度发生的,能够先被观察到;相反,如果Δ

即在此网络中将不能观察到部分同步。另一方面,为了验证(8.31)式,还可以定义另一个物理

imtv

λ=λλ>0,则说明部分同步可以在该网络中发生。量Δ-λ这样若Δ21。

下面我们以N=20;Ns=2的情况为例,将理论与数值结果相互对照一下。图8.6(b),

ε与Δλ的比较。ε总是图8.7(b)和图8.8(b)为在前述三种网络中Δ对于平行网络(图8.6),Δ

大于零的,这即是前面部分同步总是能在平行网络中发生的实验事实。而有趣的是,对于所有

λ>0,这说明了平行网络中的部分同步解总是稳定的,而这一结论与上面的实验结果的k,Δ

λ>0的区也是吻合的。图8.7和图8.8是交叉网络和八字网络的情况。我们同样能看到,在Δ

ε>0。域,相应的Δ但是,不是所有的交叉网络和八字网络都可以观察到部分同步态,因而在Δλ<0的区域,Δε=0。总的来说,对于交叉网络,两个非局域连接越远离对称轴,则这个网络

越容易同步;对八字网络,两个非局域连接越接近,则这个网络越容易同步。例如,在图8.7(b)中(即交叉网络k=9的情况),部分同步只发生在j<7的区域。而在图8.8(b)中(即八字网j=3的情况),部分同步只有k<7的区域才能观察得到。

图8.9　部分同步可以发生的区域。(a)平行网络;(b)交叉网络;(c)八字网络

进一步来说,因为不是每一个交叉网络和八字网络都可以发生部分同步,自然要问上面的

结论是否跟网络的振子数有关呢?不妨在(2j=N-2k=N)平面上讨论这个问题,实际上,注意到j、k取值范围(1>2k=N>2j=N>0),我们只需要在该平面第一象限的角分线以下讨论。在(2j=N-2k=N)平面上,每一个点就代表一组确定的(j,k)值,也就代表一种确

λ。定的网络结构,因而也就可以分别计算出其Δ图8.9给出了在这个平面上使部分同步发生

λ=0的线,不同的线代表不同振子数的情况,图中的阴影区域即代表可以的分界线,即使得Δ

同步的网络。图8.9(a)的角分线下都是阴影区,正说明平行网络都可以发生部分同步的。而

(c)的阴影区则分别表示交叉网络和八字网络可产生部分同步;图8.9中各个分界图8.9(b)、

线都重合在了一起,这说明部分同步的发生不依赖于振子数。

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

387

C.Ns=3:一个例子

我们也可以对于Ns=3的情况进行讨论。图8.10(a)是一种特殊情况(N=20),三条非局域边中的两条已经固定为(3,11)与(11,19),然后变化第三条边(11-m,11+m),这样网络就只有一条对称轴。它的上面所定义的两个量的对比见图8.10(b),可见也不是每个这样的网络都可以发生部分同步。有趣的是如果没有那两条固定的非局域连接,则网络与Ns=1的情况相同,故而对1≤m≥9部分同步都可以发生。但是,现在对于Ns=3的情况,虽然

网络中有了更多的非局域连接,可能同步的网络拓扑位形却变少了。

图8.10　Ns=3的一种部分同步情况,(a)为网络拓扑图;(b)为同步区域图

D.多对称性网络中的部分同步:斑图竞争和选择

让我们看一下在一个网络中有多个对称性的情况。我们曾经给出图8.3(c)所示网络的两个对称性的所对应的矩阵。而数值实验表明与这两个对称性相连的同步解中有一个是不稳

imtv

定的,另一个是不稳定的。试验证之,经过T1的本征矢矩阵S1的变换后,λ=-1>λ21=-3,即(8.31)成立,因而部分同步能依此对称性发生。而经过T2的本征矢矩阵S2的变换后,

imtv有λ=-3<λ21=-1即(8.31)不成立,因而部分同步不能依此对称性发生。

我们通过上述几个不同的方面讨论得到了网络的拓扑结构与其动力学性质的关系:通过给网络加入非局域连接,可以使原来不稳定的部分同步解(如在规则的环状网络中)实现稳定(可以被观察到);通过给原来的网络一个小的改变,也可以增加新的部分同步态;但不是系统中有更多的非局域连接部分同步就越容易被观察到。这些结果都说明了网络的拓扑结构(尤其是非局域连接)对网络动力学性质有巨大影响以及这个问题的复杂性。

8.5　李指数(LE)谱:部分同步的表现形式

内在分叉导致部分同步,暗示了耦合动力学系统的相变行为。现在我们观察部分同步发生时李指数谱的行为。图8.11给出了一个环状网络的最大四个李指数,该网络有6个振子,没有长程连接,其连接矩阵:

388物理学进展　　27卷

-21

-2

01

-2

001

-2

0001

-2

10001

-21

Cr=

0001

1000

102

10

1

　　从图8.11(a)可见:在吸纳共同的耦合强度下,第三和第四李指数同时穿零变负,这是环

(degeneration)。状网络李指数(LE)的“简并现象”

我们继续观察加了一个长程连接(1,4)时的李指数谱的情况。该网络的连接矩阵为-310101

1

Csr=

-2

1

-2

01-3001

-2

0001

-201011000

102

10

1

　　如图8.11(b),我们发现简并现象就消失了,而其第四大李指数在一个很小的耦合强度就

变负了,该耦合强度刚好对应于部分同步的发生。这可以理解为由于长程连接而导致对称性破缺。有趣的是在相同的耦合强度下d1,4,d2,3,d5,6同时变0,而d1,2,d4,3,d1,6,d4,5都不变0。因此,李指数简并破缺,部分同步发生。

图8.11　(a)环状网络(N=6)的最大的四个李指数;(b)Ns=1的环状网络(1,4)的最大的四

个李指数;(c)d1,4～ε;(d)部分同步流型(SM,方块线)及横向系统(TM,叉形线)的条件李指数。(c)(d)网络与(b)同

我们下面对以上的三个现象利用前面第三节关于网络邻接矩阵的讨论来解释。首先,对于一个环状的网络,其临界矩阵的本征值可以写为

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

389

πk-1σk=2-1+cos2

N

,(k=1,2,…,N)(8.33)

　　显然,存在关系:

σσσ2=σN,3=σN-1,…,

N+12

=σN2+2

(8.34)

　　环状网络的邻接矩阵的本征值一定是成对出现的。所以当增大耦合强度时,环状网络最

CσσC大的两个非零本征值会使得ε这时系统就有两个李指数同2和εN同时落入γ平面的稳定区。时变负,而考虑到矩阵的最大本征值为零,故而系统还有一个正的李指数以及一个零李指数,

所以前面的两个李指数正是第三及第四大李指数。这时如果加入一条非局域的连接,由于部分

tvimimλελλ同步的发生,ε=0与ε这也就解释了为什么部1已经进入了稳定区,还有12在稳定区外。

分同步当系统还有两个正李指数的时候发生。为了进一步验证上述结论,我们还可以计算方程(8.24)和(8.25)的条件李指数。如图8.11(d)所示,当横向系统的最大条件李指数(虚线)穿过零的时候正好就是系统发生部分同步的时候,而这时的不变同步流形的李指数(实线)正好还有两个是正的。因此部分同步的发生的确是系统在网络同步出现之前的一种合作分岔行为,这种分岔反映在(耦合网络)系统动力学和相空间结构的变化上。8.6　应用:稀疏连接的时空系统的同步上面得到的理论判据我们可以自然地应用到两个稀疏连接的时空系统中。我们考虑一个网络由两个群构成,N个混沌振子耦合成链形成独立的群,群内部耦合强度为ε,两个群之间加入了nl条稀疏连接L1,L2,…,Ln1,群之间的连接强度为r,如图8.12(a)所示。在两个群完全同步之前,当两个群对应的元素都互相达到同步时,我们称这个网络发生了群同步GrS。我们用分析部分同步的方法分析群同步。显然,整个系统具有镜像对称性,按照其对称性,我们可以把整个网络系统分成两部分(子网络):G1=X+X′代表全局同步流形上的动力学,Gr

lr

=X-X′代表群同步流形上的动力学。同样,我们可以得到群同步发生的判据:λ2>λ1。

图8.12　(a)两个存在稀疏连接的社区网络;(b)社区间强耦合极限下的网络合二为一,可以分解

为闭型小集团环状网络;(c)与开型小集团链式网络;(d)两种子结构

事实上,网络的拓扑性质对群同步的发生有很大影响。有的网络结构,无论我们如何加强群之间的联系,都不能达到群同步。这里我们不妨把群间连接强度极限化,即r→∞。这样群之间耦合的振子可以视为两群共同的振子。这样整个系统被分为三个部分,如图8.12(b),

390物理学进展　　27卷

而这三个部分属于两种结构形态:闭型小集团(大小为No,如图8.12(c))和开型小集团(Nc,如图8.12(d))。我们可以解析得到这两种结构形态的群同步流形上的本征值。对于开型小集团,群同步流形上的本征值为:

(2m-1)ππrr

λεε=-4sin2,m=1,2,…,N0,λsin21=-4

2(2N0+1)2(2N0+1)

　　对于闭型小集团,群同步流形上的本征值为:

ππmrr

λεε=-4sin2,m=1,2,…,Nc,λsin21=-4

2(Nc+1)2(2Nc+1)

对于整个系统,全局同步流形上的本征值为:

ππm1

λεε′=-4sin2,m=0,1,2,…,N-1,λsin22=-4

2N2N　　通过对本征值的分析,我们可以得到如下有用的结论:

1.耦合网络系统中存在的两种小集团的规模越大,那么越难以满足群同步发生的条件:1rλ2>λ1。

2.如果耦合网络系统中存在一个开型小集团,它的规模超过了整个系统的一半大小,即No>N/2,那么群同步发生的条件不能满足,无论如何不能达到群同步。

3.大小为n的开型小集团与大小为2n的闭型小集团群同步能力相同。4.如果要两个群达到同步,那么在它们

之间最少加入两条稀疏连接,并且最好的方式为使得耦合网络系统被分为两个相同大小的开型小集团,一个二倍大小的闭型小集

团。

进一步,我们分析网络的参数空间,举两个群为例,每个群由N=10的链构成,在两个群之间加入稀疏连接:L1=3,L2=8。如图8.13所示,我们可以看到ε-r参数空间被分成了5个部分:完全不同步态(US)、群同步态(GrS)、群内同步态(IS)、过渡区间、完全　图8.13　稀疏连接的社区网络的相图,US为非同步同步态(CS)。其中箭头显示了网络通向完全区,IS为链内同步区,GrS为链间同步区,同步的不同路径,揭示了时空网络中的一种相变行为。8.7　小结

CS为整体同步区

作为结论,本章讨论了由于网络拓扑结构的对称性破缺而导致的部分同步动力学。揭示了不同网络多种形式的部分同步斑图。部分同步斑图的存在依赖于网络拓扑的对称性质。另一方面,部分同步斑图的发生依赖于相应解的稳定性。我们给出了判断给定网络部分同步斑图稳定性的方法和判据,这是基于对部分同步斑图流形和横向流形的李指数谱的比较。这个判据已被大量的模拟证明。基于这个判据,我们也讨论了由网络的多种拓扑对称导致的部分同步斑图的竞争和选择问题。对于环形网络,我们系统地研究了当2条长程连接引入时的部分同步斑图的

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

391

动力学问题。发现:对于一些存在2条长程连接的拓扑对称环形网络,部分同步斑图的状态可能是不稳定的。本章给出了同步和异步情况下的相图。我们发现部分同步斑图依赖于网络拓扑,而与节点或连接的数量、节点上的动力学等等无关。我们对节点上的不同动力学进行了检验:洛仑兹振子、罗斯勒系统和非线性映象(如:logistic映象和圆映象)等。所有这些大量的数值实验支持我们的判据和结论。我们还讨论了部分同步斑图的李指数谱,发现长程连接的加入破坏了环形网络李指数谱的简并性质,这可以理解为网络拓扑结构的对称性破缺。

需要注意的是:上述部分同步斑图现象发生的判据具有一般性,因为我们没有限定任何具体的同步形式和耦合系统。只要网络的对称属性已知和相应的部分同步斑图解可被解出时,本章我们建议的判据都可用。而且,我们的讨论也可以延伸至由不同节点组成的网络和非线性耦合的情况。对于异质网络,我们同样可以块对角化它们的连接矩阵,比较两个矩阵块的本征值。

当前的研究揭示了在达到全局同步之前的耦合动力学系统固有的合作分叉行为。我们所揭示的现象和本章进行的理论探讨对于揭示时空系统的斑图动力学有着潜在的重要意义。它们也可以应用到对空间延伸系统和复杂网络的同步时空斑图的控制中。例如,我们可以通过边的加减或者改变耦合强度的分布来控制部分同步斑图。我们希望当前的研究能够对同步动力学、斑图动力学以及网络动力学的研究有所启发。

9　具有SW和SF拓扑的束流传输网络中的同步与控制

9.0　引　言

第728两章我们讨论了复杂网络系统中的动力学行为的完全同步,发现了在耦合混沌振子网络中的部分同步现象,并作了理论分析。作为这个令人关注的课题的继续发展和应用,本章进而讨论节点为非线性(混沌)方程的所谓“混沌复杂网络”中动力学的同步和控制。近年来,牵制控制等方法被应用于复杂网络的动力学控制[205～209],发现在全局耦合的条件下只通过对部分节点可实现复杂动力学网络的平衡点等控制[207～209],从数学上还证明了仅用单一控制器即可达到复杂网络中单一目标的控制[254]。然而,许多现实社会中的复杂网络经常由多个“社区”“群落”、或“模块”等构成,除了社会网络还有通讯网络及生物网络等都各有不同的层次和控制目标,如互联网就是一个由很多不同层次的多局域网络组成的,并不是一个单一的网络。因此,实际上需要把一个复杂网络划分为多个局域网络,而每个局域网络有各自的行为或功能,即存在不同的目标控制,人们期望实现网络分区的多目标控制,果能如此必将开辟复杂网络的多种用途。我们在第2章介绍过规则网络,其中有一种具有实际应用潜力的束流传输网络(BTN),也属于“混沌复杂网络”。不失一般性,本章将以BTN为典型实例,从理论上通过利用小世界(SW)效应和无标度(SF)特性的拓扑结构,构造SW/SF2BTN系统,进而提出驾驭束流传输网络中的束晕-混沌的新思路和新方法,达到所需的整体控制目标以及多目标的分区同步控制目的。这里讨论的BTN系统还与当前能源问题有一定的关系。由于煤和石油的储量有限,而且大量使用化石能源已经引起了严重的环境污染和导致全球气温上升,这引起了各国的极大忧虑。人类解决能源的途径颇多,但是从目前知识看,普遍的共认是,其根本途径是利用核能。如何更有效地发展洁净裂变核能仍是21世纪面临的十分重要课题。为了战胜目前常规核电的主要弊端(铀资源利用率低

392物理学进展　　27卷

(1%)、易致超临界事故和核废料后处理难等)问题,1993年西欧核子中心(CERN)诺贝尔奖获得

者C.Rubbia小组,提出强流加速器驱动的放射性洁净核能系统(ADS)[247,248],巧妙地把20世纪两大核装置———粒子加速器和核反应堆结合起来,构成了新的更安全、更干净、更便宜的洁净核能系统。但是,由于强流加速器中质子流产生了特有的束晕2混沌现象,成为解决洁净核能系统的关键问题之一:必须提出新方法和新技术来解决对束晕2混沌的有效控制问题。这就是洁净核能系统面临的挑战性课题之一[248～250]。

通常,强流加速器中的束流传输通道是由几十上百或数百上千个电磁聚焦单元所组成的规则的BTN系统。虽然已经提出一些控制束晕-混沌的方法[247,248],但是主要是非线性控制方法,一般代价高。如何利用当前复杂网络的研究进展来解决束流传输网络(BTN)中的束晕-混沌现象,是我们理论上的新思路。本章将结合复杂网络的新发现,提出具有SW和SF拓扑结构的束流输运网络,所谓SW2BTN和SF2BTN系统,首先从理论上寻找在SW2BTN和SF2BTN系统中束晕混沌的同步和控制方法,为实验和工程设计提供一定依据和参考。同时这更重要的是为其他“混沌复杂网络”系统的同步与控制提供新方法和新途径,因此本章讨论的思想和方法具有普适性。9.1　传统规则的BTN系统可控性的稳定性原理通常一个典型的传统的BTN系统是由许多螺线管或四极磁极子产生的周期性聚焦磁场组成的束流传输通道,它们可以构成直线型加速器或圆环形加速器,当一束具有园形的质子束在初始满足Kapchinsky2Vladimirki(K2V)分布时,质子束在BTN系统中水平方向运动的无量纲包络方程为[247,248]

d2rbk1-3=02+kZ(S)rb-rbdsrb

(9.1)

　　我们假设周期外场和质子束都是轴对称的,且在周期聚焦磁场通道网络中质子速度和周期长度保持不变。为了控制或抑制质子束的径向发散。令2gG为欲加的非线性反馈控制器

(函数),则受控BTN系统质子束包络方程可改写成:

d2rbk1-3=-gG2+kzrb-rbdsrb

(9.2)

在非线性常微分方程的特定条件下,我们提出和证明了非线性常微分方程(9.2)的稳定性定理如下[249]:

对于方程(9.2),以下仅列出代表性的3种非线性控制函数形式:

),G=-αsin(y-β

)2G=-α(y-β

22-)],f(x)=-21-xeG=α[f(y)-f(β

σσ

σ

x2

(9.3)

其中,α,β,σ为正的常数,只要参数α选择得合适,对于受控方程(9.2),在原点的适当小的邻

)时,则rb(t)→β(t→∞),这时β是所要控制的目标。域内,ε→0(t→∞

这个定理的重要性是:对于传统的BTN系统,必须应用非线性反馈控制方法才能实现对束晕2混沌的稳定控制,我们已经提出多种非线性反馈控制器,大量数值模拟证明了这个理论

的正确性。但是,根据目前复杂网络新发现———小世界(SW)效应和无标度(SF)特性,启发了

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

393

我们从理论上来重新研究束晕-混沌的控制问题,我们的新研究发现:在具有SW或SF拓扑结构下的束流输运网络中,大可不必应用“非线性反馈”,而仅仅采用“线性反馈控制方法”就可实现束晕-混沌的控制目标。我们已经提出束晕-混沌的控制和同步的线性耦合方法[250],实现了束流的不动点和周期稳定控制。具有小世界网络的BTN中的控制和同步效率高于随机耦合网络以及最近邻耦合网络[250～252],复杂网络的各个因子对耦合的网络同步也有一定的影响[92,93]。我们下面采用了小世界特性的两种拓扑结构:一是WS的小世界模型[1],二是我们提出的网络度不变(DS)的小世界模型[92,93],并比较了这两种情形的同步能力。我们还计算了复杂网络的拓扑特性因子:平均最短路径、平均集群系数对这两种网络模型同步能力的影响。特别是我们拓广了全局耦合和分区控制方法实现不同目标的网络分区控制[252,253],具有应用前景。

9.2　小世界拓扑结构的束流传输网络(SW2BTN)

考虑具有N个完全相同的、线性耦合振子的小世界网络,每个节点表示n维的连续时间动力学系统:

・xi=f(xi)+cn

j=1

6NaijΓxj,　　i=1,2,…,N(9.4)

其中:xi=(xil,xi2,…,xin)∈R是点i的状态变量;c为耦合强度;Γ为连接耦合变量的矩阵,即为耦合网络的矩阵A。Pecora等在这基础上提出了用主稳定函数方法确定动力学网络同步的稳定性[66,215]。我们利用第7章中网络同步判据可知,类型I网络的同步化能力可以用耦合矩阵的第二大特征值来刻画,较小的λ对于类型II网络2(负值小,绝对值大)导致较强的同步化能力。的同步化能力可以用耦合矩阵的特征值之比R=λN/λ2来刻画,随着R的减小,耦合网络的同步能力增加;随着|λ2|的增加,所需要的耦合强度c减小,耦合网络的同步能力增强。

两种小世界模型的耦合网络中的节点动力学是束晕-混沌的粒子包络方程(9.1)变换为:

dx1=x2

dt

dx2K12

(9.5)=-(a+bcos(x3))x1++3

dtx1x1

dx3=ωdt

　　我们分别计算了网络规模N=1000时,WS模型和SD模型的矩阵特征根,可以得到随着网络演化概率p增加,R值的变化情况,如图9.1所示;随着起始规则网络节点度k的增加,R值和|λ2|值的变化情况如图9.2所示,从两图中数据可以得出随着演化概率的增加,R值减小、|λ|λ2|值增加;随着起始规则网络节点度k的增加,R值减小、2|值增加;表明随着演化概率p增加,两种网络的同步能力都增强了;随着起始规则网络节点度k的增加,两种网络的同步能力也是增强了。图9.1和图9.2中的变化趋势都是在p值比较小时有一个比较急剧的变化,之后随之而来的是在演化概率p=0.4之后同步能力平缓变化,也就是说长程边的增加对网络同步能力的影响减弱了。同时,我们研究了平均最短路径距离和平均集群系数对两种网络模型同步能力的影响。图9.3示出R值随着网络的平均最短路径距离的变化情况,网络规模N=1000时,初始网络节点度k=4,随着平均最短路径距离的增加,R值增加、|λ2|值减

394物理学进展　　27卷

小。图9.4给出R值随着网络的平均集群系数变化,随着平均集群系数的增加,R值增加、|λ计算结果表明随着平均最短路径距离的增加,平均集群系数的增加,都导致两2|值增加。种网络的同步能力减弱。

图9.1　在不同模型下网络同步能力与演化　　图9.2　不同初始网络节点度下网络同步能力与概率的关系　　　　　　　　　　　　　　　　　演化概率的关系图9.3　网络同步能力与平均最短路径距离的关系　　图9.4　网络同步能力与平均集群系数关系

　　以上结果表明:尽管具有小世界特性的WS模型和SD(度不变)模型的生成机制不同,但是它们的同步能力随复杂网络特性因子的变化的趋势是基本相同的。9.3　SW2BTN中同步的线性耦合控制[250]

传统的BTN系统中束流失匹配、离子的非线性共振等作用都会导致强流离子束产生的束晕-混沌现象。我们利用WS模型[1]和SD(度不变)[92,93]二个小世界模型分别生成线性耦合网络,每个

节点为束晕-混沌振子,构成SW2BTN系统。从线性耦合网络的同步判据可知,如果取合适的耦合强度c,则可能使束流输运网络达到匹配同步。节点的束晕-混沌方程就是(9.5)形式,系统的参数

ω=2π、K=5;当方程的初始条件为(x1(0),x2(0),x3(0))=(1.0,分别为:a=1.65、b=1.25、

π)时,BTN系统处于束晕-混沌态。0.5,2数值模拟计算发现:只在每个节点束晕混沌方程的变量x2方向上引入线性耦合控制项G,才能较好地达到同步。变量x2的受控方程为:

dx2K12

(9.6)=-(a+bcos(x3))x1++3+G,

dtx1x1其中:控制器为G=c

j=1

6

N

aij(xi2-xj2),i=1,2,…N,根据前面对网络因子与同步稳定判据的

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

395

相关计算以及所拥有的计算设备的计算能力条件,不失一般性,选择耦合网络规模N=100,开始规则网络节点度k=4,并且保证演化过程中耦合网络保持连通,这样尽可能使束流传输网络达到同步,同时耦合网络中每个振子的初值各不相同:x1∈[0,2]、x2∈[0,1]、x3=2π,我们分别计算了不同耦合强度c=2、c=5和演化概率p=0.04、0.4、0.9时的WS模型和SD模型下SW2BTN束流输运网络的同步情况,表征网络的同步状态可以用最大同步误差:

Dmax(t)=maxk=1

6

n

(xik(t)-xjk(t))

2

12

,　i,j=1,2,…,N,(9.7)

　　当任意两个节点之间的最大同步误差Dmax趋于0,则SW2BTN系统达到同步态(流形),如图9.5所示。从图9.5可见对于WS模型,耦合强度c由2增加到5后,演化概率p=0.04时的两种同步误差曲线均趋于0了,演化概率p=0.4、0.9时两种同步误差曲线在不同耦合强度下都趋于0。对于SD模型其结果完全类似。研究表明:小的耦合强度不足以使所有演化概率下的耦合网络达到同步,实际的动力学网络2束流传输网络同步情况与前面我们仅计算耦合网络模型的同步能力的计算是相一致的,这个结果也

图9.5　SW模型c=2时最大同步误差有助于我们在实验或工程上找到比较经济的同步控　

Dmax随时间变化情况

制的参数。

9.4　SW2BTN中同步的噪声控制方法[251]

我们在束流输运网络中每个节点束晕2混沌动力学系统的x2变量方程上引入一个高斯型噪声驱动项,并且每个噪声驱动项不相同,加入噪声驱动项后,变量x2的受控方程变为:

dx2K12

=-(a+bcos(x3))x1++2+ξci(t)dtx1x　31

j=1

6

N

aijH(xj2),　i=1,2,…,N(9.7)

其中:ξi(t)是分布在[0,1]内的高斯噪声驱动项。

图9.6　(a)SW2BTN系统中同步误差与耦合强度c关系;(b)SW2BTN系统中同步误差与演化概率p关

系曲线图

396物理学进展　　27卷

图9.6分别示出(a)SW2BTN系统中同步误差与耦合强度c关系和(b)SW2BTN系统中同步误差与演化概率p关系曲线图。图9.6(a)结果表明:WS2BTN和SD2BTN系统在噪声驱动控制下分别存在一个耦合强度临界值Cc=0.6和7,5,只有当C≥Cc时,最大同步误差Dmax曲线才能趋于零,即实现了网络同步。图9.6(b)也表明WS2BTN和SD2BTN网络同步误差与演化概率关系,同样分别存在一个临界值pc,当p≥pc=0.4和p≥pc=0.6时,这两种拓扑的束流输运网络同步误差趋于0,达到了同步控制。

比较表明,WS2BTN和SD2BTN系统受噪声驱动达到同步所要求的二个主要参数临界值并不相同,显然,需要选择不同的网络演化参数才能实现不同拓扑下的SW/SF2BTN系统中噪声驱动控制同步,这样既能够很快达到同步控制的目的,又能够提高同步效率,节约控制成本。9.5　SF2BTN中同步的噪声驱动的同步控制效果[252]

近年来,复杂网络研究发现现实中的许多复杂网络(如:WWW,新陈代谢网络和Internet等等)节点度分布具有幂律分布的形式,我们根据无标度(SF)BA模型,构造了SF2BTN系统。

图9.7　BA模型同步能力与网络规模N变化曲线图　　图9.8　BA模型同步能力随m变化曲线图

图9.7给出SF2BTN模型同步能力随网络规模N改变而变化的情况,随着规模N的增大,R值增大,并且|λ2|值减小,可见BA模型网络的同步能力随着网络规模的增大而减弱;网络同步能力随每次引入的新边m的变化如图9.8所示,随着m增大,R值减小而|λ2|值增大,BA模型网络的同步能力随每次引入的新边条数m的增大而增强,网络SF2BTN变得容易同步。计算结果表明:对于无标度模型,网络节点数的增加会引起整个网络同步能力的减弱,但新引入边数的增加却不会引起整个网络同步能力的下降.计算发现当m=3时,不同耦合强度下噪声驱动束流输运网络都能够达到同步,但当m=4或5的时候,增加耦合强度c并不能使噪声驱动束流输运网络达到同步.因此在对BA拓扑束流输运网络实现噪声驱动同步时可选择比较小的耦合强度,这样同样有利于提高控制效率。9.6　SW2BTN和SF2BTN中的周期稳定控制[251,252]

对于束流输运网络,关心的是把束晕-混沌控制到所需的稳定的单周期态。为此,通过利用一个简单的线性控制器G,可以实现了对束流输运网络的周期轨道的稳定控制。类似地,SW2BTN和SF2BTN中受控的节点动力学方程为:

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

397

dx2K12

=-(a+bcos(x3))x1++3+caij(xj2)+G,　　i=1,2,…,Ndtx1x1j=1

6

N

(9.8)

我们设计了一个简单线线性控制器是G=s(x2+d),根据前面的计算,选取耦合强度c=1,

控制器的参数s=-5和d=1,WS拓扑演化概率p=0.4,SD拓扑演化概率p=0.6,m=3。分别在三种拓扑结构:WS拓扑、SD拓扑和BA拓扑下构造了束流输运网络SW2BTN和SF2BTN,利用上述简单线线性控制器都可以实现稳定控制到单周期态,图9.9和图9.10示出了

x2随时间变化曲线图及最大同步误差Dmax随时间变化曲线图.从图可见:不同拓扑时网络的

最大同步误差Dmax都趋于0,这表明整个动态网络在简单线性控制器控制下达到了周期同步.每个节点上的振子x1、x2变量的轨迹都是单周期态,x1的周期振幅为0.5,不超过束流匹配半径0.7。同时控制后的李雅普诺夫指数为负值(-2.494,-2.5061)证明确实实现了单周期态的稳定控制。

图9.9　SW2BTN系统中x2随时间变化曲线及最大同步误差Dmax随时间变化曲线图。N=100,K

=6,c=1

图9.10　SF2BTN系统中x2随时间变化曲线及最大同步误差Dmax随时间变化曲线图。N=100,m=3,c=1

9.7　SF2BTN中的全局耦合分区牵制控制方法[252～254]

类似于方程(7.1)(7.2)(7.3),我们可以讨论同步态分别为平衡点、周期态和混沌态情形。将连续时间耦合动态网络分为两个局域网络G1和G2,局域网络规模分别为N1和N2,并且有N=N1+N2。分别在局域网络G1和G2选择I1、I2个节点在全局耦合作用下引入牵制控制项,则式(7.1)变成:

398

・

N

物理学进展　　27卷

xi=f(xi)+cxi=f(xi)+cxi=f(xi)+cxi=f(xi)+c

j=1

6666

NNN

aijH(xj)-cd(xi-s1),i=1,…,I1aijH(xj),i=I1+1,…,N1

aijH(xj)-cd(xi-s2),i=N1+1,…,N1+I2aijH(xj),i=N1+I2+1,…N2

(9.8)

・

j=1

・

j=1

・

j=1

假设δxi(t)=xi-s,s=s1,i∈N1,s=s2,i∈N2则式(9.8)可简化为:

+c

δdxi=f(xi)-f(s)dt

j=1

6

N

～

aijδxj,其中耦合网络A的对角线元素为:aii=[a11-d,…,al1l1-d,aI1+1I1+1,…aN1N1,

～

～～

aN1+1N1+1-d,…,aN1+I2N1+I2-d,…aNN]。耦合网络矩阵A的其它元素与耦合网络矩阵A相同。

根据文献[254]中证明可得耦合网络矩阵A的特征值λ<0,即耦合网络矩阵A的所有特

δdx(t)征值为负值。假设δx(t)=[δx1(t),…,δxN(t)],对s(t)求导可得=Df(s)(t)δx(t)+

dt

c

j=1

TTT

λ1,…,N],0>λ1

6δx(t)A,设A=WJW为耦合矩阵A的特征值分解,其中J=diag[λN

～～

～～～

δdyk(t)λδλkI)δ≥…≥y(t)=δx(t)W,则可得:=Df(s(t)+cyk(t)。N,

dt

为了确保控制策略有效,文献[254]中提出了4个定理,其中主要2个定理是:

1(T

定理(1):设μDf(s(t)))+Df(s(t))的特征值,且μ(t)=i(t),(i=1,…,N)是矩阵

2λmaxi=1,…,Nμi(t)。如果t>0时有u(t)<-c1,则整个耦合系统能够局部地(指数地波)牵制控

制同步到s(t)(平衡点)。

～

λ定理(2):假设0>λ如果存在正定的对角线矩阵P=1≥…≥N是矩阵A的特征值。

T

diag[p1,…,pm],Δ=diag[ΔP(f(x,t)-f(y,t)+Δy-Δx)1,…,Δm]且有η>0,使得(x-y)η(x-y)T(x-y),且对k=1,…,m,有Δ≤c1<0,则整个耦合系统能够全局耦合指数地波k+λ牵制控制到同步s(t)(流形)。

利用与文献[252]类似的证明,我们不难把牵制控制方法推广到多分区的不同目标的同步控制。由定理(1)和定理(2)可知:只要全局耦合强度c足够大,而局域反馈控制选择满足稳定性要求,则全局耦合下分区牵制控制可以使整个混沌复杂网络系统分别控制到两个同步态上。9.8　基于SW2BTN和SF2BTN系统的牵制控制[251～253]

本节试图把上述全局耦合分区的牵制控制方法分别应用于SW2BTN系统和SF2BTN系统,实现对它们的不动点和周期态两个目标的分区稳定控制。9.8.1　SW2BTN系统中的束晕2混沌的控制[252,253]

以WS模型为基础,在网络中每个节点引入束晕-混沌振子,我们构造出基于WS模型的

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

399

SW2BTN系统。然后应用上述全局耦合分区的牵制控制方法实现多目标分区控制目的。我们对于方程(9.5)通过数值模拟计算得到束晕2混沌振子有两个平衡点分别为:xeq1=[1.41,0,0],xeq2=[-1.41,0,0]。对于整个网络G,考虑我们计算机计算能力有限,不失一般性,可

取其节点总数N=50,初始节点度k=6,并将整个网络分为两个局域网络G1和G2,如下式所示:

G1={1,2,…,N1}G2={N1+1,N1+2,…,N}

(9.9)

其中,局域网络G1节点数N1为20,局域网络G2节点数N2=N-N1为30。

牵制控制在选择牵制控制节点时有两种方法:(1)随机牵制,在网络所有节点中随机地选取若干个节点引入牵制控制;(2)择优牵制,在网络所有节点中按照每个节点的度的大小选取若干个度大的节点引入牵制控制。两者研究比较发现:采用择优牵制的同步能力更强,因此我们选用择优牵制的方法对所构造的束流输运网络进行同步控制[9]。9.8.2　不同平衡点的分区控制[252]

为此,在全局线性耦合的局域网络G1和G2中分别择优选择I1和I2个节点,并在这些节点上引入控制器u1和u2,可以得到:

・

xi=f(x)+cxi=f(xi)+c

j=1

6

N

a1jh(xi,xj)+u1　　i=1,…,I1　　　　　　　　　aijh(xi,xj)　　i=I1+1,…,N1

・

j=1

66

N

N

・

xi=f(xN1+1)+cxi=f(xi)+c

j=1

6

N

aN1+1,jh(xi,xj)=u2　　i=N1+1,…,N1+I2

(9.10)

・

j=1

aN1+1,jh(xi,xj)　　i=N1+I2+1,…,N

　　为了将这两个局域网络分别控制到两个平衡点xeq1和xeq2上,即有

x1=…=xN1=xeq1xN1+1=…=xN=xeq2

(9.11)

　　全局耦合函数h(xi,xj)设计为:

h(xi,xj)=(xj-xeq2)-(xi-xeq1)

(9.12)

　　根据数值模拟效果,控制器u1和u2分别设计为:

u1=-cd(xi-xeq1)　　(i=1,2,…,I1),

u2=-cd(xi-xeq2)　　(i=N1+1,N1+2,…,N1+I2)

(9.13)

　　通过数值模拟研究,我们得到,在局域网络G1中择优选择节点数I1=14,在局域网络G2择优选择节点数I2=24时,束流输运网络的全局耦合强度c=200,控制器的控制强度d=7.5,局域网络G1分别被控制到平衡点xeq1=[1.41,0,0],局域网络G2被控制到平衡点xeq2

=[-1.41,0,0],结果如图9.11所示。我们分别计算了控制后的Lyapunov指数,局域网

络G1的Lyapunov指数λ1=-50.067,局域网络G2的Lyapunov指数λ2=-49.985,这表明经过控制后两个局域网络分别达到了不同的稳定的平衡点。

400物理学进展　　27卷

　　　　图9.11　SW2BTN系统中分区的不同平衡点控制结果。(a)局域网络G1同步控制到衡点xeq1,c=200,d=7.5;(b)局域网络G2同步控制到平衡点xeq2,c=200,d=7.5

9.8.2　平衡点和周期态的分区控制[252,253]

由于混沌吸引子包含无数个不稳定周期态和平衡点,为了不同分区控制目标,我们通过引入控制器的方法,可以将束晕-混沌系统分别分区控制到所需的平衡点和周期态上。为此,我们采用上述提出的类似的全局耦合分区的牵制控制方法,将网络的一个局域网络G1控制到平衡点附近,另一个局域网络G2控制到周期态s上,即

x1=…=xN1=xeqxN1+1=…=xN=s

(9.14)

其中一个目标为平衡点xeq=[1.41　0　0],另一目标状态s为周期态。我们设计采用了适当的误差全局耦合h(xi,xj),即:

h(xi,xj)=(xj-s)-(xi-xeq)

(9.15)

并分别设计两区的控制器u1和u2为:

u1=-cd1(xi-xeq)　　(i=1,2,…,I1)

u2=-d2(xi2+1)　　(i=N1+1,N1+2,…,N1+I2)

(9.16)

实现控制的数值计算结果如图9.12所示,我们得到:在局域网络G1中择优选择节点数I1=14,在局域网络G2择优选择节点数I2=24时,束流输运网络的全局耦合强度c=200,控制器

的控制强度d1=7.5,d2=5,局域网络G1被控制到平衡点xeq1=[1.41,0,0],局域网络G2被控制到单周期态s。此时算得局域网络G1的Lyapunov指数λ1=-50.025,局域网络G2的Lyapunov指数λ2=-2.61,这表明局域网络G1和局域网络G2分别被控制到了稳定的平衡点

和单周期轨道上。

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

401

　　　　图9.12　SW2BTN系统中分区的平衡点和周期态控制结果。

(a)局域网络G1同步控制到平衡点xeq1,c=200,d1=7.5;(b)局域网络G2同步控制到单周期态S,d2=5

9.9　SF2BTN系统中的不同平衡点的分区控制[251～253]9.9.1　不同平衡点的分区控制我们以束晕2混沌振子为网络节点,利用BA模型构造出基于无标度拓扑的束流输运网络SF2BTN。假定,生成一个无标度网络G,网络G的节点总数N=50,初始节点和每次加入新节点所引入的边m=3。同样,我们将整个网络分为两个局域网络G1和G2,局域网络G1节点数N1为20,局域网络G2节点数N2N1为30。在局域网络G1和G2中分别按照节点度的大小择优选择I1和I2个节点,并在这些节点上引入与(9.16)相同的控制器u1和u2。数值计算得到,在局域网络G1中仅择优选择节点数I1=4,在局域网络G2仅择优选择节点数I2=4时,束流输运网络的全局耦合强度取c=100,控制器的控制强度d=20,局域网络G1被控制到一个平衡点xeq1=[1.41,0,0],局域网络G2被控制到另一个平衡点xeq2=[-1.41,0,0]。控制结果如图9.13所示。我们同样分别计算了控制后的两个局域网络的Lyapunov指数,此时局域网络G1的Lyapunov指数λ1=-61.012,局域网络G2的Lyapunov指数λ2=-62.343,这表明经过控制后两个局域网络分别达到了稳定的平衡点。

　　　　　图9.13　SF2BTN系统中分区的不同平衡点控制结果。

(a)局域网络G1同步控制到平衡点xeq1,c=100,d=20;(b)局域网络G2同步控制到平衡点xeq2,c=100,d=20

402物理学进展　　27卷

9.9.2　平衡点和周期态的分区控制

与上面类似地,以BA模型为基础,在网络中每个节点引入束晕2混沌振子,我们构造出基于BA模型的SF2BTN系统,然后应用上述类似的全局耦合分区的牵制控制方法来达到多目标分区控制目的。首先,我们将SF2BTN系统的一个局域网络G1节点控制到平衡点,另一个局域网络G2节点控制到单周期态上。数值计算得到,在局域网络G1中择优选择节点数I1=4,在局域网络G2择优选择节点数I2=4时,束流输运网络的全局耦合强度c=150,控制器的控制强度d1=10,d2=15,控制结果如图9.14所示,在局域网络G1中实现了平衡点xeq1=[1.41,0,0]的稳定控制,在局域网络G2中达到了单周期态s的稳定控制。此时局域网络G1的Lyapunov指数λ1=-47.873,局域网络G2的Lyapunov指数λ2=-7.41,这表明局域网络G1和局域网络G2分别被控制到了稳定的平衡点和周期轨上。　　　　图9.14　BA2BTN系统中分区的平衡点和周期态控制结果。

(a)局域网络G1同步控制到平衡点xeq1,c=150,d1=10;(b)局域网络G2同步控制到单周期态s,d2=15

9.10　结论和简单讨论

束晕-混沌的同步与控制是加速器驱动的放射性洁净核能系统关心的一个问题,从理论上我们构造了具有小世界拓扑和无标度拓扑的束流输运网络,不失一般性,可拓广于任意混沌复杂网络。研究发现:WS小世界模型和SD小世界模型的同步能力在同样的网络参数下基本相同;进一步研究了基于这两种小世界模型的束流输运网络的同步,采用全局线性耦合控制方法实现了束晕2混沌的同步控制。同时研究了应用噪声驱动控制方法,实现了对具有两种小世界拓扑和无标度BA拓扑的束流输运网络的束晕2混沌的周期态控制。特别是,应用全局耦合和分区的牵制控制方法能够实现在具有小世界和无标度拓扑特性的BTN系统(即SW2BTN系统和SF2BTN系统)中的多目标(平衡点和周期态)分区控制。比较SW2BTN系统和SF2BTN系统发现:为了实现复杂网络分区的多目标控制,利用小世界拓扑结构需要选取的节点

数明显多于无标度BA结构所选取的节点。这样,采用无标度束流输运网络(SF2BTN)系统相对成本较低。研究结果还表明:BA模型的择优机制能够比较快地通过部分关键节点把同步控制作用传播到网络中的每一个节点,容易实现控制。但当关键节点遇到攻击时,BA模型由于无标度特性将不如小世界模型的鲁棒性好。现实社会中的复杂网络节点众多,在同步控制

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

403

时如果采用每个节点都引入控制器的方法来实现对网络同步的控制,显而易见成本较高,实现起来难度较大。如果采用随机的或有选择的选取几个甚至一个节点加以控制,实现起来相对容易,能够提高控制效率。这里强调指出:本章研究结果具有理论意义和应用参考价值,不仅可以直接推广应用于其他复杂混沌网络的多目标同步控制,而且我们设计了节点为束晕2混沌动力学方程的网络通信线路,以便应用束晕-混沌网络来实现不同目标的分区网络的保密通信[255,256]。需要指出:如何使受控的节点最少就能达到多目标控制仍然是值得深入研究的一个问题。

10　复杂多智能体网络的协调控制

生物体(鸟群,鱼群等)的群集现象,近年来引起了不同领域科研人员越来越大的兴趣,包括生物、物理、社会和控制等领域的人们的关注。最初,科学家们发现生物界存在普遍的群集现象,即智能体之间局部相互影响的作用会使得整体涌现出一致的行为,比如说鸟群朝着同一个方向飞翔,鱼群绕着某个中心转圈等,人们试图建立理论模型来解释和揭示这个现象及其奥秘,近年提出了两个经典模型—Vicsek模型和Boid模型,并被广泛地研究。随着工业科技的发展,在某些危险环境下,多智能体系统继而取代原来传统的单个智能体,由于其灵活性和抗干扰性都有所提高,多智能体系统被越来越广泛地应用,人们期望通过模拟群集生物体之间相互协调合作的现象来控制实际群集智能体系统,以达到期望的目标,比如:群机器人,无人驾驶飞行器,自治海底探矿船和高速公路系统等的控制.本章结合国际进展情况,综述了我们近期基于Vicsek模型和Boid模型,从动态网络的拓扑结构和动力学方程以及控制策略方面研究该课题所取得的主要进展。10.0　引言

近年发现生物界普遍存在生物体(鸟群,鱼群等)的群集现象,即智能体之间的局部相互作用会使得整体涌现出一致的行为,人们试图对这一有趣的现象建立动态网络的理论模型来解释和揭示其奥秘。同时,随着科技和工业的日益发展,开展了对于群集智能体(如自治水下交通船,无人驾驶飞行器,群机器人等)的广泛应用,期望通过掌握生物体的涌现规律,进而应用到群集智能体中。群集智能体的研究已经成为一个热点.生物群集的涌现行为是自然界中难解的美丽现象[257～259]。我们经常可以看到生物体群集的现象,如鸟群、鱼群等,这些由大量个体(或者称之为粒子)组成的群生物体常常就如同单个的生物体一样运动,而且很多情况下,群集中并没有明显的领导者或者外界环境的刺激,粒子(个体)通常也不知道整体的群集信息。它们是如何形成群集并一致运动的呢?其运动具有什么特征呢?近年来,这方面大量的工作主要基于三种方法:拉格朗日方法[260～265],欧拉方法[266～269]和离散方法[193,267～277]。

1987年Reynolds提出了三条启发式规则:聚合、分离和速度匹配,并用计算机来模拟鸟群的飞行[269]。接着Vicsek等人提出了一个简化的模型,使用方向一致的最近邻规则研究涌现行为[270],Vicsek模型是最常研究的模型之一,如该群集在噪声情况下的行为[270]、间歇性行为和群聚(簇)的形成和系统初始阶段的运动特征[270]。还研究了网络拓扑连接性的稳定性分析[276～277]。具有方向一致这一重要特征的离散模型还有:Couzin模型[278],粒子之间除了具有

404物理学进展　　27卷

方向排列之外,还具吸引和排斥力,其模型可以产生更复杂的宏观斑图。Couzin等人还研究了具有领导者的模型,研究表明,领导者获取的信息可以在整个群集中传播。10.1　Vicsek模型的相关介绍及研究概况10.1.1　Vicsek模型的介绍

考虑N个粒子,均匀随机的分布在边长为L的矩形上,该矩形具有周期边界条件,假设所有的粒子都以相同的绝对速度v0运动,初始的运动方向是随机分布在[0,2π)上。在系统演化的每一步中,每一个粒子的方向都采用其半径为r的邻域内所有粒子的方向平均值θ〈i(k)〉R,R称为作用半径。

θθθ,i=1,2,…,N(10.1)i(k+1)=〈i(k)〉R+Δ

θ表示噪声。在不同的参数下,Vicsek模型具有不同的运动特征。这个方向排列的局这里Δ

部邻居规则和常速率运动的粒子模型对群集的研究具有很大的影响[270～286]。描述系统从零传输到涌现出相变广泛采用下面序参量[193,271,278,380～282]:

Φv(k)=1Nv0

i=16N

vi(k),0≤Φv(k)≤1(10.2)

Φ其中vi(k)是粒子i的速度,其方向为θi(k),v0=|vi(k)|,i=1,2,…,N是一个常量。v(k)具

Φv=0表有明确的物理意义,即它可以度量系统的平均动量。当Φv(k)µ0,系统具有涌现行为。明方向不存在极性,而Φv(k)=1意味着系统的方向收敛,即所有的粒子具有相同的运动方向。10.1.2　群集的自适应速度策略10.1.2.1　自适应速度模型

现在假设一般情况下不同的粒子具有不同的速率。v0是所有粒子的最大可几速率的平均值,即v0=

1N

6

Ni=1

vi0,vi0为粒子i的最大可几速率。一般来说,即使系统中粒子的运动呈现出

各向同性,这种情况下没有涌现行为,但是仍可能存在Φv(k)>0的情况。而且Φv(k)=1并不必然表示所有粒子的一致运动(coherence),除非要求所有的vi0都相同,即v0=vi0,i=1,2,

…,N。当粒子在稳态的情况下具有不同的速率时,则群集的一致运动是不可能达到的。所以Φv(k)不适用于描述系统的一致运动。

另一个常用的序参量是[274,283]

Φθ(k)=

1N

i=1

6e

N

θii(k)

,0≤Φθ(k)≤1

(10.3)

该序参量特别适用于表征网络耦合振子的同步性,它除了对速率的影响,只考虑方向的宏观特

征。当然对常速率Vicsek模型来说,显然有Φv=Φθ。

在本节中,我们给出自适应速度模型。首先定义粒子i在第(k+1)步的复数局部序参量为:

1θθii(k+1)i(k)

φ(10.4)=ej,i=1,2,…,N;k=0,1,2,…i(k+1)e6ni(k+1)i∈Γi(k+1)

其中eij

是单位方向矢量,Γi(k+1)表示粒子i在第(k+1)步的邻居集合,该集合含有ni(k+

1)个元素(包括该粒子自身)。复数序参量的幅值(局部极性)φi(k+1)表征粒子i在第(k+1)

θ(k)

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

405

φθ(k)的局部化形式,显然有步的邻居运动方向的一致性程度。i(k+1)可以看作是全局序参量φφ0≤幅i(k+1)≤1.φi(k+1)的值大表明粒子i的邻居具有比较高的方向一致性(如图10.1)。角θi(k+1)表示粒子i在第(k+1)步的运动方向,它是在Γi(k+1)集合中的粒子运动方向的均值。使用复数表示计算方向均值可以避免一些不期望的方向问题[275,276]。

图10.1　粒子i的局部极性φ图中箭头表示粒子i的邻居粒子的运动方向。为了简单起见,这i的图示。

些具有单位模长的矢量都画在同一个单位圆上。(a)粒子四下散开,φi≈0。(b)粒子具有一个相对强的极性,φ0。极性φi≠i=1当且仅当在集合Γi(k)中的所有粒子的运动方向一致

令Pi(k)表示粒子i在复平面上的位置,Pi(k)∈C。在自适应速度模型中,每一个粒子不仅调整它

的运动方向,而且还根据邻居方向极性的大小调整其运动速率。具体的说,粒子i在第k步的速率系数

α

φφα(c(i(k))是其局部极性的幂函数,即ci(k))Χφi(k),其中幂指数α≥0(图10.2)自适应速度模型数学形式表示如下:

αθii(k)φΔPi(k+1)=Pi(k)+v0××ti(k)e

θi(k+1)

φ(k+1)ei=

1ni(k+1)

Γ()j∈ik+1

6

e

θij(k)

(10.5)

i=1,2,…,N;k=0,1,2,…

αθii(k)

φ其中Δt=1。vi(k)≡v0×表示粒子i在第i(k)e

k步的速度,速度方向为θ因为对任何α≥0都有i(k)。α

φ0≤1,所以粒子的速率vi(k)i(k)≤

φ=v0×i(k)

α

满足约束0≤vi≤v0。

φ)局部极性φ　图10.2　速率系数c之间的函α(

幂指数α≥0反映了粒子根据其邻居的运动方向数关系。对任何α值,如果φ=1,则

φ)=1。φ)≡1;当0过局部极性如何,每一个粒子都始终以最大速率v0运动。而当α>0的情况下,只有当邻居的运动方向达到了一致,粒子才会以最大的速率运动。大的α值意味着在给定邻居运动方向的情况下,该粒子的运动速率较小。当α=1,粒子i的局

)为:部序参量就是粒子i的所有邻居的方向矢量的直接和。当α=∞,则c∞(φ

0,0≤φ<1+∞

(φ)φ(10.6)c∞==

φ=11,

这意味着粒子的局部极性达到1,即该粒子所有的邻居粒子的方向都达到一致了,该粒子才运动。

2维的离散自适应速度模型(10.5)可以推广到M维欧几里得空间。令

α<∞及0<φ<1,0406物理学进展

T

　　27卷

Pi=[pi1,pi2,…,piM]

表示粒子i的空间位置,i=1,2,…,N。粒子i的运动方向由单位向量di=[di1,di2,…,diM]T表示,‖di‖=1,-1≤dij≤1,j=1,2,…,M。粒子i和粒子j是邻居,如果满足‖pi(k)-pj(k)‖≤R。定义序参量:

1(10.7)ri(k+1)=‖6dj(j)‖ni(k+1)j∈Γi(k+1)显然0≤ri(k+1)≤1。则M维空间中群集的自适应速度模型可以表示为:α

ΔPi(k+1)=pi(k)+v0×ri(k)×di(k)×t,　　k=0,1,2,…

di(k+1)=

Γi(k+1)j∈

(10.8)(10.9)

6

dj(k)/‖

Γi(k+1)j∈

6

dj(k)‖,　　k=0,1,2,…10.1.2.2　仿真与讨论

假设粒子可以在整个复平面上运动,这样可以保证系统可以充分地在空间上演化,这比周期边界条件[271]更符合自然界的群集运动的特征。为了说明自适应速度的效果,考虑N个粒

子在复平面上运动,这N个粒子的初始位置随机分布在边长为L的矩形内。令粒子i初始分布的位置和方向分别表示为θ注意到初始分布的方向θi和Pi(0),i=1,2,…,N。i不是粒子的初始运动的方向θ下面计算粒子i的初始运动的方向θi(0)。i(0)和初始的局部极性φi(0):

1θθii(0)ij

φ=i(0)eΓ(0)ej∈6ini(0)

φ这意味着系统在一开始就采用了自适应速度策略。粒子i的初始的速率为v0×i(0)。

在仿真中,取N=300,L=5,R=2,计算都取200次结果的平均值。首先研究自适应速度模型中幂指数α对系统收敛概率p的影响。系统收敛概率p定义为:N个随机分布的粒子最终能够达到方向一致,并以相同的最大速率v0运动的概率。

图10.3　(a)收敛概率p和最大速率参数v0的函数关系[283];(b)收敛概率p和幂指数α的函数关系

图10.3(a)表明对任一给定的α值,收敛概率p是最大速率v0的递减函数,但是其降低速度随着α值的增大而减小;图10.3(b)表明对一给定的速率参数v0,收敛概率p是指数α的增函数,而且小的v0值能使系统有高的收敛概率。当常速率参数v0足够大的情况下,即使对常速率Vicsek模型(α=0)来说是不大可能收敛的,对自适应速度模型来说,充分大的α值仍然可以使系统有较高的收敛概率。特别的,当α=∞,从图10.3中可以看到,在当前参数下,收敛

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

407

概率趋近于1。该模型表明了即使在没有领导者和全局信息的情况下,粒子使用局部信息来调整速率和方向,可以达到较高的收敛概率。

注意到,不管群集是否收敛,所有粒子的速率总是可以达到最大值v0;但是方向到达一致性需要一定的条件。一般来说,在自适应速度模型中,粒子的速率动态变化的过程是比较短暂的,所有粒子的平均速率vave(k)一般来说是单调增大的,直到稳态平均速率vave(k)到达最大值v0。

令τ表示群集中所有粒子的速率平均值达到最大速率值v0的98%所需要的步数,τ服从:

τ≈β+γ(10.10)　　　　　　　　　　　　　　　log2α,α≥1图10.4(a)和(b)表明:τ与最大速度v0及密度N/L2没有明显的相关性。这里,γ≈1.4。即使

对很大的指数α=1000,大多数粒子趋近最大速率的平均步数也小于18步。图10.4(c)和(d)表明系数β和γ随感知半径R的增加而减小。当R趋近于6,粒子近似于全耦合,β趋近于1,即1步就可以达到稳态。

图10.4　(a)对不同的v0,暂态步数τ和α的函数关系;(b)对不同的N,τ和α的函数关系[283]。v0=

1.0;(c)对不同的R,τ和α的函数关系。v0=1.0.;(d)τ和α、R的函数关系。v0=1.0

为什么在自适应速度模型中,随着指数α的增大,收敛概率会增大呢?这是因为自适应速

度策略有阻滞局部粒子扩散的作用。在一个粒子的邻居方向集合呈现各向同性的情况下,φ≈0意味着邻居粒子的运动方向是四下散开的,根据自适应规则,该粒子的速率相对是很小的。即使对0νφ<1的情形,当α很大的情况下,粒子的速率仍然很小。所以这些粒子的相对位移就不会很大,邻居关系在下一步甚至以后都将继续成立,它们之间的通信仍然保持,这种通信是使粒子方向一致的作用因素。

从复杂网络的角度看[19,20],群集的拓扑可以表达成一个无向图G=(V,E):每一个粒子i表示图中的一个顶点vi;顶点i和顶点j之间无向边意味着它们是邻居,反之亦然。一个顶点属于的一个片是指能够从该顶点沿着图中的边所到达的那些顶点的集合[20].随着群集的演化,群集的拓扑图G(k)=(V,E(k))也在变化。这里关心的是群集在稳态情形下的最大连通片MCSG(maximalcomponentofswarmgraph)。最近工作表明,对一个在2维平面上运动的

408物理学进展　　27卷

群集,其收敛与否是因为群集粒子之间的连接性[275～277],而不是群集之间的长程作用[271,284]。

令S表示在稳态情形下中的粒子数量与上群集粒子总数的比值。显然0对任一给定的α值,MCSG显然是最大参数v0的减函数,当α增大,其减小的趋势变慢(图10.5(a));对任一给定的速率参数v0,MCSG是指数α的增函数,而且小的v0值将导致大的最大连通片MCSG(图10.5(b))。这样,对具有比较大的速率参数v0的群集,当α=0,即Vicsek模型不大可能产生大的簇,而对自适应速度模型仍然可以比较容易的产生大的簇,如

果α足够大的话。考虑对粒子方向作用的随机均匀分布的噪声,

θii(k+1)

φ=i(k+1)e

1ni(k+1)

e

ξi

Γi(k+1)j∈

6e

θij(k)

,k=0,1,2,…图10.6表明了具有幅值为η的噪声对群集的全局序参量Φθ(见公式(10.3))和平均速率系数

Cα的影响。其中平均速率系数Cα定义为:CαΧ1Ni=1

6NφCα(i)=

1N

i=1

α

φ6i

N

(10.11)

可以看到在图10.6(a)中,对大的噪声幅值η和大的指数α,全局序参量Φθ近似沿着同一条直

线下降,这是值得进一步研究的问题之一。

图10.5　(a)MCSG和最大速率参数v0的函数关系;(b)MCSG和幂指数α的函数关系[283]

图10.6　(a)群集的全局序参量Φ的函数关系。对大的噪声幅值η和大的指数α,全局序参量Φθ和噪声幅值ηθ沿着

同一条直线下降;(b)群集的平均速率系数C的增大而逐渐减小到0(取自文献[245])α随着噪声幅值η

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

409

10.1.2.3　异质影响网络中的一致性问题

Vicsek等人提出了一个完全均质的感知模型用于模拟二维空间中的粒子群集现象,在二

维平面上有一群初始随机分布的粒子在运动,每个粒子具有相同的绝对运动速率和不同的方向,且有一个感知半径,在每一时刻以此粒子的位置为圆心,它的感知半径作圆,分布于此圆中的粒子即为该粒子此刻的邻居。该模型假设粒子均为同质的,即每个粒子都具有相同的感知半径。每个粒子下一时刻的方向等于此刻它邻居粒子的方向与自己方向的平均值[271]。在该文中,大量的仿真结果表明,即使不存在集中式控制,只根据邻居间相互影响的规则,所有粒子的方向最终将达到一致。

完全均质的感知网络等价于完全均质的影响网络。然而,在许多社会网络中,主体的影响能力不同,是异质分布的。比如说,一群人正在协商达成一致的观点,每个人在试图劝服他人接受他的观点的同时,也在受其他人观点的影响而调整自己的观点。这群人中比较有权威的人,即影响能力较大的人,他们必定比普通人能够影响更多的人,即更加具有说服力。因此,本文提出了异质影响网络模型,期望通过研究它来更好的反映真实世界中的群集行为。同步现象是一个非常有趣的现象,目前,复杂网络中的同步化能力已得到了广泛的研究,我们在第728已经进行了介绍。文献[209,286]中针对异质网络(节点度、强度、参数分布异质)中的同步问题进行了深入的研究,并给出了一些重要的结论。上述研究都基于一个基本的假设,即网络具有固定的拓扑结构,则绝大部分网络如果异质性越强,网络的同步化能力越差[285]。文献[209]也指出,在Watts2strogatz的小世界网络中,随着网络度分布的异质性增强,网络的同步化能力也增强。

在本章中,我们基于均质Vicsek感知模型提出了异质影响网络模型,并进一步进行了研究和分析。在该模型中,每个智能体的影响能力由它的影响半径来决定,其影响半径根据幂率分布P(r)～r-r(γ∈[2,∞])随机选择,影响半径越大则影响能力越强。随着尺度参数γ的减小,影响半径分布越异质。本章旨在研究该模型的异质性对所有的智能体最终达到方向一致的影响和作用。得出的主要结论是,异质影响网络中的少数领导智能体是最终全局达到方向一致的关键。

如同Vicsek模型[270],我们考虑在一个边长为L,边界周期变化的正方形格子中,有N个运动的智能体,分别编号为1,……,n。每个智能体都有相同的绝对运动速率,但是运动方向

)之内。最初,每个智能体随机分布在格子中。在我们的模型中,假设不同,且分布在[0,2π

每个智能体可以将自己的方向信息传递给那些分布在它影响范围之内的智能体,影响范围即

以该智能体当前时刻位置为圆心,影响半径作圆所覆盖的范围。每个智能体下一时刻的方向等于此刻能够影响它的智能体的方向与自己方向的平均值,即对于智能体i,其方向演化如下:

θi(t+1)=

11+ni(t)

(θi(t))+

j∈Ni(t)

6

θj(t))

(10.12)

其中ni(t)和Nt(t)为影响智能体i的智能体的个数和集合。每个智能体i的位置更新如下:

(10.13)xi(t+1)=xi(t)+vi(t)Δt其中智能体i的速度矢量vi(t)的模值为绝对数率,方向为θi(t)。

γ

我们进一步假设每个智能体的影响半径根据幂率分布P(r)～r-(γ∈[2,∞])随机选

410物理学进展　　27卷

择[274]。当γ=∞时,影响半径的分布完全均质,即所有智能体影响半径相同。在这种情况下,我们的模型就等同于Vicsek模型。随着参数γ的减小,网络变得越来越异质,即有非常大影响半径的智能体减少,而有相对小影响半径的智能体增多。

为了衡量智能体方向达到一致的程度,我们用簇来表示达到相同方向的智能体的集合,并定义当系统演化到稳定时,即所有智能体的方向不再改变时,最大簇中智能体个数与整个系统智能体总数之比为最大簇的相对值S,则0≤S≤1。S越大,则表明方向达到一致的程度越高。当S=1时,所有的智能体均达到方向一致。

为了重点研究由参数γ决定的影响网络的异质性,我们固定系统的某些参数n,L和[288,289]v。设在边长L=50的格子中有n=1250个智能体在移动,且每个智能体都有相同的绝对速率v=0.1。如图10.7(a)所示,在不同尺度参数值γ变化时,最大簇的相对值S关于平均影响半径r的一组曲线。对于γ的某一定值,我们可以发现,S是关于r的一递增函数,这是

-),当由于增大r值就意味着网络连通性能变好。进一步,对于任意γ值,存在一个域值r(γ〈r〉

--)后,系统就能达到所有智能体方向一致,即,S=1。)≥r(γ从图10.7(a)中,可以明显看出r(γ

)≈30;然而,在网络最异质关于γ是一个递增函数。在网络最均质的情况下,即γ=∞时r(∞

--)的十分之一。时(γ=2),r(2)≈3,几乎是r(∞这表明随着网络异质性的增强,智能体更加容

易达到完全一致。从图10.7(a)中还可进一步看出,当平均影响半径〈r〉不够充分大且智能体还没有达到完全一致时,S是关于γ的一个递减函数,这也表明网络的异质性有利于智能体达

-

到方向一致。除了上述考察异质性对达到一致性程度的影响,我们还计算了多智能体系统达到稳定状态所需要的收敛时间。如图10.7(b)所示,随着γ减小或者〈r〉增加,收敛时间递减。这表明网络的异质性或者连通性增强,均有利于系统收敛,即收敛时间缩短。

图10.7　不同尺度参数γ时,最大簇的相对值S(a)和收敛时间(b)随着平均影响半径〈r〉的变化

[345]

趋势。其他系统参数取为L=50,n=1250和v=0.1。仿真取30次平均值

对于上述现象,下面我们进行一些启发式的解释。在异质性强的影响网络中,大部分的智能体具有相对较小的影响半径,但是许多中心智能体有非常大的影响半径。在比较理想的情况下,我们可以近似的认为,这些中心智能体形成了一个连通的完全子图,他们几乎不受其他智能体的影响,然而却可以影响几乎其余所有的智能体。因为网络达到一致的充分条件是,随着网络的演化网络保持强连通[289,290],则那些中心智能体将最终收敛到一个相同的方向θ。假设中心智能体的个数为nh,并且其他每个智能体i至少被一个中心智能体影响,基于这个理想的假设,得到下式:

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

411

θi(t+1)=

其中1≤ni(t)≤nh,有

θi(t+1)-θ=

11+ni(t)1(θ)i(t)+ni(t)θ(10.14)

1+ni(t)

(θ)≤1(θ)i(t)-θi(t)-θ

2

(10.15)

因此,所有智能体的方向最终都将趋于θ。

以上仿真结果和启发式解释表明:在异质网络中,少量的中心智能体在系统的演化过程中起着领导者的作用,且中心智能体之间达到一致是全体智能体收敛到完全一致的关键。

实际上任意实际系统都无法避免噪声的干扰,因此我们进一步考虑噪声给系统带来的影响,并且研究是否异质模型对于噪声干扰具有更好的鲁棒性。我们重写方向更新方程(10.12)如下:

1θ(θθ(10.16)i(t+1)=i(t)+i(t)6θj(t))+Δ1+ni(t)j∈Ni(t)θ其中Δ/2,η/2](η≥0)中均匀随机取的一个i(t)是智能体i所受的随机噪声,其值是从[-η值。在这种情况下,我们不能将所有的智能体精确的分成簇。因此,为了衡量方向一致程度,我ρ=们采用如下序参数ρ:　　　　　　　1nv

i=1

6

n

-

vi(10.17)

从图10.8中,我们可以清楚的看出,ρ是关于η的递减函数,这表明方向一致性的程度随着噪声幅值的增加而递减。另一方面,对于一个固定的η值,ρ是关于γ的递减函数,这意味着随着网络越来越异质,网络抗干扰的鲁棒性越强。

图10.8　不同尺度参数γ时,序参量ρ随着噪声η的变化曲线[345].其他系统参数取为〈r〉=3,n

=100,L=31和v=1。在仿真中,我们取足够长时间之后的ρ值(10这里取大约200个仿真步长)。仿真取100次平均

现在如果我们不仅期望所有的智能体达到方向上的一致,而且希望它们能够收敛到我们

设定的一个方向。很明显,为了达到上述目的,我们需要采用一些控制策略。而异质网络中的中心智能体发挥的领导者的作用促使我们去研究是否仅仅通过控制少数起领导作用的智能体就可以使得所有的智能体一致收敛到期望的方向上。

牵制控制在规则动态网络的时空混沌控制中使用非常普遍,并且最近被应用到无标度动态网络中[285～287]。如果我们希望通过牵制控制一小部分智能体,使得影响网络中的所有智能体一致收敛到一个期望值θ我们比较两种牵制控制策略。若采用特定牵制控制策略,则选p。

412物理学进展　　27卷

择f%具有最大影响半径的智能体进行牵制控制;若采用随机牵制策略,则随机选择f%智能体进行牵制控制。被牵制的智能体以一个固定的期望方向θp,与其余未被牵制的智能体相同的绝对速率在格子中移动。未被牵制的智能体根据式(10.12)更新方向。为了衡量所有智能体最终一致收敛到期望方向的程度,我们定义网络演化到稳定状态时,达到期望值θp的最大簇中智能体个数与整个系统智能体总数之比为最大期望簇的相对值Sd。

我们比较在影响网络中分别采用特定牵制与随机牵制控制策略的效果。仿真中取平均影响半径为〈r〉=6。如图10.9所示,在不同的尺度参数γ下,Sd是关于f%的一个函数。对任意一种牵制控制策略和每个一个γ值,存在一个关键值f,当超过它时,所有的智能体都能够一致收敛到期望的方向(即,Sd=1)。在具有同一较大尺度参数γ的非常均质影响网络中(比如,γ=5),两种牵制控制策略时的关键值f之间的差别相对较小(图10.9(c))。这是由于在均质影响网络中,所有的智能体的影响能力几乎相同。相反,当影响网络十分异质时,采用两种策略关键值f之间的差别很大。比如,在γ=2时,采用随机牵制策略时的关键值f为0.43,大约是采用特定牵制控制时关键值f为0.04的10倍(图10.9(a))。这种大幅度的差别主要是由于影响网络的异质性。正如在文中前面所解释的,如果那些拥有较大半径的中心智能体朝着同一个期望的方向移动,那么几乎所有其他的智能体都会最终朝着期望的方向移动。另一方面,由于网络中大部分的智能体具有较小的影响半径,随机选择到具有较小半径的智能体的概率较大,牵制控制这些智能体几乎不会影响其余大部分的智能体。

)或者特定(Δ)地图10.9　最大期望簇的相对值Sd与关键值f的关系。其中存在,智能体随机(□

被固定在期望方向的两种情形。(a)γ=2;(b)γ=3;(c)γ=5。其他系统参数取为〈r〉=6,n=1250,L=50与v=0.1。仿真取30次平均值[345]

在牵制控制策略中,我们假设被选择的那一小部分智能体被固定在相同的方向。现在假设被选择的智能体被固定在不同的方向,那么是否最终智能体还能一致收敛呢?为了研究这个问题,我们分别采用特定和随机牵制f智能体两种情形。这些被牵制的智能体的方向被随

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

413

机选定并且在演化过程中保持不变。其他未被牵制的智能体的方向根据式(10.12)更新。如图10.10所示,在一个非常异质(γ=2)的影响网络中,当采用特定牵制策略时,最大簇相对大小S迅速减小,从S=1随着f增加而递减:当仅仅2%含有最大影响半径的智能体被固定在不同的方向,则最终最大一致簇包含的智能体仅仅占所有智能体的10%。这就意味着中心智能体达到一致是全局最终达到一致的充分条件。另一方面,当采用随机牵制策略时,最大簇的相对大小S缓慢递减,从S=1随着f增加而递减。即使有10%的随机智能体被固定在不同方向时,最大一致簇也包含了将近占所有智能体40%的智能体。这些结果都表明,异质影响网络中所有智能体达到全局一致的能力当存在随机干扰时很鲁棒,然而在特定攻击时却很脆弱。小部分被牵制的具有较大影响半径的中心智能体不能达到一致将破坏所有智能体达到全局一致,然而小部分被牵制的具有较小影响半径的随机选择智能体不能达到一致并不会影响其余的智能体达到一致。这个结论与最近的研究发现保持一致,即异质网络的连通性能的容错性能好,但是却在恶性攻击时很脆弱[296]。

)或者特定(Δ)地图10.10　最大簇的相对大小S与关键值f的关系,其中存在f智能体随机(□

被固定于随机选择的方向的两种情形。系统其他参数取为〈r〉=6,n=1250,L=50及v=0.1。仿真取30次平均值[345]

10.2　Boid模型的相关介绍及研究10.2.1　Boid模型的介绍

1986年,Reynolds提出了一个在三维空间上用计算机

来模拟群体行为的模型[269]。这个模型要求一群智能体的群体行为满足三条基本规则:(1)分离:与邻域内的智能体避免相撞;(2)聚合:与邻域内的智能体保持紧凑;(3)速度匹配:与邻域内的智能体速度保持一致。智能体的邻域是以智能体所处位置为中心点,围绕该中心点的一定距离范围和角度来构成。图10.11所示的智能体的邻域是以智能体所处位置为圆心的一个圆。这个邻域可以认为是由于有限的知觉感应造成的。每个智能体的运动受到其邻域范围

内智能体的位置和速度的影响。下面我们对Reynolds提

图10.11　智能体i的邻域[22]

出的三条控制规则分别加以详细介绍[269]。

分离规则是使一个智能体能够和它附近的智能体保持一定的分隔距离,以避免发生碰撞,

414物理学进展　　27卷

如图10.12(a)所示。为了能够保持分离,一个智能体首先需要获得邻域内其它智能体的位置信息。对于一个智能体来说,每一个在其邻域范围内的智能体与它之间有一个排斥力,这个排斥力是与它们之间的距离成反比的。每一个智能体所受到的排斥力是它邻域范围内其它智能体对它的排斥力的累加。

聚合规则是使智能体间具有凝聚力,保持队列的紧凑,如图10.12(b)所示。为了能够产生聚合的群体行为,智能体需要获得其邻域内其它智能体的位置信息,并计算出邻域内智能体位置的平均值,由此产生一个作用于该邻域内所有智能体的吸引力,这样使得智能体向平均位置方向运动。

速度匹配规则使各个智能体与其邻域内其他智能体的速度保持一致。如图10.12(c)所示。为了能够实现速度匹配,每个智能体需要获得邻域内其它智能体速度信息,计算出邻域内智能体的平均速度。通过速度匹配,可以使得智能体速度大小和方向与这个邻域内所有智能体速度Ni的平均值保持一致。

图10.12　Reynolds模型的三条控制规则[22,269]

10.2.2　蜂拥控制算法简介

Olfati2Saber提出的第一种蜂拥算法是[263],考虑N智能体在n维的欧式空间中运动,可

以假设第i个智能体的运动方程为:

・

qi=p

・

i

p

i

=ui,i=1,…,N(10.18)

其中qi∈Rn代表智能体i的位置向量,pi∈Rn代表智能体i的速度向量,ui∈Rni代表智能体i的控制输入(加速度)向量。由于每个智能体的感知能力有限,它们只能感应到自己邻域半径内的智能体的状态。在本章中,智能体i的邻域定义为:

(10.19)Ni(t)=j|‖qi-qj‖其中我们假设每个智能体具有相同的邻域半径r。如果我们将智能体作为节点,在节点i与它的邻域中的任一节点之间用一条无向边相连,那么在时刻t由这些节点V={1,…,N}和边E(t)={(i,j)∈V×V:i～j}组成的网络可以用一个无向图G(V,E)或G(t)来表示。对于网络中某个智能体,上一时段在其邻域范围内的一些智能体在某时刻移动到了其邻域之外,或者上一时段不在其邻域范围内的一些智能体在某时刻移动到其邻域之内,因此系统具有切换拓扑结构。

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

415

Olfati2Saber提出的第二个蜂拥算法包括三项:速度一致项、人工势函数梯度项和引导反

馈项。对于智能体i,它的控制输入为[261～266]

ui=fi+fi+fi

g

d

γ

(10.20)

其中第一项fig代表人工势函数V(qi)对位置qi的梯度,用于实现分离和聚合两条规则。其中人工势函数定义为

V(qi)=

j=1,j≠1

6

N

Ψα(‖qij‖σ)=ε1j∈Ni(t),j≠i

6

Ψα(‖r‖σ)+

j∈Ni(t)

6

Ψσ)α(‖qij‖(10.21)

σ=其中qij=qi-qj,‖r‖

2

σ处处可导,我们可以看到范数‖r‖1+ε‖r‖-1和ε>0。

而‖r‖在r=0处是不可导的。这个性质是被用来构造光滑可导的人工势函数。方程中Ψα是一

个非负的光滑可导的人工势函数,它满足:1)当‖qij‖σ→0时,Ψα达到其最大值;2)Ψα在i和j之间的距离‖qij‖σ为某一个期望值时达到最小值,并且‖qij‖σ≥‖r‖σ时,Ψα恒为某个很小的正常数。第二项fid是采用速度一致算法,用于实现智能体之间的速度匹配;第三项fid是引导反馈项,用来实现跟踪虚拟领导者。其中虚拟领导者满足运动方程q・γ=pγ

p・γ=fγ(qγ,pγ)

(10.22)

其中qγ,pγ∈Rn是虚拟领导者的位置和速度向量,它的初始值为(qγ(0),pγ(0))=(qd,pd)。

控制输入方程(10.20)的具体形式为[259]

ui=-j∈Ni(t)

6

ΔqiΨα(‖qj-qi‖σ)+

j∈Ni(t)

6

aij(t)(pj-pi)

(10.23)

+c1(qγ-qi)+c2(pγ-pi),c1,c2>0

其中A(t)=(aij(t))是邻接矩阵,它满足当j=i时,aij(t)=0;当j≠i时aij(t)=ρh(‖qj-qi‖σ/‖σ),其中r‖

1,2∈[0,h)

ρ0.5[1+cos(π(z-h))/(1-h)],z∈[h,1)h(z)=

0,otherwise

和h∈(0,1)。这样我们可以发现:当(j,i)∈E(t)时,aij(t)=aji(t)>0,当(j,i)|E(t)时,aij=0。正常数c1和c2代表引导反馈项的系数。

在初始能量有限的假设情况下,Olfati2Saber证明了在控制算法(10.23)的作用下,Reynolds的三条规则均可以实现,参看文献[260～269]中关于聚合,速度匹配和避免相撞的定理2。

在控制算法(10.23)中,每个智能体都被假设知道虚拟领导者的信息。在本章中我们假设群体中只有一部分智能体知道虚拟领导者的信息。对于智能体i,它的控制输入为

ui=-j∈Ni(t)

6

󰂈qiΨα(‖qj-qi‖σ)+

j∈Ni(t)

6

aij(t)(pj-pi)+hi[c1(qγ-qi)+c2(pγ-pi)]

(10.24)

其中如果智能体i知道虚拟领导者的信息,则hi=1;如果不知道信息,则hi=0。不失一般性,

我们假设前M0(1≤M0416物理学进展　　27卷

运动,因此运动方程(10.22)可以简化为

qγ=pd,qγ=qd

・

・

(10.25)

　　本章我们对整个群体中只有一小部分个体知道信息的情况进行研究(也就是说1≤M0≤

N)。在这种情况下,直观上我们就可以得到一个知道信息的个体在运动的过程中不仅仅受到虚拟领导者的作用,而且还会受到某些不知道信息的个体的影响。因此,对于一个知道引导信息的个体来说,它就不一定会跟踪到虚拟领导者;对其它的个体来说,情况就更难确定。本章我们从理论上证明了不仅仅所有的知道引导信息的智能体能够跟踪虚拟领导者,部分没有引导信息的个体也能跟踪到虚拟领导者,而且通过仿真我们发现即使群体中有少部分个体具有引导信息,绝大多数不知道引导信息的个体也能够跟踪到虚拟领导者。10.2.3　主要的结论和理论分析

我们定义系统的总的能量为系统个体之间总的势能和个体与虚拟领导者之间相对动能之和,它的表达式如下:

1Q=2i=1

6N

(Ui+(pi-pγ)T(pi-pγ))(10.26)

其中

Ui=

j=1,j≠i

6

N

T

Ψα(‖qij‖σ)+hic1(qi-qγ)(qi-qγ)

(10.27)

我们易得Q是一个正半定的方程。我们下面介绍主要的结论:定理一:考虑一个具有N个智能

体的系统,它们的运动方程为(10.18),每个智能体的控制输入为(10.24)。假设初始的能量Q0为一个有限值,则对于M0(1≤M0i)从初始时刻开始,任意一个知道引导信息的智能体和虚拟领导者之间的距离不会超过2Q0/c1。

ii)所有知道引导信息的智能体的速度将会渐近收敛到期望的速度pγ。对群体中的某一

个不知道引导信息的智能体来说,假设在有限的时间T0(≥0)后,在它和群体中某个知道引导信息的智能体之间总可以找到一条有限长度的路径。于是,我们可以得到下面的结论:

iii)从T0时刻开始,这个不知道引导信息的智能体和虚拟领导者之间的距离将不会超过

2Q0/c1+(N-M0)r。

　　iv)这个不知道引导信息的智能体将会渐近收敛到期望的速度pγ。

下面的结论对于所有的智能体都满足:

v)如果系统的初始能量Q0小于(k+1)c3,其中c3=Ψα(0)和k∈Z+,则至多有k对智能体会相撞(k=0能保证所有的智能体都不会相撞)。(证明从略)10.2.4　仿真结果

仿真的对象是N个在平面上移动的智能体。它们的初始位置随机在一个L×L的正方形区域选取,并且满足密度ρ=N/L2=0.1;初始的速度在[0,2]2中随机选取;影响半径取为r=4,智能体之间期望的距离d=3.33,σ范数中ε的值取为0.1,方程ρh(z)中的h取为0.9,

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

417

虚拟领导者的初始位置向量为qd=[10,10]T,速度为pγ=[3,3]T,引导反馈的系数c1和c2均取为0.5。

我们在N个智能体中随机地选取比例为δ的智能体为知道引导信息的智能体。图10.13显示了在N分别取为100、300、500和1000的情况下,能够达到期望速度的智能体所占整个群体的比例p与选取为知道引导信息的智能体的比例δ之间的关系。这里仿真的结果是100次平均后得到的。我们从图10.13中可以看到:对于给定规模的群体,达到期望速度的智能体的比例p是群体中知道引导信息的智能体的比例δ的一个递增函数。而且,群体的规模越大,达到相同的期望比例p所需要的知道引导信

图10.13　知道引导信息的智能体所占的息的智能体的比例δ越小。例如,为了使80%的智　能体达到期望的速度,对于N=100的情况,我们需要使27%的智能体知道引导信息;而对于N=1000

比例和达到期望速度的智能体所占的比例的关系图

的情况,我们只需要使10%的智能体知道引导信息。因此,我们可以推断对于足够大的群体,只需要有很少量的个体知道引导信息,则绝大多数个体可以达到期望的运动。10.3　小结本章综述了我组近期关于Vicsek模型和Boid模型的控制策略的研究工作。主要包括以下三个方面:

(1)推广了经典的常速率Vicsek模型,考虑粒子的速率可以随着其所处的邻居运动方向的极性的情况自适应地进行调整快慢。每一个粒子的运动方向仍然是邻居粒子的运动方向的均值,而速率则是局部极性的一个幂函数。当幂指数α=0,则模型简化为常速率Vicsek模型。在自适应模型中,一些重要的也是困难问题值得进一步研究。例如,指定的收敛概率(或者平均的最大连通片的大小),在什么条件下,我们可以保证存在一个临界指数α,使得对大于该值的所有α,群集具有不低于该收敛概率的演化特征?关于线性化的Vicsek模型的稳定性的问题,现在已有的结果主要是考虑群集在演化过程中粒子之间的连接性特征[275～276]。对自适应速度模型(包括Vicsek模型)更为实用的稳定性判据值得进一步研究。

(2)基于经典的Vicsek模型,并构造了异质影响网络,发现再影响网络中所有智能体达到方向一致的能力随着智能体影响半径分布异质性的增加而提高。实际上,在异质影响网络中,如果可以将起领导作用的那一小部分智能体控制到期望的方向,则所有的智能体就能够全局一致收敛到期望方向。这些结论对控制社会网络和其他人造多智能体网络达到一致起到了一定的作用。

(3)将Olfati2Saber提出的一个蜂拥控制算法推广到更加一般的情况。我们理论上证明了所有知道引导信息的智能体将会按照期望的速度运动。而且,如果一个不具有引导信息的智能体在运动的过程中始终能够直接或者间接的受到具有引导信息的智能体的作用,那么这个智能体也能最终达到期望的速度。通过仿真,我们发现在群体中只需要有很少一部分智能

418物理学进展　　27卷

体知道引导信息,就能使绝大多数智能体最终都能按照期望的速度运动;并且随着群体规模的增加,按照期望的速度运动的智能体所占的比例也会增加。但是需要指出的是这里的虚拟领导者是被假设按照某个固定的速度运动。这样一来,一个未来的研究方向是针对任意运动方式的虚拟领导者和群体中只有少量知道引导信息的智能体的情况,究竟应如何设计一个有效的蜂拥控制算法。

11　复杂网络的演化博弈理论研究

11.0　引言

博弈论(GameTheory),又称对策论,是模拟和分析理性的个体在利益冲突环境下的决策问题,研究个体之间行为的相互影响和相互作用,在经济学、管理学、政治学、军事学、生物学等诸多学科领域具有广泛的实际背景和应用价值。为了刻画个体间利益的冲突对整个系统的影响,人们提出和研究了许多博弈模型。其中,著名的囚徒困境(Prisoner’sDilemma,PD)问题[295]描述了存在于个体理性与集体理性之间的一种矛盾:为了追求自身利益的最大化,参与者的利己抉择使得双方最终都受害,整个系统处于低效状态。雪堆博弈(SnowdriftGame,SG)[296])是囚徒困境的变形和拓展,丰富了个体间冲突和合作的内涵。而以酒吧问题(BarProblem)[297,298]和少数者博弈(MinorityGame,MG)[299]为代表的拥塞博弈是描述现实生活中参与者面对有限资源的竞争行为的模型。

经典的博弈模型中通常假设个体是完全理性的,不会犯任何错误,这显然与现实情况不完全相符。通过引入生物进化论中的自然选择和变异机制,演化博弈理论(EvolutionaryGameTheory)研究了有限理性的个体如何随着时间的推移在不断的重复博弈过程中,通过自适应地学习,来实现最终目标收益达到最优的问题[303～306]。它一方面有助于“精炼”经典博弈中存在的多个均衡———通过个体长期寻优得到的演化稳定策略(EvolutionaryStableStrategy,ESS)对随机扰动是稳定的;另一方面,在博弈演化过程中展现出丰富的动力学行为特性,也为理解现实生活中合作的涌现和合理分配有限的系统资源的内部驱动机制提供了一种有力的理论支持。

个体以及它们之间的关系可以抽象为具有一定拓扑结构的网络:网络中的节点代表个体,边代表两个个体之间的信息传递和作用关系。然而,这种空间结构对合作涌现和资源分配的作用一度被忽视,直到Nowak和May[301]才将其引入到重复囚徒困境博弈的研究中去,人们开始关注空间博弈,从群体组织关系角度来研究网络结构对策略演化和合作涌现的影响。早期的空间博弈都将注意力集中在规则格子与个体博弈行为的关系。随着近年来复杂网络理论的兴起[13,15,19,20,241],人们对自然界的组织关系有了更高层次的理解:我们生活的世界组织既不是完全规则的,又不是完全随机的,而是介于两者之间具有小世界[1]和无标度[3]等特性。这些网络特性的发现提供了新的理论框架和研究方法,让人们开始考虑从整个网络连接拓扑结构的系统角度来理解自然界中广泛存在的合作现象和资源分配。本章以少数者博弈、囚徒困境和雪堆博弈为例,结合我们的工作进展,综述了复杂网络对博弈演化作用方面所取得的最新的研究成果。

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

419

11.1　少数者博弈11.1.1　少数者博弈简介

少数者博弈模型是由Challet和张翼成于1997年提出的[298]。它假设在一个系统中有N(N为奇数)个参与者,在某一时刻各自独立地在两个策略中做出选择,参与人数少的策略获胜。少数者博弈模型实际上是对美国著名经济学家提出的酒吧问题[297,298]的一种抽象和简化。酒吧问题研究的是一群生活在圣塔菲的人们在周四晚上是否去附近的ElFarol酒吧的决策问题:每周四晚上这个酒吧都会有优雅的爱尔兰音乐演奏,然而如果去的人数过多,超过了酒吧所能容纳的人数(阈值c),酒吧就会变得嘈杂拥挤,人们也无法悠闲地欣赏音乐。因此人们需要根据过去的公共信息来对当晚去酒吧的人数做预测,以决定自己究竟是去酒吧还是留在家里。酒吧问题和少数者博弈模型都反映了社会经济活动中众多千差万别的参与者对有限资源竞争的基本特征,其思想是金融市场中的普遍原则———少数人获胜。

演化少数者博弈(EvolutionaryMinorityGame,EMG)将进化论与少数者博弈结合在一起,发现通过学习过去的公共历史信息,参与者的平均收益得到了提高[297]。在EMG模型中,对于某一轮博弈,参与者根据他记忆中保存的公共历史信息来独立地决策本轮自己是加入“1”组还是“0”组;当所有人都做出选择后,进入人数少的一组的人为获胜者,进入人数多的一组的人为失败者。因此,EMG模型初始化时会随机赋予每个参与者一个基因值p(0≤p≤1),在每一轮博弈中,每个参与者以概率p及以概率1-p选择失败方的策略。博弈后获胜者的收益加S,而失败者的收益减1,因此S也被称为奖惩比。随着时间的演化,参与者的累积收益将发生变化,当低于阈值d(d<0),它需要修改其策略p值来改变窘境。新的p是在以初始的p值为中心,宽度为R的区间里随机选择的,而这时该参与人的财富也重置为零。人们通过对EMG模型的研究发现一个有趣的现象:当S≥1时一个相互间竞争的人群最终总是趋向于分离成为具有两种相反的极端行为的人群[303]。如图11.1(a)所示,无论初始分布如何,参与者的策略基因值p的分布P(p)最终都会演化为一个对称的“∪”形分布。这意味着为了在竞争社会中生存,参与者的行为最终会走向极端:要么始终遵循基本策略(p=1),要么始终反其道而行之(p=0)。但是,最近Hod和Nakar研究了S<1情况下的策略分布,如图11.1(a)所示,P(p)的分布可以形成“M”形(S=0.992),也可以形成“∩”形分布(S=0.971),这说明随着奖惩比S的减小,参与者采取的策略从极端转向中庸。进一步,陈侃和汪秉宏研究发现P(p)的这种转变过程不但与奖惩比S有关,而且与参考者数目N、阈值d相关。如图11.1(b)所示,S=0.971时随着参与者人数N的增加,策略分布P(p)也会从“∪”形过渡为“∩”形,并对人群分布相变问题作了统计力学的理论解释。

在演化少数者博弈模型中,参与者仅依据公共的信息,并凭借自己的过去经验参与博弈,参与者之间并没有直接的信息交流或相互作用。随着近年来复杂网络研究的蓬勃发展,复杂网络上的少数者博弈问题的研究也开始受到关注。引入空间概念后,一个参与者i可以知道他的邻居的财富以及他们的基因值p[298]。由此进化变异规则可以为:当某一时步一个人的财富低于阈值d(d<0)时,他将财富重置为零,并按照如下规则修改p:新的概率p′是在以pbest为中心、宽度为R的区间内重新按均匀分布随机选取的,即

420物理学进展　　27卷

图11.1　(a)不同奖惩比S情况下,基本EMG模型的策略分布P(p)。其中N=10001,d=-4(取自

文献[HodS.,NakarE.,Phys.Rev.Lett.2002,88:238702]);(b)奖惩比S=0.971时策略分布P(p)随着参与者数目N的变化,其中d=-4(取自文献[ChenK,WangBH,YuanB,Phys.Rev.E.2004,69:025102];Intern,J.Modern.Phys.B,2004,18:238722393)

p→p′∈[pbest-R/2,pbest+R/2](11.1)

其中pbest是i的ki个邻居中财富最高的邻居的概率。R(0≤R≤2)可以看成是参与者的新概率p′和他的最佳邻居的概率pbest之间关联程度的度量:R越小,相关性越大;R=0时,该参与者直接用最佳邻居的概率替换自己的概率;R=2时,参与者的新概率p′和他的最佳邻居的概率pbest之间几乎没有什么关联。这里需要保证0≤p≤1。

为了衡量系统资源的使用效率,定义:

σ=

2

1T

t=1

6

T

(A(t)-A)2,

(11.2)

其中A(t)表示t时步选择“1”组的人数,A=1/T

22“1”组的平均人数。因此σ越小,资源利用率就越高。在最理想的情况下,σ等于0.25,即(N-1)/2个参与者在少数方,其余(N+1)/2个参与者在多数方。而在没有任何信息可以依据的情况下,每个参与者完全随机独立地选择自己的行动———这被称为随机选择博弈(Random

2

ChoiceGame,RCG),由于参与者选择每个组的概率都为0.5,系统的方差σ等于0.25N。

近期研究表明,在最近邻网络上的参与者通过模仿邻居的行为,可以促进参与者之间的协调性[308,309]。而对于Kauffman随机网络上的演化少数者博弈研究表明,当网络中所有节点的入度都为2时系统性能达到最优,资源利用率最高[310]。小世界网络的随机重连概率对系统性能也有重要影响,尤其是当随机重连概率等于0.3时系统协调性达到最优[311]。而对于无标度网络中参与者的度与它的成功率之间的相关性研究,也表明这二者之间存在正的相关性[310]。本节后面主要介绍在规则格子、小世界网络、无标度网络上EMG模型的动态演化行为[311],以及在随机Kauffman网络上修正演化少数者博弈模型的研究结果[309]。

6

Tt=1

A(t)～N/2表示演化时间T步内选择

11.1.2　星形网络

我们研究了星形网络上的演化少数者博弈模型,首先固定奖惩比参数S=1来观察参与者的策略分布P(p)(如图11.2(a)所示)。对于任意的初始策略分布,参与者的稳态策略分布P(p)都是关于中心(p=1/2)对称的分布。因此星型网络结构改变了EMG模型在收益

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

421

函数对称(S=1)时的策略分布形状,与图11.1(a)基本EMG模型的“自分离现象”完全相反。这是由于参与者基于星型网络进行交互作用时,所有的参与者都模仿中心参与者的策略,因此发生了群集行为。而当群集作用导致成功率下降后,人们的策略开始变得摇摆不定,会更频繁地改变自己的选择。此时采用谨慎小心的中庸策略的成功率往往高于极端策略的成功率,从而形成了如图11.2(a)所示的p=1/2处的峰值。当S远小于1时,如图11.2(b)所示,人群分布仍会出现向中庸人群靠拢的趋势,因此策略分布在p=1/2处形成了峰值———这与基本

2

EMG模型在S<1时是相似的,我们也观察到参数S对约化方差σ/N的影响:当S<1.1时,随着S的增加约化方差迅速减少,但仍然高于RCG的约化方差(等于0.25);当S>1.1后,随S的增加约化方差基本保持不变,略小于RCG的约化方差。

(a)S=1;(b)S=0.9(取自文献[298])图11.2　星型网络上参与者的策略分布P(p),N=1001,R=0.05。

2

与S对约化方差σ/N的影响相反,在R较小时,系统有较大的约化方差,并随着R的增加缓慢减少;而当R>0.5后,约化方差随着R

的增大而迅速减小,这是因为参与者的新概率与被模仿的中心点参与者的概率之间的相关性越来越小,避免了群集行为的产生。因此随着R的增大,如图11.3所示,策略分布P(p)逐渐由“∩”形分布过渡为“∪”形分布,所以星形网络的策略分布对区间宽度变化比较敏感。11.1.3　小世界网络

　图11.3　星型网络上不同的R值对应的参与者的

策略分布P(p)(取自文献[303])

我们来进一步观察小世界网络上的少数者

博弈行为。当S取值不同时WS小世界网络上稳态的策略分布情况如图11.4所示:当S=1时,对应不同的小世界网络的重连概率psw,WS小世界网络拓扑结构没有改变EMG模型的稳态策略分布,参与者都是趋向于自组织分离成极端行为的相反人群。图11.4(b)给出了S=0.9时,相应的p值分布。当S远小于1时,WS小世界网络上参与者的稳态策略分布演化为对称的、在中心峰化的“∩”形分布。这意味着在困难时期(惩罚大于奖励的时候)人们倾向

422物理学进展　　27卷

于采用和多数人一致的决定。

图11.4　WS小世界网络的重连概率为psw=0,0.01,0.1,1时,参与者的策略分布P(p),(a)S=

1;(b)S=09(取自文献[303])

当小世界网络中重连概率psw取不同值时(如图11.52

所示),通过研究约化方差σ/N随奖惩比S的变化关系可2

以发现:当S<1时,σ/N随着奖励的增加缓慢减小;当惩

2罚等于奖励(S=1)时,σ/N出现了一个最小值;当惩罚小

2

于奖励(S>1)时,σ/N随着奖励的增加缓慢增大,最后基

2

本保持不变,此时σ/N的值小于RCG的约化方差0.25。

2

而约化方差σ/N随区间宽度R的变化关系表明:随着psw

2越大,对应的σ/N也越大。因此最近邻耦合网络(psw=0)

22

上的σ/N最小。同时,小世界网络上最大的σ/N仍小于　图11.5　WS小世界网络的重连概率为

psw=0,0.01,0.1,1时,约化参与者随机选择时的约化方差0.25。

2

方差σ/N随奖惩比S的变化

WS小世界网络的重连概率psw也影响着整个网

关系(取自文献[303])

络的合作行为,如图11.6所示:基于最近邻网络作用

的人群能够相互协调彼此合作,整个社会达到了较高的

2

资源优化配置,σ/N很小。当psw较小时(022

01),小世界网络上的σ/N与最近邻网络上的σ/N差别不是很大;随着psw的进一步增大(psw≥0.01),参与

2

者之间的协调性减弱,σ/N迅速增加。因此系统的约

2

化方差σ/N的变化与网络的群聚系数相关:网络的群聚系数越大,系统的约化方差越小。11.1.4　无标度网络

2

　图11.6　约化方差σ/N随WS小世

基于无标度网络的少数者演化博弈行为与小世界网络界网络的重连概率psw的变

化关系,S=1,R=005(取时的情况相似。当奖惩比S=1时,参与者趋向于自组织分

自文献[303])离成两种极端行为的相反人群;而当奖惩比S=0.9时,P(p)

在中心(p=1/2)处形成峰值,对0和1附近的极端值p,P(p)=0。同时,我们发现无标度网络上的

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

423

2约化方差σ/N随S的变化趋势也与小世界网络上的情况相似。

在BA无标度网络中,绝大部分节点的度相对很低,但存在少量的度相对很高的节点,即节点的度分布具有异质性。参与者的度与他们的收益之间的相关性如图11.7所示,其中每个点代表一个参与者,实线为所有参与者成功率的平均值。可以发现参与者的度与他们的成功率之间存在正相关,度大的参与者的表现要好于那些度小的参与者。因此参与者获得的信息

越多,他们的收益就越高。

为了研究无标度网络上系统的动态行为与群聚系数的关系,图11.8显示了可变群聚系数

2

的无标度网络[48]上的约化方差σ/N随区间宽度R的变化曲线,可以发现,随着群聚系数的增加,约化方差越来越小,参与者之间的协调性越来越好。因此,EMG模型的动态性能与网络的群聚特性相关,网络群聚系数越大,越容易促进人群与反人群的形成,资源利用率越高。

图11.7　BA无标度网络上参与者的成功率与　图11.8　可调群聚系数的无标度网络中,约化方差

2

σ度的关系,S=1,R=005　　　　　　　　　/N随R的变化关系,S=1(取自文献[302])

11.1.5　随机Kauffman网络上修正演化少数者博弈模型

随机Kauffman网络是由入度相同的节点组成的有向网络[310],这类网络在研究自然界中的基因调控机理[311],分析社会系统中的规范形成[308]和设计存储器等工程应用[313]中具有重要的地位。随机Kauffman网络的生成过程如下:首先指定N个节点,然后每个节点i随机选择其它K个节点j(j=1,…,K)相连,连边的方向是由节点j指向节点i,表示节点i可以获

得节点j的信息。针对随机Kauffman网络,我们提出了一种改进的演化少数者博弈(Modi2fiedEvolutionaryMinorityGame,MEMG)模型[304],参与者在局域范围内模仿自己邻居的行

为,而不需要其他历史信息。

对于有两种选择的MEMG模型,假设共有N(N为奇数)个参与者,在某时刻他们以概率

p选择“1”组,而以概率1-p选择“0”组。这里p的意义与基本EMG模型不同:在基本EMG

模型中[297],参与者以概率p选择公共信息中获胜方的行为,以概率1-p选择失败方的行为。而MEMG模型在博弈过程中没有使用任何公共信息,仅仅依据局域信息进行决策。当所有人都做出选择后,进入人数少的一组的人为获胜者,每人的财富加一分;进入人数多的一组的人为失败者,每人的财富减一分。假设每个参与者知道K个邻居的财富以及他们的基因值。一个参与者在t时步的最佳邻居是t-1时步他的K个邻居中财富最多的一个,若有多个具有相同财富的最佳邻居,则随机地在其中选择一个。当被选中的最佳邻居的财富比自己的财富

424物理学进展　　27卷

多时,该参与者就会修改自己的概率(基因值)p,修改规则与EMG相同,因此该模型通过模仿和修正达到了演化的目的。进一步地,MEMG模型可以自然推广到多选择情况,即每个参与者可以从Q个组中选择其中一个加入。

对于两个选择的MEMG模型,研究表明,K取1、2、4、8时,类似于基本EMG模型S=1的情况,系统通过自组织都会形成了两种截然不同的群体,分别采用极端的概率值:p=0或

2

p=1。同时,我们发现,策略分布P(p)对K值存在依赖性。图11.9给出了约化方差σ/N与K的关系曲线,结果显示约化方差在K=2处取得了最小值,并且不依赖于N和R的改变。该结果与文献[310]的结果是一致的,即MG和MEMG模型都是在K=2处取得了最佳的合作效果。而K=2的系统的高度适应性与Kauffman网络本身的内在特性是密切相关的,网络的行为在K=2时处于有序和混沌的临界形态[306]。同时,我们发现系统约化方差的大小也与区间

2σ宽度R有关。/N与R的函数关系不是单调的:R2

很小时,σ/N随着R的增加迅速减少,在R≈0.045。2处有一个最小值;然后随着R的继续增加,σ/N也迅速增加。因此参与者经过模仿和修正p值后,新的概率和最佳邻居的概率之间的关联程度有一个最优值。我们还比较了MG、EMG和MEMG3种模型在相同参数设置下的约化方差变化,发现随机Kauffman　2

图11.9　约化方差σ/N随网络平均连接度

2

网络中MEMG模型的约化方差σ/N不但小于EMGK的变化关系(取自文献[304])模型,而且小于MG模型的情况。这说明基于

MEMG模型的模仿和修正机制,随机Kauffman网络中犹豫不决的参与者几乎消失,形成了大小

几乎相等的两组人群,从而改善了参与者之间的协作效应,使系统的整体收益得到了显著提高。而对于多选择的MEMG模型,其约化方差也远小于RCG模型的约化方差,所以模仿机制的引入可以明显增强系统的协作,从而取得了最佳的资源配置效果,提高系统的整体收益。11.2　囚徒困境1.2.1　囚徒困境简介

囚徒困境(Prisoner’sDilemma)作为一个经典的博弈模型受到广泛关注。这个博弈模型假设两个小偷合伙作案时被捕,分别关在不同的屋子里,如果双方都拒绝承认同伴的罪行,则由于证据不足两人都会被轻判R;为此,警方设计了一个机制:如果一方出卖同伴,而另一方保持忠诚,则背叛者将无罪释放,收益为T,坚持忠诚的合作方将被重判S;如果双方都背叛了对方,则双方都会被判刑P。这里假设收益参数满足下面的条件:T>R>P>S。对每个参与者来说,如果对手选择合作,则他也选择合作得到的收益R小于他选择背叛得到的收益T;如果对手选择背叛,则他选择合作得到的收益S仍小于他选择背叛得到的收益P。因此,无论对手采取哪种策略,自己的最佳策略就是背叛。双方都选择背叛是囚徒困境的唯一纳什均衡。然而,两个人同时选择背叛所取得的平均收益要低于两个人同时选择合作取得的平均收益,在这种情况下,理性参与者面临着两难的困境。但是,自然界中广泛存在的合作现

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

425

象———从单细胞生物的协同工作到人类的无私奉献的行为———说明,还有其他的动力学机制激励自私的个体认识到合作的重要性。为了揭示这种潜在的演化机制,Axelrod提出了“针锋相对(TitForTat,TFT)”演化规则[294];Nowak等人则研究了“去输存赢(Win2Stay2Lose2Shift,WSLS)”策略[301,318],以改进囚徒困境中的两难结局。然而,这些研究忽视了一个重要的方面:群体的组织结构对合作涌现和系统性能的影响。而日益深入的复杂网络研究为研究空间结构对合作现象的涌现作用提供了强有力的理论基础和实验手段[317,318,330]。这一节将介绍近年来这方面取得的一些主要研究成果。11.2.2　规则格子上囚徒困境

Nowak和May首先将空间结构引入囚徒困境[301,318],研究了二维方格网络上的重复囚徒困境博弈,他们发现,合作现象能够在规则格子上涌现。在网络上,每个个体占据一个节点,个体仅与直接连接的邻居根据囚徒困境进行博弈,并累计收益;然后在更新策略时,个体与它的邻居比较本轮的收益,取收益最高者的策略作为下一轮博弈的策略,直到网络进入稳定状态为止。通常设定囚徒困境中的收益参数为T=b>1,R=1,P=S=0。Nowak发现,引入空间结构后,通过演化,当b在一定范围之内变化(1≤bν2,如图11.10(b)所示),合作者可以通过相互结成紧密的簇来抵御背叛者的入侵。虽然从微观角度来看,这种合作簇并不固定,其形状会随时间的改变而改变,但它并不会消亡;从宏观角度观察,稳定状态时系统的合作者的比例会趋于稳定(如图11.10(a)显示),它被称为合作频率(fractionofcooperators),是衡量系统合作涌现程度的重要指标。

图11.10　在方格上进行囚徒困境博弈。(a)合作频率随着b的变化;(b)b接近于合作湮灭时的斑图[318]

11.2.3　无标度网络

Santos研究了无标度网络对合作涌现的作用[319～321]。初始化时,网络中所有节点随机地

选择合作或者背叛;在一轮博弈中,网络中每个节点x与它的邻居进行一次囚徒困境博弈收益累计为Px;当本轮博弈结束时,每个节点将更新自己的策略:一个节点x在它的kx个邻居中随机地选择y并比较二者的本轮收益,如果Px>Py,下一轮博弈中x保持自己的策略不变,反之x采用概率:

p=(Py-Px)/[kmax(max(T,R)-min(S,P)]

(3)

426物理学进展　　27卷

图11.11　网络模型对囚徒困境中合作涌现的影响。网络平均度为4,规模为10000(取自文献[316])

采取y本轮的策略作为其下一轮的策略。其中kmax=max(kx,ky)是这两节点中的最大度。可以发现,网络中具有较大度的中心节点具有较高的累计收益,因此对网络中的合作涌现具有至关重要的影响。Santos的研究表明(如图11.11所示),具有较大异质性的无标度网络比ER随机网络和规则格子更容易涌现合作。进而,Santos等发现当中心节点之间保持较好连通性时,网络进入稳定状态后,中心节点都会转变为合作者[320,321]。这是由于初始合作的中心节点会坚持自己的策略不变,并带动其周围越来越多的邻居成为合作;而初始背叛的中心节点虽然短期内会从它周围的采取合作策略的邻居处取得较高收益,然而随着邻居模仿它的行为,其收益会下降,低于与它相邻的合作中心节点的收益,所以随着时间的演化,背叛的中心节点有学习合作的中心节点行为的趋势,最终导致了中心节点都成为了合作者,进而这些中心节点带动了大量小度节点成为合作者。因此,与其它类型网络相比较,囚徒间的合作行为更容易在无标度网络中传播。

最近,Gomez2Gardenes等人[322]将进入稳定状态后节点根据状态分为三类:随时间演化保持合作的节点、保持背叛的节点和不断改变策略的节点。他们发现与ER随机图中保持合作的节点零散分布不同,在无标度网络中,保持合作的节点会结成一个连通的簇,阻断背叛者的入侵,而保持背叛的节点零散地分布在无标度网络中。别外文献针对网络的异质结构提出了一种优先选择机制[323],可以有效促进合作在网络中的涌现。这些研究加深了对现实网络中合作涌现机制的理解。11.3　雪堆博弈11.3.1　雪堆博弈简介

雪堆博弈(SnowdriftGame)又称为鹰鸽博弈(HawksandDovesGame)或者小鸡博弈(ChickenGame),是另一类两人对称博弈模型,描述了两个人相遇时是彼此合作共同受益,还是彼此欺骗来相互报复,它揭示了个体理性和群体理性的矛盾对立。可以这样来描述雪堆博弈:在一个风雪交加的夜晚,两人相向而来,被一个雪堆所阻,假设铲除这个雪堆使道路通畅需要的代价为c,如果道路通畅则带给每个人的好处量化为b>1。如果两人一齐动手铲雪,则他们的收益为R=b-c/2;如果只有一人铲雪,虽然两个人都可以回家,但是背叛者逃避了劳动,它的收益为T=b;而合作者的收益为S=b-c;如果两人都被雪堆挡住而无法回家,

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

427

他们的收益都为P=0。这里假设收益参数满足下面的条件:T>R>S>P。雪堆模型与囚徒困境不同的是,遇到背叛者时合作者的收益高于双方相互背叛的收益。因此一个人的最佳策略取决于对手的策略:如果对手选择合作,他的最佳策略是背叛;反过来,如果对手选择背叛,那么他的最佳策略是合作。这样合作在系统中不会消亡,而与囚徒困境相比,合作更容易在雪堆博弈中涌现。雪堆博弈有两个纯策略纳什均衡:(合作,背叛)和(背叛,合作)。同时存在一个演化稳定策略,根据模仿者动态[324],在进化的稳定状态中,合作者所占的比例是1-r,其中r=c/(2b-c)称为相互合作的损益比(Cost2to2benefitRatio),并且它对于随机扰动是稳定的。需要指出的是,雪堆博弈在稳定状态时的平均收益仍然低于两个人同时选择合作时的平均收益,所以它仍体现了理性参与者的两难困境[316]。11.3.2　规则网络上雪堆博弈

自从Nowak和May将网格结构引入囚徒困境,随后的众多研究表明空间结构能够促进和维系囚徒困境中合作的产生和稳定[325～327]。一直以来囚徒困境被认为是博弈论研究合作和背叛的经典模型,雪堆博弈并没有受到太多的关注,直到2004年Hauert和Doebeli在Nature上发表工作[325,326]使人们开始把目光投向了雪堆博弈。他们研究了二维规则网格上雪堆博弈的合作演化,却得出了与囚徒困境不一致的结论:如图11.12(a)所示,规则格子上雪堆博弈的合作频率低于模仿者动态下的演化稳定策略———说明空间结构抑制了合作的产生,它与人们的通常认识是截然相反的。这是因为与囚徒困境的斑图不同,在雪堆博弈中合作者更容易聚成丝状簇,如图11.12(b)所示,这导致了当损益比r较高时,背叛者更容易侵入,使系统合作频率的下降,这是雪堆博弈与囚徒困境在合作演化上的本质区别。

图11.12　在方格上进行雪堆博弈,(a)合作频率随着r的变化;(b)r接近于合作湮灭时的斑图(取自文献[316])

11.3.3　基于距离的空间小世界网络上的雪堆博弈

我们也研究了基于空间距离的小世界网络上的雪堆博弈[327]。对于一个二维方格网格,受Kleinberg小世界模型的启发[91,100,326],我们依“按距离概率”为网络加边:对于一个给定的常数q≥0,以概率psw给网络中的每个节点u随机地增加q个长程连接。选择节点v与节点u建立长程连接的概率

∏满足如下关系:

v

428物理学进展　　27卷

∏

v

=

其中指数α≥0,d(u,v)是点对(u,v)之间的欧几里德距离。因此随着α的增加,每个节点都趋向于和邻近的节点相连,两个相隔比较远的节点之间建立长程连接的可能性越来越小。基于距离的空间小世界模型反映了现实生活中朋友关系网络的一种特性:住得比较近的邻居和在一起工作的同事之间更容易成为朋友;另一方面,住的比较远的,甚至是在不同国度的两个人相遇的几率比较小,因此彼此成为朋友的可能性也较小。

首先我们研究指数α=0时与距离无关的小世界网络对博弈行为的影响。按照文献[319]提出的演化规则(公式(3)),参与者在每轮博弈结束后更新自己的策略。图11.13画出了雪堆博弈的合作频率与损益比r的关系曲线,其中小世界网络的连边概率分别取为psw=0,0.01,0.1,0.5,1,q=4。与规则网络(psw=0)相比,小世界拓扑结构下(psw>0)人群的合作频率得到了很大的提高,即使在r取值比较大的时候也能够促进合作演化。因此小世界模型引入的长程连接扩大了合作者的优势,有利于合作社会的形成与繁荣。而且由图11.13的插图可以明显看出合作频率对psw的依赖是非单调的。

进而我们研究了对应不同的r,合作频率随psw变化的关系。图11.14中显示对应不同的r,合作频率都是随着概率psw的增加而迅速增加,直至psw增加到psw≈0.2,此时与psw=0相比,系统的合作频率得到了显著的提高。这意味着整体的合作行为对底层小世界网络的平均路径长度L比较敏感。对于损益比r=0.3和r=0.4,合作频率在psw≈0.2时出现了一个峰值,表明在该点处人群取得了最佳的合作效果。随着r的增大,psw≈0.2处的峰值越来越趋于饱和,当psw>0.2时,合作频率基本保持不变,因为该阶段L的变化比较平缓。

6

-α

[d(u,v)]-α

[d(u,j)]j

(4)

图11.13　雪堆博弈的合作频率随损益比r的变化关系[328],其中指数α=0,概率psw=0,0.01,0.1,0.5,1

我们再来观察指数α>0时,基于距离的空间小世界网络对合作演化动态的影响。随着α

的增大,网络群聚特性快速增长以及平均路径长度变大。我们发现只有当r比较小时(r<0.3),合作者的比例才不会受到α变化的影响,主要是因为利益较高,损失较小,所以合作相对

背叛来说占有绝对的优势。然而,当r比较大时,具有较大α的空间结构对合作动态产生了重要的影响。合作频率随α变化的关系曲线如图11.15所示,当r<0.4时,参与者的利益虽然降低了,但相对来说还是比较高,因此不同的概率psw产生了不同的演化行为。当r≥0.4时,

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

429

图11.14　雪堆博弈的合作频率随概率psw随不同概率PSW的变化关系[328],其中指数α=0,损益比

r=0.3,0.4,0.5,0.6,0.7,0.8

利益相对来说比较小,此时网络的群聚特性使合作策略容易被背叛策略取代,因此这种情况会抑制合作,对背叛者更有利。随着α的增加,网络中长程连接越来越稀疏,短程连接越来越密集,因此长程连接模式在r较大时有利于整体合作水平的提高。在其他的网络规模下进行相同的仿真研究表明,系统性能在r=0.4处都存在这样一个转折点,因此该结论不依赖网络规模的变化。

与规则网格相比,我们的研究表明距离无关的小世界网络有利于整体合作水平的提高。然而当损益比比较大时,随着长程连接的减少和短程连接的增加,距离相关的小世界网络结构趋向于抑制雪堆博弈的合作行为演化。

430物理学进展　　27卷

图11.15　雪堆博弈的合作频率随指数α的变化关系[323],其中概率psw=0.01,0.2,1,损益

比r=0.3,0.6。(a)2(c))对应r=0.3,(d)2(f)对应r=0.6

11.3.4　标度网络

Santos将研究无标度网络上囚徒困境的方法移植到雪堆博弈上[319],观察到类似于图11.11的现象,这说明无标特性同样有利于雪堆博弈中合作的涌现。特别的,通过对小规模网

络(128个节点)进行仿真,弱化了影响合作涌现的无标度网络其它统计学特性,着重突出了节点度的异质性的因素。再次验证了关于异质因素促进合作涌现的一般性结论,指出无标度网络为研究演化博弈理论提供了统一的理论框架。

基于少数者博弈,王文旭等[323]提出了一种基于记忆的雪堆博弈模型,在无标度网络上,小度节点为了追求自身的高收益,会选择与中心大度节点相反的策略,这导致了合作频率随着损益比的变化而震荡。因此,中心节点对网络中合作涌现的程度具有至关重要的作用。11.4　结论

本章以少数者博弈、囚徒困境和雪堆博弈为例,从复杂网络的不同拓扑结构对合作涌现和资源分配的作用角度,评述了近年所取得的包括我组在内一些新的研究成果。研究表明,群聚系

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

431

数、网络异质性等结构特性对个体的博弈行为起到了至关重要的影响。Nowak近期指出[329]:自然合作(NaturalCooperation)与自然选择,与遗传变异居于同等地位,是演化理论的第三条基本规则。而博弈行为对网络结构演化的作用也同样值得关注[325～327]。我们相信,随着演化博弈动力学行为与复杂网络之间的关系的逐渐清晰,空间博弈必定会推动复杂网络中其他领域的发展。

12　网络科学面临的挑战和展望

全文通过上述11章,从不同课题和多角度评述了正在国内外迅猛发展的新兴的网络科学的若干重要研究进展,特别是,重点总结了我们正在开展的复杂网络的国家自然科学基金重点项目所取得的主要成果。由于网络科学日新月异,层出不穷,涉及非常广泛的交叉课题,与众多领域的密切相关。尽管,网络科学的诞生仅仅10个年头,但是她所取得的辉煌成果已经极大地改变和丰富了人们对复杂客观世界的认识,揭示了前所未有的理论和技术问题,展现了巨大的应用前景。网络科学不仅适用于自然网络,而且适于人造网络,进一步探索复杂网络的规律,无论对自然界,还是对人类社会,对于国防和经济都具有极端重要性和长远的意义。

鉴于我们知识和篇幅有限,无疑,还有一些重要的课题,尤其与应用有关的研究课题未能完全涵盖在内,这里无法详细介绍。为了稍微弥补不足,这里简略提及值得关心和正在研究中的几个课题[332～355]。关于复杂网络上的拥塞问题[336～345]。拥塞现象是发生在互联网、通信网络和交通网络等复杂网络上的一种典型的普遍的动态行为,与实际应用关系密切。网络科学对这个问题的研究和以往的不同之处在于,这里更关注网络的结构特性对网络的搜索、拥塞、路由优化策略和动力学行为的影响。人们寻找是否可以利用网络的局部或者全局拓扑特性来提高搜索以及路由策略的有效性?路由策略是否会导致网络拓扑结构发生变化从而反过来影响算法的效率?如何设计对网络拓扑结构变化具有自适应能力的路由策略?局部优化策略对无标度/小世界网络的结构的形成和动力学行为有何影响?等等。为此,我组基于各种复杂网络模型(小世界网络和无标度网络),分析了不同网络结构上拥塞产生的原因及其控制策略,探讨了网络结构和其上发生的动态行为之间的相互关系,提出了改进的路由策略来提高网络的吞吐量以及传输性能等动态指标。基于BA无标度网络模型,按照节点在网络中的重要程度将其分类,并通过定义相应的动态过程及性能指标,研究了网络节点自身容量或者处理速度以及网络的无标度结构特性对拥塞的影响,提出了相应的控制策略来改善网络的拥塞,研究发现控制策略是否有效和网络的拓扑结构紧密相关,相同的策略在不同网络结构中的作用会有明显不同。特别是,不需要对整个网络施加控制作用,而仅需要对一些最关键的节点加以控制,就可以得到类似控制所有节点所产生的控制效果。我组还提出了一种结合最短路径路由和局部拥塞检测的改进路由算法,并在具有可变群聚系数的无标度网络上进行了仿真。研究表明:可以得到一个最优参数使得该路由算法具有最好的性能,并且随着拥塞的加剧,最优参数值随之增加。和最短路径路由算法相比,该算法在拥塞加剧的情形下性能提升更加明显,但是最优参数值几乎不随群聚系数的变化而变化。进一步研究了可以感知局部拥塞信息的路由算法在不同网络拓扑结构下的性能,发现在各种拓扑结构下,该路由算法均可以显著的提高网络吞吐量,并且路由算法的提高程度和网络的异质性有密切关系。复杂网络上的拥塞问题的研究仍然有不少课题

432物理学进展　　27卷

值得研究。

关于复杂网络上的传播问题[346～349]。以前研究的网络传播主要基于规则网络。现在则是基于小世界和无标度特性的复杂网络的统计性质,重新审理和发现的网络上的传播/演化的规律。主要问题是:如何控制具有小世界和无标度特性的复杂网络上的传播动力学?计算机病毒和各种传染病是如何传播的?包括:艾滋病、非典和禽流感等在复杂网络上的传播问题。如何有效应对这类传染病给人类造成巨大的威胁?这些是当前及今后国内外急迫需要研究和解决的应用问题。我们通过仿真分析了网络上即时通讯病毒的传播特性,发现其传播速度非常迅速,并且在无标度网络上比在其他网络拓扑结构下传播得更快.单独增强用户的安全意识只能减少病毒的感染范围,却不能减缓病毒的传播速度。于是,我们结合复杂网络理论提出了监控和防御即时通讯病毒的一条思路,引入杀毒软件的病毒库的概念,使即时通讯软件能对已知病毒免疫。根据病毒的传播特点和即时通讯网络的无标度特性分别提出了基于IM客户端的监控方案和基于服务器端的病毒监控方案。同时研究了个体差异对均匀小世界网络中病毒传播的影响。经过理论分析和大规模的数值仿真发现:只有满足一定的临界条件后,病毒才能在均匀小世界网络中流传。进一步,复杂网络中即使存在有感染能力非常强的个体,只要这些个体的比例小于一定的临界值,病毒也无法在均匀小世界中传播开。此外,我们提出针对复杂网络上的病毒传播实行局域控制策略,在图论距离意义下研究其局域控制范围,对比完全随机网络、小世界网络和无标度网络这三种典型结构,指出了局域控制的有效性及适用范围。我们发现对于小世界网络,存在零感染的局域控制范围为3。在欧氏距离空间下,我们建立了一个欧氏空间的小世界网络模型,研究其病毒传播临界行为,并提出研究欧氏空间意义下的免疫控制范围的问题,指出了临界免疫半径的存在性,为现实中如何控制区域传播范围提供参考。针对一维小世界网络,我们提出了三种线性SIR模型,研究了在不同类型延迟响应下的病毒传播动力学行为特性,并具体分析了延迟与震荡传播行为之间的依赖关系。国内外虽然开展了不少关于网络上的传播的研究工作[346],但是理论研究仍然与实际有一定的差距,今后有待进一步深入研究。

关于复杂网络的控制问题。在第9章讨论过混沌复杂网络的牵制控制方法,特别是,如何控制最少的节点来实现整个网络中分区(不同社区、群落、子网络等)中的多目标控制,这是需要继续深入研究的兼具理论和实际意义的一个重要课题。为此,我组还提出了另外一种控制策略:首先对每个节点增加一个常数输入,这个常数依赖于这个节点和它的邻居期望达到的状态;然后随机从每个子网络中选取少量的节点,并给它们分别加上一个线性误差反馈。我们导出了网络局部或者全局稳定到不同(异质)平衡点/周期态的条件,模拟实验结果验证了控制策略的有效性。同时,针对一类有向动态网络网络模型,我组提出了一种新的基于ControlRank(CR)的牵制控制方式,得到了比单纯基于度的牵制控制策略更好的控制效果。这是基于Google使用的衡量网页重要性的PageRank(PR)值而提出复杂网络中的节点的CR值的概念的。CR与PR之间存在对偶关系:如果把一个原始网络中所有边的方向都反向,则新网络中每个节点的PR值就等于原始网络中相应节点的CR值。我们通过在一类有向无标度网络中的仿真实验验证了基于CR的牵制策略的优势。另外我们在第10章中还把牵制控制应用于带有虚拟领导者的蜂拥算法中,研究了群体中只有一部分的智能体能够感应到虚拟领导者的信息的控制问题,构造了一个新的蜂拥控制算法,得到多智能体的质心是以指数收敛虚拟领导

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

433

者的运动位置和速度,且给出了具体的指数表达式。同时提出了多虚拟领导者的蜂拥控制算法,得到了有意义的结果。该课题的实际应用仍然任重而道远。

关于混沌网络的入侵检测问题[350]。针对Internet网的复杂性和网络的入侵检测重要性,我们不仅研究了复杂网络的复杂性问题,而且探索了网络上入侵检测问题。在传统的客观距离的基础上,我们引入了主观距离的概念,从而提高了判断入侵事件的准确性。研究发现:入侵检测的遗传算法可以归结为一个复杂性网络问题,而且是一个NP困难问题。为了便于进行入侵检测,我们提出了基于SQL(StructuralQueryLanguage)语言的TTT语言(Time2To2pology2Topic),其主要目的是从时间、逻辑拓扑、关键标题三个方面着手,分析复杂网络数据库中的数据。TTT语言类似SQL语言,但具有在线分析和数据挖掘的功能,作为分析复杂网络入侵检测的分析用。关于网络的入侵检测问题尚需深入研究。

关于社会经济网络的问题[351～356]。这是复杂网络很重要的一个研究方向,研究内容十分丰富,问题相当复杂。例如,创业型经济网络,是一个以创业企业为主体,政府规制和激励为引导、投资为支撑、中介服务为纽带、大学和科研机构为创业源的“生态网络”。中小创业网络的各结点间存在着相互制约、相互依赖、相互促进的关系。如果某一节点演化/发育不全,发展滞后,就会影响整个生态网络的发展和功能的发挥,从而影响创业网络的整体活力和效率。如何运用复杂网络理论及其分析方法,试图打开创业生态网络的结构与网络功能相互作用机制的黑箱,主要关心的问题是:创业网络的总体结构属性如何具体刻画?创业网络结构属性与创业网络功能有什么内在联系?建立在网络结构基础上的创业网络功能具有哪些特征?如何通过创业网络的结构优化实现创业型经济的持续健康发展?等等。为此,我们试图应用第5章中

[49～52]

提出的“和谐统一的混合模型”以及具有变速增长的统一混合模型[369,377]来探讨高科技产业网及创业经济网络,把理论研究与实证研究结合起来,进一步改进和完善网络理论模型,使之更适用于高新科技企业网络和创业型经济网络,及其它社会经济网络的的发展需要,试图从微观和宏观相结合的多层面来揭示高新科技发展的机制、特征和规律,这方面研究已取得了一定进展,这些结果[369～375]包括:中国硅谷———中关村科技园技术产业网、全国高新园区网、全国高新技术产业网、全国电子信息技术百强企业网和世界五百强企业网,全国高校科技园网,等等。期望能为我国宏观经济调控和高新技术发展,以及提升国家自主创新能力提供支撑和一定的科学决策依据。

面临的挑战性课题。序言里已经指出,网络科学的研究面临着许多迫切需要解决的课题,无论是应用方面,还是学术方面,都极具挑战性。为了进一步深入开展网络科学及其应用研究,我们认为,还有以下九方面的学术与应用课题是值得关注和重视的[357～369,376,378]。第一,继续深入探索网络科学的理论模型。正如美国有影响的科学家E.O.Wilson指出[357]:“今天最大的挑战性,不仅是细胞生物学和生态学,而是科学的所有方面,特别是如何精确地和完全地描述复杂系统。科学家已经认识了许多类型的复杂系统。他们认为已经知道系统中大多数元素和受力情况。下一步的任务就是怎么综合起来,至少在数学模型方面必须抓住整个系综的关键性质。”因此,首先,必须从理论上继续深入探索复杂动态网络的数理模型,建立全面符合实际网络的精确理论框架.对于我们,将继续探索大统一混合网络模型,量子信息和纳米网络等微观模型,广义(部分)同步模型,各种权重演化网络、下一代互联网、社区网络和高科技企业网络等模型,都需要进一步进行精确的理论分析,以准确地反映实际网络的

434物理学进展　　27卷

全面特性。

第二,开拓复杂网络的各种拓扑构造方法。鉴于复杂网络的多样性和复杂性,需要深入开展从不同角度研究不同复杂网络的生长方式,不仅有随机方法和确定性方法,而且非常需要统一的多层次的混合方法。深入揭示和科学理解不同拓扑特性与动力学行为(功能)之间关系及互相转变规律,进一步探索复杂网络从特殊到一般的深层次的统计规律、非统计和动力学新特性,一直是物理学家对网络科学研究所追求目标之一。

第三,加强变速增长的有权网络的研究工作。这是一个需要深入开展研究的一个重要方向。在许多真实世界的形形色色网络中,不仅几乎都是有权的(包含各种相互作用),而且网络演化和增长的速度各不相同,如在因特网、万维网、科学家合作网及新陈代谢网等现实网络中,边的增长比节点增长速度要快,人们称这种现象为加速增长,实际上网络演化包含有正常、减速、加速和超加速等变速的多种复杂情况。目前,这方面变速的增长网络的演化机制与模型研究尚少。因此,需要对加速和超加速增长现象进行深入研究,分析不同变速网络形成的机理,在此基础建立符合各类型的变速有权网络的演化模型,使其适合于现实世界中网络的统计和动力学特性,并进一步分析超、加、减速增长对网络结构特性和动力学特性的影响,这将是一项很有意义课题,值得引起重视。第四,深入研究复杂网络结构与动力学及其各种同步的关系,一直是网络科学中人们普遍关注的一个焦点和热点问题。在网络动力学或网络功能的研究中非平凡结构所导致最令人关注是涌现现象,最典型的一个涌现问题就是同步,包括完全同步和广义/部分同步,这复杂网络上的典型的集体运动形式和网络结构导致的涌现现象。目前研究显示,复杂网络各个结构特征量都与之有着复杂的内在关系,目前许多工作所依据的网络结构特征还是局域的和静态的。研究进展表明,从集体运动角度考察复杂网络结构和动力学过程是非常重要的,同时复杂网络中还普遍存在着分形或自相似结构也与此相关,这些都是还没有解决的而值得深入研究的前沿课题。因此,继续探索和定量描述复杂网络全局的结构/拓扑特征和动力学行为以及功能之间的内在关系和规律,无疑是推进复杂网络动力学理论和应用的关键之一。

第五,深入研究非线性动态加权复杂网络中的动力学时空复杂性和协调控制方法。结合小世界和无标度等不同类型网络,进一步进行理论分析,深入搞清动态网络演化过程产生的分叉、混沌、阵发混沌和班图涌现的演化表现形式、特点和规律。从统一的高度来探讨动态网络上动力学或功能与拓扑结构的关系。特别关注象生物集群这样动态网络中动力学行为的协调控制策略和方法的研究,并应用于研究多智能体机器人问题。

第六,深入探索不同领域及其不同类型网络的差异、共性及其内在联系,特别是,社会网络与技术网络和生物网络等拓扑特性和动力学行为的特殊性及其关联,进一步从宏观到微观全面揭示不同类型网络的非线性动态演化的规律和物理机制。为此,随着研究的不断深入,需要提出能够全面深刻揭示更本质的复杂动态网络的新理论、新方法和新特征量,后者不仅需要与网络拓扑结构有关的几何量,而且需要反映网络动态变化特点和全局性质的物理量,还需要与动态有关的,如网络通信、传输和传播等与信息有关的特征量。第七,进一步发展网络科学的多种交叉理论和新方法,推进整体复杂动态网络研究向综深入发展,为全面深入研究和复杂系统发展提供广泛的理论基础。

第八,关注生物复杂网络研究,促进系统生物学发展。国内外正在掀起的系统生物研究的

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

435

热潮,国内有关单位,包括我们合作单位也已经开始这方面的探索工作,该课题具有深远意义,但是任务艰巨,相信只要下大功夫,不断努力,有望今后取得有意义的进展。

第九,加强与复杂动态网络的各种应用相关的课题研究,如前所述,提出了许多应用研究课题。其中,如怎么确保复杂网络可靠运转以有效应对网络灾变?如何控制计算机病毒和各种传染病传播(包括:艾滋病、非典和禽流感等)以应对给人类造成巨大的威胁?这些具有实际应用需要而紧切的课题,一直是网络科学应用研究的重要任务之一。只有结合大量的实际应用问题深入研究,从特殊到一般,才能揭开复杂网络的所有奥秘,才有希望找到解决一系列难题的正确答案和有效的途径。

人们虽然面对网络科学如此空前严峻的挑战,但是由于她那诱人而广阔的应用发展前景更加坚定了人们勇攀科学高峰的意志。迄今,网络科学发现的若干重大规律极大地改变和丰富了人们对复杂世界的认识,揭示和开拓了前所未有的理论和技术问题,因而极大地激起人们对网络科学的研究浓厚的兴趣和探索的激情,国内外掀起了网络科学研究和应用的热潮,一浪高过一浪,方兴未艾。各类论文如雨后春笋般地大量地发表,国内的研究进展同样喜人。网络科学诞生的时间虽短,已经从多角度、多领域对复杂网络展开了空前的研究,无论是理论模型探索,还是真实世界网络的实证研究,网络科学的成果必将大大有助于人类解决面临的一些重要问题,促进整个社会的物质文明和精神文明的建设。

展望未来,网络科学的研究首先对于整个国家安全具有现实的和长远的战略意义。21世纪将是互联网和网络信息的时代。互联网的发展已经和继续强劲地带动计算机、微电子、通信和软件等信息产业的高度发展,这成为21世纪全球军事和经济的主要推动力,正是这些迫切需求不断地推动因特网向宽带、高速方向发展,特别是正在向光互联网和未来的量子互联网方向发展。所谓“量子互联网”是科学家们正在探索中的一个崭新量子网络,通过这个新奇的量子网络传送宇宙间最奇特的物质,这种新物质就是所谓“缠结”信息。“缠结”信息可用于量子密码翻译、极小规模的量子计算和远距传送信息。科学家认为,建立一个产生、储存和传送“缠结”信息的“量子互联网”,是向开发远距离传物系统迈出的重大的关键的第一步,并为“量子互联网”诞生铺平道路。因为利用这种“缠结”信息能够制造超快速量子计算机,把它们连接成量子互联网,这为目前传统的互联网和正在发展中的光因特网开辟一条崭新的途径,只有这样,国家才能适应未来对于安全问题的挑战和对经济持续高度发展的需求,下一代“量子互联网”的应用前景将不可估量。

令人关注国内外正在大力研发中的特设(adhoc)通信网,已经成为国内外网络科学研究的一大热点。所谓“特设通讯网”的主要特点是,不仅具有多跳、无中心、自组织、无基站网络的特点和功能,而且又具有结点快速移动、拓扑结构频繁变化和连接短暂的特性。特设通讯网不仅大为军方所青睐,因为它可以为少数高级指挥人员所专用,在现代信息化战争中具有特殊的实战意义,而且可以为重要的政府部门和商务机密使用,成为未来最快捷、最高效的、最先进的一种通信手段。

值得强调指出,军事上对网络科学有迫切的需求和应用前景。在1991年海湾战后,美军发现在网络中心作战实践中由于对伊军因特网具有先进的动态路由选择技术一直并不知情,因而迟迟没有完全切断伊军指挥控制网络,直到最后伊军还保留一条主要干线的光纤电缆,这就是导致美军当时相当艰难才打垮了伊拉克的军事指挥控制网络的主要原因,充分暴露了美

436物理学进展　　27卷

军网络中心作战存在的一个关键问题。将美军打击伊军因特网战例与具有无标度特性的因特网的攻击实验进行分析对比便知,当时如果美军能有组织地协同攻击伊军网络中心集散节点,就必然能够加速战争的胜利。鉴于美国两次海湾战争的教训,从1997年4月美国军队开始拨巨汇,加紧推进所谓“网络中心作战研究”,从而对网络科学提出了更高的需求和投入研究。如2004年9月美国科学院国家研究委员会所属“陆军科学技术专业委员会”开始了“网络科学在未来陆军的应用”等多个与军事密切的相关项目研究,不难看出网络科学为发达国家军方高度注视的程度了。这种态势也应该引起我国有关部门的重视。

网络科学的研究兼备基础性、前瞻性、指导性和实用性,网络理论已经成为新原理和新技术的源头,是网络工程的设计、安全防护和开发应用的强有力的工具。人们从“爱虫”“、熊猫烧

(台湾地震香”等病毒在互联网上大肆传播、震惊世界的“北美大停电”和“台湾海底电缆事故”演变成史无前例的亚太区通讯网络大灾难)等一系列网络事件中吸取教训,更加关注的问题是:计算机病毒如何在互联网上和万维网上传播而导致流行?人们应该如何阻止/控制病毒在复杂网络上传播蔓延?如何有效地防止黑客侵入?怎样来设计出具有能够有效抵抗意外故障和攻击能力的复杂网络以防止网络上的一系列级联效应?怎样消除不断恶化的全球气候变化和生态环境网络,而保持气候和生态环境良性平衡?等等,这一系列棘手问题无不与与社会生活息息相关,涉及到空间大气网络、因特网、万维网、各种交通运输网、电力网、各种通信网络、卫星电视网、电子邮件网、生态环境网络和食物链网等复杂网络。因此,网络科学的深入研究和发展,必将揭开整个网络的“庐山真面目”,具有特别重要的现实的和长远的意义。

网络科学对于社会领域的的应用更具紧迫性和实际意义,特别是,近年来对人类造成巨大的威胁遍及全球的艾滋病、非典、禽流感等各种传染病等,都已成为世界各国当前最关注的焦点之一。主要问题是:在特定的社会网络中这些传染病是如何通过接触关系传播而导致流行?决策者如何控制和减少这些疾病对人类的伤害把损失降到最低限度呢?深入探索复杂网络的发展规律,可望寻找到解决这些问题的合理的答案和有效途径。类似地,人们已经应用复杂网络理论来模拟社会上谣言的传播过程,寻找控制流言飞语的扩散和降低其负面影响。

网络科学同样在政治、经济和管理等众多领域中能够发挥意想不到的作用,大有用武之地。例如,发生在美国的9・11恐怖事件,令人十分痛恨!当今恐怖主义已成为21世纪的十大危害之一,恐怖网络遍布全球各地。国际上各国家联合起来如何有效地破坏和摧毁恐怖网络和犯罪组织网络?这已成为网络科学研究的当务之急。从无标度网络的特性可知,只要通过捉拿逮捕恐怖网络中节点度最大的少数几个中心首领和主要人物,就可摧毁整个网络系统,使其网络恐怖的功能失常或崩溃,以维护人类社会政治的稳定。网络科学已用于揭示政府议员的组织、层次关系,甚至可以成功地模拟政治选举,预先分析一些因素对选举结果的影响,以避免出现突变情形和导致社会动乱。

综上可见,网络科学已经引起了国内外不同学科的高度重视和密切关注,她的研究任务任重而道远,面临着极具挑战性的理论和应用研究的双重艰巨任务。显然,这些正是网络科学深入开展研究的迫切需求和继续发展的强大推动力,网络科学已成为国内格外最广泛交叉的新兴科学,正在以雷霆万钧之势向前发展,前景十分美好。

为了实现网络科学的研究任务和目标,需要结成世界上最广泛的科学家—工程师合作网,这应该是世界上最强的、最庞大的、最广泛的、加速的、有权的网络。在这样空前的科学合作网

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

437

中,更需要坚持不懈的创新拼搏精神和严谨的科学作风。我们非常愿意与国内同行一道,竭尽全力,为推动我国网络科学在科研、国防和国民经济中的应用发展,发挥我们的积极作用,作出应有的贡献。

让我们为网络科学新时代的到来而热情欢呼吧!

我们衷心地祝愿,并寄重望于:中国科学事业的有志者,特别是我国年轻的一代骄子们,“海阔任鱼跳,天空任飞鸟”,尽情地发挥自己的聪明才智,努力为国际科学界这座最高、最雄伟、最庞大的网络科学的大厦添砖加瓦吧,贡献我国宝贵的一分力量!去铸造21世纪的科学辉煌!

结束语:注记和致谢

完成此文之际,我们不禁想起撰写这个长篇综述文的来龙去脉。2006年本文第一作者在(参看文献[362])上发表了一篇评论文““科技导报”:网络科学的理论模型探索及其进展”,该文

简要评述了复杂网络的分类和进展概况后,着重概述了我们正在进行的国家自然科学基金重点项目近年来在网络科学的理论模型方面的探索及取得的与国际上的同步进展概况。完全出乎意料之外,文章刚刚出来立即引起了《物理学进展》主编冯端院士先生的关注。正值2007年元旦来临之际,第一作者荣幸地收到冯端院士先生的亲笔约稿函,信里提出建议“:在该文的框架下进一步充实具体的科学内容,将可作为“物理进展”的一篇高质量的综述论文。特此向您征稿。望能得到您的支持和响应。”在此,我再次衷心地感谢冯老先生提出了一个很好的建议,给了我们一个进一步整理、分析和系统总结研究成果的难得机会,对我们整个联合项目组起到了莫大的激励和促进作用。特别是,令我们十分惊喜和万分感动是,冯先生作为我国物理学界的泰斗,仍然在80

(中英文版)中将阐明网络科学的前沿课题有关内多高年龄,不仅正在撰写专著《凝聚态物理学》容,而且如此敏锐地抓住当前今世界上新兴科学,并关心国内网络科学的发展。我们欣然同意和

接受了约写这篇综述文的艰巨任务,并且按照冯老信函上指明的方向“:从无规(随机)网络、小世界网络到无标度网络的基本规律谈起,以便初学者有门径可循。”本文撰写正是遵循这个基本宗旨。为此,第一作者立即主持组织了“一院两校”国家自然科学基金重点项目的三个合作单位的主要成员,共同合作来完成这篇综述文。在冯老先生的关心与指引下经过5个多月的共同努力,该文终于完成了。由于题材广泛,内容丰富,篇幅很长,征得冯先生和编辑部的同意,我们分成上下两篇。全文既反映了冯先生的主要意图,具有基础性和导引性,又反映了近年网络科学在国内外的主要进展和重要成果。同时,本文体现了最大的一个特色就是,我们用十章篇幅,以十大专题重点突出地反映和系统总结了我们在该领域在国家自然科学基金重点项目资助(70431002)下,所取得的若干重要进展和主要成果,并提出和评论了网络科学发展所面临的若干挑战性课题,展望了网络科学美好而巨大的发展应用前景。

全文撰写工作,由第一作者负责拟定整个写作提纲目录,并负责撰写了序言、第1、2、3、4(部分)、5、9、12章和结束语,最后对全文负责统一修饰工作。北京师范大学郑志刚和狄增如负责撰写第4、8章;上海交通大学汪小帆和李翔负责撰写第7、10、11章;毕桥负责撰写第6章和第2章部分内容。除了本文署名作者外,中国科学技术大学的博士后杨会杰和复旦大学的博士后章忠志分别提供了网络谱结构和确定性网络的研究资料,在此一并致谢。另外,参与有关资料整理和图文协助的还有我们项目组内三校的博士生李永、卢新彪、荣智海、冯晓琴、王

438物理学进展　　27卷

延、马晓娟、李梦辉、苏厚胜,杨文、李维以及硕士生刘强等。可见,该长篇综述文是我们联合项目组全体大力合作的产物,显示了我们团队真诚合作、开拓创新、共同奋进的精神。

我们向国家自然科学基金委员会对我们重点项目(批准号:70431002)的资助和支持,对以李正道教授为主任的中国高等科技中心多次资助和支持我们主办的“全国复杂动态网络论坛”,以及向以不同方式关心和支持我们课题的国内同行们,特别是陈关荣、汪秉宏、李德毅、邢修三、何大韧和刘增荣等诸教授,还有一些年轻的博士们,一并表示谢意。

由于网络科学是一门广泛交叉的新兴科学,课题内容浩瀚,文章很长,不当之处难免,敬请国内外同行不吝指正,我们将不尽感谢。

最后,我们对《物理学进展》主编冯端院士先生为我们提供这么难得的与国内外同行进行交流的机会,以及对我们课题的关注、鼓励和推动,再此表示真诚的感谢,同时感谢刘协祯副主编的辛勤编辑。我们将以此为动力,在今后网络科学等研究中更上一层楼。参考文献[1]　WattsDJ,StrogatzSH.Nature,1998,393:4402442.[2]　WattsDJ.SixDegrees:TheScienceofaConnectedAge,NewYork:W.W.Norton&Company,

2003.

[3]　BarabásiAL,AlbertR.Science,1999,286:5092512.

[4]　BarabásiAL.ThenewScienceofNetworks.Cambridge:Prerseus,2002.[5]　网络科学:参看http://www.newsinfo.nd.edu/

[6]　WattsDJ.Annualreviewofsociology,2004,30:2432270.

[7]　欧拉(LeonhardEuler,170721783)简况:http://www2.zzu.edu.cn/math/classes/2003/y1/oula.htm[8]　布鲁斯・谢克特著,王元、李文林译《.我的大脑敞开了———天才数学家保罗・爱多士传奇》.上海:译文

出版社,2002.

[9]　MilgramS.PsychologyToday,1967,2:60267.

[10]　DoddsPS,MuhamadR,WattsDJ.2003,301(5634):8272829.[11]　StragotzSH.Nature,2001,410(8):2682276.

[12]　Buchanan,NexusM.SmallWorldsandtheGroundbreakingScienceofNetworks.NewYork:W.W.

NortonandCompany,2002.

[13]　DorogovtsevSN,MendesJFF.AdvancesinPhysics,2002,51:107921181.

[14]　WattsDJ.Smallworlds:Thedynamicsofnetworksbetweenorderandrandomness.NewJersey:

PrincetonUniversityPress,1999.

[15]　BoccalettiS,LatoraV,MorenoY,etal.PhysicsReport,2006,424:1752308.

[16]　DorogovtsevSN,MendesJFF.EvolutionofNetworks:FromBiologicalMetstotheInternetand

WWW.Oxford:OxfordUniversityPress,2003.

[17]　Paster2SatorrrasR.AVespignaniEvolutionandStrictureoftheInternet:AStatisticalPhysicsAp2

proach,Cambridge:CambridgeUniversityPress,2004.

[18]　Editors:EliBen2Naim,HansFrauenfelder,ZoltanTorozckai.ComplexNetworks,LectureNotesin

Physicsvolume650,Amazon:Springer,2004.

[19]　AlbertR,BarabasiAL.Rev.Mod.Phys,2002,74:47298.[20]　NewmanMEJ.SIAMReview,2003,45:1672256.

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

439

[21]　汪小帆,李翔,陈关荣编著.复杂网络理论及其应用,北京:清华出版社,2006.[22]　郭雷,许晓铭主编.复杂网络.上海:上海科技教育出版社,2006.[23]　HolmeP,ParkSM,KimBJ.arXiv,2004cond2mat/0411634.

[24]　MacdonaldPJ,AlmaasE,BarabásiAL.arXiv,2004,cond2mat/0405688.[25]　TieriP,ValensinS,LatoraV,etal.Bioinformatics,2005,21(8):163921643.[26]　BarabásiAL,JeongH,NédaZ,etal.PhysicaA,2002,311:5902614.

μaR,MossaS,TurtschiA,etal.Proc.Nat.Acad.Sci.USA2005,102:779427799.[27]　Guimer

μaR,AmaralLAN.Eur.Phys.J.B2004,38:3812385.Guimer

[28]　AlmaasE,KovacsB,VicsekT,etal.Nature,2004,427:8392843.

[29]　BarratA,BarthélemyM,Pastor2SatorrasR,etal.Proc.Natl.Acad.Sci.U.S.A.2004,101(11):

374723752.

[30]　OnnlaJP,SaramAakiJJ,Kertesz,etal.Phys.Rev.E,2005,71:065103.[31]　NewmanMEJ.Phys.Rev.E,2003,67:026126.[32]　NewmanMEJ.Phys.Rev.Lett.,2002,89(20):208701.[33]　PastorSR,VázquezA,VespignaniA.Phys.Rev.Lett.2001,87:258701.[34]　RavaszE,BarabásiAL.Phys.Rev.E,2003,67:026112.[35]　FreemanL.Sociometry,1977,40:35241.[36]　BarthélemyM.TheEuropeanPhysicalJournalB,2004,38:1632168.[37]　MiloR,ItzkovitzS,KashtanN,etal.2004,303(5663):153821542[38]　GagenGM.MattickJS.Phys.Rev.E,2005,72:16123.[39]　SenP.Phys.Rev.E,2004,69:46107.

[40]　DavidMDS,JukkaPO,NeilFJ.arXiv,2007,physics/0701339.

[41]　GirvanM,NewmanMEJ.Proc.Natl.Acad.Sci.USA,2002,99(12):782127826.[42]　NewmanMEJ.TheEuropeanPhysicalJournalB,2004,38:3212330.[43]　NewmanMEJ,GirvanM.Phys.Rev.E,2004,69:026113.[44]　NewmanMEJ.Phys.Rev.E,2004,69:066133.

[45]　Erd.sP,RényiA.Publ.Math.Debrecen,1959,6:2902297.[46]　BianconiG,BarabásiAL.Europhys.Lett.,2001,54:4362442.[47]　JinEM,GirvanM,NewmanMEJ.Phys.Rev.E,2001,64:046132.[48]　HolmeP,KimBJ.Phys.Rev.E,2002,65:026107.[49]　方锦清,毕桥,李永,等.中国科学G辑,2007,3(2):2302249.[50]　FangJQ,LiangY.ChiPhysLett,2005,22(10):271922722.[51]　LuXB,WangXF,FangJQ.PhysicaA,2006,371(2):8412850.[52]　FangJQ,BiQ,LiY,etal.Chin.Phys.Lett.,2007,24(1):2792283.[53]　DorogovtsevSN,MendesJFF.Europhys.Lett.2000b,52:33239.

[54]　YookSH,JeongH,BarabásiAL,etal.Phys.Rev.Lett.,2001,86:583525838.[55]　ZhengD,TrimperS,ZhengB,etal.Phys.Rev.E,2003,67:040102.[56]　BarratA,BarthélemyM,VespignaniA.Phys.Rev.Lett.,2004,92:228701.[57]　AntalT,KrapivskyPL.Phys.Rev.E,2005,71:026103.

[58]　WangWX,WangBH,HuB,etal.Phys.Rev.Lett.,2005,94:188702.[59]　ParkJ,NewmanMEJ.Phys.Rev.E,2004,70:066117.

[60]　PrigogineI.Structure,DissipationandLife,TheoreticalPhysicsandBiology,Versailles,1967,Am2

440物理学进展　　27卷

sterdam:North2HollandPubl.Company,1969.(Itisinthiscommunicationthattheterm\"structuredissipative\"isusedforthefirsttime.)

[61]　JeansJH.Thedynamicaltheoryofgases,Inc.NewYork:DoverPublications,1954.[62]　XingXS.(inChinese)arXiv,2004,cond2mat/0505009.

[63]　BianconiG,BarabasiAL.Phys.Rev.Lett.2001,86:563225635.[64]　PecoraLM.Phys.Rev.E,1998,58:3472360.

[65]　BoccalettiS,PecoraLM,PelaezA.Phys.Rev.E,2001,63:066219.[66]　PecoraLM,CarrollTL.Phys.Rev.Lett.,1998,80:210922112.[67]　JohnS,SompolinskyH,StephenMJ.Phys.Rev.B,1983,27:559225603.[68]　SongC,HavlinS,MakseHA.Nature,2005,43:3922395.[69]　SongC,HavlinS,MakseHA.NaturePhysics2006,2:2752281.

[70]　FarkasI,DerényiJI,BarabásiAL,etal.Phys.Rev.E,2001,64:026704.[71]　YangH,ZhaoF,QiL,etal.Phys.Rev.E,2004,69:066104.[72]　ZhaoF,YangH,WangB.Phys.Rev.E,2005,72:046119.[73]　YangH,ZhaoF,WangB.PhysicaA,2006,364:5442556.[74]　YangH,ZhaoF,WangB.CHAOS,2006,16:043112.[75]　BuldyrevSV,GoldbergerAL,HavlinS,etal.Phys.Rev.E,1995,51:508425091.[76]　Bo.ekP,PoszajczakM,BotetR.Phys.Rep.1995,252:101.

[77]　YangH,ZhaoF,ZhuoY,etal.Phys.Lett.A,2002,292:349.YangH,ZhaoF,ZhuoY,etal.

PhysicaA,2002,312:23.

[78]　AndersonPW.Phys.Rev.1958,109:1492.

[79]　JohnS,SompolinskyH,StephenMJ.Phys.Rev.B,1983,28:6358.[80]　JohnS,SompolinskyH,StephenMJ.Phys.Rev.Lett.,1984,53:2169.

[81]　Erd.sP,RényiA.Ontheevolutionofrandomgraphs,PublicationoftheMathematicalInstituteofthe

HungarianAcademyofScience,1960,5(1):17261.参看:郭雷,许晓铭主编.复杂网络.2842328,上

海:上海科技教育出版社,2006.

[82]　BollobasB.RandomGraph,Secondedition,Cambridge:CambridgeUniversityPress,2001.[83]　DurrettR.RandomGraphDynamics,Cambridge:CambridgeUniversityPress,2006.[84]　BenderEA,CanfieldER,CombJ.TheoryA,1978,24:296.[85]　MolloyM,ReedB.RandomStruct.Algorithm,1995,6:161.[86]　MolloyM,ReedB.Conbin.Probab.Comput.1998,7:295.[87]　NewmanMEJ,WattsDJ.Phys.LettA,1999,263:3412346.[88]　NewmanMEJ,WattsDJ.Phys.Rev.E,1999,60:733227342.[89]　KasturiranganR.arXiv,1999,cond-mat/9904055.

[90]　Dorogovtsev,SN,MendesJFF.Euro.Lett.,2000,50:127.[91]　KleinbergJ.Nature,2000,406:8452845.

[92]　LiY,FangJQ,LiuQ,etal.CommuTheorPhys.2006,45(5):9502954.[93]　刘强,方锦清,李永,等.复杂系统和复杂性科学,2005,2(2):13219.[94]　FronczakA,Ho.ystJA,JedynakM,etal.PhysicaA,2002,316:6882694.[95]　FronczakA,FronczakP,HolystJA.Phys.Rev.E,2003,68(4):046126.[96]　KrapivskyPL,RednerS,LeyvrazF.Phys.Rev.Lett.,2000,85:462924632.[97]　DorogovtsevSN,MendesJFF.Phys.Rev.E,2000,62:184221845.

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

441

[98]　AmaralLAN,ScalaA,BarthélémyM,etal.Proc.Natl.Acad.Sci.USA,2000,97(21):111492

11152.

[99]　DorogovtsevSN,MendesJFF,SamukhinAN.Phys.Rev.E,2001,63(6):062101.

[100]　KleinbergJM,KumarR,RaghavanP,etal.ProceedingsoftheInternationalConferenceonCombi2

natoricsandComputing.LectureNotesinComputerScience,1999,1627:1217.

[101]　KumarR,RaghavanP,RajalopaganS,etal.Proceedingsofthe42stAnnualIEEESymposiumon

FoundationsofComputerScience,InstituteofElectronicandElectronicsEngineers,NewYork,2000:57265.

[102]　ChungF,LuLY,DeweyTG,etal.J.ofComput.Biology,2003,10(5):6772688.[103]　KrapivskyLP,RednerS.Phys.Rev.E,2001,63:066123.[104]　VázquezA.Europhys.Lett.2001,54:4302435.

[105]　LiuJG,DangYZ,WangZT.Chin.Phys.Lett.,2006,23(3):462749.[106]　AlbertR,BarabásiAL.Phys.Rev.Lett.2000,85:523425237.[107]　ChenQH,ShiDH.PhysicaA,2004,335:2402248.[108]　ErgünG,RodgersGJ.PhysicaA,2002,303:2612272.[109]　SzabóG,AlavaM,KertészJ.Phys.Rev.E,2003,67:056102.[110]　KlemmK,EguíluzVM.Phys.Rev.E,2002,65:036123.[111]　KlemmK,EguíluzVM.Phys.Rev.E,2002,65:057102.[112]　VázquezA.Phys.Rev.E,2003,67:056104.[113]　Saram.kiJ,KaskiK.PhysicaA,2004,341:80286.

[114]　ZhangZZ,RongLL,ComellasF.JournalofPhysicsA:MathematicalandGeneral,2006,39(13):

325323261.

[115]　RavaszE,SomeraAL,MongruDA,etal.Science,2002,297:155121555.[116]　WangSJ,ZhangCS.Phys.Rev.E,2004,70:066127.[117]　LiX,JinYY,ChenGR.PhysicaA,2003,328:2872296.[118]　LiX,ChenGR.PhysicaA,2003,328:2742286.

[119]　WangB,TangHW,ZhangZZ,etal.InternationalJournalofModernPhysicsB,2005,19(26):395123959.[120]　BarabásiAL,RavaszE,VicsekT.PhysicaA,2001,299:5592564.[121]　ComellasF,OzónJ,PetersJG.Inf.Process.Lett.2000,76:83290.[122]　ComellasF,SampelsM.PhysicaA,2002,309:2312235.[123]　XiaoWJ,ParhamiB.Inf.Process.Lett.2006,97:1152117.[124]　CorsoG.Phys.Rev.E,2004,69:036106.[125]　AchterJ.Phys.Rev.E,2004,70:058103.

[126]　ChandraAK,DasguptaS.PhysicaA,2005,357:4362446.[127]　ZhangZZ,RongLL,GuoCH.PhysicaA,2006,363(2):5672572.[128]　IguchiK,YamadaH.Phys.Rev.E,2005,71:036144.

[129]　DorogovtsevSN,GoltsevAV,MendesJFF.Phys.Rev.E,2002,65:066122.[130]　ComellasF,FertinG,RaspaudA.Phys.Rev.E,2004,69:037104.[131]　ZhangZZ,RongLL,ZhouSG.PhysicaA,2007,377:3292339.[132]　JungS,KimS,KahngB.Phys.Rev.E,2002,65:056101.[133]　NohJD.Phys.Rev.E,2003,67:045103.

[134]　NacherJC,UedaN,KanehisaM,etal.Phys.Rev.E,2005,71:036132.

442物理学进展　　27卷

[135]　ZhouTao,WangBinghong,HuiPM,etal.PhysicaA,2006,367:6132618.

[136]　AndradeJS,HerrmannHJ,AndradeRFS,etal.Phys.Rev.Lett.,2005,94:018702.[137]　DoyeJPK,MassenCP.Phys.Rev.E,2005,71:016128.[138]　ZhouT,YanG,WangBH.Phys.Rev.E,2005,71:046141.

[139]　ZhangZZ,ComellasF,FertinG,etal.JournalofPhysicsA:MathematicalandGeneral,2006,39

(8):181121818.

[140]　ZhangZZ,RongLL,ComellasF.PhysicaA,2006,364:6102618.[141]　ZhangZZ,RongLL,ZhouSG.Phys.Rev.E,2006,74:046105.[142]　HinczewskiM,BerkerAN.Phys.Rev.E,2006,73:066125.[143]　ZhangZZ,ZhouSG,ZouT.arXiv:cond2mat/0612427.[144]　ZhangZZ,ZhouSG,ZouT.(unpublished)

[145]　DorogvtsevSN,MendesJFF.AIPConf.Proc.,2005,776,29236.[146]　LiW,CaiX.Phys.Rev.E,2005,69:046106.[147]　LiCG,ChenGR.PhysicaA.2004,343:2882294.[148]　WangWX,WangBH,ZhouT,etal.Phys.Rev.E,2005,72:046140.[149]　WuZX,XuXJ,WangYH.Phys.Rev.E,2005,71:066124.

[150]　WuZX,XuXJ,WangYH.TheEuropeanPhysicalJournalB,2005,45:3852390.[151]　PanZF,LiX,WangXF.Phys.Rev.E,2006,73:056109.[152]　XieZ,LiX,WangXF.CommunTheorPhys,2006(Accepted).

[153]　方锦清,中国当代科技学术专刊,北京,发现杂志社,2006,122[详见170,171].[154]　LiMH,FanY,ChenJW,etal.PhysicaA,2005,350:6432656.

[155]　FanY,LiMH,ChenJW,etal.InternationalJournalofModernPhysicsB,2004,18(17219):

250522511.

[156]　ZhangP,LiMH,WuJS,etal.PhysicaA,2006,367:5772585.[157]　LiMH,WangDH,FanY,etal.NewJournalofPhysics,2006,8:72.[158]　LiMH,WuJS,WangDH,etal.PhysicaA,2007,375:3552364.[159]　NewmanMEJ.Phys.Rev.E,2004,70:056131.

[160]　ReichardtJ,BornholdtS.Phys.Rev.Lett.,2004,93:218701.[161]　DuchJ,ArenasA.Phys.Rev.E,2005,72:027104.

[162]　FanY,LiMH,ZhangP,etal.PhysicaA,2007,378:5832590.

[163]　LiMH,FanY,WangDH,etal.Phys.Lett.A(Inpress),2006,doi:10.1016/j.physleta.2006.

12.040.

[164]　PallaG,DerenyiI,FarkasI,etal.Nature,2005,453:8142818.[165]　PallaG,DerenyiI,VicsekT.Phys.Rev.Lett.2005,94:160202.[166]　FanY,LiMH,ZhangP,etal.PhysicaA,2007,377:3632372.[167]　NewmanMEJ.Phys.Rev.E,2001,64:016131.

[168]　FangJQ,BiQ.In:ProceedingsofTheForthInternationalConferenceonNonlinearScience,Pohang,

2006,34235.

[169]　WangWX,HuB,WangBH,etal.Phys.Rev.E,2006,73:016133.[170]　方锦清.科技导报,2006,24(12):67272.

[171]　FangJQ,BiQ.LiY.Front.Phys.China,2007,1:1092124.[172]　LiuZH,LaiYC,YeN,etal.Phys.Lett.A,2002,303:3372344.

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

443

[173]　Soon2HyungY,FilippoR,HildegardMO.Phys.Rev.E,2005,72(4):045105.[174]　BrajendraKS,NeelimaG.Phys.Rev.E,2003,68(6):066121.[175]　SumiyoshiA,NorikazuS.Phys.Rev.E,2006,74(2):026113.[176]　GrabowskiA,KosiskiRA.Phys.Rev.E,2006,73(1):016135.[177]　JulianS,JanuszAH.Phys.Rev.E,2005,72(4):046127.

[178]　MicheleC,GuidoC,LucianoP.Phys.Rev.E,2004,70(3):037101.[179]　Xulvi2BrunetR,SokolovIM.Phys.Rev.E,2004,70(6):066102.[180]　LiCG,PhilipKM.Phys.Rev.E,2005,72(4):045102.

[181]　Xulvi2BrunetR,PietschW,SokolovIM.Phys.Rev.E,2003,68(3):036119.

[182]　BreuerHP.Thetheoryofquantumopensystems,NewYork:OxfordUniversityPress,2002.[183]　HedinER,CosbyRM,JoeYS,etal.J.Appl.Phys.,2005,97:063712.[184]　SchwingerJS.Particles,SourcesandFields,MA:PerseusBooks,1998.[185]　VanKampenNG.Physica,1974,74:214.

[186]　BiQ,FangJQ,ZouQ.Chi.Phys.Lett.2006,23:1947.[187]　BiQ,RudaHE,ZhouDZ.PhysicaA,2006,364:170.[188]　HasegawaH,MurankaT,KasaiS,etal.Jpn.J.Appl.Phys.2003,42:2375.[189]　BiQ,RudaHE,ZhouDZ.PhysicaA,2006,198:210.[190]　AntoniouI,MelnikovY,BiQ.PhysicaA,1997,246:97.

[191]　CatanzaroM,CaldarelliG,PietroneroL.PhysicaA,2004,338:119.[192]　StrogatzSH,StewartI.Sci.Am.1993,269(6):1022109.[193]　NedaZ,RavaszE,VicsekT,etal.Nature,2000,403:8492850.

[194]　NédaZ,RavaszE,VicsekT,etal.Phys.Rev.E,2000,61(6):698726992.[195]　赵明,汪秉宏,将品群等.物理学进展,2005,25(3):273229.[196]　MohantyP.Nature,2005,437(7057):3252326.

[197]　KakaS,PufallMR,RippardWH,etal.Nature,2005,437(7057):3892392.[198]　LuXB,LiX,WangXF.Communicationintheoreticalphysics,2006,45:9552960.[199]　LuXB,WangXF,LiX,FangJQ.PhysicaA,2006,370:3812389.[200]　FanJ,WangXF.PhysicaA,2005,349(3):4432451.[201]　FanJ,LiX,WangXF.PhysicaA,2005,355(5):6572666.[202]　BarahonaM,PecoraLM.Phys.Rev.Lett.,2002,89(5):054101.[203]　NishikawaT,MotterAE.Phys.Rev.E,2006,73(6):065106.[204]　ZhouC,KurthsJ.Phys.Rev.Lett.2006,96(16):164102.

[205]　WangXF,ChenG.IEEETrans.CircuitsandSystemsI,2002,49(1):54262.[206]　WangXF,ChenG.J.SystemsScienceandComplexity,2003,16:1214.

[207]　WangXF,ChenGR.InternationalJournalofBifurcation&Chaos,2002,12:1872192.[207]　NishikawaT,MotterAE.Phys.Rev.Lett.,2003,91:014101.

[208]　HongH,KimBJ,ChoiMY,etal.Phys.Rev.E,2004,69:0671052067109.[209]　MotterAE,ZhouCS,andKurthsJ.Phys.Rev.E,2005,71:016116.[210]　ChavezM,HwangDU,AmannA,etal.Phys.Rev.Lett.,2005,94:218701.[211]　McGrawPN,MenzingerM.Phys.Rev.E,2005,72:01510.[212]　KimBJ.Phys.Rev.E,2004,69:045101.[213]　WuChW.Nonlinearity,2005,18:105721064.

444物理学进展　　27卷

[214]　LuWL,ChenTP.IEEEtransactionsoncircuitsandsystems2II:expressbrief,2007,54(2):1362140.[215]　PecoraLMandCarrollTL.Phys.Rev.Lett.1990,64:8212824

[216]　PikovskyA,RosenblumM,KurthsJ.Synchronization:Auniversalconceptinnonlinearsciences,

Cambridge:CambridgeUniversityPress,2001.

[217]　BoccalettiS,GrebogiC,LaiYC,etal.Phys.Rep.2000,329:103.[218]　BoccalettiS,KurthsJ,OsipovG,etal.Phys.Rep.2002,366:1.[219]　TimmeM,WolfF,GeiselT.Phys.Rev.Lett.,2004,92:074101.[220]　RosenblumM,PikovskyA,KurthsJ.Phys.Rev.Lett.,1996,76:1804.

[221]　ErnstU,PawelzikK,GeiselT.Phys.Rev.Lett.1995,74:157021573.HuG,ZhangY,Cerdeira

HA,etal.Phys.Rev.Lett.,2000,85:337723380.

[222]　ParekhN,ParthasarathyS,SinhaS.Phys.Rev.Lett.,1998,81:140121404.[223]　HuntBR,OttE,YorkeJA.Phys.Rev.E,1997,55(4):4029.[224]　LiuZH,ChenSG.Phys.Rev.E,1997,56:7297.[225]　ZhengZG,WangXG,CrossMC.Phys.Rev.E,2002,65(5):056211.[226]　OsipovG,PikovskyA,RosenblumM,etal.Phys.Rev.E,1997,55(3):2353.[227]　ParlitzU,JungeL,KocarevL.Phys.Rev.E,1996,54:2115.[228]　JalanS,AmritkarRE.Phys.Rev.Lett.,2003,90:014101.

[229]　HwangDU,ChavezM.AmannA,etal.Phys.Rev.Lett.2005,94:138701.[230]　NewmanMEJ,MooreC,WattsDJ.Phys.Rev.Lett.,2000,84:320123204.

[231]　Lago2FernandezLF,HuertaR,CorbachoF,etal.Phys.Rev.Lett.,2000,84:275822761.[232]　WeiGW,ZhanM,LaiCH.Phys.Rev.Lett.,2002,89:284103.[233]　AtayFM,JostJ.Phys.Rev.Lett.,2004,92:144101.

[234]　JiangY,Lozada2CassouM,VinetA.Phys.Rev.E,2003,68:065201.[235]　RestrepoJG,OttE,HuntBR.Phys.Rev.E,2004,69:066215.[236]　MorenoY,PachecoAF.Europhys.Lett.,2004,68:603.[237]　YonkerS,WackerbauerR.Phys.Rev.E,2006,73:026218.[238]　GrayCM,KoenigP,EngelAK,etal.Nature,1998,338:334.

[239]　PyragasK.Phys.lett.A,1992,170:421;Phys.Rev.E,1996,54:4508.

[240]　StewartI,GolubitskyM,PivatoM.SIAMJ.AppliedDynamicalSystems,2003,2(4):609.[241]　WangXF,ChenGR.IEEECircuits&SystemsMagazine,2003,3(1):6220.

[242]　ZhigangZheng,XiaoqinFeng,BinAo,andCrossM.Fronti.Phys,China,2006,1(4):4582467.[243]　BinAoandZhigangZheng.Europhys.Lett.2006,74(2):2292235.

[244]　BinAo,XiaojuanMa,YunyunLi,andZhigangZheng.Chin.Phys.Lett.2006,23(4):7862789.[245]　BinAo,XiaojuanMa,XiulinLi,andZhigangZheng.Int.J.Mod.Phys.2007,21:99521003.[246]　BinAoandZhigangZheng.Phys.Rev.E(2007,submitted).[247]　方锦清.驾驭混沌与发展高新技术,北京:原子能出版社,2002.[248]　方锦清,陈关荣.物理学进展,2003,23(3):3212388.

[249]　FangJQ,WangZS,ChenGR.Commun.Theor.Phys.2004,42(4):5572560.[250]　LiuQ,FangJQ,LiY.Commu.Theor.Phys.(Beijing,China),2007,47(4):7522758[251]　LiuQ,FangJQ,LiY.ProgressinNatureScience,accepted,2006[252]　刘强,方锦清,李永.复杂系统与复杂性科学,2007,4(1):6214.

[253]　FangJQ.TheThirdChina2SingaporeJointSymposiumonResearchFrontiersinPhysics,Invited

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

talk,25227,May2007

445

[254]　ChenTP,LiuXW,LuWL.IEEETransactionsoncircuitsandsystems2I:Regularpaper,2007,54

(6):131721326.

[255]　FangJQ,YuXH.Chin.Phys.Lett.,2004,21(8):142921433.

[256]　FangJQ,LiuQ,LiY.TheFourthAsia2PacificWorkshoponChaosControlandSynchronization(4th

Apwccs),keynotespeech,Jin2QingFang,MultigoalcontrolofHalo2ChaosinBeamTransportNet2workswithSW/SF,proceedingsofAPWCCS2007,pp123,3282338.24226August2007atHarbin,China.

[257]　ParrishJK,ViscidoSV,GrunbaumD.Biol.Bull.,2002,202:2962305.[258]　BuhlJD,SumpterJT,CouzinID,etal.Science,2006,312:140221406.[259]　LevineH,RappelWJ,CohenI.Phys.Rev.E,2001,63:0171012017104.[260]　GaziV,PassinoKM.IEEETrans.Autom.Control,2003,48:6922697.[261]　ErdmannU,EbelingW,MikhailovAS.Phys.Rev.E,2005,71:0519042051907.[262]　ChenHY,LeungKT.Phys.Rev.E,2006,73:0561072056109.[263]　Olfati2SaberR.IEEETrans.Autom.Control,2006,51:4012420.

[264]　D’OrsognaMR,ChuangYL,BertozziAL,etal.Phys.Rev.Lett.,2006,96:1043022104304.[265]　TonerJ,TuY.Phys.Rev.Lett.,1995,75:432624329.[266]　TonerJ,TuY.Phys.Rev.E,1998,58:482824858.

[267]　AlexanderMogilnerLEK.J.Math.Biol.,1999,38:5342570.

[268]　TopazCM,BertozziAL,LewisMA.Bull.Math.Biol.,2006,68:160121623.[269]　ReynoldsCW.Comput.Graph.,1987,21:25234.

[270]　VicsekT,CzirókA,Ben2JacobE,etal.Phys.Rev.Lett.,1995,75:122621229.[271]　CzirokA,VicsekT.PhysicaA,2000,281:17229.

[272]　HuepeC,AldanaM.Phys.Rev.Lett.,2004,92:1687012168704.[273]　GrégoireG,ChatéH.Phys.Rev.Lett.2004,92:0257022025705.

[274]　JadbabaieA,LinJ,MorseAS.IEEETrans.Autom.Control,2003,48:98821001.[275]　SavkinAV.IEEETrans.Autom.Control,2004,49:9812983.[276]　MoreauL.IEEETrans.Autom.Control,2005,50:1692182.

[277]　CouzinD,KrauseJ,JamesR,etal.J.Theor.Biol.2002,218:1211.[278]　CouzinID,KrauseJ,FranksNR,etal.Nature,2005,433:5132517.[279]　AldanaM,HuepeC.J.Stat.Phys.,2003,112:135.

[280]　CzirokA,StanleyHE,VicsekT.J.Phys.A,1997,30:1375.[281]　GrégoireG,ChatéH,TuY.PhysicaD,2003,181:157.[282]　StrogatzSH.PhysicaD,2000,143:1.

[283]　LiW.andWangX.Phys.Rev.E,2007,75(3):021917.

[284]　YangW,CaoL,WangXandLiX.Phys.Rev.E,2006,74(9):037101.[285]　SuH,WangXandLinZ.IEEETrans.Autom.Control,submitted.[286]　VicsekT.Nature,2001,411:421.

[287]　AssisiCG,JirsaVK,KelsoJAS.Phys.Rev.Lett.,2005,94:0181062018110.[288]　Olfati2SaberR,MurrayRM.IEEETrans.AutomaticControl,2004,49:152021533.[289]　RenW,BeardRW.IEEETrans.AutomaticControl,2005,50:6552661.[290]　HuG,XiaoJ,GaoJ,etal.Phys.Rev.E.,2000,62:304323047.

446物理学进展　　27卷

[291]　WangXF,ChenGR.Phys.A.,2002,310,5212530.

[292]　LiX,WangXF,ChenGR.IEEETrans.CircuitsandSystems.2004,51:207422087.[293]　AlbertR,JeongH,BarabásiAL.Nature,2000,406:3782381.[294]　AxelrodR,HamiltonWD.Science,1981,211:139021396.

[295]　SugdenR.Theeconomicsofrights,Co2Operationandwelfare,Oxford:Blackwell,1986.[296]　ArthurWB.AmericanEconomicReview,1994,84(2):4062411.[297]　ArthurWB.Science,1999,284:1072109.

[298]　ChalletD,ZhangYC.PhysicaA,1997,246(324):4072418.

[299]　MaynardSmithJ.Evolutionandthetheoryofgames,Cambridge:CambridgeUniversityPress,1982.[300]　WeibullJW.EvolutionaryGameTheory,Cambridge:MITPress,1995.[301]　NowakM,MayR.Nature,1992,359:8262829.

[302]　JohnsonNF,HuiPM,JonsonR,etal.Phys.Rev.Lett.,1999,82(16):336023363.[303]　QuanHJ,WangBH,HuiPM,etal.PhysicaA,2003,321(122):3002308.[304]　MoelbertS,LosRiosP.PhysicaA,2002,303(122):2172225.[305]　BastollaU,ParisiG.PhysicaD,1998,115(324):2032218.[306]　KirleyM.PhysicaA,2006,365(2):5212528.[307]　ZhouT,WangBH,ZhouPL,etal.Phys.Rev.E,2005,72:046139.[308]　ShangLH,WangXF.PhysicaA,2007,377(2):6162624.[309]　ShangLH,WangXF.PhysicaA,2006,361(2):6432650.

[310]　KauffmanSA.TheOriginsofOrder2Self2OrganizationandSelectioninEvolution.Oxford:Oxford

UniversityPress,1993.

[311]　BornholdtS,SneppenK.Phys.Rev.Lett.1998,81(1):2362239.[312]　AlexanderJM.PhilosophyOfScience.2003,70(5):128921304.

[313]　WattaPB.WangKN,HassounMH.IEEETransactionsOnNeuralNetworks.1997,8(6):126821280.[314]　GalstyanA,LermanK.Phys.Rev.E.2002.66(1):015103.[315]　NowakMKS.Nature.1993,364:56258.

[316]　DoebeliM,HauertC.EcologyLetters,2005,8(7):7482766.[317]　SzabóG,FathG.PhysicsReports.(toappear).

[318]　NowakM,MayR.Int.J.Bifurcat.Chaos.1993,3:35278.[319]　SantosFC,PachecoJM.Phys.Rev.Lett.2005,95(9):098104.

[320]　SantosFC,PachecoJM,LenaertsT.ProceedingsOfTheNationalAcademyOfSciencesOfTheU2

nitedStatesOfAmerica.2006,103(9):349023494.

[321]　SantosFC,PachecoJM.JournalOfEvolutionaryBiology.2006,19(3):7262733.

[322]　Gomez2GardenesJ,CampilloM,FloriaLM,MorenoY.Phys.Rev.Lett.s.2007,98:108103.[323]　WuZX,XuXJ,HuangZG,WangSJ,WangYH.Phys.Rev.E.2006,74:021107.

[324]　HofbauerJandSigmundK.Evolutionarygamesandpopulationdynamics,Cambridge:CambridgeU2

niversityPress,1998.

[325]　HauertC,DoebeliM.Nature.2004,428:6432646.

[326]　HauertC,SzaboG.AmericanJournalOfPhysics.2005,73(5):4052414.[327]　ShangLH,LiX,WangXF.Euro.Phys.J.B.2006,54(3):3692373.[328]　WangWX,RenJ,ChenGR,WangBH.Phys.Rev.E2006,74:056113.[329]　NowakMA.Science.2006,314:156021563.

　4期方锦清等:一门崭新的交叉科学:网络科学(下篇)

447

[330]　SantosFC,PachecoJM,LenaertsT.PlosComputationalBiology.2006,2(10):128421291.[331]　PachecoJM,TraulsenA,NowakMA.Phys.Rev.Lett.s.2006,97:258103.

[332]　FangJQ,BiQ,LiY,etal.ScienceinChinaSeriesG:Physics,Mechanics&Astronomy,2007,

50(3):3792396.

[333]　FangJQ.ProgressinNatureScience,17(7):7612774.

[334]　FangJQ.BiQ,LiY,etal.AdvanceinComplexSystem,willappear2007[335]　FangJQ.BiQ,LiY,etal.ModernPhysicsB,17(7):7612774.[336]　方锦清.自然科学进展,2007,17(7):29247.

[337]　陈振毅,汪小帆.系统工程学报,2005,20(2):1322138[338]　郭维平,汪小帆,李翔.通信学报,2006,27(10):51256.[339]　ChenZY,WangXF.PhysicaA,2006,364:5952602.[340]　ChenZY,WangXF.Phys.Rev.E,2006,73(3):036107[341]　潘灶烽,汪小帆.物理学报,2006,55(8):405824064.[342]　宋丽雅,李翔,汪小帆.计算机仿真,2006,23(10):1032108.[343]　SongLY,LiX,WangXF.Proceedingofthe6thWCICA,June21223,2006,Dalian.16220.[344]　ShiMJ,LiX,WangXF.Proceedingofthe6thWCICA,June21223,2006,Dalian.26230.[345]　MaBJ,XiaoJH,YangJZ.DYNAMCONTDISSERB,200613(324):3712377.

[346]　周涛,汪秉宏.见:复杂网络(郭雷,许晓铭主编),第6章,pp1152140,上海:上海科技教育出版社,2006.[347]　潘灶烽,汪小帆,李翔.系统仿真学报,2006,18(8):234622348[348]　许丹,李翔,汪小帆.控制与决策,2006,21(7):817282.[349]　史明江,李翔,汪小帆.计算机工程与应用,2006,11:1102115.

[350]　JianXinWang,GengZhao.Security,2004,6:3162324,No70371068,70431002

[351]　张茉楠.创业型经济的发展战略及其模式研究,博士后开题报告,国家信息中心,2007,3.[352]　姜彦福,张帷.创业管理学,北京:清华大学出版社,2005.

[353]　杨德林.中国科技型创业家行为与成长,北京:清华大学出版社,2005.[354]　解绉,汪小帆.复杂系统与复杂性科学,2005,3:1212.[355]　ZhangP,LiMH,etal.arXiv:physics/0505064,2005.

[356]　GaoL,LiMH,WuJ,DiZR.DynamicsofContinuous,DiscreteandImpulsiveSystems,SeriesB,

2006,13(3):421.

[357]　WilsonEO.Consilience2TheUnityofKnowledge[M].NewYork:Knopf,1998.[358]　AmaralLAN,BarratA,BarabasiAL,etal.Euro.Phys.J.B.2004,38(2):1432145.[359]　章忠志.复杂网络演化模型研究.大连理工大学博士学位论文.2006.[360]　方锦清.自然杂志,2005,27(5):2712276.

[361]　方锦清.汪小帆,刘曾荣.科技导报,2004,188(2):9212.[362]　方锦清.科技导报,2006,24(222):67272.

[363]　方锦清.2005年全国复杂网络学术会议,特邀报告,中国武汉.2005.4.9.

[364]　方锦清主编.第二届全国复杂动态网络学术论坛论文集,中国高等科学技术中心.2005:1412146.[365]　方锦清.2005全国复杂系统学术论坛,特邀报告,北京,中国高等科学技术中心.2005.12.[366]　方锦清.第十四届全国凝聚态理论与统计物理学术会议,邀请报告,广州,2006.11.11213.[367]　方锦清.2006全国复杂网络会议,特邀报告,武汉,2006.12.16219.[368]　方锦清.复杂系统研究与系统生物学论坛,邀请报告,上海.2006.11.21224.[369]　方锦清.2007秋季物理学术会议,邀请报告,南京,2007.9.1829.20.

448物理学进展　　27卷

[370]　方锦清.考察我国高新技术产业网的发展,第三届海峡两岸统计物理学术会议,大会特邀报告,杭州和

金华:2007.11.11211.16.

[371]　方锦清.我国高新技术网的若干特点与思考———探索“中国硅谷”之路,2007年全国复杂网络会议,大

会特邀报告,上海:2007年12.1212.3.

[372]　李永,方锦清等.我国高技术产业网络的层次结构探讨,同上会议.[373]　刘强,方锦清等.我国电子信息技术百强企业网络的初步分析,同上会议.[374]　刘强,方锦清等.考察世界五百强企业网络的若干特点,同上会议.

[375]　孙伟刚,方锦清.李常平等.国家高新技术产业开发区网络的某些特点,同上会议.[376]　方锦清著.驾驭束晕与探索网络科学,北京:原子能出版社,将出版.

[377]　FangJQ.The4thInternationalWorkshopHcongzhou2007onSimnlationPhysics,InvitedSpeaker,

Nov.10211,2007,Hangzhou,China.

[378]　方锦清.略论网络科学的理论模型研究的若干进展与挑战,共建西部构和谐,张剑主编,中国西部杂志社,2007,35243.

NEWINTERDISCIPLINARYSCIENCE:

NETWORKSCIENCE(II)

FANGJin2qing1　WANGXiao2fan2　ZHENGZi2gang3

LIXiang2　DIZeng2ru3　BIQiao1

(Unitednetworkreseachgroupof“OneInstituteandTwoUniversities”

1:ChinaInstituteofAtomicEnergy,P.O.Box275-81,Beijing102413,P.R.China;

2:ShanghaiJiaoTongUniversity,Shanghai200240,P.R.China;

3:BeijingNormalUniversity,Beiijing100875,China)

Abstract:Thisreviewarticletriestosystematicallysummarizeandreviewnewinterdisci2plinaryscience2NetworkScience,whichisjustemergingintheendof20century,thecon2tentsofthissurveryincludeberifelydevelopmenthistory,basicconcepts,classificationofnetworkscience,theoreticalmodelsandtheirproperties,mainprogressesofsomeimportantissuesinnetworkscience.Mainbodyorganizationofthisarticleconsistsofpreface,12chap2tersandafterwords,whichdescribeandcoverresearchachievementsandadvancesofnet2workscienceinrecentyears.Especially,wesummarizeourownresultsandseveralprogres2sesinthisfieldsupportedbythekeyprogramofNationalNatureSciencefoundationofChina

forourunitednetworkreseachgroupof“OneInstituteandTwoUniversities”(ChinaInsti2tuteofAtomicEnergy,ShanghaiJiaoTongUniversityandBeijingNormalUniversity).Part(II)ofthissurveryincludesfromChapter7toChapter212andafterwords.

Keywords:Networkcompletesynchronization;partsynchronization;beamtransportnetworkwithsmall2worldandscale2free;multi2goalsynchronizationcontrol;flockinginswarm,networkgame,challengibleissuesandoutlook