[2019-7-20]一种简单的级联失效模型

博主“时光荏苒心依旧”在一篇关于复杂网络MATLAB工具箱的博文中，曾介绍过一个生成随机无标度网络、随机攻击、有目的攻击的教程，复杂网络MATLAB工具箱，我在一个星期前就这个教程开始正式学习Matlab编程及数学建模，目的是为了研究复杂网络级联失效现象，网上很难找到研究级联失效的程序代码，所以得自力更生。

关于级联失效的具体定义和表现，由于时间有限，在此不做阐释，有兴趣的有需要的自然知道。

我花了一天时间百度各个Matlab函数的用法含义，虽然Matlab自带help函数，但是既看不太懂，即便翻译成中文理解起来也费劲，所以，我疯狂百度，编辑了一份属于自己的帮助文档“Matlab小贴士”，从刚开始连load、rand、round、plot这一类最简单的函数都看不懂，到逐渐能看懂简单的代码，上文的教程一中一段循环语句我就看了一整天，一句一句分析，像做英语精读似的，后面都用%标注了中文说明，刚开始我连%符号都不知道啥用法，后来知道“%8.6f\n”是啥意思。然后当了两天程序员，第一天在上文的教程一中直接修改，试图得到自己想要的功能，调试要人命，要人命啊，一直在报错和即将报错的道路上，郁闷；第二天自己尝试重新写个，考虑到能力问题，将难度标准（要实现的功能）一降再降，最终运行。最后花了一天，提高要求，重新又编写了一段程序，最终调试了几下，成功运行，基本能够实现随机攻击下，级联失效的最简单模拟。当然可能还有些小问题，可还不能让萌新小嘚瑟一下嘛。

ps:看到有好几个童鞋说，Clustering_Coefficient这个报错，我可能忘了在底下加上这个自定义函数的代码了，不好意思，我在文章最后面加上了，不知道行不行。。。。。。

%% 生成一个完全连通BA无标度网络
[A,NumberOfNodes] = BA_net();
%function [A,NumberOfNodes]=BA_net()
%input未增长前的网络节点个数m0:
%input每次引入的新节点时新生成的边数m：
%input增长后的网络规模N：
%input初始网络时m0个节点的连接情况：1表示都是孤立；2表示构成完全图；3表示随机连接一些边初始网络情况1，2或3：
%output该随机图的平均路径长度为：
%output该随机图的聚类系数为：
%output该随机图的平均度为：


[Custer,aver_Custer]=Clustering_Coefficient(A);%聚类系数;平均聚类系数
fprintf('Average Clustering Coefficient: %3.4f%%\n',aver_Custer);
[DeD,aver_DeD]=Degree_Distribution(A);%度;平均度
fprintf('Average Node Degree: %2.4f\n',aver_DeD);
[D,aver_D]=Aver_Path_Length(A);%路径长度;平均路径长度
fprintf('Average Path Length: %2.4f\n',aver_D);

负载重分配模型
% 负载重分配模型
%节点初始负载F
%节点容量C
%网络节点i的  初始负载F_i  与  节点本身的度k_i  相关：F_i = rho * k_i^tau
%网络中的节点 容量C_i  与初始负载F_i  成正比：C_i=(1+alpha)*F_i
%节点i失效后，失效节点的负载分配到完好节点j上，按一定比例
%负载更新
%剩余节点的更新后负载F_j;             if   F_j>C_j   节点j失效，否则end
%网络中所有节点的负载不超过其本身容量，连锁故障过程结束

%度量
%度量连锁故障结束后网络中失效节点个数CF_i
%CFN=sum(CF_i)/(N*(N-1))

rho=4;
tau=1.6;  %负载参数（0.1,2）
beta=1.6; %和度有关的负载分配比例参数
theta=1;   %和距离有关的负载分配比例参数
alpha=0.1;  %容量参数

F = rho * DeD.^tau;%初始负载，1*N

C = (1+alpha).*F;%容量，1*N,一行矩阵，默认为固定值



f = (D.^theta) .* ( DeD.^beta)./sum(sum((D.^theta) .*  (DeD.^beta)));%负载分配比例，方阵
% F_Temp = F_Temp + F_Remove * f;  %更新后的临时负荷矩阵，1*（N-1）
% Fail=F_Temp > C_First;   %判断剩余节点更新后的负荷是否大于节点容量，是=1，否=0，1*（N-1）,即接下来要删除的节点位置

