1 Introduction↩︎

With the rapid advancement of the Internet of Things (IoT), the number of connected wireless devices is experiencing exponential growth and is projected to reach an estimated 5 trillion by 2030 [1]. The realization of such a large-scale network necessitates the widespread deployment of IoT devices, which in turn imposes significant and unprecedented demands on energy consumption and spectrum utilization. However, the limited available spectrum cannot meet the growing demand from massive IoT deployments, particularly when each device requires a dedicated frequency band [2]. Moreover, the widespread deployment of radio frequency (RF) components in IoT devices results in high energy consumption, thereby increasing operational costs. These challenges present major obstacles to the development of next-generation communication systems, highlighting the urgent need for innovative solutions that prioritize spectrum and energy efficiency [3].

Recently, symbiotic radio (SR) has attracted increasing research interest due to its ability to enable spectrum- and energy-efficient IoT communications [4], offering a promising approach to address the aforementioned challenges. In SR, the passive secondary transmission, also referred to as IoT transmission [5], is achieved by the backscatter device (BD) acting, which operates as an IoT node by embedding its own information into the primary signal to IoT receiver without generating an RF carrier. Moreover, the primary transmission can benefit from the secondary transmission by regarding the backscatter link as a multipath component. Hence, they establish a mutually beneficial symbiotic relationship[6]. Driven by these advantages, substantial research efforts have been devoted to SR[4], [7], [8]. Based on the relationship between the symbol durations of the primary transmitter (PT) and the BD, SR can be categorized into commensal SR (CSR) and parasitic SR (PSR). In PSR, the BD and PT share the same symbol duration, whereas in CSR, the BD’s symbol duration is significantly longer than that of the PT [4], [7]. To enhance spectral efficiency, transmit power minimization has been investigated for full-duplex SR systems [8]. Although BD can provide an additional reflection path, the assistance from the reflection path alone is very weak, when the line-of-sight (LoS) path experiences severe fading or blockage. Improving wireless environments has become a key research focus, aiming to create favorable channel conditions to enhance communication performance[9].

1.1 Related Work↩︎

1.1.1 Reconfigurable Antenna aided SR↩︎

To enable efficient backscatter communication, [10] investigated the roles of reconfigurable intelligent surfaces (RIS) for assisting different scenarios of backscatter communication. By configuring RIS to act as a BD, backscatter communication can benefit from the additional spatial degrees of freedom introduced by multiple reflecting elements. Furthermore, the problem of minimizing total transmit power in RIS-assisted multiple-input multiple-output (MIMO) SR systems has been analyzed in [11]. In [12], RIS-enabled SR with orthogonal frequency division multiplexing (OFDM) transmission has been investigated under imperfect symbol synchronization to improve spectral efficiency. More recently, movable antennas (MAs) [13], also known as fluid antennas [14], have emerged as a promising technology for intelligently reconfiguring wireless environments by flexibly adjusting antenna positions. The utilization of MAs in PT can effectively improve the rate of secondary transmission by optimizing the positions of MAs to strengthen the beamforming gain at the BD [15], [16]. However, the performance improvements offered by these technologies are often limited due to severe path loss, particularly in high-frequency bands [17]. For instance, the double fading effect inherent in the cascade reflection link of RIS leads to significantly higher path loss compared to a direct LoS link [18]. Similarly, the movement range of MAs is typically limited to only a few wavelengths, which restricts their overall performance gains.

1.1.2 Pinching-Antenna Systems↩︎

To overcome these limitations, the pinching antenna system (PASS), recently proposed by DOCOMO [19], has emerged as a promising solution in the domain of flexible-antenna technologies [19]–[21]. PASS leverages a dielectric waveguide as its transmission medium to establish adjustable LoS links with users. The system enables the signal radiation from any desired radiation points that are activated by implementing dielectric particles [22]. These dielectric particles are referred to as pinching antennas (PAs), which exhibit properties similar to those of leaky-wave antennas [23]. However, in contrast to the leaky-wave based systems where the antennas are fixed in place with pre-defined locations, PAs support flexible and dynamic activation, allowing signals radiated from dielectric waveguides to adapt effectively to complex and time-varying environments. This capability enables a cost-efficient and scalable MIMO implementation through the novel concept of pinching beamforming [24], which enhances communication performance by dynamically optimizing antenna configurations.

Driven by the above promising characteristics, PASS has attracted increasing research attention, although it remains in the early stages of development. In [20], both single-waveguide and multi-waveguide scenarios were investigated, and low-complexity pinching beamforming schemes were proposed for single-user and two-user multi-input single-output (MISO) systems. In [25], the authors addressed a joint transmit and pinching beamforming optimization problem for a multi-user PASS downlink framework, introducing both an optimization-based majorization-minimization and penalty dual decomposition method and a learning-based knowledge-guided dual learning approach. Furthermore, [26] developed two efficient deep learning-driven channel estimation methods for PASS, demonstrating their superior estimated performance and low pilot overhead. The achievable array gain of PASS was analyzed in [27], [28], where [27] proposed an antenna position enhancement algorithm to approximate its performance upper bound, and [28] demonstrated that the presence of LoS blockage can enhance the performance gain of pinching antennas compared to conventional antennas. Additionally, energy-efficient resource allocations for PASS were studied in [29]and [30]. Moreover, [31] presented a comprehensive analytical framework for evaluating the performance of PASS, including closed-form expressions for the average achievable rate and outage probability.

1.2 Motivations and Contributions↩︎

Based on the above discussion, PASS has demonstrated strong capabilities in establishing robust LoS links, significantly reducing free-space propagation loss, and overcoming blockage issues. Consequently, the utilization of PAs in the SR systems is essential for achieving highly reliable and spectrally efficient primary and secondary transmissions. Moreover, different from conventional antenna systems, the flexible deployment of PAs introduces additional degrees of freedom, facilitating effective pinching beamforming tailored to the locations of the BD, the IR, and primary receivers (PR). To the best of our knowledge, the application of PASS in SR systems remains largely unexplored in the existing literature. Motivated by this gap, this paper proposes a PASS-enabled downlink SR framework and develops joint beamforming methods. The main contributions of this work are summarized as follows.

• We propose a novel PASS-assisted downlink SR framework, where the PASS BS with multiple waveguides treated as a PT serves an IR and PRs with the assistance of the BD. Within this model, a joint transmit and pinching beamforming optimization problem is formulated for maximizing the sum rate, while satisfying the constraints of the detection error probability of the IR and the feasible deployment region of PAs. To tackle this highly coupled nonconvex problem, we develop both learning-aided gradient descent (LGD) and two-stage optimization-based algorithms.

• For the LGD algorithm, we address the constraints by equivalently transforming the constrained optimization problem into a tractable form that can be directly solved using gradient descent. Specifically, the position matrix of PAs is reparameterized as non-negative offsets to satisfy the minimum spacing constraint, while the maximum transmit power and waveguide length constraints are handled via projection and normalization techniques. Conventional manual gradient derivation or symbolic differentiation often leads to expression swelling and computational inefficiency. In contrast, the proposed LGD framework leverages automatic differentiation and the Adam optimizer, allowing efficient updates of optimization variables modeled as learnable parameters updated by back-propagation.

• For the two-stage optimization-based approach, a successive convex approximation and particle swarm optimization (SCA-PSO) algorithm is exploited. The original joint optimization problem is decoupled into two subproblems. In the first stage, the subproblem with respect to the transmit beamforming matrix, which is inherently non-convex, is approximated via SCA by transforming the objective and constraint functions into concave forms, enabling efficient solution through convex optimization tools. In the second stage, the pinching beamforming design is addressed using a PSO-based algorithm, where each particle encodes a candidate PA deployment matrix, and its fitness is evaluated based on the resulting achievable sum rate.

• Finally, numerical simulations are conducted to evaluate the effectiveness of the proposed framework and algorithms. The results demonstrate that: i) The proposed PASS-enabled SR system achieves significantly higher sum rates compared to conventional antenna-based SR schemes. ii) The proposed SCA-PSO algorithm improves performance by 17.1% and 35.5% over the low-complexity LGD and fixed-PA schemes, respectively, and approaches the performance of the element-wise optimization method.

