# Dynamics of a Discrete Internet Congestion Control Model

Authors: Yingguo Li
Journal: Discrete Dynamics in Nature and Society
 (2011)

---

## Abstract

We consider a discrete Internet model with a single link accessed by a single source, which responds to congestion signals from the network. Firstly, the stability of the equilibria of the system is investigated by analyzing the characteristic equation. By choosing the time delay as a bifurcation parameter, we prove that Neimark-Sacker bifurcations occur when the delay passes a sequence of critical values. Then, the explicit algorithm for determining the direction of the Neimark-Sacker bifurcations and the stability of the bifurcating periodic solutions is derived. Finally, some numerical simulations are given to verify the theoretical analysis.


---

## Body

## 1. Introduction

Congestion control algorithms and active queue management (AQM) for Internet have been the focus of intense research since the seminal work of Kelly et al. [1]. The primal algorithms they have introduced are based on explicit feedback functions, where the sum of user utilities is maximized within the capacity constraints of the links. In primal algorithms, the users adapt the source rates dynamically based on the congestion signal generated by the link, and the links select a static law to determine the link prices directly from the arrival rates at the links [2]. In [3], a stability condition was provided for a single proportionally fair congestion controller with delayed feedback. Since then, this result was extended to networks in [4, 5] and for different classes of controllers in [6].In this paper, we consider an Internet model with a single link and single source, which can be formulated as a congestion control system with feedback delay. The model is described by [3, 7, 8]
(1.1)x′(t)=k[w-x(t-τ)p(x(t-τ))],
where x(t) is the sending rate of the source at time t, k is a positive gain parameter, τ is the sum of forward and return delays, w>0 is a target (set point), and the congestion indication function p(·) is increasing, nonnegative, and not identically zero, which can be viewed as the probability that a packet at the resource receives a feedback congestion indication signal. System (1.1) corresponds to a rate control algorithm for the user, which can be used to adjust the sending rate x(t) so that the expected number of signals received by the user will tend to target w. A great deal of research has been devoted to the asymptotic behaviour, periodic solutions, global stability, and other properties of this model, for which, we refer to [3, 7–9].But considering the need of scientific computation and real-time simulation, our interest is focused on the behaviors of discrete dynamics system corresponding to (1.1). Many authors considered the numerical approximation of a scalar delay differential equation by using different numerical methods, such as θ-method, Euler method, and Trapezoidal method (see [10–15]). In this paper, we use the forward Euler scheme to make the discretization for system (1.1).Moreover, it is also of interest to find what will happen when the congestion control system loses stability. The purpose of this paper is to discuss this version as a discrete dynamical system by using Neimark-Sacker bifurcation theory of discrete systems. We not only investigate the stability of the fixed point and the existence of the Neimark-Sacker bifurcations, but also the stability and direction of the Neimark-Sacker bifurcation of the discrete system.The paper is organized as follows: in Section2, we analyze the distribution of the characteristic equation associated with the discrete model and obtain the existence of the local Neimark-Sacker bifurcation. In Section 3, the direction and stability of closed invariant curve from the Neimark-Sacker bifurcation of the discrete delay model are determined by using the theories of discrete systems in [16]. In the final section, some computer simulations are performed to illustrate the analytical results found.


## 2. Stability Analysis

