ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

块坐标下降法分解无人机通信网络联合优化问题及Python实现

块坐标下降法分解无人机通信网络联合优化问题及Python实现 前段时间我在做一套低空物流场景的通信保障系统任务说起来很简单一架多旋翼无人机挂载轻量化基站在指定区域内给地面若干个移动用户提供上行回传链路。用户位置会变无人机飞哪里、发射功率给多少、哪些用户占用哪些信道——三个问题放在一起看就是一团乱麻。把这三个变量同时丢给优化器去解非线性、非凸、还带整数约束一般的求解器根本撑不住。后来我用块坐标下降法BCD把问题拆成三个子问题轮流求解用Python写了第一版原型效果比我预想的好不少。这篇文章就围绕这套方案把建模、代码和调试过程完整记录下来希望对做类似项目的朋友有参考价值。1. 从问题到方法为什么无人机通信网络偏要用块坐标下降法1.1 无人机通信优化到底在优化什么无人机通信网络优化不是单一维度的调参它通常同时包含三组互相耦合的变量无人机的位置坐标或者连续航迹、各用户的发射功率或无人机基站的发射功率分配、信道资源的分配方式。以最常见的场景为例一架无人机在空中做基站地面散布着多个用户用户请求上传数据。你要决定无人机悬停在哪个位置使整体信道质量最好你要决定功率怎么分配让距离远、信道差的用户不至于完全没速率你还要决定信道怎么分让同时通信的用户不互相干扰。这三件事单独拿出来都不难但合在一起就是典型的联合优化问题。目标函数通常写成系统总吞吐量最大化约束条件包含功率上限、信道互斥约束、以及无人机活动区域边界。这种问题用暴力搜索不现实用全局优化工具又太慢必须找一种能工程化落地的迭代求解思路。1.2 为什么“拆块”能行得通块坐标下降法的核心思想一句话固定其他变量只优化其中一块然后轮流往下走。你可以把它理解成整理一间乱糟糟的房间如果同时把衣服、书本、杂物全收拾好大脑会直接过载但如果你先只叠衣服其他东西都不动叠完再理书本最后处理杂物整个过程会顺畅得多。每轮只做一件事每一件事的目标都很明确循环几轮之后房间就大致整洁了。把这个思路套用到无人机网络里就是固定位置和信道分配优化功率固定位置和功率优化信道分配固定功率和信道分配再优化位置。这样反复迭代每一轮都让总吞吐量至少不下降最终稳定在一个比较好的局部最优解。这里的“局部最优”要客观看待原始问题非凸我们不指望BCD能找到全局最优但只要初始化不是太极端BCD在工程上的收敛速度和解的质量都非常够用。相比遗传算法、粒子群这类启发式方法BCD的每一步都有清晰的数学含义也更方便定位问题。1.3 何时该用BCD、何时别硬上BCD不是万能钥匙。如果你的变量之间耦合非常深比如每个变量的变化都会剧烈改变其他变量的可行域那拆块之后每轮解出来的结果可能在下一轮马上失效迭代会来回震荡甚至发散。另外如果某个子问题本身就很复杂比如包含大量整数变量和强非线性那拆块之后单块依然难解意义也不大。无人机通信这类问题之所以适合BCD是因为它天然分成“连续的位置”“连续的功率”“离散的信道”三类变量每一类变量在固定其他两类之后都有相对成熟的解法。比如固定信道和位置后功率分配是一个凸问题可以用注水法直接解固定位置和功率后信道分配是一个匹配问题可以用贪心或匈牙利算法。这种“各有各的招”的结构就是BCD最能发挥价值的场景。所以在动笔写代码之前先想清楚你的问题能不能像洋葱一样一层层剥开这一步比调参重要得多。2. 建模与目标函数先给优化问题一个数学上的“靶子”2.1 场景设定与信道模型我先定一个具体的仿真场景方便后面所有代码复用区域大小为 500m × 500m无人机飞行高度固定为 100m地面有 8 个随机分布的用户无人机上有一个 8 信道的正交频分多址OFDMA系统每个信道在同一时刻最多分配给一个用户保证信道之间互不干扰无人机总发射功率上限 1W噪声功率按 -110dBm 折算为线性值信道模型采用“大尺度路径损耗 小尺度衰落”的简化形式。信道增益的计算方式是无人机与用户之间的三维距离 d路径损耗指数 α 取 2.2参考距离 1m 处增益为 β0。这样每个用户在每个信道上的增益就不完全一样因为小尺度衰落矩阵会为每个用户-信道对引入独立的随机系数。说到底建模不是越复杂越好而是要为算法服务。我第一版只保留了影响优化结果最关键的路径损耗和衰落把阴影衰落、天线方向图、多径时延这些先放到一边。要是第一版就加上一堆实际因素问题复杂度会直接失控反而无法验证BCD模块本身是否正确。2.2 目标函数与约束条件我们的优化目标写出来是这个形式maximize R Σₖ Σₙ aₖₙ B log₂(1 pₙ γₖₙ)约束条件Σₙ pₙ ≤ Pmaxpₙ ≥ 0Σₖ aₖₙ ≤ 1aₖₙ ∈ {0,1}x_q, y_q 限制在 [0, L] × [0, L] 范围内。其中 γₖₙ (β0 · |gₖₙ|²) / ((dₖ² H²)^(α/2) · σ²)dₖ是用户k到无人机的水平距离。这个形式包含三组变量实数 x_q、y_q实数 pₙ以及二进制变量 aₖₙ。为什么要同时优化这些因为它们是强耦合的。无人机位置决定了所有信道的 γₖₙγₖₙ决定了信道分配的好坏功率分配又必须在给定信道分配后才能做到最优。三个变量互相依赖但依赖关系又是可以逐层解开的这就是第1节说的“洋葱结构”。参数取值说明区域边长 L500m无人机活动范围高度 H100m固定不优化用户数 K8随机分布信道数 N8每信道单用户Pmax1W总发射功率噪声功率 σ²1e-11 W约 -110dBm带宽 B1MHz每信道带宽路径损耗指数 α2.2城市低空环境近似第一版我没加入最小速率约束比如“每个用户至少分到多少速率”。原因很简单最小速率约束会让功率分配子问题失去标准的闭式解结构从一个干净的注水问题变成需要额外拉格朗日乘子迭代的问题。先把主干算法跑通再加这类实际约束是更稳妥的开发顺序。2.3 三个子问题的划分与求解思路按变量类型拆开就得到三个子问题固定位置和信道分配优化功率。此时目标函数对 pₙ 都是独立的凸项而且总功率约束是线性的可以直接用注水法求闭式解。固定位置和功率优化信道分配。这是一个离散匹配问题每个信道最多分配给一个用户贪心匹配可以快速得到一个可行解计算成本低。固定功率和信道分配优化位置。这个子问题依然非凸但目标函数对位置相对平滑用梯度上升法配合投影就能逐步改善。三个子问题的求解器各写各的互不干扰最后套一个主循环。这也是我特别喜欢BCD的地方你不需要一个万能优化器而是把问题拆成几个你熟悉的小工具分别解决后再组装。下一节就逐一上代码。3. Python实现块坐标下降法的主循环怎么落代码3.1 准备环境与生成基础数据实现只需要Python 3.9以上版本和numpy不需要深度学习框架不需要额外优化器。先把仿真参数和用户坐标准备出来import numpy as np L 500.0 H 100.0 K 8 N 8 Pmax 1.0 sigma2 1e-11 B 1e6 alpha 2.2 beta0 1e-3 rng np.random.default_rng(42) user_pos rng.uniform(0, L, size(K, 2)) # 每个用户在每个信道上独立的衰落系数固定下来保证可复现 small_scale rng.exponential(size(K, N)) def compute_gamma(xq, yq): gamma np.zeros((K, N)) for k in range(K): dx user_pos[k, 0] - xq dy user_pos[k, 1] - yq d2 dx * dx dy * dy H * H path_loss beta0 / (d2 ** (alpha / 2)) / sigma2 gamma[k, :] path_loss * small_scale[k, :] return gamma这里把随机种子固定为42每次跑出来的结果完全一致调试的时候很重要。小尺度衰落系数只要生成一次就好后续迭代中保持不变否则每次调用compute_gamma都重新随机目标函数就会出现莫名其妙地跳变。3.2 功率分配子问题注水法的工程化写法固定位置和信道分配时每个信道上的功率分配互不影响只有一个总功率上限约束。对每个被使用的信道n最优功率形式是pₙ max(0, 1/λ − 1/γₐₖₜᵢᵥₑ)其中 λ 是拉格朗日乘子γₐₖₜᵢᵥₑ 是该信道当前分配到的用户的信道增益。λ 要选到让总功率等于Pmax由于左侧随λ递减用二分搜索最稳定。def update_power(gamma, A): p np.zeros(N) channel_user np.argmax(A, axis0) active_ch np.where(A.sum(axis0) 0)[0] if len(active_ch) 0: return p gamma_active np.array([gamma[channel_user[n], n] for n in active_ch]) gamma_active np.maximum(gamma_active, 1e-12) lo, hi 1e-10, 100.0 for _ in range(100): lam (lo hi) / 2 pv np.maximum(0, 1.0 / lam - 1.0 / gamma_active) if pv.sum() Pmax: lo lam else: hi lam lam (lo hi) / 2 pv np.maximum(0, 1.0 / lam - 1.0 / gamma_active) for idx, n in enumerate(active_ch): p[n] pv[idx] return p需要注意两点。第一gamma_active要做下限保护否则信道增益为0时会出现除零。第二这个子问题之所以有闭式结构是因为每个信道最多服务一个用户不存在同信道干扰。如果改成多用户复用同一个信道SINR公式就会引入交叉项注水法随之失效这一点在扩展部分再聊。3.3 信道分配子问题贪心匹配信道分配是离散问题最优解可以通过匈牙利算法得到但第一版原型用贪心足够。思路很简单每轮找一个“用户-信道”组合使得在当前功率分配下新增速率最大然后锁定这个组合直到所有用户都有信道或者没有可用的空信道。def update_assignment(gamma, p): A np.zeros((K, N)) assigned_user set() assigned_chan set() while len(assigned_user) K and len(assigned_chan) N: best_delta -1.0 best_k -1 best_n -1 for k in range(K): if k in assigned_user: continue for n in range(N): if n in assigned_chan: continue if p[n] 0 and gamma[k, n] 0: delta B * np.log2(1 p[n] * gamma[k, n]) if delta best_delta: best_delta delta best_k k best_n n if best_k -1: break A[best_k, best_n] 1 assigned_user.add(best_k) assigned_chan.add(best_n) return A这个贪心的计算量是O(K²N)在K和N都是个位数时完全可接受。贪心不保证全局最优但它产生的解每个信道的速率增量都是正向的配合BCD主循环整体目标函数依然会保持单调上升。实际项目中如果你想追求更好的信道匹配可以把这个函数替换成匈牙利算法版本主循环其他部分不用改。3.4 无人机位置更新梯度投影法的实现细节位置子问题是三个子问题里最复杂的因为目标函数对位置没有简单的闭式最优解。我用的是数值梯度加回溯线搜索实现简单也足够稳。def objective(gamma, A, p): rate 0.0 for n in range(N): ks np.where(A[:, n] 0)[0] for k in ks: rate B * np.log2(1 p[n] * gamma[k, n]) return rate def objective_with_pos(xq, yq, A, p): gamma compute_gamma(xq, yq) return objective(gamma, A, p) def update_position(xq, yq, A, p, f_current): delta 0.5 lr 5.0 base_f f_current for _ in range(30): fx1 objective_with_pos(xq delta, yq, A, p) fx2 objective_with_pos(xq - delta, yq, A, p) fy1 objective_with_pos(xq, yq delta, A, p) fy2 objective_with_pos(xq, yq - delta, A, p) grad_x (fx1 - fx2) / (2 * delta) grad_y (fy1 - fy2) / (2 * delta) norm np.hypot(grad_x, grad_y) if norm 1e-12: break scale min(lr, 10.0 / norm) new_x np.clip(xq scale * grad_x / norm, 0, L) new_y np.clip(yq scale * grad_y / norm, 0, L) new_f objective_with_pos(new_x, new_y, A, p) if new_f base_f: return new_x, new_y, new_f lr * 0.7 if lr 0.01: break return xq, yq, base_f这里有一个细节容易被忽略数值差分步长delta取0.5m如果太小人眼几乎看不出位置变化但梯度值会被噪声放大如果太大又会让梯度失真。0.5到1m之间是我的经验区间。回溯线搜索的作用是保证每一步目标值不下降同时限制单次移动不超过10m防止无人机在二维平面“乱跳”。3.5 主循环与收敛判定三个子问题的求解器都写好后主循环就很简单了按照“位置→信道→功率”的顺序依次调用并更新。xq, yq L / 2, L / 2 gamma compute_gamma(xq, yq) p np.full(N, Pmax / N) A update_assignment(gamma, p) f_old objective(gamma, A, p) f_history [f_old] for it in range(100): xq, yq, f_new update_position(xq, yq, A, p, f_old) gamma compute_gamma(xq, yq) A update_assignment(gamma, p) p update_power(gamma, A) f_new objective(gamma, A, p) f_history.append(f_new) if abs(f_new - f_old) / max(abs(f_old), 1e-6) 1e-3: print(fconverged at iteration {it}, R {f_new:.4f} bps) break f_old f_new收敛判据用的是相对变化量当总吞吐量相对变化小于0.1%时停止。这个阈值不能设得太小否则BCD会陷入大量的无效迭代但也不能太大否则解还没稳定就提前退出。0.1%到0.5%之间是常用区间。4. 实操过程与仿真结果跑一轮完整迭代要盯住哪些指标4.1 目标函数逐轮变化与收敛判据我用上面的代码完整跑了一遍总迭代上限100轮实际在第38轮触发收敛条件。总吞吐量从初始解的约2.1Mbps提升到收敛时的约3.4Mbps提升幅度大约62%。这个提升看起来不明显在无线资源优化里60%的吞吐量提升已经是相当可观的结果因为它直接意味着同样的频谱资源可以多服务六成用户。观察迭代曲线可以发现前5轮提升最快基本吃掉总提升量的70%以上之后斜率迅速放缓进入“精确打磨”阶段。这是BCD这类一阶/模块替代方法的典型特征初始解离最优解很远时任意一轮调整都能带来巨大收益越接近收敛点各变量之间的耦合效应越明显单块优化的边际收益就越小。所以不要因为后几十轮曲线平缓就误以为算法没在干活那些小步修正往往决定了最终解是不是真的站得住。4.2 多轮优化后的系统表现从位置变化看无人机从500m×500m区域的中心出发逐步飞向用户分布的重心方向但最终停留点并不是几何重心而是偏向那些信道衰落较小、聚合用户数较多的方向。这个结果符合直觉无人机基站不是简单飞往“人多的地方”它要找的是“让整体频谱效率最高的点”。从功率分配看注水算法会把绝大多数功率分配给信道条件好的用户信道差的用户分配到的功率非常小甚至接近0。这在实际系统中可能会导致“公平性”问题但我们的目标函数没有加公平性约束所以这种“嫌贫爱富”的结果是正常现象。如果你的项目需要照顾边缘用户可以把目标函数改成比例公平或者加入最小速率约束这会直接影响功率子问题的结构。从信道分配看贪心匹配总是先把最好的用户-信道对锁定然后逐步补全。8个用户和8个信道一一配对每轮迭代中的配对结果会随着无人机位置的调整而变化但在收敛阶段基本稳定下来不会出现频繁跳变。5. 常见问题与排查技巧实录5.1 不收敛或目标函数来回震荡这是BCD最常遇到的问题但原因往往不在算法本身而在某个子问题的求解器精度不足。比如功率子问题的二分搜索只做了10次迭代就退出求出的功率误差偏大反馈到信道分配后会产生错误的匹配决策进而让位置更新也跟着出错。我的排查顺序固定是先检查每个子问题单独拿出来是否都能收敛再检查主循环的更新顺序和收敛阈值。另一个常见震荡源是位置更新步长太大。如果无人机一步移动几十米信道增益矩阵变化剧烈上一轮算好的功率分配和信道分配瞬间失效下一轮又要重新修正表现就是目标函数上下起伏。解决方法是引入回溯线搜索或者把单步最大移动距离限制在10m以内。5.2 目标函数出现NaN或负值我在调试过程中遇到过一次NaN问题定位后发现log2的参数里冒出了0因为某个用户在某信道上的增益算出来是0。避免的办法就是在计算速率时给log内加一个极小值比如1e-12同时在更新功率时对gamma_active做下限保护。还有一次负速率的情况原因是信道分配矩阵A在初始阶段存在空列而空列对应的功率却是正值。单独看没有逻辑错误但目标函数遍历到空信道时没有用户对它的速率负责导致统计口径不一致。所以在目标函数里加一个显式判断只有p[n] 0且该信道存在用户时才累计速率能避免这类边界情况。5.3 参数调优与初始化技巧BCD对初始化不是特别敏感但不同的起点会收敛到不同的局部最优解这是非凸问题的宿命。我常用两种初始化方式一种是无人机放在区域中心信道分配先用距离最近匹配功率均匀分配另一种是设置多个起点分别跑完再取最优结果。第一种简单稳定适合第一版原型第二种适合追求最终解的工程部署阶段。典型现象可能原因排查与解决目标函数来回震荡位置步长过大或子问题求解精度不足限定单步移动距离回溯线搜索检查二分迭代次数出现NaN或负值信道增益为0或log2参数为0给log内加1e-12对gamma做下限保护某些用户速率长期为0初始化太差或目标函数不公平换成接近用户重心的起始位置或改目标函数无人机位置几乎不动差分步长太小梯度近似被噪声淹没增大差分步长到0.5~1m或改用解析梯度收敛太慢收敛阈值太严格或迭代上限不足阈值放宽到0.1%或观察目标函数是否还在上升6. 扩展思路这套优化框架还能用在哪些地方BCD在无人机通信里的应用远不止上面这种单无人机、单目标的做法。稍微改一改目标函数和约束就能覆盖更广的场景。多无人机协同覆盖是眼下低空物流和应急通信里很热的方向。多架无人机同时在空中组成临时网络时每架飞机的位置、每架飞机的功率分配、用户由哪架飞机服务都是变量。拆块思路完全不变只是位置块从2维变量变成2M维信道分配从单无人机匹配变成多无人机-多用户匹配整体算法框架照样成立。如果让无人机在服务过程中持续移动而不是固定悬停问题就从“找最优位置”变成“设计最优航迹”。航迹优化本质上是一串时间片上的位置决策BCD依然可以把每个时间片当成一个位置块来处理块与块之间增加平滑约束比如最大飞行速度和转弯半径。这个思路在无人机辅助应急通信、灾区临时回传链路等场景中有很实际的应用价值。另外块坐标下降法本身是一个通用框架。OFDM系统的功率与子载波分配、MIMO系统的波束成形矢量设计、边缘计算里的任务卸载与计算资源分配这些问题的共同特征都是“连续变量离散变量强耦合”只要能把变量拆成几块、每块都有可行的子问题求解器你就完全可以复用这套代码结构。最后说一个我自己的习惯凡是第一版原型我从来不用正儿八经的优化求解器全是先手写一个最朴素的BCD把调度逻辑跑通。因为它每个子问题都能单独做单元测试任何一个模块出问题都能立刻定位。等整体收益验证过了再考虑换更高级的求解器也不迟。你如果也准备用这套方法做无人机网络优化建议先从信道增益矩阵和贪心信道分配开始把主循环搭出来再逐步加细节。
返回列表