1.3 Organization and Notations↩︎

The structure of the paper is as follows. Section II describes the PASS-assisted SR system and formulates the sum-rate maximization problem. Section III introduces the equivalent reformulation of the original problem and proposes a GD-based joint beamforming framework. Section IV presents the proposed SCA-PSO algorithm for joint transmit and pinching beamforming, designed to further enhance system performance. Section V provides numerical results that validate the convergence behavior and demonstrate the performance advantages of the proposed framework and algorithms. The concluding remarks are provided in Section VI.

Notation: Scalars, vectors, and matrices are represented by \(x\), \(\mathbf{x}\), and \(\boldsymbol{X}\), respectively. \({\left( \cdot \right)^T}\), \({\left( \cdot \right)^*}\), and \({\left( \cdot \right)^H}\) stand for the transpose, complex conjugate, and conjugate transpose operations, respectively. The notation \({\mathop{\rm Re}\nolimits} \left\{ \cdot \right\}\) and \({\mathop{\rm Im}\nolimits} \left\{ \cdot \right\}\) denote the real and imaginary part of a complex number, respectively. \({\rm{Tr}}\left( \cdot \right)\), \(\left| \cdot \right|\), and \(\left\| \cdot \right\|\) represent the trace, the modulus operator, and the Euclidean norm, respectively. The blkdiag \((a_1,...,a_N)\) is a block diagonal matrix with diagonal blocks \(a_1,...,a_N\). \({C^{M \times N}}\) denotes the dimension of an \(M \times N\) complex-valued matrix. \({\cal C}{\cal N}\left( {\mu ,\sigma^2} \right)\) is the circularly symmetric complex Gaussian random distribution with mean \(\mu\) and variance \(\sigma^2\).

2 System Model↩︎

As illustrated in Fig. 1, this paper considers a downlink PASS-assisted SR system, where the PASS is connected to a base station (BS) to simultaneously serve \(K\) single-antenna PRs and a single-antenna IoT receiver (IR) via a BD. The PASS consists of \(N\) dielectric waveguides, with \(N \ge K\), where each waveguide is incorporated with \(M\) pinching antennas. Activating PAs at different positions along the waveguides enables flexible control over the phases of the incident signals and the large-scale fading. The waveguides in PASS are fed by the BS, which transmits primary signals to PRs with BD-assisted backscattering. This process not only facilitates primary transmission but also embeds its messages within the primary signals intended for the IR.

2.1 Channel Model↩︎

Assume that both PAs and waveguides are located at a fixed height of \(z^{\rm{PA}}\), and a three-dimensional Cartesian coordinate system is established. The location of the \(m\)-th PA associated with the \(n\)-th waveguide is denoted by \({\mathbf{I}}_{n,m}^{{\rm{PA}}} = {[{x_{n,m}},{y_n},{z^{{\rm{PA}}}}]^T}\), where \({x_{n,m}}\) is the adjustable coordinate over the \(x\)-axis, and \(y_n\) the fixed and pre-defined coordinate over \(y\)-axis, and \({\mathbf{I}}_k^{\rm{U}} = {[x_k^{\rm{U}},y_k^{\rm{U}},0]^T}\) represent the \(k\)-th PR’s position. Given that non-line-of-sight (NLoS) paths are significantly weaker than line-of-sight (LoS) paths, we adopt a practical channel model that considers only the LoS components while ignoring the NLoS components [32], [33]. Based on the geometric free-space spherical model, the channel from the PA \({\mathbf{ I} _{n,m}^{{\rm{PA}}}}\) to the \(k\)-th PR at the \(\mathbf{I} _k^{\rm{U}}\) is given by

\[h_{k,n,m}^H({x_{n,m}}) = \frac{{\kappa {e^{ - j{\beta _h}r({x_{n,m}},{\mathbf{I}}_k^{\rm{U}})}}}}{{r({x_{n,m}},{\mathbf{I}}_k^{\rm{U}})}}, \label{eq1}\tag{1}\] where \(\kappa = c/4\pi {f_c}\) denotes the reference channel gain at a distance of 1 m, with \(c\) and \(f_c\) representing the speed of light and carrier frequency, respectively. \({\beta _h} = {{2\pi } \mathord{\left/ {\vphantom {{2\pi } {{\lambda _f}}}} \right. \kern-\nulldelimiterspace} {{\lambda _f}}}\) is the wave number in the propagation medium, and \({\lambda _f}\) is the corresponding wavelength. \(r({x_{n,m}},{\mathbf{I}}_k^{\rm{U}}) = \left\| {{\mathbf{I}}_k^{\rm{U}} - {\mathbf{I}}_{n,m}^{{\rm{PA}}}} \right\|\) represents the distance from the PA \({\mathbf{I} _{n,m}^{{\rm{PA}}}}\) to \({{\mathbf{I}}_k^{\rm{U}}}\). Similarly, the channels from the PA \({\mathbf{I} _{n,m}^{{\rm{PA}}}}\) to the IR and the BD are denoted by \(h_{{\text{\text{IR}}},n,m}^H\) and \(h_{\text{\text{BD}},n,m}^H\), respectively. Stacking the channel vectors from all the PAs to the \(k\)-th PR, \(\mathbf{h}_k^H(\boldsymbol{X}) = [\mathbf{h}_{k,1}^H({x_1}),...,\mathbf{h}_{k,N}^H({x_N})] \in {\mathbb{C}^{1 \times NM}}\) is the overall channel vector for the \(k\)-th PR. Furthermore, for in-waveguide transmission, we denote the diagonal matrix \(\boldsymbol{G}(\boldsymbol{X}) \in {\mathbb{C}^{NM \times N}}\) as representing the path response from the feed point of each waveguide to the corresponding PAs, which is given by \[\begin{align} \nonumber {\boldsymbol{G}}({\boldsymbol{X}})& = {\text{blkdiag}}({\boldsymbol{g}}\left( {{x_1}} \right),...,{\boldsymbol{g}}\left( {{x_N}} \right)) \\ \nonumber &= \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{g}}\left( {{x_1}} \right)}&0& \cdots &0\\ 0&{{\boldsymbol{g}}\left( {{x_2}} \right)}& \cdots &0\\ \vdots & \vdots & \ddots & \vdots \\ 0&0& \cdots &{{\boldsymbol{g}}\left( {{x_N}} \right)} \end{array}} \right]. \end{align}\] The response vector \(\mathbf{g}\left( {{x_1}} \right)\) is characterized by \[\mathbf{g}\left( {{x_n}} \right) = [{\upsilon _1}{e^{ - j{\beta _g}{x_{n,m}}}},...,{\upsilon _M}{e^{ - j{\beta _g}{x_{n,M}}}}], \label{eq2}\tag{2}\] where \({\upsilon _1}\) is the amplitude of the transmitted signal and \({\beta _g} = {{2\pi {n_{\text{eff}}}} \mathord{\left/ {\vphantom {{2\pi n_{\text{eff}}} {{\lambda _f}}}} \right. \kern-\nulldelimiterspace} {{\lambda _f}}}\) denotes the propagation constant of the waveguide, with \({{n_{\text{eff}}}}\) being the effective refractive index of the dielectric waveguide.

2.2 Signal Model↩︎

Denote the symbol vector transmitted from the BS to the PRs by \(\mathbf{s}(l) \in {\mathbb{C}^{K \times 1}}\) with \(\mathbb{E}[\mathbf{s}(l){\mathbf{s}^H}(l)] = {\boldsymbol{I}_K}\). Let \(\boldsymbol{W} = [{\mathbf{w}_1},...,{\mathbf{w}_K}] \in {\mathbb{C}^{N \times K}}\) represent the transmit beamforming matrix, with the total transmit power satisfying \(\operatorname{Tr}({\boldsymbol{W}}{{\boldsymbol{W}}^H}) \le {P_{\max }}\). Subsequently, the transmitted signal at the BS is given by \({\boldsymbol{W}}\mathbf{s}\left( l \right)\). Additionally, the BD transmits its own signal \(c\) to the IR by using the ON-OFF keying (OOK) modulation, i.e., \(c =\) “0” and “1” correspond to OFF and ON states, respectively. At the \(l\)-th time slot, the received signal at the \(k\)-th PR is given by

