本研究提出的考虑硬件损耗的IRS辅助MIMO的SWIPT系统如图1．假设共有K_I个信息接收用户和K_E个能量收集器，基站端的发射天线数量N_b>1，信息接收用户和能量收集器的天线数量分别为N_I>1和N_E>1，信息接收者与能量收集器的集合为K_I和K_E，并在能量收集器附近部署M个IRS ，用于扩展能量收集器的作用范围，或者增强远处信息接收者的信号强度．

图1 IRS辅助的SWIPT通信系统Fig. 1 IRS aided SWIPT communication system.

假设每个信息接收用户接收到的数据流总量为d，满足1≤d≤min (N_b，N_E），基站的发射信号为

其中，x_k∈C^d×1为第k个信息接收用户的d×1元素数据符号向量，满足期望E [ x_k x ^H_k ]=I _d，且当i≠j (1≤i≤j≤d )时，E [ x_i x ^H_j ]=0，H为求矩阵的厄米特共轭；f_k∈C^N_b×d是基站向第k个信息接收用户发射的线性预编码矩阵；m_b是基站端由于硬件损耗产生的失真噪声，每个天线处的失真噪声功率都与其发射信号功率成正比，服从高斯随机分布m _b~N ( 0，μ _b diag ( f_k f ^H_k ) )，m _b∈C^N_d×1，μ _b为基站端的硬件损耗系数．基站到IRS、基站到第k个信息接收用户、基站到第l个能量收集器、IRS到第k个信息接收用户、IRS到第l个能量收集器的信道分别建模成矩阵h_b，r∈C^M×N_b，h_b，k∈C^N_I×N_b，z_b，l∈C^N_E×N_t， h _r，k∈C ^N _I×M 和 z _r，l∈C ^N _E×M．I R S 的反射矩阵 Φ=diag ( e^iθ₁，e^iθ₂，⋯，e^iθ_m，⋯，e^iθ_M )，其中， θ_m为第m个反射元件的相移角， θ_m∈[ 0 ，2 π ]．本研究只考虑信号的1次反射，多次信号反射忽略．第k个信息接收用户的接收信号为

y _I，k=( h _b，k+h _r，kΦh _b，r ) s+n _I，k+m_i（2）

其中，n_I，k是第k个信息接收用户处的加性白高斯噪声，服从N ( 0 ，σ ²_I I _{N I} )的正态分布， σ ²_I为噪声功率， I_NI为N_I阶的单位矩阵．定义第k个信息接收用户接收到的无噪声干扰的信号矩阵为yˉ_I，k=( h_b，k+h _r，kΦh _b，r ) s+n _I，k．其中，m_i为信息接收用户i处因硬件损耗产生的失真噪声，服从 分布， μ_i为信息接收用户处的硬件损耗系数^{[2 3]}．同理，第l个能量收集器的接收信号为

y _E，l=( z _b，l+z _r，lΦh _b，r ) s+n _E，l+m _E（3）

其中，n_E，l是第l个能量收集器处的加性高斯白噪声，服从正态分布N ( 0 ， σ ²_E I_{N E} ) ， σ ²_E为噪声功率， I_NE为N_E阶的单位矩阵．第l个能量接收器接收的无噪声干扰信号矩阵为yˉ_E，l=( z _b，l+z _r，lΦh _b，r ) s+n _E，l．其中，m_E为能量收集器处因硬件损耗而产生的失真噪声，服从m_E~ 分布，μ_E为能量收集器处的硬件损耗系数．将式（1）代入式（2），可得

其中，H_k为基站发射给第k个信息接收用户的级联信道， H_k=h _b，k+h _r，kΦh _b，r．

第k个信息接收用户的可实现数据传输速率为

其中，J_k是第k个用户的干扰加噪声的协方差矩阵，

本研究采用一种线性模型^[25]作为能量接收模型，通过忽略能量收集器处的噪声干扰，使能量收集器的总收集功率与总接收功率成比例．将基站发射给第l个能量收集器的级联信道定义为Z_l=z _b，l+z _r，lΦh _b，r . 第l个能量收集器总收集功率为

其中，η为能量转换效率，且0<η<1．本研究考虑能量收集器的最小能量约束，即

其中，α_l为第l个能量收集器的能量加权因子，α_l越大第 l 个能量收集器的优先级越高；Zˉ = ；Q₀为最小能量收集阈值．

本研究在基站最大发射功率、能量收集器最小收集能量以及IRS的单位模约束下，对基站的预编码矩阵f和IRS的相移矩阵Φ联合优化，最大化信息接收用户的加权和速率（weighted sum rate， WSR），即

其中，ϕ_m=e^iθ_m，对∀m有Φ=diag (ϕ₁，ϕ₂，⋯，ϕ_M )；ω_k为加权系数， 0<ω_k<1， ω_k越大第k个信息接收用户优先级越高；R_k ( f，Φ )为第k个信息接收用户的可实现数据传输速率；P_T为基站的最大发射功率；约束（8d）为相移角的单位模约束．

为将原问题中高度耦合地优化变量解耦，采用最小均方误差的方法来简化原问题^[26]．引入信息接收用户端的解码矩阵U={U_k，∀k∈K_I}，则每个信息接收用户的接收信号矩阵为

其中， U_k为第k个信息接收用户的解码矩阵， U_k∈C^N_I×d．第k个信息接收用户的均方误差矩阵为

引入辅助矩阵W={W_k≥0，∀k∈K _I}，则问题（8）可转换^[23]为

其中，h_k ( W，U，f，Φ )=lg|W_k|-tr ( W_k E_k )+d．

转换后的问题（11）虽然增加了优化变量，但比问题（8）更易处理．利用BCD算法来解决此问题．首先，对于给定的W_k、Φ和f，使h_k(W，U，f，Φ)对U_k求一阶导数，并令其一阶导数等于0，可得U_k的最优解为

同理，固定U_k、Φ和f，使h_k(W，U，f，Φ)对W_k求一阶导数，并令其等于0，可得到W_k的最优解为

本研究通过固定W、U、Φ来优化预编码矩阵f．将式（1 0）代入问题（1 1），并将f的无关项去掉，可得

其中， ．

由于目标函数和最小能量约束都是非凸的，为将问题（14）转换成一个等价的凸问题，利用不等式^[27]CA ^H diag ( BB ^H ) A≥B ^H diag ( ACA ^H ) B ，可将目标函数转换成凸形式，即将问题（14）中的目标函数转换为

接着采用连续凸逼近（successive convex approximation，SCA）方法解决最小能量约束问题^[28]．通过一阶泰勒展开来近似能量约束，可得

其中，f ⁽ⁿ⁾k 为第n迭代所得的解．能量约束（8c）可写成

其中， ．问题（14）可转换为

经此将问题（18）转变为凸问题，再采用拉格朗日对偶乘子算法来求解^[29]．由于问题（18）是一个凸问题，满足强对偶条件，且对偶间隙为0，所以可通过解决其对偶问题来获得最优解．引入功率约束相关的拉格朗日乘子λ，则问题（18）的部分拉格朗日函数为

对偶函数F ( λ)为

F ( λ )= m_fin L ( f，λ )， s. t.式(17 )（20）

则对偶问题可表示为

在解决式（21）前，需先求解F ( λ)的表达式．通过引入能量约束（17）的拉格朗日乘子μ，可得问题（20）的拉格朗日函数为

拉格朗日函数L ( f，μ )对f_k求一阶导数，并令其等于0，可得到f_k的最优解为

其中，μ满足约束（17）的互补松弛条件，即

当μ=0时，如果满足下列条件

则式（23）的最优解为f ^⋆_k (0)，∀k∈K_I．当μ≠0时，由式（23）和式（24）求得μ的最优解为

然后，对对偶问题（21）进行求解，找到最优的λ．λ应满足约束（17）的互补松弛条件，即

当λ=0时，如果满足下列条件

则可得到对偶问题的最优解为f_k (0)．若不满足条件（28），则假设在λ处求得问题（21）的最优解，有

由于μ的表达式过于复杂，难以得到P ( λ)的表达式，本研究通过证明可得到P ( λ)是一个单调递减函数（证明过程请扫描论文末页右下角二维码查看补充材料证明S1），再使用二分法算法可找到最优的λ．二分法算法详情请扫描论文末页右下角二维码查看补充材料图S1．图S1对λ>0的情况提供了解决问题（18）的详细步骤．

图2给出了SCA算法的详细步骤，解决了问题（14），并得到了基站预编码矩阵的最优解．

图2 SCA算法伪代码Fig. 2 SCA algorithm pseudocode.

定理1由算法1产生的{ f ⁽ⁿ⁾，n=1，2，…}序列会收敛于问题（14）的卡罗需-库恩-塔克（Karush-Kuhn-Tucker，KKT）最优点．

【证】 与文献[29]证明相似，在此不再赘述．

在本小节中，将固定其他优化变量来优化相移矩阵Φ．将式（10）带入问题（11），并去除Φ的无关项，可得相移角的优化问题

其中，

将H_k=h _b，k+h _r，kΦh _b，r带入问题（3 0），可得

W_kU ^H_k H_k f_k=W_kU ^H_k h _r，kΦh _b，r f_k+W_kU ^H_k h _b，k f_k（32）

定义B_k=ω_kh ^H_r，kU_kW_kU ^H_k h _r，k，C ₁=h _b，r F͂h ^H_b，r ，D_k=ω_kh _b，r F͂h ^H_b，kU_kW_kU ^H_k h _r，k 和 T_k=ω_kh _b，r f_kW_kU ^H_k h _r，k．式（31）和式（32）等价为

其中，res₁和res₂都为与Φ无关的常数，可舍去．

将Z_l=z _b，l+z _r，lΦh _b，r带入式（8 c），可得

定 义 ，将约束（8c）转换成等价形式

定义 ，将式（33）和式（34）代入到问题（30），则问题（30）转换成等价形式

其中， ．

为更好地处理优化问题，将相移矩阵Φ的对角元素用矩阵ϕ=[ϕ₁，ϕ₂，⋯，ϕ_M ]^T表示，取矩阵恒等式^[28]，则有

这里，符号∘为哈达玛积.

将V和Z_b，r的元素对应到对角矩阵中，如v=[ V _1，1，V _2，2，⋯，V_M，M ] ^T 和 z=[ ( Z _b，r ) _1，1，( Z _b，r ) _2，2，⋯， ( Z_b，r )_M，M ]^T，可得

能量约束（36）可写成

其中， χ=Z_r∘C₂； Q₂=Q₀-tr (Z_bFˉ)；为将能量约束式（36）变成凸形式，可使用SCA算法．在ϕ⁽ⁿ⁾处将式子（36）一阶泰勒展开，可得

其中，（n）表示第n次迭代；Q₃=Q₂+[ ]ϕ⁽ⁿ⁾ H χϕ⁽ⁿ⁾．

因此，问题（37）可改写为

其中， ；v^*为v的共轭矩阵．

接下来，利用最优最小化（majorization-minimization，MM）算法^[30]求解优化问题（42）．MM算法包括2个步骤：①找到一个局部逼近目标函数的替代函数，且在当前点处两个函数的差最小．②最小化替代函数．在替代函数达到最小值处，重复步骤①和②，直到目标函数与替代函数在当前点的值低于一个阈值，就可得到目标函数的最小值．

MM算法最重要的步骤在于如何找到一个合适的替代函数，为此，本研究引入引理1^[31]．

引理1 对于给定的ϕ⁽ⁿ⁾，对于任意的ϕ都满足不等式

其中， X = λ _max I _M ， λ _max是γ的最大特征值； I_M为M阶单位矩阵．

定义y ( |ϕ ϕ⁽ⁿ⁾ )为替代函数，有

在替代函数y ( |ϕ ϕ⁽ⁿ⁾ )中，ϕ^H Xϕ=λ_max M为常量，可舍去，问题（42）可改写为

其中， p ^{( n )} =( λ _{m a x} I_M-γ )_ϕ^{( n )} -v^*，由于约束（4 1），和约束（8d）是非凸的，导致问题（45）也是非凸问题，因此不能使用的格朗日方法来求解．

文献[32]提出了一种等价方法，在约束（41）左边引入非负的辅助矢量q．

对于给定的q，该优化问题的全局最优解为

任意q都应当满足约束（41）的互补松弛条件

q(J(q)-Q₃)=0（48）

其中，J (q)=2Re éëϕ(q) ^H( z^*+χϕ⁽ⁿ⁾ ) ù 当q=0时û若满足条件

J(0)≥Q₃（49）

则有q^*=0．若不满足式（49），则q^*满足J ( q^*)=Q₃．可以证明，J (q)是一个单调递增的函数（证明与补充材料S1类似，在此不再赘述），再通过二分法算法求得q^*（详情请扫描论文末页右下角二维码查看补充材料图S2），提供了算法步骤，解决了问题（46）．

MM联合SCA算法详情如图3．

图3 MM算法伪代码Fig. 3 MM algorithm pseudocode.

图4提供了全局性BCD算法求解问题（8）的详细步骤．其中，R ( f ⁽ⁿ⁾，ϕ⁽ⁿ⁾ )为问题（8）目标函数第n次迭代的值.

定理 2 由算法 3 产生的{R ( f ⁽ⁿ⁾，ϕ⁽ⁿ⁾ )，n=1，2，⋯}序列会收敛于问题（14）的KKT最优点．

【证】 与文献[32]证明相似，在此不再赘述．

图4 BCD算法伪代码Fig. 4 BCD algorithm pseudocode.