Optimizations of Multi-hop Cooperative Molecular Communication in Cylindrical Anomalous-Diffusive Channel

Abstract

In this paper, the optimizations of multi-hop cooperative molecular communication (CMC) system in cylindrical anomalous-diffusive channel in three-dimensional enviroment are investigated. First, we derive the performance of bit error probability (BEP) of CMC system under decode-and-forward relay strategy. Then for achieving minimum average BEP, the optimization variables are detection thresholds at cooperative nodes and destination node, and the corresponding optimization problem is formulated. Furthermore, we use conjugate gradient (CG) algorithm to solve this optimization problem to search optimal detection thresholds. The numerical results show the optimal detection thresholds can be obtained by CG algorithm, which has good convergence behaviors with fewer iterations to achieve minimized average BEP compared with gradient decent algorithm and Bisection method which are used in molecular communication.

1. Introduction

The field of molecular communication (MC) [1-2] attracts more attentions of many researchers because of its application prospects in biological environments and industrial fields [3-5]. Recently, more studies have been focused on MC with diffusion in an unbounded fluid environment [6-7].

However, MC in an unbounded environment has some limitations in the practical applications due to its channel characteristics, such as in the blood vessel of human body [8]. Then MC in bounded environment has attracted a lot of attentions of many researchers. The authors [9-10] analytically derived the concentration Green’s function (CGF) of MC system in biological cylindrical channel. Based on the Poiseuille flow and the Robin boundary condition, the authors [11] utilized a Markovian-based channel model for MC system to deduce the probability density function. In 2022, Dhok et al. [12] evaluated the performance of probability error of a cooperative molecular communication (CMC) system in cylindrical environment by using a fusion center.

Usually, the molecules propagate in the channel which follows the conventional Brownian motion. However, the scenario of the anomalous diffusion is more extensive comparing with conventional Brownian motion. In 2017, Mai et al. [13] proposed an algorithm for optimizing the network throughput in MC system with anomalous diffusion. In 2019, Trinh et al. [14] derived the first passage time of the molecules by using timing modulation. In 2020, Chouhan et al. [15] calculated the bit error probability (BEP) based on deriving the expression of first passage time density of MC system. In 2021, Trinh et al. [16] derived the observed molecules and the closed-form expressions of the bit error rate. In 2022, the formulation of the first hitting time density of anomalous-diffusive MC system was deduced in [17]. The authors in [18] considered the anomalous diffusion and derived the expression of the CGF of an underlay cognitive MC system.

As so far, some traditional algorithms are used to obtain optimal decision threshold in order to achieve minimum BEP of MC system. Tavakkoli et al. [19] minimized the BEP of two-hop MC system by using bisection algorithm to obtain optimal detection threshold. Chouhan et al. [20] implemented gradient descent (GD) optimization for finding the solutions of the optimization problem which are the values of optimization variables including optimal decision threshold of MC system. Cheng et al. implemented particle swarm optimization algorithm [21] and adaptive genetic algorithm [22] to solve different optimization problems in order to achieve minimum BEP of MC system.

However, the multi-hop CMC system in cylindrical anomalous-diffusive channel in three-dimensional (3D) environment has not been studied. On one hand, because of multiple relay nodes in the CMC system, there are multiple detection thresholds needed to be computed by using the maximum a posteriori (MAP) probability detection method with multiple times. However, all the detection thresholds at cooperative nodes and destination node can be obtained simultaneously by using the optimization algorithm. On the other hand, the detection thresholds are optimized by minimizing the average BEP of the CMC system. Therefore, the optimized detection thresholds can improve the reliability of communication of the CMC system. In order to obtain minimum BEP of this system, how to optimize the detection thresholds at cooperative nodes and destination node is a challenge study. But the existing traditional algorithms needs more iterations. It is important to develop more efficient offline optimization techniques. In this paper, we use the conjugate gradient (CG) algorithm to solve the optimization problem. The main contributions of our paper are concluded as follows:

(1) The decode-and-forward (DF) relaying strategy is implemented at cooperative nodes for the multi-hop CMC system. Then the formulation of average BEP is derived.

(2) We set up an optimization problem for achieving minimum average BEP with optimization variables which are detection thresholds at cooperative nodes and destination node. Then CG algorithm is adopted to solve this optimization problem for searching the optimal detection thresholds.

