金彥亮,姚 林,王 雪,羅雪濤

(上海大學通信與信息工程學院,上海 200444)

近幾年,復雜網(wǎng)絡科學逐漸成為熱門研究領域.生活中的很多現(xiàn)象,如生物體的運作、電網(wǎng)、因特網(wǎng)系統(tǒng)的運作及失效等問題都可以抽象成復雜網(wǎng)絡和系統(tǒng)進行研究[1].復雜網(wǎng)絡模型是一種由節(jié)點和連接節(jié)點的邊所組成的網(wǎng)絡模型,現(xiàn)實生活中的許多復雜系統(tǒng)都可以抽象成復雜網(wǎng)絡研究,因此在工程學、生物學、社會學和物理學等各個交叉領域中具有重要的地位[2-3].而同步是復雜網(wǎng)絡研究中的一個重要分支,受到國內外大量學者的關注.目前,該領域提出了幾種振子模型來研究復雜網(wǎng)絡中的同步動力學行為,常見的有Kuramoto 振子模型[4]、Lorenz 振子模型[5]以及FitzHugh[6]和Nagumo[7]提出的FHN(FitzHugh-Nagumo)振子模型等.

爆發(fā)式同步(explosive synchronization)現(xiàn)象在2011 年由G′omez-Gardenes 等[8]研究無標度網(wǎng)絡同步現(xiàn)象時發(fā)現(xiàn),該現(xiàn)象指在同步化過程中,復雜系統(tǒng)從無序態(tài)到同步態(tài)的相變過程是不連續(xù)的一級相變過程.研究者們將這種動力學行為命名為爆發(fā)式同步,此后便掀起了相關研究熱潮.Li等[9]研究了網(wǎng)絡拓撲結構特征對爆發(fā)式同步產生的影響,發(fā)現(xiàn)節(jié)點的度與節(jié)點自身的本征頻率呈正相關時會產生爆發(fā)式同步;Leyva等[10]研究了帶權重的復雜網(wǎng)絡模型上的爆發(fā)式同步現(xiàn)象;Zhang等[11]提出了一種以振子本征頻率的絕對值為耦合權重的新模型,將對爆發(fā)式同步現(xiàn)象的研究從無標度拓撲結構推廣到了一般性的網(wǎng)絡拓撲結構.Zhou等[12]借助Kuramoto 頻率權重耦合模型,研究了振子本征頻率分布非對稱情況下的爆發(fā)式同步現(xiàn)象.

目前,對復雜網(wǎng)絡的研究大多是將復雜系統(tǒng)抽象成單個網(wǎng)絡,而忽略了現(xiàn)實生活中不同復雜系統(tǒng)之間并非是獨立的.復雜系統(tǒng)之間常常存在著相互作用依賴關系,如計算機網(wǎng)絡中終端系統(tǒng)與服務器系統(tǒng)相互依存[13],電力基礎設施中通信網(wǎng)與電網(wǎng)的交互控制等[14].因此,對復雜網(wǎng)絡同步的研究從之前的單一網(wǎng)絡逐漸開始轉向多層網(wǎng)絡,并逐漸成為近幾年的研究熱點.Su等[15]研究了雙層共同演化網(wǎng)絡下的爆發(fā)式同步,得出了兩層網(wǎng)絡中相互依賴關系的強弱會影響爆發(fā)式同步現(xiàn)象的產生;Zhang等[16]采用局域同步序參量作為調節(jié)振子耦合強度的控制變量,將爆發(fā)式同步擴展到了多層網(wǎng)絡,并觀察到了一級相變過渡到二級相變的現(xiàn)象;Jiang等[17]利用社區(qū)結構網(wǎng)絡作為外力作用層直接作用于另一層網(wǎng)絡來研究其對爆發(fā)式同步的影響;Kachhvah等[18]利用二階Kuramoto 慣性模型研究了不同雙層網(wǎng)絡之間相互作用時對爆發(fā)式同步的影響.