\[{y_k}(l) = \left( {\mathbf{h}_k^{eq} + c{f_{\text{BD},k}}\mathbf{h}_{\text{BD}}^{eq}} \right){\boldsymbol{W}}\mathbf{s}(l) + {n_k}(l), \label{eq3}\tag{3}\] where \(\mathbf{h}_k^{eq} = \mathbf{h}_k^H(\boldsymbol{X})\boldsymbol{G}(\boldsymbol{X})\) and \(\mathbf{h}_{\text{BD}}^{eq} = \mathbf{h}_{\text{BD}}^H(\boldsymbol{X})\boldsymbol{G}(\boldsymbol{X})\) are the equivalent channels from the BS to the \(k\)-th PR and BD, respectively. \({{f_{\text{BD},k}}}\) and \({n_k}(l) \sim {\cal C}{\cal N}(0,\delta _k^2)\) denote the channel reflective-link from the BD to the \(k\)-th PR and the additive white Gaussian noise, while \({{{P_{\max }}} \mathord{\left/ {\vphantom {{{P_{\max }}} {\delta _k^2}}} \right. \kern-\nulldelimiterspace} {\delta _k^2}}\) is the transmit signal-to-noise ratio (SNR) of the \(k\)-th user [20]. Since the communication rate of \(c\) is much lower than that of \(\mathbf{s}(l)\), we assume \(T_c = L T_s\), \(L \gg 1\), where \(T_c\) and \(T_s\) denote the symbol period of \(c\) and \(\mathbf{s}(l)\). When decoding \(\mathbf{s}(l)\), the backscatter link formed by the BD can be treated as an additional path. Since the PRs have no prior knowledge about the BD’s symbol \(c\), we assume non-coherent detection can be applied to detect \(\mathbf{s}(l)\) with partial CSI [34]. Therefore, the signal-to-interference-plus-noise ratio (SINR) at the \(k\)-th PR can be expressed as

\[{\text{SINR}_k} = {\frac{{{{\left| {\left( {\mathbf{h}_k^{eq} + c{\mathbf{f}_{b,k}}} \right){\mathbf{w}_k}} \right|}^2}}}{{\sum {_{i = 1,i \ne k}^K{{\left| {\left( {\mathbf{h}_k^{eq} + c{\mathbf{f}_{b,k}}} \right){\mathbf{w}_i}} \right|}^2} + \delta _k^2} }}} . \label{eq4}\tag{4}\] where \({\mathbf{f}_{b,k}} = {f_{\text{BD},k}}\mathbf{h}_{\text{BD}}^{eq}\) denote the backscatter cascade channel from the BS to the \(k\)-th PR via the BD. Assuming that the BS transmits the symbols ‘0’ and ‘1’ with equal probability, the average achievable rate of decoding \(\mathbf{s}(l)\) for the \(k\)-th PR is given by [11]

\[\begin{align} \nonumber {R_k}& = {\mathbb{E}_c}\left[\log(1 + {\text{SINR}_k}) \right] \\ \nonumber &= \frac{1}{2}{\log _2}(1 + \frac{{{{\left| {\mathbf{h}_k^{eq}{\mathbf{w}_k}} \right|}^2}}}{{\sum {_{i = 1,i \ne k}^K{{\left| {\mathbf{h}_k^{eq}{\mathbf{w}_i}} \right|}^2} + \delta _k^2} }}) \\ &\quad + \frac{1}{2}{\log _2}(\frac{{{{\left| {(\mathbf{h}_k^{eq} + {\mathbf{f}_{b,k}}){\mathbf{w}_k}} \right|}^2}}}{{\sum {_{i = 1,i \ne k}^K{{\left| {(\mathbf{h}_k^{eq} + {\mathbf{f}_{b,k}}){\mathbf{w}_i}} \right|}^2} + \delta _k^2} }}). \end{align}\]

At the \(l\)-th time slot within one BD symbol period, the received signal at the IR is given by

