找回密码
 注册
关于网站域名变更的通知
查看: 725|回复: 1
打印 上一主题 下一主题

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

[复制链接]

该用户从未签到

跳转到指定楼层
1#
发表于 2021-2-22 10:03 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式

EDA365欢迎您登录!

您需要 登录 才可以下载或查看,没有帐号?注册

x
事情是这样的,我最近在研究团队编组及内部模式对发挥团队能力的影响,以及如何正确编组让团队能力发挥实现最大化,别问我为什么研究这个,反正稀里糊涂就研究上了。我发现在描述团队编组间及内部同步能力的时候,人们对Kuramoto模型(藏本模型)作了大量的研究,其中包括模型达到完全相位同步的充分条件、耦合强度对于同步的影响、一定条件下振子的收敛速率等。但具体实现一般都在MATLAB中,且网上代码过于复杂(我运行了一遍一堆报错),这里我使用Python和MATLAB对Kuramoto模型简单模拟。
" y7 @7 ^& _/ Y8 b4 w1 m模拟的话还是一遍举个栗子,边分析边测试效果最好,百度学术上有一篇关于Kuramoto模型的简单论文,我们就用它来实现模拟。4 N) A0 b  H8 m
) w. A3 p6 w; q/ Q6 F5 q
好吧,那我们开始,首先是Python实现:0 p) g/ k4 d/ I+ P" A1 K! ]" N

: C+ d2 l. I8 G9 R3 m, o3 Z% W$ [3 g) W( IN维Kuramoto模型的数学描述如下:
: F) ^* N: }, N* h: h# Z; I
; r7 Y. h+ m$ T- h* S  v. |
6 g4 l# K8 t- U! s3 y; w( }没错,是讨厌的数学公式,没事,它可以改写成这样:( v$ n  _2 D9 W- R, w/ G6 ~  r6 d
4 o; c. O6 E- n% H

