2020-11-16

function [Xr, Theta] = solveP2(Qr) %Qr 为初始的轨迹
global M K Ek rk

% Equation (1)
pathloss = eq1(Qr);

% Equation (4)
R = eq4(pathloss);

cvx_begin quiet %(P2)
%variable 是待优化的变量
variables x(M,K);v %M X K
variables theta(1,1); %松弛变量

minimize(sum(theta))
subject to:
0 <= x <= 1; % Constraint(11)
sum(x.*Ek, 1) <= theta % Constraint(5)
sum(x.*R, 1) >= rk % Constraint(6)
sum(x , 2) <= 1 % Constraint(7)
cvx_end

Xr = x;

Theta = sum(theta);

end

function [pathloss] = eq1(Qr)
global K M H w %传感器 时隙数量 飞行高度 %用于记录每个时隙下每个SN的坐标状态

Qr = repelem(Qr,1,1,K); %matrix
d = (H^2 + sum((Qr - w).^2 ,1)).^(0.5) %1 X100 X 4 100个时间间隔 4个传感器
pathloss = getPathLoss(d); % eq1
pathloss = reshape(pathloss, [M,K] ); % matrix
end

function [R] = eq4(pathloss)
R = getAchievableRate(pathloss);
end

function [dB] = dec2dB(dec)
dB = 10*log10(dec);
end

function [dec] = dB2dec(dB)
dec = 10.^(dB/10);
end

function [R] = getAchievableRate(pathloss)
% input, pathloss is dB scale
% output, R’s unit is bps/hz
global F_1 Pk sigma_2 Lamda
R = log2(1+(F_1PkdB2dec(pathloss)) / (dB2dec(sigma_2)*dB2dec(Lamda)));
end

function [pathloss] = getPathLoss(d)
% input, d is decimal scale 十进制
% output pathloss is dB scale
global beta0 alpha

da = -alpha*dec2dB(d);
pathloss = beta0 + da;

end