(3) The numerical results have revealed that CG algorithm is more efficient in finding the optimal detection thresholds at cooperative nodes and the destination node with fewer iterations compared with GD algorithm and Bisection method.

The rest of our paper is organized as follows. Section 2 presents the CMC system model in 3D environment. In Section 3, the detection thresholds are optimized by CG algorithm. In Section 4, the performances of this CMC system with optimized detection thresholds are evaluated. Section 5 summarizes the paper.

2. The Multi-hop CMC System

The system model of multi-hop CMC in cylindrical channel with anomalous diffusion and drift in 3D enviroment is shown in Fig. 1. This cylindrical channel with the radius ρ_cy is a semi-infinite cylinder with boundary. This system is composed of one source node S, K cooperative nodes CN₁, CN₂, …, CN_K and one destination node D, which have fixed locations denoted by (ρ_S, ϕ_S, z_S), (ρ_CN1, ϕ_CN1, z_CN1), (ρ_CN2, ϕ_CN2, z_CN2), … , (ρ_CNK, ϕ_CNK, z_CNK), (ρ_D, ϕ_D, z_D), respectively. Here, (ρ, ϕ, z) describes the radial, azimuthal and axial coordinates of each node, respectively. In addition, we have 0 ≤ ρ ≤ ρ_cy, 0 ≤ ϕ ≤ 2π . In the axial direction, it extends from origin to infinity, then the range of z is 0 ≤ z <∞.

Fig. 1. The multi-hop CMC system.

We assume that the transmission time is time-slotted. Each time slot duration is T_s. The modulation method of ON/OFF keying [23] is adopted. The node S and K cooperative nodes instantaneously release molecules at the beginning of each time slot for transmitting bit 1, while releasing no molecules represents the transmission of bit 0. We also suppose that source node S, K cooperative nodes and node D can keep perfect synchronization in time [24-25].

After the transmission process starts, the DF relaying strategy is adopted at each relay node. Different types of molecules are used in each hop. When time slot j begins, the source node S releases molecules with type A₁ for transmitting bit 1. For the DF relaying, when time slot j ends, the cooperative node CN₁ receives the transmitted bit, which is forwarded to its adjacent node CN₂ when next time slot begins. Each node CN_k (k=1, 2, ... , K) can detect type A_k molecules from node CN_k-1 at the end of time slot (j+k-1). Then it releases A_k+1 molecule types to forward the decoded bit to node CN_k+1. CN₀ and CN_K+1 represent node S node D, respectively. The cooperative nodes and node D are the passive spherical receivers with the radius r_CNk (k=1, 2, ... , K) and r_D, respectively. In each time slot, they can count the number of molecules to make the decision that the received bit is 1 or 0 by comparing with the detection thresholds.

After node S and cooperative nodes release molecules, these molecules diffuse in cylindrical anomalous channel. This process can be modeled based on the type of diffusion phenomenon. D_α(t) represents the instantaneous diffusion coefficient which is computed by [26]

D_α(t) = αD_pt^α−1,       (1)

where D_p represents the diffusion coefficient of molecules. α∈[0,2] is the diffusion exponent which is defined by the normal diffusion with α = 1, sub-diffusion with α∈[0,1) and super-diffusion with α∈(1,2].

The CGF under the Robin’s boundary condition shows the concentration of molecules at node S at time t with location (ρ_S, ϕ_S, z_S) under given initial time t₀, which is denoted by C_S(t; t₀) [9]

\(\begin{align}\begin{aligned} C_{\mathrm{S}}\left(t ; t_{0}\right)= & \frac{\exp \left(-\frac{\left(z-z_{\mathrm{S}}-v\left(t-t_{0}\right)\right)^{2}}{4\left(t-t_{0}\right)^{\alpha} D_{p}}-\xi\left(t-t_{0}\right)\right)}{\sqrt{4 \pi\left(t-t_{0}\right)^{\alpha} D_{p}}} \\ & \times\left(\sum_{n=0}^{\infty} \sum_{m=1}^{\infty} Q_{n m} \cos \left(n\left(\phi-\phi_{\mathrm{S}}\right)\right) J_{n}\left(\lambda_{n m} \rho\right) \times \exp \left(-D_{p} \lambda_{n m}^{2}\left(t-t_{0}\right)^{\alpha}\right) \delta\left(t-t_{0}\right)\right),\end{aligned}\end{align}\),       (2)

where z_S is the coordinate value along z axis. v is the drift velocity. ξ is the degradation constant. ρ_S and ρ_cy are the radii of node S and the cylindrical channel, respectively. t₀ is the initial time instant. \(\begin{align}Q_{n m}=\frac{L_{n} J_{n}\left(\lambda_{n m} \rho_{\mathrm{s}}\right)}{N_{n m}}, N_{n m}=\frac{\rho_{c y}^{2}}{2}\left(J_{n}^{2}\left(\lambda_{n m} \rho_{c y}\right)-J_{n-1}\left(\lambda_{n m} \rho_{c y}\right) J_{n+1}\left(\lambda_{n m} \rho_{c y}\right)\right)\end{align}\).

J_n and λ_nm are the n-th order Bessel function of the first kind and the m-th eigenvalues, respectively. J_n-1 and J_n+1 are defined as J_n. n and m are integers. L_n and δ(t) are given as follows:

\(\begin{align}L_{n}=\left\{\begin{array}{l}\frac{1}{\pi}, n \geq 1, \\ \frac{1}{2 \pi}, n=0,\end{array}\right.\end{align}\),       (3)

\(\begin{align}\delta\left(t-t_{0}\right)=\left\{\begin{array}{l}1, t \geq t_{0}, \\ 0, t<t_{0} .\end{array}\right.\end{align}\).       (4)

The probability of the case that one released molecule is observed at node CN₁ at time t is obtained by [9]

p_(S,CN1)(t;t₀) = ∫∫∫_VCN1 C_S(t; t₀) ρ d ρdϕdz.       (5)

where V_CN1 is the volume of node CN₁ which a spherical region. Then (5) is approximated as p_(S,CN1)(t;t₀) = V_CN1C_S(t; t₀). When t₀ represents the l-th time slot in which one molecule emitted by node S, the probability of the case that this one molecule is observed in time slot j by node CN₁ is given by

p^lj_(S,CN1) = V_CN1C_S((j - l)T_s; t₀).       (6)

3. Derivation of average BEP of the CMC system

N_S[l] and x^l_S (1 ≤ l ≤ j) are used to denote the quantity of released molecules and the transmitted bit by node S in time slot l, respectively. In time slot j, considering the reception at node CN₁, N_(S,CN1)[j] is defined as the quantity of received molecules. Then we have

\(\begin{align}N_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}[j]=\sum_{l=1}^{j} N_{\mathrm{S}}[l] x_{\mathrm{S}}^{l} p_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{l j}+N_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{\mathrm{Noise}}\end{align}\),       (7)

where N^Noise_(S,CN1) is the noise generated for this link. It is modelled as a Normal distribution 𝒩(0,(σ^Noise_(S,CN1))²), which has mean 0 and variance (σ^Noise_(S,CN1))² [25].

When N_S[l] is large enough, according to the central limit theorem [25], N_(S,CN1)[j] is modelled as Gaussian approximation. Therefore, we have

N_(S,CN1)[j] ~ 𝒩(µ_(S,CN1)[j], σ²_(S,CN1)[j]),       (8)

where µ_(S,CN1)[j] is the mean and σ²_(S,CN1)[j] is the variance of N_(S,CN1)[j], which are obtained as follows:

\(\begin{align}\mu_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}[j]=\sum_{l=1}^{j} \pi_{1} N_{\mathrm{S}}[k] p_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{l j}\end{align}\),       (9)

\(\begin{align}\sigma_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{2}[j]=\sum_{l=1}^{j}\left[\pi_{1} N_{\mathrm{S}}[l] p_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{l j}\left(1-p_{\left(\mathrm{S}, \mathrm{CN}_{\mathrm{N}}\right)}^{l j}\right)+\left(N_{\mathrm{S}}[l]\right)^{2} \pi_{1} \pi_{0}\left(p_{\left(\mathrm{S}, \mathrm{CN}_{\mathrm{N}}\right)}^{k j}\right)^{)^{k}}\right]+\left(\sigma_{\left(\mathrm{S}, \mathrm{CN}_{\mathrm{N}}\right)}^{\mathrm{Naise}}\right)^{2}\end{align}\),       (10)