單層網(wǎng)絡上爆發(fā)式同步現(xiàn)象的研究給出了臨界耦合強度[11],但對于多層復雜網(wǎng)絡產生爆發(fā)式同步的臨界耦合強度還未有相關結論.因此,本工作提出了一種層間一對一相連的由Kuramoto 振子所組成的多層頻率權重耦合模型,并基于該模型研究了層間相互作用強度和網(wǎng)絡平均度對爆發(fā)式同步的影響,通過理論推導出了多層復雜網(wǎng)絡上爆發(fā)式同步的臨界耦合強度,最后通過數(shù)值仿真驗證了理論結果的正確性.多層頻率權重耦合Kuromoto 模型的提出,有助于分析實際中各個領域的多層網(wǎng)絡,例如由通信網(wǎng)和電網(wǎng)組成的智能電網(wǎng);服務器系統(tǒng)和終端系統(tǒng)所組成的計算機網(wǎng)以及生物系統(tǒng)中由不同分子組成的多層神經(jīng)網(wǎng)絡等,這對于研究多層復雜網(wǎng)絡具有很重要的意義.

1 Kuramoto 模型與同步判據(jù)

Kuramoto 模型是一種沒有振幅、只有相位的經(jīng)典相振子模型.網(wǎng)絡中的各個振子在未受到耦合作用時以各自的本征頻率獨立運動,當振子之間以某些方式發(fā)生耦合時,系統(tǒng)便會隨著耦合作用的增強達到同步狀態(tài).對于由N 個Kuramoto 振子組成的復雜系統(tǒng),Kuramoto 耦合模型[8]為

式中:θi和ωi分別表示振子i 的相位和本征頻率;λ 為振子之間的相互耦合強度;Aij為網(wǎng)絡的鄰接矩陣.

在復雜網(wǎng)絡同步研究中,同步序參量R 是衡量整個系統(tǒng)中振子同步化的參量.對于N 個Kuramoto 振子的網(wǎng)絡,序參量可定義為[8]

式中:ψ 為整個系統(tǒng)的平均相位.若把系統(tǒng)中所有振子的相位看作是分布在一個復平面單位圓上的點,每個振子的相位可用一個單位向量eiθj表示.當系統(tǒng)處于完全無序狀態(tài)時,振子的相位均勻地分布在單位圓上,此時所有相位會相互抵消,R=0;隨著系統(tǒng)同步程度的增加,序參量R 也隨之增大;而當系統(tǒng)達到完全同步狀態(tài)時,所有振子的單位向量均指向單位圓的某一點附近,此時序參量R=1.

為了進一步觀察同步現(xiàn)象的細節(jié),可以通過計算振子在長時間內的有效頻率來判定網(wǎng)絡同步化程度[8],

系統(tǒng)在未達到同步之前,每個振子會按照自身的本征頻率自由運動,隨著同步程度增加直至完全同步時,整個系統(tǒng)的有效頻率會變成系統(tǒng)的平均頻率.

2 多層復雜網(wǎng)絡下的頻率權重耦合模型

本工作研究了兩層復雜網(wǎng)絡上的爆發(fā)式同步.如圖1 所示,每層網(wǎng)絡由N 個節(jié)點組成;λ 表示同一層內節(jié)點之間的相互耦合作用強度;h 表示兩層網(wǎng)絡之間的相互作用強度.

圖1 由隨機網(wǎng)絡構成的雙層網(wǎng)絡拓撲結構Fig.1 Topology structure of 2-layered network composed of random networks

對于由2N 個Kuramoto 振子所組成的雙層網(wǎng)絡,本工作提出一種雙層頻率權重耦合復雜網(wǎng)絡模型,其動力學方程可表示為

式中:i=1,2,···,N;下標1 和2 分別表示兩層網(wǎng)絡;λ 為層內耦合強度;h 表示層間相互作用強度;ki,1和ki,2代表網(wǎng)絡中各個振子自身的度;ωi,1和ωi,2為振子的本征頻率;Mij表示多層網(wǎng)絡鄰接矩陣M 中的元素.M 定義為

