复杂网络：随机、蓄意打击和级联失效模型

          因为方向定的是复杂网络这块，写小论文的话最初说是做抗打击实验，就是基于度中心性和介数中心性两个指标的排序对节点进行随机打击、蓄意打击，之后又紧接着开始做级联失效这方面。刚开始也不懂代码，就在网上各种找，发现这方面的资料并不是很多，最后自己就利用暑假时间开始学习，现在把这些代码拿出来跟大家交流分享下，有问题的可以私信啊，这东西都是自己写的，肯定会存在不足啊，见谅

主要就是两部分：随机、蓄意打击和级联失效模型，用的是matlab。

1.随机、蓄意打击

% 初始网络性能
%生成随机数，以此进行随机攻击（注释掉即为蓄意攻击），随机数值改为你自己网络的节点数
%   Name_Struct.Node_Key_Degree = randperm(238);

这个是打击的时候切换随机和蓄意用的

%% 
%  按照 Degree 算法排序，删除节点
A = A_Init;          %% 网络邻接矩阵 A
B=[];                %%定义空数组，接收被删除的节点序号
for i = 1:NumDelete
    A( Name_Struct.Node_Key_Degree(i),: ) = 0;     %% 删除节点 Node_Key_Degree(i),用 0 占位，不能置空
    A( :,Name_Struct.Node_Key_Degree(i) ) = 0;
    AA = A;
    AA( sum(A)==0,: ) = [];
    AA( :,sum(A)==0 ) = [];
    B(i) = Name_Struct.Node_Key_Degree(i);
     Con_Index_NetEff = testEglob( AA );
     Eglob(i) = Con_Index_NetEff.Net_Eff_Mymod;
     Con_Index_NetEff = largestcomponent1 ( AA );
     largest(i) = Con_Index_NetEff.Net_Eff_Mymod;
     Con_Index_NetEff = largestcomponent ( AA );
   components(i) = Con_Index_NetEff.Net_Eff_Mymod;
    save  
end

%%
%接下来就是生成网络连通效率图
%Eglob存储了相应的网络效率的数值

这部分就是在网络中（一般用的是邻接矩阵，如果是对边进行打击的话是邻接表），把需要打击的节点从网络中给删除掉，之后计算指标，像网络效率、最大联通子图相对大小、联通片数量等。

function Eglob = ATestAver(str,numDelete)
%输入 ： 
%        str：意为数据文件的路径
%        numDelete:删除节点的个数

%返回值：Eglob，即网络效率值的数组

刚才上面那个是进行蓄意打击时用到的，随机打击的话就需要多次做实验取其平均值了，所以把代码改成上面的函数形式。

%定义随机攻击节点的个数，具体数值根据网络规模拟定
numDelete = 50;

%定义网络效率初始矩阵
netSum = zeros(1,numDelete);

%定义随机攻击的次数，也就是函数循环的次数
numRandom = 5;
for i=1:numRandom
    
    %把得到的网络效率数组赋给netI
    net1(i) = ATestAver('Eglob.mat',numDelete);
%     net2 = ATestAver('largest(i).mat',numDelete);
%      net3 = ATestAver('components(i).mat',numDelete);
    %累加
    netSum1(i)= netSum + net1(i);
% %     netSum2 = netSum + net2;
%     netSum3 = netSum + net3;
end

%求出平均值
netAver1 = netSum1(i)/numRandom;
% netAver2 = netSum2/numRandom;
% netAver3 = netSum3/numRandom;
save netAver

然后接着就是重复实验。

function Con_Index_NetEff = largestcomponent1 ( AA )

%A = triu(A);
n=length(AA);

for i=1:n
    mz{i}=find(AA(i,:));
end

x(1:n)=0;
z=0;
k=0;
for i=1:n
    if x(i)==0;
        z=z+1;
        clear v
        v(1)=i;
        while nnz(v)>0
                    x(v(1))=z;
                    k=union(k,v(1));
                    b=setdiff(mz{v(1)},k);
                    v=setdiff(v,v(1));
                    v=union(v,b);         
        end
    end
end

c(1:max(x))=0;
for i=1:max(x)
    c(i)=length(find(x==i)); 
end

largest = max(c)/238;
% components = size(c,2);
% zuidazitu = largest;
% result = struct('zuidazitu',zuidazitu );
%cm=find(c==max(c));
Net_Eff_Mymod = largest 
Con_Index_NetEff = struct('Net_Eff_Mymod',Net_Eff_Mymod);
%B=find(x==cm(1));
save

end

这个代码就是求最大联通子图相对大小的，可以改成求联通片数量，上面有的。

2.级联失效

这个模型有点复杂，当时想着先把最基本的ML模型给搞明白，然后再这个基础上取增加一些其他东西，像非线性负载-容量、节点过载、正常、失效和拥堵的节点状态。

