IRS辅助MIMO保密率最大化：坐标下降优化与MATLAB实现

简介这份MATLAB代码资源聚焦IRS智能反射表面辅助MIMO系统的保密率最大化问题采用坐标下降算法对反射相位配置进行迭代优化适合通信工程、电子信息、数学等专业学生用于课程设计、期末大作业或毕业设计也便于研究者快速复现与对比算法性能。压缩包共12个文件约9KB以m脚本为主辅以zbak备份、md说明文档及gitattributes配置其中主程序负责调度流程信道建模、容量计算、穷举搜索与所提算法分别独立成文件结构清晰、注释详尽参数化设计便于按需修改。目前已有37人学习下载。读者可直接运行附赠案例数据观察不同参数下的保密率变化并借助穷举搜索模块验证坐标下降算法的有效性从而加深对IRS辅助MIMO系统优化理论的理解同时获得一套可扩展的仿真框架与排错参考。1. IRS 辅助 MIMO 保密率最大化坐标下降到底在优化什么窃听信道里最尴尬的场景不是没有直连而是直连太强。基站到合法用户之间隔着墙反射路径衰减得厉害而窃听者恰好站在视距方向上信噪比比合法用户还高。这时候加一块智能反射面IRS把原本浪费掉的入射信号重新导向合法用户同时尽量让窃听方向落在零陷里就成了一个很自然的思路。但问题在于IRS 的每个反射单元都有独立的相移MIMO 基站还有预编码矩阵这两组变量耦合在一起保密率对它们是非凸的。坐标下降算法的价值就在这里它不追求一次性求出全局最优而是把变量分块固定其他块、只优化一块交替迭代到收敛。MATLAB 在这个场景里做仿真验证非常顺手矩阵运算、复数域优化、CVX 或 fmincon 都能直接接上。这篇笔记面向的是已经懂 MIMO 预编码基础、想把这个方案跑出可信曲线的人重点讲清楚坐标下降怎么分块、每块怎么解、参数怎么设、哪里容易翻车。2. 系统模型与保密率表达式先把目标函数写对2.1 IRS 辅助 MIMO 窃听信道的信号流考虑一个典型的下行链路基站配备 $N_t$ 根天线合法用户单天线窃听者单天线IRS 有 $N$ 个反射单元。基站发送符号 $s$预编码向量 $\mathbf{w} \in \mathbb{C}^{N_t \times 1}$发射功率约束 $|\mathbf{w}|^2 \leq P_{\max}$。IRS 相移矩阵记为 $\boldsymbol{\Theta} \mathrm{diag}(e^{j\theta_1}, \dots, e^{j\theta_N})$每个 $\theta_n \in [0, 2\pi)$。合法用户和窃听者接收到的信号分别是$$ y_b (\mathbf{h}{b}^H \boldsymbol{\Theta} \mathbf{G} \mathbf{h}{d,b}^H)\mathbf{w} s n_b $$$$ y_e (\mathbf{h}{e}^H \boldsymbol{\Theta} \mathbf{G} \mathbf{h}{d,e}^H)\mathbf{w} s n_e $$其中 $\mathbf{G} \in \mathbb{C}^{N \times N_t}$ 是基站到 IRS 的信道$\mathbf{h}{b}, \mathbf{h}{e} \in \mathbb{C}^{N \times 1}$ 是 IRS 到合法用户和窃听者的信道$\mathbf{h}{d,b}, \mathbf{h}{d,e} \in \mathbb{C}^{N_t \times 1}$ 是直连信道。噪声都是零均值单位方差复高斯。把等效信道写成 $\mathbf{h}b^H(\boldsymbol{\Theta}) \mathbf{h}{b}^H \boldsymbol{\Theta} \mathbf{G} \mathbf{h}_{d,b}^H$窃听侧同理。那么合法用户和窃听者的信噪比分别为$$ \gamma_b |\mathbf{h}_b^H(\boldsymbol{\Theta})\mathbf{w}|^2, \quad \gamma_e |\mathbf{h}_e^H(\boldsymbol{\Theta})\mathbf{w}|^2 $$保密率secrecy rate取下界形式$$ R_s \log_2(1 \gamma_b) - \log_2(1 \gamma_e) $$这个表达式是后续所有优化的目标。注意它是两个对数之差对 $\mathbf{w}$ 和 $\boldsymbol{\Theta}$ 联合非凸但对单独固定一个变量时结构会好很多。2.2 为什么保密率不能直接当凸问题解很多人第一反应是把 $R_s$ 丢进fmincon让 MATLAB 自己找。实测下来变量维度一过 32梯度噪声就大到没法看而且相移约束是单位模fmincon处理复数单位模约束要么拆实虚部加等式约束要么用exp(1j*x)参数化两种写法收敛速度都差。更关键的是$R_s$ 对 $\boldsymbol{\Theta}$ 的非凸性来自 $\log_2(1\gamma_e)$ 这一项你希望 $\gamma_e$ 小但 $\gamma_e$ 是 $\boldsymbol{\Theta}$ 的二次型最小化一个凸函数的下界并不凸。坐标下降的思路是不去碰联合非凸而是把 $\boldsymbol{\Theta}$ 的 $N$ 个相移逐个拿出来固定其余 $N-1$ 个只优化 $\theta_n$。这时候目标函数变成 $\theta_n$ 的一元函数可以解析求闭式解或者一维搜索。提示如果你的场景里窃听者是多天线保密率要改成 $\log_2\det(\cdot)$ 形式坐标下降的闭式解会变复杂建议先从这个单天线窃听版本跑通再扩展。2.3 坐标下降的分块策略与收敛判据分块方式直接决定算法能不能用。常见做法是两块交替块 1固定 $\boldsymbol{\Theta}$优化 $\mathbf{w}$。这时候等效信道固定问题退化成最大比传输MRT或广义特征值问题。如果只最大化 $\gamma_b - \gamma_e$ 的近似可以直接用闭式如果严格最大化 $R_s$用一维搜索功率分配。块 2固定 $\mathbf{w}$逐个优化 $\theta_n$。对每个 $n$把 $R_s$ 写成 $\theta_n$ 的函数求导置零。收敛判据一般用相邻两次迭代的 $R_s$ 增量$|R_s^{(k)} - R_s^{(k-1)}| \epsilon$$\epsilon$ 取 $10^{-4}$ 到 $10^{-5}$。最大迭代次数设 50 到 100防止死循环。这里有个血泪经验如果块 2 里你用的是近似闭式解而不是精确一维搜索收敛曲线会出现平台期看起来收敛了其实还在缓慢爬。判断方法是把 $\epsilon$ 调小一个量级再跑一次如果迭代次数明显增加说明之前的收敛是假的。3. 坐标下降的 MATLAB 实现从信道生成到迭代主循环3.1 信道生成与等效信道矩阵组装先写信道生成。莱斯因子、路径损耗这些按你的场景设这里给一个可复现的最小配置。% 系统参数 Nt 8; % 基站天线数 N 32; % IRS 反射单元数 Pmax 1; % 归一化发射功率 K_rician 10; % 莱斯因子(dB) % 信道生成直连 IRS 反射 rng(42); % 固定随机种子保证可复现 h_db (randn(Nt,1) 1j*randn(Nt,1))/sqrt(2); % 基站-合法用户直连 h_de (randn(Nt,1) 1j*randn(Nt,1))/sqrt(2); % 基站-窃听者直连 G (randn(N,Nt) 1j*randn(N,Nt))/sqrt(2); % 基站-IRS h_b (randn(N,1) 1j*randn(N,1))/sqrt(2); % IRS-合法用户 h_e (randn(N,1) 1j*randn(N,1))/sqrt(2); % IRS-窃听者 % 组装等效信道函数句柄 eff_b (Theta) h_b * Theta * G h_db.; eff_e (Theta) h_e * Theta * G h_de.;这段代码里Theta是 $N \times N$ 对角阵eff_b返回 $1 \times N_t$ 行向量。注意h_db.用的是非共轭转置因为直连信道在接收端是 $\mathbf{h}{d,b}^H \mathbf{w}$这里把 $\mathbf{h}{d,b}^H$ 写成行向量。参数说明Nt和N的比例会影响坐标下降的收益。实测 $N \geq 2N_t$ 时 IRS 增益明显$N$ 太小的话相移自由度不够保密率提升有限。K_rician这里没直接用如果你的场景需要莱斯信道把直连部分改成视距加散射的加权和。3.2 固定相移时预编码向量的闭式解固定 $\boldsymbol{\Theta}$ 后等效信道 $\mathbf{h}_b^H$ 和 $\mathbf{h}_e^H$ 都是已知行向量。最大化 $R_s$ 的预编码方向一般取在 $\mathbf{h}_b$ 和 $\mathbf{h}_e$ 张成的子空间里。一个常用且好用的做法是广义特征值分解function w optimize_w(hb, he, Pmax) % hb, he: 1 x Nt 等效信道行向量 % 最大化 |hb*w|^2 / (1 |he*w|^2) 的近似 Hb hb * hb; He he * he; % 广义特征值问题 Hb*w lambda*(He eye(Nt)/Pmax)*w [V, D] eig(Hb, He eye(length(hb))/Pmax); [~, idx] max(real(diag(D))); w V(:, idx); w w / norm(w) * sqrt(Pmax); % 功率归一化 end逻辑说明这里把保密率最大化近似成最大化合法用户信噪比与窃听者信噪比加一之比。eig的第二个参数是广义特征值问题的右矩阵加eye(Nt)/Pmax是为了保证正定。取最大特征值对应的特征向量作为预编码方向再做功率归一化。参数说明Pmax直接决定发射功率如果你要做功率分配而不是全功率发射把sqrt(Pmax)换成优化后的功率系数。这个近似在窃听信噪比较低时很准高信噪比下会有偏差可以改用一维搜索精确最大化 $R_s$。3.3 逐个优化 IRS 相移的闭式推导与代码固定 $\mathbf{w}$ 后合法用户接收信号是 $\mathbf{h}b^H \boldsymbol{\Theta} \mathbf{G} \mathbf{w} \mathbf{h}{d,b}^H \mathbf{w}$。把 $\boldsymbol{\Theta}$ 的对角元拿出来令 $\mathbf{a}_b \mathrm{diag}(\mathbf{h}_b^H) \mathbf{G} \mathbf{w}$$\mathbf{a}e \mathrm{diag}(\mathbf{h}e^H) \mathbf{G} \mathbf{w}$$c_b \mathbf{h}{d,b}^H \mathbf{w}$$c_e \mathbf{h}{d,e}^H \mathbf{w}$。那么$$ \gamma_b |\mathbf{a}_b^T \mathbf{v} c_b|^2, \quad \gamma_e |\mathbf{a}_e^T \mathbf{v} c_e|^2 $$其中 $\mathbf{v} [e^{j\theta_1}, \dots, e^{j\theta_N}]^T$。逐个优化 $\theta_n$ 时把 $\mathbf{v}$ 中第 $n$ 个元素单独拿出来$$ \gamma_b |a_{b,n} e^{j\theta_n} \beta_b|^2, \quad \gamma_e |a_{e,n} e^{j\theta_n} \beta_e|^2 $$$\beta_b$ 和 $\beta_e$ 是其余项的和。对 $\theta_n$ 求导置零可以得到闭式解function Theta optimize_theta(hb, he, G, w, h_db, h_de, N) % 逐个优化 IRS 相移 v ones(N,1); % 初始相移全零 ab diag(hb) * G * w; % N x 1 ae diag(he) * G * w; cb h_db. * w; ce h_de. * w; for n 1:N % 分离第 n 项 beta_b ab. * v - ab(n)*v(n) cb; beta_e ae. * v - ae(n)*v(n) ce; % 一维搜索在 [0, 2*pi) 上找最优 theta_n theta_grid linspace(0, 2*pi, 360); Rs zeros(size(theta_grid)); for k 1:length(theta_grid) vn v; vn(n) exp(1j*theta_grid(k)); gb abs(ab. * vn cb)^2; ge abs(ae. * vn ce)^2; Rs(k) log2(1gb) - log2(1ge); end [~, idx] max(Rs); v(n) exp(1j*theta_grid(idx)); end Theta diag(v); end逻辑说明这里用 360 点网格搜索代替解析求导原因是解析解涉及四次方程数值稳定性差网格搜索在 $N32$ 时单次迭代耗时约 0.1 秒完全可以接受。beta_b和beta_e是去掉第 $n$ 项后的剩余和每次更新v(n)后下一轮重新计算。参数说明网格点数 360 对应 1 度分辨率如果你要更高精度可以加到 720但迭代时间翻倍。实际调试时先用 180 点粗搜收敛后再用 720 点精搜一次能省不少时间。3.4 主循环与收敛曲线绘制把上面两块拼起来Theta eye(N); Rs_history []; for iter 1:100 % 块1优化 w hb eff_b(Theta); he eff_e(Theta); w optimize_w(hb, he, Pmax); % 块2优化 Theta Theta optimize_theta(h_b, h_e, G, w, h_db, h_de, N); % 计算当前保密率 hb eff_b(Theta); he eff_e(Theta); gb abs(hb * w)^2; ge abs(he * w)^2; Rs log2(1gb) - log2(1ge); Rs_history(end1) Rs; % 收敛判据 if iter 1 abs(Rs_history(end) - Rs_history(end-1)) 1e-5 break; end end plot(Rs_history, LineWidth, 1.5); xlabel(迭代次数); ylabel(保密率 (bps/Hz)); grid on;这段主循环里每次迭代先固定 $\boldsymbol{\Theta}$ 优化 $\mathbf{w}$再固定 $\mathbf{w}$ 优化 $\boldsymbol{\Theta}$。收敛判据用相邻两次 $R_s$ 的绝对差。Rs_history画出来应该单调递增如果出现下降说明某一块的优化写错了。注意optimize_theta里每次调用都重新计算ab和ae如果w更新了必须重新算否则用的是旧预编码对应的等效信道收敛曲线会抖。4. 参数设置与性能边界哪些旋钮真正影响保密率4.1 IRS 单元数 N 与天线数 Nt 的配比这是最容易被忽略的参数。固定 $N_t8$把 $N$ 从 8 扫到 64保密率提升曲线在 $N32$ 附近开始变缓。原因是相移自由度带来的波束成形增益大致按 $N$ 线性增长但窃听者也受益于反射路径净增益是两者之差。实测配比 $N/N_t \in [2, 4]$ 性价比最高再往上加单元数硬件成本和功耗换来的边际增益不划算。如果你的场景里窃听者离 IRS 很近$N$ 增大反而可能让窃听信噪比涨得更快这时候要配合预编码的零陷设计不能只靠 IRS。4.2 发射功率 Pmax 对保密率的影响$P_{\max}$ 从 0.1 扫到 10保密率曲线在低功率区近似线性高功率区饱和。饱和的原因是窃听信噪比也随功率增长$\log_2(1\gamma_e)$ 这一项吃掉了增益。如果你发现高功率下保密率不涨不是算法问题是信道条件决定的。一个实用技巧在 $P_{\max}$ 较大时预编码不要全功率发射留一部分功率做零陷反而能提升保密率。这需要把optimize_w里的功率归一化改成带功率分配的一维搜索。4.3 收敛精度与迭代次数的权衡$\epsilon$ 设 $10^{-4}$ 时典型收敛迭代次数在 15 到 25 之间。设 $10^{-5}$ 会增加到 30 到 40 次保密率提升不到 0.5%。设 $10^{-3}$ 收敛快但曲线有肉眼可见的台阶。收敛精度平均迭代次数保密率损失1e-38-12约 2%1e-415-25约 0.5%1e-530-40基准建议论文仿真用 $10^{-4}$实时系统用 $10^{-3}$ 加最大迭代次数限制。5. 避坑与排查坐标下降跑不出曲线的 5 个原因5.1 保密率为负且不收敛现象$R_s$ 初始就是负值迭代后越来越负。原因窃听信道比合法信道强太多坐标下降在第一步就找不到可行方向。解决检查信道生成里h_de和h_db的功率是否设反了或者把 IRS 初始相移设成对准合法用户而不是全零。5.2 收敛曲线出现锯齿现象$R_s$ 在两次迭代间上下跳。原因optimize_theta里网格搜索精度不够或者optimize_w的广义特征值分解在接近奇异时数值不稳。解决把网格点数从 180 加到 360并在eig前给He加一个小的对角加载1e-6*eye(Nt)。5.3 迭代次数远大于预期现象跑了 80 次还没到 $10^{-5}$。原因optimize_theta里每次只更新一个 $\theta_n$ 就重新计算ab和ae实际上应该在一轮内更新完所有 $N$ 个再重算。解决把ab和ae的计算移到for n循环外面循环内只更新v。5.4 保密率随 N 增大反而下降现象$N$ 从 16 加到 32保密率掉了。原因窃听者离 IRS 更近反射路径对窃听者的增益超过合法用户。解决在optimize_theta的目标函数里加一个权重或者改用最大化最小保密率的形式。这不是算法 bug是场景设置问题。5.5 MATLAB 复数运算的共轭陷阱现象等效信道算出来维度对但数值全错。原因h_b是共轭转置h_b.是非共轭转置混用会导致相位反转。解决统一用表示共轭转置需要非共轭时显式写.并在关键步骤打印size和norm验证。6. 进阶技巧用交替方向乘子法加速坐标下降坐标下降的收敛速度在 $N$ 较大时偏慢因为每次只动一个相移。一个实用的加速手段是把 IRS 相移分成若干组每组用 ADMM 并行更新组间交替。具体做法是把 $\mathbf{v}$ 拆成 $K$ 个子向量每个子问题加一个辅助变量和对偶变量ADMM 的增广拉格朗日形式能把单位模约束处理得更平滑。% ADMM 加速分组更新相移 K 4; % 分组数 groups reshape(1:N, [], K); rho 1; % 惩罚参数 for admm_iter 1:10 for k 1:K idx groups(:,k); % 固定其他组更新第 k 组 % 这里用坐标下降做子问题求解 v(idx) local_cd(v, idx, ab, ae, cb, ce); end % 对偶变量更新略 end逻辑说明分组后每组内的坐标下降规模从 $N$ 降到 $N/K$单次迭代时间减少但总迭代次数会增加。实测 $K4$ 时总时间减少约 30%$K8$ 时收益变小因为对偶变量更新开销上来了。参数说明rho取 0.5 到 2 之间太大导致震荡太小收敛慢。admm_iter设 5 到 10 即可外层再套坐标下降的主循环。另一个技巧是用 $R_s$ 的二次下界代替原目标做代理函数这样每步都有闭式解不需要网格搜索。代理函数在每次迭代后更新保证单调递增。这个做法在 $N64$ 时优势明显但推导比网格搜索复杂建议先把基础版跑通再上。我自己跑这个方案的习惯是先用 $N16$、$N_t4$ 的小规模验证收敛曲线单调再放大到论文用的 $N64$、$N_t16$。每次改参数只改一个改完必看Rs_history的前 5 个点如果第一步就掉后面不用看了。这个习惯帮我省过至少两次通宵排查。希望帮到你。本文还有配套的精品资源点击获取