基于频率权重耦合模型的多层复杂网络爆发式同步

金彦亮,姚 林,王 雪,罗雪涛

(上海大学通信与信息工程学院,上海 200444)

近几年,复杂网络科学逐渐成为热门研究领域.生活中的很多现象,如生物体的运作、电网、因特网系统的运作及失效等问题都可以抽象成复杂网络和系统进行研究[1].复杂网络模型是一种由节点和连接节点的边所组成的网络模型,现实生活中的许多复杂系统都可以抽象成复杂网络研究,因此在工程学、生物学、社会学和物理学等各个交叉领域中具有重要的地位[2-3].而同步是复杂网络研究中的一个重要分支,受到国内外大量学者的关注.目前,该领域提出了几种振子模型来研究复杂网络中的同步动力学行为,常见的有Kuramoto 振子模型[4]、Lorenz 振子模型[5]以及FitzHugh[6]和Nagumo[7]提出的FHN(FitzHugh-Nagumo)振子模型等.

爆发式同步(explosive synchronization)现象在2011 年由G´omez-Gardenes 等[8]研究无标度网络同步现象时发现,该现象指在同步化过程中,复杂系统从无序态到同步态的相变过程是不连续的一级相变过程.研究者们将这种动力学行为命名为爆发式同步,此后便掀起了相关研究热潮.Li等[9]研究了网络拓扑结构特征对爆发式同步产生的影响,发现节点的度与节点自身的本征频率呈正相关时会产生爆发式同步;Leyva等[10]研究了带权重的复杂网络模型上的爆发式同步现象;Zhang等[11]提出了一种以振子本征频率的绝对值为耦合权重的新模型,将对爆发式同步现象的研究从无标度拓扑结构推广到了一般性的网络拓扑结构.Zhou等[12]借助Kuramoto 频率权重耦合模型,研究了振子本征频率分布非对称情况下的爆发式同步现象.

目前,对复杂网络的研究大多是将复杂系统抽象成单个网络,而忽略了现实生活中不同复杂系统之间并非是独立的.复杂系统之间常常存在着相互作用依赖关系,如计算机网络中终端系统与服务器系统相互依存[13],电力基础设施中通信网与电网的交互控制等[14].因此,对复杂网络同步的研究从之前的单一网络逐渐开始转向多层网络,并逐渐成为近几年的研究热点.Su等[15]研究了双层共同演化网络下的爆发式同步,得出了两层网络中相互依赖关系的强弱会影响爆发式同步现象的产生;Zhang等[16]采用局域同步序参量作为调节振子耦合强度的控制变量,将爆发式同步扩展到了多层网络,并观察到了一级相变过渡到二级相变的现象;Jiang等[17]利用社区结构网络作为外力作用层直接作用于另一层网络来研究其对爆发式同步的影响;Kachhvah等[18]利用二阶Kuramoto 惯性模型研究了不同双层网络之间相互作用时对爆发式同步的影响.

单层网络上爆发式同步现象的研究给出了临界耦合强度[11],但对于多层复杂网络产生爆发式同步的临界耦合强度还未有相关结论.因此,本工作提出了一种层间一对一相连的由Kuramoto 振子所组成的多层频率权重耦合模型,并基于该模型研究了层间相互作用强度和网络平均度对爆发式同步的影响,通过理论推导出了多层复杂网络上爆发式同步的临界耦合强度,最后通过数值仿真验证了理论结果的正确性.多层频率权重耦合Kuromoto 模型的提出,有助于分析实际中各个领域的多层网络,例如由通信网和电网组成的智能电网;服务器系统和终端系统所组成的计算机网以及生物系统中由不同分子组成的多层神经网络等,这对于研究多层复杂网络具有很重要的意义.

1 Kuramoto 模型与同步判据

Kuramoto 模型是一种没有振幅、只有相位的经典相振子模型.网络中的各个振子在未受到耦合作用时以各自的本征频率独立运动,当振子之间以某些方式发生耦合时,系统便会随着耦合作用的增强达到同步状态.对于由N 个Kuramoto 振子组成的复杂系统,Kuramoto 耦合模型[8]为

