
简介这是一套面向通信工程、电子信息与数学等专业学生及研究者的Matlab复现资源聚焦无人机无线传感器网络中的节能数据采集问题。代码兼容Matlab2014、2019a与2024a采用参数化编程思路关键参数可方便调整并配有详细注释附赠案例数据可直接运行也能替换为自己的数据适合课程设计、期末大作业或毕业设计阶段快速上手。压缩包共26个文件包含9个m源码文件、11个bmp结果图、2个png示意图以及mat数据、说明文档和README等整体约420KB结构清晰便于按模块对照学习。目前已有97人学习下载。代码实现涵盖路径损耗、可达速率求解及多个优化子问题的求解可对无人机轨迹、节点唤醒与数据采集过程进行仿真帮助读者理解网络能耗机制并验证改进算法的效果。1. 无人机数据采集不省电问题不在飞在于等做无线传感器网络的人都知道地面节点电池容量有限每次采集任务最怕的不是数据多而是无人机绕着区域飞了一圈节点全在监听等待能量耗在空转上。这个Matlab复现项目解决的核心问题很直接给定地面传感器位置和能量预算无人机作为移动汇聚节点轨迹怎么规划、哪些节点在哪个时隙唤醒、数据速率怎么分配使得整个采集周期内采集到的数据量最大。代码基于Matlab 2014/2019a/2024a编写参数化程度高主程序和绘图脚本分离适合做通信方向课程设计、期末大作业和毕业设计的同学直接复现也适合刚开始接触无人机辅助WSN轨迹优化的从业者快速搭建仿真基线。文件包里solveP1.m、solveP2.m、solveP3.m三个脚本对应了三个子问题的求解这是目前该方向论文中最主流的块坐标下降分解框架值得动手拆一遍。2. 无人机辅助WSN的系统模型与能量-数据量权衡2.1 系统模型与变量定义先把场景固定下来一片区域内随机部署N个地面传感器节点每个节点有固定的坐标和有限的初始能量无人机从起始点出发以固定高度H飞行飞行周期为T将采集到的数据在飞行过程中通过无线链路传输到地面基站或直接接收节点上传的数据。这个项目里核心变量包括无人机的二维飞行轨迹、每个节点的唤醒调度变量某时刻节点是监听状态还是休眠状态、以及数据传输速率分配。这套模型里有几个先验假设需要搞清楚否则在Matlab里调参数时很容易对不上结果。路径损耗模型用的是自由空间损耗叠加视距/非视距概率模型高度固定时无人机与地面节点的距离只由水平投影距离决定所以轨迹优化本质上是在二维平面上做。唤醒调度是0-1变量睡眠状态的节点不接收也不发送数据能量消耗为0这直接影响无人机只有在节点唤醒时才能采集数据因此轨迹、调度、速率三者耦合在一起不能分开单独优化。2.2 为什么必须联合优化而不是分别贪心如果把问题拆成两段来做先飞一条经过所有节点的最短路径再让节点在无人机靠近时唤醒传输会发现采集数据量远低于理论最优。原因在于节点唤醒时间窗和无人机到达时间是由轨迹决定的而轨迹又要迁就唤醒窗口这本质上是鸡生蛋的问题。单独最优的路径可能让某些节点的唤醒窗口错开单独最优的调度又可能让节点长时间监听浪费能量。在能量约束下系统吞吐量最大化问题可以写成最大化所有节点在整个飞行周期内上传的数据量总和约束条件包括每个节点的能量上限、无人机的最大飞行速度、以及每时隙只能与一个节点通信的干扰避免约束。这个混合整数非凸问题直接求解是NP难的所以项目中把它分解为三个子问题用块坐标下降法BCD迭代求解。solveP1.m负责给定期望速率下的轨迹优化solveP2.m负责给定轨迹下的唤醒调度solveP3.m负责给定调度下的速率分配三者循环迭代直到目标值收敛。2.3 收敛判据与三个子问题的边界块坐标下降法收敛的前提是每个子问题都能找到最优解或至少单调不减的解。solveP1.m中轨迹子问题是非凸的通常用连续凸近似SCA做一阶泰勒展开每轮迭代内用CVX或内点法求一个凸近似问题solveP2.m中的唤醒调度在固定速率后实际上是一个带能量约束的0-1背包问题松弛为连续变量后结合投影即可solveP3.m在固定轨迹和调度后是凸问题因为对数速率函数是凹的可以直接求解。在MATLAB中判断迭代是否停止要看相邻两次迭代目标值的相对变化量一般阈值设为1e-3就够。文件里的Fig1.mat存的就是一组迭代过程的目标值序列可以画出来验证收敛性。这一层理解透了后面改参数、换数据、调到自己的场景时才不会盲目。3. Matlab复现的代码结构与核心脚本实现3.1 parameter_setting.m所有实验参数的入口打开压缩包后建议先看parameter_setting.m这是整个复现实验的配置中心。脚本用一个结构体把仿真参数统一管理起来结构大致的逻辑如下% parameter_setting.m 核心结构示意 param.N 8; % 传感器节点数量 param.T 100; % 无人机飞行周期秒 param.H 100; % 无人机飞行高度米 param.Vmax 30; % 无人机最大飞行速度m/s param.B 1e6; % 信道带宽Hz param.Pmax 0.1; % 无人机最大发射功率W param.noise 1e-11; % 噪声功率谱密度 param.E0 ones(1,param.N) * 100; % 每个节点的初始能量J param.dt 1; % 时隙长度秒T/dt 100个时隙所有仿真脚本通过调用这个参数文件来获取变量修改场景时只需要改动这个文件不需要在solveP1.m、solveP2.m里到处找硬编码值。param.dt决定了时隙数量时隙越细轨迹和时间调度的精度越高但每个子问题的矩阵规模也越大求解时间会明显上升。有一个值得注意的点param.dt从1改为0.5时时隙数从100变成200getPathLoss.m里生成的路径损耗矩阵维度也跟着变如果矩阵索引没对齐solveP1.m会报维度不匹配的错误。所以改时间步长时需要同时检查轨迹矩阵和调度矩阵的列数定义。3.2 getPathLoss.m与getAchievableRate.m信道和速率怎么算getPathLoss.m根据无人机位置矩阵与节点坐标计算每个时隙的路径损耗输入输出结构是标准的function [PL] getPathLoss(trajectory, node_pos, param) % trajectory: T_slots x 2 的无人机水平位置序列 % node_pos: N x 2 的地面节点坐标 % param.H: 无人机固定高度 % 输出 PL: T_slots x N 的路径损耗矩阵单位 dB num_slots size(trajectory, 1); num_nodes size(node_pos, 1); PL zeros(num_slots, num_nodes); for k 1:num_slots delta_x trajectory(k,1) - node_pos(:,1); delta_y trajectory(k,2) - node_pos(:,2); dist_3d sqrt(delta_x.^2 delta_y.^2 param.H^2); PL(k,:) 20*log10(4*pi*param.fc*dist_3d/3e8); % 自由空间路径损耗 end end这段代码里trajectory每一行是无人机在一个时隙的水平坐标node_pos是N行两列的节点坐标矩阵输出PL是时隙数乘以节点数的矩阵表示每个时隙无人机到每个节点的信道损耗。实际项目中这段代码还加了视距概率修正项但核心逻辑一致。getAchievableRate.m是在路径损耗基础上用香农公式把信噪比映射成速率% 给定路径损耗和发射功率计算可达速率bit/s/Hz function [R] getAchievableRate(PL, Ptx, param) snr Ptx ./ (10.^(PL/10) * param.noise * param.B); R param.B * log2(1 snr); end注意snr计算时PL要从dB转回线性值。Ptx可以是固定最大功率也可以在solveP3.m里优化分配这个函数的独立性设计方便了两种用法的切换。3.3 三个求解脚本的主循环装配核心的迭代逻辑是把三个子问题串起来主循环会直接在Figure1.m或Figure2a.m中运行大致结构如下% 初始化轨迹例如从圆心出发 traj repmat([500,500], T/dt, 1) ... schedule ones(num_nodes, num_slots); % 初始状态全唤醒 rates zeros(num_nodes, num_slots); for iter 1:30 % 给定调度和速率优化轨迹 traj_new solveP1(schedule, rates, param); % 给定新轨迹和速率优化唤醒调度 schedule_new solveP2(traj_new, rates, param); % 给定轨迹和调度优化速率分配 rates_new solveP3(traj_new, schedule_new, param); % 检查收敛计算目标函数值并比较 obj_old sum(rates(:)); obj_new sum(rates_new(:)); if abs(obj_new - obj_old) / obj_old 1e-3 break; end endsolveP1.m内部是连续凸近似加CVX求解的过程没有CVX工具箱的环境下可以改成用fmincon做序列二次规划近似但需要手动把凸近似矩阵写出来不建议新手在第一次复现时替换。solveP2.m里注意唤醒调度是一个二进制矩阵直接松弛后求解最后用阈值0.5投影回0/1值这个投影操作可能会带来目标函数值的轻微波动属于正常现象。solveP3.m是三个子问题中最容易求解的因为速率分配在固定信道增益后是一个标准的几何规划问题内点法迭代几次就能收敛。4. 复现实验跑通基线、改参数、替换数据集4.1 从Figure1.m到Figure2b.m的执行顺序项目文件里的脚本可以从parameter_setting.m开始依次执行。首次跑通建议按照Figure1.m、Figure2a.m、Figure2b.m的顺序来因为这三个脚本分别绘制了不同飞行周期下的结果对比。Figure1.m的完整执行步骤是% 步骤1加载参数设置 run(parameter_setting.m); % 步骤2生成节点位置并初始化轨迹 rng(0); % 固定随机种子保证可复现 node_pos 1000 * rand(param.N, 2); traj_init linspace(0, 2*pi, param.T/param.dt); traj_init [500 400*cos(traj_init), 500 400*sin(traj_init)]; % 步骤3调用主迭代循环 [traj_opt, schedule_opt, rates_opt] main_loop(node_pos, param); % 步骤4绘制与保存 figure; plot(traj_opt(:,1), traj_opt(:,2), r-, LineWidth, 1.5); hold on; scatter(node_pos(:,1), node_pos(:,2), 60, filled); saveas(gcf, fig1_result.png);注意rng(0)这一步很关键如果不固定随机种子每次运行生成不同的节点坐标轨迹图和调度结果都会变不利于对比实验。main_loop函数在压缩包中没有单独出现它实际被整合在绘图脚本里启动每个Figure脚本之前它会自动执行迭代求解所以直接点击运行Figure1.m即可看到收敛后的轨迹图。4.2 T40、T50、T100轨迹图差异的解读文件包里有T40_trajectory.bmp、T50_trajectory.bmp和T100_trajectory.bmp三张图它们对应param.T 40/50/100三个飞行周期下的最优轨迹。把param.T改到不同值然后分别运行Figure2a.m会看到两个明显的趋势飞行周期T越小无人机速度上限约束相对越紧因为飞行路程固定而时间变短轨迹会偏向直线冲过节点的最短路径牺牲部分节点的传输速率来保证全部节点至少访问一次飞行周期T越大能量预算不变的情况下可飞行路程变长无人机会绕到每个节点正上方悬停通信轨迹走向与TSP路径接近但唤醒调度优化会让它在低能量节点停留更长时间。这个对比在论文里常用来展示飞行周期对无人机辅助采集系统能效的影响分析。做课程设计写报告时这三张图可以直接作为仿真结果对比图加入章节配上每一张图对应的参数配置表和解释即可。4.3 数据替换成自己的传感器能量值如果要做自己的场景替换数据的方式很直接。假设你有本地传感器坐标与能量值存在一个CSV文件% 读取自己的节点坐标和初始能量 node_table readtable(node_data.csv); node_pos [node_table.x, node_table.y]; param.E0 node_table.initial_energy; param.N size(node_pos, 1);替换之后要检查两个地方。第一param.N是否仍然为8如果不是需要把parameter_setting.m中的param.N同步修改。第二initial_energy列的数值分布不能差异过大如果某些节点能量只有其他节点的十分之一solveP2.m中能量约束会让这些节点几乎睡满全程导致这些节点采集速率极低这不是bug而是优化结果需要在报告中解释这一现象。替换后重新运行Figure1.m如果出现Variable param has been deleted的报错通常是脚本执行过程中不小心用clear all清掉了param结构体注意在循环内只清理临时变量不要清除param。4.4 修改节点数量N时易出现的错误把节点数从8改成16是课程设计里最常见的需求。直接改param.N 16还不够需要同步处理初始轨迹和随机种子。Figure1.m中生成节点位置的代码会根据param.N生成随机坐标但如果主循环和绘图块引用了node_pos而node_pos在上一次运行时被保存为8行两列的矩阵再次运行时会报维度不匹配的错误建议在修改param.N后对整个工作区执行clear并重新运行所有脚本。另外getAchievableRate.m中的信噪比计算对节点数并不敏感复杂性在于getPathLoss.m的矩阵维度随节点数线性增长16个节点时三重循环的耗时大约是8个节点的6到8倍这是正常的不需要优化代码耐心等待即可。5. 进阶验证收敛曲线检查与常见坑的定位5.1 用wake-up.bmp反向验证调度合理性wake-up.bmp和T100_wake-up.bmp展示了不同周期下的唤醒时刻分布这是一张热力图横轴是时隙序号纵轴是节点编号颜色深浅表示唤醒深度或累计通信时长。验证结果是否合理的方法是选中一个能量非常低的节点例如第4号节点初始能量只有其他节点的20%观察它的唤醒行中是否有大段的浅色块。如果完全没有说明该节点能量预算不足以支持传输这是预期结果如果出现了深色块说明该节点在某些时隙有高功率通信此时需要检查它的剩余能量是否在约束范围内。5.2 三个子问题的收敛性观测方法在运行主循环时把每轮迭代的目标函数值存下来并绘制成曲线是比较推荐的验证手段obj_history zeros(max_iter, 1); for iter 1:max_iter % ... 三个solveP的调用 ... obj_history(iter) sum(rates_new(:)); if iter 1 abs(obj_history(iter) - obj_history(iter-1)) 1e-4 break; end end figure; plot(obj_history(1:iter), b-o); xlabel(迭代次数); ylabel(目标函数值);正常情况下曲线呈单调上升后趋平。如果曲线上下波动优先检查solveP2.m中阈值投影是否太粗糙例如0.5阈值附近的值对目标影响敏感可以将阈值调整到0.6或0.7或改用排序保留固定数量的唤醒时隙。5.3 运行报错的三类常见原因第一类是矩阵维度对不齐报错信息通常是Matrix dimensions must agree。位置在getPathLoss.m函数调用处时将size(trajectory,1)打印出来与param.T/param.dt比较确认时隙数是否一致。第二类是CVX提示Disciplined convex programming error这说明solveP3.m或者solveP1.m中的约束写法破坏了凸性检查是否有slog(1exp())或除法出现在约束不等式左侧。第三类是Assignment has more non-singleton rhs dimensions than non-singleton subscripts看到这种报错就检查solveP2.m中的调度矩阵schedule_new在投影后又按列赋值的索引顺序确认列数与时隙数相同。5.4 修改Matlab版本兼容性的处理建议代码兼容2014、2019a和2024a三个版本但不同版本间有一些小的函数差异。2014版本不支持string类型数组和contains函数代码中如果使用了contains判断文件名或参数字符串需要改成strfind或ismember的匹配方式。2024a版本对legend和colororder的默认颜色序列做了更新旧代码生成的Figure2a.m曲线图如果不额外指定颜色会变得不好区分可以在绘图前明文指定线条颜色。如果遇到graphics object array赋值失败通常是2014版本返回的是句柄而不是对象将句柄转为double后就能兼容。综合来看优先推荐在2019a或2024a环境运行2014环境适合验证关键脚本的语法兼容性。本文还有配套的精品资源点击获取