超可靠低时延通信(URLLC，ultra-reliable and low-latency communication)需要在时延受限的条件下准确地传递信息^[1]，故可采用短数据包进行传输；由于数据包长度有限，不再使用香农定理的最大互信息算法^[2].最大可达速率(MAR, maximal achievable rate)不但与信道分布有关，而且也与块错误概率(PEP, packet error probability)有关^[3]；同时，在URLLC下，收发设备大多是能量受限的，需要采用能量捕获的方式获取能量，以供数据传输.多天线技术能够利用空间资源实现复用增益，提高短数据包传输系统的MAR以及能量传输效率，进而降低传输时延；同时也可以实现分集增益，从而降低系统的PEP.

Khan和Lopez等^[3-4]推导了点对点的单天线无线能量传输(WET，wireless energy transmission)系统中短包传输的MAR和PEP；Lopez等^[5]认为能量捕获URLLC系统的时延其与系统的PEP成反比；Yang等^[6]研究了能量捕获受限条件下系统的MAR性能，认为MAR与PEP和信道散度有关.由此可以看出，PEP直接影响短数据包传输的MAR，从而进一步影响其传输时延；同时捕获能量的多少直接影响发射功率的大小，进而影响PEP和MAR.所以，PEP能够表征短数据包传输系统的可靠性、有效性和能量捕获效率，并影响整个系统的时延.然而，以上研究都是针对理想状态的点对点通信系统，未考虑多设备连接场景. Gu等^[7]分别在时分多址接入和迫零波束成形情景下研究了多用户短包传输系统的PEP，并分别对块长度分配和功率分配系数进行优化，使系统的总PEP最小，但未考虑能量捕获问题. Hu等^[8]研究了协作多终端短包传输系统的PEP，但由于能量开销较大，无法使用WET为设备提供能量以供数据传输.由以上研究发现，多用户MIMO-WET系统中短包传输的PEP性能尚无研究.

针对此问题，在多用户MIMO-WET短包传输系统中，为能量发射器(ET，energy transmitter)和源节点配置多根天线，利用源节点捕获到的能量进行信息发送；根据源节点与用户的空间分布和迫零波束成形方法得到每个用户的信噪比累积分布函数，并利用高斯函数近似的方法，推导出每个用户的块错误概率的近似闭合表达式.利用基于块坐标下降(BCD，the block coordinate descent)和凹凸过程(CCCP，concave-convex procedure)的迭代算法，对WET阶段中的时隙数和无线信息传输(WIT，wireless information transmission)阶段中的功率分配系数因子进行联合优化，实现总块错误概率的最小化.理论分析结果表明，块错误概率与信息传输的块长度、功率分配系数成反比.仿真结果验证了所推导的块错误概率闭合表达式的正确性.根据仿真数据可知，在给定信息传输的块长度和天线数量条件下，WET和WIT所用时隙数存在最优值.

1 系统模型

将多天线技术应用到无线能量捕获短包传输系统中，系统由ET、源节点S和K个用户组成，如图 1所示.由于节点能量受限，K个用户都配备单天线，并分布于不同的位置，分别用U₁, U₂, …, U_K表示；ET的天线数量为N_ET，源节点的天线数量为N.假设能量捕获节点使用存储然后传输协议进行无线通信，整个过程分为WET阶段和WIT阶段.首先，源节点S在m个时隙中通过ET发出的射频(RF，radio frequency)信号进行能量收集，并存入电池中，即WET阶段；然后，源节点利用收集的能量，通过n个时隙与目的节点进行信息传输, 即WIT阶段.假设信道为准静态衰落信道，并且在每一个传输块中其衰落系数都固定不变.另外，假设ET到S具有完美的信道状态信息(CSI，channel state information)，每个时隙的持续时间为T_c.

1.1 WET阶段

首先考虑在ET处的能量波束成形，假设在ET处一共有d个能量波束，1≤d≤N_ET. ET处传输的能量信号为$\mathit{\boldsymbol{x}} = \sum\limits_{i = 1}^d {{\mathit{\boldsymbol{\omega }}_i}{s_i}} $，其中，ω_i∈C^M_T×1为第i个传输波束成形向量，s_i为相应的能量信号.在ET处的传输协方差矩阵为

$ \mathit{\boldsymbol{S}} = E(\mathit{\boldsymbol{x}}{\mathit{\boldsymbol{x}}^{\rm{H}}}) = \sum\limits_{i = 1}^d {{\mathit{\boldsymbol{\omega }}_i}} \mathit{\boldsymbol{\omega }}_i^{\rm{H}} \ge 0 $

(1)

源节点S捕获的能量为

$ Q = \eta LE({\left\| {\mathit{\boldsymbol{Hx}}} \right\|^2}) = \eta L{\rm{tr}} (\mathit{\boldsymbol{GS}}) $

(2)

其中：η为能量捕获效率；L为块长度；H为ET到S之间的MIMO信道矩阵，G$\buildrel \Delta \over = $H^HH.

