Kuramoto模型在Python和MATLAB中的简单实现

haidaowang

事情是这样的，我最近在研究团队编组及内部模式对发挥团队能力的影响，以及如何正确编组让团队能力发挥实现最大化，别问我为什么研究这个，反正稀里糊涂就研究上了。我发现在描述团队编组间及内部同步能力的时候，人们对Kuramoto模型（藏本模型）作了大量的研究,其中包括模型达到完全相位同步的充分条件、耦合强度对于同步的影响、一定条件下振子的收敛速率等。但具体实现一般都在MATLAB中，且网上代码过于复杂（我运行了一遍一堆报错），这里我使用Python和MATLAB对Kuramoto模型简单模拟。
模拟的话还是一遍举个栗子，边分析边测试效果最好，百度学术上有一篇关于Kuramoto模型的简单论文，我们就用它来实现模拟。
. |4 c. s) i- y1 \0 D7 ]3 h
- f1 @' ?/ ^9 W7 \% Z  K9 \/ u好吧，那我们开始，首先是Python实现：
( t2 m% c: O& h3 H) H, K, v
, D( q+ O) Y1 ^- c7 a* x) lN维Kuramoto模型的数学描述如下：
' a& n7 ?! h: F7 }5 q/ J1 ^; x
" `  z1 G' D( v) p% P6 T) H
没错，是讨厌的数学公式，没事，它可以改写成这样：
- J* R+ v8 p& }& Y9 c

3 N+ c' P2 D4 Y8 p; F& p( O好像还是有点长，那我们在改写一下：
" h) q* h9 v8 [0 \! h
) t- |" M- k7 v( S2 J$ g
; h, {5 t7 Q; `0 q  A6 S看着好多了，那我就来说说式子中参数的意义，Kij为耦合矩阵，是为了便于描述不同振子间耦合程度不同的情形。最下面那个式子的r就是我们的目标，反应振子间的相关性，这个相关性就可以描述我们想要的编组内部同步能力。
2 }( j0 W" @9 f0 c' k+ Z4 f
- y9 T. }# ^$ v( s5 p9 M哎呦，这个式子看起来好简单，这里要补充一下知识点：同步能力可不是一下子各组该怎么同步直接确定的了，它是一个从开始到稳定的阶段，也就是说随时间变化，最终反映在各组的同步能力才会确定，那么最后图像是什么样子才算同步能力好呢？
7 J$ v7 l  l/ ?1 h, p) k
7 a% F. }2 Y5 R同步能力好，是指随着时间的推移，各组的同步能力r逐渐稳定，波动现象消失或固定在某一个小范围内。需要注意的是这和各组r值之间的差距没有关系，我们要的是一个平稳的状态。那怎么办找r和t的关系呢？

注意看最上面那两个式子，相位（第一项，等号左边那个）上面有个点，这样他可就不简简单单是个相位了，它代表的是相位的变化值，是一个微小的微分值，好吧具体意思就是，那个式子左边展开之后是这样的：
8 U9 I, s& U# ?+ Y6 ]

哎呀，t出现了，其实\theta 与t有关，这里你可能有点绕，因为它们之间的关系是一一对应的，就是说每个时间的t对应了一个\theta ，我下面带入具体数值的时候你就知道了。

3 K- j4 R6 B4 o( [组间同步能力与时间t的关系出现了！
! c/ n& @+ C6 M( k6 H: h) F8 g/ q
2 [+ Q6 [2 W2 }' R3 v5 V也就是说我先用上面的那个公式4计算出来\theta 的值，在带入到公式5，那么t-r关系就可以明确下来了，那现在我们再回过头来看看文章中已经给的例子，看看还有没有未知量。

栗子是这样的栗子：
. W- u2 ^% K4 z5 [
* U! ^( J# V: |+ G+ @5 U8 z% ^假设某机构内部有 4 个编组，每个编组包括 5 个节点(其中 1 个节点为领导节点) 。另外，将上级领导作为一个独立的编组，且只包含一个节点。假设在领导机关增加4名信息传递人员。当以独立编组模式编组时，指定1名信息传递人员为指挥者，其指挥关系与其他编组一致; 当分散编组时，信息传递力量节点的关系与所在编组其他节点指挥关系一 致。其中，完全分散编组模式时，各信息传递力量节点之间无信息共享通道; 不完全分散编组时，在各信息传递人员节点之间建立一条信息共享通道。各编组模式及其拓扑结构如下图所示：

, P0 c! ?( L0 V4 k独立编组模式


: }8 C! N" {6 K: C6 Y3 a完全分散编组
$ L* z& x" B! j" N) U) J

3 Y1 X! H6 M" t9 K' o& w
不完全分散编组
" {5 F# D. H0 {6 ^. l& n


% {  {: Q( _% M+ p4 H$ i参数数据

, _( o& v0 t2 e) t0 B


& p- @* S& w; g参数确定一下有没有未知量：
2 a: b' A2 Q0 W9 ^- E: b8 V0 S
首先N，数据数目已知，这个有了。
2 ?; l7 l* H0 z! w# k
K值是分组内的连接强度，这个是看实际情况，由甲方提供或者自己看着给的，这里就是甲方给的编组图，i与j点的链接强度一目了然，这个有了。
4 }2 ]; F1 o: x/ R. e( u
4 N: j* g) N+ E- R0 ], K4 f( ]# c\omega _{i} 是振子i的固有频率，也称自然频率，甲方会给，没法自己估计，这个有了。
0 T5 T3 t# N- G
4 v2 U5 Y* ^2 G- y6 f. }9 k; i\theta ，\theta 怎么办，初始的\theta 会给，自己也能测的出来，但那么多\theta {j} -\theta{i} 得多少不知道啊，这里通过翻看文章，我发现其实文章是有一个特殊条件的，不然的话是需要研究耦合因子求三种约束条件解情况的，特殊条件就在这：
+ m% N7 I6 k0 ~( D) k
' H  M+ a5 n: F. v1 w; v( W+ Q' f假设编组内节点的初始相位差为π/2，且编号最小者为0，随编号增大而增大。

# N( H/ K% C; ^) }: ~哦，初始相位差知道了，你还告诉了我各个初始相位，那么\theta _{j} -\theta _{i} 的值就在一个范围内的几个固定值里面啊！
" `: i$ ]" [% b9 [' Z' e( L4 C3 Y
" i) I; i3 W* [+ H好的，没有未知量了，就是找K的时候麻烦点，没办法，这个决定了编组的不同，写脚本算一下吧：

