非均勻Kuramoto振子系統(tǒng)的頻率同步

熊河浪，張 華

(重慶理工大學 理學院， 重慶 400054)

群體運動一直是近幾十年國內(nèi)外眾多學者研究的重要問題，其中耦合振子系統(tǒng)是研究熱點問題之一。17世紀物理學家惠更斯發(fā)現(xiàn)了兩鄰近鐘擺同步的振動現(xiàn)象。隨后，Kuramoto Y[1]提出的耦合振子模型解決了鐘擺的同步特性問題，并可以描述自然界中很多集群同步現(xiàn)象。振子集群同步現(xiàn)象在物理、生物、化學、社會等不同領域具有非常廣泛的應用，如在物理和化學中，用于自旋玻璃模型的建模和分析[2]、Josephson結陣列的同步現(xiàn)象[3]、原子反沖激光中的同步現(xiàn)象[4]、中微子的味道演化[5]、電化學振蕩器中的同步現(xiàn)象[6]、在互聯(lián)網(wǎng)連接網(wǎng)絡中各路由器同步發(fā)射信號時發(fā)生的網(wǎng)絡堵塞現(xiàn)象[7]；在生物中，心臟起搏器細胞[8]、晝夜節(jié)律[9]、神經(jīng)科學[10]、動物群集行為[11]、魚群涌現(xiàn)[12]等。

早期Kuramoto振子模型常是基于平均場耦合的，其解析處理方法一般基于平均場理論和自洽方法，如Thomas等[13]通過大型網(wǎng)絡上的廣泛數(shù)值模擬，系統(tǒng)地比較了非均勻度平均場(HMF)和淬滅平均場(QMF)方法的預測結果。Yook等[14]研究了Kuramoto振子在去同步階段的2個不同階參數(shù)的行為。Wesley等[15]分析了2個改變標準的易感-感染-易感(SIS)動力學模型。隨后有學者開始將Kuramoto振子模型轉化為梯度系統(tǒng)的形式，利用Lojasiewicz不等式理論分析該振子模型(梯度系統(tǒng)形式)的動力學行為[16]。

1階Kuramoto振子的模型還拓展到不同的連通拓撲效應、噪聲效應、混雜耦合形式的效應以及復雜網(wǎng)絡連接效應等多種復雜外部和內(nèi)部環(huán)境中考慮[17]。2階Kuramoto振子的模型可以拓展到固有頻率服從不同的分布、相位延遲等制約因素下考慮。常用于分析模擬動物群集行為的粒子模型[18]、結構保持的電力系統(tǒng)模型[19]、網(wǎng)絡簡化電力系統(tǒng)模型[20]、耦合節(jié)拍器[21]等。如Choi等[22]利用Gronwall不等式方法分析有限慣量Kuramoto振子的漸近完全相頻同步，Koditschek等[23]利用Lyapunov穩(wěn)定性理論證明了閉環(huán)機械控制系統(tǒng)的幾乎全局漸近穩(wěn)定性，D?rfler等[24]利用奇異攝動理論分析了搖擺方程與非均勻Kuramoto模型之間的關系等。然而，這些方法或者完全局限于2階Kuramoto振子，或者只適用于均勻慣性和單位阻尼，或者需要充分小的慣性系數(shù)。因此，本文給出了一個充分條件改進現(xiàn)有方法，使結果應用更加普遍。

本文提出了一個使2階Kuramoto振子達到頻率同步的充分條件；重新定義了振子的頻率同步和相位內(nèi)聚，將本文結果和奇異攝動方法得到的結果進行對比，發(fā)現(xiàn)在滿足一定的條件下慣性項對于同步的發(fā)生無影響，通過計算機數(shù)值模擬進一步驗證了充分條件的正確性。

1 預備知識

給定Rn空間中，元素全為1的向量記為1n，向量z∈Rn的2范數(shù)和∞范數(shù)分別記為‖z‖2和‖z‖∞；In表示n階單位矩陣[26]。

給定一個無向連通圖G=(V，E，A)，V= {1，2，…，n}表示所有頂點的集合，記e為總邊數(shù)， {1，2，…，e}=E?V×V表示頂點之間構成邊的集合，A=[aij]∈Rn×n表示圖G的鄰接矩陣，其中權重aij>0，且A=AΤ[26]。

如果diag({as}s∈{1，2，…，e})是邊權重的對角矩陣，那么對應的拉普拉斯矩陣為L=Bdiag({as}s∈{1，2，…，e})BΤ[26]。

Brouwer’s不動點定理：假設X?Rn是一個非空緊凸集，f∶X→X為連續(xù)映射，則存在x*∈X，使得f(x*)=x*[27]。

2 問題闡述

考慮連通圖G中具有n≥2個的2階Kuramoto振子系統(tǒng)：

(1)

記θ=[θ1，θ2，…，θn]Τ∈Rn，令Φ(θ)=BΤθ表示角的增量。同理，ω=[ω1，ω2，…，ωn]Τ∈Rn，Φ(ω)=BΤω表示固有頻率的增量。

