# A Novel Four-Dimensional Energy-Saving and Emission-Reduction System and Its Linear Feedback Control

Authors: Minggang Wang; Hua Xu
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

This paper reports a new four-dimensional energy-saving and emission-reduction chaotic system. The system is obtained in accordance with the complicated relationship between energy saving and emission reduction, carbon emission, economic growth, and new energy development. The dynamics behavior of the system will be analyzed by means of Lyapunov exponents and equilibrium points. Linear feedback control methods are used to suppress chaos to unstable equilibrium. Numerical simulations are presented to show these results.


---

## Body

## 1. Introduction

Since energy saving and emission reduction is the most effective way to control carbon emissions, how to promote energy saving and emission reduction is becoming the hot topic of academic research. Calculation and control analysis of carbon emissions have attracted a great deal of attention from various fields of researchers. Feng et al. [1] conducted a research on the long-run equilibrium relationships, temporal dynamic relationships, and causal relationships between energy consumption structure, economic structure, and energy intensity in China. In order to decrease energy intensity, the Chinese government should continue to reduce the proportion of coal in energy consumption, increase the utilization efficiency of coal, and promote the upgrade of economic structure. Amjad et al. [2] indicated that petroleum is the major energy consumption of most of the nations in the world, so taking actions to reduce the petroleum consumption, such as replacing the diesel locomotive with hybrid electrical vehicle, could reduce humans dependence on petroleum and hence decrease carbon emissions. Liao [3] analyzed the role of developing hydro-energy, wind energy, nuclear power, and so forth. Guo et al. [4] found technical innovations and standard management a decisive role in the energy use per unit of GDP, for which the ratio of oil consumption was the major limiting factor. Mendiluce et al. [5] compared the evolution of energy intensity in Spain with that in the EU15 and found the increasing of energy intensity in Spain since 1990 is mainly due to strong transport growth and the construction boom.From the above analysis, we can see that previous researches mainly explored the influence of energy intensity by energy structures, technical change and management level, energy consumption structure, economic structure, energy prices, and so on [4, 6]. Some were proceeded from the world or state perspectives [4–6], some from provincial and local governments perspectives [7, 8], which figured out the variables which influence energy intensity and came up with the corresponding measures to reduce energy intensity. The research achievements are satisfying. Compared with the previous researches, this study was undertaken from continuous dynamic equation, making clear the quantitative relationships among the concerning variables. It brought energy saving and emission reduction, carbon emissions, economic growth, and new energy development into a nonlinear dynamics system with the analysis of the relationship between the variables and their influence on energy intensity. With the aid of simulation figures, the evolution behavior and the change regularity of the four-dimension system, and their influence trends on energy intensity are shown vividly. It is clear that this paper is more vivid and more adherent to the reality.Chaos analysis and applications in dynamical systems are observed in many practical applications in engineering, biology, and economics [9–13]. Energy-saving and emission-reduction system is a complex nonlinear system, which includes energy-saving and emission-reduction, carbon emissions, economic growth, energy efficiency, carbon tax, energy intensity, and so forth [11–13]. One of the most noticeable problems is how to conduct a further research of energy saving and emission reduction through nonlinear dynamics, which is currently a method of rapid development. While most previous studies focused on scenario analysis, Fang et al. [14] established a three-dimensional system in accordance with the complicated relationship between energy saving and emission reduction, carbon emissions and economic growth. This system displays a very complex phenomenon and contains a special chaotic attractor named the energy-saving and emission-reduction attractor, which is different from the previous chaotic attractor, such as Lorenz attractor [15], Chen attractor [16], Lü attractor [17], Energy resource attractor [18–20], and so forth.In the three-dimensional energy-saving and emission-reduction system, the authors have not considered clean energy development (including wind energy, solar energy, hydropower, geothermal, biomass energy, and so forth), but most of the countries are developing and making use of new energy resources. Therefore, it is necessary to add new energy resources to the three-dimensional energy-saving and emission-reduction system.By adding a new variable (new energy development) to the three-dimensional energy-saving and emission-reduction system, a new four-dimensional energy-saving and emission-reduction system is obtained.This paper establishes a new four-dimensional energy-saving and emission-reduction system. It is organized as follows: Section2 sets up the model; Section 3 discusses basic properties of the system and gives numerical results. Simulation results show that the system can generate complex chaotic attractors when the system parameters are chosen appropriately. Linear feedback control criterions are presented in Section 4. Conclusions are finally given in Section 5.


## 2. Establishment of the Model