#codingutf-8

##ScriptName:KuramotoSimulation.py
: Y' m- W3 {/ Q  N
% E9 b7 Y- J; m( kimport matplotlib.pyplot as plt
4 r6 C* N# m: E4 X/ C, H
from pylab import

  R  l8 r7 L# O; L. y4 Ifrom sympy import

from matplotlib.ticker import MultipleLocator, FormatStRFormatter
% Y  g1 n: D' \( N
  E2 U1 Q5 U- r1 |import math

8 e5 l" s1 m: n1 I, _% _1 Zimport numpy as np
% U) A( a4 w4 f
N = 31      #总节点数
9 W0 t- I; i. u# _6 @* H6 Y% |! ~0 L
/ m# [# o) P7 s! m( w9 {2 N) Jc=[0,0,math.pi  2,math.pi,3  math.pi 2,0,0,0,math.pi  2,math.pi,3  math.pi 2,0,0,0,math.pi  2,math.pi,3  math.pi 2,0,0,0,math.pi  2,math.pi,3  math.pi 2,0,0,0,math.pi  2,math.pi,3  math.pi 2,0,0]
3 c- y5 Q- C0 a9 P
% Q8 z  J. @* G6 r: ww = [4,3,3,2,2,1,4,3,3,2,2,1,4,3,3,2,2,1,4,3,3,2,2,1,4,3,3,2,2,1,5]
; J5 V1 a0 d; M* ^; V3 w" O3 _0 A
5 b1 @* n& }0 ^1 X% Hk1 = [

) G6 ]0 |, _1 {1 I[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
3 S: K/ l5 r' i: e/ S' f" Y
[1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2],
6 A3 `3 d: ^' D6 {4 e1 q
' i7 g% u9 i# M( _' }[1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
' p4 x) o: f3 P; y0 R; w& I8 H( C1 L
  n2 s' G& N) s3 d0 P9 W[1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

, i5 |$ C# j- E. J. J5 z- r[1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
+ U- j& E; l9 M, D+ r# o3 ]; v
; ~( x( \, @8 c3 c8 ~[1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
0 h: ~1 _& H2 `
[1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
2 R% R' V3 p1 m4 @
[0,0,0,0,0,0,1,0.9,0.9,0.9,0.9,0.9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2],

[0,0,0,0,0,0,0.9,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
9 f" x1 B) m4 f% t; N$ x
1 @3 z  V& P/ v[0,0,0,0,0,0,0.9,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  G8 t) Q6 ~+ m( X# K
[0,0,0,0,0,0,0.9,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
+ C& a! v4 G, X5 S( M0 P6 b
[0,0,0,0,0,0,0.9,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

4 k5 X: G$ k) M2 u. o- G. K[0,0,0,0,0,0,0.9,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
1 c3 _' y9 \8 Q8 d8 a) e" N5 q
! Z  D7 D; j  [; N) e* C[0,0,0,0,0,0,0,0,0,0,0,0,1,0.8,0.8,0.8,0.8,0.8,0,0,0,0,0,0,0,0,0,0,0,0,2],
8 w9 q. j$ I+ x( D: [8 |- O- ^
[0,0,0,0,0,0,0,0,0,0,0,0,0.8,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
, Y: K" Y) Z% l& B9 |$ J2 q
: v# u# J- g! i[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0],

& H2 Y! m- \! L[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0.7,0.7,0.7,0.7,0.7,0,0,0,0,0,0,2],

& y' a, _  T( B2 q[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,1,0,0,0,0,0,0,0,0,0,0,0],
/ z& A& f) {2 U* l, @5 S
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,1,0,0,0,0,0,0,0,0,0,0],
6 }% u- d: n+ z8 a+ S2 `+ S" y1 L
0 G6 f3 b/ ^. S, ?& P( w( w[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,1,0,0,0,0,0,0,0,0,0],

9 A; `  ~% m3 [0 ]  F0 }' _3 Z% \[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,0,1,0,0,0,0,0,0,0,0],
- p7 H# u" @3 O! U/ r! B. V" O+ ^
; U! |" h( s2 c& b* e3 e- k[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,0,0,1,0,0,0,0,0,0,0],
: u: C3 a: o% n, u
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0.6,0.6,0.6,0.6,0.6,2],

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,1,0,0,0,0,0],

/ S1 Q8 {5 J9 p[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,1,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,1,0,0,0],
# n$ {- ^+ m# s4 C1 n2 Y
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,1,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,0,1,0],
- ]  L1 g" ]; E" d$ Y7 @
  ~9 f7 V+ g: M+ {6 y0 x8 ][2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,1]]   
  |. y% R& v4 \5 F0 g5 o0 k- g
n = [i + 1 for i in range(22)]                #目标划分,24个值,1-24

t = [j  for j in range(1000)]

1 R% r8 e! _+ c/ N, Fci = 0

: t3 t6 {) t8 NC = []
  R3 ^5 L+ x% C, w' p, ]2 i" X
( ~1 [( e) P/ a. c3 KC.append(c)
# b! a" I5 l( z  Q$ l
( b  L/ z) j: k( {) x) \. N! zfor d in range(1100)
' ^# f  z: h+ V/ e- d$ i
3 f% N5 v- G7 Y; ?; J    for i in n
+ C8 _# B7 Q7 @7 N$ m
; E: \9 T; ?- `, h7 `/ R        for j in range(22)
$ D5 }3 f, j1 f6 V  z
3 Z/ q* g" C9 B: R' `$ D            cj = c[j + 1] - c
/ E3 X  d& g" p5 ]) a
            ci += k1[j + 1]  math.sin(cj)
8 H: S% W# \0 k( _
/ F3 f) Q  U" X( Z, B        cii = ci + w
, J, ]: M2 X% u
! c7 \& Q" {; i0 q  i. H" m8 N5 B        h = [0.01  cii + z for z in c]

* E% G& Z( h0 q+ Z0 @& Y  `; ^        C.append(h)
& ~1 ?& p, S1 ~# @* L8 x
        c = h
) g  N5 P3 q0 l. I4 I& t
def r_function(u)
9 e! v2 K& ?1 [' p, o0 j  u8 Z% w
    y1 = 0
0 u3 j9 C8 }* Q3 a# V
    y2 = 0

    y12 = 0

    y22 = 0
  e1 {: d3 s' |8 E3 w* f
    r = []

    for x in range(1000)

        for ul in u
, w9 d7 e4 F4 r0 o. z6 M) h1 X6 G
            y1 += math.cos(C[x][ul])
* E* ?1 d' M4 z# e+ R
& X& \0 {3 U$ c            y2 += math.sin(C[x][ul])

        y12 = y1  2

! l" Z9 W3 h, L1 n' ]& G# ]0 Z        y22 = y2  2
+ a' N' H; V4 M! W) G. s
% y$ Q1 n1 n- F/ W- N- I        r.append((float(1)  N )  ((y12 + y22)  0.5))
6 X4 W# l3 I' s! t: V; c
0 @8 J. S8 T; m6 ?$ h    return r
# o( F2 _! L; G* o! v; {" }
' r; M4 F0 m7 s! Dr1 = r_function([1,2,3,4,5,6])

r2 = r_function([7,8,9,10,11,12])
" L) g5 q% Q% a0 @; Q2 P& c0 ?
r3 = r_function([13,14,15,16,17,18])

r4 = r_function([19,20,21,22,23,24])
% t" _' s% d2 L0 n4 c
* c+ |2 k4 F/ e' g, H0 er5 = r_function([25,26,27,28,29,30])
0 n4 ?/ D3 r8 a& D) a9 Z4 O! ?
) B; l% {0 f+ p( {" ?2 H3 j9 B#r6 = r_function([31])


% y4 n; N0 [& o2 `2 ~- K
& T- ?+ G  W1 V% m. G& x* Oax = subplot(111) #注意一般都在ax中设置,不再plot中设置
3 d( d, {% P3 v' v3 e0 D" `+ R8 r
6 @! x2 Y$ @1 S! @$ b+ R% ZymajorLocator  = MultipleLocator(0.1) #将y轴主刻度标签设置为0.5的倍数
* {4 ^6 C9 f9 l
& I0 g' W2 l  W" {ax.yaxis.set_major_locator(ymajorLocator)

plt.plot(t, r1, marker='o', color='green', label='1')

plt.plot(t, r2, marker='D',color='red', label='2')
$ i1 P( m3 z: g+ E
( E. P$ D' h" c9 l, i, Aplt.plot(t, r3, marker='+',color='skyblue', label='3')
5 Q8 h6 ^: K+ Q; u2 V8 |) m9 Q
plt.plot(t, r4, marker='h', color='blue', label='4')

$ [. z: p+ j0 x# P6 aplt.plot(t, r5, marker='_',color='yellow', label='5')

#plt.plot(t, r6, color='red', label='6')

, G/ e9 X) W2 K( r$ R; i

plt.legend() # 显示图例

plt.xlabel('iteration times')
& _; i9 H9 |$ i7 B6 \# x& P5 D9 E
plt.ylabel('r')

plt.show()
* C7 M& b- n, i6 Q
, V9 l9 o" ]: B" B
8 }+ a5 X- e0 `% \独立编组结果如图：
; r9 V" m) g% b4 C
( [* c0 [5 x5 }' v

k1独立编组
7 k: F# ]: |& v, R$ c  b6 K) g& e7 G好的，从图像我们来看看Kuramoto模型在描述这个编组的时候，5组最终稳定，我们说这个团队编组还算科学，但我们改变一下K的值，换成分散编组：
# s, P( W/ M' j: i& _
k2 = [

' E* q% }) Z' ?2 C  h[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

/ C) d3 L  [) A[1,1,0.9,0.8,0.7,0.6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2],

[1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  a: T" H6 U0 v9 V! Q' \
[0.9,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

! S6 e- L* E7 p. H+ C[0.8,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

% G2 n/ d0 ]/ L0 E4 D[0.7,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  J* ~5 `& K0 H& J' f; x% C
9 j+ d: |2 n5 N$ Q[0.6,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

. [. }# x) s" O+ J+ |9 ~[0,0,0,0,0,0,1,1,0.9,0.8,0.7,0.6,0,0,0,0,0,0,0,0,0,0,0,2],

6 i9 ]9 I6 b6 ~4 O0 o[0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0.9,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
; A# Z* K) ?$ A4 g
/ P- p/ _0 g3 S# f3 \" r* o/ ~[0,0,0,0,0,0,0.8,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],

! S5 q& M& ^. C- q[0,0,0,0,0,0,0.7,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
: C  m7 U: U! R4 C3 `; A$ O
[0,0,0,0,0,0,0.6,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0],
- Z1 o5 b) Q8 a% v
& j  B  F8 }" ?6 R/ y& U, Q7 |[0,0,0,0,0,0,0,0,0,0,0,0,1,0.9,0.8,0.7,0.6,0,0,0,0,0,0,2],
) u" N5 V) R$ w, j) Z
[0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0.9,0,1,0,0,0,0,0,0,0,0,0],
' n! w4 u1 s2 v3 V! W  C
[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,1,0,0,0,0,0,0,0,0,0],
( ?+ x' y/ q% g) H( m" J+ O
[0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,0,1,0,0,0,0,0,0,0,0],

7 s6 ~# [4 z9 r6 ?% m; j2 W[0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,0,1,0,0,0,0,0,0,0],
0 q' Z( [) n, o$ Z8 V
, J9 i% }" |/ S1 w[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0.9,0.8,0.7,0.6,2],
! ~) V3 h& @* d# t' _5 q- t
/ Z' ?) V3 L* k$ X  S! f' [[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.9,0,1,0,0,0,0],
$ |; i  {* c, Z2 T
. g+ z' q0 z6 |( f) o[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,1,0,0,0],

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,0,1,0,0],
7 g. i- [. ~* }
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,0,1,0],
* c  q- j( _2 j+ j( u8 j/ i( x
[2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,1]
' C+ {/ y6 t; |; k+ I* j
3 I5 n/ |0 I6 H3 Z* s
8 s6 D& ~; D* z" A& I
5 u, {( i, I7 V3 R# H0 Q6 U/ @9 nk2独立编组
1 r3 u* ]2 `0 _7 d5 Y: o! Y% W# d3 o- B


2 o" W$ m8 J3 f" T2 g* U这两幅图像都在开始阶段大幅波动，而后在一定范围内趋于稳定，那么到底哪个分组模式最符合实际，最能突出编组能力呢？
- z/ T1 d* s0 Q4 f
- Q$ o5 L1 s6 t4 G2 _$ y4 u这里还有一个公式，来解决这个问题，编组同步能力的量化：


+ i8 b8 i  \3 l( EM_{s} 就可以描述某编组的同步效果，r_{s} 是达到稳定状态后序参量的均值，β∈(0,1)是调节因子。我们可以用M_{s} 比较编组内部的好坏。那编组间能力的好坏怎样比较呢？
5 @, A% V) T7 o6 G
9 f) R3 E$ y, d3 `7 ~9 a这个Kuramoto模型同样有所考虑，它有一个描述整个系统编组能力的公式：
/ L/ W# j5 Q! T% o

其中，P是编组的数量，M_{i} 是第i个编组的同步能力，\omega _{i} 是编组在整个系统中的权重，\psi {e} 是各编组平均 相位的均值，\Delta \psi 是各编组平均相位的标准差。具体的计算不是这篇文章的重点，就不在计算M{s} 和M的值来比较上述例子独立编组和分散编组的好坏了，本篇文章主要是讲下Kuramoto模型的解决思路，尤其是上面解决\theta 值的方法可以套用在其他Kuramoto模型中，做一个目标估计绰绰有余的。
, w9 [$ ~# p/ g" o' L4 k
2 l9 d. y' h! U+ J( V下面是解决Kuramoto模型常用的MATLAB编程方法，具体思路与上述基本一致，这里不再赘述，K的值我们给另一种编组模式：不完全分散编组模式，也是现在实际上最长用的编组方式，直接上代码：
' i1 C0 B. [/ K. I: P; g; z( K
clc;
0 ]7 O, w5 n" R/ \
# W. j' X0 y: s4 r1 t1 hclear all;
. m8 b- r, v0 ?4 S
k=[0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0.5;

( i) E4 ?& a$ R3 W$ t  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
; {7 n& m' j- c/ @6 P
  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
6 u0 y; u8 Q) t
! p5 H6 X, S2 \& k2 c8 m; c2 ]  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0.5;
) H7 o0 k' H3 U  U/ L% x6 ?  V
  |+ t7 G: G+ a2 N  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
4 Q/ g- V8 W" t3 o$ u
  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

& ]4 O& N8 ^5 \3 {  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
! w) ^6 Q( p  c( I! ^  K  s
% D1 ^) E- e. [" f  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0.5;

2 I, q1 V! ^5 J* j* _  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
) l, S1 N  D% }
  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
- \9 y2 X' U/ N3 Y: [9 C
  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   
6 i+ F- p& l3 K+ v3 s* u
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0.5;
! e& R6 g9 _) W5 g' A# J" B
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0;   

) ]% k8 ]0 X5 h1 c9 q$ U  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0;   
/ G& N7 @8 ?. w/ L
! g( y2 f1 m; X" L% u! ^! d  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0;   
; [# `/ u4 `& n3 Z$ k6 r+ W0 r
& G5 [. T% v, ^6 ^8 k# I+ @  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  2  2  0.5;

! |7 d; @" w# P  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  0  0  0  0;   

1 H$ B* F: N+ X5 P  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  0  0  0  0;   
3 _* v! v) w" G$ |, e9 V% P
7 Q' r2 C- q2 X+ t2 s6 B* S' s  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  0  0  0  0;   

: Y- e' H  r$ F& z3 l+ ]; G- s  0.5 0 0  0  0  0.5 0 0  0  0 0.5 0  0  0  0  0.5 0 0  0  0  0.5 0 0  0  0;   

  ]     
& V+ N& N$ D' @2 d3 }
/ r7 \7 j5 r+ [N = 25;

Step= 10000;

Theta=zeros(Step,N);
5 P( U0 w/ L( B, p4 k$ \6 A1 r
Omega =zeros(Step,N);
/ x2 U; \4 i3 [
DeltaT=0.01
6 s) D: b4 ~. b( J) m2 Q, ^
Theta(1,：）=[0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0];

; _6 h) n+ `2 @' U4 ~7 P4 Y- nOmega(1,：D=2*pi*[2  3  3  4      4 2    3  3      4 4    2  3      3 4    4  2      3 3    4  4      2 2    3  4      1];

%  求解微分方程

for n = 1:Step

    for i = 1:N

: G$ f( w4 V3 s& j9 C            Sigma = 0;

        for j = 1:N

: ^6 A: }- y# i* K            %%判断k对同步效果的影响
6 d, B! j0 c; h4 M( C1 M/ N
= Sigma + k(i,j) * sin( Theta(n,i) - Theta(n,j) );
, Z: _5 w$ P- {" X/ n
# H. p0 }8 Q/ R+ T5 @/ C0 l' l: I0 Q        end

            Omega(n+1,i) = Omega(n,i) - Sigma;

, y5 L% J% ?: }            Theta(n+1,i) = Omega(n+1,i)*DeltaT + Theta(n,i);
6 Y3 C/ l5 c+ [* H- {' X
    end
6 x; a2 O5 r  b- T6 v
end
3 X* @0 `; b# q  b  {7 G
9 u( M$ s  t( f" y* T7 `$ N& {%  求解序参量r(t)

groupN = 5;
7 S7 o* h; f, Z4 c$ ^
Step = 1000;
1 e0 t! T. R7 X
- Y5 B- L2 G4 efor i = 1:Step
- ~  v  n7 j/ D% d# l3 Y* I
    sigma1 = 0;

' N7 C8 [5 u1 n5 S& j% I3 `6 Q8 d3 l    sigma2 = 0;

) ~) _- M' u5 p$ q. `& d    sumTheta = 0;

    for j =6:10  %表示相应的群组

        sigma1=sigma1+cos(Theta(i,j));

        sigma2=sigma2+sin(Theta(i,j));

. r: d: Q# `. w  K' {% F    end
$ R. b. i8 Z. P
    r(i)= sqrt( sigma1^2+sigma2^2 )/groupN;

) q3 D) r" e9 l% R# G; B9 iend

x= 0.01：DeltaT：DeltaT*Step;
) |% O( w# r4 l, G$ s) |8 |4 F! v
' Z$ n) {5 Y2 q3 fplot(x,r);
6 y) y7 x8 O9 bMATLAB版不完全分散编组结果如下图：
2 O- n* s  d7 X0 k; {+ Q

regngfpcb

Kuramoto模型在Python和MATLAB中的简单实现