开头直接进入主题。做无线通信系统级仿真的人肯定都绕不开这个问题多对收发机共用一个频段互相干扰既要让每个链路的速率尽量高又不能让发射功率和电路功耗刷上去。放在MIMO干扰信道里这个能效优化问题天生就是非凸的没有一个现成的凸优化工具箱能直接扔进去求解。这几年圈里最主流的做法就是两条线一条是WMMSE加权最小均方误差及其带半定规划形态的SDP-WMMSE另一条是逐次凸近似SCA。这篇文章就把这两条线从头到尾拆透顺带给出可以直接跑的Matlab代码框架适合正在做无线资源分配、能效优化、或者刚入门MIMO干扰信道建模的研究生和工程师参考。先说清楚我们到底在解什么问题。假设系统里有K个发送端、K个接收端每个发送端配N_t根天线每个接收端配N_r根天线。第k个发送端发给第k个接收端同时给其他K-1个接收端制造干扰。用Q_k表示第k个发送端的发射协方差矩阵维度N_t x N_t半正定。第k个用户的接收信号是y_k H_kk s_k sum_{j ! k} H_kj s_j n_k其中s_k是第k个发送端发射的信号向量协方差为Q_k。H_kk是第k个链路的信道矩阵H_kjj ! k是第j个发送端到第k个接收端的干扰信道矩阵。n_k是零均值循环对称复高斯噪声方差sigma^2。干扰加噪声协方差写出来就是R_k(Q) sigma^2 I_Nr sum_{j ! k} H_kj Q_j H_kj^H由于用户k的有用信号协方差是H_kk Q_k H_kk^H在接收端把干扰当噪声处理的话可达速率就是经典的MIMO信道容量公式R_k(Q) log2 det( I_Nr H_kk Q_k H_kk^H R_k(Q)^{-1} )注意log里的det是真分数规划的分子矩阵里既含Q_k、又含其他用户的Q_j所以R_k(Q)是这些矩阵变量的非凹函数。多个用户速率加起来再加一个分式目标问题就麻烦了。我们要优化的能效目标定义如下EE sum_k R_k(Q) / ( sum_k Tr(Q_k) P_c )其中P_c是电路功耗、基带功耗等与发射功率无关的固定开销Tr(Q_k)是第k个链路的实际发射功率。能效单位是bits/Joule含义是整个网络每消耗一焦耳能量能可靠传输多少比特信息。题目里的“能效优化”就是在给定总功率预算P_total的情况下找到一组Q_k让EE最大同时满足sum_k Tr(Q_k) P_total且Q_k半正定。如果还要考虑单用户公平性可以再加每用户功率约束Tr(Q_k) P_k这取决于你实际仿真的场景加约束不改变后续算法框架。为什么这个优化难解三件事叠在一起分式目标。分子凹非凹、分母凹线性整个目标非凹非凸。R_k(Q)里的log det项本身它不是关于Q的凹函数。把R_k拆成log det(有用干扰噪声)减log det(干扰噪声)前一项是凹函数后一项前面有负号整体丧失凹性。用户之间通过Q_j互相耦合不能逐个发射端单独优化。所以这类问题标准的处理套路是分层外层用Dinkelbach处理分式内层把“给定能效参数eta最大化sum R_k - eta * (sum Tr P_c)”这个加权和问题想办法解出来。Dinkelbach的思想其实特别朴素。我们要最大化f(x)/g(x)等价于解方程f(x) - eta * g(x) 0其中x是最优解、eta是最优能效值。算法第n轮迭代时固定当前的eta_n去求解子问题max_{Q} sum_k R_k(Q) - eta_n * ( sum_k Tr(Q_k) P_c )设解出Q_n然后更新eta_{n1} sum_k R_k(Q_n) / ( sum_k Tr(Q_n) P_c )一直迭代到eta不再上升或|eta_{n1} - eta_n|很小时就认为逼近最优能效。Dinkelbach的理论保证是线性收敛且单调增加这个性质即使在非凸内层问题下也表现出很强的鲁棒性。实测中只要内层不是解得太差外层基本10到20轮就收敛。内层的加权和问题本身还是非凸的这里就需要WMMSE或者SCA出场。先说WMMSE。这个名字全称是Weighted Minimum Mean Square Error它的核心是发现了一个非常漂亮的等价关系在很一般的条件下最大化干扰信道下的和速率等价于最小化加权MSE。为什么绕到MSE因为速率函数是log det的形式求导和迭代都别扭而MSE在给定接收矩阵时可写成关于发射预编码矩阵的二次型优化性质好得多。考虑发射预编码矩阵V_k每用户d个数据流时V_k是N_t x dQ_k V_k V_k^H。线性接收机U_kN_r x d对接收信号做估计s_hat_k U_k^H y_k。这个估计对应的MSE矩阵是E_k (U_k^H H_kk V_k - I)(U_k^H H_kk V_k - I)^H U_k^H R_k U_k给定V_k时最优接收机是MMSE接收机U_k ( H_kk V_k V_k^H H_kk^H R_k )^{-1} H_kk V_k把U_k代回去E_k在最优接收下就变成E_k I - V_k^H H_kk^H ( H_kk V_k V_k^H H_kk^H R_k )^{-1} H_kk V_k接着定义权重矩阵W_k E_k^{-1}。核心等价关系是log det( I H_kk Q_k H_kk^H R_k^{-1} ) max_{U_k, W_k} [ log det(W_k) - Tr(W_k E_k) d ]所以最大化sum R_k可以写成对一个很复杂的MSE表达式的最大化。通过交替优化U_k、W_k和V_k逼近局部最优这就是WMMSE的标准迭代格式。迭代三步走第一步固定V_k更新U_k为MMSE接收矩阵。 第二步用新U_k算E_k令W_k E_k^{-1}。 第三步固定U_k和W_k求解所有V_k。当目标函数带上能效项-eta * (sum Tr P_c)时第三步的优化变成总功率惩罚下的加权MSE最小化典型解是V_k_new ( sum_j H_jk^H U_j W_j U_j^H H_jk mu_k I )^{-1} H_kk^H U_k W_k其中mu_k是拉格朗日乘子用来满足第k个用户或者全局功率约束。mu_k没法闭式解要用二分法去搜。这个公式在实际代码里很容易写错因为求和遍历所有j不能只算本链路这个sum_j H_jk^H U_j W_j U_j^H H_jk把所有用户对第k个发射端的干扰路径都汇总进来了漏掉任何一个都会导致结果发散。那SDP-WMMSE又是什么你在很多论文里看到SDP-WMMSE这个名字其实就是把第三步或者整个WMMSE发射端更新用半定规划来表达。因为Tr(W_k E_k)关于V_k并不是凸的直接对V_k闭式更新时mu_k的搜索比较麻烦再加上你可能想加更多约束比如每用户最大功率、秩约束或者是协方差矩阵整体结构这时候直接用CVX把整个更新建模成半定规划就舒服得多。核心变量从V_k变成Q_k约束Q_k半正定目标函数里Tr(W_k E_k)的部分要做等价改写然后把所有凸项交给SDP求解器处理。这种做法的代价是每一轮迭代都要调用一次CVX仿真时间会明显变长但稳定性高尤其适合调试初始版本。我个人建议如果只是验证算法趋势先用闭式WMMSE如果要在复杂约束下对比界再上SDP-WMMSE。SCA的思路跟WMMSE完全不同。它不绕到MSE直接对速率函数做逐次凸近似。前面说了R_k(Q) log det(有用干扰噪声) - log det(干扰噪声)。记F_k(Q) log det( sigma^2 I sum_j H_kj Q_j H_kj^H ) G_k(Q) log det( sigma^2 I sum_{j ! k} H_kj Q_j H_kj^H )那么R_k(Q) F_k(Q) - G_k(Q)。问题在于F_k是凹的G_k也是凹的凹减凹整体不再是凹。SCA的处理方式是在每次迭代的当前点Q^{(t)}把G_k做一阶泰勒展开。因为G_k是凹函数一阶泰勒展开是所有点的上界G_k(Q) G_k(Q^{(t)}) sum_{j ! k} Tr( A_{k,j}^{(t)} ( Q_j - Q_j^{(t)} ) )其中A_{k,j}^{(t)} H_kj^H ( sigma^2 I sum_{l ! k} H_kl Q_l^{(t)} H_kl^H )^{-1} H_kj把这个上界代入R_k得到当前点的凹下界R_k_tilde(Q) F_k(Q) - G_k(Q^{(t)}) - sum_{j ! k} Tr( A_{k,j}^{(t)} ( Q_j - Q_j^{(t)} ) )这个式子对Q是凹的因为F_k是凹的减去的只是常数项和线性项。于是把每个用户的凹下界加起来再减去能效惩罚项就是当前SCA迭代要最大化的凸问题。用CVX直接建问题变成标准的半定规划max sum_k R_k_tilde(Q) - eta * ( sum_k Tr(Q_k) P_c ) s.t. sum_k Tr(Q_k) P_total Q_k 0每次解出来一个新的Q作为下一轮的线性化点不断迭代直到目标值不再明显上升。SCA收敛性有较完整的理论支撑只要每一步解出全局最优的凹近似子问题生成的序列会收敛到原问题的稳定点。这里需要注意一个细节SCA里的“凹下界”是逐点构造的所以每次迭代点变了A_{k,j}^{(t)}就要重算不能拿上一轮的下界直接再优化那就成了一次性近似的解不是SCA。WMMSE和SCA哪个好这个问题没有标准答案。我个人的体会是WMMSE的收敛速度通常更快特别在信噪比不算太高、约束比较简单时闭式更新的几百次迭代在Matlab里也就是毫秒量级。SCA的优势在于问题拓展性强能效、感知、安全速率、谐波干扰这些带复杂约束的目标函数往里边加约束只需要改CVX模型而且半定规划的形式让我们可以顺便观察发射协方差矩阵的结构。另外当每用户多流且发射端天线数较多时SCA对初值的敏感度通常比WMMSE更低不容易飞到奇点附近。下面进入Matlab实现部分。先说参数配置实际仿真里可以按下面这个基底来。clear; clc; rng(2024); K 3; % 用户对数 Nt 4; % 发射天线数 Nr 2; % 接收天线数 ds min(Nt,Nr); % 每用户数据流数 SNR_dB 10; sigma2 1; P_total 10^(SNR_dB/10); % 总功率约束 Pc 5; % 电路及固定功耗 outer_max 20; inner_max 100; tol_eta 1e-4;信道生成用瑞利平坦衰落每个元素独立同分布复高斯。注意这里H是一个K x K的cell矩阵H{k,j}表示第j个发射端到第k个接收端的信道尺寸Nr x Nt。多小区多用户场景里这个命名顺序最容易乱建议一开始就统一好后面写任何公式都按这个顺序。H cell(K,K); for k 1:K for j 1:K H{k,j} (randn(Nr,Nt) 1i*randn(Nr,Nt)) / sqrt(2); end end初始化V时不要所有用户都用同一个功率分配最简单可靠的是均匀分配总功率。每个V_k取随机酉矩阵乘以sqrt(P_total / (K*ds))保证初始时总功率不超过约束。WMMSE核心三步更新函数可以这么组织。先写接收矩阵和权重更新:function [U,W] update_UW(V, H, sigma2, K) d size(V{1}, 2); U cell(K,1); W cell(K,1); for k 1:K Rint sigma2 * eye(size(H{1,1},1)); for j 1:K if j k, continue; end Rint Rint H{k,j} * (V{j} * V{j}) * H{k,j}; end A H{k,k} * V{k}; U{k} (A * A Rint) \ A; E (U{k} * A - eye(d)) * (U{k} * A - eye(d)) U{k} * Rint * U{k}; E 0.5 * (E E); W{k} inv(E); end end然后发射预编码更新。这一步把能效惩罚项eta加进去。标准的闭式更新里矩阵B_k sum_j H_jk^H U_j W_j U_j^H H_jk是半正定的我在B_k上加的是mu_k * I其中mu_k是拉格朗日乘子。求解V_k后如果Tr(V_k V_k^H)超过P_total/K就增大mu_k低于则减小mu_k用二分法收敛。function V update_V(U, W, V, H, sigma2, P_per_user, eta) K length(V); Nt size(V{1},1); d size(V{1},2); V_new cell(K,1); for k 1:K Bk zeros(Nt,Nt); for j 1:K Bk Bk H{j,k} * U{j} * W{j} * U{j} * H{j,k}; end Ak H{k,k} * U{k} * W{k}; mu_l 0; mu_h 1; v_opt zeros(Nt,d); for step 1:30 mu (mu_l mu_h) / 2; v_opt (Bk (mu eta) * eye(Nt)) \ Ak; if trace(v_opt * v_opt) P_per_user mu_l mu; else mu_h mu; end end V_new{k} v_opt; end V V_new; end注意我在括号里写的是mu eta因为能效惩罚项在目标里是-eta * sum(trace(V_j V_j^H))带入KKT条件后等效于在功率矩阵上加了(eta)倍的I但这个系数不一定严格出现在该位置取决于目标函数有无xi系数。如果只要总功率约束不要每用户约束可以把P_per_user换成P_totalmu共享。为了让二分搜出来的V_k功率接近上界但不超过收敛条件不能只看最大迭代次数最好是检查power和P_per_user的相对误差小于1e-6。实际用的时候我会把二分迭代数调到50否则mu可能没收敛到位V_k功率会明显小于上界速率损失很大。外层Dinkelbach循环拼接起来V cell(K,1); for k 1:K V{k} sqrt(P_total / (K*ds)) * (randn(Nt,ds) 1i*randn(Nt,ds)) / sqrt(2); end eta 0; EE_record []; for outer 1:outer_max for inner 1:inner_max [U,W] update_UW(V, H, sigma2, K); V update_V(U, W, V, H, sigma2, P_total/K, eta); end R calc_sum_rate(V, H, sigma2); P calc_power(V) Pc; eta_new R / P; EE_record [EE_record eta_new]; if abs(eta_new - eta) tol_eta, break; end eta eta_new; endsum_rate和power函数比较基础就是按公式循环求log2 det。这里有一个工程上的隐藏坑矩阵E_k求逆前必须做对称化E (E E)/2否则由于浮点误差可能出现虚部极小但非零的情况导致inv结果略偏。虽然对最终结果影响很小但在调试时看不出问题算det时却会蹦出微小虚部画图时轻则警告重则NaN。SCA版本用CVX建。注意要用CVX的同学先确认solver装全了Solve是SeDuMi还是SDPT3都行小规模问题差别不大。SCA的Matlab核心部分如下function [Q_opt, obj] sca_inner(Q_cur, H, sigma2, P_total, eta) K size(H,1); Nt size(Q_cur{1},1); Nr size(H{1,1},1); obj_tilde 0; % 预计算线性化矩阵 A_{k,j} A_lin cell(K,K); for k 1:K Rint sigma2 * eye(Nr); for j 1:K if j k, continue; end Rint Rint H{k,j} * Q_cur{j} * H{k,j}; end Rint_inv inv(Rint); for j 1:K if j k, continue; end A_lin{k,j} H{k,j} * Rint_inv * H{k,j}; end end cvx_begin sdp quiet variable Q(Nt, Nt, K) hermitian semidefinite expr 0; for k 1:K Rtot sigma2 * eye(Nr); for j 1:K Rtot Rtot H{k,j} * Q(:,:,j) * H{k,j}; end logdet_Rtot log_det(Rtot); logdet_Rint_cur 0; Rint_cur sigma2 * eye(Nr); for j 1:K if j k, continue; end Rint_cur Rint_cur H{k,j} * Q_cur{j} * H{k,j}; end logdet_Rint_cur log_det(Rint_cur); lin_part 0; for j 1:K if j k, continue; end lin_part lin_part trace( A_lin{k,j} * ( Q(:,:,j) - Q_cur{j} ) ); end expr expr logdet_Rtot - logdet_Rint_cur - lin_part; end maximize( expr - eta * (sum_Tr_Q Pc) ) subject to sum_Tr_Q 0; for k 1:K sum_Tr_Q sum_Tr_Q trace(Q(:,:,k)); end sum_Tr_Q P_total; cvx_end for k 1:K Q_opt{k} Q(:,:,k); end obj cvx_optval; end写这段有一个点要注意代码里sum_Tr_Q这种变量必须在constraint之前用表达式定义CVX要求约束里出现的中间表达式是凸表达式。另外log_det在CVX里要求矩阵半正定Rtot包含噪声项sigma2 I理论上一定是正定的但数值上极小的特征值可能导致log_det报错。稳妥做法是在Rtot里加一个微小正则项比如1e-10 * eye(Nr)。这不会改变结果但能避免偶尔的数值崩溃。SCA的外层同样套Dinkelbach循环每次把当前Q作为下次线性化点。这里我习惯每轮SCA只内迭代5到10次因为Dinkelbach外层迭代会在不同eta下反复调用内层内层迭代太多总时间爆炸。实测下来SCA内层5次加上外层20次就已经能收敛到不错的结果继续加次数只是小数点后四五位的差别。关于仿真结果建议画三样东西能效随Dinkelbach迭代次数变化的收敛曲线、对比不同发射功率下的能效曲线、以及各用户速率组成柱状图。能效曲线如果出现先上冲后下降再回升的锯齿通常是内层没有收敛就更新了eta导致的可以加大inner_max或者把Dinkelbach外层收敛阈值放宽。如果出现单调上升但很久不收敛检查外层是否忘记重新初始化内层变量也就是说外层每轮迭代都要紧接着上次的V/Q继续迭代而不是清零重来这个细节很多人踩坑。仿真中我总结过一份高频问题清单供排查用第一个是NaN和Inf。要么是inv(P)遇上奇异矩阵要么是log_det里出现负数特征值。应对办法Rint加对角正则sigma2 * 1e-8 * eye(Nr)E_k求逆前做对称化Q更新后强制投影到半正定锥简单做法是[V,D] eig(Q); D(D0)0; Q VDV。这种投影在理论上有争议但在工程仿真里能显著提升鲁棒性作为调试手段完全够用。第二个是收敛到无意义的低能效。最常见的原因是发射预编码的功率分配在二分法下没有解到位导致实际功率远低于约束上限。解决方式是检查总功率P_actual和约束P_total的大小关系如果差了0.1以上说明mu搜索太粗糙或初始mu范围太小。我把mu_h初始值从1改到1000后这个现象基本消失。第三个是SDP-WMMSE里CVX求解时间太长。小规模问题K3, Nt4, Nr2每次CVX调用也就几十毫秒一旦K到7以上或者天线数到8CVX迭代次数指数增长。这时候建议只保留SCA里单变量交替更新的步骤或者是用MOSEK这类商用solver代替SeDuMi。如果只是做算法验证K不超过5就够了。第四个也是最后一个关于发射协方差的秩。SDR把Q的半正定约束保留、但丢掉了秩约束因此SDP解出来的Q_k很可能不是秩一也就是发射多流时实际需要多个数据流。这在某些文献里算作松弛上界不能直接当最终波束赋形使用。实际工程中要么直接把最大特征值对应的特征向量拿出来当单流预编码要么在目标里加一个小迹惩罚项推着解向低秩靠。前者简单后者在论文里更常见但调惩罚系数的过程比较费劲。我个人做这个小课题最大的感受是这类算法最大的风险不在数学推导而在Matlab数值细节。WMMSE和SCA的理论都非常成熟论文里的伪码三行就能写完可是落到代码执行时矩阵维数、共轭转置、循环索引、正则项设置每一个地方都能让结果悄悄走样。调试时不要一上来就跑完整迭代先把用户数设为2、天线数设为2、把SDR和闭式WMMSE的结果对比确认它们收敛到同一个目标值以后再加复杂度和规模。如果你后续想扩展这个框架可以直接改造成鲁棒能效优化、智能反射面辅助MIMO能效优化或者IRS与波束成形联合设计。核心变化无非是在目标或者约束里多几个矩阵变量SCA这边多几组线性化项WMMSE这边多几组更新模块整体的Dinkelbach加交替优化骨架完全不用动。最后再分享一个小技巧仿真脚本里把所有关键中间量都做成可选项输出比如每一轮的sum_rate、power、eta、U的范数、W的条件数。这样跑参数扫描的时候不用反复打断重跑就能定位是哪一步开始异常。我在初始版本里只输出最终能效值结果发现算法在某些参数下悄悄崩了查半天才靠增加中间输出日志定位到是W矩阵条件数过大。做优化仿真日志和可视化比什么都重要。