By adding a new variableu(t): the time-dependent variable of new energy development during a given economic period to the three-dimensional energy-saving and emission-reduction system, a new four-dimensional energy-saving and emission-reduction system is obtained as follows:
(2.1)x˙=a1x(yM-1)-a2y+a3z,y˙=-b1x+b2y(1-yC)+b3z(1-zE)-d4u,z˙=c1x(xN-1)-c2y-c3z+c4u(uL-1),u˙=d1y+d2z(zK-1)-d3u,
where x(t) is the time-dependent variable of energy saving and emission reduction, y(t) is the time-dependent variable of carbon emissions, and z(t) is the time-dependent variable of economic growth (GDP).  ai,  bi,  cj,  dj,  M,  N,  L,  K  are positive constants, t∈I, I is a given economic period (i=1,2,3,j=1,2,3,4, the units of M,N,L,K can be transformed into tons of standard coal). a1 is the development coefficient of x(t), a2 is the influence coefficient of y(t) to x(t), a3 is the impudence coefficient of z(t) to x(t), M is the inflexion (local maximum point) of y(t) to x(t); b1 is the influence coefficient of x(t) to y(t), b2 is the development coefficient of y(t), b3 is the influence coefficient of z(t) to y(t), C is the peak value of y(t) during a given period, E is the peak value of z(t) during a given period, d4 is the influence coefficient of u(t) to y(t), c1 is the influence coefficient of x(t) to z(t), c2 is the influence coefficient of y(t) to z(t), c3 is the influence coefficient of x(t) to z(t), N is the inflexion of x(t) to z(t), c4 is the influence coefficient of u(t) to z(t), and  L is the inflexion of u(t) to z(t). d1 is the influence coefficient of y(t) to u(t), d2 is the influence coefficient of z(t) to u(t), K is the inflexion of z(t) to u(t), and d3 is the influence coefficient of to itself.The first formula in (2.1) expresses the complicated relationship between the change rate of time-dependent energy saving and emission reduction d(x)/d(t), energy saving and emission reduction, carbon emissions and economic growth during a given period, which indicates that the change rate of time-dependent energy saving and emission reduction d(x)/d(t) is associated with energy saving and emission reduction x(t) and the share of energy saving and emission reduction potential y/M-1 simultaneously, in a positive proportion to them. As for a1x(y/M-1), when y<M, that is,  Y/M-1<0, the development trend of x(t) becomes weaker; when y>M, the development trend of x(t) becomes faster. d(x)/d(t) is inversely proportional to carbon emissions y(t), that is, the accession of y(t) will counteract the change rate of d(x)/d(t). d(x)/d(t) is positively proportional to economic growth z(t), that is,  the increasing investment in x(t) will promote the growth of d(x)/d(t).The second formula in (2.1) indicates that the change rate of time-dependent carbon emissions dy/dt is positively proportional to x(t), that is,  the development of x(t) will slowdown the pace of dy/dt. The development speed of y(t) is fast before the peak value C and will slow down after the peak value. The early stage of development of z(t) will bring about much carbon emissions, the influence of which on y(t) will become moderate after the peak value E. As for b2y(1-y/C), when y<C, that is,  1-y/C>0, the development speed of y(t) is fast; when y>C, the development trend of y(t) gets weaker. As for b3z(1-z/E), when z<E, that is,  1-z/E>0, the influence of z(t) for y(t) is positive; when z(t) arrives at the peak value E, the influence of z(t) on y(t) will be negative. dy/dt is inversely proportional to new energy development u(t), that is, the accession of u(t) will counteract the change rate of dy/dt.The third formula in (2.1) indicates that the early investment to x(t) will counteract the development of z(t). With technology progress and integrated development of x(t), x(t) will promote z(t) in turn. As for c1x(x/N-1), when x<N, that is, x/N-1<0, the influence of x(t) on z(t) is negative; when x>N, the influence of x(t) on z(t) is positive. The change rate of time-dependent economic growth dz/dt is inversely proportional to y(t), that is, the accession of y(t) will counteract the development of z(t). dz/dt is inversely proportional to investment to energy-saving and emission-reduction, that is,  the investment will counteract the development of z(t) to a certain extent. The early investment to u(t) will counteract the development of z(t). With the progress of u(t), u(t) will promote z(t) in turn. As for c4u(u/L-1), when u<L, that is,  u/L-1<0, the influence of u(t)  on z(t) is negative; when u>L, the influence of u(t) on z(t) is positive.The fourth formula in (2.1) indicates that du/dt is positively proportional to carbon emissions y(t), that is, the increasing carbon emissions will promote the growth of du/dt. The low level of z(t) will counteract the development of u(t). With the development of z(t), z(t) will promote u(t) in turn. As for d2z(z/K-1), when z<K, that is, z/K-1<0, the influence of z(t) on u(t) is negative. The change rate of time-dependent new energy development du/dt is inversely proportional to u(t), that is, du/dt will decrease with the increase of u(t).