定義1給定常數(shù)α∈[0，π/2)和上述θ，若||Φ(θ)||∞≤α成立，則稱系統(tǒng)(1)是相位內(nèi)聚的。

定義2對于任意充分小的正數(shù)ε，如果||Φ(ω)||∞≤ε成立，則稱系統(tǒng)(1)是頻率鎖相的或者是頻率同步的。

帶有慣性項的2階Kuramoto振子可以寫成如下梯度系統(tǒng)：

(2)

3 主要結果

引理1[26]無向連通圖G的拉普拉斯矩陣L是可對角化的，即L=Udiag(0，λ2，…，λn)UΤ，U∈Rn×n表示正交矩陣，令L?∈Rn×n是L的加號廣義逆，記：

L?=Udiag(0，1/λ2，…，1/λn)UΤ

(3)

證明：首先對式(1)的所有方程求和，得：

(4)

定義如下1階Kuramoto振子系統(tǒng)：

(5)

定義：

由于Φ(θ)=BΤθ，那么平衡方程可轉化為：

ω-Dωsync=Bdiag({as}s∈{1，2，…，e})sin(BTθ)=

Bdiag({as}s∈{1，2，…，e})y(BTθ)BTθ=

(6)

(7)

式(7)左邊是θ∈Ω的連續(xù)函數(shù)。若存在緊凸集Γ∞(α)={x∈Rn∶‖x‖∞≤α}，使x=BTθ∈Γ∞(α)，則式(5)存在一個平衡θ∈Ω。

||BTU(x)diag(0，γ2，…，γn)UT(x)(ω-Dωsync)||2≤

γ2||BTU(x)diag(0，1，…，1)UT(x)(ω-Dωsync)||2=

γ2||BT(ω-Dωsync)||2

(8)

對于x∈Γ2(α)，有：

(9)

當ν∈(0，π/2)時，y(ν)=sin(v)/v的值域為(0，1]，故λ2(L)≥λ2(L)y(ν)，因此1/λ2(L)≤ 1/(λ2(L)y(ν))。

又||x*||2≤α，則式(9)等價于：

(10)

即：

λ2(L)sin(α)≤λ2(L)

(11)

注1 在定理1條件下，當α∈[0，π/2)時，相位差滿足‖θi-θj‖∞<α；對于任意充分小的ε，固有頻率滿足‖ωi-ωj‖∞<ε，則式(1)處于同步平衡狀態(tài)。

奇異攝動分析[28]假設式(1)是過阻尼的，即Di>>Mi。定義擾動參數(shù)ε=Mmax/Dmin，其中：

Mmax=max{M1，M2，…，Mn}

Dmin=min{D1，D2，…，Dn}

或

在等式兩邊同時乘以ε，得：

(12)

在式(12)兩邊同時除以Mi，得：

(13)

(14)

4 數(shù)值模擬

考慮n=5的情形，設固有頻率初值ω=[0.01，0.02，0.03，0.04，0.05]T，阻尼系數(shù)初值D=[3，1，4，2，3]T，相位初值θ= [π/4，π/3，π/6,π/5,π/6]T，鄰接矩陣為：

系統(tǒng)(5)在時間[0，5]內(nèi)，利用奇異攝動分析得到振子的演化曲線，如圖1所示。

圖1 滿足條件Di>>Mi，時，振子同步曲線

系統(tǒng)(5)在時間[0，5]內(nèi)，模擬滿足定理1條件，得到振子的演化曲線，如圖2所示。

圖2 滿足條件時，振子同步曲線

系統(tǒng)(5)在區(qū)間[0，5]內(nèi)，模擬不滿足定理1條件，得到振子的演化曲線，如圖3所示。

圖3 不滿足條件時，振子不同步曲線

注2定理1的方法是通過構造緊凸集，運用不動點定理，得到慣性項，不會影響同步發(fā)生。奇異攝動理論是使阻尼系數(shù)遠大于慣性系數(shù)，通過計算找到一個充分條件，得出慣性項，會影響同步發(fā)生。因為慣性系數(shù)對于同步發(fā)生一直是有爭議的，在其他文章中也提到這個問題[29]。所以，通過2種情形的對比，在滿足一定的條件下，爭議是可以避免的。參考圖1、2，系統(tǒng)(5)最終會都會收斂到相同的數(shù)值。

下面驗證當振子數(shù)目n=2時，得出的結論是否與定理1給出的充分條件一致。模型如下：

(15)

(16)

其拉普拉斯矩陣為：

對應的特征值為λ1=0和λ2=2a。給出一個差分函數(shù)h∶R→R：

(17)

圖4 不同K值下的零點分布圖

5 結論

主要討論了2階Kuramoto振子的頻率同步問題，在滿足充分條件的情況下，所有振子之間的相位差都趨于零，即系統(tǒng)(1)達到同步。

提供的方法適用于1階和2階Kuramoto振子等價的情形，而在不等價的情形是否可以運用這個方法還需進一步研究。本文可以推廣到其他情形，如將固有頻率推廣到時變的固有頻率、在振子間的相互作用函數(shù)中加入相位延遲、慣性系數(shù)與阻尼系數(shù)的其他比率等。上述情形都是需要繼續(xù)研究的方向。