tne version of matlab is R2016a,the cvx is SDPT3.
H =100; %% height
V = 3; %% speed
noise_pw = 10^(-14);
reference = 10^(-6);
gamma=10.^8;
u0=[100;600];
uN=[100;-600];
xs=[0;0];
xe=[200;0];
N=799; % T=200s, dt=0.5s;
i=1:N;
Q1_r=zeros(2,N);
Q1_r=[u0(1).*ones(1,N);u0(2)+(uN(2)-u0(2))/(N-1).*(i-1)];
Q2_r=zeros(2,N);
Q2_r=[u0(1).*ones(1,N);u0(2)+(uN(2)-u0(2))/(N-1).*(i-1)];
phi=3; %%
avg_Ps=10^(1)/1000;
avg_Pj=10^(0)/1000; %%%%%平均功率
Ps_r=avg_Psones(1,N);
Pj_r=avg_Pjones(1,N);
Ps_max=4avg_Ps;
Pj_max=4avg_Pj;
dt=0.5;
BisecConvergenceThreshold = 10^(-4);
N_num = 0;
disp(’%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%’)
disp(['the number of timeslots is ',num2str(N)]);
t_traj = 0;
ll=0;kk = 0;
while (1)
ll=ll+1;
disp(’%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%’)
disp(['the UAV1 power control ',num2str(ll)]);
for i = 1:N
a(i) = gamma./(H^2+norm(Q1_r(:,i)-xs).^2)./(gamma.*Pj_r(i)./(H^2+norm(Q2_r(:,i)-xs).^2)+1); %%an
b(i) = gamma./(H^2+norm(Q1_r(:,i)-xe).^2)./(gamma.*Pj_r(i)./(H^2+norm(Q2_r(:,i)-xe).^2)+1);
end
%first set all power to be zero, and then only allocate power to Splus
Ps = zeros(1,N);
%find the set in which power should be positive
Splus = find(a>b);
if isempty(Splus) %if the set is empty, uniform power allocation
Ps = Ps_r;
else
%only allocate power in Splus
ap = a(Splus);
bp = b(Splus);
%bisection search for lambda
%the lower bound and upper bound of the water level lambda
lambda_lb = 0;
lambda_ub = max(a-b);
while abs(lambda_ub-lambda_lb)/abs(lambda_ub) > BisecConvergenceThreshold
lambda = ( lambda_lb + lambda_ub )/2;
Ptemp = -(2*bp).^(-1) - (2*ap).^(-1) + sqrt( ( (2*bp).^(-1) - (2*ap).^(-1) ).^2 ...,
+ (2/lambda)*( (2*bp).^(-1) - (2*ap).^(-1) ) );
Pplus = min( max( Ptemp , 0 ) , Ps_max );
%compare the sum power, change the value of upper and lower bounds
%Fact: lower lambda brings larger sum power; larger lambda brings lower power
if sum(Pplus) > N*avg_Ps
lambda_lb = lambda;
else
lambda_ub = lambda;
end
end %while
%only allocate power to the subcarriers in Splus
Ps(Splus) = Pplus;
end %if
Ps_r=Ps;
ratej = zeros(1,N);
Pj = zeros(1,N);
for i=1:N
au(i) = gamma.*Ps_r(i)./(H^2+norm(Q1_r(:,i)-xs).^2);
bu(i) = gamma./(H^2+norm(Q2_r(:,i)-xs).^2);
cu(i) = gamma.*Ps_r(i)./(H^2+norm(Q1_r(:,i)-xe).^2);
du(i) = gamma./(H^2+norm(Q2_r(:,i)-xe).^2);
end
A=zeros(1,N);
cvx_solver sdpt3
cvx_save_prefs
cvx_begin
variable Pj(1,N) %%
expression ratej(1,N)
A = au.*bu./(bu.*Pj_r+1)./(bu.*Pj_r + au + 1);
ratej = log(1-cu.*inv_pos(du.*Pj+cu+1))-A.*Pj; % -log(1+temp_aa/(temp_bb*t+temp_cc)
maximize ( sum(ratej) )
subject to
0 <= Pj; %%
Pj <= Pj_max.*ones(1,N);
sum(Pj)<= avg_Pj*N;
%Pu(Ps==0) == 0;
cvx_end
Pj_r=Pj; %%迭代
rate1 = zeros(1,N);
%disp(’%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%’)
disp(['the UAV1 trajectory ',num2str(ll)]);
for i=1:N
l_r(:,i)=H^2+norm(Q1_r(:,i)-xs).^2;
m_r(:,i)=H^2+norm(Q1_r(:,i)-xe).^2;
a1(i) = gamma.*Ps_r(i)./(gamma.*Pj_r(i)./(H^2+norm(Q2_r(:,i)-xs).^2)+1);%% gn
b1(i) = gamma.*Ps_r(i)./(gamma.*Pj_r(i)./(H^2+norm(Q2_r(:,i)-xe).^2)+1);
end
cvx_begin
variable l(1,N) %% (x1[n],y1[n]) 2D location, UAV1’s position sequence across time slots
variable m(1,N)
variable Q1(2,N)
expression rate1(1,N)
B = a1./(l_r.^2+l_r.*a1);
rate1=log(1-b1.*inv_pos(m+b1))-B.*l;
maximize sum(rate1)
subject to
%%% speed constraints
for i = 1:N-1
norm( Q1(:,i+1) - Q1(:,i) ) <= V*dt;
end
for i=1:N
%m(i)+Q1_r(1,i)^2-2*Q1_r(1,i)*Q1(1,i)+2*xe(1)*Q1(1,i)-xe(1)^2+Q1_r(2,i)^2-2*Q1_r(2,i)*Q1(2,i)-H^2<=0;
m(i)+pow_pos(norm(Q1_r(:,i)),2)-2*(Q1_r(:,i)-xe(:,1))'*Q1(:,i)-pow_pos(norm(xe(:,1)),2)-H^2<=0;
% Q1(1,i)^2+Q1(2,i)^2+H^2-l(i)<=0;
pow_pos(norm(Q1(:,i)),2)+H^2-l(i)<=0;
end
0<= m;
Q1(:,1) == u0; %% 1a
Q1(:,N) == uN;%% 1b
cvx_end
Q1_r = Q1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
the result:Status: Infeasible
Optimal value (cvx_optval): -Inf
other result: Solved
Thanks.