## 3. The Basic Properties Analysis of the Four-Dimensional Energy-Saving and Emission-Reduction System

### 3.1. Equilibrium Point

We can obtain that the system (2.1) has four equilibriums: O(0,0,0,0), S1(x1,y1,z1,u1), S2(x2,y2,z2,u2), and S3(x3,y3,z3,u3). Linearizing the system (2.1) at equilibrium O(0,0,0,0) yields the Jacobian matrix
(3.1)(a1yM-a1a1xM-a2a30-b1b2-2b2yCb3-2b3zE-d42c1xN-c1-c2-c32c4uL-c40d12d2zK-d2-d3).For simplicity, we fix the following parameters:a1=0.09, a2=0.003, a3=0.012, b1=0.0412, b2=0.08, b3=0.8, c1=0.035, c2=0.0062, c3=0.08, c4=0.02, d1=0.01, d2=0.02, d3=0.06, d4=0.03, M=0.9, C=1.6, E=2.8, N=0.35, K=2, and L=2. By calculations, we can obtain that the eigenvalues of the Jacobian matrix of the system (2.1) at O(0,0,0,0) are
(3.2)λ1=-0.0103,λ2=0.0380,λ3=-0.0899+0.0216i,λ4=-0.0899-0.0216i.
Therefore O(0,0,0,0) is an unstable saddle focus.Theorem 3.1.
(
1
)  If 5.9986<d1<9.8401, the equilibrium point O(0,0,0,0) is stable. (2) If 0<d1≤5.9986 or d1≥9.8401, the equilibrium point O(0,0,0,0) is unstable.Proof.
We fix the following parameters:a1=0.09, a2=0.003, a3=0.012, b1=0.0412, b2=0.08, b3=0.8, c1=0.035, c2=0.0062, c3=0.08, c4=0.02,  d2=0.02, d3=0.06, d4=0.03, M=0.9, C=1.6, E=2.8, N=0.35, K=2, and L=2. While we let parameter d1 be varied, and the corresponding characteristic equation of Jacobian matrix of the system (2.1) at O(0,0,0,0) is
(3.3)f(λ)=(λ+0.0489)(λ3+0.1011λ2+0.0057λ+0.00015d1-0.00089979).
Solving (3.3) gives λ1=-0.0489<0, and the following equation:
(3.4)λ3+0.1011λ2+0.0057λ+0.00015d1-0.00089979=0.
Letp1=0.1011,  p2=0.0057, and  p3=0.00015d1-0.00089979. By the Routh-Hurwitz criterion, all real eigenvalues and all real parts of complex conjugate eigenvalues of (3.4) are negative if and only if the following conditions hold:
(3.5)p1>0,p3>0,p1p2-p3>0.
That is, 5.9986<d1<9.8401. Therefore, when 5.9986<d1<9.8401, the equilibrium point O(0,0,0,0) is stable; when 0<d1≤5.9986 or d1≥9.8401, the equilibrium point O(0,0,0,0) is unstable.We fix parameters as above, and then obtain the equilibrium pointS1(1.3748,0.7691,1.6915,0.0412), S2(0.8733,0.8609,0.4999,0.0185), and S3(-1.5441,1.1453,3.4423,1.0183). By calculations, we can obtain that the eigenvalues of the Jacobian matrix of system (2.1) at S1 are λ1=-0.1516, λ2=-0.0644, λ3,4=0.0330±0.2482i; the eigenvalues of the Jacobian matrix of system (2.1) at S2 are λ1=0.0691, λ2=-0.0634, λ3,4=-0.0778±0.2044i; the eigenvalues of the Jacobian matrix of system (2.1) at S3 are λ1=-0.4812, λ2=-0.0613, λ3,4=0.1963±0.2874i. Therefore, S1, S2, and S3 are three saddle points.


### 3.2. Dissipation

Consider the following:(3.6)∇V=∂x˙∂x+∂y˙∂y+∂z˙∂z+∂u˙∂u=a1yM-a1+b2-2b2yC-c3-d3=(a1M-2b2C)y+(b2-a1-c3-d3).Ifa1/M=2b2/C and b2-a1-c3-d3<0, then the system (2.1) is a dissipative system.


### 3.3. Theoretical Proof of the Existence of Smale Horseshoes and the Horseshoes Chaos