由文献[1]得知，设一个连续时间信号的近似持续时间为t，近似带宽为B，则数据包长度为：L≈Bt.因此，L≈BmT_c.式(2)可以写为

$ Q \approx \eta Bm{T_{\rm{c}}} {\rm{tr}} (\mathit{\boldsymbol{GS}}) $

(3)

其中:S$\underline \succ $0，tr(S)≤P. Xu等^[9]通过优化的方法得到式(3)的最大值Q^*=ηBmT_cPλ₁. λ₁为矩阵G的最大特征值，即λ₁≥…≥λ_MET≥0.

1.2 WIT阶段

源节点S利用捕获的能量与K个用户进行信息传输.此阶段采用迫零波束成形方案. K个用户同时与S进行通信，数据包长度为M, M≈BnT_c，功率分配系数为β_k. P_S为源节点S的传输功率，有

$ {P_{\rm{S}}} = \frac{Q}{{n{T_{\rm{c}}}}} = \frac{{\eta Bm {\rm{tr}} (\mathit{\boldsymbol{GS}})}}{n} $

(4)

最大传输功率为

$ P_{\rm{S}}^* = \frac{{{Q^*}}}{{n{T_{\rm{c}}}}} = \frac{{\eta BmP{\lambda _1}}}{n} $

(5)

用户k的接收信号为

$ {y_k} = \sqrt {{P_{\rm{S}}}} \mathit{\boldsymbol{W}}{\mathit{\boldsymbol{h}}_k} + {n_k} $

(6)

其中：h_k~CN(0_N, Ω_kI_N), k=1, 2, …, K为S与U_k之间的N×1信道矢量；Ω_k为S与U_k的平均信道功率增益；x_k为源节点S给用户k传输的信息；n_k~CN(0，σ_k²)为加性高斯白噪声；$\mathit{\boldsymbol{W}} = \sum\limits_{k = 1}^K {\sqrt {{\beta _k}} {x_k}{\mathit{\boldsymbol{v}}_k}} $为S的1×N预编码矩阵；v_k, k=1, 2, …, K为信息x_k的一个1×N波束成形矢量；‖v_k‖=1；β_k为x_k的功率分配系数，$\sum\limits_{k = 1}^K {{\beta _k}} = 1$.由文献[5]得到U_k的信噪比(SNR，signal-to-noise ratio)为

$ {\gamma _k} = \frac{{{P_{\rm{S}}}{\beta _k}|\mathit{\boldsymbol{h}}_k^\dagger {\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}_k}{\mathit{\boldsymbol{h}}_k}|}}{{\sigma _k^2}} $

(7)

其中：Θ_k=I_N－H_k(H_k^†H_k)^－1H_k^†；H_k=[h₁, h₂, …, h_K].

1.3 性能分析

源节点S以固定的速率r进行信息传输，系统的性能主要由块错误概率和时延来描述.假设υ_k是S为用户k传输的信息比特数，则传输速率为r_k=υ_k/M.

类似文献[5]的时延定义，假设从WET阶段到WIT阶段完成整个过程中，系统的时延为δ，δ^*为满足给定错误概率的最小时延.

$ \delta = (m + n){T_{\rm{c}}} $

(8)

对于系统总块错误概率，首先需要求出用户k错误概率的近似表达式.假设ε_k为用户k的错误概率，当n>100时，在准静态衰落信道中，ε_k可以近似为

$ {\varepsilon _k} \approx E\left\{ {Q\left[ {\frac{{C({\gamma _k}) - {r_k}}}{{\sqrt {V({\gamma _k})/M} }}} \right]} \right\} $

(9)