Letu(t)=x(tτ). Then (1.1) can be rewritten as
(2.1)u′(t)=τk[w-u(t-1)p(u(t-1))].
We consider sufficient small step size of the form h=1/m, where m∈ℤ+. The Euler method applied to this equation yields the delay difference equation
(2.2)un+1=un+τkh[w-un-mp(un-m)],
where un is an approximate value to u(nh). Let u* be a positive equilibrium of (2.2), then we have(2.3)u*p(u*)=w.
Set yn=un-u*, then there follows:
(2.4)yn+1=yn+τkh[w-(yn-m+u*)p(yn-m+u*)].
By introducing a new variable Yn=(yn,yn-1,…,yn-m)T, we can rewrite (2.4) in the form
(2.5)Yn+1=F(Yn,τ),
where F=(F0,F1,…,Fm)T, and
(2.6)Fj={yn+τkh[w-(yn-m+u*)p(yn-m+u*)],j=0yn-j+1,1≤j≤m.
Clearly the origin is a fixed point of (2.5), and the linear part of (2.5) is(2.7)Yn+1=AYn,
where
(2.8)A=(10⋯0-τkh[p(u*)+u*p′(u*)]10⋯0001⋯00⋮⋮⋱⋮⋮00⋮10).
The characteristic equation of A is given by
(2.9)a(λ):=λm+1-λm+τkh[p(u*)+u*p′(u*)]=0.
It is well known that the stability of the zero equilibrium solution of (2.5) depends on the distribution of the zeros of the roots of (2.9). In order to proof the existence of the local Neimark-Sacker bifurcation at fixed point, we need some lemmas as follows.Lemma 2.1.
There exists aτ¯>0 such that for 0<τ<τ¯ all roots of (2.9) have modulus less than one.Proof.
Whenτ=0, (2.9) becomes
(2.10)λm+1-λm=0.
The equation has, at τ=0, an m-fold root λ=0 and a simple root λ=1.
Consider the rootλ(τ) such that λ(0)=1. This root depends continuously on τ and is a differential function of τ. From (2.9), we have
(2.11)dλdτ=kh[p(u*)+u*p′(u*)]mλm-1-(m+1)λm,dλ¯dτ=kh[p(u*)+u*p′(u*)]mλ¯m-1-(m+1)λ¯m.
Since p(·) is increasing and nonnegative, we have
(2.12)d|λ|2dτ|τ=0,λ=1=[λdλ¯dτ+λ¯dλdτ]|τ=0,λ=1=-2kh[p(u*)+u*p′(u*)]<0.
So with the increasing of τ>0, λ cannot cross λ=1. Consequently, all roots of (2.9) lie in the unit circle for sufficiently small positive τ>0, and the existence of the τ¯ follows.A Neimark-Sacker bifurcation occurs when a complex conjugate pair of eigenvalues ofA cross the unit circle as τ varies. We have to find values of τ such that there are roots on the unit circle. Denote the roots on the unit circle by eiω, ω∈(-π,π]. Since we are dealing with complex roots of a real polynomial, we only need to look for ω∈(0,π]. For ω∈(0,π], eiω is a root of (2.9) if and only if
(2.13)eiω-1+τkh[p(u*)+u*p′(u*)]e-imω=0.
Separating the real part and imaginary part from (2.13), there are
(2.14)cos⁡⁡ω+τkh[p(u*)+u*p′(u*)]cos⁡⁡mω=1,(2.15)sin⁡⁡ω-τkh[p(u*)+u*p′(u*)]sin⁡⁡mω=0.
So
(2.16)cos⁡⁡ω=1-12(τkh)2[p(u*)+u*p′(u*)]2.
From (2.14) and (2.15), we obtain that
(2.17)sin⁡⁡mωcos⁡⁡ω+cos⁡⁡mωsin⁡⁡ω=sin⁡⁡mω,
that is,
(2.18)sin⁡⁡(m+1)ω=sin⁡⁡mω,
then there exists a sequence of
(2.19)ωj∈(2jπm,(2j+1)πm),j=0,1,2,…,[m-12],
where [·] denotes the greatest integer function. It is clear that there exists a sequence of the time delay parameters τj satisfying (2.14) and (2.15), according to ω=ωj.From (2.15), we have
(2.20)τj=τ(ωj)=sin⁡⁡ωjkh[p(u*)+u*p′(u*)]sin⁡⁡mωj.
Differentiating (2.16) with respect to ω, we obtain that(2.21)sin⁡⁡ω=τ{kh[p(u*)+u*p′(u*)]}2dτdω.
This gives(2.22)dτdω|ω=ωj,τ=τj=sin⁡⁡ωjτj{kh[p(u*)+u*p′(u*)]}2>0.
Therefore,
(2.23)0<τ0<τ1<⋯<τ[(m-1)/2].Lemma 2.2.
Letλj(τ)=rj(τ)eiωj(τ) be a root of (2.9) near τ=τj satisfying rj(τj)=1 and ωj(τj)=ωj. Then
(2.24)drj2(τ)dτ|τ=τj,ω=ωj>0.Proof.
From (2.11), we have
(2.25)drj2dτ|τ=τj,ω=ωj=d|λj|2dτ|τ=τj,ω=ωj=[λdλ¯dτ+λ¯dλdτ]|τ=τj,ω=ωj=2kh[p(u*)+u*p′(u*)][mcos⁡⁡mωj-(m+1)cos⁡⁡(m+1)ωj]|mei(m-1)ωj-(m+1)eimωj|2.
From (2.14) and (2.15), we obtain that
(2.26)cos⁡⁡mωj=1-cos⁡⁡ωjτjkh[p(u*)+u*p′(u*)],sin⁡⁡mωj=sin⁡⁡ωjτjkh[p(u*)+u*p′(u*)].
It is easy to see that
(2.27)cos⁡⁡(m+1)ωj=cos⁡⁡mωjcos⁡⁡ωj-sin⁡⁡mωjsin⁡⁡ωj=cos⁡⁡ωj-1τkh[p(u*)+u*p′(u*)].
Using (2.25)–(2.27), we have
(2.28)d|λj|2dτ|τ=τj,ω=ωj=2(2m+1)(1-cos⁡⁡ωj)τj|mei(m-1)ωj-(m+1)eimωj|2>0.
This completes the proof.Lemma 2.3.
Equation (2.9) has a pair of simple roots e±iωj on the unit circle when τ=τj,j=0,1,2,…,[(m-1)/2]. Furthermore, if τ∈[0,τ0), then all the roots of (2.9) have modulus less than one; If τ=τ0, then all roots of (2.9) except e±iω0 have modulus less than one. But if τ∈(τj,τj+1], for j=0,1,2,…,[(m-1)/2], (2.9) has 2(j+1), roots have modulus more than one.Proof.
From (2.18)–(2.20), we have that (2.9) has toots e±iωj if and only if τ=τj and ω=ωj given in (2.19) and (2.20). By Lemmas 2.1 and 2.2, we know that if τ∈[0,τ0), then all roots of (2.9) have modulus less than one; if τ=τ0, then all roots of (2.9) except e±iωj have modulus less than one;  furthermore, by the Rouchs theorem (Dieudonné [17, Theorem 9.17.4]), the statement on the number of eigenvalues with modulus more than one is as follows.Lemma2.3 immediately leads to stability of the zero solution of (2.4), and equivalently, of the positive fixed point u=u* of (2.2). So we have the following results on stability and bifurcation in system (2.2).Theorem 2.4.
There exists a sequence of values of the time delay parameter0<τ0<τ1<⋯<τ[(m-1)/2] such that the positive fixed point u=u* of (2.2) is asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0. Equation (2.2) undergoes a Neimark-Sacker bifurcation at the u* when τ=τj,j=0,1,2,…,[(m-1)/2], where τj satisfies (2.20).


## 3. Direction and Stability of the Neimark-Sacker Bifurcation in Discrete Model

In the previous section, we obtain the conditions under which a family of periodic solutions bifurcate from the steady state at the critical valueτ=τj,j=0,1,2,…,[(m-1)/2]. Without loss of generality, denote the critical value τ=τj by τ*. For continuous systems, as pointed out in Hassard et al. [18], it is interesting and important to determine the direction, stability, and the period of these periodic solutions bifurcating from the steady state. In this section, following the idea of Hassard et al., we shall study the direction, stability, and the period of the bifurcating periodic solution when τ=τ0 in the discrete Internet congestion control model. The method we used is based on the theories of discrete system by Kuznetsov [16].Setτ=τ0+ν,ν∈ℝ. Then ν=0 is a Neimark-Sacker bifurcation value for (2.4). Rewrite (2.4) as
(3.1)yn+1=yn-τkh[(p(u*)+u*p′(u*))yn-m+12(2p′(u*)+u*p′′(u*))yn-m2+16(3p′′(u*)+u*p′′′(u*))yn-m3]+O(|yn-m|4).
So system (2.5) is turned into
(3.2)Yn+1=AYn+12B(Yn,Yn)+16C(Yn,Yn,Yn)+O(‖Yn‖4),
where(3.3)B(Yn,Yn)=(b0(Yn,Yn),0,…,0),C(Yn,Yn,Yn)=(c0(Yn,Yn,Yn),0,…,0),(3.4)b0(ϕ,ψ)=-τkh(2p′(u*)+u*p′′(u*))ϕmψm,c0(ϕ,,ψ,η)=-τkh(3p′′(u*)+u*p′′′(u*))ϕmψmηm.
Let q=q(τ0)∈ℂm+1 be an eigenvector of A corresponding to eiω0, then
(3.5)Aq=eiω0q,Aq¯=e-iω0q¯.
We also introduce an adjoint eigenvector q*=q*(τ)∈ℂm+1 having the properties
(3.6)ATq*=e-iω0q*,ATq¯*=eiω0q¯*,
and satisfying the normalization 〈q*,q〉=1, where 〈q*,q〉=∑i=0mq¯i*qi.Lemma 3.1 (see [19]).
Define a vector-valued functionq:ℂ→ℂm+1 by
(3.7)p(ξ)=(ξm,ξm-1,…,1)T.
If ξ is an eigenvalue of A, then Ap(ξ)=ξp(ξ).In view of Lemma3.1, we have
(3.8)q=p(eiω0)=(eimω0,ei(m-1)ω0,…,eiω0,1)T.Lemma 3.2.
Suppose thatq*=(q0*,q1*,…,qm*)T is the eigenvector of AT corresponding to eigenvalue e-iω0, and 〈q*,q〉=1. Then
(3.9)q*=K¯(1,a0eimω0,a0ei(m-1)ω0,…,a0ei2ω0,a0eiω0)T,
where a0=-τkh(2p′(u*)+u*p′′(u*)) and
(3.10)K=(eimω0+ma0e-iω0)-1.Proof.
Assignq* satisfies ATq*=z¯q* with z¯=e-iω0, then the following identities hold
(3.11)q0*+q1*=e-iω0q0*,qk*=e-iω0qk-1*,k=2,…,m,a0q0*=e-iω0qm*.
Let qm*=a0eiω0K¯, then
(3.12)q*=K¯(1,a0eimω0,a0ei(m-1)ω0,…,a0ei2ω0,a0eiω0)T.
From normalization 〈q*,q〉=1 and computation, (3.10) is held.Leta(λ) be characteristic polynomial of A and λ=eiω0. Following the algorithms in [16] and using a computation process similar to that in [19], we can compute an expression for the critical coefficient c1(τ0),
(3.13)c1(τ0)=g20g11(1-2λ0)2(λ02-λ0)+|g11|21-λ¯0+|g02|22(λ02-λ¯0)+g212,
where
(3.14)g20=〈q*,B(q,q)〉,g11=〈q*,B(q,q¯)〉,g02=〈q*,B(q¯,q¯)〉,g21=〈q*,B(q¯,ω20)〉+2〈q*,B(q,ω11)〉+〈q*,C(q,q,q¯)〉,ω20=b0(q,q)a(λ02)p(λ02)-〈q*,B(q,q)〉λ02-λ0q-〈q¯*,B(q,q)〉λ02-λ¯0q¯,ω11=b0(q,q¯)a(1)p(1)-〈q*,B(q,q¯)〉1-λ0q-〈q¯*,B(q,q¯)〉1-λ¯0q¯.
By (3.4), (3.8), and Lemma 3.2, we get
(3.15)b0(q¯,p(ei2ω0))=-τkh(2p′(u*)+u*p′′(u*)),b0(q,q)=-τkh(2p′(u*)+u*p′′(u*)),b0(q,q¯)=-τkh(2p′(u*)+u*p′′(u*)),c0(q,q,q¯)=-τkh(3p′′(u*)+u*p′′′(u*)),a(ei2ω0)=ei2(m+1)ω0-ei2mω0+τkh(2p′(u*)+u*p′′(u*)),a(1)=τkh(2p′(u*)+u*p′′(u*)),b0(q,p(1))=-τkh(2p′(u*)+u*p′′(u*)).
Substituting these into (3.13), we have
(3.16)c1(τ0)=K2{[τkh(2p′(u*)+u*p′′(u*))]2a(ei2ω0)+2[τkh(2p′(u*)+u*p′′(u*))]2a(1)-τkh(3p′′(u*)+u*p′′′(u*))[τkh(2p′(u*)+u*p′′(u*))]2a(ei2ω0)}.Lemma 3.3 (see [20]).
Given the map (2.5) and assume that(1)
λ(τ)=r(τ)eiω(τ), where r(τ*)=1,r′(τ*)≠0 and ω(τ*)=ω*;(2)
eikω*≠1 for k=1,2,3,4;(3)
Re[e-iω*c1(τ*)]≠0.Then an invariant closed curve, topologically equivalent to a circle, for map (2.5) exists for τ in a one-side neighborhood of τ*. The radius of the invariant curve grows like O(|τ-τ*|). One of the four cases below applies:(1)
r′(τ*)>0,Re[e-iω*c1(τ*)]<0. The origin is asymptotically stable for τ<τ* and unstable for τ>τ*. An attracting invariant closed curve exists for τ>τ*.(2)
r′(τ*)>0,Re[e-iω*c1(τ*)]>0. The origin is asymptotically stable for τ<τ* and unstable for τ>τ*. An repelling invariant closed curve exists for τ<τ*.(3)
r′(τ*)<0,Re[e-iω*c1(τ*)]<0. The origin is asymptotically stable for τ>τ* and unstable for τ<τ*. An attracting invariant closed curve exists for τ<τ*.(4)
r′(τ*)<0,Re[e-iω*c1(τ*)]>0. The origin is asymptotically stable for τ>τ* and unstable for τ<τ*. An repelling invariant closed curve exists for τ>τ*.From the discussion in Section2, we know that r′(τ*)>0; therefore, by Lemma 3.3, we have the following result.Theorem 3.4.
For (2.2), u=u* is asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0. An attracting (repelling) invariant closed curve exists for τ>τ0 if Re[e-iω0c1(τ0)]<0(>0).


## 4. Computer Simulation

In this section, we will confirm our theoretical analysis by numerical simulation. We choosep(x)=x/(20-3x) as in [8] and set k=1.5,m=4,w=1. Then τ0=1.0159 is the Neimark-Sacker bifurcation value.Figures1 and 2 are about delay difference equation (2.2) when step size h=0.25.(a) The equilibriumu* of (2.2) is asymptotically stable for τ=0.9<τ0=1.0159, (b) Phase plot of (2.2) on the plane (u(n),u(n-4)) for τ=0.9<τ0=1.0159.
(a)(b)(a) Waveform plot of (2.2) for τ=1.1>τ0=1.0159, (b) A bifurcating periodic solution appears on the plane (u(n),u(n-4)) for τ=1.1>τ0=1.0159.
(a)(b)In Figure1, we show the waveform plot and phase plot for (2.2) with initial values uj=3(j=0,1,…,4) for τ=0.9<τ0=1.0159. The equilibrium u*=3.2170 of (2.2) is asymptotically stable. In Figure 2, we show the waveform plot and phase plot for (2.2) with initial values uj=3(j=0,1,…,4). The equilibrium u*=3.2170 of (2.2) is unstable for τ=1.1>τ0=1.0159. When τ varies and passes through τ0=1.0159, the equilibrium loses its stability and a periodic solution bifurcates from the equilibrium for τ=1.1>τ0=1.0159. That is the delay difference equation (2.2) which has a Neimark-Sacker bifurcation at τ0.



---

*Source: 628369-2011-11-28.xml*