EDA365电子论坛网

标题: Kuramoto模型在Python和MATLAB中的简单实现 [打印本页]
作者: haidaowang    时间: 2021-2-22 10:03
标题: Kuramoto模型在Python和MATLAB中的简单实现
事情是这样的，我最近在研究团队编组及内部模式对发挥团队能力的影响，以及如何正确编组让团队能力发挥实现最大化，别问我为什么研究这个，反正稀里糊涂就研究上了。我发现在描述团队编组间及内部同步能力的时候，人们对Kuramoto模型（藏本模型）作了大量的研究,其中包括模型达到完全相位同步的充分条件、耦合强度对于同步的影响、一定条件下振子的收敛速率等。但具体实现一般都在MATLAB中，且网上代码过于复杂（我运行了一遍一堆报错），这里我使用Python和MATLAB对Kuramoto模型简单模拟。 ' a( X5 A  Y# j
模拟的话还是一遍举个栗子，边分析边测试效果最好，百度学术上有一篇关于Kuramoto模型的简单论文，我们就用它来实现模拟。
; W  c' L4 z( R& ^. Q5 B' v+ A" r2 I$ L; f" r% q
好吧，那我们开始，首先是Python实现：
% w: y4 j6 D$ K% K$ q6 ^( y# f. t" p% s6 c; i4 e# L" u$ F# |# e
N维Kuramoto模型的数学描述如下：
9 Z; h2 O; H1 K! k% W
8 S7 Q0 x% S6 K9 y. {! Z) X9 E# H, `0 E$ `- Z" R
没错，是讨厌的数学公式，没事，它可以改写成这样：
7 M$ D( m2 x" T# A! I 5 S2 a* C2 Y0 v( W' t$ w
; v: l7 k4 T1 k
好像还是有点长，那我们在改写一下：: c4 l: o7 u5 s, D0 y6 @/ J
% i! \3 W! U0 `7 G

* P* V7 w; p6 E/ K2 u. X- o0 e5 v+ I看着好多了，那我就来说说式子中参数的意义，Kij为耦合矩阵，是为了便于描述不同振子间耦合程度不同的情形。最下面那个式子的r就是我们的目标，反应振子间的相关性，这个相关性就可以描述我们想要的编组内部同步能力。& Y: a$ M; l3 Q/ a
6 S  A( ~$ j6 p3 v+ J# w/ O
哎呦，这个式子看起来好简单，这里要补充一下知识点：同步能力可不是一下子各组该怎么同步直接确定的了，它是一个从开始到稳定的阶段，也就是说随时间变化，最终反映在各组的同步能力才会确定，那么最后图像是什么样子才算同步能力好呢？5 ~- T1 V5 J( V: B- [" O
6 [( W& K! B* u$ n
同步能力好，是指随着时间的推移，各组的同步能力r逐渐稳定，波动现象消失或固定在某一个小范围内。需要注意的是这和各组r值之间的差距没有关系，我们要的是一个平稳的状态。那怎么办找r和t的关系呢？1 O- r3 l) Z0 y3 h; _7 O! ]

1 L* z% @) q: z0 ?3 G注意看最上面那两个式子，相位（第一项，等号左边那个）上面有个点，这样他可就不简简单单是个相位了，它代表的是相位的变化值，是一个微小的微分值，好吧具体意思就是，那个式子左边展开之后是这样的：
1 N$ W( h: Q9 e9 i* M   T* P; |, b5 Y  M2 L/ X5 z* U( y

% \- z2 x% y/ |4 E+ S+ }$ y; T哎呀，t出现了，其实\theta 与t有关，这里你可能有点绕，因为它们之间的关系是一一对应的，就是说每个时间的t对应了一个\theta ，我下面带入具体数值的时候你就知道了。& Y  l5 j: t& V, D3 t

& ^) y7 W; G7 G' J组间同步能力与时间t的关系出现了！9 [6 m: e' r- j

. Y: r3 \0 `/ s9 a也就是说我先用上面的那个公式4计算出来\theta 的值，在带入到公式5，那么t-r关系就可以明确下来了，那现在我们再回过头来看看文章中已经给的例子，看看还有没有未知量。
) X2 \5 |# d5 X$ I
3 B; z6 W) e: C栗子是这样的栗子：' z/ f8 g+ ^$ g( a
/ a- G5 }$ Q3 X; x/ `7 N
假设某机构内部有 4 个编组，每个编组包括 5 个节点(其中 1 个节点为领导节点) 。另外，将上级领导作为一个独立的编组，且只包含一个节点。假设在领导机关增加4名信息传递人员。当以独立编组模式编组时，指定1名信息传递人员为指挥者，其指挥关系与其他编组一致; 当分散编组时，信息传递力量节点的关系与所在编组其他节点指挥关系一 致。其中，完全分散编组模式时，各信息传递力量节点之间无信息共享通道; 不完全分散编组时，在各信息传递人员节点之间建立一条信息共享通道。各编组模式及其拓扑结构如下图所示：% x) L4 }& E0 g
* `3 i4 m/ g2 W8 U7 l
独立编组模式0 j+ \7 D# Y8 d; M* [1 d2 @1 \7 ]