\[\label{eq6} {y_{\text{IR}}}(l) = \left( {h_{\text{IR}}^{eq} + c{f_{\text{b},\text{IR}}}} \right){\boldsymbol{W}}\mathbf{s}(l) + {n_{\text{IR}}}(l),\tag{5}\] where \({\mathbf{f}_{\text{b},\text{IR}}} = {f_{\text{BD},{\text{IR}}}}\mathbf{h}_{\text{BD}}^{eq}\). The IR aims to recover the symbol \(c\) from the received signal by distinguishing between two hypotheses corresponding to the BD’s transmitted symbol, either “0” or “1”. Before detecting, the IR first decodes \(\mathbf{s}(l)\) and employs the SIC technique to remove the direct-link signal \(\mathbf{h}_{\text{IR}}^{eq}{\boldsymbol{W}}\mathbf{s}(l)\) [4], [35]. Subsequently, the two hypotheses can be written as \[{\bar y_{\text{IR}}}(l) \leftarrow \left\{ {\begin{array}{*{20}{c}} {{{n_{\text{IR}}}(l)},\qquad \qquad {\quad}{{H_0}}}\\ {{{\mathbf{f}_{\text{b},\text{IR}}} {\boldsymbol{W}}\mathbf{s}(l) + {n_{\text{IR}}}(l)},{} H_1} \end{array}}, \right. \label{eq7}\tag{6}\] where \(H_0\) is the null hypothesis corresponding to the OFF state in OOK, and \(H_1\) is the alternative hypothesis corresponding to the ON state. The detection performance at the IR is then evaluated in terms of the detection error probability, which is expressed as

\[\label{eq8} \xi = \Pr \left( {{{\cal B}_1}\left| {{H_0}} \right.} \right) + \Pr \left( {{{\cal B}_0}\left| {{H_1}} \right.} \right)\tag{7}\] where \(\Pr \left( {{{\cal B}_1}\left| {{H_0}} \right.} \right)\) and \(\Pr \left( {{{\cal B}_0}\left| {{H_1}} \right.} \right)\) denote the false alarm rate and miss detection rate, respectively. \({{\cal B}_1}\) and \({{\cal B}_0}\) represent the binary decisions that determine whether the backscatter link is present or not, respectively. Based on the Neyman-Pearson criterion, the likelihood ratio test is employed to minimize the detection error probability \(\xi\) [36], which can be written as

\[\label{eq9} \frac{{{P_1} \buildrel \Delta \over = \prod {_{l = 1}^Lf({{\bar y}_{\text{IR}}}(l)\left| {{H_1}} \right.)} }}{{{P_0} \buildrel \Delta \over = \prod {_{l = 1}^Lf({{\bar y}_{\text{IR}}}(l)\left| {{H_0}} \right.)} }}{\rm{ }}\frac{{\mathop{>} \limits^{{{\cal B}_0}} }}{{\mathop{<} \limits_{{{\cal B}_1}} }}1.\tag{8}\] The likelihood functions of \({\bar y_{\text{IR}}}(l)\) in \(H_0\) and \(H_1\) are denoted as \(f({\bar y_{\text{IR}}}(l)\left| {{H_0}} \right.) \sim CN(0,\delta _{{\rm{\text{IR}}}}^2)\) and \(f({\bar y_{\text{IR}}}(l)\left| {{H_1}} \right.) \sim CN(0,{\gamma _b} + \delta _{{\rm{\text{IR}}}}^2)\) with \({\gamma _b} = {\left\| {{f_{\text{BD},{\text{IR}}}}\mathbf{h}_{\text{BD}}^{eq}{\boldsymbol{W}}} \right\|^2}\), respectively. Then, the minimum detection error rate \(P_e\) can be derived from (7 ) and (8 ). However, as the resultant expression of \(P_e\) involves the incomplete Gamma function, it poses challenges to further analytical and design efforts[36]. To address this, a tractable lower bound on \(P_e\) is obtained according to [37], expressed as follows

\[\label{eq10} {P_e} \le 1 - \sqrt {\frac{1}{2}{\cal D}\left( {{P_0}\left\| {{P_1}} \right.} \right)},\tag{9}\] where \({\cal D}\left( {{P_0}\left\| {{P_1}} \right.} \right) = L\left[ {\ln (\frac{{{\gamma _b} + \delta _{{\rm{\text{IR}}}}^2}}{{\delta _{{\rm{\text{IR}}}}^2}}) + \frac{{\delta _{{\rm{\text{IR}}}}^2}}{{{\gamma _b} + \delta _{{\rm{\text{IR}}}}^2}} - 1} \right]\) denotes Kullback-Leibler (KL) divergence from \(P_0\) to \(P_1\). Hence, the detection constraint for the secondary transmission is derived as \({\cal D}\left( {{P_0}\left\| {{P_1}} \right.} \right) \ge 2{\varepsilon ^2}\), which is a more stringent constraint to guarantee \({P_e} \le 1 - \varepsilon\).

2.3 Problem Formulation↩︎

The objective of this paper is to jointly optimize the transmit beamforming at the AP and the pinching beamforming formed by the PAs, with the aim of maximizing the sum rate of PRs. The optimization problem is formulated as

\[\tag{10} \begin{align*} (\text{P1})\mathop{\max }\limits_{\boldsymbol{W}, \boldsymbol{X}}\, &\sum\limits_{k = 1}^K {{R_k}} \tag{11}\\ {\rm{s.t.}} &\operatorname{Tr}({\boldsymbol{W}}{{\boldsymbol{W}}^H}) \le {P_{\max }}, \tag{12}\\ &D({P_0}|{P_1}) \ge 2{\varepsilon ^2}, \tag{13}\\ &{x_{n,m + 1}} - {x_{n,m}} \ge {d_{\min }},\forall n,m, \tag{14}\\ &0 \le {x_{n,m}} \le {S_x},\forall n,m, \tag{15} \end{align*}\] where (12 ) represents the maximum transmit power constraint, (13 ) ensures the minimum detection error rate at the IR, (14 ) imposes a minimum antenna spacing \(d_{\text{min}}\) to avoid mutual coupling between adjacent PAs, and (15 ) guarantees that the positions of the PAs are within the maximum range of the connected waveguide.

The optimization problem (P1) is highly non-convex and intractable due to the fractional expressions and multivariable coupling in both the objective function and the constraints. In the following section, we first introduce a learning-based beamforming framework to solve problem (P1). Furthermore, an alternating optimization approach is adopted to decompose the original problem into two sub-problems, which are then solved iteratively using SCA and PSO methods.

3 Gradient Descent-based Joint Beamforming Framework↩︎

In this section, we propose an LGD method to solve the joint beamforming problem (P1). Specifically, we first address the constraints and equivalently transform the constrained optimization problem into an unconstrained form. Subsequently, we present the proposed LGD-based beamforming design algorithm, which is implemented using self-defined neural network layers, where the optimization variables are treated as learnable parameters.

3.1 GD-based Reformulation↩︎

A major challenge in designing the pinching beamforming lies in efficiently handling constraints (14 ) and (15 ). To overcome this, we reformulate the original optimization variable \(\mathbf{x}_n\) into a more tractable form by introducing an offset variable \(\Delta \mathbf{x}_n\). Specifically, inspired by the minimum distance constraints of PAs, constraint (14 ) can be reformulated as \({x_{n,m + 1}} \ge {x_{n,m}} + {d_{\min }} \to {x_{n,m + 1}} = {x_{n,m}} + {d_{\min }} + \Delta {x_{n,m + 1}}\), where \(\Delta {x_{n,m + 1}}\) denotes the non-negative offset of the \((m+1)\)-th PA relative to \(m\)-th PA at the \(n\)-th waveguide. Hence, we define the position of the first PA \({x_{n,1}} = \Delta {x_{n,1}},\Delta {x_{n,1}} \ge 0\), and constraint (14 ) can be further rewritten as

\[\label{eq12} \left\{ \begin{align} &x_{n,1} = \Delta x_{n,1}, \\ &x_{n,m} = (m - 1){d_{\min }} + \sum_{j = 1}^m \Delta x_{n,j},\quad m = 2,\ldots,M \end{align} \right.\tag{16}\] for \(\forall n\). For clarity, Fig. 2 illustrates the parameter mapping of the PA position matrix in the \(n\)-th waveguide, with the offsets \(\Delta {x_{n,m}}\), \(m=1,…,M\), being defined as the optimization variables.

However, \(\left\{ {\Delta {x_{n,m}}} \right\}_{m = 1}^M\) are also restricted by the maximum length constraint (15 ), such as \({x_{n,M}} \le {S_x}\). which can be equivalently written as

\[\label{eq13} \sum {_{j = 1}^m\Delta {x_{n,j}}} \le {\Delta _{\max }},\tag{17}\] where \({\Delta _{\max }} = {S_x} - (m - 1){d_{\min }}\). When this constraint is not satisfied, the optimized \(\Delta {\mathbf{x}_n}\) will be normalized through the softmax function, i.e., \(f_{\text{SM}}(\Delta {x_{n,m}}) = \frac{{\Delta {x_{n,m}}}}{{\sum {_{i = 1}^M\Delta {x_{n,i}}} }}{\Delta _{\max }}\). Based on the above transformation process, optimizing \(x_n\) is equivalent to optimizing the corresponding offset \(\Delta \mathbf{x}_{n}\), which efficiently handles the constraints (14 ) and (15 ).

For the transmit beamforming, the projected gradient descent method [38] is adapted to deal with the maximum transmit power constraint (12 ). Let the constraint set be denoted by \({\cal C} \buildrel \Delta \over = \left\{ {\left. {\boldsymbol{W}} \right|\operatorname{Tr}({\boldsymbol{W}}{{\boldsymbol{W}}^H}) \le {P_{\max }}} \right\}\), the projection operation \(\prod {_{\cal C}}\) can be expressed as

\[\label{eq14} \prod {_{\mathcal{C}}} \{ W\} = \begin{cases} {\boldsymbol{W}}, & \text{if } \operatorname{Tr}({\boldsymbol{W}} {\boldsymbol{W}}^H) \le P_{\max} \\ \frac{{\boldsymbol{W}}}{\| {\boldsymbol{W}} \|}\sqrt{P_{\max}}, & \text{otherwise} \end{cases}\tag{18}\] Moreover, following [39], the penalty method can be adopted to guarantee the constraint (13 ). Specifically, we introduce a penalty parameter \(\xi\) to the objective function (11 ), allowing the original problem to be reformulated as follows

\[\label{eq15} (\text{P1.1})\mathop{\min }\limits_{\boldsymbol{W}, \boldsymbol{X}} {\cal F} \buildrel \Delta \over = -\sum\limits_{k = 1}^K {{R_k}} + \xi {[\max (0,2{\varepsilon ^2} - D({P_0}\left| {{P_1}} \right.))]^2}\tag{19}\] where \(\xi\) controls the penalty magnitude. It is important to note that the objective function \({\cal F}\) is differentiable, and its gradient with respect to the transmit beamforming \({\boldsymbol{W}}\) can be denoted as

\[\label{eq16} {\nabla _{\boldsymbol{W}}}{\cal F} = \sum\limits_{k = 1}^K {{\nabla _{\boldsymbol{W}}}{R_k}} - 2\xi \frac{{\partial D({P_0}|{P_1})}}{{\partial {\gamma _b}}} \cdot \frac{{\partial {\gamma _b}}}{{\partial {\boldsymbol{W}}}}\tag{20}\] Thus, the update of \({\boldsymbol{W}}\) is obtained by

\[\label{eq17} {\widetilde{\boldsymbol{W}}^{(i)}} = {\boldsymbol{W}}^{(i-1)} - {\eta _1}{\nabla _{\boldsymbol{W}}}{\cal F}\tag{21}\] \[\label{eq18} {{\boldsymbol{W}}^{(i)}} = \prod {_{\cal C}} ({\widetilde{\boldsymbol{W}}^{(i)}})\tag{22}\] where \({\eta _1}\) is the step size. Then, the gradient vector for the objective function with respect to the pinching beamforming \(\Delta {\mathbf{x}_n}\) is given by

\[\label{eq19} {\nabla _{\Delta {\mathbf{x}_n}}}{\cal F} = {[\frac{{\partial {\cal F}}}{{\partial \Delta {x_{n,1}}}},\frac{{\partial {\cal F}}}{{\partial \Delta {x_{n,2}}}},...,\frac{{\partial {\cal F}}}{{\partial \Delta {x_{n,M}}}}]^T},\forall n.\tag{23}\] Notably, the derivation of \({\cal F}\) involves the chain rule, rendering a closed-form analytical expression intractable. Instead, by computing the gradient \({\nabla _{\Delta {\mathbf{x}_n}}}{\cal F}\), the variable \({\Delta {\mathbf{x}_n}}\) can be updated by

\[\label{blkmwcsu} \widetilde{\Delta} {\mathbf{x}_n}^{(i)} = \Delta {\mathbf{x}_n}^{(i - 1)} - {\eta _2}{\nabla _{\Delta {\mathbf{x}_n}}}{\cal F},\tag{24}\] \[\label{vyqesdlx} \Delta {\mathbf{x}_n}^{(i)} = {f_{\text{SM}}}\left( {\widetilde{\Delta} {\mathbf{x}_n}^{(i)}} \right),\tag{25}\] where \({\eta _2}\) also denote the step size. After obtaining the optimized \(\Delta {\mathbf{x}_n}\), we can recover \({{x_{n}}}\) from (16 ).

3.2 Learning-aided GD Approach↩︎

In conventional GD, computing and projecting the gradients as described in (16) and (19) is challenging and often intractable. In this paper, we propose a learning-aided GD (LGD) framework that introduces two self-defined layers wherein the optimization variables are directly modeled as learnable parameters, instead of converting the original problem into an alternating maximization framework. This effectively avoids complex matrix inversion operations and alternating optimization between variables. Furthermore, the proposed LGD framework leverages automatic differentiation, enabling efficient and broadly applicable gradient descent solutions for practical optimization tasks [40].

Fig. 3 illustrates the architecture of the proposed LGD for joint beamforming, where \(\Delta \boldsymbol{X}\) layer and \({\boldsymbol{W}}\) layer are implemented, with their parameters corresponding to the position offsets of PAs and the transmit beamforming matrix, which are also updated at each iteration by back-propagating the loss function \({\cal F}\). Note that the channel location information of the BD, IR, and PRs constitutes the known parameter set \(\mathbf{a} = \{ \mathbf{I}_k^{\text{U}},{\mathbf{I}_{\text{BD}}},{\mathbf{I}_{\text{IR}}}\}\), which is incorporated into the computation of \({\cal F}\). Moreover, the learnable parameters are not trained in advance using a dataset, but instead are adapted in a task-driven manner using gradient descent on problem (P1), without requiring labeled data. As a result, each update of the learnable parameter set \(\boldsymbol{\Theta} = \{ \Delta \boldsymbol{X},\boldsymbol{W}\}\) can be regarded as a form of implicit training embedded within the optimization process.

As can be seen from the Fig. 3, the constraints are achieved through different activation functions. For the maximum operator, i.e., \(\max ( \cdot , \cdot )\) in (19 ), the rectified linear unit (ReLU) function can be adopted to replace them to facilitate the network training, and the Softmax activation function is employed to normalize the output of \(\Delta {\boldsymbol{\widetilde{X}}^{(i)}}\). Moreover, we design a custom projection activation function that normalizes the output \({\boldsymbol{W}}\) according to the operation defined in (22 ). Consequently, (16 ) can be implemented using the cumsum function in PyTorch. Furthermore, a fixed learning rate that is too large causes oscillations, while too small slows convergence [41]. In contrast, Adam demonstrates superior convergence speed and performance, and is less sensitive to the initial learning rate due to its adaptive learning rate mechanism. The learnable parameters \(\boldsymbol{\Theta}\) are updated by the Adam optimizer as follows

\[\label{eq22} {\boldsymbol{\Theta}^{(i)}} = {\boldsymbol{\Theta}^{(i - 1)}} + \eta \cdot \text{Adam}({\nabla _\boldsymbol{\Theta}}{\cal F}),\tag{26}\] where \(\eta\) denotes the initial learning rate. Upon completion of the training phase, the optimized variables \(\boldsymbol{X}\) and \({\boldsymbol{W}}\) can be efficiently reconstructed from the learned network parameters and their associated transformations.

The conventional neural network approach is typically trained offline with a large number of samples to minimize the average loss over the entire dataset. However, its performance may be limited by insufficient training data or convergence to local minima. In contrast, the proposed LGD method can be employed as a plug-and-play tool without requiring any pre-collected training dataset. Each update of the parameter set \(\boldsymbol{\Theta}\) corresponds directly to a step in solving the original optimization problem, while the update direction and step size are still determined in accordance with the principles of GD [42]. In addition, LGD exhibits strong interpretability, as its iterative variable update process is inherently explainable and can be seamlessly integrated with expert knowledge and prior information [40]. Thanks to its lightweight network architecture, the proposed LGD method achieves low computational complexity, which can be characterized by \(\mathcal{O}\left(I_G \big(KN^2M + K^2N\big)\right)\), where \(I_G\) denotes the number of training iterations [42].

4 Proposed Two-stage Iterative Algorithm↩︎

Although the proposed LGD adopts a lightweight network architecture with low computational complexity, solving problem (P1) via gradient descent can easily lead to convergence to local optima. This is primarily due to the highly multimodal nature of the PA position optimization, which leads to a large number of local minima[22]. To further enhance performance, we propose a two-stage SCA-PSO algorithm. Given the strong coupling between the two optimization variables in (P1), the original problem is first decoupled into two subproblems, which are then solved in an alternating manner. Specifically, in the first stage, given the PA position matrix \(\boldsymbol{X}\), the transmit beamforming matrix \({\boldsymbol{W}}\) is optimized by using the SCA method. In the second stage, the PSO method is developed to optimize \(\boldsymbol{X}\) based on the optimized \({\boldsymbol{W}}\).

4.1 Transmit Beamforming Optimization↩︎

With fixed PA position matrix \(\boldsymbol{X}\), the transmit beamforming subproblem can be written as

\[\tag{27} \begin{align*} (\text{P2.1})\mathop{\max } \limits_{{\boldsymbol{W}}}\, &{\quad}\sum\limits_{k = 1}^K{\log _2}(1 + \frac{{{{\left| {\mathbf{h}_k^{eq}{\mathbf{w}_k}} \right|}^2}}}{{\sum {_{i = 1,i \ne k}^K{{\left| {\mathbf{h}_k^{eq}{\mathbf{w}_i}} \right|}^2} + \delta _k^2} }}) \\ &+ {\log _2}(1 + \frac{{{{\left| {(\mathbf{h}_k^{eq} + {\mathbf{f}_{b,k}}){\mathbf{w}_k}} \right|}^2}}}{{\sum {_{i = 1,i \ne k}^K{{\left| {(\mathbf{h}_k^{eq} + {\mathbf{f}_{b,k}}){\mathbf{w}_i}} \right|}^2} + \delta _k^2} }}). \tag{28}\\ {\rm{s.t.}} \, &(\text{\ref{11b}}), (\text{\ref{11c}}) \end{align*}\] However, it can be seen that the objective function (27 ) and constraint (13 ) are non-convex due to the quadratic terms, rendering the problem intractable to solve directly. For convenience, we introduce auxiliary variables \({t_{1,k}} = \frac{{{u_{1,k}}}^2}{{{v_{1,k}}}}\), \({t_{2,k}} = \frac{{{u_{2,k}}}^2}{{{v_{2,k}}}}\), where \({u_{1,k}} = {\mathbf{h}_{k}^{eq}}{\mathbf{w}_k}\), \({v_{1,k}} = \sum {_{i = 1,i \ne k}^K{{\left| {{ \mathbf{h}_{k}^{eq}}{\mathbf{w}_i}} \right|}^2} + \delta _k^2}\), \({v_{1,k}} ={\left| {\left( {\mathbf{h}_k^{eq} + {\mathbf{f}_{b,k}}} \right){\mathbf{w}_k}} \right|^2}\), and \({v_{2,k}} = \sum {_{i = 1,i \ne k}^K{{\left| {\left( {\mathbf{h}_k^{eq} + {\mathbf{f}_{b,k}}} \right){\mathbf{w}_i}} \right|}^2} + \delta _k^2}\). Subsequently, the first-order Taylor expansion is applied to derive a convex lower bound for \(t_{1,k}\) and \(t_{2,k}\), which can be expressed as follows

\[\label{eq24} 2\frac{{{{\bar u}_{1,k}}^*{u_{1,k}}}}{{{{\bar v}_{1,k}}}} - \frac{{{{\left| {{{\bar u}_{1,k}}} \right|}^2}}}{{{{\left( {{{\bar v}_{1,k}}} \right)}^2}}}{v_{1,k}} \ge {t_{1,k}},\tag{29}\] where \({\bar u}_{1,k}\) and \({{\bar v}_{1,k}}\) represents the value of \({u_{1,k}}\) and \({v_{1,k}}\) at the previous iteration, respectively. Similarly, \(t_{2,k}\) can be approximated by \(2\frac{{{{\bar u}_{2,k}}^*{u_{2,k}}}}{{{{\bar v}_{2,k}}}} - \frac{{{{\left| {{{\bar u}_{2,k}}} \right|}^2}}}{{{{\left( {{{\bar v}_{2,k}}} \right)}^2}}}{v_{2,k}}\). Moreover, based on (9 ), we can introduce \(D({P_0}|{P_1}) = L(\ln \rho + {1 \mathord{\left/ {\vphantom {1 \rho }} \right. \kern-\nulldelimiterspace} \rho } - 1)\), where \(\rho = {{({\gamma _b} + \delta _{{\rm{\text{IR}}}}^2)} \mathord{\left/ {\vphantom {{({\gamma _b} + \delta _{{\rm{\text{IR}}}}^2)} {\delta _{{\rm{\text{IR}}}}^2}}} \right. \kern-\nulldelimiterspace} {\delta _{{\rm{\text{IR}}}}^2}}\). Accordingly, constraint (13 ) can be simplified as

\[\label{eq225} \ln \rho + {1 \mathord{\left/ {\vphantom {1 \rho }} \right. \kern-\nulldelimiterspace} \rho } - 1 \ge {{2{\varepsilon ^2}} \mathord{\left/ {\vphantom {{2{\varepsilon ^2}} L}} \right. \kern-\nulldelimiterspace} L},\tag{30}\] which can be rewritten as \(\frac{{{{\left\| {{f_{\text{BD},{\text{IR}}}}\mathbf{h}_{\text{BD}}^{eq}{\boldsymbol{W}}} \right\|}^2} + \delta _{{\rm{\text{IR}}}}^2}}{{\delta _{{\rm{\text{IR}}}}^2}} \le {\overline{a} _0}\) or \(\frac{{{{\left\| {{f_{\text{BD},{\text{IR}}}}\mathbf{h}_{\text{BD}}^{eq}{\boldsymbol{W}}} \right\|}^2} + \delta _{{\rm{\text{IR}}}}^2}}{{\delta _{{\rm{\text{IR}}}}^2}} \ge {\overline{a} _1}\), where \({\overline{a} _0}\) and \({\overline{a} _1}\) are the two roots of the function \(\ln \rho + {1 \mathord{\left/ {\vphantom {1 \rho }} \right. \kern-\nulldelimiterspace} \rho } - 1 = {{2{\varepsilon ^2}} \mathord{\left/ {\vphantom {{2{\varepsilon ^2}} L}} \right. \kern-\nulldelimiterspace} L}\). It can be readily observed that \({\overline{a} _0} \le 1 \le \rho\), and therefore, constraint (13 ) can be reformulated as

\[\label{eq25} {\left \| {{f_{\text{BD},{\text{IR}}}}\mathbf{h}_{\text{BD}}^{eq}{\boldsymbol{W}}} \right \|^2} \ge \left( {{{\overline{a} }_1} - 1} \right)\delta _{{\rm{\text{IR}}}}^2.\tag{31}\] To handle this non-convex constraint, the quadratic term \({\left\| {{f_{\text{BD},{\text{IR}}}}\mathbf{h}_{\text{BD}}^{eq}{\boldsymbol{W}}} \right\|^2}\) is also approximated as \(\varpi \buildrel \Delta \over = 2{\rm{Re}}\left\{ \operatorname{Tr}({{{({f_{\text{BD},{\text{IR}}}}\mathbf{h}_{\text{BD}}^{eq}\bar {\boldsymbol{W}})}^H}{f_{\text{BD},{\text{IR}}}}\mathbf{h}_{\text{BD}}^{eq}{\boldsymbol{W}})} \right\} - {\left\| {{f_{\text{BD},{\text{IR}}}}\mathbf{h}_{\text{BD}}^{eq}\bar {\boldsymbol{W}}} \right\|^2}\) by using the first-order Taylor expansion, where \({\rm{Re}}\left\{ \cdot \right\}\) denotes the real part of the corresponding variable. Subsequently, the problem (P2.1) can be efficiently solved by utilizing the SCA algorithm. The approximated convex problem can be written as

\[\tag{32} \begin{align*} (\text{P2.2})\mathop{\max } \limits_{X}\, &{\quad}\sum\limits_{k = 1}^K{\log _2}(1 + {t_{1,k}}) + {\log _2}(1 + {t_{2,k}}) \tag{33}\\ {\rm{s.t.}} \,&{\left\| {{f_{\text{BD},{\text{IR}}}}\mathbf{h}_{\text{BD}}^{eq}{\boldsymbol{W}}} \right\|^2} \ge \left( {{{\overline{a} }_1} - 1} \right)\delta _{{\rm{\text{IR}}}}^2, \tag{34}\\ &(\text{\ref{11b}}), \end{align*}\] which can be directly solved by the existing convex optimization tools such as CVX [43]. The detailed SCA algorithm for solving the problem (P2.1) is summarized in , where the initial \({\boldsymbol{W}}\) is randomly initialized in the feasible region.

4.2 Pinching Beamforming Optimization↩︎

With fixed transmit beamforming matrix \({\boldsymbol{W}}\), the pinching beamforming subproblem with respect to \(\boldsymbol{X}\) can be formulated as

\[\tag{35} \begin{align*} (\text{P2.3})\mathop{\max } \limits_{\boldsymbol{X}}\, &\sum\limits_{k = 1}^K {{R_k}}. \tag{36}\\ {\rm{s.t.}} \,&(\text{\ref{11d}}), (\text{\ref{11e}}), (\text{\ref{26b}}) \end{align*}\] which is challenging due to the ill-conditioned constraints and coupled variables in the objective function. To address this issue, we adopt the PSO method to search for the optimal positions of the PAs, where PSO is a population-based stochastic optimization algorithm that mimics the social behavior of swarms to search for the optimal solution[44]. While the exhaustive search guarantees the global optimum, its computational complexity is prohibitively high. PSO efficiently explores the solution space through swarm intelligence, achieving near-optimal performance with significantly reduced computational cost. Specifically, each waveguide is associated with a swarm, where each swarm consists of \(Q\) particles. Taking a single waveguide as an example, we begin by randomly initializing \(Q\) particles with positions \(\mathbf{x}_q^{(t)} = {[x_{q,1}^{(t)},x_{q,2}^{(t)},...,x_{q,M}^{(t)}]^T}\) and velocities \(\mathbf{v}_q^{(t)} = {[v_{q,1}^{(t)},v_{q,2}^{(t)},...,v_{q,M}^{(t)}]^T}\), \(q = 1,...,Q\) within the feasible search space. Here, \(x_{q,m}^{(t)}\) and \(v_{q,m}^{(t)}\) denote the position and the update velocity of \(m\)-th PA in the \(q\)-th particle during the t-th iteration, respectively.

To satisfy constraint (15 ), all values of \(x_{q,m}^{(t)}\) are restricted within the range \([0,{S_x}]\). Based on the PSO algorithm framework [44], [45], each particle updates its position according to the current personal best position \({\widetilde{x}_{q,{\rm{pb}}}}\) and the swarm global best position \({\widetilde{x}_{g{\rm{b}}}}\). Accordingly, in \((t+1)\)-th iteration, the update process for each particle’s velocity and position is formulated as follows

\[\label{eq28} \mathbf{v}_q^{(t + 1)} = {\omega _0}\mathbf{v}_q^{(t)} + {\omega _1}{c_1}({\widetilde{x}_{q,{\rm{pb}}}} - \mathbf{x}_q^{(t)}) + {\omega _2}{c_2}({\widetilde{x}_{{\rm{gb}}}} - \mathbf{x}_q^{(t)}),\tag{37}\] \[\label{eq29} \mathbf{x}_q^{(t + 1)} = \mathbf{x}_q^{(t)} + \mathbf{v}_q^{(t + 1)},\tag{38}\] where \({\omega _0}\) is the inertia weight that regulates the momentum of the particle and is defined as \({\omega _0} = {\omega _{\max }} - ({\omega _{\max }} - {\omega _{\min }}){t \mathord{\left/ {\vphantom {t T}} \right. \kern-\nulldelimiterspace} T}\), where \({\omega _{\max }}\) and \({\omega _{\min }}\) represent the upper bound and lower bound of \({\omega _0}\), and \(T\) denotes the maximum iteration number. \({\omega _1}\) and \({\omega _2}\) are random factors uniformly distributed in the range \([0,1]\), which improve the randomness of the search process to avoid premature convergence. The parameters \(c_1\) and \(c_2\) act as the personal and global learning factors, respectively, regulating the influence of the personal best and global best positions on the velocity update.

In each iteration, the fitness value of the \(q\)-th particle is evaluated using equation (35 ), based on its current position \(\boldsymbol{X}_q\). To enforce constraints (14 ) and (34 ), an adaptive penalty method integrates the constraint violations into the objective function as penalty terms. The resulting penalized objective function is given by

\[\label{eq30} {\cal L}({\boldsymbol{X}_q}) = \sum\limits_{k = 1}^K {{R_k}({\boldsymbol{X}_q})} - \mu \left| {{\cal P}({\boldsymbol{X}_q})} \right|,\tag{39}\] where \({\cal P}({\boldsymbol{X}_q})\) denotes the set of penalty terms associated with violations of the minimum detection error rate constraint (34 ) and the minimum PA spacing constraint (14 ). Specifically, \(\mathcal{P}(\boldsymbol{X}_q)\) is defined as

\[\begin{align} \label{eq31} \nonumber {\cal P}({\boldsymbol{X}_q}) = \{ {x_{q,n,m}}&\left| {\{ {x_{q,n,m + 1}} - } \right.{x_{q,n,m}} < {d_{\min }},\forall n,m\} ,\\ &+D({P_0}|{P_1}) < 2{\varepsilon ^2}\} , \end{align}\tag{40}\] where \(\mu > 0\) is a sufficiently large penalty factor that drives particles toward feasible regions. If a particle violates either constraint, the resulting fitness value \(\mathcal{L}(\boldsymbol{X}_q)\) is penalized accordingly, potentially reducing it below zero to discourage infeasible solutions. As each particle is evaluated, its personal best and the global best positions are progressively updated until convergence is achieved. The detailed PSO algorithm for solving problem (P2.3) is summarized in Algorithm 2.

Based on the above analysis, the problem (P1) can be effectively solved using the proposed SCA-PSO algorithm, where the transmit beamforming matrix and the pinching position matrix are alternately optimized until reaching convergence or the maximum number of iterations.

4.3 Convergence and Complexity Analysis↩︎

Since the proposed SCA-PSO algorithm operates in two stages, its convergence behavior depends on the performance of the SCA-based algorithm in the first stage and the PSO-based algorithm in the second stage. Let \(f({\boldsymbol{W}^{(j)}},{\boldsymbol{X}^{(j)}})\) denote the objective function value at the \(j\)-th iteration. Then, the following inequality holds: \[\label{eq32} f({\boldsymbol{W}^{(j - 1)}},{\boldsymbol{X}^{(j - 1)}})\mathop{\le} \limits^{(a)} f({\boldsymbol{W}^{(j)}},{\boldsymbol{X}^{(j - 1)}})\mathop{\le} \limits^{(b)} f({\boldsymbol{W}^{(j)}},{\boldsymbol{X}^{(j)}}),\tag{41}\] where the inequality in (a) arises from the fact that problem (P2.2) is solved to optimality in each iteration, and its optimal objective value provides a lower bound for that of problem (P2.1). The inequality (b) holds because the fitness value of the global best position is non-decreasing throughout the iterations of Algorithm 2, i.e., \({\cal L}({\boldsymbol{X}}_q^{(t+1)}) \ge {\cal L}({\boldsymbol{X}}_q^{(t)})\). If the updated pinching-antenna position obtained from the PSO algorithm fails to yield an improvement in sum rate, the previous position is preserved [30]. Moreover, since the objective value of problem (P2.3) is bounded above, the sequence of best fitness values converges, ensuring the convergence of the PSO-based algorithm.

The computational complexity of the proposed SCA-PSO algorithm is primarily determined by the subproblems (P2.2) and (P2.3). Specifically, as indicated in [46], [47], the subproblem (P2.2) is a second-order cone programming (SOCP) problem, which can be efficiently solved via an interior-point method with a complexity of \({\log (\frac{1}{{{\varepsilon _1}}}){{({K^2}N)}^{3.5}}}\). To address the subproblem (P2.3), the PSO algorithm used for updating \(\boldsymbol{X}\) has a computational complexity of \({\mathcal{O}}\left( {TQMN} \right)\). Therefore, the overall complexity of the proposed algorithm is given by \({\mathcal{O}}\left( {TQMN + \log (\frac{1}{{{\varepsilon _1}}}){{({K^2}N)}^{3.5}}} \right)\).

5 Numerical Results↩︎

In this section, we present the simulation setup and corresponding numerical results to demonstrate the effectiveness of the proposed PASS-enabled SR design and associated algorithms. For consistency with prior work, the system parameters are adopted from [25] and [48]. Specifically, the operating frequency is set to \(f=28 \,\text{GHz}\), the ma transmit power, the noise power is \(\delta _k^2 = \delta _{\text{IR}}^2 = - 80 \, \text{dBm}\), and the effective index of the waveguide is \(n_{\text{eff}}=1.4\). We assume that the number of PRs and chains/waveguides is \(N = K = \{ 2,4\}\). Each waveguide is equipped with \(M=3\) PAs, and the heights of all PAs are fixed at \(z^{\text{PA}}=5 m\). Both the PRs and the IR are randomly deployed within a rectangular area of size \({S_x} \times {S_y} = 30 \times 4~\text{m}^2\), while the BD is placed at a fixed location of \((5, 2, 2) \, \text{m}\). The waveguides are uniformly distributed along the y-axis (vertical direction) with a consistent interval of \({{{S_{\rm{y}}}} \mathord{\left/ {\vphantom {{{S_{\rm{y}}}} N}} \right. \kern-\nulldelimiterspace} N}\) meters and the minimum spacing of PA positions is set to \(d_{\min}\) = 0.1 m. The symbol duration ratio between primary and secondary transmissions is \(L = 60\), and the detection probability threshold is \(\varepsilon = 0.95\). The convergence tolerance and penalty factor in Algorithm 1 are set to \(\varepsilon_1 = 10^{-3}\) and \(\mu=10\), respectively. Numerical results are computed as the average over 100 independently generated channel realizations.

For performance comparison, the following benchmark schemes are considered: 1) Element-wise algorithm: Each antenna position \(x_{n,m}\) is optimized individually using a one-dimensional search to obtain a near-optimal solution[22]. Specifically, denote by \({{\cal D}_x}\) a uniform sampled feasible grid over the interval \([0,{S_x}]\). The optimization with respect to \(x_{n,m}\) can be formulated as \(\mathop{\max }\limits_{{x_{n,m}} \in {{\cal D}_x}} {\rm{ }}\sum {_{k = 1}^K} {R_k}\), where the suboptimal solution \(x_{n,m}^ *\) is obtained by searching over the feasible discrete grid \({{\cal D}_x}\), subject to the constraints (13 )-(15 ). 2) Fixed-PA: The PAs are uniformly distributed along the x-axis within each waveguide, with an inter-element spacing \({{{S_x}} \mathord{\left/ {\vphantom {{{S_x}} M}} \right. \kern-\nulldelimiterspace} M}\). 3) Massive MIMO: The massive MIMO BS with the hybrid beamforming architecture is positioned at the origin point, where \(N\) RF chains are deployed, each connected to \(M\) antennas through phase shifters. 4) Conventional MIMO: The MIMO BS equipped with \(N\) antennas is placed at the origin point, each connected to a dedicated RF chain.