式中:θi和ωi分别表示振子i 的相位和本征频率;λ 为振子之间的相互耦合强度;Aij为网络的邻接矩阵.

在复杂网络同步研究中,同步序参量R 是衡量整个系统中振子同步化的参量.对于N 个Kuramoto 振子的网络,序参量可定义为[8]

式中:ψ 为整个系统的平均相位.若把系统中所有振子的相位看作是分布在一个复平面单位圆上的点,每个振子的相位可用一个单位向量eiθj表示.当系统处于完全无序状态时,振子的相位均匀地分布在单位圆上,此时所有相位会相互抵消,R=0;随着系统同步程度的增加,序参量R 也随之增大;而当系统达到完全同步状态时,所有振子的单位向量均指向单位圆的某一点附近,此时序参量R=1.

为了进一步观察同步现象的细节,可以通过计算振子在长时间内的有效频率来判定网络同步化程度[8],

系统在未达到同步之前,每个振子会按照自身的本征频率自由运动,随着同步程度增加直至完全同步时,整个系统的有效频率会变成系统的平均频率.

2 多层复杂网络下的频率权重耦合模型

本工作研究了两层复杂网络上的爆发式同步.如图1 所示,每层网络由N 个节点组成;λ 表示同一层内节点之间的相互耦合作用强度;h 表示两层网络之间的相互作用强度.

图1 由随机网络构成的双层网络拓扑结构Fig.1 Topology structure of 2-layered network composed of random networks

对于由2N 个Kuramoto 振子所组成的双层网络,本工作提出一种双层频率权重耦合复杂网络模型,其动力学方程可表示为

式中:i=1,2,···,N;下标1 和2 分别表示两层网络;λ 为层内耦合强度;h 表示层间相互作用强度;ki,1和ki,2代表网络中各个振子自身的度;ωi,1和ωi,2为振子的本征频率;Mij表示多层网络邻接矩阵M 中的元素.M 定义为

式中:A1和A2分别表示第一层网络和第二层网络的邻接矩阵,若节点i 和j 间存在相互耦合作用,则两层网络的邻接矩阵元素Aij=1,反之Aij=0;I 表示单位矩阵,表示网络层与层之间是一对一相连的相互作用.为便于理论分析,两层网络均采用相同拓扑结构的随机网络,两层网络振子的本征频率分布采用同一种对称的钟形曲线分布,这样只需要研究其中一层网络的动态演化同步行为,并且随机网络振子自身的度可以通过平均度〈k〉来代替.因此,两层网络的频率权重耦合模型方程可表示为同一方程,

式中:P(k)为网络的度分布;ρ(k;θ,t)为网络振子相位的概率密度.非全连网下的序参量可定义为

根据自洽理论分析,利用平均场思想[20]将式(8)代入式(7)中可改写成如下平均场形式:

定义振子相位与网络平均相位的相位差为∆θi=θi−ψ,则平均场方程(9)可变为

由式(11)可知,相位差的正负与振子本征频率的正负相同,且随着耦合强度的增大逐步趋向于0.由于振子本征频率分布是关于0 对称的钟形分布,因此网络振子的相位仅与本征频率的正负有关,而与具体数值无关.由此,序参量的计算方程可改写为[9]

式中:∆θ+和∆θ−分别表示本征频率为正负时的同步团相位.将式(11)代入式(12)中可得,

进一步化简可得到关于序参量R 的一元四次方程为

通过稳定性分析,可得出方程(14)存在非0 解时耦合强度λ 的范围,

由此可得产生爆发式同步时的临界耦合强度为

由式(16)可知,爆发式同步的临界耦合强度λc与网络层间的相互作用强度h 和网络的平均度〈k〉有关.当层间相互作用强度为0 时,特例λc=2 便是单层网络下本征频率权重耦合模型产生爆发式同步时后项相变的临界耦合强度.通过方程(14)可得到序参量R 与耦合强度λ的关系为

