67 views

Skip to first unread message

Jun 21, 2022, 11:13:09 AM6/21/22

to YALMIP

L0 = 1;

for iter=1:inter_iter_num

[V1_ini,~] = eig(F1_ini);

[V2_ini,~] = eig(F2_ini);

lambda_max_F1 = V1_ini(:,end);

lambda_max_F2 = V2_ini(:,end);

yalmip('clear');

options = sdpsettings('solver','mosek','verbose',0,'debug',0);

F1 = sdpvar(bs_antenna_number,bs_antenna_number,'hermitian','complex');

F2 = sdpvar(bs_antenna_number,bs_antenna_number,'hermitian','complex');

p = sdpvar(2,1);

epsilon = sdpvar(1);

eta = sdpvar(1);

theta = sdpvar(1);

v = sdpvar(1);

tau = sdpvar(1);

mu = sdpvar(1);

L = sdpvar(1);

omega = sdpvar(1);

R1 = Phi_ini*H_bs_irs*((2^rate_user_th-1)*F1-F2)*H_bs_irs'*Phi_ini';

R2 = Phi_ini*H_bs_irs*(F1-F2)*H_bs_irs'*Phi_ini';

R3 = -Phi_ini*H_bs_irs*F1*H_bs_irs'*Phi_ini';

R4 = Phi_ini*H_bs_irs*F2*H_bs_irs'*Phi_ini';

R5= Phi_ini*H_bs_irs*((2^rate_eve_th-1)*F1-F2)*H_bs_irs'*Phi_ini';