%% 原始数据，这里的负载和容量直接给定，有很多负载容量模型可以嫁接修改
Node_load = [2 4 4 3 1 3 1 2 6 2];
Node_Capacity = 1.5*Node_load;

这是最开始的线性负载-容量模型，同样还有非线性模型。

 if i==0
         [new_isolated, new_isolated_node_num ] = find_isolated(A_temp);
        for m=1:new_isolated_node_num
            A_temp(m,:)=0;
            A_temp(:,m)=0;              % 将A_temp矩阵中孤立节点的连边关系移除
            
            index=find(node_index==new_isolated(m));
            node_index(index)=[];       % 将节点索引里面的孤立节点移除
            A_change(index,:)=[];
            A_change(:,index)=[];       % 将A_change矩阵中孤立节点移除
        end
        isolated_node = [isolated_node, new_isolated];
        total_failure = [total_failure, isolated_node];
    end

在中间需要求网络中的孤立点，这个我最后时算在失效集合里面的。

for m=1:length(new_failure)
        % 先找失效节点的邻居节点的编号
        neiber=find(A_temp(new_failure(m),:)==1);

找失效节点的邻居节点。

 if length(neiber)>1 &  size(A_temp,1)>0
                temp_neiber=unique(neiber);  %  去重
                
                alltype=combntns(temp_neiber,2);  %  所有的排列组合

这里面会出现一个问题，如果失效节点过多，且他们之间还有线连接的话，在负载分配时会影响数值结果，所以需要断开。

 % 计算额外增加的负载，这里可以有很多的公式类型
        delta(neiber) = deal(Node_load_temp(new_failure(m)) .* Node_load(neiber)/sum(Node_load(neiber)));
        Node_load_temp(neiber) = deal(Node_load_temp(neiber) + delta(Neiber)); % 负载重分配

然后进行负载分配，这个就是最初基于节点负载比例分配的。之后都是需要改成基于剩余容量进行分配。

  if Node_load_temp(kk)>Node_Capacity(kk)
            fail_node_temp = [fail_node_temp, kk]; % 得到本轮引起的失效节点集合

然后就是统计失效节点数。

  
            if  [Node_load_temp(kk)-Node_Capacity(kk)]<=0.15*Node_Capacity(kk) && Node_load_temp(kk)>=Node_Capacity(kk);
                pause_node_temp = [pause_node_temp,kk];                           % 得到本轮引起的暂停节点集合
            end

也可以增加拥堵节点，之后就从网络中把失效和拥堵删除掉，继续进行循环，负载分配。

 Standby_node1 = unique(fail_node_temp);

把失效节点传给下一个失效循环。

  fail_nodes = unique([fail_nodes,Standby_node1]);                % 累计失效数，（第一次失效+本轮失效）；
        
        % 累计总的失效数（第一次失效数，本轮新增孤立节点数，本轮新失效数）
        total_failure = unique([total_failure, new_isolated_node, Standby_node1]);
        
        i=i+1;
        figure
        
        
        Z(i) =  length(find(total_failure));                     % Z（i）为每一次级联后，累计总的失效数；
        Z_e(i) =  length(find(total_failure))/238;               % Z_e（i）为每一次级联后，累计总的失效数占比；
        Z_e_e =  length(find(total_failure))/238;                % 最后一次，累计总的失效数占比；
        
        S(i) = length( fail_nodes);                              % 累计部分失效节数；
        
        E(i) = length(unique(pause_node_temp));                  % 暂停节点个数（累计每一轮）
        F(i) = length(unique(fail_node_temp));                   % 失效节点个数（累计每一轮）(没去掉重复的)
        
        G(i) = Z(i)-S(i);                                        % G(i)为累计孤立点数
        
        
        
        Con_Index_NetEff = testEglob(A_temp);               % 网络效率，导出的数据最后一组为重复的，因为要符合条件才能停止
        Eglob(i) = Con_Index_NetEff.Net_Eff_Mymod;
        Con_Index_NetEff =largestcomponent1(A_temp);        % 最大连通子图尺寸（记得修改子函数中的分母）
        largest(i) = Con_Index_NetEff.Net_Eff_Mymod;
        Con_Index_NetEff =largestcomponent(A_temp);         % 连通片数量
        components(i)= Con_Index_NetEff.Net_Eff_Mymod;

这里是统计新的指标值。

             之后就是对模型的改进，可以根据网络结果的不同，修改负载分配策略，其中还有节点过载情况，同时也可以加上过载系数和概率。我是对网络负载分配策略进行改进，考虑了两个节点间的实际距离，然后根据距离的远近，筛选出合适的节点进行负载分配。

欢迎交流分享。。。