由式(17)可以看出,R 与振子本征频率分布的具体函数没有任何关系,只要本征频率分布是关于0 对称的分布,产生爆发式同步的相变临界点就在λc=2(1+h/〈k〉)处.

3 数值模拟结果

为验证上述理论的正确性,本工作对双层频率权重耦合复杂网络进行了数值仿真验证.实验采用复杂随机网络中具有代表性的ER(Erdos-Renyi)网络作为网络模型.网络中振子的本征频率需采用单峰对称的钟形分布,因此选取洛伦兹分布(g(ω)=,ω0=0,γ=1)作为振子的本征频率分布,振子的初始相位为随机大小.当层间相互作用强度h 取不同数值时,使用四阶-龙格库塔算法分别对前后项相变的两个过程计算稳态值R(λ),其中前项相变的耦合强度λ 由λ0,λ0+∆λ,λ0+2∆λ,···,λ0+n∆λ 逐渐增加;反之,后项相变则由λ0+n∆λ 按对应耦合步长∆λ 逐渐减小到λ0.本工作参照文献[8]中的仿真参数,选取耦合步长∆λ=0.02.在对每个λ 值计算相应的稳态值时,需演化足够长的时间后抛去暂态,并对稳态的一段时间T 内取平均值,因此仿真中选取演化时间步长∆t=0.001,演化时间为5 000 时间步,稳态时间T=2 000[16].仿真实验具体步骤如下.

步骤1 选定网络振子的本征频率分布和初始相位;

步骤2 设定层内耦合强度λ 的取值范围、网络平均度〈k〉和层间相互作用强度h 的大小;

步骤3 设定λ 的步长∆λ 并进行该耦合强度值下的动态演化,待演化足够长时间后取平均时间T 记录稳态;

步骤4 计算出当前耦合强度λ 的稳态后,再将该稳态作为下一个新的耦合强度λ+∆λ情况的初态进行步骤3 的迭代;

步骤5 重复迭代步骤3 和4 并依次进行前项和后项同步相变的数值模拟,最后得出理论和仿真的拟合曲线.

图2 给出了平均度〈k〉=98.87 的ER 随机网络在不同层间相互作用强度h 下,序参量随耦合强度变化的同步相变图,图中R1F,R1B,R2F和R2B分别表示两层网络前项和后项的同步序参量.可以看出,图(a)∼(d)中均出现了明显的磁滞回线,两层网络随耦合强度的变化均产生了爆发式同步现象,且爆发式同步相变的临界点也发生了偏移,这说明层间相互作用强度h 的变化对爆发式同步造成了一定的影响.当h 取0,20,40 和60 时,后项相变的临界耦合强度值λc分别在1.96,2.36,2.76 和3.18 处,即随着h 的增大λc变大.该结果表明在多层网络的爆发式同步中,层间相互作用强度的增大会在一定程度上阻碍爆发式同步的产生.

图2 不同层间相互作用强度下序参量随耦合强度变化的同步相变图Fig.2 Synchronous phase transition diagram of order parameter variation with coupling strength under different inter-layer interaction strength

为进一步研究临界耦合强度与层间相互作用强度的关系,图3 给出了产生爆发式同步时的临界耦合强度λc和层间相互作用强度h 的关系,图中红色曲线为理论结果(λc=2(1+h/98.87)),蓝色圆圈代表层间相互作用强度h 在0∼100 间每间隔10 取值时,通过数值仿真得到的临界耦合强度λc,可以看出λc和h 呈线性递增关系.由表1 可以看到,仿真结果与理论数值之间的误差特别小,这说明实验结果与理论分析的结果相吻合,验证了理论结果的正确性.

图3 临界耦合强度随层间相互作用强度变化下理论结果与数值仿真的对比Fig.3 Comparison of theoretical results and numerical simulations that critical coupling strength variation with inter-layer interaction strength

表1 层间相互作用影响下临界耦合强度λc 理论值与实验值的误差Table 1 Deviation between the theoretical value and the experimental value of the critical coupling strength λc under the influence of inter-layer interaction strength