式中:A1和A2分別表示第一層網(wǎng)絡和第二層網(wǎng)絡的鄰接矩陣,若節(jié)點i 和j 間存在相互耦合作用,則兩層網(wǎng)絡的鄰接矩陣元素Aij=1,反之Aij=0;I 表示單位矩陣,表示網(wǎng)絡層與層之間是一對一相連的相互作用.為便于理論分析,兩層網(wǎng)絡均采用相同拓撲結構的隨機網(wǎng)絡,兩層網(wǎng)絡振子的本征頻率分布采用同一種對稱的鐘形曲線分布,這樣只需要研究其中一層網(wǎng)絡的動態(tài)演化同步行為,并且隨機網(wǎng)絡振子自身的度可以通過平均度〈k〉來代替.因此,兩層網(wǎng)絡的頻率權重耦合模型方程可表示為同一方程,

式中:P(k)為網(wǎng)絡的度分布;ρ(k;θ,t)為網(wǎng)絡振子相位的概率密度.非全連網(wǎng)下的序參量可定義為

根據(jù)自洽理論分析,利用平均場思想[20]將式(8)代入式(7)中可改寫成如下平均場形式:

定義振子相位與網(wǎng)絡平均相位的相位差為?θi=θi?ψ,則平均場方程(9)可變?yōu)?/p>

由式(11)可知,相位差的正負與振子本征頻率的正負相同,且隨著耦合強度的增大逐步趨向于0.由于振子本征頻率分布是關于0 對稱的鐘形分布,因此網(wǎng)絡振子的相位僅與本征頻率的正負有關,而與具體數(shù)值無關.由此,序參量的計算方程可改寫為[9]

式中:?θ+和?θ?分別表示本征頻率為正負時的同步團相位.將式(11)代入式(12)中可得,

進一步化簡可得到關于序參量R 的一元四次方程為

通過穩(wěn)定性分析,可得出方程(14)存在非0 解時耦合強度λ 的范圍,

由此可得產生爆發(fā)式同步時的臨界耦合強度為

由式(16)可知,爆發(fā)式同步的臨界耦合強度λc與網(wǎng)絡層間的相互作用強度h 和網(wǎng)絡的平均度〈k〉有關.當層間相互作用強度為0 時,特例λc=2 便是單層網(wǎng)絡下本征頻率權重耦合模型產生爆發(fā)式同步時后項相變的臨界耦合強度.通過方程(14)可得到序參量R 與耦合強度λ的關系為

由式(17)可以看出,R 與振子本征頻率分布的具體函數(shù)沒有任何關系,只要本征頻率分布是關于0 對稱的分布,產生爆發(fā)式同步的相變臨界點就在λc=2(1+h/〈k〉)處.

3 數(shù)值模擬結果

為驗證上述理論的正確性,本工作對雙層頻率權重耦合復雜網(wǎng)絡進行了數(shù)值仿真驗證.實驗采用復雜隨機網(wǎng)絡中具有代表性的ER(Erdos-Renyi)網(wǎng)絡作為網(wǎng)絡模型.網(wǎng)絡中振子的本征頻率需采用單峰對稱的鐘形分布,因此選取洛倫茲分布(g(ω)=,ω0=0,γ=1)作為振子的本征頻率分布,振子的初始相位為隨機大小.當層間相互作用強度h 取不同數(shù)值時,使用四階-龍格庫塔算法分別對前后項相變的兩個過程計算穩(wěn)態(tài)值R(λ),其中前項相變的耦合強度λ 由λ0,λ0+?λ,λ0+2?λ,···,λ0+n?λ 逐漸增加;反之,后項相變則由λ0+n?λ 按對應耦合步長?λ 逐漸減小到λ0.本工作參照文獻[8]中的仿真參數(shù),選取耦合步長?λ=0.02.在對每個λ 值計算相應的穩(wěn)態(tài)值時,需演化足夠長的時間后拋去暫態(tài),并對穩(wěn)態(tài)的一段時間T 內取平均值,因此仿真中選取演化時間步長?t=0.001,演化時間為5 000 時間步,穩(wěn)態(tài)時間T=2 000[16].仿真實驗具體步驟如下.

