复杂网络建模与控制研究报告——BA无标度网络的牵制控制

1、复杂网络建模与控制期末研究报告BA无标度网络的牵制控制 姓名： 学号： 专业：摘要本文首先简要介绍了复杂网络的相关知识，其次构造了一个BA无标度网络结构，并给出了网络的具体结构参数；然后以构造的BA无标度网络为研究对象，分析了耦合强度、牵制密度以及牵制强度三个参数对网络稳定性的影响；最后，分析和比较了特定牵制控制和随机牵制控制策略对网络稳定性的影响。20目录摘要1目录21 复杂网络简介31.1 复杂网络的介绍31.2 复杂网络的常见网络模型31.2.1 规则网络模型31.2.2 随机网络模型31.2.3 小世界网络模型41.2.4 无标度网络模型41.3 网络牵制控制52 牵制控制稳定性条件5

2、3 BA无标度网络模型的构造54 不同参数对网络稳定性的影响84.1 对网络施加牵制控制84.1.1 Lorenz系统84.1.2 BA无标度网络的状态方程94.2 耦合强度对网络稳定性的影响94.3 牵制密度对网络稳定性的影响124.4 牵制强度对网络稳定性的影响145 不同控制方法对网络稳定性的影响156 总结191 复杂网络简介1.1 复杂网络的介绍复杂网络是具有复杂的结构和/或具有复杂的节点行为的网络系统。网络系统的复杂性主要体现在：结构复杂性、节点复杂性、结构与节点之间的相互影响、网络之间的相互影响。人们生活在一个充满着各种各样的复杂网络的世界中，例如：生命科学领域的各种网络（如细胞

3、网络、蛋白质蛋白质作用网络、蛋白质折叠网络、神经网络、生态网络）、Internet/WWW网络、社会网络、流行性疾病的传播网络、科学家合作网络、语言学网络，等等。人类社会的网络化是一把双刃剑：它既给人类社会的生产与生活带来了极大的便利，提高了生产效率和生活水准，但也带来了一定的负面冲击，如局部动荡或传染病等更容易向全球扩散。学术界关于复杂网络的研究方兴未艾。特别是国际上有两项开创性工作掀起了一股不小的研究复杂网络的热潮。一是1998年Watts和Strogatz在Nature杂志上发表文章，引入了小世界(Small-World)网络模型，以描述从完全规则网络到完全随机网络的转变。小世界网络既具

4、有与规则网络类似的聚类特性，又具有与随机网络类似的较小的平均路径长度。二是1999年Barabási和Albert在Science上发表文章指出，许多实际的复杂网络的连接度分布具有幂律形式。由于幂律分布没有明显的特征长度，该类网络又被称为无标度(Scale-Free)网络。1.2 复杂网络的常见网络模型1.2.1 规则网络模型比较常见的规则网络有全局耦合网络、最近邻耦合网络以及星形耦合网络。如果一个网络中任意两个节点之间都有边直接相连，那么就称该网络为全局耦合网络；如果一个网络中每个节点只和它周围的邻居节点相连，那么就称该网络为最近邻耦合网络；如果一个网络中有一个中心节点，其余节点都

5、只与这个中心节点连接，而它们彼此之间不连接，则称该网络为星形耦合网络。这三种规则网络如图1-1所示。图1-1 三种规则网络（）1.2.2 随机网络模型最经典的随机网络模型是ER随机图，分为具有固定边数的ER随机图和具有固定连边概率的ER随机图。后者的构造算法如下：1) 初始化：给定个节点以及连边概率。2) 随机连边：1 选择一对没有边相连的不同的节点。2 生成一个随机数。3 如果，那么在这对节点之间添加一条边；否则就不添加边。4 重复步骤，直至所有的节点对都被选择过一次。图1-2 和时所生成的随机图的三个实例图1-2表示的是取节点数和概率时，所生成的ER随机图的三个实例。1.2.3 小世界网络

