Hello, my codes contains five agents. But I found that the rel_gap is far from 0 when I run to agent4 at iteration k=2, the control screen is as follows.
I found that the rel_gap keeps decreasing toward 0 all the time. But the time is long. So I want to know what causes the problem possibly. Please help me. The relative code is as follows.
clear all; %#ok<CLALL>
close all;
clc;
%Quaratic Function Parameter
n=4;m=24;
e_abs = 1e-4;e_rel = 1e-4; MAX_ITER = 500;%迭代次数
t =0.01;nu = 2; v=100;rho = 0.01;%change
MM=85;
X = zeros(n*m,n+1);
X_new = zeros(n*m,n+1);X_ess= zeros(MAX_ITER,n+1);lambda = zeros(n*m,2*(n+1));
h=ones(m,1);h1=zeros(m,1);
PV=[0, 0, 0, 0, 0, 0, 0, 37.46, 51.45, 56.01,58.06,85, 93, 82.98,51.59,29.70, 20.92, 0, 0, 0, 0, 0, 0, 0];
WT=[46.54,41.27,60.5,35,45, 54, 50, 45, 45.5, 40, 35, 0, 5, 4, 3, 6, 44.5, 48.04,40.5, 36.02,45.5,42.42, 44, 45.5];
d=[120.5,120, 150, 110, 130,131,130, 150.2, 280, 320.5,340,350.4,366.6,341.5, 240.2,220.6, 320.6,302.4,320.7, 208,180, 146.2, 119.6,105.3];
PW=PV+WT;Dd=d-PW;
a1=0.02*h;a2=0.0175*h;a3=0.0625*h;a4=0.01*h;
a=[a1;a2;a3];A=diag(a);
A1=diag([a1;h1;h1;h1]); A2=diag([h1;a2;h1;h1]);A3=diag([h1;h1;a3;h1]);A4=diag([h1;h1;h1;a4]);
a5=2*h;a6=1.75*h;a7=1*h;
B1=[a5;h1;h1;h1];B2=[h1;a6;h1;h1];B3=[h1;h1;a7;h1];
B0=[a5;a6;a7];
SOC_0=0.5;SOC_min=0.15;SOC_max=0.85;n_d=0.95;n_c=0.95;
delta_t=1;E_r=1000;M=tril(ones(24));
lower1=50*h;lower2=20*h;lower3=15*h;lower4=0*h;
Upper1=120*h;Upper2=80*h;Upper3=50*h;Upper4=60*h;
lower=[lower1;lower2;lower3;lower4;]; Upper=[Upper1;Upper2;Upper3;Upper4;];
r1L=-20*ones(m-1,1);r2L=-20*ones(m-1,1);r3L=-20*ones(m-1,1); r1U=20*ones(m-1,1);r2U=20*ones(m-1,1);r3U=20*ones(m-1,1);
R1L=[r1L;0];R2L=[r2L;0];R3L=[r3L;0]; R1U=[r1U;0];R2U=[r2U;0];R3U=[r3U;0];
m1=[-1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
m2=zeros(m,m);
M1=[m1,m2,m2,m2];M2=[m2,m1,m2,m2];M3=[m2,m2,m1,m2,];
H1=diag([h;h1;h1;h1]);H2=diag([h1;h;h1;h1]);H3=diag([h1;h1;h;h1]);H4=diag([h1;h1;h1;h]);
for k = 1:MAX_ITER
cvx_solver Mosek
%Agent 1
% x-update
cvx_begin
variable x_1(n*m)
obj1=x_1'*A1*x_1+B1'*x_1;
for i=1:24
cs1=inv_pos(-lower1(1)+H1(i,:)*x_1)+ inv_pos(Upper1(1)-H1(i,:)*x_1);
obj1=obj1+cs1*(1/t);
end
obj2=0;
for j=1:23
rate1=M1*x_1;
cs2=inv_pos(-r1L(1)+rate1(j,:))+ inv_pos(r1U(1)-rate1(j,:));
obj2=obj2+cs2*(1/t);
end
minimize(obj1+obj2+ (lambda(:,1)-lambda(:,10))'*x_1+(lambda(:,2)-lambda(:,3))'*x_1+...
rho* sum_square_abs(X(:,5)- x_1)+rho* sum_square_abs(X(:,2)- x_1))
cvx_end
disp(x_1);
%lambda-update
X_new(:,1) = x_1;
lambda(:,1) = lambda(:,1) -rho*(X(:,5)-X_new(:,1));
lambda(:,2) = lambda(:,2) -rho*(X(:,2)-X_new(:,1));
lambda_result1(n*m*(k-1)+1:n*m*k,1)=lambda(:,1);
lambda_result1(n*m*(k-1)+1:n*m*k,2)=lambda(:,2);
%Agent 2
% x-update
cvx_begin
variable x_2(n*m)
obj1=x_2'*A2*x_2+B2'*x_2;
for i=25:48
cs1=inv_pos(-lower2(1)+H2(i,:)*x_2)+ inv_pos(Upper2(1)-H2(i,:)*x_2);
obj1=obj1+cs1*(1/t);
end
obj2=0;
for j=1:23
rate2=M2*x_2;
cs2=inv_pos(-r2L(1)+rate2(j,:))+ inv_pos(r2U(1)-rate2(j,:));
obj2=obj2+cs2*(1/t);
end
minimize(obj1+obj2 +(lambda(:,3)-lambda(:,2))'*x_2+(lambda(:,4)-lambda(:,5))'*x_2+...
rho* sum_square_abs(X(:,1)- x_2)+rho* sum_square_abs(X(:,3)- x_2))
cvx_end
disp(x_2);
%lambda-update
X_new(:,2) = x_2;
lambda(:,3) = lambda(:,3) -rho*(X(:,1)-X_new(:,2));
lambda(:,4) = lambda(:,4) -rho*(X(:,3)-X_new(:,2));
lambda_result1(n*m*(k-1)+1:n*m*k,3)=lambda(:,3);
lambda_result1(n*m*(k-1)+1:n*m*k,4)=lambda(:,4);
%Agent 3
% x-update
cvx_begin
variable x_3(n*m)
obj1=x_3'*A3*x_3+B3'*x_3;
for i=49:72
cs1=inv_pos(-lower3(1)+H3(i,:)*x_3)+ inv_pos(Upper3(1)-H3(i,:)*x_3);
obj1=obj1+cs1*(1/t);
end
obj2=0;
for j=1:23
rate3=M3*x_3;
cs2=inv_pos(-r3L(1)+rate3(j,:))+ inv_pos(r3U(1)-rate3(j,:));
obj2=obj2+cs2*(1/t);
end
minimize(obj1+obj2 +(lambda(:,5)-lambda(:,4))'*x_3+(lambda(:,6)-lambda(:,7))'*x_3+...
rho* sum_square_abs(X(:,2)- x_3)+rho* sum_square_abs(X(:,4)- x_3))
cvx_end
disp(x_3);
%lambda-update
X_new(:,3) = x_3;
lambda(:,5) = lambda(:,5) -rho*(X(:,2)-X_new(:,3));
lambda(:,6) = lambda(:,6) -rho*(X(:,4)-X_new(:,3));
lambda_result1(n*m*(k-1)+1:n*m*k,5)=lambda(:,5);
lambda_result1(n*m*(k-1)+1:n*m*k,6)=lambda(:,6);
%Agent 4
% x-update
cvx_begin
variable xx(1,m) binary
variable yy(1,m) binary
variables x_41(n*m) x_42(n*m) x_4(n*m)
variables z1(1,m) z2(1,m)
obj3=z1*ones(m,1)+z2*ones(m,1);
obj2=0;
for j=1:24
rate5=-(delta_t/E_r)*M*(1/n_d*z1'-n_c*z2');
cs2=inv_pos(-SOC_min+SOC_0+rate5(j,:))+ inv_pos(SOC_max-SOC_0-rate5(j,:));
obj2=obj2+cs2*(1/t);
end
minimize(obj3+(lambda(:,7)-lambda(:,6))'*x_4...
+(lambda(:,8)-lambda(:,9))'*x_4+rho* sum_square_abs(X(:,3)- x_4)+rho*
sum_square_abs(X(:,5)- x_4))
subject to
for i=1:24
z1(i)<=xx(i)*Upper4(1);
z1(i)<=x_41(72+i);
z1(i)>=x_41(72+i)-(1-xx(i))*Upper4(1)
0<=z1(i)<=Upper4(1);
lower4(1)<=x_41(72+i)<=Upper4(1);
z2(i)<=yy(i)*Upper4(1);
z2(i)<=x_42(72+i);
z2(i)>=x_42(72+i)-(1-yy(i))*Upper4(1)
0<=z2(i)<=Upper4(1);
lower4(1)<=x_42(72+i)<=Upper4(1);
-MM*(1-xx(i))<=x_4(72+i)-x_41(72+i)<=MM*(1-xx(i));
-MM*(1-yy(i))<=x_4(72+i)+x_42(72+i)<=MM*(1-yy(i));
xx(i)+yy(i)==1;
end
cvx_end
disp(x_4);
%lambda-update
X_new(:,4) = x_4;
lambda(:,7) = lambda(:,7) -rho*(X(:,3)-X_new(:,4));
lambda(:,8) = lambda(:,8) -rho*(X(:,5)-X_new(:,4));
lambda_result1(n*m*(k-1)+1:n*m*k,7)=lambda(:,7);
lambda_result1(n*m*(k-1)+1:n*m*k,8)=lambda(:,8);
%Agent 5
% x-update
cvx_begin
variable x_5(n*m)
F1=x_5(1,:)+x_5(24+1,:)+x_5(2*24+1,:)+x_5(3*24+1,:);
F2=x_5(2,:)+x_5(24+2,:)+x_5(2*24+2,:)+x_5(3*24+2,:);
F3=x_5(3,:)+x_5(24+3,:)+x_5(2*24+3,:)+x_5(3*24+3,:);
F4=x_5(4,:)+x_5(24+4,:)+x_5(2*24+4,:)+x_5(3*24+4,:);
F5=x_5(5,:)+x_5(24+5,:)+x_5(2*24+5,:)+x_5(3*24+5,:);
F6=x_5(6,:)+x_5(24+6,:)+x_5(2*24+6,:)+x_5(3*24+6,:);
F7=x_5(7,:)+x_5(24+7,:)+x_5(2*24+7,:)+x_5(3*24+7,:);
F8=x_5(8,:)+x_5(24+8,:)+x_5(2*24+8,:)+x_5(3*24+8,:);
F9=x_5(9,:)+x_5(24+9,:)+x_5(2*24+9,:)+x_5(3*24+9,:);
F10=x_5(10,:)+x_5(24+10,:)+x_5(2*24+10,:)+x_5(3*24+10,:);
F11=x_5(11,:)+x_5(24+11,:)+x_5(2*24+11,:)+x_5(3*24+11,:);
F12=x_5(12,:)+x_5(24+12,:)+x_5(2*24+12,:)+x_5(3*24+12,:);
F13=x_5(13,:)+x_5(24+13,:)+x_5(2*24+13,:)+x_5(3*24+13,:);
F14=x_5(14,:)+x_5(24+14,:)+x_5(2*24+14,:)+x_5(3*24+14,:);
F15=x_5(15,:)+x_5(24+15,:)+x_5(2*24+15,:)+x_5(3*24+15,:);
F16=x_5(16,:)+x_5(24+16,:)+x_5(2*24+16,:)+x_5(3*24+16,:);
F17=x_5(17,:)+x_5(24+17,:)+x_5(2*24+17,:)+x_5(3*24+17,:);
F18=x_5(18,:)+x_5(24+18,:)+x_5(2*24+18,:)+x_5(3*24+18,:);
F19=x_5(19,:)+x_5(24+19,:)+x_5(2*24+19,:)+x_5(3*24+19,:);
F20=x_5(20,:)+x_5(24+20,:)+x_5(2*24+20,:)+x_5(3*24+20,:);
F21=x_5(21,:)+x_5(24+21,:)+x_5(2*24+21,:)+x_5(3*24+21,:);
F22=x_5(22,:)+x_5(24+22,:)+x_5(2*24+22,:)+x_5(3*24+22,:);
F23=x_5(23,:)+x_5(24+23,:)+x_5(2*24+23,:)+x_5(3*24+23,:);
F24=x_5(24,:)+x_5(24+24,:)+x_5(2*24+24,:)+x_5(3*24+24,:);
FF=[F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,
F13,F14,F15,F16,F17,F18,F19,F20,F21,F22,F23,F24];
minimize(v*sum_square_abs(FF-Dd)+(lambda(:,9)-lambda(:,8))'*x_5...
+(lambda(:,10)-lambda(:,1))'*x_5+rho* sum_square_abs(X(:,4)- x_5)+rho* sum_square_abs(X(:,1)- x_5))
cvx_end
disp(x_5);
%lambda-update
X_new(:,5) = x_5;
lambda(:,9) = lambda(:,9) -rho*(X(:,4)-X_new(:,5));
lambda(:,10) = lambda(:,10) -rho*(X(:,1)-X_new(:,5));
lambda_result1(n*m*(k-1)+1:n*m*k,9)=lambda(:,9);
lambda_result1(n*m*(k-1)+1:n*m*k,10)=lambda(:,10);
X = X_new;
for i = 1:n+1
X_result1(n*m*(k-1)+1:n*m*k,i) = X(:,i);
end
result(k)=x_1'*A1*x_1+B1'*x_1+x_2'*A2*x_2+B2'*x_2+x_3'*A3*x_3+B3'*x_3+x_4'*A4*x_4;
disp(result(k));
% Stop Criteria
if( k > 1)
for j = 1 : 5
if( j == 1 )
e_p(j) = norm( X(:,5) - X(:,1) ) + norm( X(:,2) - X(:,1) );
e_pri(j) = e_abs + ( max( norm( X(:,5) ) , max(norm( X(:,1) ) , norm(X(:,2)) ) ) ) * e_rel ;
end
if( j == 5)
e_p(j) = norm( X(:,4) - X(:,5) ) + norm( X(:,1) - X(:,5) );
e_pri(j) = e_abs + ( max( norm( X(:,4) ) , max( norm( X(:,5) ) , norm(X(:,1)) ) ) ) * e_rel ;
end
if(j~=1 && j~=5)
e_p(j) = norm( X(:,j-1) - X(:,j) ) + norm( X(:,j+1) - X(:,j) );
e_pri(j) = e_abs + ( max( norm( X(:,j-1) ) , max( norm( X(:,j) ) , norm(X(:,j+1)) ) ) ) * e_rel ;
end
e_d(j) = rho * norm( X_result1(n*m*(k-2)+1:n*m*(k-1),j) - X(:,j) );
e_dual(j) = e_abs + ( max( norm(lambda(:,2*j-1)) , norm(lambda(:,2*j)) ) ) * e_rel ;
end
if( (sum(e_p <= e_pri) ==5) && ( sum(e_d <= e_dual) == 5 ) )
break
end
end
t = nu*t;
end