其中:C(γ_k)=lb(1+γ_k)为由香农定理得到的信道容量；V(γ_k)=[1－1/(1+γ_k)²](lbe)²为信道散度；$Q\left(x \right) = \int_x^\infty {(1/\sqrt {2{\rm{ \mathsf{ π} }}})} {{\rm{e}}^{ - {t^2}/2}}{\rm{d}}t$为高斯函数；E{·}为均值.由于高斯函数极为复杂，求解较难，所以从高斯函数的近似形式着手，可求得ε_k的闭合表达式.通过文献[10]得到高斯函数近似为

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} Q\left( {\frac{{C({\gamma _k}) - {r_k}}}{{\sqrt {V({\gamma _k})/M} }}} \right) \approx \varOmega ({\gamma _k}) = \\ \left\{ {\begin{array}{*{20}{l}} {1,\quad {\gamma _k} \le {\zeta _k}}\\ {\frac{1}{2} - {v_k}\sqrt M ({\gamma _k} - {\theta _k}),\quad {\xi _k} \le {\gamma _k} \le {\zeta _k}}\\ {0,\quad {\gamma _k} \ge {\xi _k}} \end{array}} \right. \end{array} $

(10)

其中:${\zeta _k} = {\theta _k} - 1/(2{\upsilon _k}\sqrt M)$，${\xi _k} = {\theta _k} + 1/(2{\upsilon _k}\sqrt M)$，${\upsilon _k} = 1/(2{\rm{ \mathsf{ π} }}\sqrt {{2^{2{r_k}}} - 1})$，${\theta _k} = {2^{{r_k}}} - 1$.

式(7)已经得出U_k的SNR，通过文献[11]得到用户k SNR的累积分布函数(CDF，cumulative distribution function)

$ {F_{{\gamma _k}}}(x) = 1 - \sum\limits_{i = 0}^{N - K} {\frac{{{{\left( {\frac{{\sigma _k^2}}{{{P_{\rm{S}}}{\beta _k}{\varOmega _k}}}x} \right)}^i}}}{{i!}}} {\rm{exp}}\left( { - \frac{{\sigma _k^2}}{{{P_S}{\beta _k}{\varOmega _k}}}x} \right) $

(11)

因此，其概率密度函数(PDF，probability density function)为

$ \begin{array}{*{20}{c}} {{f_{{\gamma _k}}}(x) = \sum\limits_{i = 0}^{N - K} {\frac{{{\psi ^{i + 1}}{x^i}}}{{i!}}} {\rm{exp}}( - \psi x) - }\\ {\sum\limits_{i = 1}^{N - K} {\frac{{{\psi ^i}{x^{i - 1}}}}{{(i - 1)!}}} {\rm{exp}}( - \psi x)} \end{array} $

(12)

其中$\psi = {\sigma _k}^2/{P_{\rm{S}}}{\beta _k}{\mathit{\Omega }_k}$.因此，式(9)可以表示为

$ {\varepsilon _k} = \int_0^\infty {{f_{{\gamma _k}}}} (x)Q\left( {\frac{{C(x) - {r_k}}}{{\sqrt {V(x)/M} }}} \right){\rm{d}}x $

(13)

把式(9)代入式(13)中，得

$ \begin{array}{*{20}{c}} {{\varepsilon _k} \approx \int_0^\infty {{f_{{\gamma _k}}}} (x)\varOmega ({\gamma _k}){\rm{d}}x = }\\ {\int_0^{{\zeta _k}} {{f_{{\gamma _k}}}} (x){\rm{d}}x + {\omega _1}\int_{{\zeta _k}}^{{\xi _k}} {{f_{{\gamma _k}}}} (x){\rm{d}}x - }\\ {{\omega _2}\int_{{\zeta _k}}^{{\xi _k}} x {f_{{\gamma _k}}}(x){\rm{d}}x} \end{array} $

(14)

其中：${\omega _1} = 1/2 + {\upsilon _k}\sqrt M {\theta _k}, {\omega _2} = {\upsilon _k}\sqrt M $.

为进一步推导，首先要计算$\smallint {x^n}{\rm{exp}}\left({ - \mu x} \right){\rm{d}}x$.通过文献[12]中式(3.351.1)知

$ \int_0^u {{x^n}} {\rm{exp}}( - \mu x){\rm{d}}x = {\mu ^{ - n - 1}}\Gamma (n + 1,\mu u) $

其中Γ(·, ·)为低阶不完全Gamma函数.因此

$ \begin{array}{*{20}{c}} {\int_{{u_1}} {{x^n}} {\rm{exp}}( - \mu x){\rm{d}}x = }\\ {\int_0^{{u_2}} {{x^n}} {\rm{exp}}( - \mu x){\rm{d}}x - \int_0^{{u_1}} {{x^n}} {\rm{exp}}( - \mu x){\rm{d}}x = }\\ {{\mu ^{ - n - 1}}[\Gamma (n + 1,\mu {u_2}) - \Gamma (n + 1,\mu {u_1})]} \end{array} $

(15)

然后把用户k SNR的PDF和以上关于∫xⁿexp(－μx)dx的计算代入式(14)中进行化简整理，得到用户k的块错误概率近似表达式，有

$ \begin{array}{l} \begin{array}{*{20}{c}} {{\varepsilon _k} \approx \sum\limits_{i = 0}^{N - K} {\frac{{\Gamma (i + 1,\psi {\zeta _k})}}{{i!}}} - \sum\limits_{i = 1}^{N - K} {\frac{{\Gamma (i,\psi {\zeta _k})}}{{(i - 1)!}}} + }\\ {{\omega _1}\sum\limits_{i = 0}^{N - K} {\frac{{\Gamma (i + 1,\psi {\xi _k}) - \Gamma (i + 1,\psi {\zeta _k})}}{{i!}}} - } \end{array}\\ \begin{array}{*{20}{c}} {{\omega _1}\sum\limits_{i = 1}^{N - K} {\frac{{\Gamma (i,\psi {\xi _k}) - \Gamma (i,\psi {\zeta _k})}}{{(i - 1)!}}} - }\\ {{\omega _2}{\psi ^{ - 1}}\sum\limits_{i = 0}^{N - K} {\frac{{\Gamma (i + 2,\psi {\xi _k}) - \Gamma (i + 2,\psi {\zeta _k})}}{{i!}} + } } \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\omega _2}{\psi ^{ - 1}}\sum\limits_{i = 1}^{N - K} {\frac{{\Gamma (i + 1,\psi {\xi _k}) - \Gamma (i + 1,\psi {\zeta _k})}}{{(i - 1)!}}} \end{array} $

