简介面向无线通信与优化算法方向的MATLAB学习者压缩包提供了IRS辅助MIMO系统保密率最大化的坐标下降算法完整实现。代码采用参数化编程变量与参数便于调整适合电子信息工程、数学等专业学生用于课程设计、期末大作业或毕业设计也可作为智能优化与物理层安全研究的入门参考。资源共7个文件其中6个m文件覆盖主程序、坐标下降迭代、全搜索对比、IRS信道建模、固定相位下的MIMO容量计算等功能另有1个md文件作为说明文档整个压缩包仅12KB结构精简可直接在MATLAB 2014/2019a/2021a中运行。包内包含proposed_algorithm.m、full_search.m等脚本读者可对照主程序理解迭代逻辑并通过全搜索法验证坐标下降的逼近性能同时借助IRS信道与容量函数深入掌握系统模型。资源已有160人学习适合需要快速上手IRS辅助MIMO系统仿真、希望获取可修改代码框架的研究者与高年级本科生。1. 用坐标下降啃下 IRS 辅助 MIMO 保密率优化一位同事在仿真里发现32 根天线的发射端面对隔壁窃听者时单纯增加发射功率反而让保密率停滞。IRS 的出现改变了玩法用无源反射单元重塑合法用户和窃听者之间的信道差异代价是引入一块单位模约束的非凸相移使问题变成 NP-hard。坐标下降算法是最实用的切入方式把发射预编码与 IRS 相移拆成两个交替求解块相移块内再逐坐标更新。本文给出在 MATLAB 中实现整套算法的工程思路覆盖信道生成、波束成形闭式解、相移逐坐标扫描、收敛判断和验证技巧适合做物理层安全的在校生和无线算法工程师。2. 保密率目标为何非凸坐标下降的拆分思路2.1 IRS 辅助 MIMO 系统模型与保密率表达式考虑一个单小区下行链路基站配置 Nt 根天线IRS 含 N 个反射单元合法用户配置 Nu 根天线窃听者配置 Ne 根天线。为了把坐标下降算法的核心逻辑讲清楚这里先采用单流传输发射符号为 s预编码向量为 w发射功率约束为 ||w||^2 ≤ P。接收端使用最大比合并这样用户端信干噪比只取决于等效信道矩阵的 Gram 矩阵多天线带来的合并增益已经包含在二次型里。令 G 表示基站到 IRS 的 N × Nt 信道矩阵H_ru 表示 IRS 到用户的 Nu × N 信道矩阵H_bu 表示基站到用户的直达信道矩阵。IRS 的反射系数写成对角阵 Φ diag(θ_1, …, θ_N)其中 |θ_n| 1。用户侧有效信道矩阵为H_u(θ) H_bu H_ru · Φ · G窃听者侧同理为 H_e(θ) H_be H_re · Φ · G。单流场景下保密率可写为R_s(w, θ) log2(1 w^H H_u^H H_u w / σ_u^2) − log2(1 w^H H_e^H H_e w / σ_e^2)再取与 0 的较大值。严格说当用户和窃听者都支持多流时保密率应写成 log det 形式但逐坐标更新的整体框架不变最后一章会给出扩展方向。下面用一个参数对照表固定本文后面代码使用的符号避免在函数调用里来回翻找。符号含义本文演示值Nt基站天线数4NIRS 反射单元数32Nu / Ne用户 / 窃听者天线数2 / 2KRician K 因子dB10P发射功率预算1线性σ_u^2 / σ_e^2用户 / 窃听者噪声功率由 SNR 反推L每坐标相位扫描点数642.2 分块坐标下降发射波束成形与 IRS 相移交替更新直接同时优化 w 和 θ 的困难在于θ 的单位模约束不是凸集而 w 出现在信道矩阵内部导致目标函数对 θ 并没有好的凸结构。坐标下降的做法是把变量分块每一轮先固定 θ 更新 w再固定 w 更新 θ。这样每个子问题都比原问题简单且每次更新都保证保密率不降。固定 θ 时等效信道 H_u、H_e 是常数矩阵。对于单流保密率常见做法是求如下广义特征值问题的最大特征向量max_w w^H (H_u^H H_u εI) w / w^H (H_e^H H_e εI) w这里的 ε 是正则项避免 Gram 矩阵奇异。这个解不是严格的最优闭式解但实现简单实践中收敛稳定。固定 w 时θ 的更新则完全交给坐标扫描保持其他 θ_n 不变只让第 n 个反射单元在一组离散相位上取值选出使 R_s 最大的相位。整个迭代过程用伪代码可以写成下面这样这也是后面 MATLAB 主脚本的骨架。% 初始化 theta exp(1j * 2 * pi * rand(N, 1)); w randn(Nt, 1); w w / norm(w); for iter 1:max_iter % 块1固定 theta 更新 w H_u H_bu (H_ru .* theta.) * G; H_e H_be (H_re .* theta.) * G; w update_beamformer(H_u, H_e, P); % 块2固定 w对每个反射单元执行坐标下降 for n 1:N theta coordinate_descent(theta, n, w, H_ru, H_re, G, H_bu, H_be); end R_history(iter) compute_rate(w, theta); end这里theta.是数学转置而非共轭转置目的是让 H_ru 的第 n 列乘以 theta(n) 后再和 G 相乘。坐标下降法的单调性来自每次更新都会让目标值上升因此主循环里的R_history天然是单调不减序列可以直接用来画收敛曲线。为什么不直接对 θ 做梯度上升再投影到单位圆上原因是保密率函数在部分区域内呈多峰形态梯度投影法极易在反射单元数量较大时撞到平坦区。逐坐标扫描虽然每次都多花一点时间但每个维度都能得到真实的目标函数值对调试非常友好。特别是在仿真初始阶段你至少能分辨出性能差是因为算法没收敛还是因为信道模型本身设置错了。3. MATLAB 代码结构信道生成、波束成形与逐坐标相移更新3.1 信道生成与参数表信道生成是 IRS 仿真里最容易出错的地方。IRS 反射阵面一般部署在基站与用户之间的 LoS 路径上所以基站到 IRS 链路、IRS 到用户链路都应带有 Rician 衰落。下面的函数生成五组信道矩阵并直接通过pathloss_dB控制平均接收功率。function [G, H_ru, H_re, H_bu, H_be] gen_irs_channels(...) % 生成 IRS 辅助 MIMO 链路信道 % 输入: Nt, Nu, Ne, N, K_dB, pathloss_dB % 输出: % G : N x Nt 基站到 IRS % H_ru : Nu x N IRS 到用户 % H_re : Ne x N IRS 到窃听者 % H_bu : Nu x Nt 基站到用户直达 % H_be : Ne x Nt 基站到窃听者直达 K 10^(K_dB / 10); pl 10^(-pathloss_dB / 20); % 基站-IRSRicianLoS 部分用全一矩阵近似平面波前 los_G ones(N, Nt) / sqrt(Nt); nlos_G (randn(N, Nt) 1j*randn(N, Nt)) / sqrt(2); G pl * (sqrt(K/(K1))*los_G sqrt(1/(K1))*nlos_G); % IRS-用户 los_ru ones(Nu, N) / sqrt(N); nlos_ru (randn(Nu, N) 1j*randn(Nu, N)) / sqrt(2); H_ru pl * (sqrt(K/(K1))*los_ru sqrt(1/(K1))*nlos_ru); % IRS-窃听者窃听者位置不同LoS 方向角信息丢失直接采用瑞利更稳妥 nlos_re (randn(Ne, N) 1j*randn(Ne, N)) / sqrt(2); H_re pl * nlos_re; % 直达链路 H_bu pl * (randn(Nu, Nt) 1j*randn(Nu, Nt)) / sqrt(2); H_be pl * (randn(Ne, Nt) 1j*randn(Ne, Nt)) / sqrt(2); end代码里的sqrt(K/(K1))和sqrt(1/(K1))是 Rician 分解的标准写法。窃听者链路故意不给 LoS因为实际中窃听者不会把自己暴露在强视距路径上这个设置会让保密率曲线更接近真实场景。pathloss_dB在演示里取 0 表示所有信道平均增益相同实际应用时应该根据基站到 IRS、IRS 到用户、直达链路的距离分别计算。3.2 固定 IRS 相位时发射波束成形的闭式求解固定 θ 之后H_u 和 H_e 都变成了常数矩阵。下面这个函数用 MATLAB 的eig求解广义特征值问题取最大广义特征值对应的右特征向量作为波束成形向量。function w update_beamformer(H_u, H_e, P) % 固定 IRS 相位时的发射波束成形 % 最大化 (1 gamma_u) / (1 gamma_e) 的广义特征值近似解 Nt size(H_u, 2); reg 1e-6 * eye(Nt); A H_u * H_u reg; B H_e * H_e reg; [V, D] eig(A, B); % 广义特征值分解 [~, idx] max(real(diag(D))); % 广义特征值可能不排序必须手动取最大 w V(:, idx); w sqrt(P) * w / norm(w); % 满足功率约束 end这里的正则项reg很关键。当窃听者天线数小于基站天线数时B 可能是奇异矩阵直接调用eig(A, B)会得到 Inf 或 NaN。加入1e-6 * eye(Nt)后数值上稳定且对广义特征向量方向的影响可忽略。取最大广义特征值的原理是目标函数近似等于 w^H A w / w^H B w 的对数形式而广义特征向量正是使该比值最大的方向。需要提醒的是这个解是对原保密率函数的一种近似。严格的最大化 log(1γ_u) − log(1γ_e) 并不会有这么简单的闭式解。工程脚本里先用广义特征值把上半轮跑通再在调优阶段换成 MM 算法或 CVX 求解是常见的演进路径。3.3 IRS 相位的逐坐标扫描更新代码IRS 相位更新是本标题的核心。假设当前迭代已经有了一组 θ 和 w接下来要对第 n 个反射单元进行坐标更新。为了降低实现难度采用离散相位扫描在 [0, 2π) 内取 L 个候选相位对每个候选相位重新计算有效信道矩阵和保密率选出使保密率最大的相位。function theta coordinate_descent(theta, n, w, H_ru, H_re, G, H_bu, H_be, L) % 固定其他反射单元只更新第 n 个反射单元 % 返回更新后的 theta 向量 ru_n H_ru(:, n); % Nu x 1 re_n H_re(:, n); % Ne x 1 g_n G(n, :); % 1 x Nt % 先把第 n 项从有效信道中剔除 theta_off theta; theta_off(n) 0; H_u_off H_bu (H_ru .* theta_off.) * G; H_e_off H_be (H_re .* theta_off.) * G; % 候选相位 ph exp(1j * 2 * pi * (0:L-1) / L); best_R -inf; best_phase theta(n); for k 1:L H_u H_u_off (ru_n * ph(k)) * g_n; H_e H_e_off (re_n * ph(k)) * g_n; gamma_u real(w * (H_u * H_u) * w); gamma_e real(w * (H_e * H_e) * w); R log2(1 gamma_u) - log2(1 gamma_e); if R best_R best_R R; best_phase ph(k); end end theta(n) best_phase; end这段代码里H_u_off和H_e_off相当于挖去了第 n 个反射单元贡献后的信道矩阵。候选相位ph(k)是通过(ru_n * ph(k)) * g_n加回去的乘号顺序对应H_ru(:,n) * θ_n * G(n,:)的秩一结构。每个候选相位只需要两次矩阵乘法整个坐标下降一轮的复杂度大约是 O(N · L · max(Nu,Ne) · Nt)对 N64、L64 的小仿真完全可接受。如果想把扫描改成连续优化可以用 MATLAB 优化工具箱里的fminbnd对相位角直接搜索。但需要注意保密率在相位角上可能有多个峰值fminbnd只能收敛到局部极值。坐标扫描的方式天然支持多峰问题而且可以直接模拟低分辨率 IRS 单元比如每单元只有 2 或 3 bit 相移这和硬件实现更贴近。4. 仿真中的收敛判定与三个常见坑4.1 主循环、收敛曲线与 MATLAB 优化工具箱对比主脚本的职责是串起信道生成、波束成形、坐标下降和曲线绘制。下面给出核心循环的写法% 参数设置 Nt 4; Nu 2; Ne 2; N 32; P 1; sigma2 10^(-snr_dB/10); max_iter 40; L 64; % 生成信道 [G, H_ru, H_re, H_bu, H_be] gen_irs_channels(...); % 初始化 theta exp(1j * 2 * pi * rand(N, 1)); w randn(Nt, 1); w w / norm(w); R_history zeros(max_iter, 1); for iter 1:max_iter H_u H_bu (H_ru .* theta.) * G; H_e H_be (H_re .* theta.) * G; w update_beamformer(H_u, H_e, P); for n 1:N theta coordinate_descent(theta, n, w, ... H_ru, H_re, G, H_bu, H_be, L); end H_u H_bu (H_ru .* theta.) * G; H_e H_be (H_re .* theta.) * G; gamma_u real(w * (H_u * H_u) * w) / sigma2; gamma_e real(w * (H_e * H_e) * w) / sigma2; R_history(iter) max(0, log2(1 gamma_u) - log2(1 gamma_e)); end收敛曲线通常会在前 5 到 15 轮快速上升后面进入慢速微调阶段。判断收敛的标准不是保密率不再变化而是连续若干轮变化量小于某个阈值比如abs(R_history(iter) - R_history(iter-1)) 1e-4。如果曲线出现抖动下降几乎可以肯定是坐标扫描代码里把某个索引写错了因为坐标下降的单调性应该保证曲线不上扬。如果想偷懒用 MATLAB 优化工具箱里的fmincon或ga做对比建议只跑 N8 的小规模场景否则在 N64 时求解质量不如坐标扫描稳定。坐标下降的价值在于它在反射单元数量变大后仍然保持接近线性的复杂度而全局优化方法会在 N 超过 20 之后明显变慢。4.2 单位模约束、扫描点数和初始化陷阱第一个常见坑是在更新相位后忘记归一化。theta必须始终保持abs(theta) 1但浮点运算里ph(k)本身是单位模所以只要不直接对 theta 做加法就不会破坏约束。如果写的是theta(n) theta(n) delta再投影到单位圆性能会明显下降。第二个坑是扫描点数 L 的选择。L 太小时相位量化误差会直接抬高误差地板。L16 时每个相位间隔 22.5°保密率可能损失 0.3 bit/s/HzL64 时间隔约 5.6°精度已经接近连续相位下界。如果仿真目标是论文图建议 L 取 128并把坐标下降轮数降到 20总时间仍在可接受范围。第三个坑是初始化。随机相位初始化对坐标下降来说足够但为了获得更好的局部最优最好跑 3 到 5 个随机种子取保密率最高的一组作为最终结果。每次随机种子重置后信道也应重新生成否则会高估 IRS 的增益。另一个稍反直觉的结论是从所有相位都为 0 初始化往往会收敛到较差的局部极值因为 IRS 的镜面反射方向被固定坐标更新后很难越过相位势垒。随机初始化反而能覆盖更多方向。5. 验证代码正确性单步扰动检验与 MIMO 信道容量图像坐标下降代码写完后不要直接上大参数仿真先用一个“单步扰动检验”确认逐坐标更新的逻辑没写反。具体做法是在任意一次迭代后固定当前 w对第 n 个反射单元手动遍历全部 L 个候选相位画出保密率与 θ_n 角度的关系曲线并标出算法最终选择的点。如果选出的点确实落在曲线的最高峰上说明该维度的更新方向正确如果落在波谷多半是H_u_off构造时把第 n 项的错误符号留下了常见原因是没有将theta_off(n)置零。还可以进一步做穷举对比当 N 很小时比如 N4L8遍历全部 8^44096 种相位组合计算每种组合下的保密率上界。坐标下降得到的结果与穷举结果之间的差距如果小于 0.05 bit/s/Hz说明整个迭代框架没有结构性问题。这是我最常用的一招也是排查信道矩阵维度错误最有效的方法。验证通过后可以绘制固定某个反射单元相位变化时的保密率曲线俗称“相位扫描图”。代码片段如下phase_range linspace(0, 2*pi, 128); R_array zeros(size(phase_range)); for k 1:length(phase_range) theta_test theta; % 使用算法收敛后的相位 theta_test(1) exp(1j * phase_range(k)); H_u H_bu (H_ru .* theta_test.) * G; H_e H_be (H_re .* theta_test.) * G; % 固定 w 计算保密率 R_array(k) compute_rate_at_w(w, H_u, H_e, sigma2); end plot(phase_range/pi, R_array, LineWidth, 1.2); xlabel(Phase / \pi); ylabel(Secrecy rate (bit/s/Hz));这张图本质上就是你说的“MIMO 信道容量图像”的一种变体它不是原始信道容量而是保密率对相移的响应面。通过观察这条曲线的峰值位置可以反推该反射单元在级联信道中贡献的相位偏置也能解释为什么坐标下降总是先更新具有最大信道增益的反射单元。实际工程里我还会把每个单元的“相位-保密率”曲线叠加在同一张图上找出贡献最小的几个单元并在后续硬件控制中跳过它们进而降低 IRS 配置开销。这种优化技巧在反射单元数量超过 256 时尤其有价值。本文还有配套的精品资源点击获取 
长河君 · 2026/9/12 7:04:42 · 阅读 1.2k IRS辅助MIMO保密率优化：坐标下降算法与MATLAB实现
简介面向无线通信与优化算法方向的MATLAB学习者压缩包提供了IRS辅助MIMO系统保密率最大化的坐标下降算法完整实现。代码采用参数化编程变量与参数便于调整适合电子信息工程、数学等专业学生用于课程设计、期末大作业或毕业设计也可作为智能优化与物理层安全研究的入门参考。资源共7个文件其中6个m文件覆盖主程序、坐标下降迭代、全搜索对比、IRS信道建模、固定相位下的MIMO容量计算等功能另有1个md文件作为说明文档整个压缩包仅12KB结构精简可直接在MATLAB 2014/2019a/2021a中运行。包内包含proposed_algorithm.m、full_search.m等脚本读者可对照主程序理解迭代逻辑并通过全搜索法验证坐标下降的逼近性能同时借助IRS信道与容量函数深入掌握系统模型。资源已有160人学习适合需要快速上手IRS辅助MIMO系统仿真、希望获取可修改代码框架的研究者与高年级本科生。1. 用坐标下降啃下 IRS 辅助 MIMO 保密率优化一位同事在仿真里发现32 根天线的发射端面对隔壁窃听者时单纯增加发射功率反而让保密率停滞。IRS 的出现改变了玩法用无源反射单元重塑合法用户和窃听者之间的信道差异代价是引入一块单位模约束的非凸相移使问题变成 NP-hard。坐标下降算法是最实用的切入方式把发射预编码与 IRS 相移拆成两个交替求解块相移块内再逐坐标更新。本文给出在 MATLAB 中实现整套算法的工程思路覆盖信道生成、波束成形闭式解、相移逐坐标扫描、收敛判断和验证技巧适合做物理层安全的在校生和无线算法工程师。2. 保密率目标为何非凸坐标下降的拆分思路2.1 IRS 辅助 MIMO 系统模型与保密率表达式考虑一个单小区下行链路基站配置 Nt 根天线IRS 含 N 个反射单元合法用户配置 Nu 根天线窃听者配置 Ne 根天线。为了把坐标下降算法的核心逻辑讲清楚这里先采用单流传输发射符号为 s预编码向量为 w发射功率约束为 ||w||^2 ≤ P。接收端使用最大比合并这样用户端信干噪比只取决于等效信道矩阵的 Gram 矩阵多天线带来的合并增益已经包含在二次型里。令 G 表示基站到 IRS 的 N × Nt 信道矩阵H_ru 表示 IRS 到用户的 Nu × N 信道矩阵H_bu 表示基站到用户的直达信道矩阵。IRS 的反射系数写成对角阵 Φ diag(θ_1, …, θ_N)其中 |θ_n| 1。用户侧有效信道矩阵为H_u(θ) H_bu H_ru · Φ · G窃听者侧同理为 H_e(θ) H_be H_re · Φ · G。单流场景下保密率可写为R_s(w, θ) log2(1 w^H H_u^H H_u w / σ_u^2) − log2(1 w^H H_e^H H_e w / σ_e^2)再取与 0 的较大值。严格说当用户和窃听者都支持多流时保密率应写成 log det 形式但逐坐标更新的整体框架不变最后一章会给出扩展方向。下面用一个参数对照表固定本文后面代码使用的符号避免在函数调用里来回翻找。符号含义本文演示值Nt基站天线数4NIRS 反射单元数32Nu / Ne用户 / 窃听者天线数2 / 2KRician K 因子dB10P发射功率预算1线性σ_u^2 / σ_e^2用户 / 窃听者噪声功率由 SNR 反推L每坐标相位扫描点数642.2 分块坐标下降发射波束成形与 IRS 相移交替更新直接同时优化 w 和 θ 的困难在于θ 的单位模约束不是凸集而 w 出现在信道矩阵内部导致目标函数对 θ 并没有好的凸结构。坐标下降的做法是把变量分块每一轮先固定 θ 更新 w再固定 w 更新 θ。这样每个子问题都比原问题简单且每次更新都保证保密率不降。固定 θ 时等效信道 H_u、H_e 是常数矩阵。对于单流保密率常见做法是求如下广义特征值问题的最大特征向量max_w w^H (H_u^H H_u εI) w / w^H (H_e^H H_e εI) w这里的 ε 是正则项避免 Gram 矩阵奇异。这个解不是严格的最优闭式解但实现简单实践中收敛稳定。固定 w 时θ 的更新则完全交给坐标扫描保持其他 θ_n 不变只让第 n 个反射单元在一组离散相位上取值选出使 R_s 最大的相位。整个迭代过程用伪代码可以写成下面这样这也是后面 MATLAB 主脚本的骨架。% 初始化 theta exp(1j * 2 * pi * rand(N, 1)); w randn(Nt, 1); w w / norm(w); for iter 1:max_iter % 块1固定 theta 更新 w H_u H_bu (H_ru .* theta.) * G; H_e H_be (H_re .* theta.) * G; w update_beamformer(H_u, H_e, P); % 块2固定 w对每个反射单元执行坐标下降 for n 1:N theta coordinate_descent(theta, n, w, H_ru, H_re, G, H_bu, H_be); end R_history(iter) compute_rate(w, theta); end这里theta.是数学转置而非共轭转置目的是让 H_ru 的第 n 列乘以 theta(n) 后再和 G 相乘。坐标下降法的单调性来自每次更新都会让目标值上升因此主循环里的R_history天然是单调不减序列可以直接用来画收敛曲线。为什么不直接对 θ 做梯度上升再投影到单位圆上原因是保密率函数在部分区域内呈多峰形态梯度投影法极易在反射单元数量较大时撞到平坦区。逐坐标扫描虽然每次都多花一点时间但每个维度都能得到真实的目标函数值对调试非常友好。特别是在仿真初始阶段你至少能分辨出性能差是因为算法没收敛还是因为信道模型本身设置错了。3. MATLAB 代码结构信道生成、波束成形与逐坐标相移更新3.1 信道生成与参数表信道生成是 IRS 仿真里最容易出错的地方。IRS 反射阵面一般部署在基站与用户之间的 LoS 路径上所以基站到 IRS 链路、IRS 到用户链路都应带有 Rician 衰落。下面的函数生成五组信道矩阵并直接通过pathloss_dB控制平均接收功率。function [G, H_ru, H_re, H_bu, H_be] gen_irs_channels(...) % 生成 IRS 辅助 MIMO 链路信道 % 输入: Nt, Nu, Ne, N, K_dB, pathloss_dB % 输出: % G : N x Nt 基站到 IRS % H_ru : Nu x N IRS 到用户 % H_re : Ne x N IRS 到窃听者 % H_bu : Nu x Nt 基站到用户直达 % H_be : Ne x Nt 基站到窃听者直达 K 10^(K_dB / 10); pl 10^(-pathloss_dB / 20); % 基站-IRSRicianLoS 部分用全一矩阵近似平面波前 los_G ones(N, Nt) / sqrt(Nt); nlos_G (randn(N, Nt) 1j*randn(N, Nt)) / sqrt(2); G pl * (sqrt(K/(K1))*los_G sqrt(1/(K1))*nlos_G); % IRS-用户 los_ru ones(Nu, N) / sqrt(N); nlos_ru (randn(Nu, N) 1j*randn(Nu, N)) / sqrt(2); H_ru pl * (sqrt(K/(K1))*los_ru sqrt(1/(K1))*nlos_ru); % IRS-窃听者窃听者位置不同LoS 方向角信息丢失直接采用瑞利更稳妥 nlos_re (randn(Ne, N) 1j*randn(Ne, N)) / sqrt(2); H_re pl * nlos_re; % 直达链路 H_bu pl * (randn(Nu, Nt) 1j*randn(Nu, Nt)) / sqrt(2); H_be pl * (randn(Ne, Nt) 1j*randn(Ne, Nt)) / sqrt(2); end代码里的sqrt(K/(K1))和sqrt(1/(K1))是 Rician 分解的标准写法。窃听者链路故意不给 LoS因为实际中窃听者不会把自己暴露在强视距路径上这个设置会让保密率曲线更接近真实场景。pathloss_dB在演示里取 0 表示所有信道平均增益相同实际应用时应该根据基站到 IRS、IRS 到用户、直达链路的距离分别计算。3.2 固定 IRS 相位时发射波束成形的闭式求解固定 θ 之后H_u 和 H_e 都变成了常数矩阵。下面这个函数用 MATLAB 的eig求解广义特征值问题取最大广义特征值对应的右特征向量作为波束成形向量。function w update_beamformer(H_u, H_e, P) % 固定 IRS 相位时的发射波束成形 % 最大化 (1 gamma_u) / (1 gamma_e) 的广义特征值近似解 Nt size(H_u, 2); reg 1e-6 * eye(Nt); A H_u * H_u reg; B H_e * H_e reg; [V, D] eig(A, B); % 广义特征值分解 [~, idx] max(real(diag(D))); % 广义特征值可能不排序必须手动取最大 w V(:, idx); w sqrt(P) * w / norm(w); % 满足功率约束 end这里的正则项reg很关键。当窃听者天线数小于基站天线数时B 可能是奇异矩阵直接调用eig(A, B)会得到 Inf 或 NaN。加入1e-6 * eye(Nt)后数值上稳定且对广义特征向量方向的影响可忽略。取最大广义特征值的原理是目标函数近似等于 w^H A w / w^H B w 的对数形式而广义特征向量正是使该比值最大的方向。需要提醒的是这个解是对原保密率函数的一种近似。严格的最大化 log(1γ_u) − log(1γ_e) 并不会有这么简单的闭式解。工程脚本里先用广义特征值把上半轮跑通再在调优阶段换成 MM 算法或 CVX 求解是常见的演进路径。3.3 IRS 相位的逐坐标扫描更新代码IRS 相位更新是本标题的核心。假设当前迭代已经有了一组 θ 和 w接下来要对第 n 个反射单元进行坐标更新。为了降低实现难度采用离散相位扫描在 [0, 2π) 内取 L 个候选相位对每个候选相位重新计算有效信道矩阵和保密率选出使保密率最大的相位。function theta coordinate_descent(theta, n, w, H_ru, H_re, G, H_bu, H_be, L) % 固定其他反射单元只更新第 n 个反射单元 % 返回更新后的 theta 向量 ru_n H_ru(:, n); % Nu x 1 re_n H_re(:, n); % Ne x 1 g_n G(n, :); % 1 x Nt % 先把第 n 项从有效信道中剔除 theta_off theta; theta_off(n) 0; H_u_off H_bu (H_ru .* theta_off.) * G; H_e_off H_be (H_re .* theta_off.) * G; % 候选相位 ph exp(1j * 2 * pi * (0:L-1) / L); best_R -inf; best_phase theta(n); for k 1:L H_u H_u_off (ru_n * ph(k)) * g_n; H_e H_e_off (re_n * ph(k)) * g_n; gamma_u real(w * (H_u * H_u) * w); gamma_e real(w * (H_e * H_e) * w); R log2(1 gamma_u) - log2(1 gamma_e); if R best_R best_R R; best_phase ph(k); end end theta(n) best_phase; end这段代码里H_u_off和H_e_off相当于挖去了第 n 个反射单元贡献后的信道矩阵。候选相位ph(k)是通过(ru_n * ph(k)) * g_n加回去的乘号顺序对应H_ru(:,n) * θ_n * G(n,:)的秩一结构。每个候选相位只需要两次矩阵乘法整个坐标下降一轮的复杂度大约是 O(N · L · max(Nu,Ne) · Nt)对 N64、L64 的小仿真完全可接受。如果想把扫描改成连续优化可以用 MATLAB 优化工具箱里的fminbnd对相位角直接搜索。但需要注意保密率在相位角上可能有多个峰值fminbnd只能收敛到局部极值。坐标扫描的方式天然支持多峰问题而且可以直接模拟低分辨率 IRS 单元比如每单元只有 2 或 3 bit 相移这和硬件实现更贴近。4. 仿真中的收敛判定与三个常见坑4.1 主循环、收敛曲线与 MATLAB 优化工具箱对比主脚本的职责是串起信道生成、波束成形、坐标下降和曲线绘制。下面给出核心循环的写法% 参数设置 Nt 4; Nu 2; Ne 2; N 32; P 1; sigma2 10^(-snr_dB/10); max_iter 40; L 64; % 生成信道 [G, H_ru, H_re, H_bu, H_be] gen_irs_channels(...); % 初始化 theta exp(1j * 2 * pi * rand(N, 1)); w randn(Nt, 1); w w / norm(w); R_history zeros(max_iter, 1); for iter 1:max_iter H_u H_bu (H_ru .* theta.) * G; H_e H_be (H_re .* theta.) * G; w update_beamformer(H_u, H_e, P); for n 1:N theta coordinate_descent(theta, n, w, ... H_ru, H_re, G, H_bu, H_be, L); end H_u H_bu (H_ru .* theta.) * G; H_e H_be (H_re .* theta.) * G; gamma_u real(w * (H_u * H_u) * w) / sigma2; gamma_e real(w * (H_e * H_e) * w) / sigma2; R_history(iter) max(0, log2(1 gamma_u) - log2(1 gamma_e)); end收敛曲线通常会在前 5 到 15 轮快速上升后面进入慢速微调阶段。判断收敛的标准不是保密率不再变化而是连续若干轮变化量小于某个阈值比如abs(R_history(iter) - R_history(iter-1)) 1e-4。如果曲线出现抖动下降几乎可以肯定是坐标扫描代码里把某个索引写错了因为坐标下降的单调性应该保证曲线不上扬。如果想偷懒用 MATLAB 优化工具箱里的fmincon或ga做对比建议只跑 N8 的小规模场景否则在 N64 时求解质量不如坐标扫描稳定。坐标下降的价值在于它在反射单元数量变大后仍然保持接近线性的复杂度而全局优化方法会在 N 超过 20 之后明显变慢。4.2 单位模约束、扫描点数和初始化陷阱第一个常见坑是在更新相位后忘记归一化。theta必须始终保持abs(theta) 1但浮点运算里ph(k)本身是单位模所以只要不直接对 theta 做加法就不会破坏约束。如果写的是theta(n) theta(n) delta再投影到单位圆性能会明显下降。第二个坑是扫描点数 L 的选择。L 太小时相位量化误差会直接抬高误差地板。L16 时每个相位间隔 22.5°保密率可能损失 0.3 bit/s/HzL64 时间隔约 5.6°精度已经接近连续相位下界。如果仿真目标是论文图建议 L 取 128并把坐标下降轮数降到 20总时间仍在可接受范围。第三个坑是初始化。随机相位初始化对坐标下降来说足够但为了获得更好的局部最优最好跑 3 到 5 个随机种子取保密率最高的一组作为最终结果。每次随机种子重置后信道也应重新生成否则会高估 IRS 的增益。另一个稍反直觉的结论是从所有相位都为 0 初始化往往会收敛到较差的局部极值因为 IRS 的镜面反射方向被固定坐标更新后很难越过相位势垒。随机初始化反而能覆盖更多方向。5. 验证代码正确性单步扰动检验与 MIMO 信道容量图像坐标下降代码写完后不要直接上大参数仿真先用一个“单步扰动检验”确认逐坐标更新的逻辑没写反。具体做法是在任意一次迭代后固定当前 w对第 n 个反射单元手动遍历全部 L 个候选相位画出保密率与 θ_n 角度的关系曲线并标出算法最终选择的点。如果选出的点确实落在曲线的最高峰上说明该维度的更新方向正确如果落在波谷多半是H_u_off构造时把第 n 项的错误符号留下了常见原因是没有将theta_off(n)置零。还可以进一步做穷举对比当 N 很小时比如 N4L8遍历全部 8^44096 种相位组合计算每种组合下的保密率上界。坐标下降得到的结果与穷举结果之间的差距如果小于 0.05 bit/s/Hz说明整个迭代框架没有结构性问题。这是我最常用的一招也是排查信道矩阵维度错误最有效的方法。验证通过后可以绘制固定某个反射单元相位变化时的保密率曲线俗称“相位扫描图”。代码片段如下phase_range linspace(0, 2*pi, 128); R_array zeros(size(phase_range)); for k 1:length(phase_range) theta_test theta; % 使用算法收敛后的相位 theta_test(1) exp(1j * phase_range(k)); H_u H_bu (H_ru .* theta_test.) * G; H_e H_be (H_re .* theta_test.) * G; % 固定 w 计算保密率 R_array(k) compute_rate_at_w(w, H_u, H_e, sigma2); end plot(phase_range/pi, R_array, LineWidth, 1.2); xlabel(Phase / \pi); ylabel(Secrecy rate (bit/s/Hz));这张图本质上就是你说的“MIMO 信道容量图像”的一种变体它不是原始信道容量而是保密率对相移的响应面。通过观察这条曲线的峰值位置可以反推该反射单元在级联信道中贡献的相位偏置也能解释为什么坐标下降总是先更新具有最大信道增益的反射单元。实际工程里我还会把每个单元的“相位-保密率”曲线叠加在同一张图上找出贡献最小的几个单元并在后续硬件控制中跳过它们进而降低 IRS 配置开销。这种优化技巧在反射单元数量超过 256 时尤其有价值。本文还有配套的精品资源点击获取