where π₁ and π₀ are the transmission probabilities of 1 and 0 by node S when each time slot begins, respectively. We have Pr(x_S[l] = 1) = π₁ = 0.5 and Pr(x_S[l] = 0) = π₀ = 0.5.

H₀ and H₁ indicate the cases that when time slot j begins, node S transmits bits 0 and 1, respectively. Then the corresponding Normal distributions under H₀ and H₁ are formulated by

H₀ : N_(S,CN1)[j] ~ 𝒩(μ⁰_(S,CN1)[j], (σ⁰_(S,CN1)[j])²),

H₁ : N_(S,CN1)[j] ~ 𝒩(μ¹_(S,CN1)[j], (σ¹_(S,CN1)[j])²),       (11)

where µ^w_(S,CN1)[j](w=0, 1) is the mean and (σ^w_(S,CN1)[j])² (w=0, 1) is the variance of N_(S,CN1)[j] under H_w, respectively. According to (9) and (10), µ^w_(S,CN1)[j] and (σ^w_(S,CN1)[j])² (w=0, 1) can be computed by

\(\begin{align}\begin{array}{l}\mu_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{0}[j]=\sum_{l=1}^{j-1} \pi_{1} N_{\mathrm{S}}[l] p_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{l(j-1)}, \\ \left.\left(\sigma_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{0}[j]\right)^{2}=\sum_{l=1}^{j-1}\left[\pi_{1} N_{\mathrm{S}}[l] p_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{l j}\right)\left(1-p_{(\mathrm{S}, \mathrm{CN})}^{l j}\right)+\left(N_{\mathrm{S}}[l]\right)^{2} \pi_{1} \pi_{0}\left(p_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{l j}\right)^{2}\right]+\left(\sigma_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{\mathrm{Noio})^{2}},\right. \\ \mu_{(\mathrm{S}, \mathrm{CN})}^{1}[j]=N_{\mathrm{S}}[j] p_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{j j}+\mu_{(\mathrm{S}, \mathrm{CN})}^{0}[j] \\ \left(\sigma_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{1}[j]\right)^{2}=N_{\mathrm{S}}[j] p_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{j j}\left(1-p_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{j j}\right)+\left(\sigma_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{0}[j]\right)^{2} \\\end{array}\end{align}\).       (12)

The detection threshold at node CN₁ is η_CN1. Then when each time slot ends, the detection rule at node CN₁ is written as

\(\begin{align}\hat{x}_{\mathrm{CN}_{1}}[j]=\left\{\begin{array}{l}1, \text { if } N_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}[j] \geq \eta_{\mathrm{CN}_{1}}, \\ 0, \text { if } N_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}[j]<\eta_{\mathrm{CN}_{1}},\end{array}\right.\end{align}\),       (13)

where \(\begin{align}\hat{x}_{\mathrm{CN}_{1}}^{j}\end{align}\) is denoted by the bit decoded by node CN₁ when time slot j ends. When \(\begin{align}x_{\mathrm{S}}^{j} \neq \hat{x}_{\mathrm{CN}_{1}}^{j}\end{align}\), an error occurs in this time slot. \(\begin{align}\operatorname{Pr}\left(\hat{x}_{\mathrm{CN}_{1}}^{j}=0 \mid x_{\mathrm{S}}^{j}=1\right)\end{align}\) and \(\begin{align}\operatorname{Pr}\left(\hat{x}_{\mathrm{CN}_{1}}^{j}=1 \mid x_{\mathrm{S}}^{j}=0\right)\end{align}\) represent the error probabilities of bits transmission of 1 and 0, respectively. Then we have

\(\begin{align}\begin{aligned} \operatorname{Pr}\left(\hat{x}_{\mathrm{CN}_{\mathrm{1}}}^{j}=0 \mid x_{\mathrm{S}}^{j}=1\right) & =\operatorname{Pr}\left(N_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}[j]<\eta_{\mathrm{CN}_{1}} \mid x_{\mathrm{S}}^{j}=1\right) \\ & =Q\left(\frac{\eta_{\mathrm{CN}_{1}}-\mu_{\left(\mathrm{S}, \mathrm{CN} N_{1}\right)}^{0}[j]}{\sigma_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{0}[j]}\right),\end{aligned}\end{align}\),       (14)

\(\begin{align}\begin{aligned} \operatorname{Pr}\left(\hat{x}_{\mathrm{CN}_{1}}^{j}=1 \mid x_{\mathrm{S}}^{j}=0\right) & =\operatorname{Pr}\left(N_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}[j] \geq \eta_{\mathrm{CN}_{1}} \mid x_{\mathrm{S}}^{j}=0\right) \\ & =Q\left(\frac{\eta_{\mathrm{CN}_{1}}-\mu_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{1}[j]}{\sigma_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}^{1}[j]}\right),\end{aligned}\end{align}\),       (15)

where the function Q(x) is denoted by \(\begin{align}Q(x)=\frac{1}{\sqrt{2 \pi}} \int_{x}^{\infty} e^{-\frac{u^{2}}{2}} d u\end{align}\) du. Based on formulas (14) and (15), Pe_(S,CN1)[j] which represents the BEP of one bit transmission in time slot j is written as

\(\begin{align}P e_{\left(\mathrm{S}, \mathrm{CN}_{1}\right)}[j]=\pi_{1} \operatorname{Pr}\left(\hat{x}_{\mathrm{CN}_{1}}^{j}=0 \mid x_{\mathrm{S}}^{j}=1\right)+\pi_{0} \operatorname{Pr}\left(\hat{x}_{\mathrm{CN}_{1}}^{j}=1 \mid x_{\mathrm{S}}^{j}=0\right)\end{align}\).        (16)

Considering the link CN_k → CN_k+1 in the (j+k)-th time slot, N_{(CNk,CNk+1 )}[j+k] is the number of received molecules by CN_k+1. The binary hypothesis testing problem based on N_{(CNk,CNk+1 )}[j+k] is established at node CN_k+1 as follows:

H₀ : N_(CNk,CNk+1)[j+k] ~ 𝒩(μ⁰_(CNk,CNk+1)[j], (σ⁰_(CNk,CNk+1)[j])²),

H₁ : N_(CNk,CNk+1)[j+k] ~ 𝒩(μ¹_(CNk,CNk+1)[j], (σ¹_(CNk,CNk+1)[j])²),       (17)

where µ^w_(CNk,CNk+1)[j] (w=0, 1) is the mean and (σ^w_(CNk,CNk+1)[j])²(w=0, 1) is the variance of N_(CNk,CNk+1)[j+k] under H_w, respectively.

Let x^j+k_CNk = 0 and x^j+k_CNk = 1 represent the bits transmission of 0 and 1 by CN_k at the time slot (j+k), respectively. \(\begin{align}\hat{x}_{\mathrm{CN}_{k+1}}^{j+k}\end{align}\) is 1 and 0 which are decoded by node CN_k+1 when time slot (j+k) ends, respectively. Thus the corresponding error probabilities for the link CN_k → CN_k+1 are defined as \(\begin{align}\operatorname{Pr}\left(\hat{x}_{\mathrm{CN}_{k+1}}^{j+k}=1 \mid x_{\mathrm{CN}_{k}}^{j+k}=0\right)\end{align}\) and \(\begin{align}\operatorname{Pr}\left(\hat{x}_{\mathrm{CN}_{k+1}}^{j+k}=0 \mid x_{\mathrm{CN}_{k}}^{j+k}=1\right)\end{align}\), respectively, which are computed according to (17) as follows:

\(\begin{align}\operatorname{Pr}\left(\hat{x}_{\mathrm{CN}_{k+1}}^{j+k}=1 \mid x_{\mathrm{CN}_{k}}^{j+k}=0\right)=\operatorname{Pr}\left(N_{\left(\mathrm{CN}_{k}, \mathrm{CN}_{k+1}\right)}[j+k] \geq \eta_{\mathrm{CN}_{k+1}} \mid x_{\mathrm{CN}_{k}}^{j+k}=0\right)\end{align}\),       (18)

\(\begin{align}\operatorname{Pr}\left(\hat{x}_{\mathrm{CN}_{k+1}}^{j+k}=0 \mid x_{\mathrm{CN}_{k}}^{j+k}=1\right)=\operatorname{Pr}\left(N_{\left(\mathrm{CN}_{k}, \mathrm{CN}_{k+1}\right)}[j+k]<\eta_{\mathrm{CN}_{k+1}} \mid x_{\mathrm{CN}_{k}}^{j+k}=1\right)\end{align}\),       (19)

where η_CNk+1 is detection threshold at node CN_k+1.

After node CN_k+1 decodes the received bit, it forwards \(\begin{align}\hat{x}_{\mathrm{CN}_{k+1}}^{j+k}\end{align}\) to the next relay node CN_k+2. The transmitted bit by CN_k+2 at the (j+k+1)-th time slot is \(\begin{align}\hat{x}_{\mathrm{CN}_{k+1}}^{j+k+1}\end{align}\). Assume the BEP of bit 0 or 1 transmission in time slot j from S to CN_k are Pe⁰_(S,CNk)[j] and Pe¹_(S,CNk)[j], respectively. Then considering the link S → CN_k+1, the BEP of one bit transmission in time slot j are computed by

\(\begin{align}\begin{array}{l}P e_{\left(\mathrm{S}, \mathrm{CN}_{k+1}\right)}^{0}[j]=P e_{\left(\mathrm{S}, \mathrm{CN}_{k}\right)}^{0}[j] \times \operatorname{Pr}\left[\hat{x}_{\mathrm{CN}_{k+1}}^{j+k}=1 \mid x_{\mathrm{CN}_{k}}^{j+k}=1\right] \\ +\left(1-P e_{\left(\mathrm{S}, \mathrm{CN}_{k}\right)}^{0}[j]\right) \times \operatorname{Pr}\left[\hat{x}_{\mathrm{CN}_{k+1}}^{j+k}=1 \mid x_{\mathrm{CN}_{k}}^{j+k}=0\right], \\\end{array}\end{align}\),       (20)

\(\begin{align}\begin{array}{l}P e_{\left(\mathrm{S}, \mathrm{CN}_{k+1}\right)}^{1}[j]=P e_{\left(\mathrm{S}, \mathrm{CN}_{k}\right)}^{1}[j] \times \operatorname{Pr}\left[\hat{x}_{\mathrm{CN}_{k+1}}^{j+k}=0 \mid x_{\mathrm{CN}_{k}}^{j+k}=0\right] \\ +\left(1-P e_{\left(S, \mathrm{CN}_{k}\right)}^{1}[j]\right) \times \operatorname{Pr}\left[\hat{x}_{\mathrm{CN}_{k+1}}^{j+k}=0 \mid x_{\mathrm{CN}_{k}}^{j+k}=1\right] . \\\end{array}\end{align}\).       (21)

When considering the multi-hop CMC system with k=K in (20) and (21), the BEP of one bit in the j-th time slot from node S denoted by Pe_(S,D)[j] is computed by

Pe_(S,D)[j] = π₁Pe¹_(S,CNK+1)[j] + π₀Pe⁰_(S,CNK+1)[j].       (22)

4. Optimization of detection thresholds of multi-hop CMC system

The optimization problem of detection thresholds of the multi-hop CMC system in cylindrical channel is expressed as follows:

\(\begin{align}\min _{\eta_{\mathrm{CN}_{1}}, \eta_{\mathrm{CN}_{2}}, \ldots, \eta_{\mathrm{CN}_{K}}, \eta_{\mathrm{D}}} P e_{(\mathrm{S}, \mathrm{D})}[j]\end{align}\),       (23)

where η_CN1, η_CN2, …, η_CNk, are the detection thresholds of nodes CN₁, CN₂, … , CN_K, respectively. η_D is the detection threshold at destination node D.

η_u is a vector of detection thresholds at the u-th iteration which is represented by η_u = [η^u_CN1, η^u_CN2, …, η^u_CNK, η^u_D]. ∇Pe_(S,D)(η_u) is the gradient of Pe_(S,D)[j] with η_u which is computed by