(16)

2 总块错误概率最小化 2.1 问题公式化与分析

基于每个用户的块错误概率，利用文献[13]中基于BCD和CCCP的迭代算法，对WET阶段中的时隙数和WIT阶段中的功率分配系数因子进行联合优化，实现总块错误概率的最小化.优化问题可以表示为

$ {{\rm{(P1)}}:\mathop {{\rm{min}}}\limits_{m,{\beta _k} \ge 0,k = 1,2, \cdots ,K} {\varepsilon _0} = \sum\limits_{k = 1}^K {{\varepsilon _k}} } $

(17a)

$ {{\rm{ s}}{\rm{.t}}{\rm{. }}\sum\limits_{k = 1}^K {{\beta _k}} = 1} $

(17b)

$ {n \ge 100} $

(17c)

$ {0 < \delta \le {\delta _{{\rm{max}}}}} $

(17d)

式(17b)为功率分配系数的约束，式(17c)是由于式(9)仅在n大于100时才成立，式(17d)指定了最大的时延约束.由式(8)可以看出，时延δ直接受WET和WIT阶段的时隙数影响.因此，约束条件式(17d)可以写为0 < m+n≤s，s为WEI和WIT阶段最大的总时隙数.优化问题(P1)重新表示为

$ \begin{array}{*{20}{c}} {{\rm{(P2)}}:\mathop {{\rm{min}}}\limits_{m,{\beta _k} \ge 0,k = 1,2, \cdots ,K} {\varepsilon _0} = \sum\limits_{k = 1}^K {{\varepsilon _k}} }\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}\sum\limits_{k = 1}^K {{\beta _k}} = 1}\\ {n \ge 100}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0 < m + n \le s} \end{array} $

(18)

(P2) 是一个非凸问题，获得全局最优解时较为复杂.为了有效地解决问题(P2)，采用BCD算法得到局部最优解，再通过迭代算法获得全局最优解.首先，把问题(P2)分为2个子问题：

$ \begin{array}{l} \begin{array}{*{20}{c}} {({\rm{P2 - a}}):\mathop {\min }\limits_m {\varepsilon _0} = \sum\limits_{k = 1}^K {{\varepsilon _k}} }\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}n \ge 100}\\ {0 < m + n \le s} \end{array}\\ \begin{array}{*{20}{c}} {({\rm{P2 - b}}):\mathop {{\rm{min}}}\limits_{{\beta _k} \ge 0,k = 1,2, \cdots ,K} {\varepsilon _0} = \sum\limits_{k = 1}^K {{\varepsilon _k}} }\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}\sum\limits_{k = 1}^K {{\beta _k}} = 1} \end{array} \end{array} $

其中(P2-a)是给定β_k的m的优化子问题，需要在同时满足WIT阶段时隙数n和最大时延s的约束条件下，求得最优时隙数m；(P2-b)是给定m的β_k的优化子问题，在所有用户块长度都为M的条件下，进行合理的功率分配，使总块错误概率最小化.然后，迭代求解(P2-a)和(P2-b).

在第1次迭代中，初始化β_k=β_k^(t)，求解(P2-a)，并获得最优解m，作为m⁽¹⁾.然后，令m=m⁽¹⁾，求解(P2-b)并获得局部最优解β_k，作为β_k⁽¹⁾.

在第t(t>1)次迭代中，通过采用β_k=β_k^(t-1)，求解(P2-a)，并将最优m更新为m^(t)；然后，通过采用m=m^(t)，求解(P2-b)，并将β_k更新为β_k^(t).当ε₀收敛时，迭代停止.

2.2 (P2-a)的优化

通过式(5)和$\psi = {\sigma _k}^2/({P_{\rm{S}}}{\beta _k}{\mathit{\Omega }_k})$，得到ψ关于m的关系式，即$\psi = {\sigma _k}^2n/(\eta BP{\lambda _1}{\beta _k}{\mathit{\Omega }_k}m)$.对ε_k求m的二阶导函数，得

$ \frac{{{\partial ^2}{\varepsilon _k}}}{{\partial {m^2}}} = \frac{{{\partial ^2}{\varepsilon _k}}}{{\partial {\psi ^2}}}{\left( {\frac{{\partial \psi }}{{\partial m}}} \right)^2} + \frac{{\partial {\varepsilon _k}}}{{\partial \psi }}\frac{{{\partial ^2}\psi }}{{\partial {m^2}}} $

(19)