! C1 \# e% C, ~# T8 L6 V* a" W: b# C+ K3 p
完全分散编组
/ z* \) R& E: X! a2 \8 |+ l: I ' h5 D  ]# f" l0 M' Z1 Q/ L, q* n
6 n' B1 u0 q. X  ~( C4 r

: Y) p; y2 U2 Z7 ~' Y& {不完全分散编组  ~. i' v- v1 U5 U# d/ _% r

% I; ?& Z0 l8 \% @5 t8 F4 o/ n# ?  g' @- K6 j, y9 R
* H& M8 N3 ]& o
参数数据
  ]( p5 J; \0 G3 R" E
( c! z& I5 S5 N. {. D : _* e0 R+ S" c4 u/ d
4 P/ P  D* B* Y5 ~. f! d. n
; S+ _: v0 H) x9 r3 i8 e0 ^4 W
参数确定一下有没有未知量：0 t5 K2 J. r/ b. g. n1 ^; @3 V  n
6 \# g! V) p! |3 W2 [
首先N，数据数目已知，这个有了。
7 ]8 s& q- e3 d+ X
' s. l# g1 G5 f! H* y& eK值是分组内的连接强度，这个是看实际情况，由甲方提供或者自己看着给的，这里就是甲方给的编组图，i与j点的链接强度一目了然，这个有了。/ w- H* L7 y8 }5 B0 u1 z/ X

. d1 R/ l2 F$ _\omega _{i} 是振子i的固有频率，也称自然频率，甲方会给，没法自己估计，这个有了。
7 u- V# H: ?0 `) M+ J2 k8 y- a7 S% t8 W: s4 A5 M# _
\theta ，\theta 怎么办，初始的\theta 会给，自己也能测的出来，但那么多\theta {j} -\theta{i} 得多少不知道啊，这里通过翻看文章，我发现其实文章是有一个特殊条件的，不然的话是需要研究耦合因子求三种约束条件解情况的，特殊条件就在这：
: C1 _8 W1 n: \
9 r6 A# J3 m  L假设编组内节点的初始相位差为π/2，且编号最小者为0，随编号增大而增大。
. w9 G; [: }/ F, t* Q0 d! c$ r. W4 `9 w0 i9 ^
哦，初始相位差知道了，你还告诉了我各个初始相位，那么\theta _{j} -\theta _{i} 的值就在一个范围内的几个固定值里面啊！8 j; `- j7 P  F

+ H3 z7 s) W5 S" m6 t/ j5 ]4 |3 Q. x好的，没有未知量了，就是找K的时候麻烦点，没办法，这个决定了编组的不同，写脚本算一下吧：/ J0 ^) `9 Z6 b# n
4 u) @: ]- D/ \$ M( C8 q
#codingutf-8
' ^" Z, \# \: V6 l& m% `7 f+ f% z9 s) Y6 n& U
##ScriptName:KuramotoSimulation.py
) y3 J/ n  |+ x* i' E" V" w/ N' z
import matplotlib.pyplot as plt & x0 G/ \6 H) p' X- A