; j& S- g% p/ S好像还是有点长,那我们在改写一下:
8 w) M( s' E2 p5 G: i( D% L' C" {
4 Z5 a. D2 ~- ?( N) s
4 H! B* y( s9 g; p8 A" B看着好多了,那我就来说说式子中参数的意义,Kij为耦合矩阵,是为了便于描述不同振子间耦合程度不同的情形。最下面那个式子的r就是我们的目标,反应振子间的相关性,这个相关性就可以描述我们想要的编组内部同步能力。
8 X/ r- D0 L' k5 L4 M% n. A5 s$ o2 V! |/ u
哎呦,这个式子看起来好简单,这里要补充一下知识点:同步能力可不是一下子各组该怎么同步直接确定的了,它是一个从开始到稳定的阶段,也就是说随时间变化,最终反映在各组的同步能力才会确定,那么最后图像是什么样子才算同步能力好呢?5 n0 {$ x# q9 d0 `8 w8 N

. l& M4 D3 I4 [( Y; }同步能力好,是指随着时间的推移,各组的同步能力r逐渐稳定,波动现象消失或固定在某一个小范围内。需要注意的是这和各组r值之间的差距没有关系,我们要的是一个平稳的状态。那怎么办找r和t的关系呢?  A  E8 {1 h0 X& A1 Q  k- X- Q9 m

8 g0 k% q7 h- t3 N! K# Q注意看最上面那两个式子,相位(第一项,等号左边那个)上面有个点,这样他可就不简简单单是个相位了,它代表的是相位的变化值,是一个微小的微分值,好吧具体意思就是,那个式子左边展开之后是这样的:
  f6 D/ W& @1 p9 H: j$ H( B 3 m, a7 ^. R* x* n9 v

$ o" u; I$ p6 y' ^' \( ^哎呀,t出现了,其实\theta 与t有关,这里你可能有点绕,因为它们之间的关系是一一对应的,就是说每个时间的t对应了一个\theta ,我下面带入具体数值的时候你就知道了。& `: X2 x' C% o6 ~+ @

, K0 F7 U: B  `% e. p" c4 n组间同步能力与时间t的关系出现了!
3 j4 B* {# i/ A2 |: f( F, n7 a; h- j0 ~5 }. t3 e7 R
也就是说我先用上面的那个公式4计算出来\theta 的值,在带入到公式5,那么t-r关系就可以明确下来了,那现在我们再回过头来看看文章中已经给的例子,看看还有没有未知量。
* r' q2 C: y0 q7 k, p' X- e( X( ]; s7 I$ D: B3 i% F& j
栗子是这样的栗子:
% L: v7 G* g/ X2 n6 l, Q$ [! |) {6 B4 o$ r
假设某机构内部有 4 个编组,每个编组包括 5 个节点(其中 1 个节点为领导节点) 。另外,将上级领导作为一个独立的编组,且只包含一个节点。假设在领导机关增加4名信息传递人员。当以独立编组模式编组时,指定1名信息传递人员为指挥者,其指挥关系与其他编组一致; 当分散编组时,信息传递力量节点的关系与所在编组其他节点指挥关系一 致。其中,完全分散编组模式时,各信息传递力量节点之间无信息共享通道; 不完全分散编组时,在各信息传递人员节点之间建立一条信息共享通道。各编组模式及其拓扑结构如下图所示:
; a% b: |( i4 _' D# h
% m  i4 I% n0 u, ~$ w独立编组模式0 o7 {8 l& L% w" B# u8 k; a
% n/ l3 d# F& k# @8 T

+ x2 t4 d, s# m9 o3 w8 X3 E9 L完全分散编组
7 o8 _# h7 R3 P+ L8 y( q; J
+ N! a; K' _! p, A& ~4 T
. d) m, b; n5 ?; [, j1 A) Z1 v) ?* h
不完全分散编组
; J. X  v9 @/ J6 X
: G1 H& X! @$ a$ L9 G9 t- I, T* D3 J7 T
7 j4 U3 G$ r; I+ ?% H
参数数据
4 x9 n1 S6 E% R$ [2 B$ v7 [  f" Q& D; R3 [6 K

; @% S& v* m( j5 Z: t9 f* I; F7 G+ J! i/ P& `, e6 j
9 ?% S! w: I6 N2 v% T  N3 `
参数确定一下有没有未知量:
2 S# b& G9 M0 g0 h& d) Y1 m2 W9 g) b" Z8 S( F4 U9 c
首先N,数据数目已知,这个有了。6 i6 ?8 }6 d& p( G, y$ U% u! Q3 D

0 y1 w* X( T; H/ j! z# e7 y) uK值是分组内的连接强度,这个是看实际情况,由甲方提供或者自己看着给的,这里就是甲方给的编组图,i与j点的链接强度一目了然,这个有了。) o  |# H$ m" ~' j, f# R- M& ~
' l+ O9 j. B6 r* s
\omega _{i} 是振子i的固有频率,也称自然频率,甲方会给,没法自己估计,这个有了。
( ~1 S/ s8 Y. N3 k( f2 e$ @: h3 I& k2 o% Q3 c
\theta ,\theta 怎么办,初始的\theta 会给,自己也能测的出来,但那么多\theta {j} -\theta{i} 得多少不知道啊,这里通过翻看文章,我发现其实文章是有一个特殊条件的,不然的话是需要研究耦合因子求三种约束条件解情况的,特殊条件就在这:
/ x8 M( d* n2 C; r
$ Q7 X8 }( s& Z/ Y. h假设编组内节点的初始相位差为π/2,且编号最小者为0,随编号增大而增大。
# y& J7 Q$ @! [7 z4 I& r+ w
: T8 A1 O* R+ p) R% J! O0 ?4 P哦,初始相位差知道了,你还告诉了我各个初始相位,那么\theta _{j} -\theta _{i} 的值就在一个范围内的几个固定值里面啊!
+ w+ }( Z  a  Y' M( U/ T1 h2 Q* h- _2 f2 e1 ~
好的,没有未知量了,就是找K的时候麻烦点,没办法,这个决定了编组的不同,写脚本算一下吧:
( E7 B$ _( [& `, ^- y3 m) C4 X, u; X! U( i) D
#codingutf-8
8 v& N/ g# C: ~! }
5 Z2 @6 ?/ ~- Z##ScriptName:KuramotoSimulation.py 8 @/ Y: X( ?% G9 y# f, H. c) g

( |2 d& `. E) h# M& N7 wimport matplotlib.pyplot as plt
4 c) u0 X3 u/ Y7 n
) M9 G1 O0 `5 R5 u# s9 q- p! o0 zfrom pylab import " \2 n9 a. p; ?

$ h% s' O  u% i8 u, l! ]4 t9 A& `from sympy import   @5 U  C6 W0 y# l
- P! i4 |* C8 L7 l1 f# f# U
from matplotlib.ticker import MultipleLocator, FormatStRFormatter # ^  O0 \) N/ N" M
7 w2 g, @% _6 ?: B" g# ]# j
import math
& O' }( Y3 k  N) c7 J8 S) Z7 D) B4 i& C4 e, [
import numpy as np , a% n, i8 y+ w& u0 _* c  m
2 v- q" {2 l7 f/ M
N = 31      #总节点数
9 M$ x* z/ o8 S: o1 m( M/ W8 J. K! D4 F
c=[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]
- Z7 g0 q! }' \0 z4 ^: z+ z2 W5 l; r7 O) l% N0 q2 e6 J8 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]
0 D& @0 B' w( K3 \2 ~% A$ V6 }) J: _- S5 }
k1 = [
! ?! F+ ?3 f) Z( G4 C
. P' C* O; C# E* _& 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,0,0,0,0,0,0],
3 O6 |* v# I+ |8 J+ @6 ~' E1 v5 j$ I8 u! D" @5 R
[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],
1 g; S" B+ o5 J& a! v& G: O" i* j
* z  [4 `0 \) k* G# R, K[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], 3 Z& w6 @3 G! z& T; ^5 C$ }

8 b0 J8 F7 [8 t4 H+ m[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],
0 b7 L5 n: {/ g7 R$ ?" J
# x+ ~$ e0 @; y0 s9 k[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],
# z7 t' u+ ?3 U: A2 p# k& r3 t3 b" Y- N. I1 h- i: Z
[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], $ a$ W* j& |9 l; @
) H4 p! q& P3 C1 s0 T2 w; x6 K: j
[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], 5 ?$ p' m5 @' G  G  z$ N# N4 G

" N1 e. _* q1 O4 X+ }[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], " k" d$ q: C9 F9 M# K- D
/ m) w* }, I$ J. |5 _1 o( D. @
[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], 7 c! u6 K. I, v& o

4 a2 ?7 P+ b7 }5 h8 w; Y[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], + G) Z+ n1 D& o% {( W

5 i( v5 a1 n9 _, @, p: V[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],
4 W, _0 M3 l2 q( Y/ E2 ?; A$ ?& o
[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],   y6 d& o# S# B7 `$ x  _& Z  e
" r, u! W, G6 N4 U* X. S
[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], * J* X  N, H( w. k' \
0 E4 {, [. I8 A' K/ H
[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], # z+ D1 C( h0 M# e

1 u$ p; y4 e! L[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 J( o' Q! R; p6 {; x# o

7 X2 J6 F8 C& c& t& E% g[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],
# g7 t- e6 i3 ^# N4 q4 j& e8 i0 O( G1 x7 \- a
[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],
, g; ]) a( E  v. {3 \6 h, E, u# j7 n3 s6 ^
[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],
: C. o4 i, Q6 M: T
4 C! t+ m2 G: Z3 ^6 d[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],
- p4 n3 I: @2 T$ ?) {  ^4 b& a  U& B9 U0 s# i
[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],
& M( ^5 w/ c& G3 W- y; m& l! ?* g9 G- |* ^
[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], ( I( f) h3 }$ E8 G# _. J

* Q4 p/ P$ X( f' k7 y[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],
/ T  m8 b; Z* [+ f
8 K* [+ V( e: G+ P7 b8 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], 0 i/ S' w2 d5 C9 E

1 @* f& [- S& T. d, p2 @[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],
' s  o* X5 b$ v- ^% Z
5 b8 v+ I7 e0 z[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],
$ R- y/ M! T3 {# S2 x' d0 D2 ^1 R0 @7 K4 J, |
[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],
! {9 K; m1 v! z" {! Y% s) t4 X
5 n" W1 P' n  c* u1 p! M7 c9 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,1,0,0,0,0,0], $ L& ~+ _" a' m3 J
0 y" c" d8 Y7 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,1,0,0,0,0],
5 k7 m) g3 e% l% k0 u8 S# E6 @. c) W* W* e2 z0 v& u% s, A# x; 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,0,0,1,0,0,0], 6 `5 Q( H& z* p9 }

6 P3 g/ h) n# d0 U2 u# q! |# M- 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,0,0,0,1,0,0], 1 V& U& I1 g) g& {3 F- F
# X9 S, Z- S3 I  _4 t* J/ D. S: O) B
[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], ( _9 W5 b0 I* g- L, J( Y

! L% T3 }9 N! z, s  Z( u[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]]   
9 ^" R6 v8 c' H6 h& c5 {, m- r9 ?
) G3 R  E' V- an = [i + 1 for i in range(22)]                #目标划分,24个值,1-24 3 {: L8 J# W. q; U
. Z1 v1 }6 `' |9 K3 S
t = [j  for j in range(1000)]
6 s+ Z+ g# j' ~- d% a7 f5 a3 F" T1 I" ]1 _/ q5 \$ D
ci = 0
: A6 Z4 A2 c# u- Y8 S$ L+ E' Y/ K$ U% o* ^6 t
C = [] 3 F1 l1 F& N: d4 q" F& G
( h% q  F+ p% y2 m/ K$ X/ E% w7 c
C.append(c)
8 i+ x( p8 R3 U) t! d
, C! @' E* |" Sfor d in range(1100)
" y1 U7 E6 U( v0 K/ L- i
0 F. @3 M2 k; l    for i in n
# h; h( F* u. f+ F8 d( e  G: L; {: Y6 S# a6 {8 [/ C
        for j in range(22)
  C6 C: E4 c, f& I8 q
6 E& d3 y% f# g8 @+ F6 D            cj = c[j + 1] - c
2 B6 k  l6 T( k: Z( K0 a* O0 N1 F  }. Z
            ci += k1[j + 1]  math.sin(cj)
2 E: U( u  j. U$ B
7 J; \8 C0 x. B  _" X- H8 }        cii = ci + w # K. I; W5 u5 t& i: }+ E/ D

7 i1 c" E% ]& r) s        h = [0.01  cii + z for z in c] 9 p* @9 p' D* e/ I% |5 ~

7 i& B9 f1 b9 a" W9 U5 g        C.append(h) ! Y/ g% N5 q  _; _9 S9 l' S
' I* p- ?. x5 z$ a8 k3 J
        c = h 3 a2 L; J% {* ^0 Z" @) J

8 H4 k" u9 \  L5 J  _8 {def r_function(u)
! Y# Z- b: Q7 ^- U( E$ q. k3 L) p- E+ G7 u
    y1 = 0
/ v9 I& ]0 M& q2 }& V
  v+ s2 A& a- {* i/ I    y2 = 0
/ T- `, K2 ]! b7 v3 S5 D+ @; Z6 ~$ t9 T- _: `- `! @
    y12 = 0
& I8 L4 V' {0 h; ~, G$ o/ i0 @
0 |+ p: I" V8 u& m1 q3 q% p  e    y22 = 0
  m$ W8 }3 X# q) {* [) A8 i8 ~) K0 i: \
    r = []
: n" q6 Y" N) D5 f) b4 `0 l. t9 A( n1 a) O  j7 n7 q
    for x in range(1000)
' c( b- H% m, d7 G7 I  [2 h( _% C& d8 l7 ^+ p3 D9 ^
        for ul in u
4 P6 s! |1 x) ~) c
0 g- f) L$ C. u            y1 += math.cos(C[x][ul])
4 v+ ^5 U! V  [- A
3 ]6 `$ V" [' L: W3 w            y2 += math.sin(C[x][ul])
8 H/ c$ x6 M$ r  H  M- _- \# M
' L" `" ^: h1 @6 N) u        y12 = y1  2
7 e2 W- T3 I. v7 M
6 T. R0 ], a2 x' E        y22 = y2  2
$ D1 n0 U, U6 K6 G$ G1 [& A) ]$ x
        r.append((float(1)  N )  ((y12 + y22)  0.5)) + c' n) E( r1 q1 L) r

" M" ~2 f0 W% b, X  o7 E    return r
+ ?. A$ ?) j# ]$ r! n3 W" z% [1 g# q6 K1 n$ p
r1 = r_function([1,2,3,4,5,6]) / Q4 y4 u4 q, f1 t  \) {' E
5 ^% X% c3 H1 ^  e& F$ a
r2 = r_function([7,8,9,10,11,12])
) H$ C9 i' {- o0 d2 S. k3 ]9 P) l/ P' y- }- Q* K  N- B4 }
r3 = r_function([13,14,15,16,17,18]) ' d/ H2 u  T6 e$ O5 u
3 F& \$ T) l1 g3 }3 ~" C
r4 = r_function([19,20,21,22,23,24]) 9 a& ~( B# V; m3 H) X  T
1 C1 B! x, G( t+ i/ ]  {! Z2 {
r5 = r_function([25,26,27,28,29,30]) 3 Y  G! b* n, C+ P7 K( l
: l5 h5 M+ c8 Y2 M4 y( ]9 U
#r6 = r_function([31])
; E' k, K1 M) ~  Y& {) K1 O! w* u; G& @6 g3 Q2 }5 X

# {+ h; a' {' k$ V$ I& Y9 l0 e  C% }* c6 `2 o* c
ax = subplot(111) #注意一般都在ax中设置,不再plot中设置 , P. O/ e/ t9 Z

9 e! S. r1 @4 j( y7 EymajorLocator  = MultipleLocator(0.1) #将y轴主刻度标签设置为0.5的倍数 6 F6 I7 d8 {$ \0 ~( ^! t
: k& H+ M6 ^: ]  U1 f# n$ R
ax.yaxis.set_major_locator(ymajorLocator) 6 Q3 t$ h" ?0 R, h$ r, u

) N, C+ v* `# @5 j7 G' Z& Yplt.plot(t, r1, marker='o', color='green', label='1')
4 R2 K* h0 @) t+ K1 d8 y" T1 ^/ t, d& H# C5 n5 Q' D  R# D+ d
plt.plot(t, r2, marker='D',color='red', label='2')
( L2 k, u8 b# f3 h1 m  |- ?- _3 V) @( f3 E/ P/ q7 ~* Y
plt.plot(t, r3, marker='+',color='skyblue', label='3')
7 L- O  ^0 O0 k" Y$ D0 H) ~% N9 C7 F4 \
plt.plot(t, r4, marker='h', color='blue', label='4') ) h5 k7 O3 }% l8 q0 Y* @
* A1 f  C6 S3 ^- Y5 K. Y
plt.plot(t, r5, marker='_',color='yellow', label='5') 1 U% `# M, R9 b4 z0 q$ T

. a7 |' j. z/ W1 x  V#plt.plot(t, r6, color='red', label='6')
& v% i  n7 W. J) y$ S$ z" k) Q( Y" U  X* S
/ h. @; Z7 S+ r9 ^# ?
# t) c% q& n, k
plt.legend() # 显示图例
6 ^8 s/ ?$ v' l  W1 s4 V4 k5 |
! x. P. x1 s2 cplt.xlabel('iteration times')
2 M( R+ M- }9 B
* P' S# q( Z5 o4 |5 rplt.ylabel('r')
  P, ~1 W# |, v2 W1 c: K% g" N! i. D$ S) W7 j# t$ c
plt.show()
& a+ c( Q5 ]1 o, H2 M- H9 o
0 h# F$ d3 j, y5 ^0 U% ^0 h8 y9 A; x7 Q
独立编组结果如图:
1 V; h( j: r" N

7 v! O0 }3 u% g0 J- a& ^& u& m5 R6 J; O

4 ]9 C1 `+ |1 b# e- b9 \$ `( Ak1独立编组
8 ]( z0 Y! n, B  Y好的,从图像我们来看看Kuramoto模型在描述这个编组的时候,5组最终稳定,我们说这个团队编组还算科学,但我们改变一下K的值,换成分散编组:9 T' w2 {  O  U% K& c# z

% O3 j' D* ~. `& y% R1 tk2 = [
$ J6 l+ _: X, b! X# _+ P
  k1 v7 i# n8 c  {) _[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],   q9 ?) V% S% D+ w2 D3 h' M

0 v. }; o7 L- p5 Q; p. n[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], ! p' z6 K  ~2 \3 Y, E
# u3 _3 A- H: L1 z! n
[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],
3 R( P( `3 O8 F, v' D7 F9 {
! A8 t5 A7 J; ]3 O7 C! ~) U# \[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],
8 `; w; ]4 {  G4 A
" g# G6 w, L9 g% |/ x9 x[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], # B3 M# S" v4 N) X4 B0 d

9 h' L8 C+ b  N# H) _) Q8 S, M[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],
) J5 u. G# L! d7 ?( v* E. P( p( K. f
[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], " c) K1 x7 ^" G: L& r

0 x" t4 F) }% w2 `# Y0 {[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], 2 E& j6 \! ~4 A* [
$ E3 V1 ]! y( Q* `; N
[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],
, C, Z7 Z0 p0 A
/ S) p, a7 z5 K* T# L/ ^# A3 R. s[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], . C7 l! w* w, [' z0 Q: e/ [5 j, X

7 a/ E+ Q! F1 a2 M0 h3 z$ _[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],
2 s2 V4 `# N; h+ ~& F
9 e9 G- U2 B! S5 G; j+ a+ b! L[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],
# v( A+ w" B1 \  A
0 [( W5 n2 ]) c[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], 0 Q+ X) E) b0 N8 d

3 h+ D. {; n  i5 Y3 V! V[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], * i! e. e' M# M% V9 X  e
8 s; ^6 I% D4 R
[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], 2 C& Z) w0 ^: A$ [1 T
# i) ?2 f* F( J
[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], # ]( Z% ]6 T) H& o

, A& j' O  b& g  q9 N. E$ O[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],
! Q& d4 |2 ?* _9 p7 f7 \4 l' W1 r1 L( E. j) M. G% L. I0 h& l
[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],
$ Q) X- P# H! |0 [* P" d3 W  M7 Y4 r* |2 A0 d( X  l7 }- j8 t: s4 n
[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 f2 R" [2 s1 H6 B; X* F- }
# H& {8 T$ h- W5 ?
[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],
' H( P2 V) w; G2 q
  U: D, l1 g4 b) |- D1 i6 y[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], # {6 s5 n  j( G; J
# h2 @, @! s6 q7 a
[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 _- x$ b5 b2 V7 i) g2 f* Q7 W/ g: h9 r
[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],
( A. K: x( d0 T% `1 e
& l+ w7 P; P# C1 c1 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],
" b: ]; W8 Z1 F, ?  \( H0 Z
; \; c1 b& r0 v- V2 w% o+ }[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], + ?6 C8 X+ P  j$ _7 z( Z
6 a0 e+ y2 J( K0 M2 C+ x' k
[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]
& [$ }1 x& _6 V. N1 n) O
6 H# J: h6 F* o8 E1 q: Y0 a
/ c9 C  [# a7 Y7 t+ d! o$ Z5 T
3 n8 W" I1 `, q) s
k2独立编组
6 A7 A! Y: g3 e$ Z
: T; v/ C3 o/ D0 Y) g

6 a* G0 D5 J' L$ x8 B# S
* n3 p0 }7 _# {  |这两幅图像都在开始阶段大幅波动,而后在一定范围内趋于稳定,那么到底哪个分组模式最符合实际,最能突出编组能力呢?6 L5 o9 I1 z' h" h+ ]* I) S7 o
8 @; A( @  w" ]# }( [
这里还有一个公式,来解决这个问题,编组同步能力的量化:8 C$ @$ ?( D$ X5 h0 J
# n& f5 t" |# Z4 Z4 |
( |7 s. S8 @2 m5 Y
M_{s} 就可以描述某编组的同步效果,r_{s} 是达到稳定状态后序参量的均值,β∈(0,1)是调节因子。我们可以用M_{s} 比较编组内部的好坏。那编组间能力的好坏怎样比较呢?5 c' H4 D# n* G/ o/ o
4 y) n6 ^5 |2 k) p) c
这个Kuramoto模型同样有所考虑,它有一个描述整个系统编组能力的公式:
% Q, J/ Z/ n2 I7 ]# |0 s