其中：$\frac{{\partial \psi }}{{\partial m}} = - \frac{{{\sigma ^2}_kn}}{{\eta BP{\lambda _1}{\beta _k}{\mathit{\Omega }_k}}}{m^{ - 2}}$，且${\left({\frac{{\partial \psi }}{{\partial m}}} \right)^2} > 0, \frac{{{\partial ^2}\psi }}{{\partial {m^2}}} = \frac{{2{\sigma _k}^2n}}{{\eta BP{\lambda _1}{\beta _k}{\mathit{\Omega }_k}}}{m^{ - 3}} > 0$.通过计算可知ω₂ξ_k=ω₁，ω₂ζ_k=ω₁－1，因此

$ \frac{{\partial {\varepsilon _k}}}{{\partial \psi }} = \frac{{{\omega _2}{\psi ^{ - 2}}[\Gamma (N - K + 2,\psi {\xi _k}) - \Gamma (N - K + 2,\psi {\zeta _k})]}}{{(N - K)!}} $

(20)

$ \begin{array}{*{20}{c}} {\frac{{{\partial ^2}{\varepsilon _k}}}{{{\partial ^2}\psi }} = - \frac{{\zeta _k^{N - K + 1}{\psi ^{N - K - 1}}{{\rm{e}}^{ - \psi {\zeta _k}}}}}{{(N - K)!}} + }\\ {\frac{{{\omega _2}\xi _k^{N - K + 1}{\psi ^{N - K - 1}}{{\rm{e}}^{ - \psi {\xi _k}}}}}{{(N - K)!}} + ( - 2)(N - K + 1) \times }\\ {{\omega _2}{\psi ^{ - 3}}[\Gamma (N - K + 2,\psi {\xi _k}) - \Gamma (N - K + 2,\psi {\zeta _k})]} \end{array} $

(21)

由Γ(N－K+2, ψξ_k)－Γ(N－K+2, ψζ_k)>0，ω₂, ψ, ξ_k, ζ_k>0，可得$\frac{{\partial {\varepsilon _k}}}{{\partial \psi }} > 0, \frac{{{\partial ^2}{\varepsilon _k}}}{{\partial {\psi ^2}}} < 0$.进而可知$\frac{{{\partial ^2}{\varepsilon _k}}}{{\partial {m^2}}} < 0$，即ε_k是关于m的一个凹函数，进而ε₀也是关于m的一个凹函数.约束条件式(17c)和式(18)是一个仿射集.因此，(P2-a)是一个凸优化问题，可以采用拉格朗日对偶算法对问题(P2-a)进行求解.首先，构造拉格朗日函数L(m, β_k, λ)=ε_k+λ(s－n－m)，其中λ为对偶因子.根据阶优化(KKT，karush-kuhn-tucker)条件^[14]，最优解需要满足

$ {\frac{{\partial L(m,{\beta _k},\lambda )}}{{\partial m}} = \frac{{\partial {\varepsilon _k}}}{{\partial m}} - \lambda = 0} $

(22a)

$ {\lambda (s - n - m) = 0} $

(22b)

其中

$ \begin{array}{*{20}{c}} {\frac{{\partial {\varepsilon _k}}}{{\partial m}} = \frac{{\partial {\varepsilon _k}\partial \psi }}{{\partial \psi }}\frac{{{\partial _2}\eta BP{\lambda _1}{\beta _k}{\varOmega _k}}}{{\sigma _k^2n}} \times }\\ {[\Gamma (N - K + 2,\psi {\xi _k}) - \Gamma (N - K + 2,\psi {\zeta _k})]}\\ {(N - K)!} \end{array} $

(23)

通过式(22a)得到$\lambda = \frac{{\partial {\varepsilon _k}}}{{\partial m}}.\frac{{\partial {\varepsilon _k}}}{{\partial m}}$是关于m的一个减函数，并且m>0.令g(m)=$\frac{{\partial {\varepsilon _k}}}{{\partial m}}$，因此g(m)也是一个减函数，即g(0)≥g(m)≥g(s－n).通过计算发现，g(m)≤0，但是λ≥0.因此，只有当g(m)=0时，式(22a)、式(22b)成立.但是，通过g(m)≤0可知$\frac{{\partial {\varepsilon _k}}}{{\partial m}}$≤0，则ε_k是关于m的一个减函数.又因为0 < m≤s－n，所以，最优解为

$ {m^*} = s - n $

(24)

2.3 (P2-b)的优化

(P2-b)是一个非凸优化问题，通过CCCP算法^[15]，可以近似迭代的方式解决. CCCP算法的基本思想是通过迭代的方式将原问题近似为具有一阶泰勒级数展开的凸问题，并求解近似凸问题，直到收敛.

首先，令β₀为初始值，β_i=[β_{i, 1}, β_{i, 2}, …, β_{i, K}]^T(i≥1)为第i次迭代的最优值；然后，在第i次迭代时，ε_k的一阶泰勒级数展开可以表示为