6、模型作为从完全规则网络向完全随机网络的过度，只要在规则网络中引入少许的随机性就可以产生具有小世界特征的网络模型，现在常称为WS小世界模型。其构造算法如下：1) 从规则图开始：初始有数目固定的个节点，每个节点有个最近邻，构成一个规则的一维圆环。2) 随机化重连：以概率对圆环中的每一条边进行重新连接。这个过程不能自身连接和重复连接。WS小世界网络模型具有明显的聚类和小世界特征，克服了规则网络模型和ER随机网络模型的不足。WS小世界网络模型如图1-3所示。图1-3 WS小世界网络模型1.2.4 无标度网络模型构造BA无标度网络模型算法如下：1) 增长：从一个具有个节点的连通网络开始，每次引入一个新的

7、节点并且连到个已存在的节点上，这里。2) 优先连接：一个新节点与一个已经存在的节点相连接的概率正比于节点的度：经过个时间步之后，BA模型演化成一个具有个节点的网络。图1-4显示了参数为、的BA网络的演化过程。图1-4 BA模型的演化（参数）1.3 网络牵制控制所谓牵制控制，就是对网络中的少部分节点施加控制而使得整个网络达到所期望的行为。其核心思想是网络中小部分节点能够“领导”网络的其他节点逐渐实现整个网络的同步。牵制控制的优势是控制器个数少，计算量小，资源花费少。这种“牵一发而动全身”的思想已在实际众多复杂系统中得到证实。牵制控制分为特定牵制和随机牵制。特定牵制是根据节点的某些具体特性，比如节

8、点的度，有选择地选取部分节点进行控制。随机牵制是在网络中以某一概率随机选择部分节点进行控制。对于BA无标度网络而言，特定牵制控制效果和随机牵制控制效果差异较大，因此，本文主要以BA无标度网络为研究对象，来探讨不同的参数、不同的控制方法对网络稳定性的影响。2 牵制控制稳定性条件设控制网络的状态方程为 （2-1）分析可知，当存在一个常数，使得是Hurwitz稳定矩阵时，只要耦合强度满足下面的条件： （2-2）网络系统（2-1）就可以被牵制控制到平衡点，这里是在平衡点的Jacobian矩阵，是矩阵的最小特征值；矩阵，L为网络系统的Laplacian矩阵，矩阵称为控制增益矩阵，为牵制密度（受控节点个数

9、），矩阵的所有特征值均大于0；矩阵为内耦合矩阵。3 BA无标度网络模型的构造按照1.2.4节的算法构造BA无标度网络，其MATLAB程序代码如下：%m0=3;m=2;N=40;%网络结构参数randnum=rand(2,N)*100;x=randnum(1,:);y=randnum(2,:);plot(x,y,'r.','Markersize',20);%产生N=40个随机点hold on;Adjacent=zeros(N);%初始化邻接矩阵%for i=1:m0%开始三个节点两两相连 for j=i+1:m0 Adjacent(i,j)=1;Adjacent(

10、j,i)=1; endendnum=m0;for i=m0+1:N S=func(Adjacent,num); Q=zeros(1,num); for P=2:num+1 Q(1,P-1)=S(1,P)/S(1,num+1); end for j=1:m%轮盘赌 random_data=rand(1); a=find(Q>=random_data); p=a(1);%利用轮盘赌选出节点p与新的节点相连 Adjacent(i,p)=1;Adjacent(p,i)=1; end num=num+1;endfor i=1:N%根据邻接矩阵连接各个节点 for j=1:N if Adjacent

11、(i,j)=1 plot(x(i),x(j),y(i),y(j),'linewidth',1); hold on; end endendhold off%Degree=zeros(N);%初始化度矩阵for i=1:N Adjacent_du=0; for j=1:N Adjacent_du=Adjacent(i,j)+Adjacent_du; Degree(i,i)=Adjacent_du; endend%Laplacian=Degree-Adjacent;%得到拉氏矩阵save 'E:MyDiraaaa.m' Laplacian -asciiL=log(N)

