This is my code:
clear all
% Define A, B, c0, d0, D_su, D_ud, N, h_min, beta_1, beta_2, N_0, P_u, alpha_su, alpha_ud, P_su, P_ud, and f_t1, f_t2, f_alpha_su, f_alpha_ud, f_tau_su, f_tau_ud as appropriate
% 参数定义
c0 = 0.6; % 环境常数
d0 = 0.11; % 环境常数canhu
LoS = 1; % Line-of-Sight (LoS) 路径损耗
NLoS = 20; % Non-Line-of-Sight (NLoS) 路径损耗
hmin = 100; % UAV的最小高度
L = 1000; % 远程基站的位置
N = 3; % 接入点的数量
APs = [-84.5, -21.6, 0; 250.3, 5.9, 0; -331.5, 62.7, 0]; % 接入点的坐标
UAV = [0, 0, 0]; % UAV 的坐标
dud = [L, 0, 0]; % UAV到远程基站的距离
Pmax = 1; % 最大发射功率
Bt = 1; % 总可用带宽
N_0 = 1e-15; % 噪声功率谱密度
P_s=1;%固定值
beta1 = 0.5; % 接入链路带宽分配系数固定值
beta2 = 0.5; % 回传链路带宽分配系数固定值
P_u=1;%P_u=Pmax
% 给定的参数
rho0 = -80;
Omega_NLoS = 20;
Omega_LoS = 1;
% 计算A和B
A = rho0 * 10^(-0.1 * Omega_NLoS);
B = rho0 * (10^(-0.1 * Omega_LoS) - 10^(-0.1 * Omega_NLoS));
% 计算每个远程基站到 UAV 的距离
dm1 = norm(UAV - APs(1, :));
dm2 = norm(UAV - APs(2, :));
dm3 = norm(UAV - APs(3, :));
dm4 = L; %1000
%迭代
alpha_var1_l=1;alpha_var2_l=1;alpha_var3_l=1;alpha_var4_l=1;
t_l=1;
h_u_l=1;
tau_var1_l=1;tau_var2_l=1;tau_var3_l=1;
cvx_begin
%variables h_u pi_var theta t alpha_tau(length(alpha), length(tau)) s(length(s)) nonnegative;
variables pi_var theta ;
expressions alpha_var(4) tau_var(3) s(4) h_u t ;
maximize(pi_var)
subject to
tau_var(1)=1; tau_var(2)=1; tau_var(3)=1;
alpha_var(1)=1; alpha_var(2)=1; alpha_var(3)=1; alpha_var(4)=1;
s(1)=1; s(2)=1; s(3)=1; s(4)=1;
h_u=1;t=2;
% Additional constraints
h_u >= hmin; % Constraint (6g) ok
pi_var <= Ntheta; % Constraint (8c) ok
%pi_var <= Btbeta2log2(1 + (P_u * alpha_var(4)) / (Btbeta2N_0)); % Constraint (14c) ok
pi_var <= Btbeta2*-rel_entr(1,1 + (P_u * alpha_var(4)) / (Btbeta2N_0)) / log(2); % Constraint (14c) ok
inv_pos(t) <= h_u; % Constraint (18) ok
% Constraint (28c)ok
%(D_su.^2)./(A + B*c0.*(s_s_u - 15).^d0) <= f_alpha_su
(dm1^2+h_u^2)/(A + B*c0*(s(1) - 15)^d0) <= (2/alpha_var1_l-alpha_var(1)/alpha_var1_l^2);%分分开写/注意分母alpha_var应该有l
(dm2^2+h_u^2)/(A + B*c0*(s(2) - 15)^d0) <= (2/alpha_var2_l-alpha_var(2)/alpha_var2_l^2);
(dm3^2+h_u^2)/(A + B*c0*(s(3) - 15)^d0) <= (2/alpha_var3_l-alpha_var(3)/alpha_var3_l^2);
% Constraint (28d)OK
15 <= s(1) & s(1) <= (180/pi_var)*(atan(1/(dm1*t_l))-(dm1/(1+(dm1^2*(t_l)^2)))*(t-t_l));
15 <= s(2) & s(2) <= (180/pi_var)*(atan(1/(dm2*t_l))-(dm2/(1+(dm2^2*(t_l)^2)))*(t-t_l));
15 <= s(3) & s(3) <= (180/pi_var)*(atan(1/(dm3*t_l))-(dm3/(1+(dm3^2*(t_l)^2)))*(t-t_l));
% Constraint (28e)ok
(L^2+h_u^2)./(A + B*c0.*(s(4) - 15).^d0) <= (2/alpha_var4_l-alpha_var(4)/alpha_var4_l^2);
% Constraint (28f)ok
15 <= s(4) & s(4) <= (180/pi_var)*(atan(1/(L*t_l))-(L/(1+(L^2*(t_l)^2)))*(t-t_l));
% Constraint (28g) 注意i=2...N
% log(tau_var(1)) + log(dm1^2+1/(t_l)^2)-(2*(t-t_l))/(dm1^2*(t_l)^3+t_l) >= ...
% log(A+B*c0*((180/pi_var)*atan(h_u_l/dm1) - 15)^d0) + (h_u-h_u_l) * (((180/pi_var)*B*c0*d0*((180/pi_var)*atan(h_u_l/dm1)-15)^(d0-1)*(dm1/(dm1^2+h_u_l^2))) ...
% /(A+B*c0*((180/pi_var)*atan(h_u_l/dm1)-15)^d0));
-rel_entr(1,tau_var(2)) + -rel_entr(1,dm2^2+1/(t_l)^2)-(2*(t-t_l))/(dm2^2*(t_l)^3+t_l) >= ...
-rel_entr(1,A+B*c0*((180/pi_var)*atan(h_u_l/dm2) - 15)^d0) + (h_u-h_u_l) * ((180/pi_var)*B*c0*d0*((180/pi_var)*atan(h_u_l/dm2)-15)^(d0-1)*(dm2/(dm2^2+h_u_l^2))) ...
/(A+B*c0*((180/pi_var)*atan(h_u_l/dm2)-15)^d0);
-rel_entr(1,tau_var(3)) + -rel_entr(1,dm3^2+1/(t_l)^2)-(2*(t-t_l))/(dm3^2*(t_l)^3+t_l) >= ...
-rel_entr(1,A+B*c0*((180/pi_var)*atan(h_u_l/dm3) - 15)^d0) + (h_u-h_u_l) * ((180/pi_var)*B*c0*d0*((180/pi_var)*atan(h_u_l/dm3)-15)^(d0-1)*(dm3/(dm3^2+h_u_l^2))) ...
/(A+B*c0*((180/pi_var)*atan(h_u_l/dm3)-15)^d0);
% Constraint (28h)
theta <= Bt*beta1*-rel_entr(1,sum(P_s*alpha_var(1)+P_s*alpha_var(2)+P_s*alpha_var(3))+Bt*beta1*N_0)/log(2)-Bt*beta1 ...
* -rel_entr(1,sum(P_s*tau_var2_l+P_s*tau_var3_l)+Bt*beta1*N_0)/log(2)+(1/log(2))*(sum(P_s*(tau_var(2)-tau_var2_l)+P_s*(tau_var(3)-tau_var3_l))) ...
/(sum(P_s*tau_var2_l+P_s*tau_var3_l)+Bt*beta1*N_0) ; %i=1
theta <= Bt*beta1*-rel_entr(1,sum(P_s*alpha_var(2)+P_s*alpha_var(3))+Bt*beta1*N_0)/log(2)-Bt*beta1 ...
* -rel_entr(1,sum(P_s*tau_var3_l)+Bt*beta1*N_0)/log(2)+(1/log(2))*(sum(P_s*(tau_var(3)-tau_var3_l))) ...
/(sum(P_s*tau_var3_l)+Bt*beta1*N_0) ; %i=2
theta <= Bt*beta1*-rel_entr(1,sum(P_s*alpha_var(3))+Bt*beta1*N_0)/log(2)-Bt*beta1 ...
* -rel_entr(1,Bt*beta1*N_0)/log(2)+(1/log(2))* 0; %i=3
cvx_end
The result of the run is：