$ \begin{array}{l} {\varepsilon _k}({\beta _{i,k}}) = {\varepsilon _k}({\beta _{i - 1,k}}) + \nabla [{\varepsilon _k}({\beta _{i - 1,k}})]({\beta _{i,k}} - {\beta _{i - 1,k}}) + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}{\nabla ^2}[{\varepsilon _k}({\beta _{f,k}})]{({\beta _{i,k}} - {\beta _{i - 1,k}})^2} \le \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\varepsilon _k}({\beta _{i - 1,k}}) + \nabla [{\varepsilon _k}({\beta _{i - 1,k}})]({\beta _{i,k}} - {\beta _{i - 1,k}}) \buildrel \Delta \over = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\varepsilon _k}({\beta _{i,k}};{\beta _{i - 1,k}}) \end{array} $

(25)

其中：β_{f, k}表示β_{i－1, k}和β_{i, k}之间的一个数；$\nabla $[ε_k(β_{i－1, k})]为β_k=β_{i－1, k}时ε_k(β_k)的梯度；$\nabla [{\varepsilon _k}({\beta _{i -1, k}})] = \frac{{\partial {\varepsilon _k}({\beta _k})}}{{\partial {\beta _k}}}\left| {_{{\beta _k} = {\beta _{i - 1, k}}}} \right.$；${\nabla ^2}[{\varepsilon _k}({\beta _{f, k}})]$为β_k=β_{f, k}时ε_k(β_k)的二阶偏导数，${\nabla ^2}[{\varepsilon _k}({\beta _{f, k}})]$≤0.通过式(25)得到ε₀的上界：

$ {\varepsilon _0}(m,\mathit{\boldsymbol{\beta }}) \le \sum\limits_{k = 1}^K {{\varepsilon _k}} ({\beta _{i,k}};{\beta _{i - 1,k}}) \buildrel \Delta \over = {\varepsilon _{hb}}(m,\mathit{\boldsymbol{\beta }};{\mathit{\boldsymbol{\beta }}_{i - 1}}) $

(26)

进而，通过优化β，使上界函数ε_hb(m, β; β_i－1)最小化，代替了最小化ε₀(m, β).最小化上界函数ε_hb(m, β; β_i－1)可以表示为

$ ({\rm{P2}} - {\rm{c}}):\mathop {{\rm{min}}}\limits_{{\beta _k} \ge 0,k = 1,2, \cdots ,K} {\varepsilon _{hb}}(m,\mathit{\boldsymbol{\beta }};{\mathit{\boldsymbol{\beta }}_{i - 1}}){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{s}}{\rm{.t}}{\rm{.}} $

ε₀(m, β; β_i－1)也是一个凹函数，同时，约束条件(17b)是一个仿射集.因此，(P2-c)是一个凸优化问题.

在使系统的总块错误概率最小化时，需要在每个用户满足各自可靠性需求的前提下，对功率分配系数进行优化. γ_min^k为用户k所需的最小SNR，因此，每个用户的SNR需要满足

$ {\gamma _k} \ge \gamma _{{\rm{min}}}^k,k = 1,2, \cdots ,K $

(27)

根据式(7)，上面的约束条件可以转换为

$ {\beta _k} \ge \frac{{\gamma _{\min }^k\sigma _k^2}}{{{P_{\rm{S}}}|h_k^\dagger \prod\nolimits_k {{\mathit{\boldsymbol{h}}_k}} |}} $

(28)

因此，优化问题(P2-c)可以重新表示为