! {  G: V3 L: O1 [% c3 A* Gfrom pylab import
  g' g( a6 k2 H: n; X
/ V+ @* V$ {  i0 E9 L. c- Pfrom sympy import 9 m! A& _# b3 O% p" B
3 O1 k% x+ G4 G8 y
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
2 Y0 P: X3 y- Z2 @7 N
% z2 {" i6 Y7 l! _: O2 gimport math
7 k" Z* Y& u2 k1 i" d  p
' Z: t# x& j' d$ D8 kimport numpy as np
" B, Q" K' B, W: B1 t, z: M' ]6 Y0 T+ P$ g4 y* c7 C* b% P
N = 31      #总节点数 ' @9 F" `4 h" S+ R' B8 c" ~

; `! y; E  |) p3 ic=[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] 2 d( A# J: a, _! Y  J: Z, N
# P. V9 z% T5 F
w = [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]
% t' I& s' i, C; |# A+ i0 ?' W/ m
2 D; Z& L. o* f- d6 P. s  D! w' v9 Uk1 = [ 3 H9 D9 R$ C! x% y

5 v$ c2 L1 z: P4 c; u7 u6 ^[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],
- x4 M  a: O% O4 c4 n
( Z( Z3 A, h. v* G[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], % v  R5 [% o6 [( S3 y
) Y6 v6 K9 z  W8 O& T3 m7 [
[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],
+ T$ V1 H0 t. }7 X' K$ R* G6 E. N& Z: r9 C/ E, w' _5 N
[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], ' i* y1 n5 s( Y( C* Q

7 s. S. _( r. g# Z9 W[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], 0 j$ R; ]( V& b7 A

. }  |& e9 Z1 |9 o) i/ W/ ^0 b[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], " Z* P5 z) ~. x# |2 D

8 q6 O- g$ _; W5 x[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],
4 `0 ?/ ~0 Z5 `! `3 _/ ~& C3 K5 F1 |
) p5 c: b& s' y/ T8 y$ u7 }% s[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],
: F* \/ e( u: r  l; i5 g% d. S; m8 i# F# e3 C" h- K& l
[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],
) U, `' E& h  Q8 b. T. c2 ]5 w. ]5 w: S1 D8 D
[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],
  x/ _+ R, O$ N" L& `
8 Q4 R1 p0 }! h6 a[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],
2 ]' |' }) [! [5 m5 m0 P% e& G8 H+ v  ]5 A
[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], 9 q* d/ U8 t: D& x2 H1 D
7 i# W5 [" p* k1 {. |! S8 A
[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], . d; N: p. S3 K9 T0 I' s

  ~1 K4 l9 }$ T& R- h& ~) P[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], ; M) _6 ?* v  ]2 l* R. d$ {

, M. o" c0 G! D2 {( @# a9 j/ m6 }[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],
+ G* G( ~3 W, Q2 x' T! E9 f- V. N' r0 K8 e( p+ s3 \
[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],
3 i% {# T& a/ o; F* m
' G! ^) u& m* m/ D* G! o) l' t) X[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], 7 w2 @: J  _3 b0 q0 K

( ?; ^0 G% l' a2 H[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], " [8 r( \# E* v" G8 c; z9 g) r

3 Y; T, e/ R$ w. |[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],
9 p( t; A; k, ]$ z) @
% ?3 q" d! K: m# s0 |& i1 r[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], : u. \. L7 A5 M/ v

' v; ]: |/ k" \+ D+ [; z[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],
# V% g6 @8 x: A; h5 ?: i- u
& ?  I. a/ a; c+ f[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],
# C2 V+ ]; S9 l. E2 u. K
+ ^% t7 z3 `" W2 ~& x( v[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],
; R2 C+ r. n+ t! U+ |; p( x2 `( S) e/ ]$ l
[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],
8 {$ |' N" X0 `+ }7 `. o7 X$ Q( F) ?7 P7 u* C$ q- Q8 ^1 B7 s$ q* `$ S7 a
[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],
+ E" _) r) I- I3 F* x2 |% ~) f+ ^1 ]6 g
[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],
% S+ g9 n5 T1 d  y; @. o! m" o4 e  _7 ^( M2 |  }- L7 n( \
[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],
+ t( J- o4 ?) J% i% b- C
$ I2 g8 I! [% j5 J* r5 Y% 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.6,0,1,0,0,0,0],
; X/ |9 c' e3 L5 H  \" h) r9 X  V5 c, m# 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,0.6,0,0,1,0,0,0],
9 S/ o8 [+ R( D" o! ~
4 T  g, w5 G6 Y% E5 k3 i) i0 {[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],
7 l% ?* s3 x9 g5 y) ~. S7 L. H+ u5 d' I& V; w- r5 e$ v' W3 g% |
[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],
. Y2 U3 o, ?: ^. j' B0 c
# H$ y9 k$ e/ R% r! {[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]]   
+ _; Y1 M( m2 W
# G4 O  M4 p& H" Tn = [i + 1 for i in range(22)]                #目标划分,24个值,1-24
+ H$ v/ D  w' ]
% D4 ]! l$ [) G2 V/ A+ ht = [j  for j in range(1000)]
& v+ r  X# L5 ~
! @- \2 m" U# {" m; A% bci = 0
- i7 K6 P' Q* w( L3 h  b- s9 t2 i( v
C = [] ( _. f% Y1 q( Q" h2 o4 {+ B

1 L" @9 q, f# h. u+ R- MC.append(c)
$ `7 G" Y8 ?8 e1 b; c3 q4 k  z$ q4 D3 B
for d in range(1100)
. z- v( x  n1 N' a' v+ }, Y
. C$ d4 c" i4 q% F    for i in n
7 }4 h8 v8 w$ N* i0 n+ L2 u0 @! i9 }  Z* [
        for j in range(22)
5 D: A/ p, s. _, k2 c! B3 a5 p4 B. x5 O. j
            cj = c[j + 1] - c
' r! P2 s: O3 @8 k+ Z' ]0 C
3 f6 Q- h4 c$ x1 Q4 H* q& n" Z            ci += k1[j + 1]  math.sin(cj) 9 k6 [/ \* {) {* C" ]' K" M
7 r; R' e! U" @/ S
        cii = ci + w 7 @' ?* L1 N0 c3 r* r
, S8 U% G( j  @1 R% l
        h = [0.01  cii + z for z in c]
9 V- ~# {6 X1 ?1 [( M* V: K4 K/ e: y6 Q5 V* ^# s8 ?1 t! P
        C.append(h) 6 V5 |5 D" b3 N) C

3 i" ]- o7 a. p        c = h
/ G4 o8 _0 M- O1 e) i% ]. T1 B# }* U% l3 u) m
def r_function(u)
6 G, X& Q+ w( w, D: [. P, _4 i' ~1 t5 A1 O3 M+ f* A# U
    y1 = 0
) ~: s1 j6 t- l5 q* ?5 b7 v+ [  \  \2 X) s8 E% b
    y2 = 0 $ ~9 S- l; f0 l/ Q' H5 {3 G6 R
2 V! ?. P1 l8 h, ]0 ?
    y12 = 0
$ _5 l: H5 E. Q; m2 c, p2 R
6 r' G1 o' p' ], K# `    y22 = 0
7 E+ f! H' g0 e1 D+ a" `3 \, J2 p4 j8 C! W- `
    r = [] " _9 Y8 b$ Z, P# G2 J
3 N2 B) E3 \0 e2 s: a
    for x in range(1000) " v" t1 \! Z" B; M* l$ j( s" @

  b' a2 n2 U7 y, X+ K5 \9 l        for ul in u
8 _- @: ~, d- P9 p' Q& V# ?$ S" p6 v" j1 `- x
            y1 += math.cos(C[x][ul])
/ N9 }+ c$ ^9 @% ~. @: O1 B0 T0 r% S2 }' {  `$ _/ k+ ~% b
            y2 += math.sin(C[x][ul])
1 Y5 E) X6 u1 y, s8 ]* W. }) I
' {  {7 w9 H1 k        y12 = y1  2 4 I: k( [+ u! `9 W* d7 R
& j4 h3 a8 q% @$ A" W  L
        y22 = y2  2
