当前位置：首页 > news >正文

多用户情况下的无人机通信轨迹和调度联合优化开源代码

news 2026/3/26 17:28:59

前言

知乎专栏已更新15篇博文及开源代码，欢迎关注

通信、雷达算法仿真 - 知乎https://www.zhihu.com/column/c_1990131735602697785

代码介绍

两用户情况下的无人机通信轨迹和调度联合优化开源代码

公式推导

开源代码

%匀加速模型，两个用户 %匀加速模型，两个用户，联合优化轨迹和功率分配 clc,clear; close all; s=rng(1);%random seed %% paramters %飞行场（200m*50m）的经纬度,格式 %左上经纬度 %右上经纬度 %左下经纬度 %右下经纬度 pos = [109.14190184, 34.04978701; 109.14383636, 34.05046442; 109.14222334,34.04948425; 109.14420065, 34.04949029]; qs = [0,0]'; % start position,还得调 qe = [160,50]';% end position % s1 = [10,30]'; % BS station % s2 = [130,20]'; % BS station s1 = [0,20]'; % BS station s2 = [140,30]'; % BS station dis = norm(qe-qs); % distance between qs and qe T = 50; % times slots delta_T = 0.5;% delta T,s vmax =dis/(T*delta_T)*8; % maximial velocity vmin = 3;% minimial velocity amax = 5;% max accelerated velocity noise_power = 1e-7;%50dbm beta0 = 1e-3;% L0 pathloss iter_num = 20; % number of iterations P1 = zeros(T,iter_num); P2 = zeros(T,iter_num); P_avg = 2;%average power P_max = 4;%maximial power P1(:,1) = P_avg/2;%power allocation P2(:,1) = P_avg/2;%power allocation q = zeros(T,2,iter_num+1);% positions variables v = zeros(T,2,iter_num+1);% velocity variables v(:,:,1) = sqrt(1/2*dis/(T*delta_T));% a = zeros(T,2,iter_num+1);% accelerated velocity epsilon = 1e-3; % stop criterion H =20;% fixed hegith u = zeros(T,iter_num+1);% slack variables w = zeros(T,iter_num+1);% slack variables % straight line flying q(:,1,1)=linspace(qs(1),qe(1),T);%initial point q(:,2,1)=linspace(qs(2),qe(2),T); val1 = zeros(1,iter_num+1); val2 = zeros(1,iter_num+1); for ii = 1:T u(ii,1)= norm(q(ii,:,1)-s1)^2; end for ii = 1:T w(ii,1)= norm(q(ii,:,1)-s2)^2; end %% alternaating optimization for iter=1:iter_num % trajectory optimization subproblem cvx_begin quiet variable q1(T,2) variable u1(T,1) nonnegative variable w1(T,1) nonnegative variable v1(T,2,1) variable a1(T,2,1) rate = 0; % calculate objective function for ii = 1:T beta1_hat=beta0*P1(ii,iter)/noise_power; beta2_hat=beta0*P2(ii,iter)/noise_power; rate = rate + log2(1+beta1_hat/(u(ii,iter)+H^2))-1/log(2)*beta1_hat/((u(ii,iter)+H^2)^2+beta1_hat*(u(ii,iter)+H^2))*(u1(ii)-u(ii,iter))... +log2(1+beta2_hat/(w(ii,iter)+H^2))-1/log(2)*beta2_hat/((w(ii,iter)+H^2)^2+beta2_hat*(w(ii,iter)+H^2))*(w1(ii)-w(ii,iter)); end maximize (rate) subject to % start and end potisiton constraints q1(1,1)==qs(1,1); q1(1,2)==qs(2,1); q1(T,1)==qe(1,1); q1(T,2)==qe(2,1); % start and end velocity constraints v1(1,:) == [vmin*2,vmin*2] v1(T,:) == [vmin*2,vmin*2] % accelerated velocity constraint norm(a1(1,:))<=amax % slack variables constraints norm([1-u1(1),2*(q1(1,:)-s1')])<=1+u1(1) norm([1-w1(1),2*(q1(1,:)-s2')])<=1+w1(1) for kk=2:T % maximal and minimal velocity constraints norm(v1(kk,:))<=vmax % v(kk,:,iter)*v(kk,:,iter)'+2*v(kk,:,iter)*(v1(kk,:)'-v(kk,:,iter)')>=vmin^2 % % maximial accelerated constrains norm(a1(kk,:))<=amax % position constrains q1(kk,:)==q1(kk-1,:)+v1(kk,:)*delta_T+1/2*a1(kk,:)*delta_T^2 v1(kk,:)==v1(kk-1,:)+a1(kk,:)*delta_T % slack constrains norm([1-u1(kk),2*(q1(kk,:)-s1')])<=1+u1(kk) norm([1-w1(kk),2*(q1(kk,:)-s2')])<=1+w1(kk) end cvx_end % save optimized results q(:,:,iter+1)=q1; v(:,:,iter+1)=v1; a(:,:,iter+1)=a1; u(:,iter+1)=u1; w(:,iter+1)=w1; % val(iter+1)=cvx_optval; val1(iter+1)=0; for ii = 1:T beta1_hat=beta0*P1(ii,iter)/noise_power; beta2_hat=beta0*P2(ii,iter)/noise_power; val1(iter+1) = val1(iter+1) + log(1+beta1_hat/(norm(q(ii,:,iter+1)'-s1)^2+H^2))/log(2)+ log(1+beta2_hat/(norm(q(ii,:,iter+1)'-s2)^2+H^2))/log(2); end % power optimization subproblem cvx_begin quiet variable p1(T,1) nonnegative variable p2(T,1) nonnegative rate1=0; for ii = 1:T beta_hat1=beta0*p1(ii)/noise_power; beta_hat2=beta0*p2(ii)/noise_power; rate1 = rate1 + log(1+beta_hat1/(norm(q(ii,:,iter+1)'-s1)^2+H^2))/log(2)+ log(1+beta_hat2/(norm(q(ii,:,iter+1)'-s2)^2+H^2))/log(2); end maximize (rate1) subject to %maimium power constraint for kk=1:T p1(kk)<=P_max; p2(kk)<=P_max; end %average power constraint 1/T*(sum(p1+p2))<=P_avg;%可以注释掉 cvx_end val2(iter+1)=0; P1(:,iter+1)=p1; P2(:,iter+1)=p2; for ii = 1:T beta_hat1=beta0*p1(ii)/noise_power; beta_hat2=beta0*p2(ii)/noise_power; val2(iter+1) = val2(iter+1) + log(1+beta_hat1/(norm(q(ii,:,iter+1)'-s1)^2+H^2))/log(2)+ log(1+beta_hat2/(norm(q(ii,:,iter+1)'-s2)^2+H^2))/log(2); end disp([iter,val1(iter+1),val2(iter+1)]); if abs(val2(iter+1)-val2(iter))<=epsilon break end end %% plot trajectory figure hold on %plot trajectory plot(q(:,1,1),q(:,2,1),'r-.o','linewidth',1.5); plot(q(:,1,iter+1),q(:,2,iter+1),'b-d','linewidth',1.5); plot(s1(1),s1(2),'s','MarkerFaceColor', 'g','markersize',10) plot(s2(1),s2(2),'s','MarkerFaceColor', 'g','markersize',10) % plot positions plot(qs(1),qs(2),'s','MarkerFaceColor', 'g','markersize',10) plot(qe(1),qe(2),'s','MarkerFaceColor', 'g','markersize',10) % add text text(s1(1)+5,s1(2)+5,'User1') text(s2(1)-17,s2(2)-5,'User2') text(qs(1)-5,qs(2)-5,'Inital') text(qe(1)-5,qe(2)+5,'Final') grid on % lim & label & box xlim([qs(1)-15,qe(1)+25]) ylim([qs(2)-10,qe(2)+10]) legend('straight flight','Optimized trajectory','Location','southeast') box on xlabel('x(m)') ylabel('y(m)') title(['T=',num2str(T*delta_T),'s']) %% plot power figure hold on plot((1:T)*delta_T,P1(:,iter+1),'b-o','linewidth',1.5); plot((1:T)*delta_T,P2(:,iter+1),'r-v','linewidth',1.5); plot((1:T)*delta_T,P1(:,iter+1)+P2(:,iter+1),'k-d','linewidth',1.5); xlabel('Times (s)'); ylabel('Power (W)'); grid on box on legend('user 1','user 2','sum power'); %% plot convergence figure plot(val2(1:iter+1),'b-o','linewidth',1.5); grid on xlabel('iteration number') ylabel('objective value') box on set(gca,'fontsize',10) title('Convergence curve') % %% % figure % hold on % for ii = 1:iter % plot(q(:,1,ii),q(:,2,ii),'linewidth',1.5); % % pause(1) % end %% XY->latitude and longitude q_opt = q(:,1,iter+1); q2 = zeros(size(q_opt)); q2(:,1) = 109.14383636 + q_opt * 0.00001141; q2(:,2) = 34.04948425 + q_opt * 0.00000899; q3=roundn(q2,-5);%限定位数。小数点后5位 writematrix(q3,'Trajectory.csv') res = [q3,roundn(P1(:,iter+1),-5),roundn(P2(:,iter+1),-5)]; writematrix(res,'Trajectory_Power.csv'); %用户位置转经纬度 s1_new = zeros(size(s1)); s2_new = zeros(size(s2)); s1_new(1) = 109.14383636 + s1(1)* 0.00001141; s1_new(2) = 34.04948425 + s1(2)* 0.00000899; s2_new(1) = 109.14383636 + s2(1)* 0.00001141; s2_new(2) = 34.04948425 + s2(2)* 0.00000899;