@[TOC]（复杂⽹络级联失效的代码）
复杂⽹络级联失效的代码
function [CFN,NodesFailure,A_temp]=cascading_failure(A)
%%A为复杂⽹络的邻接矩阵
NumberOfNodes=100;
%%%%%%删除部分⽹络节点，进⾏⽹络级联失效%%%%%%% [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;%初始负载，1N
C = (1+alpha).F;%容量，1N,⼀⾏矩阵，默认为固定值
beta)./sum(sum((D.
f = (D.^theta) . ( DeD.theta) .* (DeD.^beta)));%负载分配⽐例，⽅阵,DeD是节点度分布
% 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 ~impty(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 具体位置.
el
%%移除节点带来的负载影响
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 ~impty(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 具体位置.
el
%%移除节点带来的负载影响
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);