: S3 G. y' G2 r* P' `2 P% u% c( Z! L
4 ^1 H/ L0 K( J5 V! N% ^其中,P是编组的数量,M_{i} 是第i个编组的同步能力,\omega _{i} 是编组在整个系统中的权重,\psi {e} 是各编组平均 相位的均值,\Delta \psi 是各编组平均相位的标准差。具体的计算不是这篇文章的重点,就不在计算M{s} 和M的值来比较上述例子独立编组和分散编组的好坏了,本篇文章主要是讲下Kuramoto模型的解决思路,尤其是上面解决\theta 值的方法可以套用在其他Kuramoto模型中,做一个目标估计绰绰有余的。
9 s9 N$ @( Q8 b/ X
/ {  B/ K: {2 N% z4 @5 E% ]下面是解决Kuramoto模型常用的MATLAB编程方法,具体思路与上述基本一致,这里不再赘述,K的值我们给另一种编组模式:不完全分散编组模式,也是现在实际上最长用的编组方式,直接上代码:
( M5 v  h+ H4 c3 J2 \& p) J  h9 n' v5 `) F
clc;
9 W0 Q( ~" f+ ~% M; m, N8 q5 f4 u/ Q' V
clear all; ( E+ y& R+ v7 w
. W% w6 G% K+ @0 ~) p* p! \
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;
$ D) v  r6 H5 s' r3 a/ r9 W3 q( F1 g
  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;   " p/ o4 w. l2 f9 t+ k* |
