How to find BUG in CVX?

I try to reproduce a paper which used SCA(successive convex approximation) algorithm to solved a jointly optimizing UAV position, communication and computing resource allocation, and task splitting decisions problem. But after I imitate their mathematical formulas I can not get the number which the paper solved and the number is too large(The optimal solution obtained by the first iteration is not on the same level as the final optimal solution in the article
), I‘m sure that the value shrinkage in my formula was in the range, so I decided to check my code right or not, then I have a question that when I use the function pow_pos(,2) in variables (in my code like f_UAV) to replace(.^2),the result are different between them, but I don’t know why, all my variables are positive, can you help me?
Here my mathematical formulas and code

**clc
clear all
close all
cvx_clear

k=1e-28;
L=[3;5;2;3;5;1;1;5;4;5].*1e6; %
C=[200;200;200;200;200;200;200;200;200;200];
lamda=[30;30;30;30;30;30;30;30;30;30];
rou=5;
tao_beitaij=1;
tao_faij=1;
tao_lamda=1;
tao_BUL=1;
tao_fUAV=1;
tao_beitai0=1;
P_UAVTX=1;
P_UAVRX=0.1;
P_MU=[0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;];
theta=10^((-100/10)-3);
B_DL=[0.5 0.5 0.5 0.5].*1e6;
B_ULtotal=1e7;
F_UAV=3e9;
F_EC=[8 9 6 7].*1e9;
a0=10^((-50)/10);
beita_i00=[0.2;0.2;0.2;0.2;0.2;0.2;0.2;0.2;0.2;0.2];
f_UAV0=[0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1].*1e9;
beita_ij0=[0.2 0.2 0.2 0.2;0.2 0.2 0.2 0.2;0.2 0.2 0.2 0.2;0.2 0.2 0.2 0.2;…
0.2 0.2 0.2 0.2;0.2 0.2 0.2 0.2;0.2 0.2 0.2 0.2;0.2 0.2 0.2 0.2;0.2 0.2 0.2 0.2;0.2 0.2 0.2 0.2];
f_EC0=[0.1 0.1 0.1 0.1;0.1 0.1 0.1 0.1;0.1 0.1 0.1 0.1;…
0.1 0.1 0.1 0.1;0.1 0.1 0.1 0.1;0.1 0.1 0.1 0.1;…
0.1 0.1 0.1 0.1;0.1 0.1 0.1 0.1;0.1 0.1 0.1 0.1;0.1 0.1 0.1 0.1].*1e9;
langda0=[0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1];
fai_j0=[0.1 0.1 0.1 0.1];
B_UL0=[0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5].*1e6;
Q_UAVx0=5;
Q_UAVy0=7;
Q_UAV0=[Q_UAVx0,Q_UAVy0,100];

Q_EC=[0 0 0;1000 0 0;1000 1000 0;0 1000 0]

r =round( 0 + (500+500)*rand(1,20));
Q_MU=[r(1) r(2) 0;r(3) r(4) 0;r(5) r(6) 0;r(7) r(8) 0;r(9) r(10) 0;…
r(11) r(12) 0;r(13) r(14) 0;r(15) r(16) 0;r(17) r(18) 0;r(19) r(20) 0]

primal_obj_old=-inf;
for jjj=1:100
cvx_begin
cvx_solver mosek
variables beita_i0(10,1)
variables f_UAV(10,1)
variables beita_ij(10,4)
variables fai_j(1,4)
variables f_EC(10,4)
variables B_UL(10,1)
variables langda(10,1)
variables Q_UAVx Q_UAVy
variables z(10,1)
expression loc1(1,4)
loc1=[];
for i=1:4
loc1=[loc1 pow_pos(norm([Q_UAVx Q_UAVy 100]-Q_EC(i,:)),2)];
end
expression loc2(1,4)
loc2=[];
for i=1:4
loc2=[loc2 pow_pos(norm([Q_UAVx0 Q_UAVy0 100]-Q_EC(i,:)),2)];
end
expression loc3(10,1)
loc3=[];
for i=1:10
loc3=[loc3;pow_pos(norm(Q_MU(i,:)-[Q_UAVx Q_UAVy 100]),2)];
end
expression loc4(10,1)
loc4=[];
for i=1:10
loc4=[loc4;pow_pos(norm(Q_MU(i,:)-[Q_UAVx0 Q_UAVy0 100]),2)];
end
minimize (sum(lamda.*(((k.L.C.((beita_i0.(f_UAV0.^2))+beita_i00.*pow_pos(f_UAV,2)))+((tao_beitai0./2).pow_pos(beita_i0-beita_i00,2))+((tao_fUAV./2).pow_pos(f_UAV-f_UAV0,2)))…
+(P_UAVRX.
((L.
((1.*inv_pos(B_UL.*langda0))+(1.*inv_pos(B_UL0.*langda))))+((tao_BUL./2).*pow_pos(B_UL-B_UL0,2))+((tao_lamda./2).pow_pos(langda-langda0,2))))…
+sum((repmat(L,1,4).P_UAVTX.((beita_ij./repmat(fai_j0,10,1))+(beita_ij0.inv_pos(repmat(fai_j,10,1)))))+((tao_beitaij./2).(pow_pos(beita_ij-beita_ij0,2)))+repmat(((tao_faij./2).pow_pos(fai_j-fai_j0,2)),10,1),2)))…
+(rou.
(sum(((L.
((1.*inv_pos(B_UL.*langda0))+(1.*inv_pos(B_UL0.*langda))))+((tao_BUL./2).*pow_pos(B_UL-B_UL0,2))+((tao_lamda./2).*pow_pos(langda-langda0,2)))+z))))
subject to

