when i use yalmip to solve this convex problem, i meet a mosek warning, can you help me?

[jifa zhang]

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

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?

[Johan Löfberg]

you would have to supply a minimal reproducible example

[jifa zhang]

can i upload .m file or .mat file?

[Johan Löfberg]

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:

[jifa zhang]

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?

[Johan Löfberg]

You haven't changed anything

>> x = sdpvar(2,1);
>> x'*x-norm(x)^2

ans =

     0

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

[jifa zhang]

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? 

[Johan Löfberg]

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)

[jifa zhang]

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