MATLAB code is as follows:
clear;clc;close all;cvx_expert true;
%--------------------------------------------------------------------------
%% Sim. Parameter Setting(Scene settings)
qs = [ 0 100 ;
0 0 ;
0 0]; % meters, location of sensors
qd = [1000; 0; 0]; % location of destination
qu = [-800; 0;100]; % initial location of UAV
qv = [+800; 0;100]; % final location of UAV
h = 100; % the flight altitude of UAV
Vmax = 25; % the maximum speed of UAV
N = 3; % Number of transmissions
K = 2; % sensor number
S_nk = [1e6 1e6 ;
2e6 1e6 ;
1e6 1e6 ]; % bit, bit number per packet
B = 1e6; % Hz, total system bandwidth
gamma0 = 10^7.3; % Reference SNR gap (73 dB)
IterMax = 3; % max number of SCA iterations
% Teps = 0.001; % precision of Tmax
Ek = [1 1 ]; % Joule, energy available at each sensor
Eu = 2; % Joule, energy available at UAV
% eps = 0.1;
% Nset = 20 : 10 : 100;
%--------------------------------------------------------------------------
%% Feasibility Check
D_nk = S_nk/Blog(2);
for n=1:3
for k=1:2
if (sum(D_nk)>Ekgamma0/(K*(h^2)))|(sum(sum(D_nk))>Eugamma0/(h^2))
fprintf(‘Infeasible\n’);
return;
end
end
end
%% Initialization if feasible
% step 1: Initialize UAV trajectory
% q1 = zeros(3,NK);
% q1(3,
= h;
% q1(1:2,1:NK) = repmat(q_sk(1:2,:),1,NK);
q1 = zeros(3,NK);
q1(3, = h;
q1(1:2,1:NK) = repmat(qs(1:2,:),1,N);
% q1 = zeros(2,N,K);
% q2 = zeros(2,N);
% for n = 1 : N
% for m = 1 : K
% q1(1:2,n,m) = qs(1:2,m);
% end
% % q2(1:2,n) = qd(1:2,1);
% end
q2 = zeros(3,N);
q2(3, = h;
q2(1:2,1:N) = repmat(qd(1:2),1,N);
% step 2: Initialize energy allocation of Sensor and UAV
Bnk = zeros(N,K);
Bnk(:, = B/K; %Initialize evenly divided bandwidth
E1 = zeros(N,K,K);
% E_nk1l(:, = Ek/N;
for k = 1:K
for n = 1:N
for m = 1:K
E1(n,m,k) = Ek(k)/NK ;
end
end
end
E2 = zeros(N,1);
for n = 1:N
E2(n) = Eu/N;
end
x1 = zeros(N,K,K);
x1 = gamma0E1B/(h^2);
x2 = zeros(N,1);
x2 = gamma0E2/(h^2);
cnk1m = zeros(N,K,K);
cnk1mold = sqrt(gamma0E1B)/(h^2);
c2 = zeros(N,1);
c2old = sqrt(gamma0*E2)/(h^2);
cnk2m = zeros(N,K,K);
b = zeros(N,K,K);
T1 = [65 2;
10 2;
10 2];
% T2 = [70 70 70];
for k = 1:K
for n = 1:N
for m = 1:K
b(n,m,k) = Bnk(n,k)*T1(n,m);
end
end
end
for k = 1:K
for n = 1:N
for m = 1:K
cnk2mold(n,m,k) = sqrt(T1(n,m))/b(n,m,k);
% cnk2mold(n,m,k) = 10;
end
end
end
for k = 1:K
plot(qs(1,k),qs(2,k),'rs','MarkerFaceColor','r','MarkerSize',10);hold on;grid on;
end
plot(qd(1),qd(2),'bs','MarkerFaceColor','b','MarkerSize',8);
plot(qu(1),qu(2),'ko','MarkerFaceColor','k','MarkerSize',8);
plot(qv(1),qv(2),'ko','MarkerFaceColor','k','MarkerSize',8);
%% SCA Iteration
Tmaxold = inf;
iter = 0;
while (iter<=IterMax)
iter = iter + 1;
cvx_begin
variable Tmax
variable T1(N,K)
variable T2(N,1)
variable E1(N,K,K) nonnegative
variable E2(N,1) nonnegative
variable q1(2,N*K)
variable q2(2,N)
variable x1(N,K,K)
variable x2(N,1)
variable Bnk(N,K)
% variable c_nk1l(N,K,K)
% variable c_nk2l(N,K,K)
% variable c2(N,1)
variable b(N,K,K)
minimize(Tmax)
subject to
for k =1 : K
for n = 1 : N
for m =1 : K
cnk1mold(n,m,k)^2*norm([q1(1:2,n*m); h]-qs(k))-2*cnk1mold(n,m,k)*sqrt(gamma0*E1(n,m,k)*B)<=-x1(n,m,k);
cnk2mold(n,m,k)^2*b(n,m,k)-2*cnk2mold(n,m,k)*sqrt(T1(n,m))<=-inv_pos(Bnk(n,k));
end
end
end
for k = 1:K
for n = 1:N
sum(-rel_entr(b(n,:,k), b(n,:,k)+x1(n,:,k)))>=S_nk(n,k)*log(2); %%%1/8 11:50
end
end
for n = 1:N
B=sum(Bnk(n,:));
norm([q2(1:2,n); h]-[q1(1:2,n*K); h])<=Vmax/2*(T1(n,K)+T2(n));
c2old(n)^2*([q2(1:2,n); h]-qd)-2*c2old(n)*sqrt(gamma0*E2(n))<=-x2(n);
-rel_entr(T2(n), T2(n)+x2(n))>=sum(D_nk(n,:));
end
for n =1:N-1
sum(T1(n,:))+T2(n) + sum(T1(n+1,:))+T2(n+1)<=Tmax;
norm([q1(1:2,(n+1)); h]-[q2(1:2,n); h])<=Vmax/2*(T1((n+1),1)+T2(n));
end
norm([q1(1:2,1); h]-qu)<=Vmax/2*T1(1,1);
for n = 1:N
for m = 2:K
norm([q1(1:2,n*m); h]-[q1(1:2,n*(m-1)); h])<=Vmax/2*(T1(n,m)+T1(n,m-1));
end
end
for k = 1:K
sum(E1(:,:,k))<= Ek(k);
end
sum(E2(:))<=Eu;
norm([q2(1:2,N); h]-qv)<=Vmax/2*T2(N);
% T1(n)+T2(n) + T1(n+1)+T2(n+1)<=Tmax;
% Tmax==(T1(1)+T2(1) + T1(N)+T2(N) + 2*(sum(T1(2:N-1))+sum(T2(2:N-1))))/(N-1);
% q1(:,1) == qu(1:2);
% q2(:,N) == qv(1:2);
% sum(Bn) = B;
% sum(sum(E_nk1l))<=Ek;
% sum(E2)<=Eu;
cvx_end
% cnk1mold = sqrt(gamma0E1B)/(h^2);
% for k = 1:K
% for n = 1:N
% for m = 1:K
% b(n,m,k) = Bnk(n,k)T1(n,m);
% end
% end
% end
% for k = 1:K
% for n = 1:N
% for m = 1:K
% cnk2mold(n,m,k) = sqrt(T1(n,m))/b(n,m,k);
% end
% end
% end
% c2old = sqrt(gamma0E2)/(h^2);
% PAI(iter)=Tmax;
% y(iter) = abs((Tmax-Tmaxold)/Tmaxold)+1e-100;
fprintf('Eu = %1.2fJ , SCA iteration: %1.0f: Tmax = %2.5f\n',Eu, iter, Tmax);
% x1old = x_nk1l + eps;
% x2old = x2 + eps;
if Tmaxold>Tmax
% break;
Tmaxold = Tmax;
end
% if iter>1
% Thandle.LineStyle=’:’;
% Thandle.Color=‘k’;
% Thandle.Marker=‘none’;
% end
g = reshape(q1(1,:),K,N);
w = reshape([g(:,:);q2(1,:)],(K+1)*N,1);
Trajx = [qu(1);w;qv(1)];
f = reshape(q1(2,:),K,N);
q = reshape([f(:,:);q2(2,:)],(K+1)*N,1);
Trajy = [qu(2);q;qv(2)];
Thandle = plot(Trajx, Trajy,'-r.');drawnow;
end
CVX output:
Successive approximation method to be employed.
SDPT3 will be called several times to refine the solution.
Original size: 347 variables, 168 equality constraints
15 exponentials add 120 variables, 75 equality constraints
Cones | Errors |
Mov/Act | Centering Exp cone Poly cone | Status
--------±--------------------------------±--------
0/ 0 | 0.000e+00 0.000e+00 0.000e+00 | Infeasible
Status: Infeasible
Optimal value (cvx_optval): +Inf