# Stability and performance analysis of Compound TCP with the Exponential-RED and the Drop-Tail queue policies

Stability and performance analysis of Compound TCP with the Exponential-RED and the Drop-Tail... Abstract The analysis of transport protocols, along with queue management policies, forms an important aspect of performance evaluation for the Internet. In this article, we study Compound TCP (C-TCP), the default TCP in the Windows operating system, along with the Exponential-RED (E-RED) queue policy and the widely used Drop-Tail queue policy. We consider queuing delay, link utilization and the stability of the queue size as the performance metrics. We first analyse the stability properties of a nonlinear model for C-TCP coupled with the E-RED queue policy. We observe that this system, in its current form, may be difficult to stabilize as the feedback delay gets large. Further, using an exogenous and non-dimensional parameter, we show that the system loses local stability via a Hopf bifurcation, which gives rise to limit cycles. Employing Poincaré normal forms and the center manifold theory, we outline an analytical framework to characterize the type of the Hopf bifurcation and to determine the orbital stability of the emerging limit cycles. Numerical examples, stability charts and bifurcation diagrams complement our analysis. We also conduct packet-level simulations, with E-RED and Drop-Tail, in small and large buffer-sizing regimes. With large buffers, E-RED can achieve small queue sizes compared with Drop-Tail. However, it is difficult to maintain the stability of the E-RED policy as the feedback delay gets large. On the other hand, with small buffers, E-RED offers no clear advantage over the simple Drop-Tail queue policy. Our work can offer insights for the design of queue policies that can ensure low latency and stability. 1. Introduction It is well known that the design of transport protocols, embedded in end systems and queue management policies, implemented in Internet routers, can impact network performance. The algorithms in transport protocols update their sending window based on feedback, about congestion in the network, from routers. However, feedback from routers is rarely instantaneous. Thus, with time-delayed feedback, stability and bifurcation phenomenon become important considerations for performance evaluation. Fluid models, for congestion control algorithms, are quite useful to investigate the stability and bifurcation properties of the underlying system. Such analysis can help understand trade-offs among the system parameters and can also guide design principles. For surveys on models for Internet congestion control, see Low & Srikant (2004) and Srikant (2004). Compound TCP (C-TCP) is widely implemented in the Internet today, as it is the default transport protocol in the Windows operating system (Tan et al., 2006). Moreover, recent studies (Yang et al., 2014) have shown that C-TCP contributes nearly $$25\%$$ of the total traffic in the Internet. Hence, its evaluation can help us understand the performance of transport protocols that are used in real networks. For an analytical evaluation, it is useful to have a model. The model we use for Compound is the fluid model presented in Raja & Raina (2012). An important metric for the performance analysis of communication networks is queuing delays, which are on the rise (Nichols & Jacobson, 2012). Excessive queuing delays are normally attributed to large, and unmanaged, router buffers (Gettys & Nichols, 2012; Wischik & McKeown, 2005). Such delays are undesirable as they can negatively impact performance, especially for delay sensitive applications. One way to address the concern of rising queuing delays is to employ effective queue policies. We also note that, in addition to the design of TCP and queue management, buffering dynamics do also play a crucial role in ensuring stable operation of communication networks (Sojoudi et al., 2014). Queuing delays could be significant in the scenario where the router buffers are sized as per the bandwidth-delay rule and would be negligible with small buffers (Raina et al., 2005). Drop-Tail is a simple queue policy, which drops all incoming packets when the router buffer is full. Another popular queue policy is random early detection (RED) (Floyd & Jacobson, 1993). However, the RED policy is not used in practice as there are no universally accepted guidelines on how to tune the various RED parameters. For some work on C-TCP, and TCP Reno, with Drop-Tail and the RED policy, see Han et al. (2005), Malangadan et al. (2013) and Patil et al. (2011). These articles studied the stability properties of the underlying system, and exhibited the emergence of limit cycles as system parameters vary, using both theory and packet-level simulations. Random exponential marking (REM) queue policy, one of the early proposals of active queue management (AQM), aims to achieve a target queue size and a target link utilization at the same time (Athuraliya et al., 2001). A detailed study on the stability and performance of REM policy along with C-TCP was carried out by Raina et al. (2016). Thus, for the performance evaluation of TCP and queue management policies, stability is an important consideration for performance evaluation. This area continues to be an active area of research (see Gong et al., 2014; Hong et al., 2015; Briscoe et al., 2016; Showail et al., 2016 and references therein). In this article, we study C-TCP with the Exponential-RED (E-RED) queue policy (Liu et al., 2005) which aims for high link utilization, while maintaining small queue sizes. As the model is nonlinear, it is useful to study the local stability and the local bifurcation properties. In our previous work (Prasad & Raina, 2014), we showed that C-TCP coupled with the E-RED queue policy undergoes a Hopf bifurcation (Hassard et al., 1981; Kuznetsov, 2004), when the condition for local stability is violated. We also provided, by employing Poincaré normal forms and the center manifold theory, an analytical framework to characterize the type of the Hopf bifurcation and the asymptotic orbital stability of limit cycles. We employ an exogenous and non-dimensional parameter as the bifurcation parameter for the Hopf bifurcation analysis. In our current work, we also derive a sufficient condition for local stability of C-TCP with E-RED, which suggests the need for the joint design of transport protocols and queue management strategies. We also conduct packet-level simulations using the network simulator (NS2) (Issariyakul & Hossain, 2011) to highlight the existence of limit cycles in the queue size dynamics. Numerical examples, stability charts and bifurcation diagrams also complement our analysis. The rest of this article is organized as follows. In Section 2, we outline models of C-TCP and the E-RED queue policy. In Section 3, we conduct a local stability analysis for the coupled system. In Section 4, we perform some packet-level simulations and finally summarize our contributions in Section 5. The detailed Hopf bifurcation analysis, for the nonlinear model of C-TCP and E-RED queue policy, is presented in the Appendix. 2. Models In this section, we first outline nonlinear, fluid models of C-TCP and the E-RED queue policy. We then conduct local stability analysis for the coupled system of C-TCP and the E-RED policy. The stability results of C-TCP with the Drop-Tail policy obtained in Raina et al. (2016) are summarized in Section 3, as we are interested in evaluating the performance of E-RED and the Drop-Tail policy with C-TCP. 2.1. C-TCP in a small buffer regime We consider the many flows, nonlinear, fluid model for C-TCP proposed in Raja & Raina (2012). In this model, $$\mathrm{w}(t)$$ denotes the average window size, $$i\big(\mathrm{w}(t)\big)$$ represents the increase of the window per positive acknowledgement and $$d\big(\mathrm{w}(t)\big)$$ represents the decrease of the current window size per negative acknowledgement. In small buffer regime, the congestion avoidance algorithm of C-TCP updates its current window size, $$\mathrm{w}(t)$$ as follows (Raja & Raina, 2012): $\mathrm{w}(t+1) = \left\{ \begin{array}{ll} \mathrm{w}(t)\Big(1+\alpha \mathrm{w}(t)^{k-1}\Big) & \mbox{if no loss},\\ \mathrm{w}(t)\big(1-\beta\big)& \mbox{if loss}, \end{array} \right.$ where $$\alpha$$, $$k$$, $$\beta$$ are protocol parameters. The functional forms for $$i\big(\mathrm{w}(t)\big)$$ and $$d\big(\mathrm{w}(t)\big)$$ are $$\label{funcforms} i\big(\mathrm{w}(t)\big) = \alpha \mathrm{w}(t)^{k-1} \quad \mbox{and} \quad \textit{d}\big(\mathrm{w}(t)\big) = \beta \mathrm{w}(t).$$ (2.1) A generalized nonlinear model for the congestion avoidance phase of transport protocols is given by Raina et al. (2005) $$\frac{\textit{d}\mathrm{w}(t)}{dt} = \biggl(\! i\big(\mathrm{w}(t)\big)\big( 1 - p(t-\tau) \big) - d\big(\mathrm{w}(t)\big) p\big(t-\tau\big) \!\biggr) \frac{\mathrm{w}(t-\tau)}{\tau},\label{eq:compoundwindow}$$ (2.2) with equilibrium \begin{equation*} i(\mathrm{w}_s)\big( 1 - p(\cdot) \big) = d(\mathrm{w}_s) p(\cdot), \end{equation*} where $$p(\cdot)$$ is the packet drop probability, $$\tau$$ is the round-trip time and $$\mathrm{w}_s$$ is the equilibrium window size. If the average window size of all flows is $$\mathrm{w}(t)$$, then the average rate at which packets are sent is assumed to be approximated by $$x(t) = \mathrm{w}(t)/\tau$$. Different forms of $$p(\cdot)$$ would be appropriate for different queue management policies. Substituting the functional forms (2.1) into (2.2) we get the following nonlinear model for C-TCP $$\frac{d\mathrm{w}(t)}{dt} = \bigg(\! \alpha \mathrm{w}(t)^{k-1}\big( 1 - p(t-\tau) \big) - \beta \mathrm{w}(t) p\big(t-\tau\big) \!\bigg) \frac{\mathrm{w}(t-\tau)}{\tau},\label{eq:compoundwindownew}$$ (2.3) where $$p(\cdot)$$ is a function of the arrival rate at the queue. A fluid model for small buffer Drop-Tail is given by Raina et al. (2005) \begin{align} p(y_l) = (y_l/C_l)^B, \notag \end{align} where $$y_l$$ is the total arrival rate, $$C_l$$ is the link service rate and $$B$$ is the router buffer size. This fluid model combined with (2.3) yields a model for C-TCP with small buffers using a Drop-Tail policy. 2.2. E-RED queue policy The algorithm for E-RED has been proposed to stabilize TCP Reno by slightly modifying the original RED algorithm (Liu et al., 2005). The packet marking/dropping probability $$p(t)$$ in E-RED queue policy is set to be an exponential function of virtual queue length whose capacity is slightly smaller than the link capacity. At the packet level, $$p(t)$$ can be calculated as (Liu et al., 2005): \begin{align}\label{equation:ered algorithm} p(t) = \begin{cases} \begin{array}{l} 0 \\ p_{\rm min} e^{\frac{k_l}{\gamma{}C_l}\big(b(t)-th_{\rm min}\big)} \\ 1\end{array} \begin{array}{l} \text{if} \,\, 0 \leq b(t) < th_{\rm min}, \\ \text{if} \,\, th_{\rm min} < b(t) < th_{\rm max}, \\ \text{if} \,\, b(t) \geq th_{\rm max}, \end{array} \end{cases} \end{align} (2.4) where $$b(t)$$ is virtual queue length, $$p_{\rm min}$$ is the minimum probability, $$b(t) = th_{\rm min}$$, $$k_l>0$$ is the link gain, $$0<\gamma{}\leq1$$ is the target link utilization, $$C_l$$ is link capacity and $$th_{\rm min}, th_{\rm max}$$ are queue length thresholds. 2.3. Nonlinear model for E-RED We now outline the model that describe the dynamics of E-RED queue policy (Liu et al., 2005). The E-RED queue policy aims to achieve high link utilization by dynamically adapting the packet loss probability at the queues. E-RED adapts the packet loss probability $$p(t)$$ based on the aggregate rate at the link. The link dynamics of this queue policy is given by Liu et al. (2005) $$\label{equation:ered} \frac{dp(t)}{dt}=k_l\,\frac{p(t)}{\gamma{}C}\Big(x(t)-\gamma{}C\Big),$$ (2.5) where $$x(t)=\mathrm{w}(t)/\tau$$ is the average arrival rate at the queue, and $$C$$ is the per-flow service rate of the bottleneck link. 3. Local stability analysis We analyse the nonlinear model for C-TCP (2.3) combined with the E-RED queue policy (2.5). First, we derive the Hopf condition for C-TCP with the E-RED queue policy. There are numerous parameters that could act as the bifurcation parameter in the coupled system of C-TCP and the E-RED queue policy. For example, the Compound parameter $$\alpha$$, the link gain $$k_l$$ or the feedback delay $$\tau$$. Instead, we use an exogenous and non-dimensional parameter $$\kappa$$, as the bifurcation parameter. Introducing this parameter into the system, we have the following nonlinear equations \begin{align} \frac{d\mathrm{w}(t)}{dt} &= \kappa{}\bigg( \alpha \mathrm{w}(t)^{k-1}\big( 1 - p(t-\tau) \big) - \beta \mathrm{w}(t) p\big(t-\tau\big) \bigg) \frac{\mathrm{w}(t-\tau)}{\tau}, \nonumber \\ \frac{dp(t)}{dt} &= \kappa{}\left( k_l\frac{p(t)}{\gamma{}C}\big(x(t)-\gamma{}C \big) \right)\!. \label{equation:kappa} \end{align} (3.1) Note that the parameter $$\kappa$$ does not affect the system equilibrium. We assume that the model parameters are chosen such that the system is at the edge of local stability. At this condition, the value of $$\kappa$$ is $$1$$. Now, instead of varying the model parameters, we gradually vary $$\kappa$$ beyond the edge of local stability. Thus varying $$\kappa$$, beyond $$1$$, will render the nonlinear system in a locally unstable state. Let $$(\mathrm{w}_s, p_s)$$ be the equilibrium of system (3.1), then at equilibrium, we get the following conditions \label{eqpoint} \begin{aligned} \mathrm{w}_s = \gamma{}C\tau \,\,\,\, \text{and} \,\,\,\, p_s = \frac{\alpha \mathrm{w}_s^{k-1}}{\alpha \mathrm{w}_s^{k-1} + \beta \mathrm{w}_s}. \end{aligned} (3.2) Let us consider the perturbations, $$u_1(t) = \mathrm{w}(t)-\mathrm{w}_s$$ and $$u_2(t) = p(t)-p_s$$. Linearizing system (3.1), about the equilibrium (3.2), we obtain \begin{align} \frac{du_1(t)}{dt} &= \kappa{}\bigg((k-2)(1-p_s)u_1(t)-\frac{\mathrm{w}_s}{p_s}u_2(t-\tau) \bigg)\frac{i(\mathrm{w}_s)}{\tau}, \notag \\ \frac{du_2(t)}{dt} &= \kappa{}\frac{k_lp_s }{\mathrm{w}_s}u_1(t). \label{lin1} \end{align} (3.3) Looking for exponential solutions, we get the following characteristic equation $$\label{ceq} \lambda^2+a\kappa{}\lambda+\kappa{}^2be^{-\lambda{}\tau}=0,$$ (3.4) where \begin{aligned} a =\frac{-i(\mathrm{w}_s)}{\tau}(k-2)(1-p_s) \,\, \text{and}\,\, b = \frac{i(\mathrm{w}_s)}{\tau}k_l. \end{aligned} Let $$\lambda$$ be a root of the characteristic equation, (3.4), and we determine the critical value $$\kappa = \kappa_c$$ at which $$\lambda$$ is purely imaginary. Substituting $$\lambda = j\omega_0$$ in the characteristic equation (3.4), we get $$\label{lsa1} -\omega_0{}^2+j\kappa{}_c a\omega_0+\kappa{}^2_cb\,\big(\cos{}\omega_0{}\tau-j\sin{}\omega_0{}\tau\big)=0.$$ (3.5) From (3.5), we have \begin{equation*} \begin{aligned} \cos{}\omega_0{}\tau =\frac{\omega_0{}^2}{\kappa{}^2_cb} \,\,\,\, \text{and} \,\,\,\, \sin{}\omega_0{}\tau =\frac{\omega_0 a{}}{\kappa{}_cb}. \end{aligned} \end{equation*} Let $$\tilde{\omega}_0=\omega_0/\kappa_c$$. Squaring and adding the above equations, we obtain \begin{align} &\tilde{\omega}_0 = \sqrt{\frac{-a^2+\sqrt{a^4+4b^2}}{2}}. \label{tilome} \end{align} (3.6) We now check whether the roots of (3.4) cross imaginary axis with a non-zero velocity i.e., $$\text{Re}\big(d\lambda/d\kappa\big)|_{\kappa_c} \neq 0$$. Differentiating (3.4) with respect to $$\kappa$$, we have $\frac{d\lambda}{d\kappa}=\frac{2\lambda^2+a\kappa\lambda}{\lambda^2\kappa\tau+\lambda(a\kappa^2\tau+2\kappa)+a\kappa^2}.$ Evaluating the above expression at $$\kappa_c$$, we get \begin{align*} \text{Re}\bigg(\frac{d\lambda}{d\kappa}\bigg)_{\kappa = \kappa_c} = \frac{\tau\big(2\omega_0^4+\kappa^2_c{} a^2 \omega_0^2 \big)}{ \kappa_c\big((\kappa_c{}a -\omega_0^2\tau)^2+\omega_0^2(\kappa_c{}a \tau{}+2)^2\big)}>0, \end{align*} which satisfies the transversality condition of the Hopf spectrum (Hassard et al., 1981). Further, the Hopf bifurcation occurs at $$\kappa{}_c = \frac{1}{\tilde{\omega}_0\tau}\tan^{-1}\left(\frac{a}{\tilde{\omega}_0}\right) \label{hopfcondition}.$$ (3.7) Hence, the nonlinear system (3.1) is locally stable for $$\kappa < \kappa_c$$, and undergoes a local Hopf bifurcation at $$\kappa= \kappa_c$$. We first make an observation on the choice of the Compound parameter $$k$$, whose default value is $$0.75$$. Observe that for $$k=1$$, the function $$i\big(\mathrm{w}(t)\big)$$ becomes independent of the window size $$\mathrm{w}(t)$$, and the equilibrium probability is $$p_s=\alpha/(\alpha+\beta{}\gamma{}C\tau)$$. Thus, for a high bandwidth-delay environment $$(1-p_s) \approx 1$$. In Fig. 1, we plot the local stability chart of system (3.1). This plot shows the relationship between the Compound parameter $$\alpha$$ and the link gain $$k_l$$, if local stability is to be maintained. Figure 2 shows the relationship between the link gain $$k_l$$ and the feedback delay $$\tau$$, if local stability is to be ensured. Fig. 1. View largeDownload slide Local stability chart of system (3.1), capturing the relationship between the Compound parameter $$\alpha$$ and the link gain $$k_l$$. The other system parameter values are $$\beta=0.5$$, $$k=1$$ and $$\gamma=1$$, and we assume $$(1-p_s)\approx1$$. The region to the left of the curves is stable. Fig. 1. View largeDownload slide Local stability chart of system (3.1), capturing the relationship between the Compound parameter $$\alpha$$ and the link gain $$k_l$$. The other system parameter values are $$\beta=0.5$$, $$k=1$$ and $$\gamma=1$$, and we assume $$(1-p_s)\approx1$$. The region to the left of the curves is stable. Fig. 2. View largeDownload slide Local stability chart of system (3.1) showing variation of link gain $$k_l$$ with respect to feedback delay $$\tau$$. The other system parameter values are $$\alpha=0.125$$, $$\beta=0.5$$, $$k=1$$, $$\gamma=1$$ and we assume $$(1-p_s)\approx1$$. The region under the Hopf condition curve is stable. Fig. 2. View largeDownload slide Local stability chart of system (3.1) showing variation of link gain $$k_l$$ with respect to feedback delay $$\tau$$. The other system parameter values are $$\alpha=0.125$$, $$\beta=0.5$$, $$k=1$$, $$\gamma=1$$ and we assume $$(1-p_s)\approx1$$. The region under the Hopf condition curve is stable. 3.1. Sufficient condition From the characteristic equation, (3.4) and for $$\kappa = 1$$, we obtain the loop transfer function of linear system (3.3) as \begin{align} L(j\omega) = \frac{b e^{-j\omega{}\tau}}{j\omega(a + j\omega )}. \end{align} (3.8) Let $$\omega_{gc}$$ denote the gain crossover frequency of $$L(j\omega)$$, which can be obtained by solving $$| L(j\omega_{gc})| = 1$$. By employing Nyquist stability criterion, the coupled system of C-TCP and E-RED is locally asymptotically stable, if phase margin defined as $$\text{PM} = \pi + \angle L(j\omega_{gc})$$ is greater than zero. \begin{align}\label{equation:phase margin} \text{PM} &= \pi - \pi/2 -\omega_{gc}\tau - \tan^{-1}(\omega_{gc}/a) > 0, \nonumber \\ & \pi/2 -\omega_{gc}\tau - \tan^{-1}(\omega_{gc}/a) > 0, \nonumber \\ & \omega_{gc}\tau + \tan^{-1}(\omega_{gc}/a) < \pi/2, \end{align} (3.9) where $$\quad \omega_{gc} = \sqrt{(-a^2 + \sqrt{a^4 + 4b^2})/{2}}$$. Note that the gain crossover frequency, $$\omega_{gc}$$, increases with increase in $$b$$, and if the value of $$b$$ gets large, then the stability condition given by (3.9) cannot be satisfied. Thus, we consider the values of $$b$$ to be smaller to obtain a sufficient condition for local stability of (3.1). For values of $$a,\,b$$ satisfying $$2b/a^2 << 1$$, $$\omega_{gc}$$ can be approximated to $$b/a$$. Now by substituting $$\omega_{gc} = b/a$$ in condition (3.9), we get $\frac{b\tau }{a}+ \tan^{-1}(b/a^2) < \pi/2.$ As $$(b/a^2) < 1$$, and using the expressions for $$a$$, $$b$$, we obtain the sufficient condition for stability as $$\label{equation:sufficient} \frac{k_l \tau }{(2-k)(1-p_s)} < \pi/4.$$ (3.10) For the choice of Compound parameter $$k=0$$, we get $k_l \tau \left(1 + \frac{\alpha}{\beta \mathrm{w}^2_s} \right) < \pi/2,$ where $$\mathrm{w}_s$$, $$p_s$$ are equilibrium window size and equilibrium packet-loss probability respectively. Sufficient condition (3.10) highlights that the choice of system parameter values is crucial in ensuring the stability of overall system (3.1). The stability analysis of C-TCP in a small buffer regime with Drop-Tail queues was conducted in Raina et al. (2016). The stability results of (i) C-TCP with E-RED and (ii) C-TCP with Drop-Tail are summarized in Table 1. Table 1 Local stability results of C-TCP with E-RED and Drop-Tail queue policies Queue policy Sufficient condition Hopf condition E-RED $$\displaystyle \frac{k_l \tau }{(2-k)(1-p_s)} < \pi/4$$ $$\displaystyle \kappa{}_c = \frac{1}{\tilde{\omega}_0\tau}\tan^{-1}\left(\frac{a}{\tilde{\omega}_0}\right)\!, \quad \tilde{\omega}_0 = \sqrt{\frac{-a^2+\sqrt{a^4+4b^2}}{2}}$$ Drop-Tail $$\displaystyle \alpha B \mathrm{w}^{k-1}_s < \frac{\pi}{2}$$ $$\displaystyle \alpha \mathrm{w}^{k-1}_s \sqrt{B^2 - \big(k-2\big)^2\big(1-p_s\big)^2} = \cos^{-1}\Big( \big(k-2\big)\big(1-p_s\big)\big/B \Big)$$ Queue policy Sufficient condition Hopf condition E-RED $$\displaystyle \frac{k_l \tau }{(2-k)(1-p_s)} < \pi/4$$ $$\displaystyle \kappa{}_c = \frac{1}{\tilde{\omega}_0\tau}\tan^{-1}\left(\frac{a}{\tilde{\omega}_0}\right)\!, \quad \tilde{\omega}_0 = \sqrt{\frac{-a^2+\sqrt{a^4+4b^2}}{2}}$$ Drop-Tail $$\displaystyle \alpha B \mathrm{w}^{k-1}_s < \frac{\pi}{2}$$ $$\displaystyle \alpha \mathrm{w}^{k-1}_s \sqrt{B^2 - \big(k-2\big)^2\big(1-p_s\big)^2} = \cos^{-1}\Big( \big(k-2\big)\big(1-p_s\big)\big/B \Big)$$ Table 1 Local stability results of C-TCP with E-RED and Drop-Tail queue policies Queue policy Sufficient condition Hopf condition E-RED $$\displaystyle \frac{k_l \tau }{(2-k)(1-p_s)} < \pi/4$$ $$\displaystyle \kappa{}_c = \frac{1}{\tilde{\omega}_0\tau}\tan^{-1}\left(\frac{a}{\tilde{\omega}_0}\right)\!, \quad \tilde{\omega}_0 = \sqrt{\frac{-a^2+\sqrt{a^4+4b^2}}{2}}$$ Drop-Tail $$\displaystyle \alpha B \mathrm{w}^{k-1}_s < \frac{\pi}{2}$$ $$\displaystyle \alpha \mathrm{w}^{k-1}_s \sqrt{B^2 - \big(k-2\big)^2\big(1-p_s\big)^2} = \cos^{-1}\Big( \big(k-2\big)\big(1-p_s\big)\big/B \Big)$$ Queue policy Sufficient condition Hopf condition E-RED $$\displaystyle \frac{k_l \tau }{(2-k)(1-p_s)} < \pi/4$$ $$\displaystyle \kappa{}_c = \frac{1}{\tilde{\omega}_0\tau}\tan^{-1}\left(\frac{a}{\tilde{\omega}_0}\right)\!, \quad \tilde{\omega}_0 = \sqrt{\frac{-a^2+\sqrt{a^4+4b^2}}{2}}$$ Drop-Tail $$\displaystyle \alpha B \mathrm{w}^{k-1}_s < \frac{\pi}{2}$$ $$\displaystyle \alpha \mathrm{w}^{k-1}_s \sqrt{B^2 - \big(k-2\big)^2\big(1-p_s\big)^2} = \cos^{-1}\Big( \big(k-2\big)\big(1-p_s\big)\big/B \Big)$$ 4. Packet-level simulations We conduct some packet-level simulations using the NS2 (Issariyakul & Hossain, 2011). The topology we consider is a single bottleneck link with a capacity of $$100$$ Mbps. We use the default parameter values of Compound (Tan et al., 2006): $$\alpha = 0.125$$, $$\beta = 0.5$$ and $$k = 0.75$$ and implemented the E-RED policy, as described in Liu et al. (2005), in NS2. The values of E-RED parameters chosen are: $$p_{\rm min} = 0.000001$$, $$th_{\rm min} = 1/5 \times$$ buffer size, $$th_{\rm max} = 3/5 \times$$ buffer size. The value of link gain parameter is chosen to satisfy the condition $$k_l*\tau_{\rm max} = 1$$. With the above choice of parameter values, we can determine the maximum drop probability, $$p_{\rm max}$$ from the definition of $$p(t)$$ (2.4). We consider two different traffic mixes: (i) long-lived flows and (ii) mix of long-lived and short-lived flows. The simulation run time is 225 s, and all the long-lived flows start randomly within the first 10 s. We illustrate the traces of the various quantities of interest over the last 25 s as the transient effects would settle down by then, and the system will be in steady state. We first conduct packet-level simulations of C-TCP, with Drop-Tail and E-RED, in a large and a small buffer regime in Section 4.1. The bottleneck link capacity is $$100$$ Mbps, and the large buffer regime (bandwidth-delay rule) translates into $$2084$$ packets. This is because a value of $$250$$ ms is used in practice, irrespective of the actual delays of the TCP flows (Wischik & McKeown, 2005). For the small buffer regime, parameters like the link capacity, the feedback delay, and the number of users do not influence the choice of buffer size. In this regime, we consider buffer sizes to be in the range of $$15$$–$$100$$ packets, as guided by the stability analysis (Raina et al., 2016). The traffic consists of $$60$$ long-lived C-TCP flows, each with a $$2$$ Mbps link feeding into the bottleneck. In Section 4.2, we explore the ability of E-RED to achieve desired link utilization and maintain small queue sizes with C-TCP. We consider the large buffer regime as it is currently widely deployed. In this set of simulations, the traffic mix is a combination of FTP, UDP and HTTP flows. The traffic consists of $$55$$ FTP flows of C-TCP, each with a $$2$$ Mbps link, and $$5$$ UDP flows contributing a total of $$5$$ Mbps. For short-lived flows, we use the PackMime-HTTP package in NS2. We generate $$10$$ new HTTP connections every second, on an average. The HTTP traffic model is described, in detail, in Cao et al. (2004). For E-RED we choose target link utilization of $$98 \%$$. The average round-trip time of the FTP flows present is $$10$$ ms and $$200$$ ms. For completeness, we also conduct simulations with the Drop-Tail policy. For all our simulations, we consider a fixed packet size of $$1500$$ bytes. 4.1. C-TCP, E-RED and Drop-Tail, Small and Large buffers With small buffers, as router buffer sizes vary from $$15$$ to $$100$$ packets, with both Drop-Tail and E-RED we note the emergence of stable limit cycles in the queue size, see Fig. 3(a) (buffer size = $$15$$ packets) and Fig. 3(b) (buffer size = $$100$$ packets). The analytical results of Compound with Drop-Tail predicted the loss of local stability, via a Hopf, as buffer sizes got larger, see Table 1. In case of E-RED queue policy, with buffer size of $$15$$ packets, observe that the virtual queue size oscillates as the round-trip time gets large. The analysis, for C-TCP with E-RED, carried out in Section 3 predicted the emergence of limit cycles as feedback delay increases. For the buffer size of $$15$$ packets, the transport layer protocol acts to control the distribution of the queue size, and it is possible that E-RED effectively degenerates into Drop-Tail. Fig. 3. View largeDownload slide Packet-level simulations of C-TCP with (i) Drop-Tail and (ii) E-RED queue policies employed individually at the bottleneck link. For the E-RED policy, we have plotted the real queue length (solid line) and the virtual queue length (dotted line). (a) With 15 packet buffers, the queue sizes with Drop-Tail and with E-RED policy do not exhibit limit cycles. However, the loss rate of C-TCP and E-RED is larger than that of C-TCP and Drop-Tail. (b)With 100 packet buffers, with Drop-Tail, note the emergence of limit cycles in queue size with larger feedback delay, which can hurt the link utilization. With E-RED queue policy, we notice limit cycles in virtual queue size (dotted line) for large RTT, which results in the loss of link utilization. (c)With large buffers, Drop-Tail achieves full utilization, but at the cost of large queuing delays. With ERED, the oscillations in queue size (real and virtual) for RTT = 200 ms are very clear, these oscillations result in loss of link utilization. Fig. 3. View largeDownload slide Packet-level simulations of C-TCP with (i) Drop-Tail and (ii) E-RED queue policies employed individually at the bottleneck link. For the E-RED policy, we have plotted the real queue length (solid line) and the virtual queue length (dotted line). (a) With 15 packet buffers, the queue sizes with Drop-Tail and with E-RED policy do not exhibit limit cycles. However, the loss rate of C-TCP and E-RED is larger than that of C-TCP and Drop-Tail. (b)With 100 packet buffers, with Drop-Tail, note the emergence of limit cycles in queue size with larger feedback delay, which can hurt the link utilization. With E-RED queue policy, we notice limit cycles in virtual queue size (dotted line) for large RTT, which results in the loss of link utilization. (c)With large buffers, Drop-Tail achieves full utilization, but at the cost of large queuing delays. With ERED, the oscillations in queue size (real and virtual) for RTT = 200 ms are very clear, these oscillations result in loss of link utilization. With large buffers and large RTTs, E-RED exhibits oscillations in queue size which leads to loss in link utilization, see Fig. 3(c). Drop-Tail, with large buffers, is expected to have full link utilization but this comes at the cost of enormous queuing delays. 4.2. C-TCP, E-RED and large buffers: short-lived and long-lived flows Large buffers are widely deployed in the Internet routers, and a good queue management policy should be able to ensure low latency, stability of queue sizes and offer a decent link utilization. We consider a traffic mix which consists of combination of $$55$$ FTP flows of C-TCP, each with a $$2$$ Mbps link, $$5$$ UDP flows (contributing a total of $$5$$ Mbps) and an average of $$10$$ new HTTP flows per second. With smaller round-trip times for the FTP flows, we observe that E-RED achieves the desired link utilization. With larger round-trip times, the queue sizes are small but the target link utilization is not achieved. For the traffic mix of FTP flows of C-TCP, UDP and HTTP flows, E-RED policy fails to achieve target link utilization as the feedback delay gets large, see Fig. 4(b). Thus, while latency could be reduced, the link utilization would be adversely impacted. Note that the nonlinear model of E-RED, analysed in Section 3, did alert us to the possibility of large delays destabilizing the system. Drop-Tail, with the same traffic mix, with large buffers has high utilization but excessive queuing delays, see Fig. 4(a). Fig. 4. View largeDownload slide (a) Drop-Tail and (b) E-RED queue policies with large buffers, with a traffic mix of FTP, UDP and HTTP flows. For E-RED policy, we have plotted the real queue length (solid line) and the virtual queue length (dotted line). Observe that the Drop-Tail policy achieves full link utilization but at the cost of large queuing delays. In case of E-RED policy, loss in link utilization occurs as the feedback delay increases. Thus, with higher RTTs, E-RED queue policy fails to achieve target link utilization. Fig. 4. View largeDownload slide (a) Drop-Tail and (b) E-RED queue policies with large buffers, with a traffic mix of FTP, UDP and HTTP flows. For E-RED policy, we have plotted the real queue length (solid line) and the virtual queue length (dotted line). Observe that the Drop-Tail policy achieves full link utilization but at the cost of large queuing delays. In case of E-RED policy, loss in link utilization occurs as the feedback delay increases. Thus, with higher RTTs, E-RED queue policy fails to achieve target link utilization. Clearly, the combination of C-TCP, queue management schemes and buffer size impacts network performance. The size of router buffers plays a central role in the end-to-end performance, as variations in buffer sizes influence both latency and queue stability. From the analysis and the packet-level simulations conducted so far, the E-RED queue policy does not seem to offer any clear advantage over Drop-Tail in terms of performance and stability of the coupled system (3.1). One can perform a variety of simulations considering different topologies and a wider range of traffic mixes, but it is out of the scope of this article. 5. Contributions Queuing delays in the Internet are on the rise (Nichols & Jacobson, 2012), and large queuing delays would adversely impact the Quality of Service of delay-sensitive web applications. One way to address the issue of queuing delays is to deploy effective AQM policies in Internet routers (Gettys & Nichols, 2012; Nichols & Jacobson, 2012). There have been many proposals for AQMs in the literature; however, there is no consensus on the optimal choice of an AQM policy. For the evaluation of queue management policies, it is also important to consider TCP variants that are currently implemented in the Internet. In this work, we focused on C-TCP, the default TCP in the Windows operating system. Among various queue management policies, we considered E-RED as it aims to offer high link utilization and achieve small queue sizes. We conducted a performance evaluation of the E-RED queue policy along with C-TCP flows operating over two different buffer-sizing regimes; Large (bandwidth-delay product) buffers and small buffers. Our study was based on the stability analysis of fluid models and packet-level simulations using the NS2 (Issariyakul & Hossain, 2011). We considered queuing delay, link utilization and the stability of the queue size as the performance metrics. We first conducted a local stability analysis of C-TCP along with the E-RED queue policy, and derived a sufficient, and the necessary and sufficient conditions for local stability. From the stability conditions, we highlighted that there exists a relation between C-TCP parameters, E-RED parameters and the feedback delay to ensure system stability. Such conditions on the parameters of TCP and AQM highlight the need for the joint design of transport protocols and queue management strategies. Using an exogenous and non-dimensional parameter, we also showed that the underlying system would undergo a Hopf type bifurcation. The non-dimensional parameter is also used to investigate the nature of the local instability. Further, we outlined a detailed local Hopf bifurcation analysis, using Poincaré normal forms and the center manifold theory (Hassard et al., 1981), to characterize the type of the Hopf bifurcation and to determine the orbital stability of the emerging limit cycles. Packet-level simulations conducted for different buffer-sizing regimes highlighted the existence of limit cycles in the queue size dynamics. These limit cycles in the queue size are undesirable and could lead to a loss in the link utilization. In case of large buffers, Drop-Tail can achieve full link utilization but at the cost of large queuing delays. With the same network setting, E-RED can maintain small queues but can fail to achieve the desired link utilization as the feedback delay gets large. Moreover, with E-RED, the stability of the overall system is also difficult to maintain with large delays. With small buffers and Drop-Tail, the system has low latency and stability can be ensured; however, the buffer size has to be dimensioned carefully. In such a regime, the E-RED policy does not offer any clear advantage over the simple Drop-Tail policy. The insights developed in this article could help in the design of queue management policies to ensure stability and low latency. 5.1. Avenues for further research To address the issues related to latency in the Internet, it is essential to analyse currently implemented transport protocols, of end systems, and queue management policies of the Internet routers. Compound TCP (Tan et al., 2006), the default protocol in Windows operating system, and Cubic TCP (Ha et al., 2008), protocol in Linux operating system, are widely used in the present Internet. To the best of our knowledge, fluid models for Cubic TCP have not been proposed in the literature. As fluid models are often used to conduct control-theoretic analysis, it would be useful to devise models for Cubic TCP, which may then aid performance evaluation. For queue management, it would be important to compare the performance of E-RED and Drop-Tail with other queue policies as well, especially in different buffer-sizing regimes. Adding to analysis, and packet-level simulations, experiments over a hardware platform would certainly aid the performance analysis of the Internet. Appendix In Section 3, we showed that the coupled system of C-TCP and E-RED (3.1) would exhibit a Hopf type bifurcation, as the stability conditions are violated. This results in the emergence of limit cycles. We conduct a local Hopf bifurcation analysis to determine the stability of these emerging limit cycles. In our analysis, we employ Poincaré normal forms and center manifold theory to determine the type of the Hopf bifurcation and orbital stability of the limit cycles (Hassard et al., 1981). Appendix A. Analysis For conducting local Hopf bifurcation analysis, we choose the exogenous non-dimensional parameter $$\kappa$$, to drive the system into a locally unstable region. System (3.1) may now be represented by the following set of nonlinear equations The nonlinear, delay differential equations for the underlying system are of the form \begin{equation*} \frac{d}{dt}\left[ \begin{array}{c} \mathrm{w}(t) \\ p(t) \end{array} \right] = \left[ \begin{array}{c} f(x,x_d,y_d) \\ g(x,y) \end{array} \right]\!, \end{equation*} where $$x, \,y, \, x_d \,\,\text{and}\,\, y_d$$ correspond to $$\mathrm{w}(t),\, p(t),\,\mathrm{w}(t-\tau) \,\,\text{and}\,\, p(t-\tau),$$ respectively. We perform a Taylor series expansion of (3.1), about its equilibrium (3.2), to get \begin{align} \dot{u}_{1}(t)&= \kappa{}\Big(\xi_x u_1(t) + \xi_{y_d} u_2(t-\tau) + \xi_{xx} u_1^2(t)\nonumber \\ & \quad+ \xi_{xx_d} u_1(t)u_1(t-\tau) + \xi_{xy_d} u_1(t)u_2(t-\tau) \nonumber \\ & \quad+ \xi_{x_dy_d} u_1(t-\tau)u_2(t-\tau) + \xi_{xxx} u_1^3(t) \label{eq:e21}\\ & \quad+ \xi_{xxx_d} u_1^2(t)u_1(t-\tau) + \xi_{xxy_d}u_1^2(t)u_2(t-\tau) \nonumber\\ & \quad+ \xi_{xx_dy_d} u_1(t)u_1(t-\tau)u_2(t-\tau) + \cdots\Big) \nonumber \\ \dot{u}_{2}(t) &= \kappa{} \bigg(\frac{k_lp_s}{\mathrm{w}_s} u_1(t) + \frac{k_l}{\mathrm{w}_s} u_1(t)u_2(t)\bigg), \nonumber \end{align} (A.1) where \begin{equation*} \xi_{x^i x_d^j y_d^k} =\frac{1}{i!\,j!\,k!}\frac{\partial^{i+j+k}}{\partial x^i \partial x_d^j \partial y_d^k}f\,\bigg|_{(\mathrm{w}_s,\,p_s)}. \end{equation*} The linear, quadratic and cubic terms in (A.1) are listed in Table A2. Recall that the system satisfies the Hopf condition at $$\kappa=\kappa_c$$, and as $$\kappa$$ goes beyond this critical value the nonlinear system will be locally unstable. Let $$\kappa{} = \kappa_{c} + \mu$$ and so the Hopf bifurcation takes place at $$\mu=0$$. Consider the following autonomous delay differential system: Table A2 Coefficients of the linear, quadratic and cubic terms in the Taylor series expansion of system (3.1) at equilibrium $$\xi_x$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\tau}\,(k-2)(1-p_s)$$ $$\xi_{y_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\frac{\mathrm{w}_s}{p_s}$$ $$\xi_{xx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\,(k-2)(1-p_s)$$ $$\xi_{xy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\,\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\xi_{x_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{p_s\tau}$$ $$\xi_{xxx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{6\mathrm{w}^2_s\tau}\,(k-1)(k-2)(k-3)(1-p_s)$$ $$\xi_{xxx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}^2_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xxy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)$$ $$\xi_{xx_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\xi_x$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\tau}\,(k-2)(1-p_s)$$ $$\xi_{y_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\frac{\mathrm{w}_s}{p_s}$$ $$\xi_{xx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\,(k-2)(1-p_s)$$ $$\xi_{xy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\,\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\xi_{x_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{p_s\tau}$$ $$\xi_{xxx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{6\mathrm{w}^2_s\tau}\,(k-1)(k-2)(k-3)(1-p_s)$$ $$\xi_{xxx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}^2_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xxy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)$$ $$\xi_{xx_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\bigg(k-2+\frac{1}{p_s}\bigg)$$ Table A2 Coefficients of the linear, quadratic and cubic terms in the Taylor series expansion of system (3.1) at equilibrium $$\xi_x$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\tau}\,(k-2)(1-p_s)$$ $$\xi_{y_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\frac{\mathrm{w}_s}{p_s}$$ $$\xi_{xx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\,(k-2)(1-p_s)$$ $$\xi_{xy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\,\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\xi_{x_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{p_s\tau}$$ $$\xi_{xxx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{6\mathrm{w}^2_s\tau}\,(k-1)(k-2)(k-3)(1-p_s)$$ $$\xi_{xxx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}^2_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xxy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)$$ $$\xi_{xx_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\xi_x$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\tau}\,(k-2)(1-p_s)$$ $$\xi_{y_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\frac{\mathrm{w}_s}{p_s}$$ $$\xi_{xx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\,(k-2)(1-p_s)$$ $$\xi_{xy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\,\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\xi_{x_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{p_s\tau}$$ $$\xi_{xxx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{6\mathrm{w}^2_s\tau}\,(k-1)(k-2)(k-3)(1-p_s)$$ $$\xi_{xxx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}^2_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xxy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)$$ $$\xi_{xx_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\label{eq:e22} \frac{d}{dt}\boldsymbol{u}(t) = \mathcal{L_{\mu}}\boldsymbol{u}_{t}+\mathcal{F}(\boldsymbol{u}_{t},\mu)$$ (A.2)$$t>0$$, $$\mu$$$$\in$$$$\mathbb{R}$$, where for $$\tau>0$$ $\boldsymbol{u}_{t}(\theta) = \boldsymbol{u}(t+\theta) \qquad \boldsymbol{u}:\big[-\tau,0\big] \rightarrow \mathbb{R}^{2}, \qquad \theta \in [-\tau,0].$ $$\mathcal{L_{\mu}}$$ is a one-parameter family of continuous (bounded) linear operators. The operator $$\mathcal{F}(\boldsymbol{u}_t,\mu)$$ contains the nonlinear terms. Further assume that $$\mathcal{F}$$ is analytic and that $$\mathcal{F}$$ and $$\mathcal{L_{\mu}}$$ depend analytically on the bifurcation parameter $$\mu$$ for small $$|\mu|$$. Note that (A.1) is of the form (A.2), where $$\boldsymbol{u} = \big[u_1\,\, u_2\big]^{T}$$. The objective now is to cast (A.2) into the form: $$\label{eq:e23} \frac{d}{dt}\boldsymbol{u}_{t} = \mathcal{A}(\mu)\boldsymbol{u}_{t}+\mathcal{R}\boldsymbol{u}_{t}$$ (A.3) which has $$\boldsymbol{u}_{t}$$ rather than both $$\boldsymbol{u}$$ and $$\boldsymbol{u}_{t}$$. First, we transform the linear problem $$d\boldsymbol{u}(t)/dt = \mathcal{L_{\mu}}\boldsymbol{u}_{t}$$. By the Riesz representation theorem, there exists a $$2$$$$\times$$$$2$$ matrix-valued function $$\eta(\cdot,\mu):[-\tau,0] \rightarrow \mathbb{R}^{2\times 2}$$, such that the components of $$\eta$$ has bounded variation and for all $$\boldsymbol{\phi}$$$$\in$$$$C[-\tau,0]$$ \begin{align} \mathcal{L_{\mu}}\boldsymbol{\phi} &= \int^{0}_{-\tau}d\eta(\theta,\mu)\boldsymbol{\phi}(\theta). \notag \end{align} In particular, \begin{align} \mathcal{L_{\mu}}\boldsymbol{u}_t &= \int^{0}_{-\tau}d\eta(\theta,\mu) \boldsymbol{u} (t+\theta). \label{eq:e24} \end{align} (A.4) Observe that \begin{align} \mathrm{\textit{d}}\eta(\theta,\mu) &= \kappa{} \left[ \begin{array}{@{}cc@{}} \xi_x\delta(\theta) & \xi_{y_d}\delta(\theta+\tau) \\ \frac{k_lp_s}{\mathrm{w}_s}\delta(\theta) & 0 \end{array} \right]\mathrm{\textit{d}}\theta, \end{align} (A.5) where $$\delta(\theta)$$ is the Dirac delta function, satisfies (A.4). For $$\boldsymbol{\phi}$$$$\in$$$$C^1\big[-\tau,0\big]$$, define \begin{align} \mathcal{A}(\mu)\boldsymbol{\phi}(\theta) &= \begin{cases}\hspace{-1mm} \begin{array}{c} \frac{d\boldsymbol{\phi}(\theta)}{d\theta},\\ \int^{0}_{-\tau}\mathrm{\textit{d}}\eta(s,\mu)\boldsymbol{\phi}(s) \equiv \mathcal{L_{\mu}}\boldsymbol{\phi}, \end{array}\!\!\! \begin{array}{l} \theta \in [-\tau,0) \\ \theta = 0 \end{array}\!\!\! \end{cases} \label{eq:e25} \end{align} (A.6) and \begin{align} \mathcal{R}\boldsymbol{\phi}(\theta) &= \begin{cases} \begin{array}{l} 0,\\ \mathcal{F}(\boldsymbol{\phi},\mu),\end{array} \begin{array}{l} \theta \in [-\tau,0) \\ \theta = 0. \end{array} \end{cases} \label{eq:e26} \end{align} (A.7) Then, as $$d\boldsymbol{u}_{t}/d\theta \equiv d\boldsymbol{u}_{t}/d{t}$$, (A.2) becomes (A.3). We now proceed to determine the coefficients required for the Hopf bifurcation analysis. Note that $$\kappa = \kappa_c +\mu$$ is the bifurcation parameter under consideration. As $$\kappa = \kappa_c$$ is the point of bifurcation, we set $$\mu=0$$ in order to compute the required terms. Let $$\boldsymbol{q}(\theta)$$ be the eigenfunction for $$\mathcal{A}(0)$$ corresponding to $$\lambda(0)$$, namely $\mathcal{A}(0)\boldsymbol{q}(\theta)= j\omega_{0}\boldsymbol{q}(\theta).$ To find $$\boldsymbol{q}(\theta)$$, let $$\boldsymbol{q}(\theta)=\boldsymbol{q}_0e^{j\omega_{0}\theta}$$, where $$\boldsymbol{q}_0=[q_{01}\,\,\,q_{02}]^{T}=[1\,\,\,\phi_1{}]^{T}$$. Substituting this in the above equation and using the expression for $$\mathcal{A}$$ in (A.6), we get \label{eq:e27} \begin{aligned} \phi_1 = \! \frac{\kappa{}}{j\omega_{0}}\frac{k_lp_s}{\mathrm{w}_s}, \,\, \omega_{0} = \kappa_c\sqrt{\frac{- \xi^2_{x}+\sqrt{\xi^4_{x}+4 \xi^2_{y_d} \big(k_lp_s/\mathrm{w}_s\big)^2}}{2}}. \end{aligned} (A.8) Define the adjoint operator $$\mathcal{A^{*}}(0)$$ as $\mathcal{A^{*}}(0)\boldsymbol{\alpha}(s) = \begin{cases} \begin{array}{c} -\frac{\mathrm{\textit{d}}\boldsymbol{\alpha}(s)}{d{s}},\\ \int^{0}_{-\tau}\mathrm{\textit{d}}\eta^{T}(t,0)\boldsymbol{\alpha}(-t),\end{array} & \begin{array}{l} s \in (0,\tau]\\ s = 0. \end{array}\end{cases}$ Note that the domains of $$\mathcal{A}$$ and $$\mathcal{A^{*}}$$ are $$C^1[-\tau,0]$$ and $$C^1[0,\tau]$$, respectively. By definition $\mathcal{A}(0)\boldsymbol{q}(\theta)= \lambda(0)\boldsymbol{q}(\theta),$ $$\overline{\lambda}(0)$$ is an eigenvalue for $$\mathcal{A^{*}}$$ and $\mathcal{A^{*}}(0)\boldsymbol{q}^*=-j\omega_{0}\boldsymbol{q}^*,$ for some non-zero vector $$\boldsymbol{q}^*$$. Let $$\boldsymbol{q}^*(s)=\boldsymbol{B}e^{j\omega_{0}s}$$ be an eigenvector of $$\mathcal{A^{*}}$$ corresponding to eigenvalue $$-j\omega_{0}$$, where $$\boldsymbol{B}=[B_{1}\,\,\,B_{2}]^{T}$$. We obtain $$\boldsymbol{B}=B[\phi_2\,\,\,1]^{T}, \,\,\,\, \text{where} \,\,\, \phi_2 = \frac{-\kappa{}}{\kappa{}\xi_x+j\omega_{0}}\frac{k_lp_s}{\mathrm{w}_s}.$$ For $$\boldsymbol{\phi}\in{}C[-\tau,0]$$ and $$\boldsymbol{\psi}\in{}C[0,\tau]$$, define an inner product $$\label{eq:e28} \langle\boldsymbol{\psi},\boldsymbol{\phi}\rangle=\overline{\boldsymbol{\psi}}^{T}(0)\boldsymbol{\phi}(0)-\underset{\theta=-\tau}{\int^{0}}\underset{\zeta=0}{\int^{\theta}}\overline{\boldsymbol{\psi}}^{T}(\zeta-\theta)d\eta(\theta)\boldsymbol{\phi}(\zeta)d\zeta.$$ (A.9) Then, $$\langle\boldsymbol{\psi},\mathcal{A}\boldsymbol{\phi}\rangle = \langle\mathcal{A^{*}}\boldsymbol{\psi},\boldsymbol{\phi}\rangle$$ for $$\boldsymbol{\phi} \in \mathrm{Dom}(\mathcal{A})$$ and $$\boldsymbol{\psi}\in{}\mathrm{Dom}(\mathcal{A}^*)$$. We will now find $$\boldsymbol{B}$$ such that the eigenvectors $$\boldsymbol{q}$$ and $$\boldsymbol{q}^*$$ satisfy the conditions $$\langle \boldsymbol{q}^*,\boldsymbol{q} \rangle=1$$ and $$\langle \boldsymbol{q}^*,\overline{\boldsymbol{q}} \rangle=0$$. Using the definition of the inner product defined in (A.9), we obtain $$1=\overline{B}\big(\phi_1+\overline{\phi}_{2}+\xi_{y_d}\kappa{}\tau{}\phi_{1}\overline{\phi}_{2} e^{-j\omega_{0}\tau}\big),$$ which gives $$B=\Big[\overline{\phi}_1+\phi_{2}+\xi_{y_d}\kappa{}\tau{}\overline{\phi}_1\phi_{2} e^{j\omega_{0}\tau}\Big]^{-1}.$$ Using the expression for $$B$$, we can verify that $$\langle \boldsymbol{q}^*,\overline{\boldsymbol{q}} \rangle=0$$. For $$\boldsymbol{u}_t$$, a solution of (A.3) at $$\mu=0$$, define \begin{aligned} z(t) &= \langle{}\boldsymbol{q}^*,\boldsymbol{u}_{t}\rangle, \\ \mathrm{\boldsymbol{w}}(t,\theta) &= \boldsymbol{u}_t(\theta)-2\text{Re}\big( z(t)\boldsymbol{q}(\theta)\big). \end{aligned} Then, on the manifold, $$C_{0}$$, $$\mathrm{\boldsymbol{w}}(t,\theta)=\mathrm{\boldsymbol{w}}\big(z(t),\overline{z}(t),\theta\big)$$, where $$\label{eq:e31} \mathrm{\boldsymbol{w}}\big(z,\overline{z},\theta\big)=\mathrm{\boldsymbol{w}}_{20}(\theta)\frac{z^{2}}{2}+\mathrm{\boldsymbol{w}}_{11}(\theta)z\overline{z}+\mathrm{\boldsymbol{w}}_{02}(\theta)\frac{\overline{z}^{2}}{2}+\cdots.$$ (A.10) In essence, $$z$$ and $$\overline{z}$$ are local coordinates for $$C_{0}$$ in $$C$$ in the directions of $$\boldsymbol{q}^*$$ and $$\overline{\boldsymbol{q}}^{*}$$, respectively. The existence of the center manifold $$C_{0}$$ enables the reduction of (A.3) to an ordinary differential equation for a single complex variable on $$C_{0}$$. At $$\mu=0$$, this is \begin{align} z'(t) &= \langle \boldsymbol{q}^*,\mathcal{A}\boldsymbol{u}_t+\mathcal{R}\boldsymbol{u}_t\rangle \notag \\ &= j\omega_{0}z(t)+\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}\Bigl(\!\mathrm{\boldsymbol{w}}(z,\overline{z},\theta)+2\text{Re}\big( z(t)\boldsymbol{q}(\theta)\big)\!\Bigr) \notag \\ &= j\omega_{0}z(t)+\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}(z,\overline{z}), \label{eq:e32} \end{align} (A.11) where $$a \cdot{} b$$ means $$\sum_{i=1}^{n}a_{i}b_{i}$$ and $$z'(t)$$ can be abbreviated as $$\label{eq:e33} z'(t) = j\omega_{0}z(t)+g\big(z,\overline{z}\big).$$ (A.12) Now the objective is to expand $$g$$ in powers of $$z$$ and $$\overline{z}$$ and determine the coefficients $$\mathrm{\boldsymbol{w}}_{ij}(\theta)$$ in (A.10). The differential (A.11) for $$z$$ would be explicit when we determine $$\mathrm{\boldsymbol{w}}_{ij}$$. Expanding $$g(z,\overline{z})$$ in powers of $$z$$ and $$\overline{z}$$ we have \begin{aligned} g(z,\overline{z}) &= \overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}\big(z,\overline{z}\big)\\ &= g_{20}\frac{z^{2}}{2}+g_{11}z\overline{z}+g_{02}\frac{\overline{z}^{2}}{2}+g_{21}\frac{z^{2}\overline{z}}{2}\cdots, \end{aligned} Following Hassard et al. (1981), we write $\mathrm{\boldsymbol{w}}' = \boldsymbol{u}'_{t}-z'\boldsymbol{q}-\overline{z}'\overline{\boldsymbol{q}},$ and using (A.3) and (A.12), we obtain $\mathrm{\boldsymbol{w}}' = \begin{cases} \hspace{-2mm} \begin{array}{l}\mathcal{A}\mathrm{\boldsymbol{w}}-2\text{Re}\big(\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}\boldsymbol{q}(\theta)\big),\\ \mathcal{A}\mathrm{\boldsymbol{w}}-2\text{Re}\big(\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}q(0)\big)+\mathcal{F}_{0},\end{array} \!\! \begin{array}{l}\theta \in [-\tau,0)\\\theta = 0\end{array}\!\!\!\!\! \end{cases}$ which is rewritten as $$\label{eq:e34} \mathrm{\boldsymbol{w}}' = \mathcal{A}\mathrm{\boldsymbol{w}}+\boldsymbol{h}\big(z,\overline{z},\theta\big)$$ (A.13) using (A.10), where $$\label{eq:e35} \boldsymbol{h}(z,\overline{z},\theta)=\boldsymbol{h}_{20}(\theta)\frac{z^{2}}{2}+\boldsymbol{h}_{11}(\theta)z\overline{z}+\boldsymbol{h}_{02}(\theta)\frac{\overline{z}^{2}}{2}+\cdots.$$ (A.14) Note that, on $$C_{0}$$, near the origin $\mathrm{\boldsymbol{w}}' = \mathrm{\boldsymbol{w}}_{z}z'+\mathrm{\boldsymbol{w}}_{\overline{z}}\overline{z}'.$ We use (A.10) and (A.12) to replace $$\mathrm{\boldsymbol{w}}_{z}$$ and $$z'$$, and equating this with (A.13) we obtain (Hassard et al., 1981) \label{eq:e36} \begin{aligned} \big(2j\omega_{0}-\mathcal{A}\big)\mathrm{\boldsymbol{w}}_{20}(\theta) &= \boldsymbol{h}_{20}(\theta)\\ -\mathcal{A}\mathrm{\boldsymbol{w}}_{11}(\theta) &= \boldsymbol{h}_{11}(\theta)\\ -\big(2j\omega_{0}+\mathcal{A}\big)\mathrm{\boldsymbol{w}}_{02}(\theta) &= \boldsymbol{h}_{02}(\theta). \end{aligned} (A.15) We have \begin{aligned} \boldsymbol{u}_t(\theta)&=\mathrm{\boldsymbol{w}}\big(z,\overline{z},\theta\big)+z\boldsymbol{q}(\theta)+\overline{z}\,\overline{\boldsymbol{q}}(\theta)\\ &=\mathrm{\boldsymbol{w}}_{20}(\theta)\frac{z^{2}}{2}+\mathrm{\boldsymbol{w}}_{11}(\theta)z\overline{z}+\mathrm{\boldsymbol{w}}_{02}(\theta)\frac{\overline{z}^{2}}{2} \\ &+z\boldsymbol{q}_0e^{j\omega_{0}\theta}+\overline{z}\,\overline{\boldsymbol{q}}_{0}e^{-j\omega_{0}\theta}+\cdots, \end{aligned} from which we obtain $$\boldsymbol{u}_t(0)$$ and $$\boldsymbol{u}_t(-\tau)$$. As we require only the coefficients of $$z^{2}$$, $$z\overline{z}$$, $$\overline{z}^{2}$$ and $$z^{2}\overline{z}$$, we only retain these terms in the expansions that follow: \begin{align*} u^2_{1,t}(0) &= z^{2}+\overline{z}^{2}+2z\overline{z} +\big({w}_{201}(0)+2{w}_{111}(0)\big)z^{2}\overline{z}+\cdots, \\ u_{1,t}(0)u_{2,t}(0) &= \phi_{1}z^{2}+\overline{\phi}_{1}\overline{z}^{2}+\left(\phi_{1}+\overline{\phi}_{1}\right)z\overline{z} \\ &\quad+\bigg( \frac{{w}_{201}(0)}{2}\overline{\phi}_{1}+{w}_{111}(0)\phi_{1} +{w}_{112}(0)+\frac{{w}_{202}(0)}{2}\bigg)z^{2}\overline{z}+\cdots, \\ u_{1,t}(0)u_{1,t}(-\tau) &=e^{-j\omega_{0}\tau}z^{2}+e^{j\omega_{0}\tau}\overline{z}^{2}+\left(e^{j\omega_{0}\tau}+e^{-j\omega_{0}\tau}\right)z\overline{z} \\ &\quad+\bigg( \frac{{w}_{201}(0)}{2}e^{j\omega_{0}\tau}+{w}_{111}(0)e^{-j\omega_{0}\tau} +{w}_{111}(-\tau)+\frac{{w}_{201}(-\tau)}{2} \bigg)z^{2}\overline{z}+\cdots,\\ u_{1,t}(0)u_{2,t}(-\tau) &=\phi_1e^{-j\omega_{0}\tau}z^{2}+\overline{\phi}_1e^{j\omega_{0}\tau}\overline{z}^{2}+\left(\phi_1e^{-j\omega_{0}\tau}+\overline{\phi}_1e^{j\omega_{0}\tau}\right)z\overline{z} \\ &\quad+\bigg( \frac{{w}_{201}(0)}{2}\overline{\phi}_1e^{j\omega_{0}\tau}+{w}_{111}(0)\phi_1e^{-j\omega_{0}\tau}+{w}_{112}(-\tau)+\frac{{w}_{202}(-\tau)}{2} \bigg)z^{2}\overline{z}+\cdots, \\ u_{1,t}(-\tau)u_{2,t}(-\tau) &=\phi_1e^{-2j\omega_{0}\tau}z^{2}+\overline{\phi}_1e^{2j\omega_{0}\tau}\overline{z}^{2}+\left(\phi_1+\overline{\phi}_1\right)z\overline{z}+\bigg( \frac{{w}_{201}(-\tau)}{2}\overline{\phi}_1e^{j\omega_{0}\tau} \\ &\quad+{w}_{111}(-\tau)\phi_1e^{-j\omega_{0}\tau}+{w}_{112}(-\tau)e^{-j\omega_{0}\tau}+\frac{{w}_{202}(-\tau)}{2}e^{j\omega_{0}\tau} \bigg)z^{2}\overline{z}+\cdots, \\ u^3_{1,t}(0) &= 3z^{2}\overline{z}+\cdots, \\ u^2_{1,t}(0)u_{1,t}(-\tau) &= \big(e^{j\omega_{0}\tau}+2e^{-j\omega_{0}\tau}\big)z^{2}\overline{z}+\cdots, \\ u^2_{1,t}(0)u_{2,t}(-\tau)&=\big(2\phi_1e^{-j\omega_{0}\tau}+\overline{\phi}_{1}e^{j\omega_{0}\tau}\big)z^{2}\overline{z}+\cdots, \\ u_{1,t}(0)u_{1,t}(-\tau)u_{2,t}(-\tau)&=\big(\phi_1+ \overline{\phi}_{1}+\phi_1{}e^{-2j\omega_{0}\tau}\big)z^{2}\overline{z}+\cdots, \end{align*} where $$[{w}_{ij1}\,\,{w}_{ij2}]^{T}=\mathrm{\boldsymbol{w}}_{ij}$$. Recall that, $g(z,\overline{z})=\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}(z,\overline{z})\equiv \overline{B}\big(\overline{\phi}_2\mathcal{F}_{01}+\mathcal{F}_{02}\big),$ where $$\mathcal{F}_{0}=[\mathcal{F}_{01}\,\,\mathcal{F}_{02}]^{T}\,$$, and $g(z,\overline{z}) = g_{20}\frac{z^{2}}{2}+g_{11}z\overline{z}+g_{02}\frac{\overline{z}^{2}}{2}+g_{21}\frac{z^{2}\overline{z}}{2}\cdots.$ By comparing the coefficients of $$z^{2}$$, $$z\overline{z}$$, $$\overline{z}^{2}$$ and $$z^{2}\overline{z}$$, we get $$g_{20}$$, $$g_{11}$$, $$g_{02}$$ and $$g_{21}$$ as \begin{align} g_{20} &= 2\kappa{}\overline{B}\bigg(\overline{\phi}_2\Big(\xi_{xx}+\xi_{xx_d}e^{-j\omega_{0}\tau}+\xi_{xy_d}\phi_1e^{-j\omega_{0}\tau} + \xi_{x_dy_d}\phi_1e^{-2j\omega_{0}\tau}\Big)+ \frac{k_l}{\mathrm{w}_s}\phi_1\bigg), \label{eq:e37} \end{align} (A.16) \begin{align} g_{11} &= \kappa{}\overline{B}\biggl(\overline{\phi}_2\Big(2\xi_{xx}+\xi_{xx_d}\big(e^{j\omega_{0}\tau}+e^{-j\omega_{0}\tau}\big) \nonumber\\ &\quad +\xi_{xy_d}\big(\,\overline{\phi}_1e^{j\omega_{0}\tau}+\phi_1e^{-j\omega_{0}\tau}\big)+\xi_{x_dy_d}\big(\phi_1+\overline{\phi}_1\big)\Big) + \frac{k_l}{\mathrm{w}_s}\big(\phi_1+\overline{\phi}_1\big)\biggr), \label{eq:e38} \end{align} (A.17) \begin{align} g_{02} &= 2\kappa{}\overline{B}\bigg(\overline{\phi}_2\Big( \xi_{xx}+\xi_{xx_d}e^{j\omega_{0}\tau}+\xi_{xy_d}\overline{\phi}_1e^{j\omega_{0}\tau} +\xi_{x_dy_d}\overline{\phi}_1e^{2j\omega_{0}\tau}\Big)+\frac{k_l}{\mathrm{w}_s}\overline{\phi}_1 \bigg), \label{eq:e39} \end{align} (A.18) \begin{align} g_{21} &= \kappa{}\overline{B}\Bigg(\overline{\phi}_2\bigg[2\xi_{xx}\Big({w}_{201}(0)+2{w}_{111}(0)\Big) + \xi_{xx_d}\Big(2{w}_{111}(-\tau) + {w}_{201}(0)e^{j\omega_{0}\tau} + {w}_{201}(-\tau) \nonumber \\ &\quad + 2e^{-j\omega_{0}\tau}{w}_{111}(0) \Big) + \xi_{xy_d}\Big(2{w}_{111}(0)\phi_{1}e^{-j\omega_{0}\tau} + {w}_{202}(-\tau) + 2{w}_{112}(-\tau) + {w}_{201}(0)\overline{\phi}_{1}e^{j\omega_{0}\tau} \Big) \nonumber \\ &\quad + \xi_{x_dy_d}\Big( {w}_{201}(-\tau)\overline{\phi}_{1}e^{j\omega_{0}\tau}+2{w}_{111}(-\tau)\phi_{1}e^{-j\omega_{0}\tau} + 2{w}_{112}(-\tau)e^{-j\omega_{0}\tau}+{w}_{202}(-\tau)e^{j\omega_{0}\tau} \Big) \nonumber \\ &\quad +2\bigg(3\xi_{xxx} + \xi_{xxx_d}\Big(e^{j\omega_{0}\tau}+2e^{-j\omega_{0}\tau}\Big)+\xi_{xxy_d}\Big(2\phi_1e^{-j\omega_{0}\tau} + \overline{\phi}_{1}e^{j\omega_{0}\tau}\Big) \nonumber \\ &\quad + \xi_{xx_dy_d}\Big(2\text{Re}\big(\phi_1\big)+\phi_1{}e^{-2j\omega_{0}\tau}\Big)\bigg)\bigg] + \frac{k_l}{\mathrm{w}_s} \Big( {w}_{201}(0)\overline{\phi}_{1} + 2{w}_{111}(0)\phi_{1} + 2{w}_{112}(0)+{w}_{202}(0)\Big)\Biggr). \label{eq:e40} \end{align} (A.19) For the expression of $$g_{21}$$, we have to evaluate $$\mathrm{\boldsymbol{w}}_{11}(0)$$, $$\mathrm{\boldsymbol{w}}_{11}(-\tau)$$, $$\mathrm{\boldsymbol{w}}_{20}(0)$$ and $$\mathrm{\boldsymbol{w}}_{20}(-\tau)$$. Now for $$\theta$$$$\in$$$$[-\tau,0)$$ \begin{aligned} \boldsymbol{h}(z,\overline{z},\theta) &= -2\text{Re}\big(\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}\boldsymbol{q}(\theta)\big)\\ &= -2\text{Re}\big(g(z,\overline{z})\boldsymbol{q}(\theta)\big)\\ &= -\left(g_{20}\frac{z^{2}}{2}+g_{11}z\overline{z}+g_{02}\frac{\overline{z}^{2}}{2}+\cdots \right)\boldsymbol{q}(\theta) -\left(\overline{g}_{20}\frac{\overline{z}^{2}}{2}+\overline{g}_{11}z\overline{z}+\overline{g}_{02}\frac{z^{2}}{2}+\cdots\right)\overline{\boldsymbol{q}}(\theta), \end{aligned} which when compared with (A.14), yields \begin{aligned} \boldsymbol{h}_{20}(\theta)&=-g_{20}\boldsymbol{q}(\theta)-\overline{g}_{02}\overline{\boldsymbol{q}}(\theta)\\ \boldsymbol{h}_{11}(\theta)&=-g_{11}\boldsymbol{q}(\theta)-\overline{g}_{11}\overline{\boldsymbol{q}}(\theta). \end{aligned} From (A.6) and (A.15), we get \begin{align} \mathrm{\boldsymbol{w}}'_{20}(\theta) &= 2j\omega_{0}\mathrm{\boldsymbol{w}}_{20}(\theta) + g_{20}\boldsymbol{q}(\theta)+\overline{g}_{02}\overline{\boldsymbol{q}}(\theta) \label{eq:e41} \end{align} (A.20) \begin{align} \mathrm{\boldsymbol{w}}'_{11}(\theta) &= g_{11}\boldsymbol{q}(\theta)+\overline{g}_{11}\overline{\boldsymbol{q}}(\theta) \label{eq:e42}. \end{align} (A.21) Solving the differential (A.20) and (A.21), we obtain \begin{align} \mathrm{\boldsymbol{w}}_{20}(\theta) &= -\frac{g_{20}}{j\omega_{0}}\boldsymbol{q}_0e^{j\omega_{0}\theta}-\frac{\overline{g}_{02}}{3j\omega_{0}}\overline{\boldsymbol{q}}_{0}e^{-j\omega_{0}\theta} +\boldsymbol{E}e^{2j\omega_{0}\theta}\!\!\!\! \label{eq:e43} \end{align} (A.22) \begin{align} \mathrm{\boldsymbol{w}}_{11}(\theta) &= \frac{g_{11}}{j\omega_{0}}\boldsymbol{q}_0e^{j\omega_{0}\theta}-\frac{\overline{g}_{11}}{j\omega_{0}}\overline{\boldsymbol{q}}_{0}e^{-j\omega_{0}\theta} +\boldsymbol{F}, \label{eq:e44} \end{align} (A.23) for some $$\boldsymbol{E}=[e_{1}\,\,e_{2}]^{T}$$ and $$\boldsymbol{F}=[f_{1} \,\, f_{2}]^{T}$$, which we will now determine. For $$\boldsymbol{h}(z,\overline{z},0) = -2\text{Re}\big(\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}q(0)\big)+\mathcal{F}_{0}$$, \begin{aligned} \boldsymbol{h}_{20}(0)=& -g_{20}\boldsymbol{q}(0)-\overline{g}_{02}\overline{\boldsymbol{q}}(0) +\left[\begin{array}{c}A_1\\A_2\end{array}\right]\\ \boldsymbol{h}_{11}(0)=& -g_{11}\boldsymbol{q}(0)-\overline{g}_{11}\overline{\boldsymbol{q}}(0) +\left[\begin{array}{c}B_1\\B_2\end{array}\right]\!, \end{aligned} where \begin{aligned} A_1 &= 2\kappa{}\Big(\xi_{xx}+\xi_{xx_d}e^{-j\omega_{0}\tau}+\xi_{xy_d}\phi_1e^{-j\omega_{0}\tau} + \xi_{x_dy_d}\phi_1e^{-2j\omega_{0}\tau}\Big), \\ B_1 &= 2\kappa{}\Big(\xi_{xx}+\text{Re}\big(\xi_{xx_d}e^{j\omega_{0}\tau}+\xi_{xy_d}\phi_1e^{-j\omega_{0}\tau} +\xi_{x_dy_d}\phi_1\big)\Big), \\ A_2 &= 2\kappa{}\phi_1 \frac{k_l}{\mathrm{w}_s}\,\,\, \text{and}\,\,\, B_2 = \kappa{} \big(\phi_1+\overline{\phi}_1\big) \frac{k_l}{\mathrm{w}_s}. \end{aligned} From (A.6) and (A.15), we get \begin{align} &g_{20}\boldsymbol{q}(0)+\overline{g}_{02}\overline{\boldsymbol{q}}(0) = \left[\begin{array}{c}A_1\\A_2\end{array}\right] + \left[ \begin{array}{c} \big(\kappa{} \xi_x-2j\omega_{0}\big){w}_{201}(0)+\kappa{}\xi_{y_d}{w}_{202}(-\tau) \\ \kappa{} \big(k_lp_s/\mathrm{w}_s\big) {w}_{201}(0)-2j\omega_{0}{w}_{202}(0) \end{array} \right]\!,\label{eq:e45} \end{align} (A.24) \begin{align} &g_{11}\boldsymbol{q}(0)+\overline{g}_{11}\overline{\boldsymbol{q}}(0) = \left[\begin{array}{c}B_1\\B_2\end{array}\right] +\left[ \begin{array}{c} \kappa \big(\xi_x {w}_{111}(0)+\xi_{y_d}{w}_{112}(-\tau)\big) \\ \kappa \big(k_lp_s/\mathrm{w}_s\big) {w}_{111}(0) \end{array} \right]\!. \label{eq:e46} \end{align} (A.25) Substituting the expression for $$\mathrm{\boldsymbol{w}}_{20}(0)$$, $$\mathrm{\boldsymbol{w}}_{20}(-\tau)$$, $$\mathrm{\boldsymbol{w}}_{11}(0)$$ and $$\mathrm{\boldsymbol{w}}_{11}(-\tau)$$ from (A.22) and (A.23) into (A.24) and (A.25), respectively, and solving for $$e_{1}$$, $$e_{2}$$, $$f_{1}$$ and $$f_{2}$$, we obtain \begin{align} e_1 &= \frac{1}{D}\Big(2j\omega_{0}A_1 + \kappa{}\xi_{y_d}A_2e^{-2j\omega_{0}\tau}\Big), \nonumber \\ e_2 &= \frac{1}{D}\Big(\kappa{} \frac{k_lp_s}{\mathrm{w}_s} A_1 - \big(\kappa{} \xi_x-2j\omega_{0}\big)A_2 \Big), \nonumber \\ f_1 &= \frac{-1}{p_s} \big(\phi_1 + \overline{\phi}_1 \big), \nonumber \\ f_2 &= \frac{-1}{\kappa{}\xi_{y_d}}\Big(\kappa{}\xi_x f_1 +B_1 \Big), \nonumber \end{align} where \begin{aligned} D &= -\bigg(2j\omega_{0}\big(\kappa{} \xi_x-2j\omega_{0}\big)+\kappa{}^2\xi_{y_d}\frac{k_lp_s}{\mathrm{w}_s} e^{-2j\omega_{0}\tau}\bigg). \end{aligned} Using the values of $$\boldsymbol{E}$$ and $$\boldsymbol{F}$$ in (A.22) and (A.23), followed by substituting $$\theta=0\,\, \text{and}\,-\tau$$, we can obtain the expressions for $$\mathrm{\boldsymbol{w}}_{20}(0)$$, $$\mathrm{\boldsymbol{w}}_{20}(-\tau)$$, $$\mathrm{\boldsymbol{w}}_{11}(0)$$ and $$\mathrm{\boldsymbol{w}}_{11}(-\tau)$$. With these expressions, we can now compute the term $$g_{21}$$. We now have the expressions for $$g_{20}$$, $$g_{11}$$, $$g_{02}$$ and $$g_{21}$$. All the necessary quantities that are essential for the Hopf bifurcation analysis have been computed. We now use the following expressions to analyse the type of the Hopf bifurcation (Hassard et al., 1981) $$\label{eq:e51} c_{1}(0) = \frac{j}{2\omega_{0}}\left(g_{20}g_{11}-2|g_{11}|^{2}-\frac{1}{3}|g_{02}|^{2}\right)+\frac{g_{21}}{2},$$ (A.26) \label{eq:e50} \begin{aligned} \mu_{2} &= \frac{ -\text{Re}\big(c_{1}(0)\big) }{ \alpha'(0) }, & \mathcal{B}_{2} &= 2\text{Re} \big( c_{1}(0) \big), \end{aligned} (A.27) where $$c_{1}(0)$$ is the lyapunov coefficient, $$\mathcal{B}_{2}$$ is the Floquet exponent, and $$\alpha'(0) = \text{Re}\big(d\lambda/d\kappa\big)\big|_{\kappa = \kappa_c}$$. We now state the conditions that will enable us to characterize the type of the Hopf bifurcation and determine the orbital stability of the bifurcating periodic solutions. The sign of $$\mu_2$$ determines the type of Hopf bifurcation. The Hopf bifurcation is super-critical if $$\mu_2 > 0$$ and sub-critical if $$\mu_2 < 0$$. The sign of (the Floquet exponent) $$\mathcal{B}_{2}$$ determines the stability of the bifurcating periodic solutions. The periodic solutions are asymptotically orbitally stable when $$\mathcal{B}_{2} < 0$$ and unstable if $$\mathcal{B}_{2} >0$$. Appendix B. Computations We conduct some numerical computations to demonstrate the Hopf bifurcation analysis. For a choice of system parameter values, using local stability analysis, one can determine the critical value of the bifurcation parameter at which the coupled system of C-TCP and E-RED (3.1) undergoes a local Hopf bifurcation. As an example, consider the following parameter values: $$\beta = 0.78$$, $$k=0.01$$, $$\gamma = 1$$, $$C = 5$$ units and $$\tau = 1$$ unit. For these parameter values, with $$\alpha = 1$$ the Hopf bifurcation occurs at $$k_l = 2.0029$$ and with $$\alpha = 4$$, we get $$k_l = 1.9248$$. Numerous system parameters, including the feedback delay, can act as the bifurcation parameter. However, we consider an exogenous parameter to be the bifurcation parameter. Our approach is to determine the Hopf condition using system parameters and then use the exogenous parameter to push the system just beyond the Hopf boundary. We do this as we are interested to know the impact of model parameters on the system behaviour in a local neighbourhood of the Hopf condition. At the Hopf boundary, by construction, the value of the exogenous non-dimensional parameter $$\kappa$$ is 1. An increase in $$\kappa$$, just beyond the Hopf condition, would render the system into a locally unstable regime. Using the Hopf bifurcation analysis, we compute the values of $$\mu_2$$ and $$\mathcal{B}_2,$$ which are presented in Fig. B.1. It can be seen that $$\mu_2 > 0$$ and $$\mathcal{B}_2 < 0$$, which implies that the Hopf bifurcation is supercritical, and the limit cycles are orbitally stable. Fig. B.1. View largeDownload slide Plots of $$\mu_2$$ and $$\beta_2$$ as the bifurcation parameter $$\kappa$$ varies. For the choice of $$\beta = 0.78$$, $$k=0.01$$, $$\gamma = 1$$, $$C = 5$$ units and $$\tau = 1$$ unit, system (3.1) undergoes a local Hopf bifurcation at $$\kappa = 1$$, for $$\alpha = 1, k_l = 2.0029$$ and $$\alpha = 4, k_l = 1.9248$$. Recall that $$\mu_2 > 0$$ implies that the Hopf bifurcation is supercritical, and $$\mathcal{B}_2 < 0$$ implies that the limit cycles are asymptotically orbitally stable. Fig. B.1. View largeDownload slide Plots of $$\mu_2$$ and $$\beta_2$$ as the bifurcation parameter $$\kappa$$ varies. For the choice of $$\beta = 0.78$$, $$k=0.01$$, $$\gamma = 1$$, $$C = 5$$ units and $$\tau = 1$$ unit, system (3.1) undergoes a local Hopf bifurcation at $$\kappa = 1$$, for $$\alpha = 1, k_l = 2.0029$$ and $$\alpha = 4, k_l = 1.9248$$. Recall that $$\mu_2 > 0$$ implies that the Hopf bifurcation is supercritical, and $$\mathcal{B}_2 < 0$$ implies that the limit cycles are asymptotically orbitally stable. The associated bifurcation diagrams and phase portraits are presented in Fig. B.2. The bifurcation diagrams, Fig. B.2(a and b), show the emergence of limit cycles in the average sending rate and the drop probability, at $$\kappa = 1$$. The phase portraits in Fig. B.2(c and d), for $$\kappa = 1.05 > \kappa_c$$, highlight the existence of stable limit cycles. The computations were performed using the software MATLAB. Using the analytical framework, one can perform similar computations for any desired set of parameter values. Fig. B.2. View largeDownload slide Bifurcation diagrams (a and b) and phase portraits (c and d), for the coupled system of C-TCP and E-RED (3.1). It can be observed from the bifurcation diagrams that the system transits from a stable equilibrium to a stable limit cycle. The phase portraits, plotted for $$\kappa = 1.05 > \kappa_c$$, establish the presence of stable limit cycles around the equilibrium ($$\mathrm{w}_s = 5$$ and $$p_s = 0.0495$$). Fig. B.2. View largeDownload slide Bifurcation diagrams (a and b) and phase portraits (c and d), for the coupled system of C-TCP and E-RED (3.1). It can be observed from the bifurcation diagrams that the system transits from a stable equilibrium to a stable limit cycle. The phase portraits, plotted for $$\kappa = 1.05 > \kappa_c$$, establish the presence of stable limit cycles around the equilibrium ($$\mathrm{w}_s = 5$$ and $$p_s = 0.0495$$). Acknowledgements The authors thank Nizar Malangadan and Santosh Chavan for their comments and feedback on this article. References Athuraliya, S., Low, S. H., Li, V. H. & Yin, Q. ( 2001 ) REM: active queue management. IEEE Netw., 15 , 48 – 53 . Google Scholar CrossRef Search ADS Briscoe, B., Brunstrom, A., Petlund, A., Hayes, D., Ros, D., Tsang, J., Gjessing, S., Fairhurst, G., Griwodz, C. & Welzl, M. ( 2016 ) Reducing internet latency: a survey of techniques and their merits. IEEE Commun. Surveys Tuts., 18 , 2149 – 2196 . Google Scholar CrossRef Search ADS Cao, J., Cleveland, W. S., Gao, Y., Jeffay, K. Smith, F. D. & Weigle, M. ( 2004 ) Stochastic models for generating synthetic HTTP source traffic. Proceedings of the IEEE INFOCOM. Hong Kong, China , pp. 1546 – 1557 . Floyd, S. & Jacobson, V. ( 1993 ) Random early detection gateways for congestion avoidance. IEEE/ACM Trans. Netw., 1 , 397 – 413 . Google Scholar CrossRef Search ADS Gettys, J. & Nichols, K. ( 2012 ) Bufferbloat: dark buffers in the internet. Commun. ACM, 55 , 57 – 65 . Google Scholar CrossRef Search ADS Gong, Y., Rossi, D., Testa, C., Valenti, S. & Täht, M. D. ( 2014 ) Fighting the bufferbloat: on the coexistence of AQM and low priority congestion control. Comput. Netw., 65 , 255 – 267 . Google Scholar CrossRef Search ADS Ha, S., Rhee, I. & Xu, L. ( 2008 ) CUBIC: a new TCP-friendly high-speed TCP variant. ACM SIGOPS Operating Systems Review, 42 , 64 – 74 . Google Scholar CrossRef Search ADS Han, H., Hollot, C. V., Towsley, D. & Chait, Y. ( 2005 ) Synchronization of TCP flows in networks with small Droptail buffers. Proceedings of the IEEE Conference on Decision and Control and European Control Conference. Seville, Spain , pp. 6762 – 6767 . Hassard, B. D., Kazarinoff, N. D. & Wan, Y.-H. ( 1981 ) Theory and Applications of Hopf Bifurcation. Cambridge, MA : Cambridge University Press . Hong, G., Martin, J. & Westall, J. M. ( 2015 ) On fairness and application performance of active queue management in broadband cable networks. Comput. Netw., 91 , 390 – 406 . Google Scholar CrossRef Search ADS Issariyakul, T. & Hossain, E. ( 2011 ) Introduction to Network Simulator NS2. New York : Springer . Kuznetsov, Y. A. ( 2004 ) Elements of Applied Bifurcation Theory. New York : Springer . Google Scholar CrossRef Search ADS Liu, S., Basar, T. & Srikant, R. ( 2005 ) Exponential-RED: a stabilizing AQM scheme for low and high-speed TCP protocols. IEEE/ACM Trans. Netw., 13 , 1068 – 1081 . Google Scholar CrossRef Search ADS Low, S. H. & Srikant, R. ( 2004 ) A mathematical framework for designing a low-loss, low-delay Internet. Netw. Spat. Econ., 4 , 75 – 101 . Google Scholar CrossRef Search ADS Malangadan, N., Rahman, H. & Raina, G. ( 2013 ) Non-linear oscillations in TCP networks with Drop-Tail buffers. Proceedings of the Chinese Control and Decision Conference, IEEE. Guiyang, China , pp. 188 – 194 . Nichols, K. & Jacobson, V. ( 2012 ) Controlling queue delay. Commun. ACM, 55 , 42 – 50 . Google Scholar CrossRef Search ADS Patil, G., McClean, S. & Raina, G. ( 2011 ) Drop tail and RED queue management with small buffers: stability and Hopf bifurcation. ICTACT J. Commun. Technol., 2 , 339 – 344 . Google Scholar CrossRef Search ADS Prasad, S. & Raina, G. ( 2014 ) Local Hopf bifurcation analysis of Compound TCP with an Exponential-RED queue management policy. Proceedings of the Chinese Control and Decision Conference, IEEE. Changsha, China , pp. 2588 – 2594 . Raja, P. & Raina, G. ( 2012 ) Delay and loss-based transport protocols: buffer-sizing and stability. Proceedings of the International Conference on Communication Systems and Networks, IEEE. Bangalore, India , pp. 1 – 10 . Raina, G., Towsley, D. & Wischik, D. ( 2005 ) Part II: control theory for buffer sizing. ACM SIGCOMM Comput. Commun. Rev., 35 , 79 – 82 . Google Scholar CrossRef Search ADS Raina, G., Manjunath, S., Prasad, S. & Giridhar, K. ( 2016 ) Stability and performance analysis of compound TCP With REM and drop-tail queue management. IEEE/ACM Trans. Netw., 24 , 1961 – 1974 . Google Scholar CrossRef Search ADS Showail, A., Jamshaid, K. & Shihada, B. ( 2016 ) Buffer sizing in wireless networks: challenges, solutions, and opportunities. IEEE Commun. Mag., 54 , 130 – 137 . Google Scholar CrossRef Search ADS Sojoudi, S., Low, S. H. & Doyle, J. C. ( 2014 ) Buffering dynamics and stability of internet congestion controllers. IEEE/ACM Trans. Netw., 22 , 1808 – 1818 . Google Scholar CrossRef Search ADS Srikant, R. ( 2004 ) The Mathematics of Internet Congestion Control. Cambridge, MA : Birkhauser . Google Scholar CrossRef Search ADS Tan, K., Song, J., Zhang, Q. & Sridharan, M. ( 2006 ) A Compound TCP approach for high-speed and long distance networks. Proceedings of the IEEE INFOCOM, IEEE. Barcelona, Spain , pp. 1 – 12 . Wischik, D. & McKeown, N. ( 2005 ) Part I: buffer sizes for core routers. ACM SIGCOMM Comput. Commun. Rev., 35 , 75 – 78 . Google Scholar CrossRef Search ADS Yang, P., Shao, J., Luo, W., Xu, L., Deogun, J. & Lu, Y. ( 2014 ) TCP congestion avoidance algorithm identification. IEEE/ACM Trans. Netw., 22 , 1311 – 1324 . Google Scholar CrossRef Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Mathematical Control and Information Oxford University Press

# Stability and performance analysis of Compound TCP with the Exponential-RED and the Drop-Tail queue policies

, Volume Advance Article – Nov 25, 2017
23 pages

/lp/ou_press/stability-and-performance-analysis-of-compound-tcp-with-the-mGgdcuaJhE
Publisher
Oxford University Press
ISSN
0265-0754
eISSN
1471-6887
D.O.I.
10.1093/imamci/dnx049
Publisher site
See Article on Publisher Site

### Abstract

Abstract The analysis of transport protocols, along with queue management policies, forms an important aspect of performance evaluation for the Internet. In this article, we study Compound TCP (C-TCP), the default TCP in the Windows operating system, along with the Exponential-RED (E-RED) queue policy and the widely used Drop-Tail queue policy. We consider queuing delay, link utilization and the stability of the queue size as the performance metrics. We first analyse the stability properties of a nonlinear model for C-TCP coupled with the E-RED queue policy. We observe that this system, in its current form, may be difficult to stabilize as the feedback delay gets large. Further, using an exogenous and non-dimensional parameter, we show that the system loses local stability via a Hopf bifurcation, which gives rise to limit cycles. Employing Poincaré normal forms and the center manifold theory, we outline an analytical framework to characterize the type of the Hopf bifurcation and to determine the orbital stability of the emerging limit cycles. Numerical examples, stability charts and bifurcation diagrams complement our analysis. We also conduct packet-level simulations, with E-RED and Drop-Tail, in small and large buffer-sizing regimes. With large buffers, E-RED can achieve small queue sizes compared with Drop-Tail. However, it is difficult to maintain the stability of the E-RED policy as the feedback delay gets large. On the other hand, with small buffers, E-RED offers no clear advantage over the simple Drop-Tail queue policy. Our work can offer insights for the design of queue policies that can ensure low latency and stability. 1. Introduction It is well known that the design of transport protocols, embedded in end systems and queue management policies, implemented in Internet routers, can impact network performance. The algorithms in transport protocols update their sending window based on feedback, about congestion in the network, from routers. However, feedback from routers is rarely instantaneous. Thus, with time-delayed feedback, stability and bifurcation phenomenon become important considerations for performance evaluation. Fluid models, for congestion control algorithms, are quite useful to investigate the stability and bifurcation properties of the underlying system. Such analysis can help understand trade-offs among the system parameters and can also guide design principles. For surveys on models for Internet congestion control, see Low & Srikant (2004) and Srikant (2004). Compound TCP (C-TCP) is widely implemented in the Internet today, as it is the default transport protocol in the Windows operating system (Tan et al., 2006). Moreover, recent studies (Yang et al., 2014) have shown that C-TCP contributes nearly $$25\%$$ of the total traffic in the Internet. Hence, its evaluation can help us understand the performance of transport protocols that are used in real networks. For an analytical evaluation, it is useful to have a model. The model we use for Compound is the fluid model presented in Raja & Raina (2012). An important metric for the performance analysis of communication networks is queuing delays, which are on the rise (Nichols & Jacobson, 2012). Excessive queuing delays are normally attributed to large, and unmanaged, router buffers (Gettys & Nichols, 2012; Wischik & McKeown, 2005). Such delays are undesirable as they can negatively impact performance, especially for delay sensitive applications. One way to address the concern of rising queuing delays is to employ effective queue policies. We also note that, in addition to the design of TCP and queue management, buffering dynamics do also play a crucial role in ensuring stable operation of communication networks (Sojoudi et al., 2014). Queuing delays could be significant in the scenario where the router buffers are sized as per the bandwidth-delay rule and would be negligible with small buffers (Raina et al., 2005). Drop-Tail is a simple queue policy, which drops all incoming packets when the router buffer is full. Another popular queue policy is random early detection (RED) (Floyd & Jacobson, 1993). However, the RED policy is not used in practice as there are no universally accepted guidelines on how to tune the various RED parameters. For some work on C-TCP, and TCP Reno, with Drop-Tail and the RED policy, see Han et al. (2005), Malangadan et al. (2013) and Patil et al. (2011). These articles studied the stability properties of the underlying system, and exhibited the emergence of limit cycles as system parameters vary, using both theory and packet-level simulations. Random exponential marking (REM) queue policy, one of the early proposals of active queue management (AQM), aims to achieve a target queue size and a target link utilization at the same time (Athuraliya et al., 2001). A detailed study on the stability and performance of REM policy along with C-TCP was carried out by Raina et al. (2016). Thus, for the performance evaluation of TCP and queue management policies, stability is an important consideration for performance evaluation. This area continues to be an active area of research (see Gong et al., 2014; Hong et al., 2015; Briscoe et al., 2016; Showail et al., 2016 and references therein). In this article, we study C-TCP with the Exponential-RED (E-RED) queue policy (Liu et al., 2005) which aims for high link utilization, while maintaining small queue sizes. As the model is nonlinear, it is useful to study the local stability and the local bifurcation properties. In our previous work (Prasad & Raina, 2014), we showed that C-TCP coupled with the E-RED queue policy undergoes a Hopf bifurcation (Hassard et al., 1981; Kuznetsov, 2004), when the condition for local stability is violated. We also provided, by employing Poincaré normal forms and the center manifold theory, an analytical framework to characterize the type of the Hopf bifurcation and the asymptotic orbital stability of limit cycles. We employ an exogenous and non-dimensional parameter as the bifurcation parameter for the Hopf bifurcation analysis. In our current work, we also derive a sufficient condition for local stability of C-TCP with E-RED, which suggests the need for the joint design of transport protocols and queue management strategies. We also conduct packet-level simulations using the network simulator (NS2) (Issariyakul & Hossain, 2011) to highlight the existence of limit cycles in the queue size dynamics. Numerical examples, stability charts and bifurcation diagrams also complement our analysis. The rest of this article is organized as follows. In Section 2, we outline models of C-TCP and the E-RED queue policy. In Section 3, we conduct a local stability analysis for the coupled system. In Section 4, we perform some packet-level simulations and finally summarize our contributions in Section 5. The detailed Hopf bifurcation analysis, for the nonlinear model of C-TCP and E-RED queue policy, is presented in the Appendix. 2. Models In this section, we first outline nonlinear, fluid models of C-TCP and the E-RED queue policy. We then conduct local stability analysis for the coupled system of C-TCP and the E-RED policy. The stability results of C-TCP with the Drop-Tail policy obtained in Raina et al. (2016) are summarized in Section 3, as we are interested in evaluating the performance of E-RED and the Drop-Tail policy with C-TCP. 2.1. C-TCP in a small buffer regime We consider the many flows, nonlinear, fluid model for C-TCP proposed in Raja & Raina (2012). In this model, $$\mathrm{w}(t)$$ denotes the average window size, $$i\big(\mathrm{w}(t)\big)$$ represents the increase of the window per positive acknowledgement and $$d\big(\mathrm{w}(t)\big)$$ represents the decrease of the current window size per negative acknowledgement. In small buffer regime, the congestion avoidance algorithm of C-TCP updates its current window size, $$\mathrm{w}(t)$$ as follows (Raja & Raina, 2012): $\mathrm{w}(t+1) = \left\{ \begin{array}{ll} \mathrm{w}(t)\Big(1+\alpha \mathrm{w}(t)^{k-1}\Big) & \mbox{if no loss},\\ \mathrm{w}(t)\big(1-\beta\big)& \mbox{if loss}, \end{array} \right.$ where $$\alpha$$, $$k$$, $$\beta$$ are protocol parameters. The functional forms for $$i\big(\mathrm{w}(t)\big)$$ and $$d\big(\mathrm{w}(t)\big)$$ are $$\label{funcforms} i\big(\mathrm{w}(t)\big) = \alpha \mathrm{w}(t)^{k-1} \quad \mbox{and} \quad \textit{d}\big(\mathrm{w}(t)\big) = \beta \mathrm{w}(t).$$ (2.1) A generalized nonlinear model for the congestion avoidance phase of transport protocols is given by Raina et al. (2005) $$\frac{\textit{d}\mathrm{w}(t)}{dt} = \biggl(\! i\big(\mathrm{w}(t)\big)\big( 1 - p(t-\tau) \big) - d\big(\mathrm{w}(t)\big) p\big(t-\tau\big) \!\biggr) \frac{\mathrm{w}(t-\tau)}{\tau},\label{eq:compoundwindow}$$ (2.2) with equilibrium \begin{equation*} i(\mathrm{w}_s)\big( 1 - p(\cdot) \big) = d(\mathrm{w}_s) p(\cdot), \end{equation*} where $$p(\cdot)$$ is the packet drop probability, $$\tau$$ is the round-trip time and $$\mathrm{w}_s$$ is the equilibrium window size. If the average window size of all flows is $$\mathrm{w}(t)$$, then the average rate at which packets are sent is assumed to be approximated by $$x(t) = \mathrm{w}(t)/\tau$$. Different forms of $$p(\cdot)$$ would be appropriate for different queue management policies. Substituting the functional forms (2.1) into (2.2) we get the following nonlinear model for C-TCP $$\frac{d\mathrm{w}(t)}{dt} = \bigg(\! \alpha \mathrm{w}(t)^{k-1}\big( 1 - p(t-\tau) \big) - \beta \mathrm{w}(t) p\big(t-\tau\big) \!\bigg) \frac{\mathrm{w}(t-\tau)}{\tau},\label{eq:compoundwindownew}$$ (2.3) where $$p(\cdot)$$ is a function of the arrival rate at the queue. A fluid model for small buffer Drop-Tail is given by Raina et al. (2005) \begin{align} p(y_l) = (y_l/C_l)^B, \notag \end{align} where $$y_l$$ is the total arrival rate, $$C_l$$ is the link service rate and $$B$$ is the router buffer size. This fluid model combined with (2.3) yields a model for C-TCP with small buffers using a Drop-Tail policy. 2.2. E-RED queue policy The algorithm for E-RED has been proposed to stabilize TCP Reno by slightly modifying the original RED algorithm (Liu et al., 2005). The packet marking/dropping probability $$p(t)$$ in E-RED queue policy is set to be an exponential function of virtual queue length whose capacity is slightly smaller than the link capacity. At the packet level, $$p(t)$$ can be calculated as (Liu et al., 2005): \begin{align}\label{equation:ered algorithm} p(t) = \begin{cases} \begin{array}{l} 0 \\ p_{\rm min} e^{\frac{k_l}{\gamma{}C_l}\big(b(t)-th_{\rm min}\big)} \\ 1\end{array} \begin{array}{l} \text{if} \,\, 0 \leq b(t) < th_{\rm min}, \\ \text{if} \,\, th_{\rm min} < b(t) < th_{\rm max}, \\ \text{if} \,\, b(t) \geq th_{\rm max}, \end{array} \end{cases} \end{align} (2.4) where $$b(t)$$ is virtual queue length, $$p_{\rm min}$$ is the minimum probability, $$b(t) = th_{\rm min}$$, $$k_l>0$$ is the link gain, $$0<\gamma{}\leq1$$ is the target link utilization, $$C_l$$ is link capacity and $$th_{\rm min}, th_{\rm max}$$ are queue length thresholds. 2.3. Nonlinear model for E-RED We now outline the model that describe the dynamics of E-RED queue policy (Liu et al., 2005). The E-RED queue policy aims to achieve high link utilization by dynamically adapting the packet loss probability at the queues. E-RED adapts the packet loss probability $$p(t)$$ based on the aggregate rate at the link. The link dynamics of this queue policy is given by Liu et al. (2005) $$\label{equation:ered} \frac{dp(t)}{dt}=k_l\,\frac{p(t)}{\gamma{}C}\Big(x(t)-\gamma{}C\Big),$$ (2.5) where $$x(t)=\mathrm{w}(t)/\tau$$ is the average arrival rate at the queue, and $$C$$ is the per-flow service rate of the bottleneck link. 3. Local stability analysis We analyse the nonlinear model for C-TCP (2.3) combined with the E-RED queue policy (2.5). First, we derive the Hopf condition for C-TCP with the E-RED queue policy. There are numerous parameters that could act as the bifurcation parameter in the coupled system of C-TCP and the E-RED queue policy. For example, the Compound parameter $$\alpha$$, the link gain $$k_l$$ or the feedback delay $$\tau$$. Instead, we use an exogenous and non-dimensional parameter $$\kappa$$, as the bifurcation parameter. Introducing this parameter into the system, we have the following nonlinear equations \begin{align} \frac{d\mathrm{w}(t)}{dt} &= \kappa{}\bigg( \alpha \mathrm{w}(t)^{k-1}\big( 1 - p(t-\tau) \big) - \beta \mathrm{w}(t) p\big(t-\tau\big) \bigg) \frac{\mathrm{w}(t-\tau)}{\tau}, \nonumber \\ \frac{dp(t)}{dt} &= \kappa{}\left( k_l\frac{p(t)}{\gamma{}C}\big(x(t)-\gamma{}C \big) \right)\!. \label{equation:kappa} \end{align} (3.1) Note that the parameter $$\kappa$$ does not affect the system equilibrium. We assume that the model parameters are chosen such that the system is at the edge of local stability. At this condition, the value of $$\kappa$$ is $$1$$. Now, instead of varying the model parameters, we gradually vary $$\kappa$$ beyond the edge of local stability. Thus varying $$\kappa$$, beyond $$1$$, will render the nonlinear system in a locally unstable state. Let $$(\mathrm{w}_s, p_s)$$ be the equilibrium of system (3.1), then at equilibrium, we get the following conditions \label{eqpoint} \begin{aligned} \mathrm{w}_s = \gamma{}C\tau \,\,\,\, \text{and} \,\,\,\, p_s = \frac{\alpha \mathrm{w}_s^{k-1}}{\alpha \mathrm{w}_s^{k-1} + \beta \mathrm{w}_s}. \end{aligned} (3.2) Let us consider the perturbations, $$u_1(t) = \mathrm{w}(t)-\mathrm{w}_s$$ and $$u_2(t) = p(t)-p_s$$. Linearizing system (3.1), about the equilibrium (3.2), we obtain \begin{align} \frac{du_1(t)}{dt} &= \kappa{}\bigg((k-2)(1-p_s)u_1(t)-\frac{\mathrm{w}_s}{p_s}u_2(t-\tau) \bigg)\frac{i(\mathrm{w}_s)}{\tau}, \notag \\ \frac{du_2(t)}{dt} &= \kappa{}\frac{k_lp_s }{\mathrm{w}_s}u_1(t). \label{lin1} \end{align} (3.3) Looking for exponential solutions, we get the following characteristic equation $$\label{ceq} \lambda^2+a\kappa{}\lambda+\kappa{}^2be^{-\lambda{}\tau}=0,$$ (3.4) where \begin{aligned} a =\frac{-i(\mathrm{w}_s)}{\tau}(k-2)(1-p_s) \,\, \text{and}\,\, b = \frac{i(\mathrm{w}_s)}{\tau}k_l. \end{aligned} Let $$\lambda$$ be a root of the characteristic equation, (3.4), and we determine the critical value $$\kappa = \kappa_c$$ at which $$\lambda$$ is purely imaginary. Substituting $$\lambda = j\omega_0$$ in the characteristic equation (3.4), we get $$\label{lsa1} -\omega_0{}^2+j\kappa{}_c a\omega_0+\kappa{}^2_cb\,\big(\cos{}\omega_0{}\tau-j\sin{}\omega_0{}\tau\big)=0.$$ (3.5) From (3.5), we have \begin{equation*} \begin{aligned} \cos{}\omega_0{}\tau =\frac{\omega_0{}^2}{\kappa{}^2_cb} \,\,\,\, \text{and} \,\,\,\, \sin{}\omega_0{}\tau =\frac{\omega_0 a{}}{\kappa{}_cb}. \end{aligned} \end{equation*} Let $$\tilde{\omega}_0=\omega_0/\kappa_c$$. Squaring and adding the above equations, we obtain \begin{align} &\tilde{\omega}_0 = \sqrt{\frac{-a^2+\sqrt{a^4+4b^2}}{2}}. \label{tilome} \end{align} (3.6) We now check whether the roots of (3.4) cross imaginary axis with a non-zero velocity i.e., $$\text{Re}\big(d\lambda/d\kappa\big)|_{\kappa_c} \neq 0$$. Differentiating (3.4) with respect to $$\kappa$$, we have $\frac{d\lambda}{d\kappa}=\frac{2\lambda^2+a\kappa\lambda}{\lambda^2\kappa\tau+\lambda(a\kappa^2\tau+2\kappa)+a\kappa^2}.$ Evaluating the above expression at $$\kappa_c$$, we get \begin{align*} \text{Re}\bigg(\frac{d\lambda}{d\kappa}\bigg)_{\kappa = \kappa_c} = \frac{\tau\big(2\omega_0^4+\kappa^2_c{} a^2 \omega_0^2 \big)}{ \kappa_c\big((\kappa_c{}a -\omega_0^2\tau)^2+\omega_0^2(\kappa_c{}a \tau{}+2)^2\big)}>0, \end{align*} which satisfies the transversality condition of the Hopf spectrum (Hassard et al., 1981). Further, the Hopf bifurcation occurs at $$\kappa{}_c = \frac{1}{\tilde{\omega}_0\tau}\tan^{-1}\left(\frac{a}{\tilde{\omega}_0}\right) \label{hopfcondition}.$$ (3.7) Hence, the nonlinear system (3.1) is locally stable for $$\kappa < \kappa_c$$, and undergoes a local Hopf bifurcation at $$\kappa= \kappa_c$$. We first make an observation on the choice of the Compound parameter $$k$$, whose default value is $$0.75$$. Observe that for $$k=1$$, the function $$i\big(\mathrm{w}(t)\big)$$ becomes independent of the window size $$\mathrm{w}(t)$$, and the equilibrium probability is $$p_s=\alpha/(\alpha+\beta{}\gamma{}C\tau)$$. Thus, for a high bandwidth-delay environment $$(1-p_s) \approx 1$$. In Fig. 1, we plot the local stability chart of system (3.1). This plot shows the relationship between the Compound parameter $$\alpha$$ and the link gain $$k_l$$, if local stability is to be maintained. Figure 2 shows the relationship between the link gain $$k_l$$ and the feedback delay $$\tau$$, if local stability is to be ensured. Fig. 1. View largeDownload slide Local stability chart of system (3.1), capturing the relationship between the Compound parameter $$\alpha$$ and the link gain $$k_l$$. The other system parameter values are $$\beta=0.5$$, $$k=1$$ and $$\gamma=1$$, and we assume $$(1-p_s)\approx1$$. The region to the left of the curves is stable. Fig. 1. View largeDownload slide Local stability chart of system (3.1), capturing the relationship between the Compound parameter $$\alpha$$ and the link gain $$k_l$$. The other system parameter values are $$\beta=0.5$$, $$k=1$$ and $$\gamma=1$$, and we assume $$(1-p_s)\approx1$$. The region to the left of the curves is stable. Fig. 2. View largeDownload slide Local stability chart of system (3.1) showing variation of link gain $$k_l$$ with respect to feedback delay $$\tau$$. The other system parameter values are $$\alpha=0.125$$, $$\beta=0.5$$, $$k=1$$, $$\gamma=1$$ and we assume $$(1-p_s)\approx1$$. The region under the Hopf condition curve is stable. Fig. 2. View largeDownload slide Local stability chart of system (3.1) showing variation of link gain $$k_l$$ with respect to feedback delay $$\tau$$. The other system parameter values are $$\alpha=0.125$$, $$\beta=0.5$$, $$k=1$$, $$\gamma=1$$ and we assume $$(1-p_s)\approx1$$. The region under the Hopf condition curve is stable. 3.1. Sufficient condition From the characteristic equation, (3.4) and for $$\kappa = 1$$, we obtain the loop transfer function of linear system (3.3) as \begin{align} L(j\omega) = \frac{b e^{-j\omega{}\tau}}{j\omega(a + j\omega )}. \end{align} (3.8) Let $$\omega_{gc}$$ denote the gain crossover frequency of $$L(j\omega)$$, which can be obtained by solving $$| L(j\omega_{gc})| = 1$$. By employing Nyquist stability criterion, the coupled system of C-TCP and E-RED is locally asymptotically stable, if phase margin defined as $$\text{PM} = \pi + \angle L(j\omega_{gc})$$ is greater than zero. \begin{align}\label{equation:phase margin} \text{PM} &= \pi - \pi/2 -\omega_{gc}\tau - \tan^{-1}(\omega_{gc}/a) > 0, \nonumber \\ & \pi/2 -\omega_{gc}\tau - \tan^{-1}(\omega_{gc}/a) > 0, \nonumber \\ & \omega_{gc}\tau + \tan^{-1}(\omega_{gc}/a) < \pi/2, \end{align} (3.9) where $$\quad \omega_{gc} = \sqrt{(-a^2 + \sqrt{a^4 + 4b^2})/{2}}$$. Note that the gain crossover frequency, $$\omega_{gc}$$, increases with increase in $$b$$, and if the value of $$b$$ gets large, then the stability condition given by (3.9) cannot be satisfied. Thus, we consider the values of $$b$$ to be smaller to obtain a sufficient condition for local stability of (3.1). For values of $$a,\,b$$ satisfying $$2b/a^2 << 1$$, $$\omega_{gc}$$ can be approximated to $$b/a$$. Now by substituting $$\omega_{gc} = b/a$$ in condition (3.9), we get $\frac{b\tau }{a}+ \tan^{-1}(b/a^2) < \pi/2.$ As $$(b/a^2) < 1$$, and using the expressions for $$a$$, $$b$$, we obtain the sufficient condition for stability as $$\label{equation:sufficient} \frac{k_l \tau }{(2-k)(1-p_s)} < \pi/4.$$ (3.10) For the choice of Compound parameter $$k=0$$, we get $k_l \tau \left(1 + \frac{\alpha}{\beta \mathrm{w}^2_s} \right) < \pi/2,$ where $$\mathrm{w}_s$$, $$p_s$$ are equilibrium window size and equilibrium packet-loss probability respectively. Sufficient condition (3.10) highlights that the choice of system parameter values is crucial in ensuring the stability of overall system (3.1). The stability analysis of C-TCP in a small buffer regime with Drop-Tail queues was conducted in Raina et al. (2016). The stability results of (i) C-TCP with E-RED and (ii) C-TCP with Drop-Tail are summarized in Table 1. Table 1 Local stability results of C-TCP with E-RED and Drop-Tail queue policies Queue policy Sufficient condition Hopf condition E-RED $$\displaystyle \frac{k_l \tau }{(2-k)(1-p_s)} < \pi/4$$ $$\displaystyle \kappa{}_c = \frac{1}{\tilde{\omega}_0\tau}\tan^{-1}\left(\frac{a}{\tilde{\omega}_0}\right)\!, \quad \tilde{\omega}_0 = \sqrt{\frac{-a^2+\sqrt{a^4+4b^2}}{2}}$$ Drop-Tail $$\displaystyle \alpha B \mathrm{w}^{k-1}_s < \frac{\pi}{2}$$ $$\displaystyle \alpha \mathrm{w}^{k-1}_s \sqrt{B^2 - \big(k-2\big)^2\big(1-p_s\big)^2} = \cos^{-1}\Big( \big(k-2\big)\big(1-p_s\big)\big/B \Big)$$ Queue policy Sufficient condition Hopf condition E-RED $$\displaystyle \frac{k_l \tau }{(2-k)(1-p_s)} < \pi/4$$ $$\displaystyle \kappa{}_c = \frac{1}{\tilde{\omega}_0\tau}\tan^{-1}\left(\frac{a}{\tilde{\omega}_0}\right)\!, \quad \tilde{\omega}_0 = \sqrt{\frac{-a^2+\sqrt{a^4+4b^2}}{2}}$$ Drop-Tail $$\displaystyle \alpha B \mathrm{w}^{k-1}_s < \frac{\pi}{2}$$ $$\displaystyle \alpha \mathrm{w}^{k-1}_s \sqrt{B^2 - \big(k-2\big)^2\big(1-p_s\big)^2} = \cos^{-1}\Big( \big(k-2\big)\big(1-p_s\big)\big/B \Big)$$ Table 1 Local stability results of C-TCP with E-RED and Drop-Tail queue policies Queue policy Sufficient condition Hopf condition E-RED $$\displaystyle \frac{k_l \tau }{(2-k)(1-p_s)} < \pi/4$$ $$\displaystyle \kappa{}_c = \frac{1}{\tilde{\omega}_0\tau}\tan^{-1}\left(\frac{a}{\tilde{\omega}_0}\right)\!, \quad \tilde{\omega}_0 = \sqrt{\frac{-a^2+\sqrt{a^4+4b^2}}{2}}$$ Drop-Tail $$\displaystyle \alpha B \mathrm{w}^{k-1}_s < \frac{\pi}{2}$$ $$\displaystyle \alpha \mathrm{w}^{k-1}_s \sqrt{B^2 - \big(k-2\big)^2\big(1-p_s\big)^2} = \cos^{-1}\Big( \big(k-2\big)\big(1-p_s\big)\big/B \Big)$$ Queue policy Sufficient condition Hopf condition E-RED $$\displaystyle \frac{k_l \tau }{(2-k)(1-p_s)} < \pi/4$$ $$\displaystyle \kappa{}_c = \frac{1}{\tilde{\omega}_0\tau}\tan^{-1}\left(\frac{a}{\tilde{\omega}_0}\right)\!, \quad \tilde{\omega}_0 = \sqrt{\frac{-a^2+\sqrt{a^4+4b^2}}{2}}$$ Drop-Tail $$\displaystyle \alpha B \mathrm{w}^{k-1}_s < \frac{\pi}{2}$$ $$\displaystyle \alpha \mathrm{w}^{k-1}_s \sqrt{B^2 - \big(k-2\big)^2\big(1-p_s\big)^2} = \cos^{-1}\Big( \big(k-2\big)\big(1-p_s\big)\big/B \Big)$$ 4. Packet-level simulations We conduct some packet-level simulations using the NS2 (Issariyakul & Hossain, 2011). The topology we consider is a single bottleneck link with a capacity of $$100$$ Mbps. We use the default parameter values of Compound (Tan et al., 2006): $$\alpha = 0.125$$, $$\beta = 0.5$$ and $$k = 0.75$$ and implemented the E-RED policy, as described in Liu et al. (2005), in NS2. The values of E-RED parameters chosen are: $$p_{\rm min} = 0.000001$$, $$th_{\rm min} = 1/5 \times$$ buffer size, $$th_{\rm max} = 3/5 \times$$ buffer size. The value of link gain parameter is chosen to satisfy the condition $$k_l*\tau_{\rm max} = 1$$. With the above choice of parameter values, we can determine the maximum drop probability, $$p_{\rm max}$$ from the definition of $$p(t)$$ (2.4). We consider two different traffic mixes: (i) long-lived flows and (ii) mix of long-lived and short-lived flows. The simulation run time is 225 s, and all the long-lived flows start randomly within the first 10 s. We illustrate the traces of the various quantities of interest over the last 25 s as the transient effects would settle down by then, and the system will be in steady state. We first conduct packet-level simulations of C-TCP, with Drop-Tail and E-RED, in a large and a small buffer regime in Section 4.1. The bottleneck link capacity is $$100$$ Mbps, and the large buffer regime (bandwidth-delay rule) translates into $$2084$$ packets. This is because a value of $$250$$ ms is used in practice, irrespective of the actual delays of the TCP flows (Wischik & McKeown, 2005). For the small buffer regime, parameters like the link capacity, the feedback delay, and the number of users do not influence the choice of buffer size. In this regime, we consider buffer sizes to be in the range of $$15$$–$$100$$ packets, as guided by the stability analysis (Raina et al., 2016). The traffic consists of $$60$$ long-lived C-TCP flows, each with a $$2$$ Mbps link feeding into the bottleneck. In Section 4.2, we explore the ability of E-RED to achieve desired link utilization and maintain small queue sizes with C-TCP. We consider the large buffer regime as it is currently widely deployed. In this set of simulations, the traffic mix is a combination of FTP, UDP and HTTP flows. The traffic consists of $$55$$ FTP flows of C-TCP, each with a $$2$$ Mbps link, and $$5$$ UDP flows contributing a total of $$5$$ Mbps. For short-lived flows, we use the PackMime-HTTP package in NS2. We generate $$10$$ new HTTP connections every second, on an average. The HTTP traffic model is described, in detail, in Cao et al. (2004). For E-RED we choose target link utilization of $$98 \%$$. The average round-trip time of the FTP flows present is $$10$$ ms and $$200$$ ms. For completeness, we also conduct simulations with the Drop-Tail policy. For all our simulations, we consider a fixed packet size of $$1500$$ bytes. 4.1. C-TCP, E-RED and Drop-Tail, Small and Large buffers With small buffers, as router buffer sizes vary from $$15$$ to $$100$$ packets, with both Drop-Tail and E-RED we note the emergence of stable limit cycles in the queue size, see Fig. 3(a) (buffer size = $$15$$ packets) and Fig. 3(b) (buffer size = $$100$$ packets). The analytical results of Compound with Drop-Tail predicted the loss of local stability, via a Hopf, as buffer sizes got larger, see Table 1. In case of E-RED queue policy, with buffer size of $$15$$ packets, observe that the virtual queue size oscillates as the round-trip time gets large. The analysis, for C-TCP with E-RED, carried out in Section 3 predicted the emergence of limit cycles as feedback delay increases. For the buffer size of $$15$$ packets, the transport layer protocol acts to control the distribution of the queue size, and it is possible that E-RED effectively degenerates into Drop-Tail. Fig. 3. View largeDownload slide Packet-level simulations of C-TCP with (i) Drop-Tail and (ii) E-RED queue policies employed individually at the bottleneck link. For the E-RED policy, we have plotted the real queue length (solid line) and the virtual queue length (dotted line). (a) With 15 packet buffers, the queue sizes with Drop-Tail and with E-RED policy do not exhibit limit cycles. However, the loss rate of C-TCP and E-RED is larger than that of C-TCP and Drop-Tail. (b)With 100 packet buffers, with Drop-Tail, note the emergence of limit cycles in queue size with larger feedback delay, which can hurt the link utilization. With E-RED queue policy, we notice limit cycles in virtual queue size (dotted line) for large RTT, which results in the loss of link utilization. (c)With large buffers, Drop-Tail achieves full utilization, but at the cost of large queuing delays. With ERED, the oscillations in queue size (real and virtual) for RTT = 200 ms are very clear, these oscillations result in loss of link utilization. Fig. 3. View largeDownload slide Packet-level simulations of C-TCP with (i) Drop-Tail and (ii) E-RED queue policies employed individually at the bottleneck link. For the E-RED policy, we have plotted the real queue length (solid line) and the virtual queue length (dotted line). (a) With 15 packet buffers, the queue sizes with Drop-Tail and with E-RED policy do not exhibit limit cycles. However, the loss rate of C-TCP and E-RED is larger than that of C-TCP and Drop-Tail. (b)With 100 packet buffers, with Drop-Tail, note the emergence of limit cycles in queue size with larger feedback delay, which can hurt the link utilization. With E-RED queue policy, we notice limit cycles in virtual queue size (dotted line) for large RTT, which results in the loss of link utilization. (c)With large buffers, Drop-Tail achieves full utilization, but at the cost of large queuing delays. With ERED, the oscillations in queue size (real and virtual) for RTT = 200 ms are very clear, these oscillations result in loss of link utilization. With large buffers and large RTTs, E-RED exhibits oscillations in queue size which leads to loss in link utilization, see Fig. 3(c). Drop-Tail, with large buffers, is expected to have full link utilization but this comes at the cost of enormous queuing delays. 4.2. C-TCP, E-RED and large buffers: short-lived and long-lived flows Large buffers are widely deployed in the Internet routers, and a good queue management policy should be able to ensure low latency, stability of queue sizes and offer a decent link utilization. We consider a traffic mix which consists of combination of $$55$$ FTP flows of C-TCP, each with a $$2$$ Mbps link, $$5$$ UDP flows (contributing a total of $$5$$ Mbps) and an average of $$10$$ new HTTP flows per second. With smaller round-trip times for the FTP flows, we observe that E-RED achieves the desired link utilization. With larger round-trip times, the queue sizes are small but the target link utilization is not achieved. For the traffic mix of FTP flows of C-TCP, UDP and HTTP flows, E-RED policy fails to achieve target link utilization as the feedback delay gets large, see Fig. 4(b). Thus, while latency could be reduced, the link utilization would be adversely impacted. Note that the nonlinear model of E-RED, analysed in Section 3, did alert us to the possibility of large delays destabilizing the system. Drop-Tail, with the same traffic mix, with large buffers has high utilization but excessive queuing delays, see Fig. 4(a). Fig. 4. View largeDownload slide (a) Drop-Tail and (b) E-RED queue policies with large buffers, with a traffic mix of FTP, UDP and HTTP flows. For E-RED policy, we have plotted the real queue length (solid line) and the virtual queue length (dotted line). Observe that the Drop-Tail policy achieves full link utilization but at the cost of large queuing delays. In case of E-RED policy, loss in link utilization occurs as the feedback delay increases. Thus, with higher RTTs, E-RED queue policy fails to achieve target link utilization. Fig. 4. View largeDownload slide (a) Drop-Tail and (b) E-RED queue policies with large buffers, with a traffic mix of FTP, UDP and HTTP flows. For E-RED policy, we have plotted the real queue length (solid line) and the virtual queue length (dotted line). Observe that the Drop-Tail policy achieves full link utilization but at the cost of large queuing delays. In case of E-RED policy, loss in link utilization occurs as the feedback delay increases. Thus, with higher RTTs, E-RED queue policy fails to achieve target link utilization. Clearly, the combination of C-TCP, queue management schemes and buffer size impacts network performance. The size of router buffers plays a central role in the end-to-end performance, as variations in buffer sizes influence both latency and queue stability. From the analysis and the packet-level simulations conducted so far, the E-RED queue policy does not seem to offer any clear advantage over Drop-Tail in terms of performance and stability of the coupled system (3.1). One can perform a variety of simulations considering different topologies and a wider range of traffic mixes, but it is out of the scope of this article. 5. Contributions Queuing delays in the Internet are on the rise (Nichols & Jacobson, 2012), and large queuing delays would adversely impact the Quality of Service of delay-sensitive web applications. One way to address the issue of queuing delays is to deploy effective AQM policies in Internet routers (Gettys & Nichols, 2012; Nichols & Jacobson, 2012). There have been many proposals for AQMs in the literature; however, there is no consensus on the optimal choice of an AQM policy. For the evaluation of queue management policies, it is also important to consider TCP variants that are currently implemented in the Internet. In this work, we focused on C-TCP, the default TCP in the Windows operating system. Among various queue management policies, we considered E-RED as it aims to offer high link utilization and achieve small queue sizes. We conducted a performance evaluation of the E-RED queue policy along with C-TCP flows operating over two different buffer-sizing regimes; Large (bandwidth-delay product) buffers and small buffers. Our study was based on the stability analysis of fluid models and packet-level simulations using the NS2 (Issariyakul & Hossain, 2011). We considered queuing delay, link utilization and the stability of the queue size as the performance metrics. We first conducted a local stability analysis of C-TCP along with the E-RED queue policy, and derived a sufficient, and the necessary and sufficient conditions for local stability. From the stability conditions, we highlighted that there exists a relation between C-TCP parameters, E-RED parameters and the feedback delay to ensure system stability. Such conditions on the parameters of TCP and AQM highlight the need for the joint design of transport protocols and queue management strategies. Using an exogenous and non-dimensional parameter, we also showed that the underlying system would undergo a Hopf type bifurcation. The non-dimensional parameter is also used to investigate the nature of the local instability. Further, we outlined a detailed local Hopf bifurcation analysis, using Poincaré normal forms and the center manifold theory (Hassard et al., 1981), to characterize the type of the Hopf bifurcation and to determine the orbital stability of the emerging limit cycles. Packet-level simulations conducted for different buffer-sizing regimes highlighted the existence of limit cycles in the queue size dynamics. These limit cycles in the queue size are undesirable and could lead to a loss in the link utilization. In case of large buffers, Drop-Tail can achieve full link utilization but at the cost of large queuing delays. With the same network setting, E-RED can maintain small queues but can fail to achieve the desired link utilization as the feedback delay gets large. Moreover, with E-RED, the stability of the overall system is also difficult to maintain with large delays. With small buffers and Drop-Tail, the system has low latency and stability can be ensured; however, the buffer size has to be dimensioned carefully. In such a regime, the E-RED policy does not offer any clear advantage over the simple Drop-Tail policy. The insights developed in this article could help in the design of queue management policies to ensure stability and low latency. 5.1. Avenues for further research To address the issues related to latency in the Internet, it is essential to analyse currently implemented transport protocols, of end systems, and queue management policies of the Internet routers. Compound TCP (Tan et al., 2006), the default protocol in Windows operating system, and Cubic TCP (Ha et al., 2008), protocol in Linux operating system, are widely used in the present Internet. To the best of our knowledge, fluid models for Cubic TCP have not been proposed in the literature. As fluid models are often used to conduct control-theoretic analysis, it would be useful to devise models for Cubic TCP, which may then aid performance evaluation. For queue management, it would be important to compare the performance of E-RED and Drop-Tail with other queue policies as well, especially in different buffer-sizing regimes. Adding to analysis, and packet-level simulations, experiments over a hardware platform would certainly aid the performance analysis of the Internet. Appendix In Section 3, we showed that the coupled system of C-TCP and E-RED (3.1) would exhibit a Hopf type bifurcation, as the stability conditions are violated. This results in the emergence of limit cycles. We conduct a local Hopf bifurcation analysis to determine the stability of these emerging limit cycles. In our analysis, we employ Poincaré normal forms and center manifold theory to determine the type of the Hopf bifurcation and orbital stability of the limit cycles (Hassard et al., 1981). Appendix A. Analysis For conducting local Hopf bifurcation analysis, we choose the exogenous non-dimensional parameter $$\kappa$$, to drive the system into a locally unstable region. System (3.1) may now be represented by the following set of nonlinear equations The nonlinear, delay differential equations for the underlying system are of the form \begin{equation*} \frac{d}{dt}\left[ \begin{array}{c} \mathrm{w}(t) \\ p(t) \end{array} \right] = \left[ \begin{array}{c} f(x,x_d,y_d) \\ g(x,y) \end{array} \right]\!, \end{equation*} where $$x, \,y, \, x_d \,\,\text{and}\,\, y_d$$ correspond to $$\mathrm{w}(t),\, p(t),\,\mathrm{w}(t-\tau) \,\,\text{and}\,\, p(t-\tau),$$ respectively. We perform a Taylor series expansion of (3.1), about its equilibrium (3.2), to get \begin{align} \dot{u}_{1}(t)&= \kappa{}\Big(\xi_x u_1(t) + \xi_{y_d} u_2(t-\tau) + \xi_{xx} u_1^2(t)\nonumber \\ & \quad+ \xi_{xx_d} u_1(t)u_1(t-\tau) + \xi_{xy_d} u_1(t)u_2(t-\tau) \nonumber \\ & \quad+ \xi_{x_dy_d} u_1(t-\tau)u_2(t-\tau) + \xi_{xxx} u_1^3(t) \label{eq:e21}\\ & \quad+ \xi_{xxx_d} u_1^2(t)u_1(t-\tau) + \xi_{xxy_d}u_1^2(t)u_2(t-\tau) \nonumber\\ & \quad+ \xi_{xx_dy_d} u_1(t)u_1(t-\tau)u_2(t-\tau) + \cdots\Big) \nonumber \\ \dot{u}_{2}(t) &= \kappa{} \bigg(\frac{k_lp_s}{\mathrm{w}_s} u_1(t) + \frac{k_l}{\mathrm{w}_s} u_1(t)u_2(t)\bigg), \nonumber \end{align} (A.1) where \begin{equation*} \xi_{x^i x_d^j y_d^k} =\frac{1}{i!\,j!\,k!}\frac{\partial^{i+j+k}}{\partial x^i \partial x_d^j \partial y_d^k}f\,\bigg|_{(\mathrm{w}_s,\,p_s)}. \end{equation*} The linear, quadratic and cubic terms in (A.1) are listed in Table A2. Recall that the system satisfies the Hopf condition at $$\kappa=\kappa_c$$, and as $$\kappa$$ goes beyond this critical value the nonlinear system will be locally unstable. Let $$\kappa{} = \kappa_{c} + \mu$$ and so the Hopf bifurcation takes place at $$\mu=0$$. Consider the following autonomous delay differential system: Table A2 Coefficients of the linear, quadratic and cubic terms in the Taylor series expansion of system (3.1) at equilibrium $$\xi_x$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\tau}\,(k-2)(1-p_s)$$ $$\xi_{y_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\frac{\mathrm{w}_s}{p_s}$$ $$\xi_{xx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\,(k-2)(1-p_s)$$ $$\xi_{xy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\,\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\xi_{x_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{p_s\tau}$$ $$\xi_{xxx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{6\mathrm{w}^2_s\tau}\,(k-1)(k-2)(k-3)(1-p_s)$$ $$\xi_{xxx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}^2_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xxy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)$$ $$\xi_{xx_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\xi_x$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\tau}\,(k-2)(1-p_s)$$ $$\xi_{y_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\frac{\mathrm{w}_s}{p_s}$$ $$\xi_{xx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\,(k-2)(1-p_s)$$ $$\xi_{xy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\,\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\xi_{x_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{p_s\tau}$$ $$\xi_{xxx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{6\mathrm{w}^2_s\tau}\,(k-1)(k-2)(k-3)(1-p_s)$$ $$\xi_{xxx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}^2_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xxy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)$$ $$\xi_{xx_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\bigg(k-2+\frac{1}{p_s}\bigg)$$ Table A2 Coefficients of the linear, quadratic and cubic terms in the Taylor series expansion of system (3.1) at equilibrium $$\xi_x$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\tau}\,(k-2)(1-p_s)$$ $$\xi_{y_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\frac{\mathrm{w}_s}{p_s}$$ $$\xi_{xx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\,(k-2)(1-p_s)$$ $$\xi_{xy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\,\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\xi_{x_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{p_s\tau}$$ $$\xi_{xxx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{6\mathrm{w}^2_s\tau}\,(k-1)(k-2)(k-3)(1-p_s)$$ $$\xi_{xxx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}^2_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xxy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)$$ $$\xi_{xx_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\xi_x$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\tau}\,(k-2)(1-p_s)$$ $$\xi_{y_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\frac{\mathrm{w}_s}{p_s}$$ $$\xi_{xx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\,(k-2)(1-p_s)$$ $$\xi_{xy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\tau}\,\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\xi_{x_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{p_s\tau}$$ $$\xi_{xxx}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{6\mathrm{w}^2_s\tau}\,(k-1)(k-2)(k-3)(1-p_s)$$ $$\xi_{xxx_d}$$ $$\displaystyle \frac{i(\mathrm{w}_s)}{2\mathrm{w}^2_s\tau}\,(k-1)(k-2)(1-p_s)$$ $$\xi_{xxy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{2\mathrm{w}_s\tau}\,(k-1)(k-2)$$ $$\xi_{xx_dy_d}$$ $$\displaystyle \frac{-i(\mathrm{w}_s)}{\mathrm{w}_s\tau}\bigg(k-2+\frac{1}{p_s}\bigg)$$ $$\label{eq:e22} \frac{d}{dt}\boldsymbol{u}(t) = \mathcal{L_{\mu}}\boldsymbol{u}_{t}+\mathcal{F}(\boldsymbol{u}_{t},\mu)$$ (A.2)$$t>0$$, $$\mu$$$$\in$$$$\mathbb{R}$$, where for $$\tau>0$$ $\boldsymbol{u}_{t}(\theta) = \boldsymbol{u}(t+\theta) \qquad \boldsymbol{u}:\big[-\tau,0\big] \rightarrow \mathbb{R}^{2}, \qquad \theta \in [-\tau,0].$ $$\mathcal{L_{\mu}}$$ is a one-parameter family of continuous (bounded) linear operators. The operator $$\mathcal{F}(\boldsymbol{u}_t,\mu)$$ contains the nonlinear terms. Further assume that $$\mathcal{F}$$ is analytic and that $$\mathcal{F}$$ and $$\mathcal{L_{\mu}}$$ depend analytically on the bifurcation parameter $$\mu$$ for small $$|\mu|$$. Note that (A.1) is of the form (A.2), where $$\boldsymbol{u} = \big[u_1\,\, u_2\big]^{T}$$. The objective now is to cast (A.2) into the form: $$\label{eq:e23} \frac{d}{dt}\boldsymbol{u}_{t} = \mathcal{A}(\mu)\boldsymbol{u}_{t}+\mathcal{R}\boldsymbol{u}_{t}$$ (A.3) which has $$\boldsymbol{u}_{t}$$ rather than both $$\boldsymbol{u}$$ and $$\boldsymbol{u}_{t}$$. First, we transform the linear problem $$d\boldsymbol{u}(t)/dt = \mathcal{L_{\mu}}\boldsymbol{u}_{t}$$. By the Riesz representation theorem, there exists a $$2$$$$\times$$$$2$$ matrix-valued function $$\eta(\cdot,\mu):[-\tau,0] \rightarrow \mathbb{R}^{2\times 2}$$, such that the components of $$\eta$$ has bounded variation and for all $$\boldsymbol{\phi}$$$$\in$$$$C[-\tau,0]$$ \begin{align} \mathcal{L_{\mu}}\boldsymbol{\phi} &= \int^{0}_{-\tau}d\eta(\theta,\mu)\boldsymbol{\phi}(\theta). \notag \end{align} In particular, \begin{align} \mathcal{L_{\mu}}\boldsymbol{u}_t &= \int^{0}_{-\tau}d\eta(\theta,\mu) \boldsymbol{u} (t+\theta). \label{eq:e24} \end{align} (A.4) Observe that \begin{align} \mathrm{\textit{d}}\eta(\theta,\mu) &= \kappa{} \left[ \begin{array}{@{}cc@{}} \xi_x\delta(\theta) & \xi_{y_d}\delta(\theta+\tau) \\ \frac{k_lp_s}{\mathrm{w}_s}\delta(\theta) & 0 \end{array} \right]\mathrm{\textit{d}}\theta, \end{align} (A.5) where $$\delta(\theta)$$ is the Dirac delta function, satisfies (A.4). For $$\boldsymbol{\phi}$$$$\in$$$$C^1\big[-\tau,0\big]$$, define \begin{align} \mathcal{A}(\mu)\boldsymbol{\phi}(\theta) &= \begin{cases}\hspace{-1mm} \begin{array}{c} \frac{d\boldsymbol{\phi}(\theta)}{d\theta},\\ \int^{0}_{-\tau}\mathrm{\textit{d}}\eta(s,\mu)\boldsymbol{\phi}(s) \equiv \mathcal{L_{\mu}}\boldsymbol{\phi}, \end{array}\!\!\! \begin{array}{l} \theta \in [-\tau,0) \\ \theta = 0 \end{array}\!\!\! \end{cases} \label{eq:e25} \end{align} (A.6) and \begin{align} \mathcal{R}\boldsymbol{\phi}(\theta) &= \begin{cases} \begin{array}{l} 0,\\ \mathcal{F}(\boldsymbol{\phi},\mu),\end{array} \begin{array}{l} \theta \in [-\tau,0) \\ \theta = 0. \end{array} \end{cases} \label{eq:e26} \end{align} (A.7) Then, as $$d\boldsymbol{u}_{t}/d\theta \equiv d\boldsymbol{u}_{t}/d{t}$$, (A.2) becomes (A.3). We now proceed to determine the coefficients required for the Hopf bifurcation analysis. Note that $$\kappa = \kappa_c +\mu$$ is the bifurcation parameter under consideration. As $$\kappa = \kappa_c$$ is the point of bifurcation, we set $$\mu=0$$ in order to compute the required terms. Let $$\boldsymbol{q}(\theta)$$ be the eigenfunction for $$\mathcal{A}(0)$$ corresponding to $$\lambda(0)$$, namely $\mathcal{A}(0)\boldsymbol{q}(\theta)= j\omega_{0}\boldsymbol{q}(\theta).$ To find $$\boldsymbol{q}(\theta)$$, let $$\boldsymbol{q}(\theta)=\boldsymbol{q}_0e^{j\omega_{0}\theta}$$, where $$\boldsymbol{q}_0=[q_{01}\,\,\,q_{02}]^{T}=[1\,\,\,\phi_1{}]^{T}$$. Substituting this in the above equation and using the expression for $$\mathcal{A}$$ in (A.6), we get \label{eq:e27} \begin{aligned} \phi_1 = \! \frac{\kappa{}}{j\omega_{0}}\frac{k_lp_s}{\mathrm{w}_s}, \,\, \omega_{0} = \kappa_c\sqrt{\frac{- \xi^2_{x}+\sqrt{\xi^4_{x}+4 \xi^2_{y_d} \big(k_lp_s/\mathrm{w}_s\big)^2}}{2}}. \end{aligned} (A.8) Define the adjoint operator $$\mathcal{A^{*}}(0)$$ as $\mathcal{A^{*}}(0)\boldsymbol{\alpha}(s) = \begin{cases} \begin{array}{c} -\frac{\mathrm{\textit{d}}\boldsymbol{\alpha}(s)}{d{s}},\\ \int^{0}_{-\tau}\mathrm{\textit{d}}\eta^{T}(t,0)\boldsymbol{\alpha}(-t),\end{array} & \begin{array}{l} s \in (0,\tau]\\ s = 0. \end{array}\end{cases}$ Note that the domains of $$\mathcal{A}$$ and $$\mathcal{A^{*}}$$ are $$C^1[-\tau,0]$$ and $$C^1[0,\tau]$$, respectively. By definition $\mathcal{A}(0)\boldsymbol{q}(\theta)= \lambda(0)\boldsymbol{q}(\theta),$ $$\overline{\lambda}(0)$$ is an eigenvalue for $$\mathcal{A^{*}}$$ and $\mathcal{A^{*}}(0)\boldsymbol{q}^*=-j\omega_{0}\boldsymbol{q}^*,$ for some non-zero vector $$\boldsymbol{q}^*$$. Let $$\boldsymbol{q}^*(s)=\boldsymbol{B}e^{j\omega_{0}s}$$ be an eigenvector of $$\mathcal{A^{*}}$$ corresponding to eigenvalue $$-j\omega_{0}$$, where $$\boldsymbol{B}=[B_{1}\,\,\,B_{2}]^{T}$$. We obtain $$\boldsymbol{B}=B[\phi_2\,\,\,1]^{T}, \,\,\,\, \text{where} \,\,\, \phi_2 = \frac{-\kappa{}}{\kappa{}\xi_x+j\omega_{0}}\frac{k_lp_s}{\mathrm{w}_s}.$$ For $$\boldsymbol{\phi}\in{}C[-\tau,0]$$ and $$\boldsymbol{\psi}\in{}C[0,\tau]$$, define an inner product $$\label{eq:e28} \langle\boldsymbol{\psi},\boldsymbol{\phi}\rangle=\overline{\boldsymbol{\psi}}^{T}(0)\boldsymbol{\phi}(0)-\underset{\theta=-\tau}{\int^{0}}\underset{\zeta=0}{\int^{\theta}}\overline{\boldsymbol{\psi}}^{T}(\zeta-\theta)d\eta(\theta)\boldsymbol{\phi}(\zeta)d\zeta.$$ (A.9) Then, $$\langle\boldsymbol{\psi},\mathcal{A}\boldsymbol{\phi}\rangle = \langle\mathcal{A^{*}}\boldsymbol{\psi},\boldsymbol{\phi}\rangle$$ for $$\boldsymbol{\phi} \in \mathrm{Dom}(\mathcal{A})$$ and $$\boldsymbol{\psi}\in{}\mathrm{Dom}(\mathcal{A}^*)$$. We will now find $$\boldsymbol{B}$$ such that the eigenvectors $$\boldsymbol{q}$$ and $$\boldsymbol{q}^*$$ satisfy the conditions $$\langle \boldsymbol{q}^*,\boldsymbol{q} \rangle=1$$ and $$\langle \boldsymbol{q}^*,\overline{\boldsymbol{q}} \rangle=0$$. Using the definition of the inner product defined in (A.9), we obtain $$1=\overline{B}\big(\phi_1+\overline{\phi}_{2}+\xi_{y_d}\kappa{}\tau{}\phi_{1}\overline{\phi}_{2} e^{-j\omega_{0}\tau}\big),$$ which gives $$B=\Big[\overline{\phi}_1+\phi_{2}+\xi_{y_d}\kappa{}\tau{}\overline{\phi}_1\phi_{2} e^{j\omega_{0}\tau}\Big]^{-1}.$$ Using the expression for $$B$$, we can verify that $$\langle \boldsymbol{q}^*,\overline{\boldsymbol{q}} \rangle=0$$. For $$\boldsymbol{u}_t$$, a solution of (A.3) at $$\mu=0$$, define \begin{aligned} z(t) &= \langle{}\boldsymbol{q}^*,\boldsymbol{u}_{t}\rangle, \\ \mathrm{\boldsymbol{w}}(t,\theta) &= \boldsymbol{u}_t(\theta)-2\text{Re}\big( z(t)\boldsymbol{q}(\theta)\big). \end{aligned} Then, on the manifold, $$C_{0}$$, $$\mathrm{\boldsymbol{w}}(t,\theta)=\mathrm{\boldsymbol{w}}\big(z(t),\overline{z}(t),\theta\big)$$, where $$\label{eq:e31} \mathrm{\boldsymbol{w}}\big(z,\overline{z},\theta\big)=\mathrm{\boldsymbol{w}}_{20}(\theta)\frac{z^{2}}{2}+\mathrm{\boldsymbol{w}}_{11}(\theta)z\overline{z}+\mathrm{\boldsymbol{w}}_{02}(\theta)\frac{\overline{z}^{2}}{2}+\cdots.$$ (A.10) In essence, $$z$$ and $$\overline{z}$$ are local coordinates for $$C_{0}$$ in $$C$$ in the directions of $$\boldsymbol{q}^*$$ and $$\overline{\boldsymbol{q}}^{*}$$, respectively. The existence of the center manifold $$C_{0}$$ enables the reduction of (A.3) to an ordinary differential equation for a single complex variable on $$C_{0}$$. At $$\mu=0$$, this is \begin{align} z'(t) &= \langle \boldsymbol{q}^*,\mathcal{A}\boldsymbol{u}_t+\mathcal{R}\boldsymbol{u}_t\rangle \notag \\ &= j\omega_{0}z(t)+\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}\Bigl(\!\mathrm{\boldsymbol{w}}(z,\overline{z},\theta)+2\text{Re}\big( z(t)\boldsymbol{q}(\theta)\big)\!\Bigr) \notag \\ &= j\omega_{0}z(t)+\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}(z,\overline{z}), \label{eq:e32} \end{align} (A.11) where $$a \cdot{} b$$ means $$\sum_{i=1}^{n}a_{i}b_{i}$$ and $$z'(t)$$ can be abbreviated as $$\label{eq:e33} z'(t) = j\omega_{0}z(t)+g\big(z,\overline{z}\big).$$ (A.12) Now the objective is to expand $$g$$ in powers of $$z$$ and $$\overline{z}$$ and determine the coefficients $$\mathrm{\boldsymbol{w}}_{ij}(\theta)$$ in (A.10). The differential (A.11) for $$z$$ would be explicit when we determine $$\mathrm{\boldsymbol{w}}_{ij}$$. Expanding $$g(z,\overline{z})$$ in powers of $$z$$ and $$\overline{z}$$ we have \begin{aligned} g(z,\overline{z}) &= \overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}\big(z,\overline{z}\big)\\ &= g_{20}\frac{z^{2}}{2}+g_{11}z\overline{z}+g_{02}\frac{\overline{z}^{2}}{2}+g_{21}\frac{z^{2}\overline{z}}{2}\cdots, \end{aligned} Following Hassard et al. (1981), we write $\mathrm{\boldsymbol{w}}' = \boldsymbol{u}'_{t}-z'\boldsymbol{q}-\overline{z}'\overline{\boldsymbol{q}},$ and using (A.3) and (A.12), we obtain $\mathrm{\boldsymbol{w}}' = \begin{cases} \hspace{-2mm} \begin{array}{l}\mathcal{A}\mathrm{\boldsymbol{w}}-2\text{Re}\big(\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}\boldsymbol{q}(\theta)\big),\\ \mathcal{A}\mathrm{\boldsymbol{w}}-2\text{Re}\big(\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}q(0)\big)+\mathcal{F}_{0},\end{array} \!\! \begin{array}{l}\theta \in [-\tau,0)\\\theta = 0\end{array}\!\!\!\!\! \end{cases}$ which is rewritten as $$\label{eq:e34} \mathrm{\boldsymbol{w}}' = \mathcal{A}\mathrm{\boldsymbol{w}}+\boldsymbol{h}\big(z,\overline{z},\theta\big)$$ (A.13) using (A.10), where $$\label{eq:e35} \boldsymbol{h}(z,\overline{z},\theta)=\boldsymbol{h}_{20}(\theta)\frac{z^{2}}{2}+\boldsymbol{h}_{11}(\theta)z\overline{z}+\boldsymbol{h}_{02}(\theta)\frac{\overline{z}^{2}}{2}+\cdots.$$ (A.14) Note that, on $$C_{0}$$, near the origin $\mathrm{\boldsymbol{w}}' = \mathrm{\boldsymbol{w}}_{z}z'+\mathrm{\boldsymbol{w}}_{\overline{z}}\overline{z}'.$ We use (A.10) and (A.12) to replace $$\mathrm{\boldsymbol{w}}_{z}$$ and $$z'$$, and equating this with (A.13) we obtain (Hassard et al., 1981) \label{eq:e36} \begin{aligned} \big(2j\omega_{0}-\mathcal{A}\big)\mathrm{\boldsymbol{w}}_{20}(\theta) &= \boldsymbol{h}_{20}(\theta)\\ -\mathcal{A}\mathrm{\boldsymbol{w}}_{11}(\theta) &= \boldsymbol{h}_{11}(\theta)\\ -\big(2j\omega_{0}+\mathcal{A}\big)\mathrm{\boldsymbol{w}}_{02}(\theta) &= \boldsymbol{h}_{02}(\theta). \end{aligned} (A.15) We have \begin{aligned} \boldsymbol{u}_t(\theta)&=\mathrm{\boldsymbol{w}}\big(z,\overline{z},\theta\big)+z\boldsymbol{q}(\theta)+\overline{z}\,\overline{\boldsymbol{q}}(\theta)\\ &=\mathrm{\boldsymbol{w}}_{20}(\theta)\frac{z^{2}}{2}+\mathrm{\boldsymbol{w}}_{11}(\theta)z\overline{z}+\mathrm{\boldsymbol{w}}_{02}(\theta)\frac{\overline{z}^{2}}{2} \\ &+z\boldsymbol{q}_0e^{j\omega_{0}\theta}+\overline{z}\,\overline{\boldsymbol{q}}_{0}e^{-j\omega_{0}\theta}+\cdots, \end{aligned} from which we obtain $$\boldsymbol{u}_t(0)$$ and $$\boldsymbol{u}_t(-\tau)$$. As we require only the coefficients of $$z^{2}$$, $$z\overline{z}$$, $$\overline{z}^{2}$$ and $$z^{2}\overline{z}$$, we only retain these terms in the expansions that follow: \begin{align*} u^2_{1,t}(0) &= z^{2}+\overline{z}^{2}+2z\overline{z} +\big({w}_{201}(0)+2{w}_{111}(0)\big)z^{2}\overline{z}+\cdots, \\ u_{1,t}(0)u_{2,t}(0) &= \phi_{1}z^{2}+\overline{\phi}_{1}\overline{z}^{2}+\left(\phi_{1}+\overline{\phi}_{1}\right)z\overline{z} \\ &\quad+\bigg( \frac{{w}_{201}(0)}{2}\overline{\phi}_{1}+{w}_{111}(0)\phi_{1} +{w}_{112}(0)+\frac{{w}_{202}(0)}{2}\bigg)z^{2}\overline{z}+\cdots, \\ u_{1,t}(0)u_{1,t}(-\tau) &=e^{-j\omega_{0}\tau}z^{2}+e^{j\omega_{0}\tau}\overline{z}^{2}+\left(e^{j\omega_{0}\tau}+e^{-j\omega_{0}\tau}\right)z\overline{z} \\ &\quad+\bigg( \frac{{w}_{201}(0)}{2}e^{j\omega_{0}\tau}+{w}_{111}(0)e^{-j\omega_{0}\tau} +{w}_{111}(-\tau)+\frac{{w}_{201}(-\tau)}{2} \bigg)z^{2}\overline{z}+\cdots,\\ u_{1,t}(0)u_{2,t}(-\tau) &=\phi_1e^{-j\omega_{0}\tau}z^{2}+\overline{\phi}_1e^{j\omega_{0}\tau}\overline{z}^{2}+\left(\phi_1e^{-j\omega_{0}\tau}+\overline{\phi}_1e^{j\omega_{0}\tau}\right)z\overline{z} \\ &\quad+\bigg( \frac{{w}_{201}(0)}{2}\overline{\phi}_1e^{j\omega_{0}\tau}+{w}_{111}(0)\phi_1e^{-j\omega_{0}\tau}+{w}_{112}(-\tau)+\frac{{w}_{202}(-\tau)}{2} \bigg)z^{2}\overline{z}+\cdots, \\ u_{1,t}(-\tau)u_{2,t}(-\tau) &=\phi_1e^{-2j\omega_{0}\tau}z^{2}+\overline{\phi}_1e^{2j\omega_{0}\tau}\overline{z}^{2}+\left(\phi_1+\overline{\phi}_1\right)z\overline{z}+\bigg( \frac{{w}_{201}(-\tau)}{2}\overline{\phi}_1e^{j\omega_{0}\tau} \\ &\quad+{w}_{111}(-\tau)\phi_1e^{-j\omega_{0}\tau}+{w}_{112}(-\tau)e^{-j\omega_{0}\tau}+\frac{{w}_{202}(-\tau)}{2}e^{j\omega_{0}\tau} \bigg)z^{2}\overline{z}+\cdots, \\ u^3_{1,t}(0) &= 3z^{2}\overline{z}+\cdots, \\ u^2_{1,t}(0)u_{1,t}(-\tau) &= \big(e^{j\omega_{0}\tau}+2e^{-j\omega_{0}\tau}\big)z^{2}\overline{z}+\cdots, \\ u^2_{1,t}(0)u_{2,t}(-\tau)&=\big(2\phi_1e^{-j\omega_{0}\tau}+\overline{\phi}_{1}e^{j\omega_{0}\tau}\big)z^{2}\overline{z}+\cdots, \\ u_{1,t}(0)u_{1,t}(-\tau)u_{2,t}(-\tau)&=\big(\phi_1+ \overline{\phi}_{1}+\phi_1{}e^{-2j\omega_{0}\tau}\big)z^{2}\overline{z}+\cdots, \end{align*} where $$[{w}_{ij1}\,\,{w}_{ij2}]^{T}=\mathrm{\boldsymbol{w}}_{ij}$$. Recall that, $g(z,\overline{z})=\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}(z,\overline{z})\equiv \overline{B}\big(\overline{\phi}_2\mathcal{F}_{01}+\mathcal{F}_{02}\big),$ where $$\mathcal{F}_{0}=[\mathcal{F}_{01}\,\,\mathcal{F}_{02}]^{T}\,$$, and $g(z,\overline{z}) = g_{20}\frac{z^{2}}{2}+g_{11}z\overline{z}+g_{02}\frac{\overline{z}^{2}}{2}+g_{21}\frac{z^{2}\overline{z}}{2}\cdots.$ By comparing the coefficients of $$z^{2}$$, $$z\overline{z}$$, $$\overline{z}^{2}$$ and $$z^{2}\overline{z}$$, we get $$g_{20}$$, $$g_{11}$$, $$g_{02}$$ and $$g_{21}$$ as \begin{align} g_{20} &= 2\kappa{}\overline{B}\bigg(\overline{\phi}_2\Big(\xi_{xx}+\xi_{xx_d}e^{-j\omega_{0}\tau}+\xi_{xy_d}\phi_1e^{-j\omega_{0}\tau} + \xi_{x_dy_d}\phi_1e^{-2j\omega_{0}\tau}\Big)+ \frac{k_l}{\mathrm{w}_s}\phi_1\bigg), \label{eq:e37} \end{align} (A.16) \begin{align} g_{11} &= \kappa{}\overline{B}\biggl(\overline{\phi}_2\Big(2\xi_{xx}+\xi_{xx_d}\big(e^{j\omega_{0}\tau}+e^{-j\omega_{0}\tau}\big) \nonumber\\ &\quad +\xi_{xy_d}\big(\,\overline{\phi}_1e^{j\omega_{0}\tau}+\phi_1e^{-j\omega_{0}\tau}\big)+\xi_{x_dy_d}\big(\phi_1+\overline{\phi}_1\big)\Big) + \frac{k_l}{\mathrm{w}_s}\big(\phi_1+\overline{\phi}_1\big)\biggr), \label{eq:e38} \end{align} (A.17) \begin{align} g_{02} &= 2\kappa{}\overline{B}\bigg(\overline{\phi}_2\Big( \xi_{xx}+\xi_{xx_d}e^{j\omega_{0}\tau}+\xi_{xy_d}\overline{\phi}_1e^{j\omega_{0}\tau} +\xi_{x_dy_d}\overline{\phi}_1e^{2j\omega_{0}\tau}\Big)+\frac{k_l}{\mathrm{w}_s}\overline{\phi}_1 \bigg), \label{eq:e39} \end{align} (A.18) \begin{align} g_{21} &= \kappa{}\overline{B}\Bigg(\overline{\phi}_2\bigg[2\xi_{xx}\Big({w}_{201}(0)+2{w}_{111}(0)\Big) + \xi_{xx_d}\Big(2{w}_{111}(-\tau) + {w}_{201}(0)e^{j\omega_{0}\tau} + {w}_{201}(-\tau) \nonumber \\ &\quad + 2e^{-j\omega_{0}\tau}{w}_{111}(0) \Big) + \xi_{xy_d}\Big(2{w}_{111}(0)\phi_{1}e^{-j\omega_{0}\tau} + {w}_{202}(-\tau) + 2{w}_{112}(-\tau) + {w}_{201}(0)\overline{\phi}_{1}e^{j\omega_{0}\tau} \Big) \nonumber \\ &\quad + \xi_{x_dy_d}\Big( {w}_{201}(-\tau)\overline{\phi}_{1}e^{j\omega_{0}\tau}+2{w}_{111}(-\tau)\phi_{1}e^{-j\omega_{0}\tau} + 2{w}_{112}(-\tau)e^{-j\omega_{0}\tau}+{w}_{202}(-\tau)e^{j\omega_{0}\tau} \Big) \nonumber \\ &\quad +2\bigg(3\xi_{xxx} + \xi_{xxx_d}\Big(e^{j\omega_{0}\tau}+2e^{-j\omega_{0}\tau}\Big)+\xi_{xxy_d}\Big(2\phi_1e^{-j\omega_{0}\tau} + \overline{\phi}_{1}e^{j\omega_{0}\tau}\Big) \nonumber \\ &\quad + \xi_{xx_dy_d}\Big(2\text{Re}\big(\phi_1\big)+\phi_1{}e^{-2j\omega_{0}\tau}\Big)\bigg)\bigg] + \frac{k_l}{\mathrm{w}_s} \Big( {w}_{201}(0)\overline{\phi}_{1} + 2{w}_{111}(0)\phi_{1} + 2{w}_{112}(0)+{w}_{202}(0)\Big)\Biggr). \label{eq:e40} \end{align} (A.19) For the expression of $$g_{21}$$, we have to evaluate $$\mathrm{\boldsymbol{w}}_{11}(0)$$, $$\mathrm{\boldsymbol{w}}_{11}(-\tau)$$, $$\mathrm{\boldsymbol{w}}_{20}(0)$$ and $$\mathrm{\boldsymbol{w}}_{20}(-\tau)$$. Now for $$\theta$$$$\in$$$$[-\tau,0)$$ \begin{aligned} \boldsymbol{h}(z,\overline{z},\theta) &= -2\text{Re}\big(\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}\boldsymbol{q}(\theta)\big)\\ &= -2\text{Re}\big(g(z,\overline{z})\boldsymbol{q}(\theta)\big)\\ &= -\left(g_{20}\frac{z^{2}}{2}+g_{11}z\overline{z}+g_{02}\frac{\overline{z}^{2}}{2}+\cdots \right)\boldsymbol{q}(\theta) -\left(\overline{g}_{20}\frac{\overline{z}^{2}}{2}+\overline{g}_{11}z\overline{z}+\overline{g}_{02}\frac{z^{2}}{2}+\cdots\right)\overline{\boldsymbol{q}}(\theta), \end{aligned} which when compared with (A.14), yields \begin{aligned} \boldsymbol{h}_{20}(\theta)&=-g_{20}\boldsymbol{q}(\theta)-\overline{g}_{02}\overline{\boldsymbol{q}}(\theta)\\ \boldsymbol{h}_{11}(\theta)&=-g_{11}\boldsymbol{q}(\theta)-\overline{g}_{11}\overline{\boldsymbol{q}}(\theta). \end{aligned} From (A.6) and (A.15), we get \begin{align} \mathrm{\boldsymbol{w}}'_{20}(\theta) &= 2j\omega_{0}\mathrm{\boldsymbol{w}}_{20}(\theta) + g_{20}\boldsymbol{q}(\theta)+\overline{g}_{02}\overline{\boldsymbol{q}}(\theta) \label{eq:e41} \end{align} (A.20) \begin{align} \mathrm{\boldsymbol{w}}'_{11}(\theta) &= g_{11}\boldsymbol{q}(\theta)+\overline{g}_{11}\overline{\boldsymbol{q}}(\theta) \label{eq:e42}. \end{align} (A.21) Solving the differential (A.20) and (A.21), we obtain \begin{align} \mathrm{\boldsymbol{w}}_{20}(\theta) &= -\frac{g_{20}}{j\omega_{0}}\boldsymbol{q}_0e^{j\omega_{0}\theta}-\frac{\overline{g}_{02}}{3j\omega_{0}}\overline{\boldsymbol{q}}_{0}e^{-j\omega_{0}\theta} +\boldsymbol{E}e^{2j\omega_{0}\theta}\!\!\!\! \label{eq:e43} \end{align} (A.22) \begin{align} \mathrm{\boldsymbol{w}}_{11}(\theta) &= \frac{g_{11}}{j\omega_{0}}\boldsymbol{q}_0e^{j\omega_{0}\theta}-\frac{\overline{g}_{11}}{j\omega_{0}}\overline{\boldsymbol{q}}_{0}e^{-j\omega_{0}\theta} +\boldsymbol{F}, \label{eq:e44} \end{align} (A.23) for some $$\boldsymbol{E}=[e_{1}\,\,e_{2}]^{T}$$ and $$\boldsymbol{F}=[f_{1} \,\, f_{2}]^{T}$$, which we will now determine. For $$\boldsymbol{h}(z,\overline{z},0) = -2\text{Re}\big(\overline{\boldsymbol{q}}^{*}(0)\cdot\mathcal{F}_{0}q(0)\big)+\mathcal{F}_{0}$$, \begin{aligned} \boldsymbol{h}_{20}(0)=& -g_{20}\boldsymbol{q}(0)-\overline{g}_{02}\overline{\boldsymbol{q}}(0) +\left[\begin{array}{c}A_1\\A_2\end{array}\right]\\ \boldsymbol{h}_{11}(0)=& -g_{11}\boldsymbol{q}(0)-\overline{g}_{11}\overline{\boldsymbol{q}}(0) +\left[\begin{array}{c}B_1\\B_2\end{array}\right]\!, \end{aligned} where \begin{aligned} A_1 &= 2\kappa{}\Big(\xi_{xx}+\xi_{xx_d}e^{-j\omega_{0}\tau}+\xi_{xy_d}\phi_1e^{-j\omega_{0}\tau} + \xi_{x_dy_d}\phi_1e^{-2j\omega_{0}\tau}\Big), \\ B_1 &= 2\kappa{}\Big(\xi_{xx}+\text{Re}\big(\xi_{xx_d}e^{j\omega_{0}\tau}+\xi_{xy_d}\phi_1e^{-j\omega_{0}\tau} +\xi_{x_dy_d}\phi_1\big)\Big), \\ A_2 &= 2\kappa{}\phi_1 \frac{k_l}{\mathrm{w}_s}\,\,\, \text{and}\,\,\, B_2 = \kappa{} \big(\phi_1+\overline{\phi}_1\big) \frac{k_l}{\mathrm{w}_s}. \end{aligned} From (A.6) and (A.15), we get \begin{align} &g_{20}\boldsymbol{q}(0)+\overline{g}_{02}\overline{\boldsymbol{q}}(0) = \left[\begin{array}{c}A_1\\A_2\end{array}\right] + \left[ \begin{array}{c} \big(\kappa{} \xi_x-2j\omega_{0}\big){w}_{201}(0)+\kappa{}\xi_{y_d}{w}_{202}(-\tau) \\ \kappa{} \big(k_lp_s/\mathrm{w}_s\big) {w}_{201}(0)-2j\omega_{0}{w}_{202}(0) \end{array} \right]\!,\label{eq:e45} \end{align} (A.24) \begin{align} &g_{11}\boldsymbol{q}(0)+\overline{g}_{11}\overline{\boldsymbol{q}}(0) = \left[\begin{array}{c}B_1\\B_2\end{array}\right] +\left[ \begin{array}{c} \kappa \big(\xi_x {w}_{111}(0)+\xi_{y_d}{w}_{112}(-\tau)\big) \\ \kappa \big(k_lp_s/\mathrm{w}_s\big) {w}_{111}(0) \end{array} \right]\!. \label{eq:e46} \end{align} (A.25) Substituting the expression for $$\mathrm{\boldsymbol{w}}_{20}(0)$$, $$\mathrm{\boldsymbol{w}}_{20}(-\tau)$$, $$\mathrm{\boldsymbol{w}}_{11}(0)$$ and $$\mathrm{\boldsymbol{w}}_{11}(-\tau)$$ from (A.22) and (A.23) into (A.24) and (A.25), respectively, and solving for $$e_{1}$$, $$e_{2}$$, $$f_{1}$$ and $$f_{2}$$, we obtain \begin{align} e_1 &= \frac{1}{D}\Big(2j\omega_{0}A_1 + \kappa{}\xi_{y_d}A_2e^{-2j\omega_{0}\tau}\Big), \nonumber \\ e_2 &= \frac{1}{D}\Big(\kappa{} \frac{k_lp_s}{\mathrm{w}_s} A_1 - \big(\kappa{} \xi_x-2j\omega_{0}\big)A_2 \Big), \nonumber \\ f_1 &= \frac{-1}{p_s} \big(\phi_1 + \overline{\phi}_1 \big), \nonumber \\ f_2 &= \frac{-1}{\kappa{}\xi_{y_d}}\Big(\kappa{}\xi_x f_1 +B_1 \Big), \nonumber \end{align} where \begin{aligned} D &= -\bigg(2j\omega_{0}\big(\kappa{} \xi_x-2j\omega_{0}\big)+\kappa{}^2\xi_{y_d}\frac{k_lp_s}{\mathrm{w}_s} e^{-2j\omega_{0}\tau}\bigg). \end{aligned} Using the values of $$\boldsymbol{E}$$ and $$\boldsymbol{F}$$ in (A.22) and (A.23), followed by substituting $$\theta=0\,\, \text{and}\,-\tau$$, we can obtain the expressions for $$\mathrm{\boldsymbol{w}}_{20}(0)$$, $$\mathrm{\boldsymbol{w}}_{20}(-\tau)$$, $$\mathrm{\boldsymbol{w}}_{11}(0)$$ and $$\mathrm{\boldsymbol{w}}_{11}(-\tau)$$. With these expressions, we can now compute the term $$g_{21}$$. We now have the expressions for $$g_{20}$$, $$g_{11}$$, $$g_{02}$$ and $$g_{21}$$. All the necessary quantities that are essential for the Hopf bifurcation analysis have been computed. We now use the following expressions to analyse the type of the Hopf bifurcation (Hassard et al., 1981) $$\label{eq:e51} c_{1}(0) = \frac{j}{2\omega_{0}}\left(g_{20}g_{11}-2|g_{11}|^{2}-\frac{1}{3}|g_{02}|^{2}\right)+\frac{g_{21}}{2},$$ (A.26) \label{eq:e50} \begin{aligned} \mu_{2} &= \frac{ -\text{Re}\big(c_{1}(0)\big) }{ \alpha'(0) }, & \mathcal{B}_{2} &= 2\text{Re} \big( c_{1}(0) \big), \end{aligned} (A.27) where $$c_{1}(0)$$ is the lyapunov coefficient, $$\mathcal{B}_{2}$$ is the Floquet exponent, and $$\alpha'(0) = \text{Re}\big(d\lambda/d\kappa\big)\big|_{\kappa = \kappa_c}$$. We now state the conditions that will enable us to characterize the type of the Hopf bifurcation and determine the orbital stability of the bifurcating periodic solutions. The sign of $$\mu_2$$ determines the type of Hopf bifurcation. The Hopf bifurcation is super-critical if $$\mu_2 > 0$$ and sub-critical if $$\mu_2 < 0$$. The sign of (the Floquet exponent) $$\mathcal{B}_{2}$$ determines the stability of the bifurcating periodic solutions. The periodic solutions are asymptotically orbitally stable when $$\mathcal{B}_{2} < 0$$ and unstable if $$\mathcal{B}_{2} >0$$. Appendix B. Computations We conduct some numerical computations to demonstrate the Hopf bifurcation analysis. For a choice of system parameter values, using local stability analysis, one can determine the critical value of the bifurcation parameter at which the coupled system of C-TCP and E-RED (3.1) undergoes a local Hopf bifurcation. As an example, consider the following parameter values: $$\beta = 0.78$$, $$k=0.01$$, $$\gamma = 1$$, $$C = 5$$ units and $$\tau = 1$$ unit. For these parameter values, with $$\alpha = 1$$ the Hopf bifurcation occurs at $$k_l = 2.0029$$ and with $$\alpha = 4$$, we get $$k_l = 1.9248$$. Numerous system parameters, including the feedback delay, can act as the bifurcation parameter. However, we consider an exogenous parameter to be the bifurcation parameter. Our approach is to determine the Hopf condition using system parameters and then use the exogenous parameter to push the system just beyond the Hopf boundary. We do this as we are interested to know the impact of model parameters on the system behaviour in a local neighbourhood of the Hopf condition. At the Hopf boundary, by construction, the value of the exogenous non-dimensional parameter $$\kappa$$ is 1. An increase in $$\kappa$$, just beyond the Hopf condition, would render the system into a locally unstable regime. Using the Hopf bifurcation analysis, we compute the values of $$\mu_2$$ and $$\mathcal{B}_2,$$ which are presented in Fig. B.1. It can be seen that $$\mu_2 > 0$$ and $$\mathcal{B}_2 < 0$$, which implies that the Hopf bifurcation is supercritical, and the limit cycles are orbitally stable. Fig. B.1. View largeDownload slide Plots of $$\mu_2$$ and $$\beta_2$$ as the bifurcation parameter $$\kappa$$ varies. For the choice of $$\beta = 0.78$$, $$k=0.01$$, $$\gamma = 1$$, $$C = 5$$ units and $$\tau = 1$$ unit, system (3.1) undergoes a local Hopf bifurcation at $$\kappa = 1$$, for $$\alpha = 1, k_l = 2.0029$$ and $$\alpha = 4, k_l = 1.9248$$. Recall that $$\mu_2 > 0$$ implies that the Hopf bifurcation is supercritical, and $$\mathcal{B}_2 < 0$$ implies that the limit cycles are asymptotically orbitally stable. Fig. B.1. View largeDownload slide Plots of $$\mu_2$$ and $$\beta_2$$ as the bifurcation parameter $$\kappa$$ varies. For the choice of $$\beta = 0.78$$, $$k=0.01$$, $$\gamma = 1$$, $$C = 5$$ units and $$\tau = 1$$ unit, system (3.1) undergoes a local Hopf bifurcation at $$\kappa = 1$$, for $$\alpha = 1, k_l = 2.0029$$ and $$\alpha = 4, k_l = 1.9248$$. Recall that $$\mu_2 > 0$$ implies that the Hopf bifurcation is supercritical, and $$\mathcal{B}_2 < 0$$ implies that the limit cycles are asymptotically orbitally stable. The associated bifurcation diagrams and phase portraits are presented in Fig. B.2. The bifurcation diagrams, Fig. B.2(a and b), show the emergence of limit cycles in the average sending rate and the drop probability, at $$\kappa = 1$$. The phase portraits in Fig. B.2(c and d), for $$\kappa = 1.05 > \kappa_c$$, highlight the existence of stable limit cycles. The computations were performed using the software MATLAB. Using the analytical framework, one can perform similar computations for any desired set of parameter values. Fig. B.2. View largeDownload slide Bifurcation diagrams (a and b) and phase portraits (c and d), for the coupled system of C-TCP and E-RED (3.1). It can be observed from the bifurcation diagrams that the system transits from a stable equilibrium to a stable limit cycle. The phase portraits, plotted for $$\kappa = 1.05 > \kappa_c$$, establish the presence of stable limit cycles around the equilibrium ($$\mathrm{w}_s = 5$$ and $$p_s = 0.0495$$). Fig. B.2. View largeDownload slide Bifurcation diagrams (a and b) and phase portraits (c and d), for the coupled system of C-TCP and E-RED (3.1). It can be observed from the bifurcation diagrams that the system transits from a stable equilibrium to a stable limit cycle. The phase portraits, plotted for $$\kappa = 1.05 > \kappa_c$$, establish the presence of stable limit cycles around the equilibrium ($$\mathrm{w}_s = 5$$ and $$p_s = 0.0495$$). Acknowledgements The authors thank Nizar Malangadan and Santosh Chavan for their comments and feedback on this article. References Athuraliya, S., Low, S. H., Li, V. H. & Yin, Q. ( 2001 ) REM: active queue management. IEEE Netw., 15 , 48 – 53 . Google Scholar CrossRef Search ADS Briscoe, B., Brunstrom, A., Petlund, A., Hayes, D., Ros, D., Tsang, J., Gjessing, S., Fairhurst, G., Griwodz, C. & Welzl, M. ( 2016 ) Reducing internet latency: a survey of techniques and their merits. IEEE Commun. Surveys Tuts., 18 , 2149 – 2196 . Google Scholar CrossRef Search ADS Cao, J., Cleveland, W. S., Gao, Y., Jeffay, K. Smith, F. D. & Weigle, M. ( 2004 ) Stochastic models for generating synthetic HTTP source traffic. Proceedings of the IEEE INFOCOM. Hong Kong, China , pp. 1546 – 1557 . Floyd, S. & Jacobson, V. ( 1993 ) Random early detection gateways for congestion avoidance. IEEE/ACM Trans. Netw., 1 , 397 – 413 . Google Scholar CrossRef Search ADS Gettys, J. & Nichols, K. ( 2012 ) Bufferbloat: dark buffers in the internet. Commun. ACM, 55 , 57 – 65 . Google Scholar CrossRef Search ADS Gong, Y., Rossi, D., Testa, C., Valenti, S. & Täht, M. D. ( 2014 ) Fighting the bufferbloat: on the coexistence of AQM and low priority congestion control. Comput. Netw., 65 , 255 – 267 . Google Scholar CrossRef Search ADS Ha, S., Rhee, I. & Xu, L. ( 2008 ) CUBIC: a new TCP-friendly high-speed TCP variant. ACM SIGOPS Operating Systems Review, 42 , 64 – 74 . Google Scholar CrossRef Search ADS Han, H., Hollot, C. V., Towsley, D. & Chait, Y. ( 2005 ) Synchronization of TCP flows in networks with small Droptail buffers. Proceedings of the IEEE Conference on Decision and Control and European Control Conference. Seville, Spain , pp. 6762 – 6767 . Hassard, B. D., Kazarinoff, N. D. & Wan, Y.-H. ( 1981 ) Theory and Applications of Hopf Bifurcation. Cambridge, MA : Cambridge University Press . Hong, G., Martin, J. & Westall, J. M. ( 2015 ) On fairness and application performance of active queue management in broadband cable networks. Comput. Netw., 91 , 390 – 406 . Google Scholar CrossRef Search ADS Issariyakul, T. & Hossain, E. ( 2011 ) Introduction to Network Simulator NS2. New York : Springer . Kuznetsov, Y. A. ( 2004 ) Elements of Applied Bifurcation Theory. New York : Springer . Google Scholar CrossRef Search ADS Liu, S., Basar, T. & Srikant, R. ( 2005 ) Exponential-RED: a stabilizing AQM scheme for low and high-speed TCP protocols. IEEE/ACM Trans. Netw., 13 , 1068 – 1081 . Google Scholar CrossRef Search ADS Low, S. H. & Srikant, R. ( 2004 ) A mathematical framework for designing a low-loss, low-delay Internet. Netw. Spat. Econ., 4 , 75 – 101 . Google Scholar CrossRef Search ADS Malangadan, N., Rahman, H. & Raina, G. ( 2013 ) Non-linear oscillations in TCP networks with Drop-Tail buffers. Proceedings of the Chinese Control and Decision Conference, IEEE. Guiyang, China , pp. 188 – 194 . Nichols, K. & Jacobson, V. ( 2012 ) Controlling queue delay. Commun. ACM, 55 , 42 – 50 . Google Scholar CrossRef Search ADS Patil, G., McClean, S. & Raina, G. ( 2011 ) Drop tail and RED queue management with small buffers: stability and Hopf bifurcation. ICTACT J. Commun. Technol., 2 , 339 – 344 . Google Scholar CrossRef Search ADS Prasad, S. & Raina, G. ( 2014 ) Local Hopf bifurcation analysis of Compound TCP with an Exponential-RED queue management policy. Proceedings of the Chinese Control and Decision Conference, IEEE. Changsha, China , pp. 2588 – 2594 . Raja, P. & Raina, G. ( 2012 ) Delay and loss-based transport protocols: buffer-sizing and stability. Proceedings of the International Conference on Communication Systems and Networks, IEEE. Bangalore, India , pp. 1 – 10 . Raina, G., Towsley, D. & Wischik, D. ( 2005 ) Part II: control theory for buffer sizing. ACM SIGCOMM Comput. Commun. Rev., 35 , 79 – 82 . Google Scholar CrossRef Search ADS Raina, G., Manjunath, S., Prasad, S. & Giridhar, K. ( 2016 ) Stability and performance analysis of compound TCP With REM and drop-tail queue management. IEEE/ACM Trans. Netw., 24 , 1961 – 1974 . Google Scholar CrossRef Search ADS Showail, A., Jamshaid, K. & Shihada, B. ( 2016 ) Buffer sizing in wireless networks: challenges, solutions, and opportunities. IEEE Commun. Mag., 54 , 130 – 137 . Google Scholar CrossRef Search ADS Sojoudi, S., Low, S. H. & Doyle, J. C. ( 2014 ) Buffering dynamics and stability of internet congestion controllers. IEEE/ACM Trans. Netw., 22 , 1808 – 1818 . Google Scholar CrossRef Search ADS Srikant, R. ( 2004 ) The Mathematics of Internet Congestion Control. Cambridge, MA : Birkhauser . Google Scholar CrossRef Search ADS Tan, K., Song, J., Zhang, Q. & Sridharan, M. ( 2006 ) A Compound TCP approach for high-speed and long distance networks. Proceedings of the IEEE INFOCOM, IEEE. Barcelona, Spain , pp. 1 – 12 . Wischik, D. & McKeown, N. ( 2005 ) Part I: buffer sizes for core routers. ACM SIGCOMM Comput. Commun. Rev., 35 , 75 – 78 . Google Scholar CrossRef Search ADS Yang, P., Shao, J., Luo, W., Xu, L., Deogun, J. & Lu, Y. ( 2014 ) TCP congestion avoidance algorithm identification. IEEE/ACM Trans. Netw., 22 , 1311 – 1324 . Google Scholar CrossRef Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Mathematical Control and InformationOxford University Press

Published: Nov 25, 2017

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations