## 3 Performance Analysis Framework
### Calculation of the MGF of the HSTS  
The MGF of the considered system model can be written, by using Equation (8), as:  
\[ \begin{array}{ll}M_{\gamma}(s)&=E_{\gamma}\left[e^{-s\frac{|h_{1}|^{2}|h_{2 }|^{2}}{|h_{2}|^{2}+\frac{E_{i}}{\sigma_{1}^{2}}|h_{2}|^{2}|h_{3}|^{2}+C} \frac{E_{s}}{\sigma_{1}^{2}}}\right]\\ &=\int\limits_{0}^{\infty}\int\limits_{0}^{\infty}\int\limits_{0}^{\infty}e^{ -s\frac{xy}{y\left(1+\frac{E_{i}}{\sigma_{1}^{2}}\right)}+\frac{E_{s}}{\sigma _{1}^{2}}}f_{|h_{1}|^{2}}(x)f_{|h_{2}|^{2}}(y)f_{|h_{3}|^{2}}(z)dxdydz\end{array}  \tag{9}\]  
We define the following integral from the above triple-integral:  
\[I_{1}\triangleq\int\limits_{0}^{\infty}e^{-s\frac{xy}{y\left(1+\frac{E_{i}}{ \sigma_{1}^{2}}\right)}+\frac{E_{s}}{\sigma_{1}^{2}}}f_{|h_{1}|^{2}}(x)dx \tag{10}\]  
By using identities from [23] (Equations (7.621.4, 9.121.1)), we obtain the following form for \(I_{1}\) as:  
\[ \begin{array}{ll}I_{1}&=\left(\alpha_{1}\left(1+z\frac{E_{i}}{\sigma_{1}^{ 2}}\right)y+\alpha_{1}C\right)\left(\left(s\frac{E_{s}}{\sigma_{1}^{2}}+\beta _{1}\left(1+z\frac{E_{i}}{\sigma_{1}^{2}}\right)\right)y+C\beta_{1}\right)^{m_ {1}-1}\\ &\times\ \left(\left(s\frac{E_{s}}{\sigma_{1}^{2}}+\left(\beta_{1}-\delta _{1}\right)\left(1+z\frac{E_{i}}{\sigma_{1}^{2}}\right)\right)y+C\left(\beta_ {1}-\delta_{1}\right)\right)^{-m_{1}}\end{array}  \tag{11}\]Let us now define the following integral from Equation (9) by using \(I_{1}\):  
\[I_{2}\triangleq\int\limits_{0}^{\infty}I_{1}f_{|h_{2}|^{2}}(y)dy \tag{12}\]  
We substitute Equations (4) and (11) in Equation (12) and reach at:  
\[ \begin{array}{rcl}I_{2}&=&\int\limits_{0}^{\infty}\lambda_{2}y^{m_{2}-1}e^{- \epsilon_{2}y}\left(\alpha_{1}\left(1+z\frac{E_{i}}{\sigma_{1}^{2}}\right)y+ \alpha_{1}C\right)\left(\left(s\frac{E_{i}}{\sigma_{1}^{2}}+\beta_{1}\left(1+ z\frac{E_{i}}{\sigma_{1}^{2}}\right)\right)y+C\beta_{1}\right)^{m_{1}-1}\\ &&\times\ \left(\left(s\frac{E_{i}}{\sigma_{1}^{2}}+(\beta_{1}-\delta_{1})\left(1+ z\frac{E_{i}}{\sigma_{1}^{2}}\right)\right)y+C\left(\beta_{1}-\delta_{1}\right) \right)^{-m_{1}}dy\end{array}  \tag{13}\]  
By employing the method outlined in [6], Equation (13) becomes:  
\[ \begin{array}{rcl}I_{2}&=&\lambda_{2}\int\limits_{0}^{\infty}\left(\alpha_{ 1}\left(1+z\frac{E_{i}}{\sigma_{1}^{2}}\right)y+\alpha_{1}C\right)\left(\left( s\frac{E_{i}}{\sigma_{1}^{2}}+\beta_{1}\left(1+z\frac{E_{i}}{\sigma_{1}^{2}} \right)\right)y+C\beta_{1}\right)^{c_{1}}\\ &&\times\ \left(\left(s\frac{E_{i}}{\sigma_{1}^{2}}+(\beta_{1}-\delta_{1}) \left(1+z\frac{E_{i}}{\sigma_{1}^{2}}\right)\right)y+C\left(\beta_{1}-\delta_{ 1}\right)\right)^{-m_{1}}\\ &&\times\ \left(\left(s\frac{E_{i}}{\sigma_{1}^{2}}+\beta_{1}\left(1+z \frac{E_{i}}{\sigma_{1}^{2}}\right)\right)y+C\beta_{1}\right)^{e}y^{m_{2}-1}e ^{-\epsilon_{2}y}dy\end{array}  \tag{14}\]  
where \(c_{1}=\lfloor m_{1}\rfloor-1\) and \(e=m_{1}-\lfloor m_{1}\rfloor\) for \(m_{1}>1\); \(c_{1}=0\) and \(e=m_{1}-1\) for \(m_{1}\leq 1\); and \(\lfloor x\rfloor\) denotes the largest integer not greater than \(x\). By the use of binomial expansion in Equation (14), we rewrite Equation (14) as:  
\[ \begin{array}{rcl}I_{2}&=&\lambda_{2}\sum\limits_{l=0}^{c_{1}}\left(c_{l}^{ 1}\right)\left(s\frac{E_{i}}{\sigma_{1}^{2}}+\beta_{1}\left(1+z\frac{E_{i}}{ \sigma_{1}^{2}}\right)\right)^{l}\left(C\beta_{1}\right)^{c_{1}-1}\int\limits_ {0}^{\infty}\left(\alpha_{1}\left(1+z\frac{E_{i}}{\sigma_{1}^{2}}\right)y+ \alpha_{1}C\right)\\ &&\times\ \left(\left(s\frac{E_{i}}{\sigma_{1}^{2}}+(\beta_{1}-\delta_{1}) \left(1+z\frac{E_{i}}{\sigma_{1}^{2}}\right)\right)y+C\left(\beta_{1}-\delta_ {1}\right)\right)^{-m_{1}}\\ &&\times\ \left(\left(s\frac{E_{i}}{\sigma_{1}^{2}}+\beta_{1}\left(1+z \frac{E_{i}}{\sigma_{1}^{2}}\right)\right)y+C\beta_{1}\right)^{e}y^{m_{2}+l-1} e^{-\epsilon_{2}y}dy\end{array}