``````    z>=(L.*C.*((0.5.*(pow_pos((beita_i0+(1.*inv_pos(f_UAV))),2)-(beita_i00.^2)-(1./f_UAV0).^2))-(beita_i00...  %?????
.*(beita_i0-beita_i00))+(pow_pos((1./f_UAV0),3).*((1.*inv_pos(f_UAV))-(1./f_UAV0)))));

repmat(z,1,4)>=(repmat(L,1,4).*(((0.5.*(pow_pos(beita_ij+repmat((1.*inv_pos(fai_j)),10,1),2)-pow_pos(beita_ij0,2)-repmat(pow_pos(1./fai_j0,2),10,1)))-...
(beita_ij0.*(beita_ij-beita_ij0))+(repmat(pow_p(1./fai_j0,3),10,1).*repmat((1.*inv_pos(fai_j))-(1./fai_j0),10,1)))))+...
(repmat(L,1,4).*repmat(C,1,4).*((0.5.*(pow_pos(beita_ij+(1.*inv_pos(f_EC)),2)-pow_pos(beita_ij0,2)-pow_pos(1./f_EC0,2)))...
-(beita_ij0.*(beita_ij-beita_ij0))+((pow_pos(1./f_EC0,3)).*((1.*inv_pos(f_EC))-(1./f_EC0)))));

0<=fai_j<=(B_DL.*log2(1+((a0./loc2).*P_UAVTX./theta))-...
((B_DL.*(a0.*P_UAVTX./theta).*(loc1-loc2))./...
(log(2).*loc2.*((a0.*P_UAVTX./theta)+loc2))));

0<=langda<=B_UL.*log2(1+((a0./loc4).*P_MU./theta))-((...
(a0.*P_MU./theta).*(loc3-loc4))./...
(log(2).*loc4.*((a0.*P_MU./theta)+loc4)));

sum(B_UL,2)<=B_ULtotal;

beita_i0+sum(beita_ij,2)==1;

sum(f_UAV)<=F_UAV;

sum(f_EC)<=F_EC;

0<=beita_ij<=1;

0<=beita_i0<=1;

B_UL>=0;

f_UAV>=0;

f_EC>=0;

0<=Q_UAVx<=1000;

0<=Q_UAVy<=1000;

cvx_end

beita_i00=beita_i00+(beita_i0-beita_i00)./jjj;
f_UAV0=f_UAV0+(f_UAV-f_UAV0)./jjj;
beita_ij0=beita_ij0+(beita_ij-beita_ij0)./jjj;
fai_j0=fai_j0+(fai_j-fai_j0)./jjj;
f_EC0=f_EC0+(f_EC-f_EC0)./jjj;
langda0=langda0+(langda-langda0)./jjj;
B_UL0=B_UL0+(B_UL-B_UL0)./jjj;
Q_UAV0=Q_UAV0+([Q_UAVx Q_UAVy 100]-Q_UAV0)./jjj;

primal_obj = cvx_optval;
``````

cur_err= abs(primal_obj -primal_obj_old);
if primal_obj>primal_obj_old
primal_obj_old=primal_obj;
end

``````if cur_err<10^-2  %terminate condition
break
end
Q_UAV0
``````

end

**

Here is the running result

Q_EC =

``````       0           0           0
1000           0           0
1000        1000           0
0        1000           0
``````

Q_MU =

789 368 0
206 87 0
772 206 0
388 552 0
229 642 0
484 152 0
782 101 0
294 237 0
531 91 0
405 105 0

Calling SDPT3 4.0: 1863 variables, 878 equality constraints For improved efficiency, SDPT3 is solving the dual problem.

num. of constraints = 878
dim. of sdp var = 812, num. of sdp blk = 406
dim. of socp var = 56, num. of socp blk = 14
dim. of linear var = 579
dim. of free var = 10 *** convert ublk to lblk
number of nearly dependent constraints = 10
To remove these constraints, re-run sqlp.m with OPTIONS.rmdepconstr = 1.

