基于頻率權重耦合模型的多層復雜網絡爆發式同步

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

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