$ ({\rm{P2}} - {\rm{d}}):\mathop {{\rm{min}}}\limits_{{\beta _k} \ge 0,k = 1,2, \cdots ,K} {\varepsilon _0}(m,\mathit{\boldsymbol{\beta }};{\mathit{\boldsymbol{\beta }}_{i - 1}}){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{s}}{\rm{.t}}{\rm{.}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({\rm{17b}}),({\rm{28}}) $

ε₀(m, β; β_i－1)的最小化即等同于ε_k(β_{i, k}; β_{i－1, k})的最小化.

通过计算得$\frac{{\partial {\varepsilon _k}}}{{\partial {\beta _k}}} < 0$，因此，ε_k是关于β_k的一个减函数. ε_k(β_{i, k}; β_{i－1, k})与ε_k具有相同的单调性，ε_k(β_{i, k}; β_{i－1, k})也是关于β_k的一个减函数.为此，ε_k(β_{i, k}; β_{i－1, k})的最小化就等同于寻找β_k的最小值.

通过拉格朗日对偶算法对问题(P2-d)进行求解.首先，构造拉格朗日函数L(m, β_k, λ)=ε_k+$\alpha \left({{\beta _k} - \frac{{{\gamma _{{\rm{min}}}}^k{\sigma _k}^2}}{{{P_{\rm{S}}}\left| {{\mathit{\boldsymbol{h}}_k}^{†} \mathop \prod \limits_k {\mathit{\boldsymbol{h}}_k}} \right|}}} \right) + {\rm{ }}\chi \left({1 - \sum\limits_{k = 1}^K {{\beta _k}} } \right)$，其中α、χ为拉格朗日乘子.根据KKT条件，最优解需要满足

$ {\frac{{\partial {\rm{L}}(m,{\beta _k},\lambda )}}{{\partial {\beta _k}}} = \frac{{\partial {\varepsilon _k}}}{{\partial {\beta _k}}} + \alpha - \chi = 0} $

(29)

$ {\alpha \left( {{\beta _k} - \frac{{\gamma _{{\rm{min}}}^k\sigma _k^2}}{{{P_{\rm{S}}}|\mathit{\boldsymbol{h}}_k^\dagger \prod\limits_k {{\mathit{\boldsymbol{h}}_k}} |}}} \right) = 0} $

(30)

$ {\alpha \ge 0} $

(31)

$ {\chi \left( {1 - \sum\limits_{k = 1}^K {{\beta _k}} } \right) = 0} $

(32)

$ {\chi \ge 0} $

(33)

类似于式(29)、式(30)和式(32)的推导，可以证明，为了导出问题(P2-d)的最优解，约束式(28)应该是等式.基于此，可以得到β_k的闭型表达式为

$ {{\beta _k} = \frac{{\gamma _{{\rm{min}}}^k\sigma _k^2}}{{{P_{\rm{S}}}|\mathit{\boldsymbol{h}}_k^\dagger \prod\limits_k {{\mathit{\boldsymbol{h}}_k}} |}}} $

(34)

最后问题(P2-d)的最优解为

$ {\beta _k} = \left\{ {\begin{array}{*{20}{l}} {\frac{{\gamma _{{\rm{min}}}^k\sigma _k^2}}{{{P_{\rm{S}}}\left| {\mathit{\boldsymbol{h}}_k^\dagger \prod\limits_k {{\mathit{\boldsymbol{h}}_k}} } \right|}},\quad 1 \le k \le K - 1}\\ {1 - \sum\limits_{i = 1}^{K - 1} {{\beta _i}} ,\quad k = K} \end{array}} \right. $

(35)

2.4 联合优化的迭代算法

算法1为解决问题(P1)的具体过程，步骤5~7是解决问题(P2-b)的迭代过程，步骤2~9是解决问题(P2-a)的迭代过程.

算法1  优化问题(P1)的迭代算法

1  初始化β=β⁽⁰⁾，t=1.

2  repeat

3  设置β=β^(t－1)，利用式(24)计算m^(t).

4  设置m=m⁽⁰⁾，β₀=β^(t－1)，i=1.

5  repeat

6  给定β_i，利用式(33)计算β_i+1，并且i←i+1.

7  until β收敛.

8  β^(t)=β_i^*，并且t←t+1.

9   until

$\mathit{\Gamma } = |{\varepsilon _0}({m^{(t)}}, {\mathit{\boldsymbol{\beta }}^{(t)}}) - {\varepsilon _0}({m^{(t - 1)}}, {\mathit{\boldsymbol{\beta }}^{(t - 1)}})|/{\varepsilon _0}({m^{(t)}}, {\mathit{\boldsymbol{\beta }}^{(t)}})$收敛.

为了保证算法1的收敛性，在这一部分需要证明ε₀(m^(t+1), β^(t+1))≤ε₀(m^(t), β^(t)).在BCD算法的每次迭代中都存在CCCP算法，因此将首先证明CCCP算法在BCD算法第t次迭代中的收敛.令β_i+1为CCCP算法第i+1次迭代时，问题(P2-b)的最优解，可以得到

$ \begin{array}{*{20}{c}} {{\varepsilon _{hb}}({m^{(t)}},{\mathit{\boldsymbol{\beta }}_{i + 1}};{\mathit{\boldsymbol{\beta }}_i}) = }\\ {\left\{ {\mathop {{\rm{min}}}\limits_{{\beta _k} \ge 0,k = 1,2, \cdots ,K} {\varepsilon _{hb}}({m^{(t)}},\mathit{\boldsymbol{\beta }};{\mathit{\boldsymbol{\beta }}_{i - 1}})} \right\} \le }\\ {{\varepsilon _{hb}}({m^{(t)}},{\mathit{\boldsymbol{\beta }}_i};{\mathit{\boldsymbol{\beta }}_i})} \end{array} $

(36)

由于ε_hb(m, β; β_i－1)是ε₀(m, β)的上界，所以

$ {\varepsilon _{hb}}({m^{(t)}},{\mathit{\boldsymbol{\beta }}_i};{\mathit{\boldsymbol{\beta }}_i}) = {\varepsilon _0}({m^{(t)}},{\mathit{\boldsymbol{\beta }}_i}) \le {\varepsilon _{hb}}({m^{(t)}},{\mathit{\boldsymbol{\beta }}_i};{\mathit{\boldsymbol{\beta }}_{i - 1}}) $

(37)

通过式(24)和式(25)可以得到

$ {\varepsilon _{hb}}({m^{(t)}},{\mathit{\boldsymbol{\beta }}_{i + 1}};{\mathit{\boldsymbol{\beta }}_i}) \le {\varepsilon _{hb}}({m^{(t)}},{\mathit{\boldsymbol{\beta }}_i};{\mathit{\boldsymbol{\beta }}_{i - 1}}) $

(38)

保证了在BCD算法的每次迭代中CCCP算法的收敛.

接下来证明BCD算法的收敛性.令m_i+1为CCCP算法第i+1次迭代时，问题(P2-a)的最优解.因此，可以得到

$ \begin{array}{*{20}{c}} {{\varepsilon _0}({m^{(t + 1)}},{\mathit{\boldsymbol{\beta }}^{(t)}}) = \{ \mathop {{\rm{min}}}\limits_m {\varepsilon _0}(m,{\mathit{\boldsymbol{\beta }}^{(t)}})\} \le }\\ {{\varepsilon _0}({m^{(t)}},{\mathit{\boldsymbol{\beta }}^{(t)}})} \end{array} $

(39)

通过式(24)~式(27)可以得到ε₀(m^(t+1), β^(t+1))≤ε₀(m^(t), β^(t))，即证明了算法1的收敛性.

3 仿真结果

为了描述路径衰落的影响，信道衰落模型为Ω_k=10^－3/(1+d_k^ω)，其中Ω_k为源节点S到用户k之间的平均信道功率增益，d_k为源结点S到用户k之间的距离，ω为路径损耗系数.假设ET-S的距离设定为12 m，带宽为1 MHz，每个时隙的持续时间T_c=1 μs.在WET阶段，能量发射器ET的发射功率P_ET=30 dBm.能量捕获效率η=0.5，路径衰落因子ω₁=2.在WIT阶段，总用户数K=3，路径损耗系数ω₂=3、噪声功率σ_k²=－80 dBm，源节点S与每个用户传输的信息量分别为υ₁=600，υ₂=800，υ₃=1 000，传输距离d₁=100，d₂=80，d₃=130.

图 2所示为源节点S的传输功率P_S与WET阶段的时隙数m的关系.可以看出，P_S除了受时隙数m的影响外，还会受到发射天线N_ET和接收天数N的影响.在相同时隙下，传输功率随着天线数量的增加而增大，同时也发现，当发射天线N_ET和接收天线N的数值进行互换时，两条区间几乎重合，说明RF信号收发端的天线数对P_S具有相同的影响.

图 3所示为不同传输功率P_S所对应的块错误概率ε_k.其中，功率分配系数β_k=1/3, ∀k，传输距离d_k=100，传输的信息量υ=600 bit.可以发现，随着源节点传输功率P_S的增大，用户k的块错误概率ε_k越小.同时，发现源节点S的天线数量N和数据包的长度M都会影响到ε_k.当传输功率P_S、数据包长度M一定时，天线数量N越多，块错误概率ε_k越小；当传输功率P_S、天线数量N一定时，数据包长度M越长，块错误概率ε_k越小.因此，可以通过增加传输功率P_S、天线数量N和数据包长度M来减少块错误概率ε_k，提高系统的可靠性.

图 4所示为用户k的错误概率ε_k与WIT时隙数n的关系.从仿真图中发现，刚开始随着时隙数n的增大，用户k的错误概率ε_k急剧减少，但是当n≥1 000时，随着n的增大，ε_k的变化速率越来越缓慢.从仿真图中还发现，源节点S的天线数相比ET的电线数对用户的块错误概率影响较大，在ET天线数一定时，增加S端的天线数量，用户的块错误概率显著减小.

从图 4中发现，WIT时隙数n为1 000时，信息传输的可靠性较好.因此，令n=1 000，其他参数保持不变. 图 5所示为总块错误概率ε₀与WEI和WIT阶段总时隙数s的关系.从图 5中可以明确的观察到，随着总时隙数s的增加，总块错误概率不断减小，这意味着系统的时延越大，错误概率越小，信息传输的可靠性越高.在时延一定时，可以通过增加WIT阶段的时隙数，来增加数据包的长度，降低用户的块错误概率，进而提高系统的可靠性.可以根据实际情况的需求，使系统的时延和可靠性进行权衡，得到最佳状态.同时，通过仿真图发现，在对系统整体进行分析时，源节点的天线数量对系统的可靠性影响比较大，这与图 2、图 4对天线数量的分析结果一致.因此，可以通过增加源节点的天线数量来减少系统的块错误概率，提高系统的可靠性.

4 结束语

研究了短包通信场景下多用户无线能量传输系统的块错误概率.推导出每个用户块错误率的近似闭合形式表达式，并验证了其准确性.根据闭合表达式，分析影响块错误概率的因素，得到可以通过优化发射天线数量、数据包长度以及发射功率提高用户的块错误概率，同时也可以优化系统的功率分配系数，提高系统的块错误概率.此外，通过仿真验证了块错误概率和时延的关系.结果表明，可靠性越严格，需要WET的时隙数越多，进而延迟越高.同时发现，当时延一定时，可以通过增加数据包长度减少系统的块错误概率，为系统的时延和可靠性的权衡提供了理论依据.