: M9 D0 W- z+ L8 {! M" v
  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 G% i. {6 A6 B' q( y. n

# X9 e; g! C# 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;   7 _! P# ~' M" k

! r+ a1 G. A& ~. ~+ l  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;   
! l, t' K; y# A! Q# }9 V6 C) w& v* Q1 r
  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; 0 {+ w$ f; i# h  }
3 \8 W5 ]" ^0 D' @5 I% ~
  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;   . {# t3 a# p) `8 `+ q6 b& W& M

  Z$ @5 \- A& y. G  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;   
; I3 }- `; M, q. q7 b0 ^* l' e
6 G5 j$ L8 q( g  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;   
+ z0 h+ N' E# B: m5 b& L+ }$ U  J  y/ q
  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;   9 c* I" Y  J: a+ ~, Y
  F* f- ~1 l: Q
  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;
9 s2 S/ d; M! ~' e/ J, |$ t  ]$ M& F3 n7 }
  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;   
* J8 h8 N8 P. H3 l* ?3 J5 [8 o% r2 ]' A2 ~0 Z2 X
  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;   
: p7 D$ F7 [6 M: m# [
# J7 |  w! y7 O+ 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;   
% v# W' m0 o# n1 h; ^" E' [/ c
' q9 P0 ]0 N7 S! [  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;   / D  d6 h- t2 {9 O+ p6 G0 u
( X: O, f! O* H) z0 q* g$ \
  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;
0 ?3 y. g" r% U7 J- [6 M4 k! k
4 s5 ^4 Y( m& Q8 j2 c  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 l+ J# H/ Y) l; g6 K

9 d, G* M1 M0 w  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;     T) l# u2 B" d6 j$ K- l5 A

8 C! @  H0 T/ B$ U$ 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;   , f' g9 b  R: o( S* Y* O

, O. a3 w' q: c1 n# P" 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;   
& p  r  ]2 s2 ^9 u- i  J* G9 x" {+ a( b7 M1 D. Y/ L0 D$ K6 c
  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; + ~( z; X' ^! `! P, K) g

% z, L/ l$ a0 i, g  ^3 ~  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;   
6 I, ~6 n; g7 K) u, w- E: }0 L  }, f8 v. M) c' I
  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;   ; a  W4 r3 Z4 n) x% d
