简介面向计算机、电子信息工程与数学专业学生这套MATLAB代码给出了最大化智能反射面IRS辅助MIMO系统保密率的坐标下降算法实现适用于课程设计、期末大作业与毕业设计等场景。资源为9KB的zip压缩包共12个文件核心为6个m脚本覆盖主程序、IRS信道建模、固定相位与一般特征值分解下的容量计算、坐标下降迭代更新及穷举搜索对比等模块同时配有说明文档、补充案例数据压缩包和备份文件可直接在MATLAB 2014/2019a/2021a中运行。代码采用参数化编程关键参数便于修改注释清晰既能帮助理解交替优化思路也能通过调整参数观察不同配置下的保密率差异加深对物理层安全与MIMO系统优化的认识。目前已有37人学习下载适合需要将理论推导与代码实现相结合的高年级本科生及研究生参考学习。1. 从IRS到坐标下降为什么保密率优化绕不开逐坐标调整IRS辅助的MIMO系统里反射单元动辄几十上百个每个单元都要调制一个单位模量的相移发射端还有多根天线要分配功率。如果直接把所有变量堆在一起上内点法非凸、变量多非常容易陷进某个糟糕的局部解。坐标下降算法是这类问题的天然选择——固定其他变量只更新当前一个变量把高维问题拆成多个一维优化。最大化IRS辅助MIMO系统保密率的坐标下降算法在MATLAB里用几十行脚本就能验证反射单元从8个加到64个保密率随信噪比上升的趋势一目了然。下面把信道建模、目标函数、迭代框架和MATLAB实现完整拆开照着一行行敲就能跑。2. 先把信道和保密率公式写进MATLAB从三跳信道到可分式目标2.1 IRS-MIMO的三段信道基站、反射面和两个终端完整的IRS辅助MIMO链路可以切分成三段。第一段是基站到IRS的MIMO信道维度是N乘M记作G_BIN是反射单元数M是基站天线数。第二段是IRS到合法用户和窃听者的两个信道分别记作h_RU和h_RE维度都是1乘N。第三段是基站直接到两个终端的直射信道分别记作h_DU和h_DE维度都是1乘M。直射信道一定要保留。IRS不会把来波全部吸收它只是叠加了一个受控的反射分量直射分量在实际场景里永远存在。很多刚接触这个方向的人为了简化仿真直接把直射信道设为零这样在高信噪比下得到的保密率会虚高尤其在有视距的场景里直射链路往往比反射链路更强。代码里我习惯这样生成信道% 信道参数与随机生成 M 4; % 基站天线数 N 16; % IRS反射单元数 pl_ru 2; % 用户IRS链路相对增益系数 pl_re 0.5; % 窃听IRS链路相对增益系数 % 复高斯随机信道每个系数 ~ CN(0,1) G_BI (randn(N,M) 1i*randn(N,M))/sqrt(2); h_RU sqrt(pl_ru) * (randn(1,N) 1i*randn(1,N))/sqrt(2); h_RE sqrt(pl_re) * (randn(1,N) 1i*randn(1,N))/sqrt(2); h_DU sqrt(0.8) * (randn(1,M) 1i*randn(1,M))/sqrt(2); h_DE sqrt(0.3) * (randn(1,M) 1i*randn(1,M))/sqrt(2);这里有一个容易被忽略的约定所有到终端的信道我都用行向量表示发射预编码w用列向量表示这样h_u * w直接是一个标量增益不用在代码里到处做共轭转置。涉及多天线的矩阵运算时这个约定能省掉一大堆.*维度匹配的报错。复高斯随机抽样后除sqrt(2)是为了让每个系数的平均功率归一化为1。pl_ru和pl_re是人为加上的链路增益系数相当于把用户路径损耗设成比窃听者低6dB左右。为什么要这样做如果用户链路和窃听链路的平均功率相同保密率最优解可能长期落在零附近坐标下降起点太低不容易观察到算法效果。仿真初期不妨让用户链路略占优势先把算法跑通再去调整成信道对称的极端场景。2.2 保密率目标函数为什么不能只把信噪比拉满保密率的定义是合法链路可达速率减去窃听链路可达速率外部取正。写成MATLAB能直接计算的标量形式R_s log2(1 p |h_u * w|^2 / sigma_u2) - log2(1 p |h_e * w|^2 / sigma_e2)其中等效信道h_u和h_e分别代表包含直射和IRS反射在内的复合信道h_u h_DU h_RU * diag(theta) * G_BI h_e h_DE h_RE * diag(theta) * G_BI这两个式子看完就能明白为什么坐标下降可行theta的第n个元素只出现在diag(theta)的第n个对角项上而diag(theta)*G_BI的每一行又只和h_RU、h_RE的第n个元素相乘。换句话说固定其他相移后当前第n个相移对保密率的影响是独立的一块可以单独做一维优化。注意保密率公式里我用了原始的差值而不是先取正。优化阶段如果直接写成max(0, R_s)当R_s小于零时梯度为零坐标下降会停在一个毫无意义的平台再也走不出来。实际代码里优化目标保留差值显示性能时再取正这一点非常重要。为什么不直接最大化合法用户信噪比因为那样做很可能让窃听者也获得更高信噪比差值反而变小。保密率的目标是拉开两条链路的差距不是单点最优。这也决定了功率分配不是简单的注水发射预编码也不是单纯的匹配滤波器。2.3 向量化计算保密率MATLAB里别写愚蠢的累加循环写保密率函数时新手容易这样先把h_u置零然后写一个for n 1:N的循环把每个反射单元的贡献累加进去。这样写没有错但是代码又臭又长而且后面做坐标下降时你想把某个相位拿出来修改变量这个累加逻辑会逼着你重复写很多遍等效信道构造。我一般把保密率计算封装成一个函数参数就是theta、w、p和各种信道矩阵function R calc_secrecy(theta, w, p, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2) h_u h_DU h_RU * diag(theta) * G_BI; h_e h_DE h_RE * diag(theta) * G_BI; S_u abs(h_u * w)^2 * p; S_e abs(h_e * w)^2 * p; R log2(1 S_u/sigma_u2) - log2(1 S_e/sigma_e2); end这段代码最大的好处是矩阵维度清晰。h_RU是1乘Ndiag(theta)是N乘NG_BI是N乘M乘出来正好是1乘M再和h_DU相加得到合法用户等效信道。如果不想用diag构造稀疏大矩阵也可以用h_RU .* theta.先对行向量按元素乘再去乘G_BI。效果一样N大时省一点内存但可读性不如diag。在线跑小例子用diag没问题真正做256个反射单元的蒙特卡洛时再考虑优化。把这个函数放到单独m文件里后面的坐标下降函数和主脚本都能反复调用。有一次我在别人代码里看到保密率计算直接嵌在主循环里坐标下降每试探一个候选相位就要重算一次h_u和h_e代码复制了三大段改参数时漏改一处就出bug。封装之后至少后面几个坑能少踩一半。3. 坐标下降怎么落地交替更新反射相位与发射预编码3.1 交替迭代框架外层换w内层扫theta坐标下降不是只能处理单个变量而是可以把变量分成块块之间交替优化。我们这里把发射预编码w、功率p和IRS相位theta当成两大块。外层先固定theta更新w和p内层固定w和p对theta逐坐标更新。这样整体还是坐标下降的骨架只不过内层又加了一层逐坐标扫描。给定theta后等效信道h_u被完全确定。这时候发射方向我一般选MRT也就是把w设成h_u的共轭方向w conj(h_u) / norm(h_u)为什么不是直接做广义特征分解因为保密率目标不是简单的信噪比最大化广义特征分解只能让两个增益比值最大和带两个log的保密率并不严格等价。MRT闭式、稳定、跑得快作为坐标下降的第一版最合适等算法整体收敛后再换更好的方向。这属于工程上的渐进式改进不算偷懒。功率p不要和w一起线性缩放否则w归一化就没有意义。固定w方向后p在0到P之间一维搜索。MATLAB自带的fminbnd够用function [w, p] update_wp(theta, P, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2) h_u h_DU h_RU * diag(theta) * G_BI; h_e h_DE h_RE * diag(theta) * G_BI; w h_u / norm(h_u); obj (pval) -calc_secrecy(theta, w, pval, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); [p, ~] fminbnd(obj, 1e-6, P); endfminbnd的搜索区间下限我给了1e-6而不是0。理论上p0也合法但p0时保密率恒为零任何小的试探都会让搜索器认为找到了边界返回一个没有意义的功率。加一个极小值下限后搜索至少在一个非退化的区间里进行。上限P就是系统总功率约束。3.2 固定(w,p)后逐坐标更新反射相移把0到2π打散成一维搜索这是整个坐标下降算法里最有意思的部分。固定w和p之后保密率R_s只与theta有关。固定其他theta_j不变改变theta_n时h_u和h_e都只有一行中的一项在动所以保密率关于theta_n是一个一维函数。理论上可以对这个一维函数做解析优化把theta_n写成atan2的形式。我试过只考虑用户链路的做法确实能化简出一个相位差公式但把窃听链路也放进去后两个log分子分母都有theta_natan2不再成立强行求导也得不到闭式解。保号后的一维函数还可能出现多峰梯度法容易陷在局部峰里。所以工程上我采用网格搜索把0到2π均匀离散成L个点逐一计算保密率取最大值对应的相位作为当前坐标的更新。这其实就是坐标下降里最常见的一种精确线搜索实现只是把连续区间离散化了。function theta coord_desc_theta(theta, w, p, L, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2) N size(G_BI, 1); phi_grid (0:L-1) * (2*pi/L); for n 1:N best_R -inf; best_phi angle(theta(n)); for idx 1:L theta_probe theta; theta_probe(n) exp(1i * phi_grid(idx)); R_probe calc_secrecy(theta_probe, w, p, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); if R_probe best_R best_R R_probe; best_phi phi_grid(idx); end end theta(n) exp(1i * best_phi); end end这段代码里有两个关键点。第一每次试探都用theta_probe theta而theta变量在循环中已经被前面坐标更新过所以这是严格的顺序坐标下降不是并行的同时更新。第二候选相位来自预先定义的phi_grid搜索到的最优相位立即写入theta(n)下一轮第n1个坐标看到的已经是新值。复杂度上N个反射单元每个单元L次保密率计算总共N乘L次。N16、L120时一次内层扫描大约1920次保密率计算。每次保密率里又包含两次diag矩阵乘法在普通笔记本上大概零点几秒。如果你用在线版MATLAB或者把N放大到64以上就可以考虑把diag(theta)*G_BI替换成bsxfun(times, G_BI, theta)或者theta .* G_BI这样的按行广播能省掉对角矩阵的乘法开销。3.3 停止准则与迭代上限保密率增量、相位漂移和最大循环外层交替迭代至少要设三个停止条件。第一是保密率增量也就是本轮和上一轮的R_s差值小于某个容忍度比如1e-3 bit/s/Hz。第二是相位向量的变化量norm(theta_new - theta_old)小于1e-2。第三是最大迭代上限maxIter防止算法在两个候选相位之间反复横跳。关于相位漂移有一个容易踩的坑theta是复数角度0和2π在数值上代表同一个相位但直接做差会得到2π的虚假漂移。好在我的坐标更新把候选角都限制在[0, 2π)区间内更新后angle函数也落在同一区间所以实际没遇到这个问题。如果你自己写连续角度优化就要把相位差值映射回(-π, π]再比较。外层框架我一般这么写theta exp(1i * 2*pi*rand(N,1)); w zeros(M,1); p 0; R_hist zeros(maxIter,1); for iter 1:maxIter [w, p] update_wp(theta, P, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); for inner 1:5 theta coord_desc_theta(theta, w, p, L, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); end R_hist(iter) calc_secrecy(theta, w, p, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); if iter 1 abs(R_hist(iter) - R_hist(iter-1)) tol break; end end内层固定w和p的时候多跑几轮theta扫描能保证theta在本次外层迭代中尽量收敛。否则只扫描一轮theta还没稳定就去更新w整体收敛曲线会出现毛刺。我试过内层只跑1轮最后的保密率序列会有明显波动但最终结果相差不大。想要文章里那种漂亮平滑的收敛曲线内层5轮就够。这套算法每轮目标函数都不会下降因为w、p和theta的每一步更新都是朝着至少不坏的方向走。所以理论上R_hist应该是单调不减的。如果看到下降多半是记录时机错了下一章会专门讲这个问题。4. 复现一份完整MATLAB脚本参数表、主循环和曲线输出4.1 参数表天线数、反射单元数、信噪比区域怎么选参数不是随便拍的。反射单元太少比如N4曲线太平坦看不出IRS优势N太多比如N128坐标下降每个内层要扫1万多次保密率教学实验会很慢。我建议第一次复现时用N16迭代30轮整个脚本运行时间控制在十秒以内。以下是一份自测过的参数表可以直接照抄参数推荐值作用说明M4基站天线数越大阵列增益越高N16IRS反射单元数越多相位自由度越大P_dB0 到 20发射功率范围dB刻度sigma_u21合法用户噪声功率sigma_e21窃听者噪声功率L120每个坐标一维搜索的离散点数maxIter30外层交替迭代上限tol1e-3保密率增量收敛门限pl_ru / pl_re2 / 0.5反射链路相对增益调节用户优势噪声功率都取1发射功率P_dB从0开始这意味着初始信噪比从0dB附近起步保密率曲线可以看得很清楚。如果你想让合法用户链路优势更大可以把pl_ru提高到3pl_re压到0.3。但注意不要一个参数试到底不同信道增益下坐标下降的收敛速度和最终保密率都差不少。4.2 主脚本框架先随机信道再交替求解把第二章和第三章的函数组合起来就是一个完整的主脚本。这个脚本可以直接新建脚本文件按F5运行clear; clc; rng(2025); % 固定随机种子保证可复现 M 4; N 16; P_dB 10; P 10^(P_dB/10); sigma_u2 1; sigma_e2 1; maxIter 30; tol 1e-3; L 120; pl_ru 2; pl_re 0.5; G_BI (randn(N,M) 1i*randn(N,M))/sqrt(2); h_RU sqrt(pl_ru) * (randn(1,N) 1i*randn(1,N))/sqrt(2); h_RE sqrt(pl_re) * (randn(1,N) 1i*randn(1,N))/sqrt(2); h_DU sqrt(0.8) * (randn(1,M) 1i*randn(1,M))/sqrt(2); h_DE sqrt(0.3) * (randn(1,M) 1i*randn(1,M))/sqrt(2); theta exp(1i * 2*pi*rand(N,1)); w zeros(M,1); p 0; R_hist zeros(maxIter,1); for iter 1:maxIter [w, p] update_wp(theta, P, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); for inner 1:5 theta coord_desc_theta(theta, w, p, L, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); end R_hist(iter) calc_secrecy(theta, w, p, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); if iter 1 abs(R_hist(iter) - R_hist(iter-1)) tol fprintf(收敛于第%d轮\n, iter); break; end end figure; plot(R_hist(1:iter), -o); xlabel(迭代次数); ylabel(保密率 (bit/s/Hz)); title(坐标下降收敛曲线); grid on;运行完你会看到一条单调上升的曲线前几轮涨得很快后面逐渐趋于平缓。这就是坐标下降在两个变量块之间交替前进的特征。如果曲线出现了突然下降优先检查R_hist记录的位置是不是在更新w之后、更新theta之前——那个位置R值一定会变乱。4.3 后处理与曲线复现保密率随噪声功率和反射单元数量变化手工复现一张论文级曲线时需要扫两个维度一是发射功率P_dB二是反射单元数量N。把P_dB从0取到20步长2dBN取8、16、32三条线。每条线都先跑上面的主脚本记录收敛后的R_s最后用loglog或semilogx画出来。注意每次改变N都要重新生成信道rng种子也要重新设置。为了对比公平可以在外层加一个蒙特卡洛循环num_realizations 100; R_avg zeros(length(P_dB_set), length(N_set)); for i 1:length(N_set) N N_set(i); for j 1:length(P_dB_set) P 10^(P_dB_set(j)/10); sum_R 0; for r 1:num_realizations rng(r * 100 i * 10 j); % 重新生成信道和初始化theta ... sum_R sum_R final_R; end R_avg(i, j) sum_R / num_realizations; end end这里随机种子的编号规则不是唯一的但目的是让每个(N, P)组合在每次蒙特卡洛中都有一组独立信道同时不同组合之间可以重复对照。不要偷懒用同一个rng(0)跑完全部否则所有曲线都来自同一组信道N增加纯粹是相位自由度变多信道随机性完全丢失曲线会异常平滑不代表真实性能。5. 避坑指南这5个问题最容易让保密率迭代翻车5.1 MATLAB中文注释乱码脚本一开就废现象用MATLAB 2023或2022版本打开别人分享的脚本中文注释全部变成方框或乱码代码本身还能运行但完全没法排查参数。原因MATLAB从R2020b之后在Windows上默认文件编码从GBK切换到UTF-8而早期脚本保存为GBK反过来用新版本保存UTF-8脚本到老版本环境也可能乱码。这个问题在中文技术群里几乎每周都有人问。解决统一在MATLAB预设项里找到“常规—MAT文件—字符编码设置”把默认编码改成UTF-8。如果你只负责运行脚本直接删除注释里的中文全部改成英文注释最省事。我在分享可复现代码时函数文档都写英文只有在文章正文里详细解释中文含义避免编码问题干扰复用。5.2 保密率曲线明明单调递增某轮突然下跌现象主脚本跑完R_hist曲线在第5轮以后出现一个明显的向下尖峰然后又继续上升。原因记录时机不对。先执行update_wp(theta,...)再计算R_hist(iter)等于用新w和新p计算保密率而上一轮的R_hist用的是旧w和旧p两个不在同一组变量上自然可能下降。解决把保密率记录放到内层theta收敛之后且下一次迭代更新w之前。也就是先记录再更新w。或者干脆把记录放在每个外层迭代末尾代码顺序改成先计算R_hist再更新w。这样每一步坐标下降都至少不会变差曲线单调性就回来了。5.3 一维功率搜索返回极端值保密率直接为负现象fminbnd搜索p时返回结果要么是1e-6要么是P上限整条曲线保密率接近零。原因p在低信噪比区域目标函数对p不敏感fminbnd的黄金分割搜索把p推到边界。下限为0时更加严重p0导致目标函数恒零搜索器判定边界就是最优。解决下限改1e-6上限保持P。如果发现搜索结果总是贴着上下限说明w方向给得不合适先检查MRT方向是不是对着h_u而不是h_e。建议在update_wp函数里临时打印h_u和h_e的范数肉眼确认合法通道更强。很多所谓保密率翻车其实是信道增益配置反了。5.4 网格搜索相位太粗最优解差了0.5dB现象L20时保密率曲线比文献参考值低0.3到0.5dB调大L又觉得慢。原因L20对应18度一个步长最优相位落在两个网格点之间时误差最大可达9度这对IRS波束成形是明显损失尤其是N较大时各个反射单元的小误差会累加。解决第一轮用L120粗筛锁定最优相位所在区间第二轮在最优相位附近加减一个步长用更细的网格或直接调用fminbnd在这个小区间内细搜。下面是细搜的替代写法phi_best best_phi; phi_range [phi_best - 2*pi/L, phi_best 2*pi/L]; phi_refined fminbnd((phi_probe) -calc_secrecy_with_angle(phi_probe), phi_range(1), phi_range(2));注意calc_secrecy_with_angle要把临时相位写回theta的第n个位置再调用保密率函数。这个是纯网格搜索的外部包装仍然符合坐标下降的框架只是每个坐标的线搜索更精确了一点。5.5 不同N的保密率曲线交叉得乱七八糟现象画N8、16、32三条曲线时N8的个别点反而比N32还高曲线交叉图表没法用。原因每次循环跑完一组信道后没有重置随机种子或者用同一个种子导致N8和N32时信道样本不同比较的是不同信道下的点自然没有单调性。解决把所有待比较场景的随机种子固定在同一组编号规则下比如rng(100*N P_dB)这种方式确保每次重新生成信道都得到相同随机序列。更严谨的做法是蒙特卡洛求平均每个N和P组合都各自跑100次独立信道取平均再用平均保密率画曲线。这样单次随机性被抹平曲线才光滑可复现。6. 让坐标下降再快一截随机坐标、冷启动与离散化技巧到手可用之后可以试着把算法改得更快。第一个技巧是随机坐标选择。当N很大顺序扫描所有反射单元耗时太长每一轮随机选择一部分坐标来更新比如每次只更新三分之一其他坐标保持不动。这在统计意义上仍然能收敛而且对变量之间的强相关性有一定的“破坏”作用有时反而跳出局部解。不过随机选择会让保密率曲线带上噪声收敛判据建议改成对最近5轮取滑动平均而不是看单轮增量。第二个技巧是冷启动。把坐标下降里的theta初始相位从0到2π随机取这没错但随机种子对最终结果影响很大。我一般跑3到5个随机初始theta取收敛后保密率最高的那一个作为最终答案。运气好的时候一个糟糕的初始相位会让坐标下降卡在局部解多试几次成本很低收益却很实在。第三个技巧是支持离散相位约束。真实IRS相移通常只有1位或2位精度比如只能取{0, π}或者{0, π/2, π, 3π/2}。你不需要改坐标下降框架只要把网格搜索的候选点从连续均匀网格换成离散相位集合代码几乎不用动。离散化之后每个坐标的搜索范围变小运行速度更快但保密率会有损失典型值在0.5到1.5dB之间。如果你的目标是工程落地这一步必须做。验证算法是否真的work我有个习惯先画一条随机相位下的保密率作为baseline再画坐标下降收敛后的点对比拉升了多少。如果一个系统中随机相位和优化后相位几乎没有差别说明问题不在算法而在信道模型——可能窃听者信道太强或者相位自由度根本没有生效。这时候不要急着调算法回头检查信道增益系数。我在跑这类仿真时习惯先把N设成8跑通再逐步加大否则收敛慢的因素和信道因素混在一起定位问题很痛苦。希望帮到你。本文还有配套的精品资源点击获取 IRS辅助MIMO保密率优化：坐标下降算法原理与MATLAB实战
简介面向计算机、电子信息工程与数学专业学生这套MATLAB代码给出了最大化智能反射面IRS辅助MIMO系统保密率的坐标下降算法实现适用于课程设计、期末大作业与毕业设计等场景。资源为9KB的zip压缩包共12个文件核心为6个m脚本覆盖主程序、IRS信道建模、固定相位与一般特征值分解下的容量计算、坐标下降迭代更新及穷举搜索对比等模块同时配有说明文档、补充案例数据压缩包和备份文件可直接在MATLAB 2014/2019a/2021a中运行。代码采用参数化编程关键参数便于修改注释清晰既能帮助理解交替优化思路也能通过调整参数观察不同配置下的保密率差异加深对物理层安全与MIMO系统优化的认识。目前已有37人学习下载适合需要将理论推导与代码实现相结合的高年级本科生及研究生参考学习。1. 从IRS到坐标下降为什么保密率优化绕不开逐坐标调整IRS辅助的MIMO系统里反射单元动辄几十上百个每个单元都要调制一个单位模量的相移发射端还有多根天线要分配功率。如果直接把所有变量堆在一起上内点法非凸、变量多非常容易陷进某个糟糕的局部解。坐标下降算法是这类问题的天然选择——固定其他变量只更新当前一个变量把高维问题拆成多个一维优化。最大化IRS辅助MIMO系统保密率的坐标下降算法在MATLAB里用几十行脚本就能验证反射单元从8个加到64个保密率随信噪比上升的趋势一目了然。下面把信道建模、目标函数、迭代框架和MATLAB实现完整拆开照着一行行敲就能跑。2. 先把信道和保密率公式写进MATLAB从三跳信道到可分式目标2.1 IRS-MIMO的三段信道基站、反射面和两个终端完整的IRS辅助MIMO链路可以切分成三段。第一段是基站到IRS的MIMO信道维度是N乘M记作G_BIN是反射单元数M是基站天线数。第二段是IRS到合法用户和窃听者的两个信道分别记作h_RU和h_RE维度都是1乘N。第三段是基站直接到两个终端的直射信道分别记作h_DU和h_DE维度都是1乘M。直射信道一定要保留。IRS不会把来波全部吸收它只是叠加了一个受控的反射分量直射分量在实际场景里永远存在。很多刚接触这个方向的人为了简化仿真直接把直射信道设为零这样在高信噪比下得到的保密率会虚高尤其在有视距的场景里直射链路往往比反射链路更强。代码里我习惯这样生成信道% 信道参数与随机生成 M 4; % 基站天线数 N 16; % IRS反射单元数 pl_ru 2; % 用户IRS链路相对增益系数 pl_re 0.5; % 窃听IRS链路相对增益系数 % 复高斯随机信道每个系数 ~ CN(0,1) G_BI (randn(N,M) 1i*randn(N,M))/sqrt(2); h_RU sqrt(pl_ru) * (randn(1,N) 1i*randn(1,N))/sqrt(2); h_RE sqrt(pl_re) * (randn(1,N) 1i*randn(1,N))/sqrt(2); h_DU sqrt(0.8) * (randn(1,M) 1i*randn(1,M))/sqrt(2); h_DE sqrt(0.3) * (randn(1,M) 1i*randn(1,M))/sqrt(2);这里有一个容易被忽略的约定所有到终端的信道我都用行向量表示发射预编码w用列向量表示这样h_u * w直接是一个标量增益不用在代码里到处做共轭转置。涉及多天线的矩阵运算时这个约定能省掉一大堆.*维度匹配的报错。复高斯随机抽样后除sqrt(2)是为了让每个系数的平均功率归一化为1。pl_ru和pl_re是人为加上的链路增益系数相当于把用户路径损耗设成比窃听者低6dB左右。为什么要这样做如果用户链路和窃听链路的平均功率相同保密率最优解可能长期落在零附近坐标下降起点太低不容易观察到算法效果。仿真初期不妨让用户链路略占优势先把算法跑通再去调整成信道对称的极端场景。2.2 保密率目标函数为什么不能只把信噪比拉满保密率的定义是合法链路可达速率减去窃听链路可达速率外部取正。写成MATLAB能直接计算的标量形式R_s log2(1 p |h_u * w|^2 / sigma_u2) - log2(1 p |h_e * w|^2 / sigma_e2)其中等效信道h_u和h_e分别代表包含直射和IRS反射在内的复合信道h_u h_DU h_RU * diag(theta) * G_BI h_e h_DE h_RE * diag(theta) * G_BI这两个式子看完就能明白为什么坐标下降可行theta的第n个元素只出现在diag(theta)的第n个对角项上而diag(theta)*G_BI的每一行又只和h_RU、h_RE的第n个元素相乘。换句话说固定其他相移后当前第n个相移对保密率的影响是独立的一块可以单独做一维优化。注意保密率公式里我用了原始的差值而不是先取正。优化阶段如果直接写成max(0, R_s)当R_s小于零时梯度为零坐标下降会停在一个毫无意义的平台再也走不出来。实际代码里优化目标保留差值显示性能时再取正这一点非常重要。为什么不直接最大化合法用户信噪比因为那样做很可能让窃听者也获得更高信噪比差值反而变小。保密率的目标是拉开两条链路的差距不是单点最优。这也决定了功率分配不是简单的注水发射预编码也不是单纯的匹配滤波器。2.3 向量化计算保密率MATLAB里别写愚蠢的累加循环写保密率函数时新手容易这样先把h_u置零然后写一个for n 1:N的循环把每个反射单元的贡献累加进去。这样写没有错但是代码又臭又长而且后面做坐标下降时你想把某个相位拿出来修改变量这个累加逻辑会逼着你重复写很多遍等效信道构造。我一般把保密率计算封装成一个函数参数就是theta、w、p和各种信道矩阵function R calc_secrecy(theta, w, p, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2) h_u h_DU h_RU * diag(theta) * G_BI; h_e h_DE h_RE * diag(theta) * G_BI; S_u abs(h_u * w)^2 * p; S_e abs(h_e * w)^2 * p; R log2(1 S_u/sigma_u2) - log2(1 S_e/sigma_e2); end这段代码最大的好处是矩阵维度清晰。h_RU是1乘Ndiag(theta)是N乘NG_BI是N乘M乘出来正好是1乘M再和h_DU相加得到合法用户等效信道。如果不想用diag构造稀疏大矩阵也可以用h_RU .* theta.先对行向量按元素乘再去乘G_BI。效果一样N大时省一点内存但可读性不如diag。在线跑小例子用diag没问题真正做256个反射单元的蒙特卡洛时再考虑优化。把这个函数放到单独m文件里后面的坐标下降函数和主脚本都能反复调用。有一次我在别人代码里看到保密率计算直接嵌在主循环里坐标下降每试探一个候选相位就要重算一次h_u和h_e代码复制了三大段改参数时漏改一处就出bug。封装之后至少后面几个坑能少踩一半。3. 坐标下降怎么落地交替更新反射相位与发射预编码3.1 交替迭代框架外层换w内层扫theta坐标下降不是只能处理单个变量而是可以把变量分成块块之间交替优化。我们这里把发射预编码w、功率p和IRS相位theta当成两大块。外层先固定theta更新w和p内层固定w和p对theta逐坐标更新。这样整体还是坐标下降的骨架只不过内层又加了一层逐坐标扫描。给定theta后等效信道h_u被完全确定。这时候发射方向我一般选MRT也就是把w设成h_u的共轭方向w conj(h_u) / norm(h_u)为什么不是直接做广义特征分解因为保密率目标不是简单的信噪比最大化广义特征分解只能让两个增益比值最大和带两个log的保密率并不严格等价。MRT闭式、稳定、跑得快作为坐标下降的第一版最合适等算法整体收敛后再换更好的方向。这属于工程上的渐进式改进不算偷懒。功率p不要和w一起线性缩放否则w归一化就没有意义。固定w方向后p在0到P之间一维搜索。MATLAB自带的fminbnd够用function [w, p] update_wp(theta, P, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2) h_u h_DU h_RU * diag(theta) * G_BI; h_e h_DE h_RE * diag(theta) * G_BI; w h_u / norm(h_u); obj (pval) -calc_secrecy(theta, w, pval, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); [p, ~] fminbnd(obj, 1e-6, P); endfminbnd的搜索区间下限我给了1e-6而不是0。理论上p0也合法但p0时保密率恒为零任何小的试探都会让搜索器认为找到了边界返回一个没有意义的功率。加一个极小值下限后搜索至少在一个非退化的区间里进行。上限P就是系统总功率约束。3.2 固定(w,p)后逐坐标更新反射相移把0到2π打散成一维搜索这是整个坐标下降算法里最有意思的部分。固定w和p之后保密率R_s只与theta有关。固定其他theta_j不变改变theta_n时h_u和h_e都只有一行中的一项在动所以保密率关于theta_n是一个一维函数。理论上可以对这个一维函数做解析优化把theta_n写成atan2的形式。我试过只考虑用户链路的做法确实能化简出一个相位差公式但把窃听链路也放进去后两个log分子分母都有theta_natan2不再成立强行求导也得不到闭式解。保号后的一维函数还可能出现多峰梯度法容易陷在局部峰里。所以工程上我采用网格搜索把0到2π均匀离散成L个点逐一计算保密率取最大值对应的相位作为当前坐标的更新。这其实就是坐标下降里最常见的一种精确线搜索实现只是把连续区间离散化了。function theta coord_desc_theta(theta, w, p, L, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2) N size(G_BI, 1); phi_grid (0:L-1) * (2*pi/L); for n 1:N best_R -inf; best_phi angle(theta(n)); for idx 1:L theta_probe theta; theta_probe(n) exp(1i * phi_grid(idx)); R_probe calc_secrecy(theta_probe, w, p, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); if R_probe best_R best_R R_probe; best_phi phi_grid(idx); end end theta(n) exp(1i * best_phi); end end这段代码里有两个关键点。第一每次试探都用theta_probe theta而theta变量在循环中已经被前面坐标更新过所以这是严格的顺序坐标下降不是并行的同时更新。第二候选相位来自预先定义的phi_grid搜索到的最优相位立即写入theta(n)下一轮第n1个坐标看到的已经是新值。复杂度上N个反射单元每个单元L次保密率计算总共N乘L次。N16、L120时一次内层扫描大约1920次保密率计算。每次保密率里又包含两次diag矩阵乘法在普通笔记本上大概零点几秒。如果你用在线版MATLAB或者把N放大到64以上就可以考虑把diag(theta)*G_BI替换成bsxfun(times, G_BI, theta)或者theta .* G_BI这样的按行广播能省掉对角矩阵的乘法开销。3.3 停止准则与迭代上限保密率增量、相位漂移和最大循环外层交替迭代至少要设三个停止条件。第一是保密率增量也就是本轮和上一轮的R_s差值小于某个容忍度比如1e-3 bit/s/Hz。第二是相位向量的变化量norm(theta_new - theta_old)小于1e-2。第三是最大迭代上限maxIter防止算法在两个候选相位之间反复横跳。关于相位漂移有一个容易踩的坑theta是复数角度0和2π在数值上代表同一个相位但直接做差会得到2π的虚假漂移。好在我的坐标更新把候选角都限制在[0, 2π)区间内更新后angle函数也落在同一区间所以实际没遇到这个问题。如果你自己写连续角度优化就要把相位差值映射回(-π, π]再比较。外层框架我一般这么写theta exp(1i * 2*pi*rand(N,1)); w zeros(M,1); p 0; R_hist zeros(maxIter,1); for iter 1:maxIter [w, p] update_wp(theta, P, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); for inner 1:5 theta coord_desc_theta(theta, w, p, L, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); end R_hist(iter) calc_secrecy(theta, w, p, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); if iter 1 abs(R_hist(iter) - R_hist(iter-1)) tol break; end end内层固定w和p的时候多跑几轮theta扫描能保证theta在本次外层迭代中尽量收敛。否则只扫描一轮theta还没稳定就去更新w整体收敛曲线会出现毛刺。我试过内层只跑1轮最后的保密率序列会有明显波动但最终结果相差不大。想要文章里那种漂亮平滑的收敛曲线内层5轮就够。这套算法每轮目标函数都不会下降因为w、p和theta的每一步更新都是朝着至少不坏的方向走。所以理论上R_hist应该是单调不减的。如果看到下降多半是记录时机错了下一章会专门讲这个问题。4. 复现一份完整MATLAB脚本参数表、主循环和曲线输出4.1 参数表天线数、反射单元数、信噪比区域怎么选参数不是随便拍的。反射单元太少比如N4曲线太平坦看不出IRS优势N太多比如N128坐标下降每个内层要扫1万多次保密率教学实验会很慢。我建议第一次复现时用N16迭代30轮整个脚本运行时间控制在十秒以内。以下是一份自测过的参数表可以直接照抄参数推荐值作用说明M4基站天线数越大阵列增益越高N16IRS反射单元数越多相位自由度越大P_dB0 到 20发射功率范围dB刻度sigma_u21合法用户噪声功率sigma_e21窃听者噪声功率L120每个坐标一维搜索的离散点数maxIter30外层交替迭代上限tol1e-3保密率增量收敛门限pl_ru / pl_re2 / 0.5反射链路相对增益调节用户优势噪声功率都取1发射功率P_dB从0开始这意味着初始信噪比从0dB附近起步保密率曲线可以看得很清楚。如果你想让合法用户链路优势更大可以把pl_ru提高到3pl_re压到0.3。但注意不要一个参数试到底不同信道增益下坐标下降的收敛速度和最终保密率都差不少。4.2 主脚本框架先随机信道再交替求解把第二章和第三章的函数组合起来就是一个完整的主脚本。这个脚本可以直接新建脚本文件按F5运行clear; clc; rng(2025); % 固定随机种子保证可复现 M 4; N 16; P_dB 10; P 10^(P_dB/10); sigma_u2 1; sigma_e2 1; maxIter 30; tol 1e-3; L 120; pl_ru 2; pl_re 0.5; G_BI (randn(N,M) 1i*randn(N,M))/sqrt(2); h_RU sqrt(pl_ru) * (randn(1,N) 1i*randn(1,N))/sqrt(2); h_RE sqrt(pl_re) * (randn(1,N) 1i*randn(1,N))/sqrt(2); h_DU sqrt(0.8) * (randn(1,M) 1i*randn(1,M))/sqrt(2); h_DE sqrt(0.3) * (randn(1,M) 1i*randn(1,M))/sqrt(2); theta exp(1i * 2*pi*rand(N,1)); w zeros(M,1); p 0; R_hist zeros(maxIter,1); for iter 1:maxIter [w, p] update_wp(theta, P, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); for inner 1:5 theta coord_desc_theta(theta, w, p, L, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); end R_hist(iter) calc_secrecy(theta, w, p, G_BI, h_RU, h_RE, h_DU, h_DE, sigma_u2, sigma_e2); if iter 1 abs(R_hist(iter) - R_hist(iter-1)) tol fprintf(收敛于第%d轮\n, iter); break; end end figure; plot(R_hist(1:iter), -o); xlabel(迭代次数); ylabel(保密率 (bit/s/Hz)); title(坐标下降收敛曲线); grid on;运行完你会看到一条单调上升的曲线前几轮涨得很快后面逐渐趋于平缓。这就是坐标下降在两个变量块之间交替前进的特征。如果曲线出现了突然下降优先检查R_hist记录的位置是不是在更新w之后、更新theta之前——那个位置R值一定会变乱。4.3 后处理与曲线复现保密率随噪声功率和反射单元数量变化手工复现一张论文级曲线时需要扫两个维度一是发射功率P_dB二是反射单元数量N。把P_dB从0取到20步长2dBN取8、16、32三条线。每条线都先跑上面的主脚本记录收敛后的R_s最后用loglog或semilogx画出来。注意每次改变N都要重新生成信道rng种子也要重新设置。为了对比公平可以在外层加一个蒙特卡洛循环num_realizations 100; R_avg zeros(length(P_dB_set), length(N_set)); for i 1:length(N_set) N N_set(i); for j 1:length(P_dB_set) P 10^(P_dB_set(j)/10); sum_R 0; for r 1:num_realizations rng(r * 100 i * 10 j); % 重新生成信道和初始化theta ... sum_R sum_R final_R; end R_avg(i, j) sum_R / num_realizations; end end这里随机种子的编号规则不是唯一的但目的是让每个(N, P)组合在每次蒙特卡洛中都有一组独立信道同时不同组合之间可以重复对照。不要偷懒用同一个rng(0)跑完全部否则所有曲线都来自同一组信道N增加纯粹是相位自由度变多信道随机性完全丢失曲线会异常平滑不代表真实性能。5. 避坑指南这5个问题最容易让保密率迭代翻车5.1 MATLAB中文注释乱码脚本一开就废现象用MATLAB 2023或2022版本打开别人分享的脚本中文注释全部变成方框或乱码代码本身还能运行但完全没法排查参数。原因MATLAB从R2020b之后在Windows上默认文件编码从GBK切换到UTF-8而早期脚本保存为GBK反过来用新版本保存UTF-8脚本到老版本环境也可能乱码。这个问题在中文技术群里几乎每周都有人问。解决统一在MATLAB预设项里找到“常规—MAT文件—字符编码设置”把默认编码改成UTF-8。如果你只负责运行脚本直接删除注释里的中文全部改成英文注释最省事。我在分享可复现代码时函数文档都写英文只有在文章正文里详细解释中文含义避免编码问题干扰复用。5.2 保密率曲线明明单调递增某轮突然下跌现象主脚本跑完R_hist曲线在第5轮以后出现一个明显的向下尖峰然后又继续上升。原因记录时机不对。先执行update_wp(theta,...)再计算R_hist(iter)等于用新w和新p计算保密率而上一轮的R_hist用的是旧w和旧p两个不在同一组变量上自然可能下降。解决把保密率记录放到内层theta收敛之后且下一次迭代更新w之前。也就是先记录再更新w。或者干脆把记录放在每个外层迭代末尾代码顺序改成先计算R_hist再更新w。这样每一步坐标下降都至少不会变差曲线单调性就回来了。5.3 一维功率搜索返回极端值保密率直接为负现象fminbnd搜索p时返回结果要么是1e-6要么是P上限整条曲线保密率接近零。原因p在低信噪比区域目标函数对p不敏感fminbnd的黄金分割搜索把p推到边界。下限为0时更加严重p0导致目标函数恒零搜索器判定边界就是最优。解决下限改1e-6上限保持P。如果发现搜索结果总是贴着上下限说明w方向给得不合适先检查MRT方向是不是对着h_u而不是h_e。建议在update_wp函数里临时打印h_u和h_e的范数肉眼确认合法通道更强。很多所谓保密率翻车其实是信道增益配置反了。5.4 网格搜索相位太粗最优解差了0.5dB现象L20时保密率曲线比文献参考值低0.3到0.5dB调大L又觉得慢。原因L20对应18度一个步长最优相位落在两个网格点之间时误差最大可达9度这对IRS波束成形是明显损失尤其是N较大时各个反射单元的小误差会累加。解决第一轮用L120粗筛锁定最优相位所在区间第二轮在最优相位附近加减一个步长用更细的网格或直接调用fminbnd在这个小区间内细搜。下面是细搜的替代写法phi_best best_phi; phi_range [phi_best - 2*pi/L, phi_best 2*pi/L]; phi_refined fminbnd((phi_probe) -calc_secrecy_with_angle(phi_probe), phi_range(1), phi_range(2));注意calc_secrecy_with_angle要把临时相位写回theta的第n个位置再调用保密率函数。这个是纯网格搜索的外部包装仍然符合坐标下降的框架只是每个坐标的线搜索更精确了一点。5.5 不同N的保密率曲线交叉得乱七八糟现象画N8、16、32三条曲线时N8的个别点反而比N32还高曲线交叉图表没法用。原因每次循环跑完一组信道后没有重置随机种子或者用同一个种子导致N8和N32时信道样本不同比较的是不同信道下的点自然没有单调性。解决把所有待比较场景的随机种子固定在同一组编号规则下比如rng(100*N P_dB)这种方式确保每次重新生成信道都得到相同随机序列。更严谨的做法是蒙特卡洛求平均每个N和P组合都各自跑100次独立信道取平均再用平均保密率画曲线。这样单次随机性被抹平曲线才光滑可复现。6. 让坐标下降再快一截随机坐标、冷启动与离散化技巧到手可用之后可以试着把算法改得更快。第一个技巧是随机坐标选择。当N很大顺序扫描所有反射单元耗时太长每一轮随机选择一部分坐标来更新比如每次只更新三分之一其他坐标保持不动。这在统计意义上仍然能收敛而且对变量之间的强相关性有一定的“破坏”作用有时反而跳出局部解。不过随机选择会让保密率曲线带上噪声收敛判据建议改成对最近5轮取滑动平均而不是看单轮增量。第二个技巧是冷启动。把坐标下降里的theta初始相位从0到2π随机取这没错但随机种子对最终结果影响很大。我一般跑3到5个随机初始theta取收敛后保密率最高的那一个作为最终答案。运气好的时候一个糟糕的初始相位会让坐标下降卡在局部解多试几次成本很低收益却很实在。第三个技巧是支持离散相位约束。真实IRS相移通常只有1位或2位精度比如只能取{0, π}或者{0, π/2, π, 3π/2}。你不需要改坐标下降框架只要把网格搜索的候选点从连续均匀网格换成离散相位集合代码几乎不用动。离散化之后每个坐标的搜索范围变小运行速度更快但保密率会有损失典型值在0.5到1.5dB之间。如果你的目标是工程落地这一步必须做。验证算法是否真的work我有个习惯先画一条随机相位下的保密率作为baseline再画坐标下降收敛后的点对比拉升了多少。如果一个系统中随机相位和优化后相位几乎没有差别说明问题不在算法而在信道模型——可能窃听者信道太强或者相位自由度根本没有生效。这时候不要急着调算法回头检查信道增益系数。我在跑这类仿真时习惯先把N设成8跑通再逐步加大否则收敛慢的因素和信道因素混在一起定位问题很痛苦。希望帮到你。本文还有配套的精品资源点击获取