【复杂网络】—— 论文阅读及复现《Catastrophic cascade of failures in interdependent networks》

1 相互依赖网络的级联失效

这篇论文研究了相互依赖网络的级联失效,主要有以下几个要点:

1.1研究背景和目标

• 大多数以前的网络研究集中在单一网络,而现实中的许多网络都是相互依赖的。
• 相互依赖网络中的节点故障可能会引发级联故障,最终导致整个网络的碎裂。
• 文章提出一个相互依赖网络的模型,目标是分析网络碎裂的临界点。

1.2 研究方法

• 提出了一个包含两个网络A和B的相互依赖网络模型。
• 使用生成函数方法求解模型,得到相互巨分量的大小。
• 推导出求解相互依赖ER网络和无尺度网络的碎裂临界点的公式。

1.3主要结果

• 对于相互依赖的ER网络,求解得到碎裂的临界平均度为2.455,大于单一ER网络的1。
• 无尺度网络的碎裂临界点确实大于0,不再像单网络那样趋向于0。
• 度分布更宽泛的网络更脆弱,与单网络结果相反。

1.4讨论和启发

• 相互依赖网络中,高度节点变得脆弱,拓扑结构对稳定性的影响与单网络相反。
• 需要重新考虑设计稳健网络的方式,将相互依赖性纳入考量。
• 模型可推广到多个相互依赖网络。

1.5总结

• 文章开创性地研究了相互依赖网络的碎裂行为,得到许多与单网络不同的重要结论,对复杂网络研究具有重要启发。

2 论文复现部分

笔者复现部分主要用MATLAB(理论)以及python（模拟）。

2.1 理论图FIG1

function sum= GA0(x)
syms k;
s=((2/k).^2-(2/(k+1)).^2)*x.^k;
sum=double(symsum(s,k,2,inf));
end

function output = GA1(x)
syms k;
s=(1/k-1/(k+1)+(1/(k+1)).^2)*x.^(k-1);
sum=symsum(s,k,2,inf);
output=double(4*sum/(2*pi*pi/3-3));
end

function output = g(x)
opt = optimset('Display','off');
z=fsolve(@(z)GA1(z)-1+1/x-z/x,0.5,opt);
g=1-GA0(z);
output=g;
end

A=0:0.05:1;%横坐标
% p的值
p = [0.7,0.752,0.80];
 re = [zeros(1,21);zeros(1,21);zeros(1,21)];
 re(1,1) = 0;
 re(2,1) = 0;
 re(3,1) = 0;
 for i = 1:3
    p1 = p(i);
    for j  = 1:20
        x=j*0.05;
        re(i,j+1) = p1.*g(p1.*g(x));
        fprintf('p=(%.2f),pga(pgb(%.2f))= %12.5f\n',p1,x,re(i,j+1));
    end
 end
plot(A,re(1,:),'r',A,A,'-k')
hold on;
plot(A,re(2,:),'b')
plot(A,re(3,:),'g')
xlabel('x');
ylabel('pga(pgb(x))');
legend({'p=0.7','y=x','p=0.752','p=0.80'},'Location','best')

结果如下：

2.2 理论图FIG2

function re = rc(a , b)
opt = optimset('Display','off');
re = fsolve(@(r)exp((r.^a-1)*(1-r.^b)/(a*r.^a+b*r.^b-(a+b)*r.^(a+b)))-r,0.5,opt);
end

function re = pc(a,b,r)
opt = optimset('Display','off');
re = fsolve(@(p)(a*r.^a+b*r.^b-(a+b)*r.^(a+b)).^(-1)-p,1,opt);
re =re*a;
end

k = 0:0.05:1;
b = 3;
%pc
re = zeros(1,21);
re(1)=1;
%P∞
re2 = zeros(1,21);
re2(1)=0.01;
for i = 2:21
    a = b*k(i);
    r = rc(a,b);
    re(i) = pc(a,b,r); 
    %P∞
    re2(i)= re(i)*(1-r.^a)*(1-r.^b);
end
plot(k,re,'--k',k,re2,'-k','LineWidth',2)
ylim([0 2.5])
legend({'p_c','p_{inf}'},'Location','best','FontSize',18)
xlabel('a/b','FontSize',18)
ylabel('ap_{inf},ap_c','FontSize',18)

结果如下：

2.3 理论图FIG3

function output= ERGA0(x)

output = exp(8*(x-1));
end

function output= ERGA1(x)

output = dERGA0(x)/dERGA0(1);
end

function output= dERGA0(x)
output = 8*exp(8*(x-1));
end

function output = ERPA(x)
syms z;
opt = optimset('Display','off');
a=fsolve(@(z)ERGA1(z)-1+1/x-z/x,0.5,opt);
output = 1-ERGA0(a);
end

p =2.45/8;
a = 8;
b = 8;
re = zeros(1,80);
%中间变量
p1 = zeros(1,80);

re(1) = p;
p1(1) = p;

for i = 2:80
    re(i) = p1(i-1)*ERPA(p1(i-1));
    p1(i) = p*ERPA(p1(i-1));
    if re(i)<=0
        break
    end
end
data1 = importdata("Sheet1.csv");
data2 = importdata("Sheet2.csv");
data3 = importdata("Sheet3.csv");
data4 = importdata("Sheet4.csv");
data5 = importdata("Sheet.csv");
plot(1:80,re/p,'o-r',1:80,data1,'-k')
hold on
plot(1:80,data2,'-k')
plot(1:80,data3,'-k')
plot(1:80,data4,'-k')
plot(1:80,data5,'-k')
hold off
xlabel('n',FontSize=18);
ylabel('p_{n}/p',FontSize=18);
legend({'theory'},'Location','best')

结果如下：

2.4 模拟部分

Github链接：https://github.com/lingjun127/pythonProject2