Hi,mcg.
I have written programs using CVX ans YALMIP, respectively, to solve the SDP. Runing the programs with the same parametres, I got different optimal objective results.
Especially, sometimes it is not feasible in YALMIP, but it is feasible in CVX, and the optimal value if strange, i.e., too large.
I commented the codes that generate the matrices, because I loaded the variable use the command ‘load HG.mat’, so they have the same H and G.
I try the case where H(:,:,134) and G(:,:,134) are used, the objective value gained using CVX is very large compared with other H and G, but it is not feasible if I use YALMIIP.
Why? Thanks. Now I don’t know which solver interface to use.
the CVX version
%clear;
num_CSI = 200;
N = 4;
M = 3;
K = 5;
gamma_b = 8;
gamma_b = 10 ^ (gamma_b / 10);
gamma_e = -3;
gamma_e = 10 ^ (gamma_e / 10);
rho_b = 0.1;
rho_e = 0.01;
alpha_b = -log(rho_b);
alpha_e = -log(rho_e);
sigma2 = 0.1;
epsilon2 = 0.002;
load HG.mat % load the variables H and G
% j = 1;
% while (j <= num_CSI)
% H(:,:,j) = (randn(N,M) + sqrt(-1) * randn(N,M)) / sqrt(2); %
% G(:,:,j) = (randn(N,K) + sqrt(-1) * randn(N,K)) / sqrt(2); %
% j = j + 1;
% end
cvx_solver sedumi;%sdpt3;%
for j = 134:134 %1:num_CSI %
cvx_begin sdp quiet
variable W(N,N,M) hermitian;
variable P_z;
variables mu_b(M) nu_b(M);
variables mu_e(M,K) nu_e(M,K);
expression sum_W;
expression U(N,N);
expression D_b(N,N);
expression d_b(N);
expression c_b;
expression D_e(N,N);
expression d_e(N);
expression c_e;
sum_W = 0;
for m = 1:M
sum_W = sum_W + W(:,:,m);
end
V = null(H(:,:,j)');
U = (P_z / (N - M)) * V * V';
minimize(trace(sum_W) + P_z); % object function
subject to
for m = 1:M
% the constraints at the legitimate users
D_b = epsilon2 * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W - U);
D_b = 0.5 * ( D_b + D_b' );%sym(D_b);
d_b = sqrt(epsilon2) * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W) * H(:,m,j);
c_b = sigma2 - real(H(:,m,j)' * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W) * H(:,m,j));
real(trace(D_b)) - sqrt(2 * alpha_b) * mu_b(m) - alpha_b * nu_b(m) - c_b >= 0;
norm(vec([D_b , sqrt(2) * d_b])) <= mu_b(m);
nu_b(m) * eye(N) + D_b >= 0;
nu_b(m) >= 0;
W(:,:,m) >= 0;
% the constraints at the eavesdroppers
for k = 1:K
D_e = epsilon2 * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U);
D_e = 0.5 * ( D_e + D_e' );%sym(D_e);
d_e = sqrt(epsilon2) * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U) * G(:,k,j);
c_e = sigma2 - real(G(:,k,j)' * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U) * G(:,k,j));
real(trace(D_e)) + sqrt(2 * alpha_e) * mu_e(m,k) + alpha_e * nu_e(m,k) - c_e <= 0;
norm(vec([D_e , sqrt(2) * d_e])) <= mu_e(m,k);
nu_e(m,k) * eye(N) - D_e >= 0;
nu_e(m,k) >= 0;
end
end
P_z >= 0;
cvx_end
if strfind(cvx_status,'Solved')
% get the power value
cvx_optval
else
unfeasible = 1
end
end
the YALMIP version
%clear;
num_CSI = 2000; % number of channel realization
N = 4;
M = 3;
K = 5;
gamma_b = 8;
gamma_b = 10 ^ (gamma_b / 10);
gamma_e = -3;
gamma_e = 10 ^ (gamma_e / 10);
rho_b = 0.1;
alpha_b = -log(rho_b);
rho_e = 0.01;
alpha_e = -log(rho_e);
sigma2 = 0.1;
epsilon2 = 0.002;
%generate the estimated channel vectors
%random and independent,standard complex Gaussian distribution CN(0,I)
load HG.mat % load the variables H and G
% j = 1;
% while (j <= num_CSI)
% H(:,:,j) = (randn(N,M) + sqrt(-1) * randn(N,M)) / sqrt(2); %
% %G(:,:,j) = (randn(N,K) + sqrt(-1) * randn(N,K)) / sqrt(2); %
% j = j + 1;
% end
for j = 134:134 %1:num_CSI
mu_b = sdpvar(M,1);
nu_b = sdpvar(M,1);
mu_e = sdpvar(M,K);
nu_e = sdpvar(M,K);
P_z = sdpvar(1,1);
W = sdpvar(N,N,M,'hermitian','complex');
sum_W = zeros(N,N);
for m = 1:M
sum_W = sum_W + W(:,:,m);
end
V = null(H(:,:,j)');
U = (P_z / (N - M)) * V * V';
C = [ ];
for m = 1:M
D_b = epsilon2 * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W - U);
d_b = sqrt(epsilon2) * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W) * H(:,m,j);
c_b = sigma2 - real(H(:,m,j)' * ((1 + 1 / gamma_b) * W(:,:,m) - sum_W) * H(:,m,j));
C = C + [real(trace(D_b)) - sqrt(2 * alpha_b) * mu_b(m) - alpha_b * nu_b(m) - c_b >= 0];
C = C + [norm(vec([D_b , sqrt(2) * d_b])) <= mu_b(m)];
C = C + [nu_b(m) * eye(N) + D_b >= 0];
C = C + [nu_b(m) >= 0];
C = C + [W(:,:,m) >= 0];
for k = 1:K
D_e = epsilon2 * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U);
d_e = sqrt(epsilon2) * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U) * G(:,k,j);
c_e = sigma2 - real(G(:,k,j)' * ((1 + 1 / gamma_e) * W(:,:,m) - sum_W - U) * G(:,k,j));
C = C + [real(trace(D_e)) + sqrt(2 * alpha_e) * mu_e(m,k) + alpha_e * nu_e(m,k) - c_e <= 0];
C = C + [norm(vec([D_e , sqrt(2) * d_e])) <= mu_e(m,k)];
C = C + [nu_e(m,k) * eye(N) - D_e >= 0];
C = C + [nu_e(m,k) >= 0];
end
end
C = C + [P_z >= 0];
Obj = trace(sum_W + P_z);
ops = sdpsettings('solver','sedumi','verbose',0);%,'debug',1sdpt3
diagnostics = solvesdp(C,Obj,ops) ;
if (diagnostics.problem == 0) % if feasible
% get the power value
pow = trace(double(Obj))
else
unfeasible = 1
end
end