步驟1 選定網(wǎng)絡振子的本征頻率分布和初始相位;

步驟2 設定層內耦合強度λ 的取值范圍、網(wǎng)絡平均度〈k〉和層間相互作用強度h 的大小;

步驟3 設定λ 的步長?λ 并進行該耦合強度值下的動態(tài)演化,待演化足夠長時間后取平均時間T 記錄穩(wěn)態(tài);

步驟4 計算出當前耦合強度λ 的穩(wěn)態(tài)后,再將該穩(wěn)態(tài)作為下一個新的耦合強度λ+?λ情況的初態(tài)進行步驟3 的迭代;

步驟5 重復迭代步驟3 和4 并依次進行前項和后項同步相變的數(shù)值模擬,最后得出理論和仿真的擬合曲線.

圖2 給出了平均度〈k〉=98.87 的ER 隨機網(wǎng)絡在不同層間相互作用強度h 下,序參量隨耦合強度變化的同步相變圖,圖中R1F,R1B,R2F和R2B分別表示兩層網(wǎng)絡前項和后項的同步序參量.可以看出,圖(a)～(d)中均出現(xiàn)了明顯的磁滯回線,兩層網(wǎng)絡隨耦合強度的變化均產生了爆發(fā)式同步現(xiàn)象,且爆發(fā)式同步相變的臨界點也發(fā)生了偏移,這說明層間相互作用強度h 的變化對爆發(fā)式同步造成了一定的影響.當h 取0,20,40 和60 時,后項相變的臨界耦合強度值λc分別在1.96,2.36,2.76 和3.18 處,即隨著h 的增大λc變大.該結果表明在多層網(wǎng)絡的爆發(fā)式同步中,層間相互作用強度的增大會在一定程度上阻礙爆發(fā)式同步的產生.

圖2 不同層間相互作用強度下序參量隨耦合強度變化的同步相變圖Fig.2 Synchronous phase transition diagram of order parameter variation with coupling strength under different inter-layer interaction strength

為進一步研究臨界耦合強度與層間相互作用強度的關系,圖3 給出了產生爆發(fā)式同步時的臨界耦合強度λc和層間相互作用強度h 的關系,圖中紅色曲線為理論結果(λc=2(1+h/98.87)),藍色圓圈代表層間相互作用強度h 在0～100 間每間隔10 取值時,通過數(shù)值仿真得到的臨界耦合強度λc,可以看出λc和h 呈線性遞增關系.由表1 可以看到,仿真結果與理論數(shù)值之間的誤差特別小,這說明實驗結果與理論分析的結果相吻合,驗證了理論結果的正確性.

圖3 臨界耦合強度隨層間相互作用強度變化下理論結果與數(shù)值仿真的對比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

為進一步觀察爆發(fā)式同步動力學行為的細節(jié),本工作采用其中一層的有效頻率對后項相變的同步相變情況進行表征分析.圖4 給出了層間相互作用強度不變時不同平均度下有效頻率隨耦合強度λ 的后項同步相變情況.圖4(a)～(d)中h=50,平均度〈k〉分別為39.92,60.03,79.94 和100.11.可以看出,當耦合強度大于臨界值時,所有振子的有效頻率演變?yōu)檎麄€系統(tǒng)的平均頻率〈ω〉,因網(wǎng)絡本征頻率采用的是對稱的鐘形分布,故圖中〈ω〉 ≈0.當耦合強度逐漸減小到同步臨界值時,網(wǎng)絡中每個振子的有效頻率突然從整個系統(tǒng)的平均頻率轉變?yōu)樽陨淼谋菊黝l率自由運動,整個系統(tǒng)從同步鎖頻態(tài)瞬間變成混亂無序態(tài),并且可以看出是一個不連續(xù)的變化狀態(tài),這說明產生了爆發(fā)式同步.由圖4(a)～(d)還可以看出,臨界耦合強度λc分別在4.62,3.64,3.20 和2.98 處,與平均度〈k〉呈遞減關系.這說明在層間相互作用強度一定的情況下,網(wǎng)絡平均度〈k〉越大越能促進爆發(fā)式同步的產生.