SDPT3: Infeasible path-following algorithms

number of iterations = 94 primal objective value = -6.64391289e+09 dual objective value = -6.64651288e+09 gap := trace(XZ) = 1.83e+03 relative gap = 1.37e-07 actual relative gap = 1.96e-04 rel. primal infeas (scaled problem) = 3.65e-14 rel. dual " " " = 2.31e-16 rel. primal infeas (unscaled problem) = 0.00e+00 rel. dual " " " = 0.00e+00 norm(X), norm(y), norm(Z) = 1.3e+10, 6.2e+20, 6.2e+20 norm(A), norm(b), norm© = 1.3e+04, 7.1e+09, 1.5e+10 Total CPU time (secs) = 2.94 CPU time per iteration = 0.03 termination code = -7 DIMACS: 1.7e-13 0.0e+00 4.0e-16 0.0e+00 2.0e-04 1.4e-07

Status: Inaccurate/Solved
Optimal value (cvx_optval): +6.64651e+09

Q_UAV0 =

90.2948 90.3252 100.0000

Calling SDPT3 4.0: 1863 variables, 878 equality constraints For improved efficiency, SDPT3 is solving the dual problem.

num. of constraints = 878
dim. of sdp var = 812, num. of sdp blk = 406
dim. of socp var = 56, num. of socp blk = 14
dim. of linear var = 579
dim. of free var = 10 *** convert ublk to lblk
number of nearly dependent constraints = 10
To remove these constraints, re-run sqlp.m with OPTIONS.rmdepconstr = 1.

SDPT3: Infeasible path-following algorithms

number of iterations = 51 primal objective value = 2.78655377e+16 dual objective value = -4.60049295e+16 gap := trace(XZ) = 1.73e+17 relative gap = 2.34e+00 actual relative gap = 1.00e+00 rel. primal infeas (scaled problem) = 2.39e+00 rel. dual " " " = 6.56e-05 rel. primal infeas (unscaled problem) = 0.00e+00 rel. dual " " " = 0.00e+00 norm(X), norm(y), norm(Z) = 1.3e+15, 8.8e+19, 8.8e+19 norm(A), norm(b), norm© = 1.3e+03, 1.3e+08, 1.5e+10 Total CPU time (secs) = 1.48 CPU time per iteration = 0.03 termination code = -5 DIMACS: 1.3e+01 0.0e+00 1.1e-04 0.0e+00 1.0e+00 2.3e+00

Status: Failed
Optimal value (cvx_optval): NaN

Q_UAV0 =

291.4287 290.9253 100.0000

Disciplined convex programming error:
Illegal operation: {concave} + {convex}

+sum((repmat(L,1,4).P_UAVTX.((beita_ij./repmat(fai_j0,10,1))+(beita_ij0.inv_pos(repmat(fai_j,10,1)))))+((tao_beitaij./2).(pow_pos(beita_ij-beita_ij0,2)))+repmat(((tao_faij./2).*pow_pos(fai_j-fai_j0,2)),10,1),2)))…

If your code is 'correct", which I have no idea whether it is,

Having non-zero input data many orders of magnitude different than 1, as does this input data, can make solver performance very precarious. The solver had difficulty on the 1st iteration and failed on the next.

Crude SCA is inherently unstable, all the more compounded by the poor input data. Its performance and path can strongly depend on the starting point (is yours the same as the paper? does the paper even provide its starting point?), and on solver performance. Even slight differences in the solution in one iteration can cause a totally different path to be taken on subsequent iterations. To expect your SCA solution to match the paper’s might be akin to expecting a miracle.

Thanks for your reply，I went through some posts on the forums that you’ve replied to before and benefit me a lot.
The paper give a part of starting point and set the tone for the entire scale, that is why I think the mathematical formulas has no problem.
And my question is: when I set a constraint that the variables(like f_UAV in my code) are positive, Logically I want get the square of f_UAV by using pow_pos(f_UAV,2) or f_UAV.^2 will get the same result, but in fact it’s not and I don’t know why.

Given the constraint f_UAV >= 0, then pow_pos(f_UAV,2) and f_UAV.^2 should be the same. Nevertheless, the formulation provided to the solver is not exactly the same, even though it would be mathematically equivalent if everything were calculated in exact arithmetic with zero infeasibility tolerance (which it is not). pow_pos uses max(0,f_UAV) which is not equal to f_UAV if f_UAV is even slightly negative. This could cause the result of the solver to differ, even if only slightly, in one iteration, and cause SCA to take a different path. I have no idea if that’s what’s going on here.

Alright, thanks for your patient answer, I’ll go check my code

Please use English on this forum. Thanks.

Sorry, I didn’t notice this. I will pay attention next time.

I am also very interested in the reproduction of this article, can you give me the contact information so I can ask you about it?