Fig. 6 presents the convergence characteristics of the proposed algorithms. As shown in Fig. 6 (a), the LGD algorithm demonstrates different behaviors under varying initial learning rates. The results show that although a larger initial learning rate facilitates faster initial ascent during the early iterations, it often suffers from pronounced fluctuations and is more prone to convergence to suboptimal solutions. In contrast, the LGD algorithm with an initial learning rate \(\eta=10^{-4}\) leads to a more stable convergence trajectory and ultimately achieves a higher sum rate. Fig. 6 (b) compares the convergence performance of the proposed SCA-PSO algorithm against the element-wise benchmark. It can be observed that both approaches exhibit comparable convergence speeds. While the element-wise method achieves slightly higher performance, the improvement is marginal. In contrast, the SCA-PSO algorithm achieves a sum rate performance close to that of the element-wise method with significantly fewer iterations and much lower computational complexity. The element-wise method requires exhaustive updates for each antenna element and involves a larger solution space, resulting in prohibitively high computational complexity. Therefore, the proposed SCA-PSO approach provides a more practical and efficient alternative with competitive performance.

Fig. 7 compares the achievable sum rates of different approaches as the number of active pinches/antennas per waveguide (or RF chain) increases. For all considered schemes, the sum rate improves with the number of active antennas, confirming the benefits of increased spatial degrees of freedom for beamforming and interference mitigation. In both the \(K=2\) and \(K=4\) user cases, the sum rate also increases with \(K\) due to spatial multiplexing gains. Leveraging PASS, the proposed SCA-PSO and LGD schemes achieve substantial sum-rate improvements over the baselines, primarily by reconfiguring large-scale path loss through flexible position optimization. When the maximum range of each waveguide is fixed, the performance gap between the LGD method and the massive MIMO scheme gradually narrows as the number of antennas increases.

Fig. 8 presents the achievable sum rates of different algorithms under varying maximum transmit power \(P_{\max}\) . For all considered schemes, the achievable rate increases with \(P_{\max}\), as higher transmit power enhances the received signal strength and improves the SINR. The proposed SCA-PSO and LGD algorithms consistently outperform the conventional MIMO and massive MIMO baselines across the entire transmit power range. This confirms the effectiveness of the proposed optimization framework in alleviating large-scale path loss significantly. Among them, SCA-PSO achieves higher performance than LGD by leveraging its stochastic–deterministic search strategy to obtain a high-quality suboptimal solution. In contrast, LGD is more prone to becoming trapped in local optima near the initial point, which can lead to performance degradation. While the element-wise algorithm attains the highest sum rate overall, it comes with significantly higher computational complexity due to its exhaustive search over individual antenna elements, which may limit its practical deployment. In contrast, the proposed SCA-PSO strikes a favorable balance between performance and computational efficiency.

In Fig. 9, we compare the achievable sum rate of different schemes as the distance \(S_x\) increases from 10 m to 30 m. It can be observed that all schemes experience performance degradation with the increase of distance. However, the sum rate achieved by PASS exhibits only a marginal decrease in the sum rate. This degradation is caused by increasing channel attenuation on the BD-to-PRs/IR links, while the path loss on the dominant primary transmission paths from the BS is effectively mitigated through flexible placement of the PAs. In contrast, the performance of MIMO and Massive MIMO decreases rapidly as the distance increases from 10 m to 30 m, with a sum rate reduction of approximately 21%, mainly due to the free-space path loss effect. Massive MIMO achieves competitive performance at shorter ranges due to its hybrid beamforming architecture, and the Fixed-PA scheme exhibits fluctuating performance as its fixed antenna positions cannot adapt to varying link conditions. Thanks to the flexibility of PA deployment, PASS not only minimizes the long-distance path loss on the primary transmission links but also alleviates the adverse effects of double fading in backscattering links.

We investigate the impact of the symbol period ratio \(L\), as illustrated in Fig. 10. It can be observed that the performance of all schemes degrades as \(L\) decreases. This is because a smaller \(L\) implies fewer information symbols \(s(l)\) within one \(c\)’s symbol period, which in turn deteriorates the decoding performance of the secondary signal and ultimately limits the overall system performance. As observed, the proposed schemes achieve a higher sum rate compared to the fixed-PA scheme, demonstrating the effectiveness of optimizing the positions of the PAs. Moreover, in this context, PASS demonstrates a stronger capability in mitigating this issue, resulting in less performance degradation compared to conventional antenna systems. This also indicates the potential of PASS as an enabling technology for PSR systems, which is worth further exploration in future work.

In Fig. 9, we evaluate the impact of transmit SNR on the achievable sum rate. As expected, the sum rate of all schemes increases with higher SNR, owing to the improved signal strength relative to noise. Specifically, for the PASS-based system with the proposed schemes, when the transmit SNR is 65 dBm, the SCA-PSO and LGD algorithms achieve a 65.6% and 55.9% gain in sum rate, respectively, compared to conventional massive MIMO. This performance improvement is attributed to the flexibility of PASSs in dynamically reconfiguring pinching positions along the waveguide, thereby enabling the formation of strong and reliable LoS links. Furthermore, the performance gap between PASS-based schemes and massive MIMO widens as SNR increases. This is because, at low SNR, noise interference not only degrades the achievable data rate but also increases the detection error probability of the secondary signal.

6 Conclusion↩︎

We proposed a PASS-enabled downlink SR framework that jointly optimizes transmit and pinching beamforming to enhance both primary and secondary transmissions, where a PASS-based BS equipped with multiple waveguides, acting as the PT, serves both the IR and PRs with the aid of the BD. A joint optimization problem of transmit and pinching beamforming was formulated with the objective of maximizing the sum rate, subject to the detection error probability constraint at the IR and the feasible deployment region of the PAs. To address the highly coupled and nonconvex nature of this problem, we developed two solution strategies. The first is a low-complexity LGD method, which leverages end-to-end learning to transform parameters and incorporate constraints, thereby efficiently solving the constrained problem based on the principle of gradient descent. To further enhance performance, we introduced the SCA-PSO-based approach, where the transmit beamforming was optimized via SCA, followed by pinching beamforming optimization through a PSO-based search over feasible PA positions. Simulation results verified that the proposed PASS framework achieves notable performance gains over conventional fixed-antenna and massive MIMO schemes under strict detection error constraints. In particular, SCA-PSO attained performance close to the element-wise benchmark while significantly reducing computational complexity.

References↩︎