12、/log(log(N);%网络参数平均路径长度Lc=(log(N-m0)2/(N-m0);%网络参数聚类系数cPP=tabulate(sum(Degree,2);%网络参数度分布通过MATLAB运行得到一个参数为的BA无标度网络结构图，如图3-1所示。MATLAB返回的数据还有Laplacian矩阵、平均路径长度、聚类系数、度分布等网络结构参数。图3-1 MATLAB返回的BA无标度网络结构图平均路径长度聚类系数度分布（第一列为节点的度，第二列为节点个数，第三列为所占百分比）：4 不同参数对网络稳定性的影响4.1 对网络施加牵制控制4.1.1 Lorenz系统对于如下状态方程描述的Lorenz

13、系统： （4-1）当系统参数为时，它是一个混沌吸引子。由所求的为平衡状态，那么，由可知，该系统具有以下三个不稳定平衡点：求Lorenz系统的Jacobian矩阵的MATLAB代码如下：syms x1 x2 x3;b1=0;0;0;b2=6*sqrt(2);6*sqrt(2);27;b3=-6*sqrt(2);-6*sqrt(2);27;%b1,b2,b3为系统的三个不稳定平衡点Jcb=jacobian(10*x2-10*x1;28*x1-x2-x1*x3;x1*x2-8/3*x3,x1 x2 x3);Jcb1=subs(Jcb,'x1',b1(1);%第一个平衡点的Jacobi

14、an矩阵Jcb1=subs(Jcb1,'x2',b1(2);Jcb1=subs(Jcb1,'x3',b1(3);lambda1=eig(Jcb1);%特征根lambda1=double(lambda1);Jcb2=subs(Jcb,'x1',b2(1);%第二个平衡点的Jacobian矩阵Jcb2=subs(Jcb2,'x2',b2(2);Jcb2=subs(Jcb2,'x3',b2(3);lambda2=eig(Jcb2);lambda2=double(lambda2);Jcb3=subs(Jcb,'x1

15、',b3(1);%第三个平衡点的Jacobian矩阵Jcb3=subs(Jcb3,'x2',b3(2);Jcb3=subs(Jcb3,'x3',b3(3);lambda3=eig(Jcb3);lambda3=double(lambda3);求出的结果为（lambda=）：， ， ， 4.1.2 BA无标度网络的状态方程现在假设网络中每个节点都是Lorenz系统。设受控网络的状态方程为 （4-2）每个节点都有3个分量，则受控网络的状态方程为 （4-3）4.2 耦合强度对网络稳定性的影响从拉氏矩阵中可看到度最大的节点为第1个节点（度为14），下面首先以第1个

16、节点作为牵制控制对象，来研究耦合强度c对网络稳定性的影响。为简单起见，取。由Jcb1和lambda1可知，当=11.8277时，即可保证为Hurwitz稳定矩阵。为了研究耦合强度对网络稳定性的影响，分别取固定不变，选择控制目标为。求取的临界值的MATLAB程序代码如下：%求取耦合强度ccN=40;l=1;d=10;D=zeros(N);for i=1:l D(i,i)=d;endB=Laplacian+D;cc=max(lambda1)/min(eig(B);求得的临界值。下面分别令取不同的值来研究耦合强度对网络稳定性的影响。受控网络的MATLAB程序代码如下：function f=zhuan

17、gtaifangcheng(t,x)%N=40;l=1;d=10;cc=50;%sigma=zeros(3*N,1);f=zeros(3*N,1);Laplacian=load ('E:MyDiraaaa.m');balance=0;0;0;%for i=1:l %施加控制节点的动力学方程 for j=1:N sigma(3*(i-1)+1)=Laplacian(i,j)*x(3*(j-1)+1)+sigma(3*(i-1)+1); sigma(3*(i-1)+2)=Laplacian(i,j)*x(3*(j-1)+2)+sigma(3*(i-1)+2); sigma(3*(i

18、-1)+3)=Laplacian(i,j)*x(3*(j-1)+3)+sigma(3*(i-1)+3); end f(3*(i-1)+1)=10*x(3*(i-1)+2)-10*x(3*(i-1)+1)-cc*sigma(3*(i-1)+1)-cc*d*(x(3*(i-1)+1)-balance(1); f(3*(i-1)+2)=28*x(3*(i-1)+1)-x(3*(i-1)+2)-x(3*(i-1)+1)*x(3*(i-1)+3)-cc*sigma(3*(i-1)+2)-cc*d*(x(3*(i-1)+2)-balance(2); f(3*(i-1)+3)=x(3*(i-1)+1)*x(

19、3*(i-1)+2)-(8/3)*x(3*(i-1)+3)-cc*sigma(3*(i-1)+3)-cc*d*(x(3*(i-1)+3)-balance(3);endfor i=l+1:N %其余节点的动力学方程 for j=1:N sigma(3*(i-1)+1)=Laplacian(i,j)*x(3*(j-1)+1)+ sigma(3*(i-1)+1); sigma(3*(i-1)+2)=Laplacian(i,j)*x(3*(j-1)+2)+ sigma(3*(i-1)+2); sigma(3*(i-1)+3)=Laplacian(i,j)*x(3*(j-1)+3)+ sigma(3*(

20、i-1)+3); end f(3*(i-1)+1)=10*x(3*(i-1)+2)-10*x(3*(i-1)+1)-cc*sigma(3*(i-1)+1); f(3*(i-1)+2)=28*x(3*(i-1)+1)-x(3*(i-1)+2)-x(3*(i-1)+1)*x(3*(i-1)+3)-cc*sigma(3*(i-1)+2); f(3*(i-1)+3)=x(3*(i-1)+1)*x(3*(i-1)+2)-(8/3)*x(3*(i-1)+3)-cc*sigma(3*(i-1)+3);endend对受控系统进行求解，并画出相应的状态图，MATLAB程序如下：%N=40;l=1;d=10;cc

21、=50;x0=3*ones(3*N,1); %微分方程的初值t,x=ode45(zhuangtaifangcheng,0 5,x0);%四五阶龙格-库塔解常微分方程hold onfor i=1:3:(3*N-2) plot(t,x(:,i),'r')endfor i=2:3:(3*N-1) plot(t,x(:,i),'b')endfor i=3:3:(3*N) plot(t,x(:,i),'g')endxlabel('t')ylabel('xi(i=1,.,N)')分别取三个值（临界值为96.7822）时，相应的

22、状态图如图4-1到图4-3所示。图4-1 c=50,l=1,d=10控制到的仿真图图4-2 c=100,l=1,d=10控制到的仿真图图4-3 c=110,l=1,d=10控制到的仿真图总结：由图4-1可以看到，当时，无法将网络节点控制到平衡点；图4-2说明当稍大于96.7822时，几乎等于临界值，要较长的时间才能把节点状态控制到平衡点；图4-3说明当时，可以很快将网络所有节点控制到到平衡点。4.3 牵制密度对网络稳定性的影响将和的取值固定，选择控制目标为，分别取不同的值来研究牵制密度对网络稳定性的影响。由于控制目标改变，需要对程序稍加改动。在计算的临界值的程序里需将改为，值改变的时候矩阵也要

23、改变。计算的临界值的程序如下：%N=40;l=1;d=100;D=zeros(N);for i=1:l D(i,i)=d;endB=Laplacian+D;cc=max(real(lambda2)/min(eig(B);取固定不变，MATLAB返回的的临界值为0.4529。取不同值时，MATLAB返回的仿真图如图4-4到图4-6所示。图4-4 c=2,l=1,d=100控制到的仿真图图4-5 c=2,l=5,d=100控制到的仿真图图4-6 c=2,l=10,d=100控制到的仿真图图4-4、图4-5和图4-6可以说明，牵制密度越大，即受控节点个数越多，网络中的节点就越快的收敛到平衡状态。实际

24、上，当时，的临界值为0.4529，当取时，矩阵B的最小特征值已改变，此时的临界值变为0.1531；当取时，的临界值变为0.1052。随着的增大，对的要求越来越低，如果保持不变，则必然会更快到达平衡点。4.4 牵制强度对网络稳定性的影响将和固定取值为，控制目标仍然为。分别取不同的值来研究牵制强度对网络稳定性的影响。当取不同值时的仿真图如图4-7到4-9所示。图4-7 c=5,l=1,d=10控制到的仿真图图4-8 c=5,l=1,d=30控制到的仿真图图4-9 c=5,l=1,d=50控制到的仿真图总结：图4-7、图4-8和图4-9可以说明，牵制强度d越大，网络中节点就越快被控制到平衡点。实际上

25、，当时，的临界值为0.7688，而时，的临界值降低为0.5317，时，c的临界值降低为0.4862。随着的增大，对的要求越来越低，在保持不变的情况下，网络节点必然会更快的收敛到平衡点。5 不同控制方法对网络稳定性的影响度最大的节点是第1个节点，度为14，因此选用节点1作为特定牵制的控制对象，控制目标仍然是，系统参数为。仿真图如图5-1所示。图5-1 c=5,l=1,d=30控制到的仿真图（特定牵制）为了比较随机牵制和特定牵制对网络稳定性的影响，需要进行随机牵制控制的仿真。随机找到一个节点作为控制对象，网络节点的状态方程的MATLAB代码要稍微进行改动，改动后如下：function f=zhua

26、ngtaifangcheng1(t,x)%N=40;l=1;d=30;cc=5;%sigma=zeros(3*N,1);f=zeros(3*N,1);Laplacian=load ('E:MyDiraaaa.m');balance=6*sqrt(2);6*sqrt(2);27;%rr=load ('E:MyDirbbbb.m')for i=1:N if i=rr %施加控制节点的动力学方程 for j=1:N sigma(3*(i-1)+1)=Laplacian(i,j)*x(3*(j-1)+1)+sigma(3*(i-1)+1); sigma(3*(i-1)+

27、2)=Laplacian(i,j)*x(3*(j-1)+2)+sigma(3*(i-1)+2); sigma(3*(i-1)+3)=Laplacian(i,j)*x(3*(j-1)+3)+sigma(3*(i-1)+3); end f(3*(i-1)+1)=10*x(3*(i-1)+2)-10*x(3*(i-1)+1)-cc*sigma(3*(i-1)+1)-cc*d*(x(3*(i-1)+1)-balance(1); f(3*(i-1)+2)=28*x(3*(i-1)+1)-x(3*(i-1)+2)-x(3*(i-1)+1)*x(3*(i-1)+3)-cc*sigma(3*(i-1)+2)-

28、cc*d*(x(3*(i-1)+2)-balance(2); f(3*(i-1)+3)=x(3*(i-1)+1)*x(3*(i-1)+2)-(8/3)*x(3*(i-1)+3)-cc*sigma(3*(i-1)+3)-cc*d*(x(3*(i-1)+3)-balance(3); else%其余节点的动力学方程 for j=1:N sigma(3*(i-1)+1)=Laplacian(i,j)*x(3*(j-1)+1)+ sigma(3*(i-1)+1); sigma(3*(i-1)+2)=Laplacian(i,j)*x(3*(j-1)+2)+ sigma(3*(i-1)+2); sigma(3*(i-1)+3)=Laplacian(i,j)*x(3*(j-1)+3)+ sigma(3*(i-1)+3); end f(3*(i-1)+1)=10*x(3*(i-1)+2)-10*x(3*(i-1)+1)-cc*sigma(