Hello everyone. When I use cvx tool to solve the following convex problem, I found the values of expressions T_f_auxi, T_w, T_tr and Y_time are not equal to the corresponding assignments, which made the optimal value of objective function seem to be wrong. I want to know if the following code is written incorrectly, which leads to this result. I would appreciate it if you could help me.
function [E, Q, d_tol, z, actT_h, T_h, Y_time, T_f, T_f_auxi, T_w, T_tr, T_com] = optHori1(M,q,alpha,Data,Freq,X,Y,B,p,h_0,sigma,gamma_th,mu,nu,Phi,varphi,iniQ,H,d_th,v,E_max,N,UAV_USVlist,x0,y0)
q0 = [x0 y0 0];
gamma_0 = p*h_0/sigma;
v_in = 4.03; % mean rotor induced velocity in hovering
v_tip = 120; % tip speed of rotor blade
d0 = 0.6; % fuselage drag ratio
varrho = 1.225; % air density
beta = 0.05; % rotor solidity
A = 0.503; % rotor disc area
p0 = 2592varrhobetaA;
p1 = 1.120^(3/2)/sqrt(2varrhoA);
p_h = p0+p1;
p_f = zeros(1,N);
for k=1:N
p_f(k) = p0*(sqrt(1+v(k)^4/(4v_in^4))-v(k)^2/(2v_in^2))^(1/2)+p1*(1+3v(k)^2/(v_tip^2))+1/2d0varrhobetaAv(k)^3;
end
% Square of distance between USV and UAV
d_hat = zeros(M,N);
for i=1:M
for k=1:N
d_hat(i,k) = norm(iniQ(k,:)-q(i,:))^2 + H^2;
end
end
d_th_bar = repmat(d_th^2, N, N); %Minimum safe distance between drones(UAVs)
% computation time cost of UAVs
UAV_f = zeros(1,N); % computation capacity of UAVs
for k=1:N
runTaskNum = length(find(alpha(:,k)));
UAV_f(k) = mu^(Phi - runTaskNum)/nu;
end
T_com = zeros(M,N); % the computing time of tasks
for i=1:M
for k=1:N
T_com(i,k) = Freq(i)/UAV_f(k);
end
end
Y_bar = repmat(Y,1,N);
% computation energy consumption of UAvs
E_com = zeros(1,N);
for k=1:N
[~,~,curUAV_USVlist] = find(UAV_USVlist(k,:));
% USV_num = length(curUAV_USVlist);
% curUAV_f = mu^(Phi-USV_num)/nu*1e9;
USV_freq = 0;
for j = 1:length(curUAV_USVlist)
USV_freq = USV_freq + Freq(curUAV_USVlist(j));
end
E_com(k) = varphi * UAV_f(k)^2 * USV_freq;
end
cvx_begin quiet
variable Q(N,2) nonnegative % deployment of UAVS
variable d_tol(1,N) nonnegative % flight distance of UAV from original point to deployment point
variable z(M,N) nonnegative
variable T_h(1,N) nonnegative; %hovering time of UAVs
expression d(M,N);
for i=1:M
for k=1:N
d(i,k) = square_pos(norm(Q(k,:) - q(i,:))) + H^2;
end
end
expression r_lb(M,N);
for i=1:M
for k=1:N
r_lb(i,k) = B*log2(1+gamma_0/d_hat(i,k))-...
(B * gamma_0 / (log(2)*(d_hat(i,k))*(d_hat(i,k)+gamma_0)) * (d(i,k) - d_hat(i,k)));
end
end
% lower bound of distance between UAVs
expression d_UU_lb(N,N);
for k=1:N
for j=1:N
if j == k
d_UU_lb(k,j) = d_th ^ 3;
else
d_UU_lb(k,j) = norm(iniQ(k,:)-iniQ(j,:))^2 + 2* (iniQ(k,:) - iniQ(j,:)) * ((Q(k,:) - Q(j,:)) - (iniQ(k,:)-iniQ(j,:))).';
end
end
end
expression Y_time(M,N); % time windows
expression T_f(1,N);
expression T_f_auxi(1,N);
expression T_w(M,N);
expression T_tr(M,N);
expression d_tol_auxi(1,N);
for k = 1:N
T_f(k) = norm([Q(k,:),H] - q0)/v(k);
d_tol_auxi(k) = power(d_tol(k),1/2);
T_f_auxi(k) = d_tol_auxi(k)/v(k);
end
for i = 1:M
for k = 1:N
T_w(i,k) = pos(X(i) - T_f_auxi(k));
T_tr(i,k) = Data(i)* inv_pos(z(i,k));
Y_time(i,k) = T_f(k) + T_w(i,k) + T_tr(i,k) + T_com(i,k);
% Y_time(i,k) = norm([Q(k,:),H] - q0)/v(k) + pos(X(i) - sqrt(d_tol(k))/v(k)) + Data(i)* inv_pos(z(i,k)) + T_com(i,k);
% Y_time(i,k) = norm([Q(k,:),H] - q0)/v(k) + pos(X(i) - T_a(k)) + Data(i)* inv_pos(z(i,k)) + T_com(i,k);
end
end
expression actT_h(M,N);
T_h_bar = repmat(T_h,M,1); % 将盘旋时间扩充成 M*N,便于约束条件书写
for i=1:M
for k=1:N
actT_h(i,k) = pos(X(i) - sqrt(d_tol(k))/v(k)) + Data(i)*inv_pos(z(i,k)) + T_com(i,k);
end
end
% energy constraint
expression E(1,N);
expression E_h(1,N); % hovering energy consumption of UAVs
expression E_f(1,N); % flying energy consumption of UAVs
expression d_tol_lb(1,N);
for k=1:N
E_h(k) = p_h * T_h(k);
E_f(k) = p_f(k) * T_f(k);
E(k) = E_h(k) + E_f(k) + E_com(k);
d_tol_lb(k) = norm(iniQ(k,:) - [x0, y0])^2 + 2* (iniQ(k,:)-[x0 y0])* (Q(k,:)-iniQ(k,:)).' + H^2;
end
minimize(sum(E)); % 目标函数
subject to
alpha.*d <= gamma_0/gamma_th;
d_UU_lb >= d_th_bar; % safe distance constraint of UAVs
alpha.*Y_time <= Y_bar; % time windows constraint
d_tol <= d_tol_lb;
alpha.*actT_h <= T_h_bar; % hovering time constraint
E <= E_max; %maximum energy constraint
z <= r_lb;
cvx_end
end