1.算法仿真效果
matlab2022a仿真结果如下：
2.算法涉及理论知识概要
负环的定义：负环是指权值和为负数的环。负环会使图的最短路径计算陷入死循环，因此，存在负环的图不存在最短路。
负环的计算方法：
负环有两种计算方法，都是基于Bellman-Ford算法或者SPFA算法。
第一种算法是：统计每个点的入队次数，如果某个点入队大于等于n次，则说明有负环
第二种算法是：统计到某个点的最短路所经过点的个数，如果经过n个点，则说明存在负环。
(一般情况下，我们使用第二种算法)
由于当负环存在时，SPFA会陷入死循环，且n是非死循环的最坏情况。所以以上两种算法是正确的。
求负环算法的编程实现
首先将所有点的距离都赋值为0
然后将所有的点入队。
在100m*100m随机构建12个节点，节点距离在30m内通过有向线相连，仿真独立进行1000循环。其中N_j可设置为500~1000随机数。
通过使用N算法进行求解，具体N算法流程已在文中介绍，主要针对111.pdf中2，3部分进行仿真，即单用户，具有优先级的多用户，没有区别的多用户，其中针对多用户，源节点目的节点可按：s_1=1,d_1=12;s_2=2,d_2=10;s_3=3,d_3=9;s_4=5,d_5=11进行设置。
3.MATLAB核心程序
```USER = [2,3,5];
for jks = 1:length(USER);
Traffic_VolumeS = [20];
SR = [10:10:60];
Es3 = [];
for jjj = 1:length(SR)
Traffic_Volume = Traffic_VolumeS;
Loop = 9;
Es = zeros(MTKL,Loop);
rng(5*jjj);
for i = 1:MTKL
i
%随机12个节点
X = SCALE*rand(1,Note);
Y = SCALE*rand(1,Note);
%30以内有相线相连接
xk_ij = zeros(Note,Note);
w_ij = zeros(Note,Note);
b_ij = zeros(Note,Note);
for j1 = 1:Note
for j2 = 1:Note
dist = sqrt((X(j1)-X(j2))^2 + (Y(j1)-Y(j2))^2);
if dist <= Dis_R & j1~=j2
%Select a feasible route with f = fij,这里构造一个连接矩阵，将相互连接的用1表示
xk_ij(j1,j2) = 1;
%随机分配cost
b_ij(j1,j2) = 50*randn-30;
end
end
end
EZ = 0;
for js2 = 1:USER(jks)
%初始值
fij = (1+rand/5)*Traffic_Volume + 2*randn(Note,Note);
fijmax0 = max(max(fij));
%选择一个路径
[f,flag] = func_sel_route(xk_ij,X,Y,fij,b_ij,Dis_R);
%开始迭代
E_1 = zeros(1,Loop);
for j = 1:Loop
if flag == 1
%补图
Nj = 500 + 500*rand;
Ys = func_complementary_graph(xk_ij,f);
[R,C] = size(Ys);
a = randperm(10);
SR_(jjj) = SR(jjj);%模拟in的Rate
Ri_in = 1e6*SR_(jjj);%转换为M，
Ri_out = 1e6*SR(jjj);%转换为M，
if j == 1
E_1(1,j) = Traffic_Volume;%初始值
else
Ej_proc = Nj/2*ECMP_max + ECMP_min;
Ei_out = frp(Ri_out) * (MOEC_min+(MOEC_max-MOEC_min)*rand);
Ej_in = frp(Ri_in) * (MIEC_min+(MIEC_max-MIEC_min)*rand);
E_1(1,j) = Ei_out + Ej_in + Ej_proc;
end
for k1 = 1:R
for k2 = 1:C
E(k1,k2) =E_1(1,j)*(1+randn/5);
end
end
%再将M除掉
E=E/1e6;
fijmax_ = zeros(R,C);
for k1 = 1:R
for k2 = 1:C
if Ys(k1,k2) == 1
fijmax_(k1,k2) = fijmax0 - fij(k1,k2);
fij(k2,k1) = fijmax_(k1,k2);
E_(k1,k2) = E(k1,k2);
E(k2,k1) = -E_(k1,k2);
end
end
end
E_=-E;
E = 7*E/sqrt(SR(jjj));
%negative cost cycle
[flags2] = func_negative_cost_cycle(Ys,xk_ij,X,Y,fij,Dis_R,E,E_,f);
if flags2 == 1%算法结束
break;
end
f_=f.*E;
for k1 = 1:R
for k2 = 1:C
tmps = fijmax0 - fij;
[xc,yc] = find(tmps==0);
tmps(xc,yc)=fijmax0;
delta = min([min(min(tmps)),min(min(fij))]);
if E_(k1,k2) > 0
fij(k2,k1) = fij(k2,k1)+delta;
else
fij(k2,k1) = fij(k2,k1)-delta;
end
end
end
E_1(1,j) = sum(sum(f_.*fij));
else
E_1(1,j) = 0;
end
end
EZ = EZ + E_1(end);
end
Es(i,:) = EZ;
end
%对1000迭代进行平均
Es2 = [];
for i = 1:Loop
tmps = Es(:,i);
index= find(tmps==0);
tmps(index)=[];
Es2 = [Es2,mean(tmps)];
end
Es3 = [Es3,Es2(end)];
end
if jks == 1
Eas = Es3;
save R1.mat Eas SR
end
if jks == 2
Ebs = Es3;
save R2.mat Ebs SR
end
if jks == 3
Ecs = Es3;
save R3.mat Ecs SR
end
end
```