%%随机删除某一节点，删除节点所在行和列，生成新的邻边矩阵A_Temp
A_Temp=A;
% aver_D_Temp=aver_D;
% aver_DeD_Temp=aver_DeD;

NodesToRemove=[];
% DeD_average=aver_DeD;
% D_average=aver_D;

Fail_All=[];    %失效节点集合，1*50矩阵.
NodesToRemove_Free=[];

step1=0;    %监测for循环次数.
step2=0;    %监测第一个if循环次数.
step3=0;    %监测第二个if循环次数.
step4=0;    %监测while中的if循环次数.

NodesToRemovePerStep =1;    %每步移除节点数.
RemainingNodes=1:NumberOfNodes; %剩余节点数，1行N列向量.
NumberOfNodes_Temp=NumberOfNodes;   %临时节点数，数值.
for i=1:50
        F_Temp=F;   %重置F_Temp;
        C_Temp=C;   %重置C_Temp；
        f_Temp=f;   %重置f_Temp；
        
        NodeIndecesToRemove=randperm(numel(RemainingNodes),NodesToRemovePerStep); %随机抽取移除节点，序号.
        NodesToRemove_Temp1 = RemainingNodes(:,NodeIndecesToRemove);  %移除节点的初始序号，1*1.
        NodesToRemove_Free=[NodesToRemove_Free,NodesToRemove_Temp1];    %随机抽取的50个节点的初始序号.
        NodesToRemove=[NodesToRemove,NodesToRemove_Temp1];   %所有移除节点的初始序号，矩阵.
        step1=step1+1;

        %移除节点后，更新网络
        RemainingNodes(:,NodeIndecesToRemove)=[];       %总的节点中的剩余节点，将每次随机抽样的节点剔除，避免重复抽取.
        
        RemainingNodes_Temp=1:NumberOfNodes;            %临时剩余节点重置为1*N.
        RemainingNodes_Temp(:,NodesToRemove_Temp1)=[];   %更新临时剩余节点，1*（N-1）.
        NodesToRemove_If=[];        %重置每一步if循环后的移除节点集合.
        NodesToRemove_If=[NodesToRemove_If,NodesToRemove_Temp1]; %更新，1*1.
        A_Temp=A;       %重置临时网络矩阵A.
        A_Temp=A_Temp(RemainingNodes_Temp,RemainingNodes_Temp);   %更新临时网络矩阵A，(N-1)*(N-1)；

        [DeD,aver_DeD_Temp,Fail_Temp,Fail_Num]=Degree_Distribution_Nofigure(A_Temp);    %更新节点度，平均度，孤立点的位置，孤立点点数量
        
        Fail_Sum=0;     %重置
        
        
        if ~isempty(Fail_Temp)   %判断是否存在孤立点，有的话，进行if
            step2=step2+1;
            %%算上孤立点，孤立点的效果等同于移除节点，两者一前一后连续发生
            
            Fail_Sum=Fail_Sum+Fail_Num; %总的孤立点数量
            NodesToRemove_Temp2=RemainingNodes_Temp(:,Fail_Temp);     %孤立点的初始序号，矩阵.
            NodesToRemove=[NodesToRemove,NodesToRemove_Temp2];
            NodesToRemove_If=[NodesToRemove_If,NodesToRemove_Temp2];   %1个移除节点+孤立点的初始序号，1*（Fail_Num+1）矩阵.
            RemainingNodes_Temp(:,Fail_Temp)=[]; %剩余的节点，初始序号，矩阵1*（N-1-Fail_Num）.
            
            
            %%直接算（移除节点+孤立点）带来的负载影响
            F_Remove=F_Temp(:,NodesToRemove_If); %（移除节点+孤立点）的负载，负载是一直变化的.
            %f = (D.^theta) .* ( DeD.^beta)./sum(sum((D.^theta) .*  (DeD.^beta)));
            %%负载分配比例，方阵，简单起见，认为f固定不变
            
            f_Temp=f(NodesToRemove_If,:);    %选择（移除节点+孤立点）所在行.
            f_Temp(:,NodesToRemove_If)=[];         %删除（移除节点+孤立点）的列.
            F_Temp(:,NodesToRemove_If)=[];   %临时负荷中，（移除节点+孤立点）所在列.
            F_Temp = F_Temp + F_Remove * f_Temp;  %更新剩余剩余负荷，即临时负荷 ，一行矩阵.
            
            %%判断临时负荷F_Temp 和 对应容量C_Temp 大小关系，if F_Temp>C_Temp，节点（级联）失效
            %%节点（级联）失效，效果等同于 孤立点，移除节点
            C_Temp = C(:,RemainingNodes_Temp);   %剩余节点对应的容量.
            
            NodesFailure=F_Temp>C_Temp;     %一行逻辑矩阵,找到F_Temp>C_Temp 具体位置.
            
        else
            %%移除节点带来的负载影响
            F_Remove=F_Temp(:,NodesToRemove_Temp1); %移除节点的负载，负载是一直变化的,1*1.
            
            f_Temp=f(NodesToRemove_Temp1,:);    %选择移除节点所在行,1*N.
            f_Temp(:,NodesToRemove_Temp1)=[];         %删除所有已经移除节点的列,1*(N-1).
            F_Temp(:,NodesToRemove_Temp1)=[];   %临时负荷中，删除移除节点所在列，1*(N-1).
            F_Temp = F_Temp + F_Remove * f_Temp;  %更新剩余剩余负荷，即临时负荷 ，一行矩阵.
            C_Temp = C(:,RemainingNodes_Temp);   %剩余节点对应的容量，1*(N-1).
            NodesFailure=F_Temp>C_Temp;     %一行逻辑矩阵,找到F_Temp>C_Temp 具体位置.
            
        end
        
        while sum(NodesFailure)>0   %级联失效，引发进一步的节点失效，即节点删除
            step3=step3+1;
            NodeIndecesToRemove=find(NodesFailure);  %找到逻辑矩阵中，逻辑值1所在位置，即为级联失效节点位置.
            %F_Remove=F_Temp(:,NodeIndecesToRemove);  %移除节点的负载，负载是一直变化的.
            Fail_Sum=Fail_Sum+numel(NodeIndecesToRemove);
            
            NodesToRemove_Temp = RemainingNodes_Temp(:,NodeIndecesToRemove);  %级联失效节点的初始序号，矩阵.
            NodesToRemove=[NodesToRemove,NodesToRemove_Temp];
            NodesToRemove_If=[NodesToRemove_If,NodesToRemove_Temp];   %所有(移除节点+孤立点+级联失效点）的初始序号，矩阵.
            
            %移除节点后，更新网络
            RemainingNodes_Temp(:,NodeIndecesToRemove)=[];   %剩余节点的初始序号，矩阵.
            RemainingNodes_Temp1=1:length(RemainingNodes_Temp);
            
            %C_Temp=C;   %重置C_Temp；
            %f_Temp=f;   %重置f_Temp；
            A_Temp=A;
            A_Temp=A_Temp(RemainingNodes_Temp,RemainingNodes_Temp);   %剩余节点的行和列；

            [DeD,aver_DeD_Temp,Fail_Temp,Fail_Num]=Degree_Distribution_Nofigure(A_Temp);    %更新节点度，平均度，孤立点的位置，孤立点点数量
         
            %%判断级联失效，是否会带来新的  孤立点， 以及  新一轮的  级联失效.
         if ~isempty(Fail_Temp)   %判断是否存在孤立点，有的话，进行if
            step4=step4+1;
            
            %%算上孤立点
            Fail_Sum=Fail_Sum+Fail_Num; %总的孤立点数量
            NodesToRemove_Temp2=RemainingNodes_Temp(:,Fail_Temp);     %孤立点的初始序号，矩阵.
            NodesToRemove_Temp3=RemainingNodes_Temp1(:,Fail_Temp);      %孤立点的前一次序号.
            NodesToRemove=[NodesToRemove,NodesToRemove_Temp2];
            NodesToRemove_If=[NodesToRemove_If,NodesToRemove_Temp2];   
            RemainingNodes_Temp(:,Fail_Temp)=[]; %剩余的节点，初始序号.
            
            
            %%直接算（移除节点+孤立点）带来的负载影响
            F_Remove=F_Temp(:,[NodeIndecesToRemove,NodesToRemove_Temp3]); %（级联失效节点+孤立点）的负载，负载是一直变化的.
            %f = (D.^theta) .* ( DeD.^beta)./sum(sum((D.^theta) .*  (DeD.^beta)));
            %%负载分配比例，方阵，简单起见，认为f固定不变
            
            f_Temp=f([NodesToRemove_Temp,NodesToRemove_Temp2],:);    %选择（级联失效节点+孤立点）所在行.
            f_Temp(:,NodesToRemove_If)=[];         %删除（移除节点+孤立点+级联失效）的列.
            F_Temp(:,[NodeIndecesToRemove,NodesToRemove_Temp3])=[];   %临时负荷中，删除（级联失效节点+孤立点）所在列.
            F_Temp = F_Temp + F_Remove * f_Temp;  %更新剩余剩余负荷，即临时负荷 ，一行矩阵.
            
            %%判断临时负荷F_Temp 和 对应容量C_Temp 大小关系，if F_Temp>C_Temp，节点（级联）失效
            %%节点（级联）失效，效果等同于 孤立点，移除节点
            C_Temp = C(:,RemainingNodes_Temp);   %剩余节点对应的容量.
            
            NodesFailure=F_Temp>C_Temp;     %一行逻辑矩阵,找到F_Temp>C_Temp 具体位置.

        else
            %%移除节点带来的负载影响
            F_Remove=F_Temp(:,NodeIndecesToRemove); %级联失效节点的负载，负载是一直变化的,1*1.
            
            f_Temp=f(NodesToRemove_Temp,:);    %选择级联失效节点所在行,1*N.
            f_Temp(:,NodesToRemove_If)=[];         %删除所有已经移除节点的列,1*(N-1).
            F_Temp(:,NodeIndecesToRemove)=[];   %临时负荷中，删除移除节点所在列，1*(N-1).
            F_Temp = F_Temp + F_Remove * f_Temp;  %更新剩余剩余负荷，即临时负荷 ，一行矩阵.
            C_Temp = C(:,RemainingNodes_Temp);   %剩余节点对应的容量，1*(N-1).
            NodesFailure=F_Temp>C_Temp;     %一行逻辑矩阵,找到F_Temp>C_Temp 具体位置.

        end
            
   
        end
        Fail_All=[Fail_All,Fail_Sum];       %更新失效节点集合,1*i.
        

end
CFN=sum(Fail_All)/(50*(50-1));
fprintf('平均失效规模: %8.6f%\n',CFN);

以上是一种简单的级联失效模型，我目前也才编出随机攻击策略的实现过程，至于蓄意攻击的还没写，想来改改就差不多了，以后再说。

以下是其中要用到的函数，计算节点度和平均度，计算节点距离，生成BA无标度网络，我都做了些许修改，总的功能没变，无非是原先有画图，我给改成输出值了，还有为了达到一定功能，我增加了些值，比如孤立点的位置。

函数：节点的度

function [DeD,aver_DeD,Fail,Fail_Num]=Degree_Distribution_Nofigure(A)
%% 求网络图中各节点的度及度的分布曲线
%% 求解算法：求解每个节点的度，再按发生频率即为概率，求P(k) 
%A————————网络图的邻接矩阵
%DeD————————网络图各节点的度分布
%aver_DeD———————网络图的平均度
%Fail-------------孤立点的位置
%Fail_Num---------孤立点的个数
N=size(A,2);
DeD=zeros(1,N);
Fail=[];

for i=1:N
   % DeD(i)=length(find((A(i,:)==1)));
   DeD(i)=sum(A(i,:));
end
aver_DeD=mean(DeD);
for i=1:N
    if DeD(i)==0         %判断是否为孤立点.
       Fail(end+1)=i; 
       
    else
    end
end
Fail_Num=length(Fail);
if sum(DeD)==0
    disp('该网络图全部由孤立点组成');
    return;
else 
end

函数：节点的距离以及平均路径长度

function [D,aver_D,disconnected]=Aver_Path_Length_Simple(A)
%% 求复杂网络中两节点的距离以及平均路径长度
%% 求解算法：首先利用Floyd算法求解出任意两节点的距离，再求距离的平均值得平均路径长度
%  A————————网络图的邻接矩阵
%  D————————返回值：网络图的距离矩阵
%  aver_D———————返回值：网络图的平均路径长度
%  disconnected--------返回值：1，不连通，存在孤立点；2，连通，不存在孤立点； 
 N=size(A,2);
 D=A;
 D(D==0)=inf;    %将邻接矩阵变为邻接距离矩阵，两点无边相连时赋值为inf，自身到自身的距离为0.
 for i=1:N           
     D(i,i)=0;   %对角线都为0.      
 end   
 for k=1:N            %Floyd算法求解任意两点的最短距离
     for i=1:N
         for j=1:N
             if D(i,j)>D(i,k)+D(k,j)
                D(i,j)=D(i,k)+D(k,j);
             end
         end
     end
 end
 aver_D=sum(sum(D))/(N*(N-1));  %平均路径长度
 if aver_D==inf
     disconnected=1;
 else
     disconnected=0;
 end

生成想要的BA无标度网络，略有修改

function [A,NumberOfNodes]=BA_net()
%%% 从已有的m0个节点的网络开始，采用增长机制与优先连接的机制生成BA无标度网络
%% A ——————返回生成网络的邻接矩阵
m0=input('未增长前的网络节点个数m0:  ');
m=input(' 每次引入的新节点时新生成的边数m： ');
N=input('增长后的网络规模N： ');
disp('初始网络时m0个节点的连接情况：1表示都是孤立；2表示构成完全图；3表示随机连接一些边');
pp=input('初始网络情况1，2或3： ');
if m>m0
    disp('输入参数m不合法');
    return;
end
x=100*rand(1,m0);
y=100*rand(1,m0);

switch  pp
    case 1
        A=zeros(m0);
    case 2
        A=ones(m0);
        for i=1:m0
            A(i,i)=0;
        end
    case 3
        for i=1:m0
            for j=i+1:m0
                p1=rand(1,1);
                if p1>0.5 
                    A(i,j)=1;A(j,i)=0;
                end
            end
        end
    otherwise
        disp('输入参数pp不合法');
        return;          
end 

for k=m0+1:N
    M=size(A,1);
    p=zeros(1,M);
    x0=100*rand(1,1);y0=100*rand(1,1);
    x(k)=x0;y(k)=y0;
    if isempty(find(A==1, 1))
        p(:)=1/M;
    else
         for i=1:M
             p(i)=length(find(A(i,:)==1))/length(find(A==1));
         end
    end
    pp=cumsum(p);          %求累计概率
    for i=1:m              %利用赌轮法从已有的节点中随机选择m个节点与新加入的节点相连
        random_data=rand(1,1);
        aa=find(pp>=random_data);jj=aa(1); % 节点jj即为用赌轮法选择的节点
        A(k,jj)=1;A(jj,k)=1;
    end
end
plot(x,y,'ro','MarkerEdgeColor','g','MarkerFaceColor','r','markersize',8);
hold on;

for i=1:N 
    for j=i+1:N
        if A(i,j)~=0
            plot([x(i),x(j)],[y(i),y(j)],'linewidth',1.2); 
            hold on;      %% 画出BA无标度网络图
        end
    end
end
axis equal;
hold off  

[C,aver_C]=Clustering_Coefficient(A);
[DeD,aver_DeD]=Degree_Distribution(A);
[D,aver_D]=Aver_Path_Length(A); 
NumberOfNodes=N;
 disp(['该随机图的平均路径长度为：',num2str(aver_D)]);  %%输出该网络的特征参数
 disp(['该随机图的聚类系数为：',num2str(aver_C)]);
 disp(['该随机图的平均度为：',num2str(aver_DeD)]);   

任重而道远，不过总体是快乐和充实的，因为既学到了新东西，又能取得可见到的成果。比较起来，感觉以前大学本科四年完全是在睡梦中打酱油——不值一提。一个星期学到的东西比以前大半年学到的还多。反正比实习打工要有趣得多。

以上就是过去一个星期的基本总结，喝水不忘挖井人，希望也能帮到别的手足无措的小伙伴，薪火相传，再接再厉。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
以下是我之前自学做的分析笔记，不算是代码吧，里面很多都没加%符号，是博主“时光荏苒心依旧”在一篇关于复杂网络MATLAB工具箱的博文中教程一的一小段，纯属晒晒笔记，嘚瑟一下。刚开始看的时候，惊为天书，后来越看越困惑，接着又开始嫌弃了，模型确实太简单了点，比我的简单模型还要简单，简单的平方。

%% remove random node untill the graph breaks apart:随机删除节点直到图分裂
TempGraph=Graph;   %图
NodesToRemovePerStep =1;   %每步移除的节点数
NumbersOfNodes = []; 			%节点数量
NumbersOfLinks = [];			%连边的数量
NumbersOfComponents = [];		%集群数量
LargestClusterSizes = [];			%最大聚类规模
SecondLargestClusterSizes = [];			%第二大聚类规模
 
RemainingNodes = 1:NumberOfNodes;    %剩余节点
 
while ~isempty(RemainingNodes)		%判断（）是否为空，空为0，否则为1
NodeIndecesToRemove = unique(round(rand(NodesToRemovePerStep,1)*(numel(RemainingNodes)-1))+1);
%% 随机移除节点
%% rand(NodesToRemovePerStep,1) 产生由在(0, 1)之间均匀分布的随机数组成的数组，每步移除的节点数
numel(RemainingNodes)  计数，剩余节点数
rand(NodesToRemovePerStep,1)*(numel(RemainingNodes)-1)
round(                     )  四舍五入
unique(round(                     )+1)		去掉矩阵中重复的元素
%%  重新计算移除节点后的图像及结构特点
NodesToRemove = RemainingNodes(NodeIndecesToRemove);
%%	移除的节点    = 	剩余节点数（NodeIndecesToRemove）

RemainingNodes = setdiff(RemainingNodes,NodesToRemove);
%%  剩余节点   =   剩余节点 中有而 移除的节点  中没有的值

TempGraph = mexGraphNodeRemove(TempGraph,NodesToRemove);
%%  图 	=	从原 图 中移除节点以及节点的边

NumbersOfNodes(end+1) = GraphCountNumberOfNodes(TempGraph);
%%  节点数量   (end+1)  =   计算	  图   中  节点数量

NumbersOfLinks(end+1) = GraphCountNumberOfLinks(TempGraph);
%%  连边的数量   增加一个长度，存放     =    计算   图   中  连边的数量

    if NumbersOfLinks(end)>0			%如果>0，否则=0   最后一个 连边的数量  元素  
        Components = mexGraphConnectedComponents(TempGraph);
		%% 集群   =    计算  图 中最大连通集群
        NumbersOfComponents(end+1) = numel(Components);
		%%  集群数量  增加一个长度，存放    =   统计 集群 数量
        ComponentsSizes = sort(cellfun('length',Components),'descend');  
		%%  集群大小      =	按照所有集群的长度    降序  排序
		cellfun('length',Components)  计算所有集群的长度，
		
		
        if ~isempty(ComponentsSizes)		%  	如果=1，否则=0 	判断（集群大小）是否为空，空为0，否则为1
            LargestClusterSizes(end+1) = ComponentsSizes(1);		%  最大聚类大小   增加一个长度，存放 =最大的集群大小
        else
            LargestClusterSizes(end+1) = 0;		%  最大聚类大小   增加一个长度，存放=0
        end
        if numel(ComponentsSizes)>1		%  如果>1，否则 
            SecondLargestClusterSizes(end+1) = ComponentsSizes(2);	%  第二大聚类规模    = 最大的集群大小
        else
            SecondLargestClusterSizes(end+1) = 0;
        end
    else
        NumbersOfComponents(end+1) = 0;
        LargestClusterSizes(end+1) = 0;
        SecondLargestClusterSizes(end+1) = 0;
    end
end
h4 = figure;
plot(NumbersOfComponents,'r');
hold on;				%画多重线
h5 = figure;
plot(NumbersOfNodes,'r');
hold on;
h6 = figure;
plot(NumbersOfLinks,'r');
hold on;
h7 = figure;
plot(SecondLargestClusterSizes,'r');
hold on;
h8=figure;
plot(LargestClusterSizes,'r');
hold on;
 
%% remove most connected nodes (with heighest outgoing degree):删除出度最大的节点
TempGraph=Graph;											%图
NodesToRemovePerStep =1;									%每步移除的节点数
NumbersOfNodes = [];										%节点数量
NumbersOfLinks = [];										%连边的数量
NumbersOfComponents = [];								%集群数量
LargestClusterSizes = [];								%最大聚类规模
SecondLargestClusterSizes = [];							%第二大聚类规模
 
RemainingNodes = 1:NumberOfNodes;						%剩余节点数
 
while ~isempty(TempGraph.Data)							%判断（图中数据）是否为空，空为0，否则为1
    Degrees = GraphCountNodesDegree(TempGraph);		%  节点度
    [OutDegrees, SortOrder]=sort( Degrees(:,3),'descend');		% [出度，位置序号]=按出度降序排列
NodesToRemove = Degrees(SortOrder(1:min([numel(SortOrder) NodesToRemovePerStep])));
%%  根据节点度大小，移除出度最大的节点
%%  移除的节点   =  根据 节点度（出度降序排列的位置序号（））
min([numel(SortOrder) NodesToRemovePerStep])   两者中最小的
numel(SortOrder)  计数  位置序号  有度的节点的个数
NodesToRemovePerStep   每步移除的节点数
%%  重新计算移除节点后的图像及结构特点
    TempGraph = mexGraphNodeRemove(TempGraph,NodesToRemove);				%  移除的节点
    NumbersOfNodes(end+1) = GraphCountNumberOfNodes(TempGraph);			%  节点数量
    NumbersOfLinks(end+1) = GraphCountNumberOfLinks(TempGraph);			%  连边的数量
    if NumbersOfLinks(end)>0
        Components = mexGraphConnectedComponents(TempGraph);
        NumbersOfComponents(end+1) = numel(Components);
        ComponentsSizes = sort(cellfun('length',Components),'descend');
        if ~isempty(ComponentsSizes)
            LargestClusterSizes(end+1) = ComponentsSizes(1);
        else
            LargestClusterSizes(end+1) = 0;
        end
        if numel(ComponentsSizes)>1
            SecondLargestClusterSizes(end+1) = ComponentsSizes(2);
        else
            SecondLargestClusterSizes(end+1) = 0;
        end
    else
        NumbersOfComponents(end+1) = 0;
        LargestClusterSizes(end+1) = 0;
        SecondLargestClusterSizes(end+1) = 0;
    end
end

最后真诚感谢博主博主时光荏苒心依旧，半师之谊，上海的奶茶给您满上。

Ps**:Clustering_Coefficient**

help Clustering_Coefficient
求网络图中各节点的聚类系数及整个网络的聚类系数
% 求解算法：求解每个节点的聚类系数，找某节点的所有邻居，这些邻居节点构成一个子图
% 从A中抽出该子图的邻接矩阵，计算子图的边数，再根据聚类系数的定义，即可算出该节点的聚类系数
A————————网络图的邻接矩阵
C————————网络图各节点的聚类系数
aver———————整个网络图的聚类系数

function  [C,aver_C]=Clustering_Coefficient(A)
%% 求网络图中各节点的聚类系数及整个网络的聚类系数
%% 求解算法：求解每个节点的聚类系数，找某节点的所有邻居，这些邻居节点构成一个子图
%% 从A中抽出该子图的邻接矩阵，计算子图的边数，再根据聚类系数的定义，即可算出该节点的聚类系数
%A————————网络图的邻接矩阵
%C————————网络图各节点的聚类系数
%aver———————整个网络图的聚类系数
N=size(A,2);
C=zeros(1,N);
for i=1:N
    aa=find(A(i,:)==1);  %寻找子图的邻居节点
    if isempty(aa)
        disp(['节点',int2str(i),'为孤立节点，其聚类系数赋值为0']);
        C(i)=0;
    else        
        m=length(aa); 
        if m==1
            disp(['节点',int2str(i),'只有一个邻居节点，其聚类系数赋值为0']);
            C(i)=0;
        else
           B=A(aa,aa);          % 抽取子图的邻接矩阵
           C(i)=length(find(B==1))/(m*(m-1));
        end
    end
end
aver_C=mean(C);

    
    
    

最后，附上后来建的一个小的级联失效扣扣交流群：#530#762#903#，不怎么活跃，人也不多，上传了一点点文件在群里，希望能起到一点帮助。