& f: U% ?0 X* D
  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;   
# Z; U9 P- ~  z1 N5 P1 D
+ s$ a% O+ T) `# Y6 C' K3 \2 H- Y  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;   
7 |  j# T. q( u+ S7 }7 D& w: }9 U
  ]     , Y7 s5 J# s- w, [+ j3 o% u
- P3 s6 _  k8 k9 j( }
N = 25;
7 `% |" B$ ]/ P+ M8 T! Z. Y! K' R  ~' ?
$ b% K& }  _3 l. C8 KStep= 10000; 7 F& V* b8 z1 O

3 B! f3 v! A) ^, s1 ]9 bTheta=zeros(Step,N);
/ j* f% t: R3 I2 d3 ?& G+ S: F- c- v# @
Omega =zeros(Step,N); , {( Q' d) E: F# T: n( r

2 B. Q0 |1 F6 K2 a8 nDeltaT=0.01 " r$ Z3 v( B( `( i) H+ |4 Q  R

  N: e* x0 U$ M2 W* v$ }% x: nTheta(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];
5 S  P- R4 ~# O  n5 q% v
5 X7 q4 M5 w  U9 p" J- XOmega(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]; . O9 X' |4 P, D' g4 t- c$ Y

0 U; v. C# b/ L7 [! w, x( f- c%  求解微分方程 2 ^( n& D+ y! u1 D( d
9 }2 M- \2 w( N# L' N, ^. Q! a
for n = 1:Step 0 l1 }  f/ @$ @3 b7 }1 P
) W. {+ b) u0 C4 x; M+ w0 u
    for i = 1:N 0 P+ C$ L6 x6 Y2 W2 y