P = [((2^rate_eve_th-1)*noise_power-epsilon)*eye(eve_antenna_num),zeros(eve_antenna_num,irs_antenna_num);

zeros(eve_antenna_num,irs_antenna_num)', epsilon*epsilon_e^(-2)*eye(irs_antenna_num)];

S = [H_irs_e1.',eye(irs_antenna_num)].';

A=[eye(irs_antenna_num), zeros(irs_antenna_num,1);

zeros(1,irs_antenna_num), -epsilon_u^2];

B1=[R1,R1*h_irs_u2; (R1*h_irs_u2)',h_irs_u2'*R1*h_irs_u2+(2^rate_user_th-1)*noise_power];

B2 = [R2,R2*h_irs_u2; (R2*h_irs_u2)',h_irs_u2'*R2*h_irs_u2];

B3 = [R3,R3*h_irs_u2; (R3*h_irs_u2)',h_irs_u2'*R3*h_irs_u2+mu-1];

B4 = [R4,R4*h_irs_u2; (R4*h_irs_u2)',h_irs_u2'*R4*h_irs_u2-2*L0*L+L0^2];

C=[mu,L;L,omega];

Constraints = [F1>=0,F2>=0,p(:)>=0,sum(p) <= Pmax, bs_antenna_number*diag(F1) == p(1), bs_antenna_number*diag(F2) == p(2),...

epsilon>=0,eta>=0,theta>=0,v>=0,tau>=0,mu>=0, omega>=0,...

2*real(trace(M1*F2))>=(real(trace(M1*F1))+noise_power)^2+omega^2, real(trace(M1*F2))>=real(trace(M1*F1)),...

C>=0, eta*A-B1>=0, theta*A-B2>=0,v*A-B3>=0,tau*A-B4>=0, P + S*R5*S'>=0];

N1 = -log(noise_power+trace(M1*F1))/log(2)+1/2/penalty*norm(F1,'nuclear')+1/2/penalty*norm(F2,'nuclear');

D0 = 1/2/penalty*(norm(F1_ini,2)+norm(F2_ini,2));

G_D1_F1 = 1/2/penalty*(lambda_max_F1*lambda_max_F1'); % gradient of D1 with respect to F1

G_D1_F2 = 1/2/penalty*(lambda_max_F2*lambda_max_F2'); % gradient of D1 with respect to F2

Objective =N1-real(D0)-real(trace(G_D1_F1'*(F1-F1_ini)))-real(trace(G_D1_F2'*(F2-F2_ini)));

sol = optimize(Constraints, Objective,options);

if sol.problem==0

F1= double(F1);

F2 =double(F2);

p=double(p);

[V1,D1] = eig(F1);

[V2,D2] = eig(F2);

if (1-max(diag(D1))/p(1)<=delta)&&((1- max(diag(D2))/p(2))*(p(2)>1e-6)<=delta)

flag=1;

f1_hat = V1(:,end)*D1(end,end)^(1/2);

f2_hat = V2(:,end)*D2(end,end)^(1/2);

u1_d = abs(h_irs_u1'*Phi_ini*H_bs_irs*f1_hat)^2;

rate_u1 = log2(1+u1_d/noise_power);

rate=max(real(rate_u1),0);

rate_obj = -value(Objective);

break;

end

F1_ini = F1;

F2_ini = F2;

L0=value(L);

else

%cvx failed

break;

end

end

for iter=1:inter_iter_num

[V1_ini,~] = eig(F1_ini);

[V2_ini,~] = eig(F2_ini);

lambda_max_F1 = V1_ini(:,end);

lambda_max_F2 = V2_ini(:,end);

yalmip('clear');

options = sdpsettings('solver','mosek','verbose',0,'debug',0);

F1 = sdpvar(bs_antenna_number,bs_antenna_number,'hermitian','complex');

F2 = sdpvar(bs_antenna_number,bs_antenna_number,'hermitian','complex');

p = sdpvar(2,1);

epsilon = sdpvar(1);

eta = sdpvar(1);

theta = sdpvar(1);

v = sdpvar(1);

tau = sdpvar(1);

mu = sdpvar(1);

L = sdpvar(1);

omega = sdpvar(1);

R1 = Phi_ini*H_bs_irs*((2^rate_user_th-1)*F1-F2)*H_bs_irs'*Phi_ini';

R2 = Phi_ini*H_bs_irs*(F1-F2)*H_bs_irs'*Phi_ini';

R3 = -Phi_ini*H_bs_irs*F1*H_bs_irs'*Phi_ini';

R4 = Phi_ini*H_bs_irs*F2*H_bs_irs'*Phi_ini';

R5= Phi_ini*H_bs_irs*((2^rate_eve_th-1)*F1-F2)*H_bs_irs'*Phi_ini';

P = [((2^rate_eve_th-1)*noise_power-epsilon)*eye(eve_antenna_num),zeros(eve_antenna_num,irs_antenna_num);

zeros(eve_antenna_num,irs_antenna_num)', epsilon*epsilon_e^(-2)*eye(irs_antenna_num)];

S = [H_irs_e1.',eye(irs_antenna_num)].';

A=[eye(irs_antenna_num), zeros(irs_antenna_num,1);

zeros(1,irs_antenna_num), -epsilon_u^2];

B1=[R1,R1*h_irs_u2; (R1*h_irs_u2)',h_irs_u2'*R1*h_irs_u2+(2^rate_user_th-1)*noise_power];

B2 = [R2,R2*h_irs_u2; (R2*h_irs_u2)',h_irs_u2'*R2*h_irs_u2];

B3 = [R3,R3*h_irs_u2; (R3*h_irs_u2)',h_irs_u2'*R3*h_irs_u2+mu-1];

B4 = [R4,R4*h_irs_u2; (R4*h_irs_u2)',h_irs_u2'*R4*h_irs_u2-2*L0*L+L0^2];

C=[mu,L;L,omega];

Constraints = [F1>=0,F2>=0,p(:)>=0,sum(p) <= Pmax, bs_antenna_number*diag(F1) == p(1), bs_antenna_number*diag(F2) == p(2),...

epsilon>=0,eta>=0,theta>=0,v>=0,tau>=0,mu>=0, omega>=0,...

2*real(trace(M1*F2))>=(real(trace(M1*F1))+noise_power)^2+omega^2, real(trace(M1*F2))>=real(trace(M1*F1)),...

C>=0, eta*A-B1>=0, theta*A-B2>=0,v*A-B3>=0,tau*A-B4>=0, P + S*R5*S'>=0];

N1 = -log(noise_power+trace(M1*F1))/log(2)+1/2/penalty*norm(F1,'nuclear')+1/2/penalty*norm(F2,'nuclear');

D0 = 1/2/penalty*(norm(F1_ini,2)+norm(F2_ini,2));

G_D1_F1 = 1/2/penalty*(lambda_max_F1*lambda_max_F1'); % gradient of D1 with respect to F1

G_D1_F2 = 1/2/penalty*(lambda_max_F2*lambda_max_F2'); % gradient of D1 with respect to F2

Objective =N1-real(D0)-real(trace(G_D1_F1'*(F1-F1_ini)))-real(trace(G_D1_F2'*(F2-F2_ini)));

sol = optimize(Constraints, Objective,options);

if sol.problem==0

F1= double(F1);

F2 =double(F2);

p=double(p);

[V1,D1] = eig(F1);

[V2,D2] = eig(F2);

if (1-max(diag(D1))/p(1)<=delta)&&((1- max(diag(D2))/p(2))*(p(2)>1e-6)<=delta)

flag=1;

f1_hat = V1(:,end)*D1(end,end)^(1/2);

f2_hat = V2(:,end)*D2(end,end)^(1/2);

u1_d = abs(h_irs_u1'*Phi_ini*H_bs_irs*f1_hat)^2;

rate_u1 = log2(1+u1_d/noise_power);

rate=max(real(rate_u1),0);

rate_obj = -value(Objective);

break;

end

F1_ini = F1;

F2_ini = F2;

L0=value(L);

else

%cvx failed

break;

end

end

i have found the warning is due to the constraint *2*real(trace(M1*F2))>=(real(trace(M1*F1))+noise_power)^2+omega^2*, however, the constraint is convex.

As you can see from the figure, there is no error in the first several iterations. Is there exists numerical problem?

Jun 21, 2022, 11:20:09 AM6/21/22

to YALMIP

you would have to supply a minimal reproducible example

Jun 21, 2022, 11:21:04 AM6/21/22

to YALMIP

can i upload .m file or .mat file?

Message has been deleted

Jun 21, 2022, 1:34:14 PM6/21/22

to YALMIP

the model is nasty as the quadratic is low-rank with somewhat bad data, meaning very small floating-point errors cause the matri to be indefinite despite being trivially psd in theory

best is to manually rewrite the quadratic expression to an expression involving a norm instead, i.e. basically derive the socp

K>> eig(quaddecomp((real(trace(M1*F1))+noise_power)^2+omega^2))

tisdag 21 juni 2022 kl. 17:29:37 UTC+2 skrev jifa...@yeah.net:

This is the minimal reproducible example, can you help me?

On Tuesday, June 21, 2022 at 11:20:09 PM UTC+8 Johan Löfberg wrote:

Message has been deleted

Jun 21, 2022, 10:32:07 PM6/21/22

to YALMIP

i have replace the** SOCP** with **2*real(trace(M1*F2))>=norm([(real(trace(M1*F1))+noise_power), omega],2)^2,**but the non-convex quadratic constraint warning still occurs. can you help me?

Jun 22, 2022, 1:18:22 AM6/22/22

to YALMIP

You haven't changed anything

>> x = sdpvar(2,1);

>> x'*x-norm(x)^2

ans =

0

>> x'*x-norm(x)^2

ans =

0

You have to write x'*x<=t as norm([2*x;1-t])<=1+t

Message has been deleted

Jun 22, 2022, 3:02:36 AM6/22/22

to YALMIP

ok, thank you. i have got it. I have replaced the quadratic constraint as** E = [2*real(trace(M1*F1)); 2*omega; sqrt(2)*(1-real(trace(M1*F2)))]; norm(E,2)<=sqrt(2)*(1+real(trace(M1*F2)))**. there is no mosek non-convex quadratic constraint warning when i run my code. however, the yalmip always fail to find a solution, can you help me?

Jun 22, 2022, 3:26:26 AM6/22/22

to YALMIP

Don't see where that sqrt(2) came from. x'*x <= 2y is norm([2*x;1-2y])<=1+2*y

But YALMIP never fails to find a solution as YALMIP never solves any problem. If mosek fails to find a solution, well then you are setting up an infeasible problem.

It looks like you are implementing some kind of iterative scheme, and these can sometimes be very fragile as you assume stuff about the solution iterates (such as the previous solution being returned by the solver without any issues implying that the solution actually is feasible. Since solvers work with infeasible methods, solutions can be slightly infeasible, and if that is a problem for the algorithm, you have to have strategies to combat this)

Jun 22, 2022, 3:55:42 AM6/22/22

to YALMIP

ok, i have got it. thank you very much.

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu