## 2 System Model and Formulation
### Formulations  
To keep the following simple, the study assumed that variables in Algorithm 1 were assigned initial values as \(time\_period=1\) and \(selected\_relay=R_{1}\) without losing general properties of Algorithm 1. It is worth noticing that we have sequential block time period \(T=\left\{T^{(1)},T^{(2)},\ldots,+\infty\right\}\in\mathbb{Z}^{+}\). Since assuming \(time\_period=1\) and \(selected\_relay=R_{1}\)we split the set \(T\) into two subsets, where \(T^{(odd)}=\left\{T^{(1)},T^{(3)},\dots\right\}\in\mathbb{Z}^{(odd)}\) for odd transmission block and \(T^{(even)}=\left\{T^{(2)},T^{(4)},\dots\right\}\in\mathbb{Z}^{(even)}\) for even transmission block. The EH at the receiver was provided as [36] (Equation (6)) for the PS protocol. From the aforementioned expressions, the EH at the coupled IoT relays in the first phase of the first time slot \(T_{1}^{(odd)}\in T^{(odd)}\) or \(T_{1}^{(even)}\in T^{(even)}\), is expressed as follows:  
\[EH_{S,R_{1}}^{T_{(odd)}^{(odd)}} =\eta\lambda_{R_{1}}P\sigma_{S,R_{1}}, \tag{1}\] \[EH_{S,R_{2}}^{T_{(odd)}^{(odd)}} =\eta\lambda_{R_{1}}P\sigma_{S,R_{2}},\] (2) \[EH_{S,R_{1}}^{T_{(even)}^{(even)}} =\eta\lambda_{R_{2}}P\sigma_{S,R_{1}},\] (3) \[EH_{S,R_{2}}^{T_{(even)}^{(even)}} =\eta\lambda_{R_{2}}P\sigma_{S,R_{2}}, \tag{4}\]  
where \(\lambda_{R_{1}}\) and \(\lambda_{R_{2}}\) are the PS factors (\(0<\lambda_{R_{1}}<1\) in the odd transmission block and \(0<\lambda_{R_{2}}<1\) in the even transmission block), \(P\) is power domain at IoT hub \(S\), \(\eta\) is collect factor for \(0\leq\eta\leq 1\), and expected channel gains \(\sigma_{S,R_{1}}\) and \(\sigma_{S,R_{2}}\) are given by \(\sigma_{S,R_{1}}=E\left\{\left|h_{S,R_{1}}\right|^{2}\right\}\) and \(\sigma_{S,R_{2}}=E\left\{\left|h_{S,R_{2}}\right|^{2}\right\}\), respectively. Furthermore, in the second phase of the first time slot \(T_{1}^{(odd)}\in T^{(odd)}\) or \(T_{1}^{(even)}\in T^{(even)}\), in terms of NOMA benefit, the IoT hub \(S\) encodes the messages \(x_{i}\) for \(i=\left\{N,\dots,1\right\}\) of IoT devices and superimposes the messages into the signal by superposition coding and sharing the remaining power domain \((1-\lambda_{R_{1}})P\) for \(T_{1}^{(odd)}\) or \((1-\lambda_{R_{2}})P\) for \(T_{1}^{(even)}\) with different PA factors \(\alpha_{i}\). Then, the received signals at the IoT relays \(R_{1}\) and \(R_{2}\) express as follows:  
\[y_{R_{1}}^{T_{(odd)}^{(odd)}} =h_{S,R_{1}}\sqrt{(1-\lambda_{R_{1}})P}\sum_{i=N}^{1}\sqrt{\alpha _{i}}x_{i}+n_{R_{1}}, \tag{5}\] \[y_{R_{2}}^{T_{(odd)}^{(odd)}} =h_{S,R_{2}}\sqrt{(1-\lambda_{R_{1}})P}\sum_{i=N}^{1}\sqrt{\alpha _{i}}x_{i}+n_{R_{2}},\] (6) \[y_{R_{1}}^{T_{(even)}^{(even)}} =h_{S,R_{1}}\sqrt{(1-\lambda_{R_{2}})P}\sum_{i=N}^{1}\sqrt{\alpha _{i}}x_{i}+n_{R_{1}},\] (7) \[y_{R_{2}}^{T_{(even)}^{(even)}} =h_{S,R_{2}}\sqrt{(1-\lambda_{R_{2}})P}\sum_{i=N}^{1}\sqrt{\alpha _{i}}x_{i}+n_{R_{2}}, \tag{8}\]  
where \(n_{R_{1}}\) and \(n_{R_{2}}\) are AWGN, i.e., \(n_{R_{1}}\sim CN(0,N_{0})\) and \(n_{R_{2}}\sim CN(0,N_{0})\), at IoT relays \(R_{1}\) and \(R_{2}\), respectively. Moreover, \(h_{S,R_{r}}=gd_{S,R_{r}}^{-\frac{5}{2}}\) is channel coefficient from IoT hub \(S\) to IoT relay \(R_{r}\) (\(r=\left\{1,2\right\}\)), \(g\) denoted Rayleigh distribution from IoT hub to coupled IoT relays with \(g\sim CN(0,1)\), \(e\) is path-loss exponent, \(d_{S,R_{1}}\) and \(d_{S,R_{2}}\) are distances from IoT hub \(S\) to IoT relays \(R_{1}\) and \(R_{2}\), respectively. Assuming IoT relays are closed to each other (\(d_{S,R_{1}}=d_{S,R_{2}}\)). From the system model as shown in Figure 1, the devices \(D_{1}\) and \(D_{N}\) are nearest and farthest from the IoT hub \(S\), respectively. In the concept of NOMA, the farthest device \(D_{N}\) was allocated the biggest PA factor since the farthest device \(D_{N}\) owned poorest CSI. Inspired by studies [21, 37], the PA factor \(\alpha_{i}\) for \(i=\left\{N,\dots,1\right\}\) carried the message \(x_{i}\) is given as follows:\[\alpha_{i}=i\left/\sum_{n=1}^{N}n,\right.\] (9) s.t.