\(\begin{align}\nabla P e_{(S, D)}[j]\left(\eta_{u}\right)=\left[\frac{\partial P e_{(S, D)}[j]\left(\eta_{u}\right)}{\partial \eta_{\left(N_{1}\right.}}, \frac{\partial P e_{(S, D)}[j]\left(\eta_{u}\right)}{\partial \eta_{\mathrm{CN}_{2}}}, \ldots, \frac{\partial P e_{(S, D)}[j]\left(\eta_{u}\right)}{\partial \eta_{\mathrm{CN}_{K}}}, \frac{\partial P e_{(S, D)}[j]\left(\eta_{u}\right)}{\partial \eta_{\mathrm{D}}}\right]\end{align}\),       (24)

where \(\begin{align}\frac{\partial P e_{(S, \mathrm{D})}[j]\left(\eta_{u}\right)}{\partial \eta_{\mathrm{CN}_{1}}}, \frac{\partial P e_{(\mathrm{S}, \mathrm{D})}[j]\left(\eta_{u}\right)}{\partial \eta_{\mathrm{CN}_{2}}}\end{align}\) , ..., \(\begin{align}\frac{\partial P e_{(\mathrm{S}, \mathrm{D})}[j]\left(\eta_{u}\right)}{\partial \eta_{\mathrm{CN}_{K}}}, \frac{\partial P e_{(\mathrm{S}, \mathrm{D})}[j]\left(\eta_{u}\right)}{\partial \eta_{\mathrm{D}}}\end{align}\) are the first derivative of Pe_(S,D)[j](η_u) with respect to η_CN1, η_CN2, …, η_CNK, η_D, respectively. The CG algorithm is described as follows:

Algorithm 1. CG Algorithm for Optimizing Detection Thresholds in CMC System

It is noted that when ||∇Pe_(S,D)[j](η₀)|| > λ(λ = 10^-6) and the corresponding Pe_(S,D)[j](η₀) at η₀ is smaller than 0.5, this algorithm can converge fast, and it performs exceptionally well. Then we can get optimized results of detection thresholds at cooperative nodes and destination node. When ||∇Pe_(S,D)[j](η₀)|| ≤ λ and the corresponding Pe_(S,D)[j](η₀) at η₀ is equal to 0.5, this algorithm is difficult to achieve convergence. Under such a case, this algorithm will face challenges. In order to solve this problem, we search for a better η'₀ at the neighborhood of η₀ which can satisfy the condition that ||∇Pe_(S,D)[j](η₀)|| > λ. Then η₀ is updated as η₀ = η'₀.

5. Numerical results

The numerical results are shown in this section and the default numerical parameters are given in Table 1.

Table 1. The numerical parameters

η_CN1 and η_D are the thresholds at node CN₁ and node D, respectively. Fig. 2 shows the average BEP is changing with η_CN1 and η_D when there is one cooperation node CN₁. When η_CN1 takes the same value, the average BEP decreases first and then increases with the increasing value of η_D. The change trend under fixed value of η_D is the same as that under fixed value of η_CN1. We can see that it is a convex problem when two detection thresholds with η_CN1 and η_D are optimized. Therefore, CG algorithm can be adopted to solve this optimization problem efficiently. The value of the lowest point with red dots is represented by (4800,7150, 0.0762), which means that the optimized results are η_CN1 =4800, η_D =7150 and the corresponding average BEP is 0.0762.

Fig. 2. The average BEP vs η_CN1 and η_D.

In Fig. 3, we give the convergence analysis of CG, GD and Bisection algorithms. Under these three algorithms, when the number of iterations is increasing, the performances of average BEP are both decreasing. Finally, they can achieve convergence. We can see that CG algorithm converges fastest than GD algorithm and Bisection algorithm. We can see that for the same parameter setting, the number of iterations of CG algorithm, GD algorithm and Bisection algorithm are 5, 11 and 7, respectively. Therefore, we use CG algorithm to solve this optimization problem with fewer iterations and obtain optimal detection thresholds and minimum average BEP.

Fig. 3. The comparison results of convergence with CG, GD and Bisection algorithms.

When α and ρ_cy take different values, the convergence results are shown in Fig. 4(a) and Fig. 4(b), respectively. Here α is the diffusion exponent which has corresponding ranges. The CG algorithm has good convergence which shows that the average BEP can converge with fewer iterations. In Fig. 4(a), when α = 0.5 and α =1.5, the corresponding number of iterations are 6 and 5, respectively. In Fig. 4(b), when ρ_cy = 5μm and ρ_cy = 7μm, the average BEP achieves convergence within 3 iterations and 5 iterations, respectively. Second, when the values of α are smaller and the values of ρ_cy are larger, the average BEP are both larger. This result is explained by the facts: on one hand, smaller values of α results in slower diffusion. On the other hand, larger values of ρ_cy means the same number of molecules can diffuse in much larger space. These two cases both lead to a decrease of the number of molecules received by cooperative nodes CN₁, CN₂ and node D, finally there is a decrease in BEP for each link.

Fig. 4. Convergence analysis under different values of (a) α ; (b) ρ_cy.

The average BEP is changing with η_CN1 under different values of α and ρ_cy in Fig. 5. Different ranges in Fig. 5(a) including α < 1, α = 1 and α > 1 represent sub-diffusion, normal diffusion and super-diffusion, respectively. We can see that when α is increasing, the average BEP is found to decrease. This is explained as follows: larger value of α will accelerate the diffusion of molecules. Then the number of molecules arriving at nodes CN₁, CN₂ and node D will increase. Finally, this reduces the average BEP. In addition, the optimal thresholds at CN₁, CN₂ and node D are given at the lowest point with minimum value of average BEP. When α = 1.5, η_CN1 = 9100, η_CN2 = 10970, η_D = 11082 and the corresponding average BEP is 0.097. In Fig. 5(b), ρ_cy ={5μm, 7μm, 9μm}. The performance of average BEP is decreasing with η_CN1. When η_CN1 achieves some value, the average BEP achieves its minimum value and then increase. Moreover, the number of received molecules is larger when ρ_cy takes smaller value, and the corresponding detection thresholds are also larger to obtain minimum value of average BEP. For the case ρ_cy = 5μm, η_CN1 = 5100, η_CN2 = 9200, and η_D = 12228, the average BEP is larger than those under the cases ρ_cy = 7μm and ρ_cy = 9μm.

Fig. 5. The average BEP vs η_CN1 under different values of (a)α ; (b) ρ_cy.

We use d₀ to represent the distance between two adjacent nodes along z-axial direction. Then consider three-hop CMC system with two cooperation nodes CN₁ and CN₂, the transmission distance is 3d₀. Fig. 6 gives the result the average BEP is changing with detection thresholds when d₀ takes different values. The change trend of average BEP is similar as in Fig. 5. When d₀ is with larger value which results in lower receiving probability, the average BEP is larger and the detection thresholds at lowest point of average BEP are smaller. When d₀=5μm, the average BEP is 0.386 which is larger than that when d₀=3μm, and the corresponding values of detection thresholds at lowest point of BEP are η_CN1 = 3300, η_CN2 = 3707, η_D = 3721 which are smaller than those when d₀=3μm.

Fig. 6. The average BEP vs η_CN1 under different distances along z-axial direction for each hop.

In order to show the differences between DF relay strategy and amplify-and-forward (AF) relay strategy, we give the comparison results in Fig. 7. The parameters are set the same for the two relay strategies. Especially, for the AF relay strategy, the amplification factor is set as 5. According to Fig. 7, we can see that the average BEP under DF relay strategy decreases faster and is smaller than that under AF relay strategy for the same value of η_CN1. In particular, the detection thresholds at lowest point are the same for these two relay strategies. Therefore, for the CMC system, we choose DF relay strategy.

Fig. 7. The average BEP vs η_CN1 with DF and AF relay strategies.

z_(S,D) is used to represent the distance between nodes S and D along z-axial direction. In Fig. 8, the average BEP is decreasing with the value of K which is the number of cooperative nodes. When the total distance along z-axial direction z_(S,D) is fixed, the distance along z-axial direction for each hop is decreasing with increasing value of K. The receiving probabilities at each cooperative nodes CN₁, CN₂ and node D also increase. Finally, the average BEP is decreasing. For each same value of K, the average BEP of this CMC system under z_(S,D) = 20μm is the lowest than those under z_(S,D) = 10μm and z(S,D) = 15μm.

Fig. 8. The average BEP versus K under different distances between nodes S and D along z-axial direction.

6. Conclusion

This paper studied the optimizations of multi-hop CMC system in cylindrical anomalous-diffusive channel in a 3D environment. We have derived the average BEP which is optimization objective function. Then the optimization problem for minimizing the average BEP was effectively solved by CG algorithm. The numerical results have shown that CG algorithm can converge by using fewer iterations in contrast to GD algorithm and Bisection algorithm. In addition, we also have shown that some main parameters have impacts on the optimization results. The lager value of α , smaller value of ρ_cy and initial distance between two adjacent nodes along z-axial direction, then the smaller value of average BEP. In future work, we intend to explore other traditional and efficient optimzaiton methods for solving the optimization problem under the scenario with mobile nodes in CMC system. In additon, we plan to use data driven methods, such as the deep learning methods to optimize optimal detection thresholds under different system parameters for CMC system.