The main code is as follows, and the commented code is the feasible solution that I have tested.
M = 50;
K = 7;
H_bs = 30;
H_uav = 100;
P_bs_TX = 0.2;
P_bs_sta = 3000;
P_bs_sleep = 600;
P_uav_max_sum = 2;
noise = 10^-14;
N = 120;
S_max = 100;
R_th = 1.5;
Rho = 10^-6;
error = 100;
P_on=P_bs_sta+P_bs_TX*M/K;
P_off=P_bs_sleep;
a=15;
f_a=138.5433;
beta=[0.0087,0.4759,5.4100];
UE_position=[-2074,2004,391,492,-3614,-2618,-2992,554,3003,1132,4206,2994,1473,-2460,429,-2008,1216,4536,-2630,2632,1814,4652,-2218,3681,-4335,-4428,-2148,-1219,3559,4486,877,1179,1958,-4595,3969,-157,-1699,521,-1252,-3157,-1335,4078,2578,1103,-1207,98,4574,3210,2480,-110;-3546,-3926,990,2091,2302,-1249,2807,11,-3040,992,-936,3788,1882,-3566,4700,-4454,4379,-179,839,-1222,-3434,-806,4076,1539,-256,3235,-162,2395,1466,-3236,-4095,-3434,1273,4435,3136,-4072,1702,-2161,1796,2437,2838,-1482,1570,-2726,-1833,-1291,4010,-523,-1844,1282];
BS_position=[0,4000,-4000,2000,-2000,-2000,2000;0,0,0,3464,3464,-3464,-3464];
constraint2=[1,0.5,0.5,0.5,0.5,0.5,0.5];
ua_array=[6,7,1,4,5,3,5,1,7,1,2,4,4,6,4,6,4,2,3,2,7,2,5,2,3,5,3,5,2,7,7,7,4,5,4,6,5,7,5,5,5,2,4,7,6,1,4,2,7,1];
L=[114,463,1065,2039,1989,1863,1191,555,1089,1505,959,1046,1667,472,1999,990,1205,566,1607,1834,191,1037,650,1572,423,2439,1860,1324,1531,2497,1289,822,2192,2771,1997,1940,1788,1971,1829,1548,914,1484,1980,1162,1814,1295,2632,948,1690,1287];
theta = [0,2*pi/(N-1):2*pi/(N-1):2*pi/(N-1)*(N-2),0];
cons2_Circle1 = S_max*(N-1)/3/pi*cos(theta);
cons2_Circle2 = S_max*(N-1)/3/pi*sin(theta);
q_nr=[cons2_Circle1;cons2_Circle2];
%% Feasible solution
% P_uav=P_uav_max_sum/M*ones(M,N);
% q_n=q_nr;
% S_str=ones(1,K);
%% CVX Mosek 9.1.9
cvx_begin
% 变量
variable S_str(1,K) nonnegative
variable q_n(2,N)
variable P_uav(M,N) nonnegative
expression P_tot(1,1)
expression P_tot_1(1,1)
expression P_tot_2(1,1)
expression R_bs(1,M)
expression r_uu_n(2,N)
expression f0_n(1,N)
expression f1(1,N)
expression Taylor(M,N)
expression R_uav(1,M)
% variable V(1,N-1)
% 优化目标
P_tot_1=sum(S_str*P_on+(1-S_str)*P_off)+sum(sum(P_uav))/N;
P_tot_2 = 0;
for i=1:N-1
V = norm((q_n(:,i+1)-q_n(:,i)),2);
P_tot_2 = P_tot_2...
+f_a...
+beta(1)*pow_pos(V-a,3)+beta(2)*square_pos(V-a)+beta(3)*max([V-a,0]);
end
P_tot=P_tot_1+1/(N-1)*P_tot_2;
R_bs=S_str(ua_array).*log(1+P_bs_TX*Rho./L./L/noise);
for i=1:M
r_uu=q_nr-UE_position(:,i);
for j=1:N
r_uu_n(:,j)=q_n(:,j)-UE_position(:,i);
end
f0=r_uu(1,:).^2+r_uu(2,:).^2;
f0_n=r_uu_n(1,:).^2+r_uu_n(2,:).^2;
f1=2*r_uu.*(q_n-q_nr);
Taylor(i,:)=f0+f1(1,:)+f1(2,:);
R_uav(1,i) = 1/N*sum(...
-log(f0+H_uav^2)...
-(f0_n-f0)./(f0+H_uav^2)...
+log(Taylor(i,:)+H_uav^2+P_uav(i,:)*Rho/noise)...
);
end
minimize(P_tot)
subject to
for i = 1:K
if (constraint2(i) == 0)
S_str(i) == 0;
elseif (constraint2(i) == 1)
S_str(i) == 1;
else
0 <= S_str(i) <= 1;
end
end
for i=1:N
sum(P_uav(:,i)) <= P_uav_max_sum;
end
for i=1:2
q_n(i,1) == q_n(i,N);
end
for i = 1:N-1
norm((q_n(:,i+1)-q_n(:,i)),2) <= S_max;
end
R_th*log(2) <= R_bs+R_uav;
cvx_end