confilter116.m

本程序时使用小世界模型对多个智能体寻优结果的仿真

function [e2]=confilter116(duration, dt)
% Kalman filter simulation for consensus
% INPUTS
% duration = length of simulation (seconds)
% dt = step size (seconds)
N=1000;%% 总人数 agent number
flag_WS_average=0; %平均speed
measnoise=0.1; % measurement nosie Q[k].
communoise=0.01; %communication noise v_ij[k], omiga_ij=E{v_ij*(v_ij)T}=0.1.
%load v_ij_data v_ij;%communication noiseload v_ij_initial v_ij;%communication noise
%v_ij=zeros(1000);
v_ij=rand(100)*0.05;
Q_ij=0.1;
omiga=8.33*10e-5;
xhat=zeros(N,duration/dt+1);
%xhat(:,1)=randn(N,1);
load confilter_initial initial_data;
xhat(:,1)=initial_data;
xhat(20,1)=100;
P=zeros(N,duration/dt+1);
%P(:,1)=rand(N,1)*0.5;
load initialp P_initial;
P(:,1)=P_initial;
P(1,1)=0.001;
% network topology g_ii=1,g_ij=1 if information from j to i.
%N=1001;%% 总人数
%	WS-model generaton algorithm 
%
%1）从规则图开始：考虑一个含有N个点的最近邻耦合网络，每个节点都与它左右相邻的2k个节点相连。
%2）随机化：以概率p随机地重新连接网络中的每个边。p=0对应于原来的规则网络，p=1则对应于完全随机的网络，通过调节p的值就可以控制从完全规则网络到完全随机网络的过渡。
%Node----	the matrix of the network
%k----	the degree of the node is 2*k	
%N----	the number of the nodes in the network
%	生成网络的对角线为点的度的负值。

%function WS_small_world_Model()
%clear
%N=500;
k=4;
%p=0.9; %重联概率
loop=10;% 循环次数
times=1 % 次数初始化
flag=1;% 特征值对应不同重联概率标志。
%flag_WS=zeros(1,loop);
while(times<=loop)  %取平均值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	construct a nearest neighbored network model	%
Node=zeros(N);
for j=1:N
    for K=1:k  
        if j+K>N
            Node(j,j+K-N)=1;Node(j+K-N,j)=1;
        else            
            Node(j,j+K)=1;Node(j+K,j)=1;
        end
    end
end
for i=1:N
    Node(i,i)=-2*k;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nearest neighbored %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   change the edge
Nodes=Node;
for p=0:0.1:1   %   概率p在0-1之间改变
    Node=Nodes;
    for K=1:k
        for i=1:N         
            q=rand(1);
            w=1;             
            if q<p           %               
                while w<2
                    r=round((N+1)*rand(1));                  
                    if r>0&r<(N+1)                        
                        if Node(i,r)==0     %   avoid multi-edge
                            Node(i,r)=1;Node(r,i)=1;Node(r,r)=Node(r,r)-1;
                            KK=i+K;
                            if KK>N
                                KK=KK-N; 
                            end
                            Node(i,KK)=0;                                                     
                            Node(KK,i)=0;
                            Node(KK,KK)=Node(KK,KK)+1;
                            w=2;                       
                        end
                    end
                end     % end for 'while'               
            end  % end for 'if'
        end  % end for 'for'                   
    end

%---------------------------produce WS
G_ij=zeros(N,N);
%L_ij=Node;
G_ij=Node;
ones_v=ones(1,1000);
for aa=1:N
    G_ij(aa,aa)=1;
end 
%---------------------------------------------
xhat(:,1)=initial_data;
xhat(20,1)=100;
P=zeros(N,duration/dt+1);
P(:,1)=P_initial;
P(1,1)=0.001;
%-----------------------------------------------
for k = 1 : dt: duration/dt
    P(:,k+1)=((P(:,k)+Q_ij).^(-1)+G_ij*((P(:,k)+omiga).^(-1))).^-1;
    xhat(:,k+1)=xhat(:,k)-P(:,k+1).*((G_ij*((P(:,k)+omiga).^(-1))*ones_v.*eye(1000)-G_ij.*(((P(:,k)+omiga).^(-1))*ones_v)')*xhat(:,k))+P(:,k+1).*(G_ij.*(v_ij)*((P(:,k)+0.1).^(-1)));
   %----------------------------------根据偏离平均值偏差判断收敛时间。
    %if  mad(xhat(:,k))<=10e-3   % 若偏离平均值10e-3 
        %flag_WS(1,times)=k;
        %break;
        %end
    %----------------------------------
end
%---------------------------------------------求除1之外最大特征值
L=eye(1000)-P(:,duration+1)*ones_v.*(G_ij*((P(:,duration)+omiga).^(-1))*ones_v.*eye(1000)-G_ij.*(((P(:,duration)+omiga).^(-1))*ones_v)');
e=eig(L);
[s,v]=max(e);
e(v)=-100;                
e2(flag,times)=max(e);
flag=flag+1;
    if flag>11
        flag=1;
    end
end % end for the last 'for'
times=times+1;
%---------------------------------------------
end %--------for 'while'
%flag_WS_average=mean(flag_WS);
y1=P(1,:);
y2=P(90,:);
y3=P(450,:);
y4=P(20,:);
y5=P(500,:);
y11=xhat(1,:);
y22=xhat(90,:);
y33=xhat(450,:);
y44=xhat(20,:);
y55=xhat(500,:);
t=1:dt:duration+1;
%subplot(1,2,1),
plot(t,y1,t,y2,t,y3,t,y4,t,y5);
%subplot(1,2,2),
plot(t,y11,t,y22,t,y33,t,y44,t,y55);
grid