/ a- V6 F: E( J7 z. \- L2 r, a- l
, e, s% G8 g( q+ ?( e5 N, w. s        r.append((float(1)  N )  ((y12 + y22)  0.5)) 8 U' C: x# M3 D6 ^: @4 G. ~

% b2 O) o' O# {    return r
2 R& N( V8 p* d" a' p6 y+ }. z0 Q$ x/ ?: M2 F& e4 M7 z
r1 = r_function([1,2,3,4,5,6])
% l5 c8 P0 G4 i4 G! [+ v5 h" n1 a1 A+ f+ F! ]
r2 = r_function([7,8,9,10,11,12])
3 M: N3 e. [; X* `. X7 E6 I2 @
3 `( n. ^' u; c2 ^+ H8 Wr3 = r_function([13,14,15,16,17,18]) " x: _. g" V5 X* P1 @! ?- P4 @

0 X/ t7 O8 c. @0 r2 b* nr4 = r_function([19,20,21,22,23,24])
* W* q$ d0 Z- a* W5 b! G' S% N- G. Y
r5 = r_function([25,26,27,28,29,30])
0 d5 E  X: d6 ?; ]
: F3 g) S4 o0 u3 g6 {& B#r6 = r_function([31]) - T' q3 [/ j3 A: c$ d

7 ?6 @7 ]5 A, Y( x6 g3 k3 N! I+ I4 x) G9 C6 ], J

8 k9 h  x5 \- k5 o) _9 Rax = subplot(111) #注意一般都在ax中设置,不再plot中设置
. O' @7 w* ?: R9 a7 T% v6 P4 X( `
ymajorLocator  = MultipleLocator(0.1) #将y轴主刻度标签设置为0.5的倍数
/ b$ q) |) ^+ C4 B
  K& {* z) z6 lax.yaxis.set_major_locator(ymajorLocator) 4 }8 t( X& y. j9 _; i$ t; U
: D/ D& X' x% t- S
plt.plot(t, r1, marker='o', color='green', label='1')
( U5 x: y) M, @9 `3 ?
6 s; x+ {  K8 o  ]+ Gplt.plot(t, r2, marker='D',color='red', label='2')
; M, x8 T. O) d& Y/ ^5 r4 S5 D3 p
plt.plot(t, r3, marker='+',color='skyblue', label='3')
8 U  Z$ m8 d8 ~; ~! Q+ o& R9 U. g8 K9 m
plt.plot(t, r4, marker='h', color='blue', label='4')
& Z4 w9 Z3 X4 M( v# B! y! c
: `( [6 }- u3 y9 L5 zplt.plot(t, r5, marker='_',color='yellow', label='5')
% D, D# S  ?7 l; U  U3 G' j$ M4 D, T5 c4 P3 `1 s$ ^
#plt.plot(t, r6, color='red', label='6')
* A, ?! X, a7 `
2 _$ K+ P1 D+ O2 e( U; J; T( m* U- i3 I0 K# A- Q& J" E5 P; I* g

! r; h3 w) t5 a1 f/ U$ n$ Q( z& ]plt.legend() # 显示图例
  v/ v5 g1 i; D; M
& o; b) |0 m5 L2 k2 D& J! ^0 Aplt.xlabel('iteration times') # |1 x- {' c5 f, t
% ~, z" J9 @0 G- o' I) M. H# g
plt.ylabel('r') 5 o/ y* G  G; o0 h8 l$ u4 Z
/ N2 n) z+ n% {) z* J! x; @
plt.show() ( r4 \  X* [) B
( |  P4 @# _# S1 Z# V
  ^$ ]) L: m, W9 d2 ~: n0 u& x$ y3 j
独立编组结果如图：
2 L4 V& \, C7 e
) \7 o5 x% W8 A1 \3 k( N, O4 {) j- n) I3 J

0 W8 f: p+ ~  o) z, W- Q. ik1独立编组
1 [1 Q/ `; Y1 [1 C. @" j好的，从图像我们来看看Kuramoto模型在描述这个编组的时候，5组最终稳定，我们说这个团队编组还算科学，但我们改变一下K的值，换成分散编组：8 n$ e+ z/ f$ w( r% @% a1 I# k
' T. m$ {$ A8 [7 M+ e$ Z+ a
k2 = [2 B% a& D4 H+ h6 V% A
' U( n1 E& y& E8 Z2 C' i; 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],
. H5 `3 H5 k$ N4 Z6 E" B
. H" |$ w, @+ O2 i7 p% ^5 @" @[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],
( n9 O; D4 q! @0 F8 c6 J1 G# ^
% H& @) C% @* {; G! ?[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], , }' n, ?) m' ^) \

3 H) n. e2 {1 e7 b[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], 7 M9 d4 @; _7 [1 }" f  ^
% Z+ `. r* m7 `* |
[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], & A" m' Z5 g, O, u: V9 J. x5 Z3 B
1 @& z* f8 y( I! y
[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], 4 Q: Z' k0 Q# P$ q* F/ g! r! z' e

  r3 T! c7 O4 f6 y[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],
! X8 ?8 W4 {+ {9 ]- c+ R1 N3 u! Z* G( }) f- Q% N
[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], 0 k$ F/ ~# X9 c/ v2 S& q) y9 J

, t; ?7 I  }# \0 U$ r- Y* d. \[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],
$ ]8 ^/ ~: i  Z8 |  I
7 h+ n% u  X' e; K$ J9 J: h[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],
+ B2 g) J4 L" G* C$ v. _& d$ ?
- @$ ]7 V1 _/ n9 b' [! W[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], # B" |7 p$ B! Y2 Q/ N# b# c
' B  U5 h1 P) k
[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],
9 ~: R7 X1 p* Z' `3 P# \/ I' ~, ^$ |- h0 D
[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],
; \, W+ L0 R- R& F3 c5 A( L* l! d# g. m  Y( f0 s3 e4 d( T
[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],
8 e" J; c) G5 Z2 y& c# K; x: [5 Q7 Z/ m' O0 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],
* ?) v; F0 A; j* ^* z: t3 M7 {% ]0 `+ {- n! E7 k  f% e
[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],
3 }4 J0 \) f1 l- w) T; w  D9 Y! v. ^" M' p2 k
[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], , L) H# @' y% e1 V; t
/ B: C+ R3 d2 r% a% p6 h
[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],
  P8 U3 t# z0 Q' a5 R4 h' R+ `" @( R$ Y+ i
[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],
5 \# s, m: W9 _4 U* F. n
5 j4 d5 K# H( L0 H7 P[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], 5 ~5 I! y5 A* k6 M, Q% S

9 I* X: C( k6 g8 X7 M7 C1 g4 G[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],
+ A0 J" t0 ]0 I5 M0 u, [
7 T( f2 T5 u& t+ J$ F. F  y[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],
8 w4 b4 U  F" H2 B: g+ o1 u4 S3 {2 v  z$ G, `- I  Y: W
[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], : x  d  X3 p% ^. d' P
2 F0 D1 ?1 G& g0 L/ R. I! _
[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],
! U9 R% [2 x( |7 L6 C' N; U5 r- h" p1 i" r) E3 y) [
[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],
8 M# h% G  J# `6 x
7 z; L3 o5 ]) `6 z$ }- c7 V3 Z[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] ( T* W+ }% g+ ]3 `0 z

; S) b6 C) s% j' j( p9 V7 N
: J: }! G8 y9 }5 q5 G$ G
# F$ r* t% W. O4 `" {k2独立编组% x, q; X1 Y, f- f8 k+ e
 7 K* g" {1 h9 W% }& X

8 P+ \1 X! a5 C% y- j) S2 K! B$ l+ g+ u2 O( D
这两幅图像都在开始阶段大幅波动，而后在一定范围内趋于稳定，那么到底哪个分组模式最符合实际，最能突出编组能力呢？) w7 M; N+ R- r6 ~, _2 {9 n
/ {# f3 ]$ n7 D8 H5 x
这里还有一个公式，来解决这个问题，编组同步能力的量化：& P/ U7 ~8 k: i! q7 Z. y6 R2 y
5 z& ?, ~5 f( `! I* }) u9 X* W, J: z# O* m; N* J
M_{s} 就可以描述某编组的同步效果，r_{s} 是达到稳定状态后序参量的均值，β∈(0,1)是调节因子。我们可以用M_{s} 比较编组内部的好坏。那编组间能力的好坏怎样比较呢？% H( R# ^5 _, h5 a! ^, S
1 M$ ]+ X2 q1 C3 D" E1 O/ D" ?
这个Kuramoto模型同样有所考虑，它有一个描述整个系统编组能力的公式：2 @0 ~/ ~' n2 L
 8 J: J/ \! p& z

* S% t+ ]% @" q$ S, R其中，P是编组的数量，M_{i} 是第i个编组的同步能力，\omega _{i} 是编组在整个系统中的权重，\psi {e} 是各编组平均 相位的均值，\Delta \psi 是各编组平均相位的标准差。具体的计算不是这篇文章的重点，就不在计算M{s} 和M的值来比较上述例子独立编组和分散编组的好坏了，本篇文章主要是讲下Kuramoto模型的解决思路，尤其是上面解决\theta 值的方法可以套用在其他Kuramoto模型中，做一个目标估计绰绰有余的。
' Q4 D% U, j$ G" x" l$ d- o' ~6 ?! C4 P2 s' ?5 H" u, O4 z" I
下面是解决Kuramoto模型常用的MATLAB编程方法，具体思路与上述基本一致，这里不再赘述，K的值我们给另一种编组模式：不完全分散编组模式，也是现在实际上最长用的编组方式，直接上代码：; g, d# V2 m9 I6 L6 v
0 b0 a( U7 N# w' T/ ?2 P: a* m+ p, V
clc;
, C4 C- V" i- N" @4 X0 F8 Y; m+ n3 h( O! p) t  g8 [
clear all;
/ z% ?$ t0 p" }: t" t: g, [1 W# T' ^4 B; g" r, t7 V
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; % k; H; x( X7 @" |; ~) `
. [/ w0 p$ v/ z' [
  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;   / s3 \* B( u6 c- q, z

. C6 L. N4 w) p7 q' g2 M: h( o) j  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;   % r6 i3 I  D& x# X* G

  Q& G6 C  x: Q! m' N  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 W. d4 L6 G2 g( b% }
+ c- r* n" m; o# `) |5 P/ Z
  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;   
+ g8 a, ~# x5 b1 M2 e8 v# N1 r- H% f( Z2 B; W
  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;   Q, u7 E6 W( Y1 Z8 S6 F# K

8 h" _- X3 p! H: A. d9 t5 ]& D: \" B  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;   
( Q& c1 C. O* l" x' g- o9 k* N! [. F1 \1 B
  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;   
6 G% j  h! M0 s3 N/ i* n% P1 \, `4 D2 c/ Z5 F1 N& d
  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;   
+ g. ~. l$ w4 Q! j" t  F) L* t; ]# m6 B6 [+ r: |
  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;   
/ ]7 f" [7 R; c: [7 n. E& u) [) t7 t! G
  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;
0 y1 t3 D+ I1 C
* _% n) V1 H: [2 _4 O6 e1 P" ~" M  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;   + U2 M$ \9 U3 @! M  i

/ {/ g% q* o- A5 l  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;   - h3 a) B* Q+ s% \/ d

8 z1 h' S2 V5 d7 s9 @" c' u4 y  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;   # F& ?* _, k( i' l6 u5 q

! u/ h6 W" n5 o0 N" A  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;   
  n+ x% @/ w5 Z  e8 F5 ~
6 k( g! _% {) I  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; 9 H! [2 e+ A6 ^+ |& b; V8 o

  Z( r' D1 g" k$ }# ^4 N& x. 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;   ; u- G' l% f# Z  i
) \1 T' J' e/ Y$ b1 i4 ~+ R
  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;   
4 q) c) t* a# P6 }4 q5 n0 r  w
- l% O* U4 R: C# r  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;   
8 D4 F7 u2 D, I3 g  y8 B9 V  b# r0 S
  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;   / m4 {4 }8 ]0 f
6 ^2 q6 }  ]0 l* p+ M/ [# N
  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; 0 Q0 e: B3 h: T, Z
, O# |  o3 M& G2 ?1 a
  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;   ' j" H, @, A/ X% r0 [' S

  T) |# u3 O' x& e. C2 U1 ~  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;   
: Y6 ~% W. |4 I8 t6 o. u
, U) h. @; e6 Q3 x  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;   
9 _4 ]0 z$ b" u( m, c5 s, {3 m& d" J) y6 g5 p$ b
  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;   " T  F4 A  `9 H( \
" c! @" y7 L2 x5 H6 ?; p: \( h
  ]     7 b; q2 m8 K  M6 s2 j& v9 h- R+ |: j
  U! U1 j; f3 C( R0 T5 J9 E" E* m
N = 25;
9 @$ _7 [8 o& j: ^/ i  I
7 T# m! g' P% q' V5 A# ?) ?$ CStep= 10000; ; h; Z6 f0 Q6 v8 j

; M3 U3 W0 d5 g* r" t5 tTheta=zeros(Step,N); 8 `" H# p. u6 x) c4 X2 l. `

$ \  j: J5 W1 V: sOmega =zeros(Step,N); - L/ `* a$ F1 |" @

7 g- q! `! p4 I: z, b& nDeltaT=0.01 4 D- L7 ~$ `& c) g7 J0 S& E

) U+ L3 W* L0 n, W) i0 ]. n6 j/ {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];
# t8 t2 i( G, R6 S' A. H: x3 v2 P' f8 r* I+ _% e
Omega(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];
1 ?; f% }" s0 Z) R/ U
/ y8 _; ]3 D. F. ]%  求解微分方程
& m: Y, e0 @% l4 W0 o$ }# n
8 o& m* B- R/ \7 M; \; U: e- @for n = 1:Step - d& @! ~/ x' z& @7 Z+ j7 X0 G

$ i* @3 c' h  r6 C5 _! r+ a1 z    for i = 1:N
4 X  o. H+ f3 j
& s# u+ X7 K% T3 j( |. o            Sigma = 0; " Q, m1 u/ f# g$ j3 Y+ h7 R

) N1 p$ ^2 p' S9 G( `% G. S# w        for j = 1:N
2 V+ n7 b+ X5 {2 r: q" X) K$ ]6 s0 Y, B0 s. e
            %%判断k对同步效果的影响
7 Q( Y' g, ]1 O+ A! x7 b
; _/ {4 A2 ]0 K: y= Sigma + k(i,j) * sin( Theta(n,i) - Theta(n,j) );
, X: ?2 d! b* d: ^* D6 B
# p3 w8 z' H. E( @% N        end . i4 ~9 C; I; \3 S+ f, i+ p
8 i  G# i3 t. E  a/ y  j
            Omega(n+1,i) = Omega(n,i) - Sigma;
0 L, r: H2 l- y6 E6 |, G, }; V" U# J- r+ `
            Theta(n+1,i) = Omega(n+1,i)*DeltaT + Theta(n,i);
( \& ^8 B. G3 P) p, g
# \, r  T4 o5 L" q. ^    end   n. E4 ]  a" C: w+ c

( ?$ S/ @; G2 r" f- L1 N5 Zend
" h: X/ N( l# ^/ F  i6 `5 c9 H% }5 T7 s5 _8 Z9 B, N
%  求解序参量r(t)
+ ~8 g7 F+ R8 S. v1 ?6 F( v% M5 R# o* k# e' x! b% P, m6 h
groupN = 5; ! W7 |" R9 D% V, y3 s# a

7 R  C+ F+ r5 Y1 P! \Step = 1000; ! o$ ~) @# c/ n) Q

" ^. n2 L  m" |. x8 S, m8 {$ dfor i = 1:Step 7 E6 C. C+ k8 O) J

5 J/ M2 |* n- O7 E, ^    sigma1 = 0;
" S# Z) m3 {) \
. e, Z; Q7 O. O- E    sigma2 = 0;   J3 H) P( e5 H7 L( \
6 K0 I+ s% R7 w- O3 y5 w1 B* ~
    sumTheta = 0;
4 V6 e! y" |; v5 W. M; x# G. h- i
    for j =6:10  %表示相应的群组 ) b) b* J% G# s2 d7 _

; b9 I* ^9 g! t1 P: R        sigma1=sigma1+cos(Theta(i,j)); ; z* q- q7 ^( R. L( {, I1 |" N+ O

& F  ?$ \' j' F5 S# Z        sigma2=sigma2+sin(Theta(i,j));
& o* \0 k3 G! f, P- q2 p# U5 d
    end ) S: _7 W' R4 A- e6 {
; E  A0 X+ n) `+ I8 X
    r(i)= sqrt( sigma1^2+sigma2^2 )/groupN;   b2 u* u. A4 g& u
# y$ U& |9 C" r% R5 Y  @
end   c$ O! V: F! i% h

+ R5 O8 b5 A2 Y8 v) G* xx= 0.01：DeltaT：DeltaT*Step; . d7 l4 z" G5 H0 e# @6 Y, c
) N. d+ v. W/ W4 t9 u9 n
plot(x,r); 8 t% W, s! n. v! J3 \) H4 M) W
MATLAB版不完全分散编组结果如下图：
+ F, A: k; l+ o5 Q" B 1 j/ F9 p  n6 r" \
作者: regngfpcb    时间: 2021-2-22 10:42
Kuramoto模型在Python和MATLAB中的简单实现



欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/) Powered by Discuz! X3.2