3 X5 Y8 ^& B3 V
            Sigma = 0;
3 Q) M( b2 ~' x+ _: g1 V4 n: X, q6 I6 T+ y
        for j = 1:N
# d/ |( M  s9 C; e* n" k: ?* y" C
" f4 U- e% ]% l% V            %%判断k对同步效果的影响
( v3 B; P7 E' l# u, W; H; S
) I% q& `. l  S! z= Sigma + k(i,j) * sin( Theta(n,i) - Theta(n,j) );
6 c- h4 t7 g% w& a) a# R, `! S! ]! g( T
        end
3 B" {& Z! S8 ?: i4 F) M& f5 F# G4 v$ k: @
            Omega(n+1,i) = Omega(n,i) - Sigma;
; T9 R; g! E/ B. v2 I, {, J. G. N) e# Y9 ~  B8 t% b
            Theta(n+1,i) = Omega(n+1,i)*DeltaT + Theta(n,i); 7 T4 v/ C. _5 z8 Z
  P1 ]5 S, {2 P; T; _+ a
    end   L8 s9 L: D1 G; N1 T2 o
' M5 a. e/ V% ]0 M  e
end
, X- {! i0 ~3 R. \( n5 z% Y, k$ ]7 E# W  P4 H: {
%  求解序参量r(t)
- Y; U9 f/ n, G
; ~6 |4 f( T8 G! v3 UgroupN = 5; " e* }6 r) C  e) A1 x9 y

3 g( T. Z% k: D# V' \Step = 1000; 1 g2 w5 i& T( `& F

- r. ~! C  p, r+ H# x5 dfor i = 1:Step . |! v- q  c5 y- @

* S3 h1 i7 w/ f    sigma1 = 0; # p  @. d* m3 z, B+ i' F

! K& S0 f; L( ~+ ~% {1 j& I- d2 u    sigma2 = 0;
( u: Q& A2 n. g/ s) g: T; {" Y7 ^0 i* |# l( J* z3 n; i
    sumTheta = 0;
  i$ l  u0 A6 r# p1 O0 N; g% K7 A* v8 D) N6 m/ a  \' g# g( ^0 r
    for j =6:10  %表示相应的群组
+ |$ C$ W. {+ _4 W. M6 v, h6 w# w# L8 u9 t! [7 M
        sigma1=sigma1+cos(Theta(i,j)); % B# ]6 `2 _! u

9 m7 V# ^) Z( H4 j; m! c6 ~        sigma2=sigma2+sin(Theta(i,j)); ( E. G$ E3 u: n2 p: k* o$ X5 [. Y; l
3 U: e3 J5 U8 }& x+ K! Z
    end , V. a6 s- E7 C* ^
' ^& o3 _% l& e8 c( `
    r(i)= sqrt( sigma1^2+sigma2^2 )/groupN;
% o' I' }0 x1 q/ `
& D& G3 X& |' ]9 M. Y  Xend ' n: ^" x) v$ f* L) I% f0 Y
, y* [, I( h5 \- J4 I: w1 N* p2 a* V
x= 0.01:DeltaT
:DeltaT*Step;
" q4 H- S' p1 g+ q/ I) q. B$ x2 p  m
plot(x,r); + R. ?  z# D* {! w9 i! d1 b8 P0 A
MATLAB版不完全分散编组结果如下图:5 P* g, @7 q3 g. l

2 a' i8 {* r7 P

该用户从未签到

2#
发表于 2021-2-22 10:42 | 只看该作者
Kuramoto模型在Python和MATLAB中的简单实现
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

推荐内容上一条 /1 下一条

EDA365公众号

关于我们|手机版|EDA365电子论坛网 ( 粤ICP备18020198号-1 )

GMT+8, 2025-11-1 05:20 , Processed in 0.187500 second(s), 26 queries , Gzip On.

深圳市墨知创新科技有限公司

地址:深圳市南山区科技生态园2栋A座805 电话:19926409050

快速回复 返回顶部 返回列表