为进一步观察爆发式同步动力学行为的细节,本工作采用其中一层的有效频率对后项相变的同步相变情况进行表征分析.图4 给出了层间相互作用强度不变时不同平均度下有效频率随耦合强度λ 的后项同步相变情况.图4(a)∼(d)中h=50,平均度〈k〉分别为39.92,60.03,79.94 和100.11.可以看出,当耦合强度大于临界值时,所有振子的有效频率演变为整个系统的平均频率〈ω〉,因网络本征频率采用的是对称的钟形分布,故图中〈ω〉 ≈0.当耦合强度逐渐减小到同步临界值时,网络中每个振子的有效频率突然从整个系统的平均频率转变为自身的本征频率自由运动,整个系统从同步锁频态瞬间变成混乱无序态,并且可以看出是一个不连续的变化状态,这说明产生了爆发式同步.由图4(a)∼(d)还可以看出,临界耦合强度λc分别在4.62,3.64,3.20 和2.98 处,与平均度〈k〉呈递减关系.这说明在层间相互作用强度一定的情况下,网络平均度〈k〉越大越能促进爆发式同步的产生.

图4 不同网络平均度下有效频率随耦合强度向后项相变变化的同步相变图Fig.4 Synchronous phase transition diagram that effective frequency variation with coupling strength of backward phase transition in different average degree of networks

同样,为了体现临界耦合强度随网络平均度变化的关系,图5 给出了临界耦合强度λc与网络平均度〈k〉的关系,红色曲线为理论结果(λc=2(1+50/〈k〉)),蓝色星号代表网络平均度〈k〉在20∼120 间大致每间隔10 取值时数值仿真得到的λc.从图中可以明显观察到,临界耦合强度随着网络平均度〈k〉的增大而减小.表2 给出了不同网络平均度下理论值与实验值的绝对误差,可以看出随着网络平均度的增大,实验值逐渐接近理论值,这是由于系统涨落与网络平均度呈负相关,当网络平均度越大时,仿真结果会越来越接近理论结果[19].这说明通过仿真结果验证了理论中临界耦合强度与网络平均度呈反比关系的结论.

表2 不同网络平均度下临界耦合强度λc 理论值与实验值的误差Table 2 Deviation between the theoretical value and the experimental value of the critical coupling strength λc under the different average degree of networks

图5 临界耦合强度随网络平均度变化时理论结果和数值仿真的对比Fig.5 Comparison of theoretical results and numerical simulations that critical coupling strength variation with average degree of networks

在上述仿真分析中,只考虑了两层随机网络的情况.为验证模型在多层网络上的普适性,可以将双层模型推广到多层网络,并通过类似理论分析得出多层网络上的临界耦合强度λc=2(1+).结果发现与两层网络情况类似,临界耦合强度只由层间相互作用强度和网络平均度共同决定,并不受网络层数的影响.图6 展示了实验仿真中临界耦合强度λc随网络层数m 的变化情况,其中红色实线代表h=40,〈k〉=98.87 时的理论曲线,蓝色圆点代表不同层数下的临界耦合强度值.可以发现,随着网络层数m 的增加,临界耦合强度λc基本维持在理论值附近,并不随网络层数的增加而变化.这说明实验结果验证了所提出模型在多层网络上的普适性.

图6 临界耦合强度随网络层数变化时理论结果和数值仿真的对比Fig.6 Comparison of theoretical results and numerical simulations that critical coupling strength variation with number of network layers

4 结束语

本工作利用Kuramoto 振子所组成的多层频率权重耦合模型研究了多层网络上的爆发式同步这一集体动力学行为,并得出了多层网络上产生爆发式同步的临界耦合强度.通过理论分析,发现临界耦合强度由层间相互作用强度和网络平均度共同决定.在网络平均度固定的情况下,增大层间相互作用强度会在一定程度上阻碍爆发式同步的产生;同理,保持层间相互作用强度不变,网络平均度越大越容易产生爆发式同步;而增加网络层数并不影响爆发式同步的产生.上述结论对于研究现实生活中多层复杂网络同步的产生和控制具有重要的指导意义.