To study the long-term dynamical behavior of the system (2.1), the system is divided into subsystems. Let u=0; then the first subsystem is obtained:
(3.7)x˙=a1x(yM-1)-a2y+a3z,y˙=-b1x+b2y(1-yC)+b3z(1-zE),z˙=c1x(xN-1)-c2y-c3z.Letz is a known function of the time t; then the second subsystem is obtained:
(3.8)x˙=a1x(yM-1)-a2y+a3z,y˙=-b1x+b2y(1-yC)+b3z(1-zE)-d4u,u˙=d1y+d2z(zK-1)-d3u.
When t=t0, z is a constant number, then the system (3.8) is a three-dimensional nonlinear system.Theorem 3.2.
The subsystem presented in (3.7) exhibits horseshoe chaos.Proof.
Equation (3.7) has four equilibrium points:
(3.9)O(0,0,0,0),S1(x1,y1,z1,u1),S2(x2,y2,z2,u2),S3(x3,y3,z3,u3).
Step  1. S1 is saddle foci, that is, the eigenvalues of the real matrix A=Df(S1) and the Jacobin derivative of f at S1 are the forms: λ1=r,  λ2,3=σ±iω,  r<0,  σ>0,  |r|>σ, where r,  σ,  ω are real. Here the equilibrium point S1 is discussed with the Jacobin matrix
(3.10)J1=(a1yM-a1a1xM-a2a3-b1b2-2b2yCb3-2b3zE2c1xN-c1-c2-c3).
The characteristic polynomial is obtained as(3.11)det(λI-J1)=λ3+q1λ2+q2λ+q3=0.
Supposeλ=η-q1/3, (3.11) produces η3+qη+p=0, where p=q3-q1q2/3+2q13/27, q=q2-q12/3. Then turn to the Cardano formula, ensure the roots of  (3.11) meet the requirement of  Šilnikov theorem. It yields Δ=(p/2)2+(q/3)2<min(-2q1/3,q1/3).
When parameters are fixed as above atS1(-1.519,1.151,3.460), the following can be obtained: λ1=-0.4786,  λ2,3=0.1943±0.2848i, which satisfy r,  σ±iω,  r<0,  σ>0,  |r|>σ.
Step  2. There exists a homoclinic orbit τ1 at S1.
For the discussion of the homoclinic orbit ofS1, it can be supposed as series form like
(3.12)x(t)=l0+∑k=1+∞lkekαt,y(t)=m0+∑k=1+∞mkekαt,z(t)=n0+∑k=1+∞nkekαt,
where lk,mk,nk(k≥1) are undetermined coefficients, α is attenuation index, when t→+∞, (x(t),y(t),z(t))→(x1,y1,z1).
Next, substitute (3.12) into (2.1), and match constants of items coefficient
(3.13)a1l0m0M-a1l0-a2m0+a3n0=0,-b1l0+b2m0-b2m02C+b3n0-b3n02E=0,c1l02N-c1l0-c2m0-c3n0=0.
The following equation can be obtained:(l0,m0,n0)=(x1,y1,z1). Comparing coefficients of ekαt of the same power terms, the following is obtained:
(3.14)(kαI-J1)(lkmknk)=(φk(1)(ai,bi,ci,α,ξ)φk(2)(ai,bi,ci,α,ξ)φk(3)(ai,bi,ci,α,ξ)).
If(l1,m1,n1)=(0,0,0), then (lk,mk,nk)=(0,0,0), (k>1); therefore, (l1,m1,n1)≠(0,0,0). J1 has the only negative eigenvalues, and then the only α(α<0) satisfies det(αI-J1)=0. Note that det(kαI-J1)≠0 (α is the only negative real eigenvalues of J1, so kα is not the eigenvalues of J1); therefore, lk,mk,nk can be identified uniquely and x(t),  y(t),  z(t) for t>0. For the opposite time symmetric track x(t),  y(t),  z(t), the linear transform can be adopted τ=-t,  t>0. The proof is similar to the procedure as t>0.
According to the above analysis, homoclinic orbitτ1 based on equilibrium point S1 is formally obtained as
(3.15)x(t)={∑k=1+∞lk(ai,bi,ci,α,ξ)ekαt,t>0,∑k=1+∞lk(ai,bi,ci,-α,ξ)e-kαt,t<0,y(t)={∑k=1+∞mk(ai,bi,ci,α,ξ)ekαt,t>0,∑k=1+∞mk(ai,bi,ci,-α,ξ)e-kαt,t<0,z(t)={∑k=1+∞nk(ai,bi,ci,α,ξ)ekαt,t>0,∑k=1+∞nk(ai,bi,ci,-α,ξ)e-kαt,t<0.lk,mk,nk,(k>1) are determined by (3.14), while α is determined by det(αI-J1)=0, and ξ is determined by ∑k=1+∞lk(ai,bi,ciα,ξ)=∑k=1+∞lk(ai,bi,ci,-α,ξ).
Step  3. S3 is saddle foci, and there exists a homoclinic orbit τ2 based at S3. The proof is similar to Steps 1 and 2. So, by the Šilnikov theorem, the horseshoe chaos may expect to occur in a reasonable regime.Theorem 3.3.
The subsystem presented in (3.8) exhibits horseshoe chaos.Proof.
The proof is similar to the proof of Theorem3.2.Remark 3.4.
Theorems3.2 and 3.3 show that the energy-saving and emission-reduction system exhibits horseshoe chaos. This is not to say that the system is chaotic all the time, which means that the system is chaotic under appropriate conditions, and stable under other appropriate conditions.


### 3.4. Numerical Results

We chose a set of parameters as follows:a1=0.09, a2=0.003, a3=0.012, b1=0.0412, b2=0.08, b3=0.8, c1=0.035, c2=0.0062, c3=0.08, c4=0.02, d1=0.01, d2=0.02, d3=0.06, d4=0.03, M=0.9, C=1.6, E=2.8, N=0.35, K=2, and L=2. Let initial condition be (0.015, 0.785, 1.83, 0.01), and the corresponding Lyapunov exponents are L1=0.0272>0, L2=-0.0011<0, L3=-0.0839<0, and L4=-0.0772<0. Therefore, the Lyapunov dimension of this system is
(3.16)Dl=j+1|Lj+1|∑i=1jLi=3+L1+L2+L3|L4|=2.2522,
which means that the Lyapunov dimension is fractional under the same condition. The system has a chaotic attractor, as shown in Figures 1(a)–1(c), the time series of x(t),y(t),z(t),u(t) as shown in Figure 1(d). Let d1=0.094, fixed other parameters and initial condition, and then produce a limit cycle as shown in Figure 2. By calculations, the Lyapunov exponent spectrum with respect to parameter d1 is shown in Figure 3. According to Figure 3, the system (2.1) has very rich dynamical behaviors, which are summarized as follows. When d1∈(0.0935,0.095), the system (2.1) is chaotic with a positive Lyapunov exponent (e.g., with d1=0.01, the phase portrait is shown in Figure 1); while for d1∈(0.0935,0.095), the maximum Lyapunov exponent equals zero, implying that the system has a periodic orbit (Figure 2 shows the periodic orbit when d1=0.094).A four-dimensional energy-saving and emission-reduction chaotic attractor.
(a)
(b)
(c)
(d)A limit cycle.
(a)
(b)
(c)
(d)Figure 3
Lyapunov exponent spectrum.


## 3.1. Equilibrium Point

We can obtain that the system (2.1) has four equilibriums: O(0,0,0,0), S1(x1,y1,z1,u1), S2(x2,y2,z2,u2), and S3(x3,y3,z3,u3). Linearizing the system (2.1) at equilibrium O(0,0,0,0) yields the Jacobian matrix
(3.1)(a1yM-a1a1xM-a2a30-b1b2-2b2yCb3-2b3zE-d42c1xN-c1-c2-c32c4uL-c40d12d2zK-d2-d3).For simplicity, we fix the following parameters:a1=0.09, a2=0.003, a3=0.012, b1=0.0412, b2=0.08, b3=0.8, c1=0.035, c2=0.0062, c3=0.08, c4=0.02, d1=0.01, d2=0.02, d3=0.06, d4=0.03, M=0.9, C=1.6, E=2.8, N=0.35, K=2, and L=2. By calculations, we can obtain that the eigenvalues of the Jacobian matrix of the system (2.1) at O(0,0,0,0) are
(3.2)λ1=-0.0103,λ2=0.0380,λ3=-0.0899+0.0216i,λ4=-0.0899-0.0216i.
Therefore O(0,0,0,0) is an unstable saddle focus.Theorem 3.1.
(
1
)  If 5.9986<d1<9.8401, the equilibrium point O(0,0,0,0) is stable. (2) If 0<d1≤5.9986 or d1≥9.8401, the equilibrium point O(0,0,0,0) is unstable.Proof.
We fix the following parameters:a1=0.09, a2=0.003, a3=0.012, b1=0.0412, b2=0.08, b3=0.8, c1=0.035, c2=0.0062, c3=0.08, c4=0.02,  d2=0.02, d3=0.06, d4=0.03, M=0.9, C=1.6, E=2.8, N=0.35, K=2, and L=2. While we let parameter d1 be varied, and the corresponding characteristic equation of Jacobian matrix of the system (2.1) at O(0,0,0,0) is
(3.3)f(λ)=(λ+0.0489)(λ3+0.1011λ2+0.0057λ+0.00015d1-0.00089979).
Solving (3.3) gives λ1=-0.0489<0, and the following equation:
(3.4)λ3+0.1011λ2+0.0057λ+0.00015d1-0.00089979=0.
Letp1=0.1011,  p2=0.0057, and  p3=0.00015d1-0.00089979. By the Routh-Hurwitz criterion, all real eigenvalues and all real parts of complex conjugate eigenvalues of (3.4) are negative if and only if the following conditions hold:
(3.5)p1>0,p3>0,p1p2-p3>0.
That is, 5.9986<d1<9.8401. Therefore, when 5.9986<d1<9.8401, the equilibrium point O(0,0,0,0) is stable; when 0<d1≤5.9986 or d1≥9.8401, the equilibrium point O(0,0,0,0) is unstable.We fix parameters as above, and then obtain the equilibrium pointS1(1.3748,0.7691,1.6915,0.0412), S2(0.8733,0.8609,0.4999,0.0185), and S3(-1.5441,1.1453,3.4423,1.0183). By calculations, we can obtain that the eigenvalues of the Jacobian matrix of system (2.1) at S1 are λ1=-0.1516, λ2=-0.0644, λ3,4=0.0330±0.2482i; the eigenvalues of the Jacobian matrix of system (2.1) at S2 are λ1=0.0691, λ2=-0.0634, λ3,4=-0.0778±0.2044i; the eigenvalues of the Jacobian matrix of system (2.1) at S3 are λ1=-0.4812, λ2=-0.0613, λ3,4=0.1963±0.2874i. Therefore, S1, S2, and S3 are three saddle points.


## 3.2. Dissipation

Consider the following:(3.6)∇V=∂x˙∂x+∂y˙∂y+∂z˙∂z+∂u˙∂u=a1yM-a1+b2-2b2yC-c3-d3=(a1M-2b2C)y+(b2-a1-c3-d3).Ifa1/M=2b2/C and b2-a1-c3-d3<0, then the system (2.1) is a dissipative system.


## 3.3. Theoretical Proof of the Existence of Smale Horseshoes and the Horseshoes Chaos

To study the long-term dynamical behavior of the system (2.1), the system is divided into subsystems. Let u=0; then the first subsystem is obtained:
(3.7)x˙=a1x(yM-1)-a2y+a3z,y˙=-b1x+b2y(1-yC)+b3z(1-zE),z˙=c1x(xN-1)-c2y-c3z.Letz is a known function of the time t; then the second subsystem is obtained:
(3.8)x˙=a1x(yM-1)-a2y+a3z,y˙=-b1x+b2y(1-yC)+b3z(1-zE)-d4u,u˙=d1y+d2z(zK-1)-d3u.
When t=t0, z is a constant number, then the system (3.8) is a three-dimensional nonlinear system.Theorem 3.2.
The subsystem presented in (3.7) exhibits horseshoe chaos.Proof.
Equation (3.7) has four equilibrium points:
(3.9)O(0,0,0,0),S1(x1,y1,z1,u1),S2(x2,y2,z2,u2),S3(x3,y3,z3,u3).
Step  1. S1 is saddle foci, that is, the eigenvalues of the real matrix A=Df(S1) and the Jacobin derivative of f at S1 are the forms: λ1=r,  λ2,3=σ±iω,  r<0,  σ>0,  |r|>σ, where r,  σ,  ω are real. Here the equilibrium point S1 is discussed with the Jacobin matrix
(3.10)J1=(a1yM-a1a1xM-a2a3-b1b2-2b2yCb3-2b3zE2c1xN-c1-c2-c3).
The characteristic polynomial is obtained as(3.11)det(λI-J1)=λ3+q1λ2+q2λ+q3=0.
Supposeλ=η-q1/3, (3.11) produces η3+qη+p=0, where p=q3-q1q2/3+2q13/27, q=q2-q12/3. Then turn to the Cardano formula, ensure the roots of  (3.11) meet the requirement of  Šilnikov theorem. It yields Δ=(p/2)2+(q/3)2<min(-2q1/3,q1/3).
When parameters are fixed as above atS1(-1.519,1.151,3.460), the following can be obtained: λ1=-0.4786,  λ2,3=0.1943±0.2848i, which satisfy r,  σ±iω,  r<0,  σ>0,  |r|>σ.
Step  2. There exists a homoclinic orbit τ1 at S1.
For the discussion of the homoclinic orbit ofS1, it can be supposed as series form like
(3.12)x(t)=l0+∑k=1+∞lkekαt,y(t)=m0+∑k=1+∞mkekαt,z(t)=n0+∑k=1+∞nkekαt,
where lk,mk,nk(k≥1) are undetermined coefficients, α is attenuation index, when t→+∞, (x(t),y(t),z(t))→(x1,y1,z1).
Next, substitute (3.12) into (2.1), and match constants of items coefficient
(3.13)a1l0m0M-a1l0-a2m0+a3n0=0,-b1l0+b2m0-b2m02C+b3n0-b3n02E=0,c1l02N-c1l0-c2m0-c3n0=0.
The following equation can be obtained:(l0,m0,n0)=(x1,y1,z1). Comparing coefficients of ekαt of the same power terms, the following is obtained:
(3.14)(kαI-J1)(lkmknk)=(φk(1)(ai,bi,ci,α,ξ)φk(2)(ai,bi,ci,α,ξ)φk(3)(ai,bi,ci,α,ξ)).
If(l1,m1,n1)=(0,0,0), then (lk,mk,nk)=(0,0,0), (k>1); therefore, (l1,m1,n1)≠(0,0,0). J1 has the only negative eigenvalues, and then the only α(α<0) satisfies det(αI-J1)=0. Note that det(kαI-J1)≠0 (α is the only negative real eigenvalues of J1, so kα is not the eigenvalues of J1); therefore, lk,mk,nk can be identified uniquely and x(t),  y(t),  z(t) for t>0. For the opposite time symmetric track x(t),  y(t),  z(t), the linear transform can be adopted τ=-t,  t>0. The proof is similar to the procedure as t>0.
According to the above analysis, homoclinic orbitτ1 based on equilibrium point S1 is formally obtained as
(3.15)x(t)={∑k=1+∞lk(ai,bi,ci,α,ξ)ekαt,t>0,∑k=1+∞lk(ai,bi,ci,-α,ξ)e-kαt,t<0,y(t)={∑k=1+∞mk(ai,bi,ci,α,ξ)ekαt,t>0,∑k=1+∞mk(ai,bi,ci,-α,ξ)e-kαt,t<0,z(t)={∑k=1+∞nk(ai,bi,ci,α,ξ)ekαt,t>0,∑k=1+∞nk(ai,bi,ci,-α,ξ)e-kαt,t<0.lk,mk,nk,(k>1) are determined by (3.14), while α is determined by det(αI-J1)=0, and ξ is determined by ∑k=1+∞lk(ai,bi,ciα,ξ)=∑k=1+∞lk(ai,bi,ci,-α,ξ).
Step  3. S3 is saddle foci, and there exists a homoclinic orbit τ2 based at S3. The proof is similar to Steps 1 and 2. So, by the Šilnikov theorem, the horseshoe chaos may expect to occur in a reasonable regime.Theorem 3.3.
The subsystem presented in (3.8) exhibits horseshoe chaos.Proof.
The proof is similar to the proof of Theorem3.2.Remark 3.4.
Theorems3.2 and 3.3 show that the energy-saving and emission-reduction system exhibits horseshoe chaos. This is not to say that the system is chaotic all the time, which means that the system is chaotic under appropriate conditions, and stable under other appropriate conditions.


## 3.4. Numerical Results

We chose a set of parameters as follows:a1=0.09, a2=0.003, a3=0.012, b1=0.0412, b2=0.08, b3=0.8, c1=0.035, c2=0.0062, c3=0.08, c4=0.02, d1=0.01, d2=0.02, d3=0.06, d4=0.03, M=0.9, C=1.6, E=2.8, N=0.35, K=2, and L=2. Let initial condition be (0.015, 0.785, 1.83, 0.01), and the corresponding Lyapunov exponents are L1=0.0272>0, L2=-0.0011<0, L3=-0.0839<0, and L4=-0.0772<0. Therefore, the Lyapunov dimension of this system is
(3.16)Dl=j+1|Lj+1|∑i=1jLi=3+L1+L2+L3|L4|=2.2522,
which means that the Lyapunov dimension is fractional under the same condition. The system has a chaotic attractor, as shown in Figures 1(a)–1(c), the time series of x(t),y(t),z(t),u(t) as shown in Figure 1(d). Let d1=0.094, fixed other parameters and initial condition, and then produce a limit cycle as shown in Figure 2. By calculations, the Lyapunov exponent spectrum with respect to parameter d1 is shown in Figure 3. According to Figure 3, the system (2.1) has very rich dynamical behaviors, which are summarized as follows. When d1∈(0.0935,0.095), the system (2.1) is chaotic with a positive Lyapunov exponent (e.g., with d1=0.01, the phase portrait is shown in Figure 1); while for d1∈(0.0935,0.095), the maximum Lyapunov exponent equals zero, implying that the system has a periodic orbit (Figure 2 shows the periodic orbit when d1=0.094).A four-dimensional energy-saving and emission-reduction chaotic attractor.
(a)
(b)
(c)
(d)A limit cycle.
(a)
(b)
(c)
(d)Figure 3
Lyapunov exponent spectrum.


## 4. Linear Feedback Control

Because an energy resource system in the chaotic state is very sensitive to its initial condition and chaos often causes irregular behavior, chaos is undesirable. In this section, linear feedback methods [21] are applied to control chaos of the energy resource system (2.1). Firstly, we prove this chaos can be controlled to equilibrium point O(0,0,0,0).We guide the chaotic trajectory(x(t),y(t),z(t),u(t)) to equilibrium point O(0,0,0,0). Let the system (2.1) be controlled by a linear feedback control of the form:
(4.1)x˙=a1x(yM-1)-a2y+a3z-F11x,y˙=-b1x+b2y(1-yC)+b3z(1-zE)-d4u-F22y,z˙=c1x(xN-1)-c2y-c3z+c4u(uL-1)-F33z,u˙=d1y+d2z(zK-1)-d3u-F44u,
where F11, F22, F33, F44 are the positive feedback gains, which are needed to be chosen such that the trajectory of the system (2.1) is stabilized to equilibrium point O(0,0,0,0).The Jacobian matrix of the system (4.1) is
(4.2)J0=(-a1-F11-a2a30-b1b2-F22b3-d4-c1-c2-c3-F33-c40d1-d2-d3-F44),
where a1=0.09, a2=0.003, a3=0.012, b1=0.0412, b2=0.08, b3=0.8, c1=0.035, c2=0.0062, c3=0.08, c4=0.02, d1=0.01, d2=0.02, d3=0.06, and d4=0.03. The Jacobian matrix (4.2) is
(4.3)J0=(-0.09-F11-0.030.0120-0.4120.08-F220.8-0.03-0.035-0.0062-0.08-F33-0.0200.01-0.02-0.06-F44).LetF11=F33=0; it has the characteristic equation
(4.4)f(λ)=(λ+0.0489)[λ3+(F22+F44+0.0511)λ2+(F22F44+0.1311F22-0.0089F44-0.0073)λ+0.0711F22F44+0.0044F22-0.0014F44+0.00049796λ3]=0.According to Routh-Hurwitz criteria, if(4.5)F22+F44+0.0511>0,0.0711F22F44+0.0044F22-0.0014F44+0.00049796>0,(F22+F44+0.0511)(F22F44+0.1311F22-0.0089F44-0.0073)>0.0711F22F44+0.0044F22-0.0014F44+0.00049796,
then we know that the Jacobian matrix J0 has four negative real part eigenvalues. When F22 and F44 satisfy (4.5), the controlled system (4.1) is asymptotically stable at the equilibrium O(0,0,0,0).Numerical experiments are carried out to integrate the controlled system (3.7) by the MATLAB. The parameters are chosen as a1=0.09, a2=0.003,  a3=0.012, b1=0.0412, b2=0.08, b3=0.8, c1=0.035, c2=0.006, c3=0.08, c4=0.02, d1=0.01, d2=0.02, d3=0.06, d4=0.03, M=0.9, C=1.6,  E=2.8, N=0.35, K=2, and L=2 to ensure the existence of chaos in the absence of control. Let initial states be (0.015, 0.785, 1.83, 0.01); when F11=F33=0, F22=F44=0.07, the equilibrium point O(0,0,0,0) of system (2.1) is stabilized as shown Figure 4.The stable equilibrium pointO(0,0,0,0) of the controlled system (4.1).
(a)
(b)
(c)
(d)Similarly, we can prove that another three equilibrium pointsS1, S2, and S3 of the system (2.1) can be stabilized. We fix parameters as above, let initial states be (1.5, 0.785, 1.5, 0.5). When F11=F33=0,F22=F44=0, the equilibrium point S1(1.3748,0.7691,1.6915,0.0412) of system (2.1) is stabilized as shown Figure 5.The stable equilibrium pointS1 of the controlled system (4.1).
(a)
(b)
(c)
(d)


## 5. Conclusion

We have established a four-dimensional nonlinear dynamics model for the energy-saving and emission-reduction system and have analyzed the dynamics behavior of the system. When some parameters are adjusted, the dynamic behavior of energy-saving and emission-reduction, carbon emissions, economic growth, and new energy development displays some regulated phenomena. By observing these phenomena, we can figure out the affecting factors for energy intensity and grasp the statistical results which meet the real situation. This four-dimensional energy-saving and emission-reduction system will be more satisfactory for actual energy saving and emission reduction and instructive for the energy saving and emission reduction of China. The research results provide a key to energy saving and emission reduction, that is, to develop energy-saving and emission reduction as soon as possible with proper strategies rather than simply increasing investment. The theoretical proof and the empirical study ensure the necessity and significance to carry out comprehensive energy saving and emission reduction.



---

*Source: 842716-2012-12-11.xml*