@[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);
