%%% 将WS小世界网络和Tunable无标度网络合并形成路网和出行网自动生成机制；
N=81;m=1;IAM_iteration=10; %K为每个节点邻居节点个数的一半；
%% 给定方格形道路网络邻接矩阵，道路网络在上层网络变化过程中保持不变；
consequence_efficiency_sourcepair1=zeros(12,11); %存放源节点对效率的数组；
consequence_efficiency_sourcepair2=zeros(12,11);
consequence_congestion=zeros(12,11); %存放最终网络拥堵程度的数组；
act_err_Array=zeros(12,11); %存放流量重分配误差的数组；
consequence_iteration_time=zeros(12,11); %存放流量重分配总循环次数；
para_wgt_Array=[1,1,1,0.833,0.833,0.833,0.625,0.625,0.625,0.5,0.5,0.5]; %无标度网络PA参数待选数组；
%% 生成上层无标度网络前的准备程序；
wgt_Array=zeros(1,N);
for wsi=1:N
wgt_Array(wsi)=wsi;
end
%% 在保证节点初始权重不变的情况下演化生成不同度分布的上层无标度出行网络，并进行后续的一系列计算；
for para_wgt_i=1:12
travel_net=zeros(N,N);
vertex_wgt=zeros(1,N);
net_flag=zeros(N,N);
wgtsum=0;
for vwi=1:N
vertex_wgt(vwi)=1/(wgt_Array(vwi)^para_wgt_Array(para_wgt_i));
wgtsum=wgtsum+vertex_wgt(vwi);
end
connection_p=vertex_wgt/wgtsum;
accumulate_p(2:N+1)=cumsum(connection_p);
for iteration_time=1:N*m
random_data1=rand(1,1);
for rdj=1:N
if (random_data1>=accumulate_p(rdj) && random_data1<accumulate_p(rdj+1))
vertex_select1=rdj;
break
end
end
while length(find(net_flag(vertex_select1,:)==1))==N-1
random_data1=rand(1,1);
for rdj=1:N
if (random_data1>=accumulate_p(rdj) && random_data1<accumulate_p(rdj+1))
vertex_select1=rdj;
break
end
end
end
vertex_select2=vertex_select1;
while((vertex_select2==vertex_select1)||(net_flag(vertex_select1,vertex_select2)==1))
random_data2=rand(1,1);
for rdj=1:N
if (random_data2>=accumulate_p(rdj) && random_data2<accumulate_p(rdj+1))
vertex_select2=rdj;
break
end
end
end
travel_net(vertex_select1,vertex_select2)=1;travel_net(vertex_select2,vertex_select1)=1;
net_flag(vertex_select1,vertex_select2)=1;net_flag(vertex_select2,vertex_select1)=1;
end
%% 以下步骤用于计算无标度网络各节点的度；
degree_statistic_sf=zeros(1,N); %计算每个节点的度；
for ds_sf_i=1:N
for ds_sf_j=1:N
degree_statistic_sf(ds_sf_i)=degree_statistic_sf(ds_sf_i)+travel_net(ds_sf_j,ds_sf_i);
end
end
%% 定义出行网络中源节点并使用引力模型计算交通需求的分布情况；
travel_edge=length(find(travel_net~=0));
srk=1;
source_vertex_Array1=zeros(1,length(find(travel_net~=0))); %为存储交通源节点的向量s_point1和s_point2预留内存空间；
source_vertex_Array2=zeros(1,length(find(travel_net~=0)));
source_pair_volume=zeros(1,length(find(travel_net~=0))); %为存储交通源节点间交通量的向量q预留内存空间；
for srki=1:size(travel_net,1)
for srkj=1:size(travel_net,2)
if travel_net(srki,srkj)~=0
source_vertex_Array1(srk)=srki; %将交通发生源节点存入向量source_vertex_Array1；
source_vertex_Array2(srk)=srkj; %将交通吸引源节点存入向量source_vertex_Array2；
source_pair_volume(srk)=1000/travel_edge;
srk=srk+1;
end
end
end
%% 按边介数定义初始容限；
Traffic_impedance_betweeness=zeros(N,N); %不妨认为所有路段的阻抗值均相等且都设为单位1，定义路段初始阻抗；
Traffic_impedance_betweeness(road_net==1)=1; %道路网络矩阵中元素取1也即两节点之间相互连通存在路段，阻抗元素值即取1；
Traffic_impedance_betweeness(road_net==0)=inf; %若两节点间不存在路段则阻抗元素值取inf；
for Tib_i=1:N
Traffic_impedance_betweeness(Tib_i,Tib_i)=0; %网络中不允许有自连边，相应阻抗元素值为0；
end
road_betweeness=zeros(N,N); %边介数初始化
source_vertex_point=1; %source_vertex_point为存储源终节点对向量组指针，用于存取源终节点；
source_vertex2_sum=0;
while source_vertex_point<=length(source_vertex_Array1)
source_vertex1=source_vertex_Array1(source_vertex_point); %记录该指针所对应的源节点向量组中的节点标号；
source_vertex2_num=length(find(source_vertex_Array1==source_vertex1)); %记录该源节点所对应的全部终节点的数量；
source_vertex2_sum=source_vertex2_sum+source_vertex2_num;
origin_vertex=source_vertex1;
shortest_path=zeros(N,N);
for pti=1:N %按Floyd方法初始化存放路径节点的临时矩阵；
for ptj=1:N
if (pti~=ptj && Traffic_impedance_betweeness(pti,ptj)<inf)
shortest_path(pti,ptj)=pti;
else
shortest_path(pti,ptj)=-1;
end
end
end
Temp_ODshortest=Traffic_impedance_betweeness; %Temp_ODshortest矩阵存放节点间的最短出行时间；
Temp_m=1; %采用Floyd方法通过更新Temp_ODshortest矩阵求解节点间最短距离；
while Temp_m<=N
for Temp_mi=1:N
for Temp_mj=1:N
if(Temp_m==Temp_mi||Temp_m==Temp_mj)
continue;
elseif Temp_ODshortest(Temp_mi,Temp_mj)>Temp_ODshortest(Temp_mi,Temp_m)+Temp_ODshortest(Temp_m,Temp_mj)
Temp_ODshortest(Temp_mi,Temp_mj)=Temp_ODshortest(Temp_mi,Temp_m)+Temp_ODshortest(Temp_m,Temp_mj);
end
end
end
Temp_m=Temp_m+1;
end
Traffic_impedance_origin=zeros(N,N); %求源节点到其他节点的可达矩阵；
for obi=1:N
for obj=1:N
if Traffic_impedance_betweeness(obi,obj)~=0 && Traffic_impedance_betweeness(obi,obj)~=inf
Traffic_impedance_origin(obi,obj)=Traffic_impedance_betweeness(obi,obj)+Temp_ODshortest(origin_vertex,obi);
else
Traffic_impedance_origin(obi,obj)=Traffic_impedance_betweeness(obi,obj);
end
end
end
ODshortest_matrix=zeros(N,N); %求后继节点矩阵；
for omi=1:N
for omj=1:N
if omi~=omj && Traffic_impedance_origin(omi,omj)~=inf && Traffic_impedance_origin(omi,omj)==Temp_ODshortest(origin_vertex,omj)
ODshortest_matrix(omi,omj)=1;
end
end
end
%% 在取定某一源节点的情况下依次改变其对应终节点，计算各节点对间的最短路径及各边上的介数；
while source_vertex_point<=source_vertex2_sum
destination_vertex=source_vertex_Array2(source_vertex_point); %在与源节点相对应的终节点组中取出一个节点；
path_pyramid=zeros(1,1); %预定义存储最短路径节点的金字塔矩阵；
path_pyramid(1,1)=destination_vertex; %最短路径金字塔中的第一个节点为终节点；
omj=destination_vertex;
link=2;colume=1; %指定金字塔矩阵的行与列；
omi=1;
while(omi<=N && colume<=length(find(ODshortest_matrix(:,omj)==1))) %将终节点对应的紧前节点赋予金字塔矩阵的第二行；
if ODshortest_matrix(omi,omj)==1 %至此金字塔矩阵的第一、二行完成，作为以后最短路径计算的基础；
path_pyramid(link,colume)=omi;
colume=colume+1;
end
omi=omi+1;
end
upper=zeros(colume-1,1);
store_Array_pre=zeros(link,colu