圖4 不同網(wǎng)絡平均度下有效頻率隨耦合強度向后項相變變化的同步相變圖Fig.4 Synchronous phase transition diagram that effective frequency variation with coupling strength of backward phase transition in different average degree of networks

同樣,為了體現(xiàn)臨界耦合強度隨網(wǎng)絡平均度變化的關系,圖5 給出了臨界耦合強度λc與網(wǎng)絡平均度〈k〉的關系,紅色曲線為理論結果(λc=2(1+50/〈k〉)),藍色星號代表網(wǎng)絡平均度〈k〉在20～120 間大致每間隔10 取值時數(shù)值仿真得到的λc.從圖中可以明顯觀察到,臨界耦合強度隨著網(wǎng)絡平均度〈k〉的增大而減小.表2 給出了不同網(wǎng)絡平均度下理論值與實驗值的絕對誤差,可以看出隨著網(wǎng)絡平均度的增大,實驗值逐漸接近理論值,這是由于系統(tǒng)漲落與網(wǎng)絡平均度呈負相關,當網(wǎng)絡平均度越大時,仿真結果會越來越接近理論結果[19].這說明通過仿真結果驗證了理論中臨界耦合強度與網(wǎng)絡平均度呈反比關系的結論.

表2 不同網(wǎng)絡平均度下臨界耦合強度λ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 臨界耦合強度隨網(wǎng)絡平均度變化時理論結果和數(shù)值仿真的對比Fig.5 Comparison of theoretical results and numerical simulations that critical coupling strength variation with average degree of networks

在上述仿真分析中,只考慮了兩層隨機網(wǎng)絡的情況.為驗證模型在多層網(wǎng)絡上的普適性,可以將雙層模型推廣到多層網(wǎng)絡,并通過類似理論分析得出多層網(wǎng)絡上的臨界耦合強度λc=2(1+).結果發(fā)現(xiàn)與兩層網(wǎng)絡情況類似,臨界耦合強度只由層間相互作用強度和網(wǎng)絡平均度共同決定,并不受網(wǎng)絡層數(shù)的影響.圖6 展示了實驗仿真中臨界耦合強度λc隨網(wǎng)絡層數(shù)m 的變化情況,其中紅色實線代表h=40,〈k〉=98.87 時的理論曲線,藍色圓點代表不同層數(shù)下的臨界耦合強度值.可以發(fā)現(xiàn),隨著網(wǎng)絡層數(shù)m 的增加,臨界耦合強度λc基本維持在理論值附近,并不隨網(wǎng)絡層數(shù)的增加而變化.這說明實驗結果驗證了所提出模型在多層網(wǎng)絡上的普適性.

圖6 臨界耦合強度隨網(wǎng)絡層數(shù)變化時理論結果和數(shù)值仿真的對比Fig.6 Comparison of theoretical results and numerical simulations that critical coupling strength variation with number of network layers

4 結束語

本工作利用Kuramoto 振子所組成的多層頻率權重耦合模型研究了多層網(wǎng)絡上的爆發(fā)式同步這一集體動力學行為,并得出了多層網(wǎng)絡上產生爆發(fā)式同步的臨界耦合強度.通過理論分析,發(fā)現(xiàn)臨界耦合強度由層間相互作用強度和網(wǎng)絡平均度共同決定.在網(wǎng)絡平均度固定的情況下,增大層間相互作用強度會在一定程度上阻礙爆發(fā)式同步的產生;同理,保持層間相互作用強度不變,網(wǎng)絡平均度越大越容易產生爆發(fā)式同步;而增加網(wǎng)絡層數(shù)并不影響爆發(fā)式同步的產生.上述結論對于研究現(xiàn)實生活中多層復雜網(wǎng)絡同步的產生和控制具有重要的指導意義.