# Research on Third-Party Collecting Game Model with Competition in Closed-Loop Supply Chain Based on Complex Systems Theory

Authors: Junhai Ma; Yuehong Guo
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This paper studied system dynamics characteristics of closed-loop supply chain using repeated game theory and complex system theory. It established decentralized decision-making game model and centralized decision-making game model and then established and analyzed the corresponding continuity system. Drew the region local stability of Nash equilibrium and Stackelberg equilibrium, and a series of chaotic system characteristics, have an detail analysis of the Lyapunov index which is under the condition of different parameter combination. According to the limited rational expectations theory, it established repeated game model based on collection price and marginal profits. Further, this paper analyzed the influence of the parameters by numerical simulations and concluded three conclusions. First, when the collection price is to a critical value, the system will be into chaos state. Second, when the sale price of remanufacturing products is more than a critical value, the system will be in chaos state. Last, the initial value of  the collection price  is sensitive, small changes may cause  fluctuations of market price. These conclusions guide enterprises in making the best decisions in each phase to achieve maximize profits.


---

## Body

## 1. Introduction

The product life cycle becomes shorter with the rapid development of market economy. A lot of products have been washed out before life end. Such condition not only creates a tremendous waste of resources but also brings great harm to people. Therefore, the “resources-products-waste products-remanufacturing product” closed-loop type of economic growth mode appears. It realizes the economic development and utilization of resources and environmental protection in coordination strategy of sustainable development goals. Enterprise began to take a positive attitude to collect products from customers. The research on collecting control problem becomes inevitable.This paper draws on and contributes to several streams of literature, each of whom we review below. A growing body of literature in operations management addresses reverse logistics management issues for remanufacturable products. References [1, 2] defined as the supply chain to reverse supply chain which formed a complete closed-loop system (closed-loop supply chain, hereinafter referred to as CLSC.) Reference [3] studied the pricing methods under three collecting channels which manufacturer, retailer, and third-party collector collected products, respectively. Reference [4] studied manufacturers should be responsible for the management of that collecting based on Xerox corporation and Kodak corporation. Reference [5] studied three kinds of collecting channels of manufacturers collecting, retailer collecting, and the third-party collecting in closed-loop supply chain. They found that the colleting distance was closer from consumers and the collecting efforts were more effective. References [6–8] described reverse logistics management mode of the third party. In this collecting mode, manufacturers entrust a third party to perform extended producer responsibility and manage used products.References [9, 10] studied pricing analysis of reverse supply chain using the game theory method under retailers collecting mode and concluded the optimal strategy of single stage. Reference [11] established the demand model of new products and remanufacturing products from social environmental protection consciousness and consumers’ utility function to new products and remanufacturing products and drew the conclusions. The first is the optimal pricing strategy in different collecting channels in decentralized decision-making. The second is the influence of social environmental protection consciousness. Reference [12] research shows that the channel of manufacturers collecting directly is more favorable to consumers and social environment. In the autoindustry, independent third parties are handling used-product collection activities for the original equipment manufacturers (OEMs). The “big three” automanufacturers in the United States started to invest in joint research and remanufacturing partnerships with dismantling centers to benefit from their scale economies and experience. Third parties such as GENCO Distribution System are also preferred by some consumer goods manufacturers for their experience in used-product collection. In references [13–15] Analyzed the complexity of model for a class of delay complex dynamics, and some valuable conclusions are obtained.The existing literatures studied major in a single phase of the mathematical model using game theory under the assumption that the closed-loop supply chain contains a manufacturer and a retailer and concludes equilibrium in perfectly rational state. However, market uncertainty caused enterprise not to be able to make decisions in perfectly rational. At present, many manufacturers devote themselves to the core competitiveness; they contract the collection of used products to a third party. So this paper sets up multistage game model of closed-loop supply chain which consists of a manufacturer and two competition collectors business based on limited rationality with China’s remanufacturing development and studies its complex dynamic characteristics.


## 2. Model Assumptions and Notation

### 2.1. Assumption

First, this paper considers only remanufacturing products market. That means new products and remanufacturing product do not form a competitive market. It is more in line with the actual conditions of China. The manufacturer manufactures new products and remanufacturing products.Second, the closed-loop supply chain includes a manufacturer and two competitive collecting corporations, the collecting corporations collect waste products from consumers and return the manufacturer. The manufacturer transfers payments to the collectors, and she sales directly to consumers. The manufacturers and two collectors are independent decision makers, and their strategic space is to choose the best collecting price. Their goal is to maximize returns in discrete time period ast
=
0,1
,
2
,
….Third, the number of the collection is increasing function of collecting price. The collecting capability and manufacturing capability are unlimited. In order to simplify the problem and emphasize main parameters influence of the system, all the collecting products can be manufactured, as shown in the MCTM collecting mode (Figure1).Figure 1
MCTM collecting mode.


### 2.2. Notation

p
0 denotes the retail price of the product, a constant, and c
r the unit cost of remanufacturing a returned product into a new one, a constant. p
m
(
t
) is the collection price of the manufacturer transfer payment to collectors in the t period, a decision variable of the manufacturer. p
r
i
(
t
) is the collection price of collectors in the t period. r
i
(
t
) is a decision variable of the collectors; that is, p
r
i
=
p
m
(
1
-
r
i
), i
=
1,2. A
i is the collecting cost of collectors including logistics cost and so forth, a constant. G
i is the collecting numbers for the waste product in the market as a function of collecting price, G
i
=
k
i
+
β
i
p
r
i
(
t
)
-
δ
i
p
r
j
(
t
), i
=
1,2, j
=
2,1, and subscript i, j will take values two collectors, where k is the waste products number of return voluntarily when the collectors’ collecting price is zero; it denotes consumers environmental awareness; β denotes consumers sensitive degree to collecting price; δ is competition coefficient in the two collectors, which satisfies β
≥
δ. The collecting function of the two collectors is, respectively, as follows:
(1)
G
1
=
k
1
+
β
1
p
r
1
(
t
)
-
δ
1
p
r
2
(
t
)
,
G
2
=
k
2
+
β
2
p
r
2
(
t
)
-
δ
2
p
r
1
(
t
)
.Π
m is the profits function of the manufacturer; its marginal profits are concluded by the first-order conditions, which are ∂
Π
m
/
∂
p
m
(
t
). Π
m can be stated as
(2)
∏
m
=
(
p
0
-
p
m
-
c
r
)
×
(
k
1
+
β
1
p
r
1
(
t
)
-
δ
1
p
r
2
(
t
)
h
h
h
h
+
k
2
+
β
2
p
r
2
(
t
)
-
δ
2
p
r
1
(
t
)
)
.Π
r
i is the profits function of the two collectors; marginal profits are concluded by the first-order conditions, which are ∂
Π
r
i
/
∂
r
i
(
t
). Π
r
i can be stated as
(3)
∏
r
1
=
(
p
m
-
p
r
1
(
t
)
-
A
1
)
(
k
1
+
β
1
p
r
1
(
t
)
-
δ
1
p
r
2
(
t
)
)
,
∏
r
2
=
(
p
m
-
p
r
2
(
t
)
-
A
2
)
(
k
2
+
β
2
p
r
2
(
t
)
-
δ
2
p
r
1
(
t
)
)
.Π is the profits function of the closed-loop supply chain, which can be stated as
(4)
∏
=
∏
m
+
∏
r
1
+
∏
r
2
.(
p
m
,
r
1
,
r
2
) is a pricing strategy, which can be obtained by (2) and (3). Furthermore they are satisfied with ∂
2
Π
m
/
∂
p
m
2
=
-
2
(
1
-
r
2
)
(
β
2
-
δ
1
)
<
0, ∂
2
Π
r
1
/
∂
r
1
2
=
-
2
β
1
p
m
2
(
t
)
<
0, and ∂
2
Π
r
2
/
∂
r
2
2
=
-
2
β
2
p
m
2
(
t
)
<
0. It can be proved that Π
i is concave function for p
i.


## 2.1. Assumption

First, this paper considers only remanufacturing products market. That means new products and remanufacturing product do not form a competitive market. It is more in line with the actual conditions of China. The manufacturer manufactures new products and remanufacturing products.Second, the closed-loop supply chain includes a manufacturer and two competitive collecting corporations, the collecting corporations collect waste products from consumers and return the manufacturer. The manufacturer transfers payments to the collectors, and she sales directly to consumers. The manufacturers and two collectors are independent decision makers, and their strategic space is to choose the best collecting price. Their goal is to maximize returns in discrete time period ast
=
0,1
,
2
,
….Third, the number of the collection is increasing function of collecting price. The collecting capability and manufacturing capability are unlimited. In order to simplify the problem and emphasize main parameters influence of the system, all the collecting products can be manufactured, as shown in the MCTM collecting mode (Figure1).Figure 1
MCTM collecting mode.


## 2.2. Notation

p
0 denotes the retail price of the product, a constant, and c
r the unit cost of remanufacturing a returned product into a new one, a constant. p
m
(
t
) is the collection price of the manufacturer transfer payment to collectors in the t period, a decision variable of the manufacturer. p
r
i
(
t
) is the collection price of collectors in the t period. r
i
(
t
) is a decision variable of the collectors; that is, p
r
i
=
p
m
(
1
-
r
i
), i
=
1,2. A
i is the collecting cost of collectors including logistics cost and so forth, a constant. G
i is the collecting numbers for the waste product in the market as a function of collecting price, G
i
=
k
i
+
β
i
p
r
i
(
t
)
-
δ
i
p
r
j
(
t
), i
=
1,2, j
=
2,1, and subscript i, j will take values two collectors, where k is the waste products number of return voluntarily when the collectors’ collecting price is zero; it denotes consumers environmental awareness; β denotes consumers sensitive degree to collecting price; δ is competition coefficient in the two collectors, which satisfies β
≥
δ. The collecting function of the two collectors is, respectively, as follows:
(1)
G
1
=
k
1
+
β
1
p
r
1
(
t
)
-
δ
1
p
r
2
(
t
)
,
G
2
=
k
2
+
β
2
p
r
2
(
t
)
-
δ
2
p
r
1
(
t
)
.Π
m is the profits function of the manufacturer; its marginal profits are concluded by the first-order conditions, which are ∂
Π
m
/
∂
p
m
(
t
). Π
m can be stated as
(2)
∏
m
=
(
p
0
-
p
m
-
c
r
)
×
(
k
1
+
β
1
p
r
1
(
t
)
-
δ
1
p
r
2
(
t
)
h
h
h
h
+
k
2
+
β
2
p
r
2
(
t
)
-
δ
2
p
r
1
(
t
)
)
.Π
r
i is the profits function of the two collectors; marginal profits are concluded by the first-order conditions, which are ∂
Π
r
i
/
∂
r
i
(
t
). Π
r
i can be stated as
(3)
∏
r
1
=
(
p
m
-
p
r
1
(
t
)
-
A
1
)
(
k
1
+
β
1
p
r
1
(
t
)
-
δ
1
p
r
2
(
t
)
)
,
∏
r
2
=
(
p
m
-
p
r
2
(
t
)
-
A
2
)
(
k
2
+
β
2
p
r
2
(
t
)
-
δ
2
p
r
1
(
t
)
)
.Π is the profits function of the closed-loop supply chain, which can be stated as
(4)
∏
=
∏
m
+
∏
r
1
+
∏
r
2
.(
p
m
,
r
1
,
r
2
) is a pricing strategy, which can be obtained by (2) and (3). Furthermore they are satisfied with ∂
2
Π
m
/
∂
p
m
2
=
-
2
(
1
-
r
2
)
(
β
2
-
δ
1
)
<
0, ∂
2
Π
r
1
/
∂
r
1
2
=
-
2
β
1
p
m
2
(
t
)
<
0, and ∂
2
Π
r
2
/
∂
r
2
2
=
-
2
β
2
p
m
2
(
t
)
<
0. It can be proved that Π
i is concave function for p
i.


## 3. Decentralized Control Decision-Making Model

### 3.1. Nash Equilibrium

#### 3.1.1. Model and Analysis

(
1) Model. The manufacturer and collectors are equal in closed-loop supply chain; that is, the two sides make decision at the same time: manufacturer’s decision is to choose the collection price p
i to maximize returns, and collectors’ decision is to choose p
r
i to maximize their returns.Nash equilibrium can be concluded by the first-order conditions of three reaction functions, which meet∂
Π
m
/
∂
p
m
=
0, ∂
Π
r
i
/
∂
r
i
=
0 simultaneously. The Nash equilibrium is the optimal decision of a game in all kinds of may guess, thus make it to gain maximum benefit. The Nash equilibrium is
(5)
(
p
m
*
,
r
1
*
,
r
2
*
)
=
(
-
B
+
B
2
-
4
A
C
2
A
,
t
1
+
t
2
p
m
*
,
t
3
+
t
4
p
m
*
)
,
where A, B, and C are
(6)
A
=
2
(
β
1
+
β
2
-
δ
1
-
δ
2
)
(
p
0
-
c
r
)
h
h
h
-
2
t
1
(
β
1
-
δ
2
)
-
2
t
3
(
β
2
-
δ
1
)
,
B
=
k
1
+
k
2
-
(
β
1
+
β
2
-
δ
1
-
δ
2
)
(
p
0
-
c
r
)
h
h
h
+
(
p
0
-
c
r
)
(
t
1
(
β
1
-
δ
2
)
+
t
3
(
β
2
-
δ
1
)
)
h
h
h
-
2
t
2
(
β
1
-
δ
2
)
-
2
t
4
(
β
2
-
δ
1
)
,
C
=
(
p
0
-
c
r
)
t
2
(
β
1
-
δ
2
)
+
t
4
(
β
2
-
δ
1
)
.In formula (6), t
1, t
2, t
3, t
4 can be expressed as
(7)
t
1
=
2
β
1
β
2
-
2
β
2
δ
1
+
δ
1
β
2
-
δ
1
δ
2
4
β
1
β
2
-
δ
1
δ
2
,
t
2
=
2
β
2
k
1
+
2
β
2
β
1
A
1
+
δ
1
k
2
-
δ
1
β
2
A
2
4
β
1
β
2
-
δ
1
δ
2
,
t
3
=
2
β
1
β
2
-
2
β
1
δ
2
+
δ
2
β
1
-
δ
1
δ
2
4
β
1
β
2
-
δ
1
δ
2
,
t
4
=
2
β
1
k
2
+
2
β
2
β
1
A
2
+
δ
2
k
1
-
δ
2
β
1
A
1
4
β
1
β
2
-
δ
1
δ
2
.In reality, the game between node enterprises in closed-loop supply chain is continuous, enterprises’ decision-making is a long-term repeated process, and its action has long-term memory. And each node enterprise does not completely control the market information and also cannot fully expect future market changes, so based on limited rational expectations decision we adjust process with marginal gains. They can make the next-period price decision on the basis of the local estimate to his marginal profits in current period. Their price adjustment processes are(8)
p
m
(
t
+
1
)
=
p
m
(
t
)
+
α
1
p
m
(
t
)
×
(
(
(
1
-
r
1
(
t
)
)
(
β
1
-
δ
2
)
+
(
1
-
r
2
(
t
)
)
(
β
2
-
δ
1
)
)
×
(
p
0
-
2
p
m
(
t
)
-
c
r
)
-
k
1
-
k
2
(
(
1
-
r
1
(
t
)
)
(
β
1
-
δ
2
)
+
(
1
-
r
2
(
t
)
)
(
β
2
-
δ
1
)
)
)
,
r
1
(
t
+
1
)
=
r
1
(
t
)
+
α
2
r
1
(
t
)
×
(
-
2
β
1
r
1
(
t
)
p
m
2
(
t
)
+
p
m
(
t
)
h
h
h
×
(
k
1
+
(
β
1
-
δ
1
+
δ
1
r
2
(
t
)
)
h
h
h
×
p
m
(
t
)
+
β
1
A
1
)
(
-
2
β
1
r
1
(
t
)
p
m
2
(
t
)
+
p
m
(
t
)
)
,
r
2
(
t
+
1
)
=
r
2
(
t
)
+
α
3
r
2
(
t
)
×
(
-
2
β
2
r
2
(
t
)
p
m
2
(
t
)
+
p
m
(
t
)
h
h
×
(
k
2
+
(
β
2
-
δ
2
+
δ
2
r
1
(
t
)
)
p
m
(
t
)
+
β
2
A
2
)
(
-
2
β
2
r
2
(
t
)
p
m
2
(
t
)
+
p
m
(
t
)
)
,
where p
1
(
t
+
1
),  r
1
(
t
+
1
), r
2
(
t
+
1
), and p
1
(
t
), r
1
(
t
), r
2
(
t
) are used to present the price in the (
t
+
1
)th period and in the (
t
)th period. Where α
i are positive parameters, i
=
1,2
,
3, which denotes adjustment speed, respectively, for the manufacturer and two collectors.From the system (8) we can be conclude that the optimal collecting price of the manufacturer is related to the collection price adjustment coefficient, the two collectors’ collecting price, sale price, and consumers’ environmental protection awareness. Similarly, the optimal collecting price of collectors is related to collecting price adjustment parameters, collecting price of the manufacturer, collecting function, consumers’s environmental protection awareness, and collecting costs.(2) Model Analysis. In closed-loop supply chain, any enterprise decision-making is all according to the maximum profits, the equilibrium less than zero is not practical significance, such as manufacturers to make power mainly because of benefits, collecting business is also benefit with collecting. So only research system (8) is the equilibrium; system (8) of the eight fixed points, respectively, is E
0
(
0,0
,
0
),  E
1
(
p
m
1
,
0,0
),  E
2
(
0
,
r
1
2
,
0
),  E
3
(
0
,
0
,
r
2
3
),  E
4
(
p
m
4
,
r
1
4
,
0
), E
5
(
p
m
5
,
0
,
r
2
5
), E
6
(
0
,
r
1
6
,
r
2
6
), E
7
(
p
m
*
,
r
1
*
,
r
2
*
). The analysis shows that by E
0, E
1, E
2, E
3, E
4, E
5 and E
6, is no real significance. Therefore, we discuss Nash equilibrium which is E
7
(
p
m
*
,
r
1
*
,
r
2
*
); E
7 is the reaction function intersection of the manufacturer and collectors, which means marginal profits of both sides are zero, but this does not mean that the result of the game will tend to equilibrium. Instead, a party rational behavior change may cause game process to occur very complex phenomenon. E
7 has local stability, this stability region by α
1, α
2, α
3 decided. The stability of E
7 will be studied. The first we calculate Jacobi matrix of system (8), means
(9)
J
=
(
J
11
J
12
J
13
J
21
J
22
J
23
J
31
J
32
J
33
)
.We putE
7
(
p
m
*
,
r
1
*
,
r
2
*
) into formula (9), according to the Jury stability criterion, the sufficient and necessary conditions of E
7 asymptotic stability are that the zeros of characteristic polynomial F
(
λ
)
=
λ
3
+
A
λ
2
+
B
λ
+
C are in a circular area. Moreover, the characteristic polynomial will satisfy the conditions as follows:
(10)
F
(
1
)
=
1
+
A
+
B
+
C
>
0
,
(
-
1
)
3
F
(
-
1
)
=
1
-
A
+
B
-
C
>
0
,
1
>
|
C
|
,
1
-
C
2
>
|
B
-
A
C
|
.At present, only a few simple dynamic systems are analyzed with analytical method, but complex dynamic system mainly uses the numerical analysis method. This paper processes numerical simulation on system (8) by Matlab to expresses the dynamic characteristics. Because remanufacturing is still initial stage in china, consumers environmental protection awareness is lower, and manufacturing cost (including the fixed and variable cost) is higher; for autoparts collecting, the parameters of system (8) can be defined as p
0
=
70, c
r
=
30, A
1
=
1, A
2
=
0.8, k
1
=
1, k
2
=
1, β
1
=
0.6, β
2
=
0.55, δ
1
=
0.2, and δ
2
=
0.2 to study the local stability of equilibrium point. Get the parameter value into E
7
(
p
m
*
,
r
1
*
,
r
2
*
) and formula (10) to conclude the value for (17.385, 0.491, 0.489).Change the discrete system (8) into a continuous system and the points of Jacobi matrix for
(11)
J
=
[
-
13.317
α
1
-
17.755
α
1
-
15.562
α
1
-
0.775
α
2
-
177.888
α
2
29.68
α
2
-
0.821
α
3
29.599
α
3
-
164.65
α
3
]
.Consider the continuous-time nonlinear dynamical system(12)
x
˙
=
A
x
+
N
(
x
)
,
A
=
(
a
i
j
)
n
×
n
,
x
∈
R
n
.Let the functionN
(
x
) be written as
(13)
N
(
x
)
=
1
2
B
(
x
,
x
)
+
1
6
C
(
x
,
x
,
x
)
+
O
(
∥
x
∥
4
)
,
where B
(
x
,
y
) and C
(
x
,
y
,
z
) are bilinear and trilinear functions, respectively. In coordinates, we have
(14)
B
i
(
x
,
y
)
=
∑
j
,
k
=
1
n
∂
2
N
i
(
ξ
)
∂
ξ
j
∂
ξ
k
|
ξ
=
0
x
j
y
k
,
C
i
(
x
,
y
,
z
)
=
∑
j
,
k
,
l
=
1
n
∂
3
N
i
(
ξ
)
∂
ξ
j
∂
ξ
k
∂
ξ
l
|
ξ
=
0
x
j
y
k
z
l
.Supposing thatA has a pair of complex eigenvalues on the imaginary axis, they are λ
1,2
=
±
i
ω
(
ω
>
0
), and these eigenvalues are the only eigenvalues with Re
λ
=
0. Let q
∈
C
n be a complex eigenvector to λ
1
=
i
ω:
(15)
A
q
=
i
ω
q
,
A
q
¯
=
-
i
ω
q
¯
.And the adjoint eigen vectorp
∈
C
n admits the properties
(16)
A
T
p
=
-
i
ω
p
,
A
T
p
¯
=
i
ω
p
¯and satisfies the normalization〈
q
,
p
〉
=
1.The first Lyapunov coefficient at the origin is defined by(17)
l
1
=
1
2
ω
Re
[
〈
p
,
C
(
q
,
q
,
q
¯
)
〉
-
2
〈
p
,
B
(
q
,
A
-
1
B
(
q
,
q
¯
)
)
〉
hhhhhhhh
+
〈
p
,
B
(
q
¯
,
(
2
i
ω
E
-
A
)
-
1
B
(
q
,
q
)
)
〉
]
.Next, we calculatei
w
q, J
q
-
=
-
i
w
q
-, J
T
p
=
-
i
w
p, and  J
T
p
-
=
i
w
p
-; If q
=
(
q
1
,
q
2
,
q
3
)
T
≠
0, calculate q
,
q
-
,
p
,
p
-:
(18)
[
-
13.317
α
1
-
i
w
-
17.755
α
1
-
15.562
α
1
-
0.775
α
2
-
177.888
α
2
-
i
w
29.68
α
2
-
0.821
α
3
29.599
α
3
-
164.65
α
3
-
i
w
]
×
[
q
1
q
2
q
3
]
=
[
0
0
0
]
having nonzero solutions.It means(19)
|
-
13.317
α
1
-
i
w
-
17.755
α
1
-
15.562
α
1
-
0.775
α
2
-
177.888
α
2
-
i
w
29.68
α
2
-
0.821
α
3
29.599
α
3
-
164.65
α
3
-
i
w
|
=
0
.It can be got from calculating(20)
f
1
(
w
)
=
-
375716.354
α
1
α
2
α
3
+
117.77755
α
2
2
α
3
+
13.317
α
1
w
2
+
177.888
α
2
w
2
+
164.65
α
3
w
2
=
0
,
f
2
(
w
)
=
-
2368.9345
α
1
α
2
w
+
0.6
α
2
2
w
-
28410.761
α
2
α
3
w
-
2179.87
α
1
α
3
w
+
w
3
=
0
.(1) When{
w
=
α
1
f
1
(
w
)
=
0
f
2
(
w
)
=
0, matrix J has pure imaginary eigenvalues ±
i
α
1. Calculate q, q
-, p, p
-.The result is as follows.(21)
L
2
(
α
1
,
α
2
,
α
3
)
=
1
2
ω
[
〈
p
,
C
(
q
,
q
,
q
-
)
〉
-
2
〈
p
,
B
(
q
,
A
-
1
B
(
q
,
q
-
)
)
〉
h
h
h
h
h
h
+
〈
p
,
B
(
q
-
,
(
2
i
ω
E
-
A
)
-
1
B
(
q
,
q
)
)
〉
]
=
(
408.75
-
29.6
i
)
α
1
α
2
(
0.82
i
α
1
-
168.9853
α
2
)
α
2
-
2
[
(
408.75
-
29.6
i
)
α
1
α
2
(
0.82
i
α
1
-
168.9853
α
2
)
α
2
-
2179.87
α
1
375716.354
α
1
α
2
-
117.7775
α
2
2
+
(
+
α
1
α
2
(
1.2
i
α
2
+
(
375716.354
+
56821.5
i
)
α
3
)
}
)
-
1
(
61899.185
-
4510.5
i
)
α
2
×
(
5846897.898
i
α
1
2
-
(
1.24
×
1
0
9
+
1832.854
i
)
α
1
α
2
+
388107.8
α
2
2
)
-
1
(
(
61899.185
-
4510.5
i
)
α
2
)
(
408.75
-
29.6
i
)
α
1
α
2
(
0.82
i
α
1
-
168.9853
α
2
)
α
2
]
+
(
+
α
1
α
2
(
1.2
i
α
2
+
(
375716.354
+
56821.5
i
)
α
3
)
}
)
-
1
{
α
1
[
(
-
4
+
26.634
i
)
α
1
h
h
h
h
h
h
h
h
+
(
2179.87
+
329.3
i
)
α
3
]
}
×
(
{
(
-
53
-
8
i
)
α
1
3
+
α
1
2
[
(
-
711.55
+
4737.8
i
)
α
2
-
(
658.6
-
4359.7
i
)
α
3
]
-
117.8
α
2
2
α
3
+
α
1
α
2
(
1.2
i
α
2
+
(
375716.354
+
56821.5
i
)
α
3
)
(
-
53
-
8
i
)
α
1
3
}
)
-
1
)
+
{
(
-
(
407.31
+
29.7
i
)
α
2
)
×
(
-
1.55
i
α
1
α
2
-
151.97
α
2
α
3
)
}
×
(
{
(
15.6
i
α
1
+
3295.3
α
2
)
(
-
53
-
8
i
)
α
1
3
×
{
(
-
53.3
-
8
i
)
α
1
3
(
-
53
-
8
i
)
α
1
3
+
α
1
2
[
(
-
711.552
+
4737.9
i
)
α
2
-
(
658.6
-
4359.74
i
)
α
3
]
-
117.77755
α
2
2
α
3
+
α
1
α
2
×
[
1.2
i
α
2
+
(
375716.354
+
56821.52
i
)
α
3
]
(
-
53
-
8
i
)
α
1
3
}
}
)
-
1
=
l
1
(
α
1
,
α
2
,
α
3
)
+
m
1
(
α
1
,
α
2
,
α
3
)
,
where L
1
(
α
1
,
α
2
,
α
3
) is the Lyapunov exponent function about α
1, α
2, and α
3 and l
1
(
α
1
,
α
2
,
α
3
) is the real part. The following is the figure of l
1
(
α
1
,
α
2
,
α
3
).(1) Fix the value ofα
1 and α
2
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
1
=
0.02, α
2
=
0.048, and α
3
∈
(
0,0.05
) (see Figure 3).
(ii)
α
1
=
0.02, α
2
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 2).Figure 2
The change ofl
1
(
α
1
,
α
2
,
α
3
) with α
1
,
α
2, and α
3.Figure 3
The change ofl
1
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.048.We can know from Figure3 that when α
1 =  0.02, α
2 = 0.048, and α
3 is between zero and 0.05, l
1
(
α
3
) is always equal to 1.375.(iii)
α
1
=
0.02, α
2
=
0.02, and α
3
∈
(
0,0.05
) (see Figure 4).
(iv)
α
1
=
0.02, α
2
=
0.004, and α
3
∈
(
0,0.05
) (see Figure 5).Figure 4
The change ofl
1
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.02.Figure 5
The change ofl
1
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.004.From Figures4 and 5 we can observe the same change of l
1
(
α
3
) as in Figures 2 and 3.(2) Fix the value ofα
2 and α
1
∈
(
0,0.05
), α
3
∈
(
0,0.05
). α
2
=
0.004, α
1
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 6). α
2
=
0.004, α
3
=
0.01, and α
1
∈
(
0,0.05
) (see Figure 7). α
2
=
0.004, α
3
=
0.03, and α
1
∈
(
0,0.05
) (see Figure 8). α
2
=
0.004, α
3
=
0.05, and α
1
∈
(
0,0.05
) (see Figure 9).Figure 6
The change ofl
1
(
α
1
,
α
2
,
α
3
) with α
2, α
1, and α
3.Figure 7
The change ofl
1
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.01.Figure 8
The change ofl
1
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.03.Figure 9
The change ofl
1
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.05.(3) Fix the value ofα
3 and α
1
∈
(
0,0.05
), α
2
∈
(
0,0.05
).(i)
α
3
=
0.03, α
1
∈
(
0,0.05
), and α
2
∈
(
0,0.05
) (see Figure 10).
(ii)
α
3
=
0.03, α
1
=
0.01, and α
2
∈
(
0,0.05
) (see Figure 11).
(iii)
α
3
=
0.03, α
1
=
0.02, and α
2
∈
(
0,0.05
) (see Figure 12).
(iv)
α
3
=
0.03, α
1
=
0.03, and α
2
∈
(
0,0.05
) (see Figure 13).Figure 10
The change ofl
1
(
α
1
,
α
2
,
α
3
) with α
3, α
1, and α
2.Figure 11
The change ofl
1
(
α
2
) with α
2 when α
1
=
0.01, α
3
=
0.03.Figure 12
The change ofl
1
(
α
2
) with α
2 when α
1
=
0.02, α
3
=
0.03.Figure 13
The change ofl
1
(
α
2
) with α
2 when α
1
=
0.03, α
3
=
0.03.(2) When{
w
=
α
2
f
1
(
α
2
)
=
0
f
2
(
α
2
)
=
0, matrix J has pure imaginary eigenvalues ±
i
α
2. Calculate q, q
-, p, p
-.The result is as follows:(22)
L
2
(
α
1
,
α
2
,
α
3
)
=
1
2
ω
[
〈
p
,
C
(
q
,
q
,
q
-
)
〉
-
2
〈
p
,
B
(
q
,
A
-
1
B
(
q
,
q
-
)
)
〉
+
〈
p
,
B
(
q
-
,
(
2
i
ω
E
-
A
)
-
1
B
(
q
,
q
)
)
〉
]
=
(
0.000851
-
0.175
i
)
+
37.57
+
0.18
i
375716.354
α
1
-
117.8
α
2
+
(
2.42
+
0.012
i
)
α
1
α
2
+
(
0.01293
-
2.74
i
)
α
2
375716.354
α
1
2
-
117.78
α
1
α
2
+
4359.74
α
1
375716.354
α
1
α
2
-
117.77755
α
2
2
+
(
{
α
2
2
[
(
-
0.014
+
0.000066
i
)
α
2
+
(
0.0065
+
1.37
i
)
α
3
]
+
α
1
α
2
[
(
-
4
+
0.19
i
)
α
2
+
(
18.784
+
329.2
i
)
α
3
]
+
α
1
2
(
26.6
i
α
2
+
2179.87
α
3
)
}
)
×
{
α
1
α
2
[
(
-
53.268
+
4737.87
i
)
α
1
α
2
-
(
711.552
+
9.2
i
)
α
2
2
+
(
375716.354
+
4359.74
i
)
α
1
α
3
-
(
776.38
-
56821.52
i
)
α
2
α
3
]
}
-
1
=
l
2
(
α
1
,
α
2
,
α
3
)
+
m
2
(
α
1
,
α
2
,
α
3
)
,
where L
2
(
α
1
,
α
2
,
α
3
) is the Lyapunov exponent function about α
1, α
2, and α
3 and l
2
(
α
1
,
α
2
,
α
3
) is the real part. The following is the figure of l
2
(
α
1
,
α
2
,
α
3
).(1) Fix the value ofα
1 and α
2
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
1
=
0.02, α
2
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 14).
(ii)
α
1
=
0.02, α
2
=
0.048, and α
3
∈
(
0,0.05
) (see Figure 15).
(iii)
α
1
=
0.02, α
2
=
0.02, and α
3
∈
(
0,0.05
) (see Figure 16).
(iv)
α
1
=
0.02, α
2
=
0.004, and α
3
∈
(
0,0.05
) (see Figure 17).Figure 14
The change ofl
2
(
α
1
,
α
2
,
α
3
) with α
3, α
1, and α
2.Figure 15
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.048.Figure 16
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.02.Figure 17
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.004.(2) Fix the value ofα
2 and α
1
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
2
=
0.004, α
1
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 18).
(ii)
α
2
=
0.004, α
3
=
0.01, and α
1
∈
(
0,0.05
) (see Figure 19).
(iii)
α
2
=
0.004, α
1
∈
(
0,0.05
), and α
3
=
0.03 (see Figure 20).
(iv)
α
2
=
0.004, α
3
=
0.05, and α
1
∈
(
0,0.05
) (see Figure 21).Figure 18
The change ofl
2
(
α
1
,
α
2
,
α
3
) with α
2, α
1, and α
3.Figure 19
The change ofl
2
(
α
1
) with α
1 when α
2
=
0.004
,
α
2
=
0.01.Figure 20
The change ofl
2
(
α
1
) with α
1 when α
2
=
0.004, α
2
=
0.03.Figure 21
The change ofl
2
(
α
1
) with α
1 when α
2
=
0.004, α
2
=
0.05.(3) Fix the value ofα
3 and α
1
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
3
=
0.03, α
1
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 22).
(ii)
α
3
=
0.03, α
1
=
0.01, and α
2
∈
(
0,0.05
) (see Figure 23).
(iii)
α
3
=
0.03, α
1
=
0.02, and α
2
∈
(
0,0.05
) (see Figure 24).
(iv)
α
3
=
0.03, α
1
=
0.03, and α
2
∈
(
0,0.05
) (see Figure 25).Figure 22
The change ofl
2
(
α
1
,
α
2
,
α
3
) with α
2, α
1, and α
3.Figure 23
The change ofl
2
(
α
2
) with α
2 when α
1
=
0.01, α
3
=
0.03.Figure 24
The change ofl
2
(
α
2
) with α
2 when α
1
=
0.02, α
3
=
0.03.Figure 25
The change ofl
2
(
α
2
) with α
2 when α
1
=
0.03, α
3
=
0.03.(
3
) When {
w
=
α
3
f
1
(
α
3
)
=
0
f
2
(
α
3
)
=
0, matrix J has pure imaginary eigenvalues ±
i
α
3. Calculate q, q
-, p, p
-.The result is as follows:(23)
L
3
(
α
1
,
α
2
,
α
3
)
=
1
2
ω
[
〈
p
,
C
(
q
,
q
,
q
-
)
〉
-
2
〈
p
,
B
(
q
,
A
-
1
B
(
q
,
q
-
)
)
〉
+
〈
p
,
B
(
q
-
,
(
2
i
ω
E
-
A
)
-
1
B
(
q
,
q
)
)
〉
]
=
408.75
α
1
α
2
-
29.6
i
α
2
α
3
α
2
(
-
168.9853
α
2
+
0.821
i
α
3
)
+
{
(
7183233.84
+
87766
i
)
α
1
2
α
2
+
(
61899.185
+
631.33
i
)
α
1
α
2
2
-
(
414.48
-
33923.1
i
)
α
1
2
α
3
-
(
13181
-
1085129.6
i
)
α
1
α
2
α
3
-
(
46
-
4510.5
i
)
α
2
2
α
3
-
(
5124.6
+
62.25
i
)
α
1
α
3
2
}
×
{
-
(
658.6
+
8
i
)
α
3
2
]
(
3295.26
α
1
α
2
+
15.562
i
α
1
α
3
)
×
[
-
(
658.6
+
8
i
)
α
3
2
(
375716.354
+
4737.87
i
)
α
1
α
2
-
(
117.78
+
1.2
i
)
α
2
2
-
(
53.3
-
4359.74
i
)
α
1
α
3
-
(
711.552
-
56821.5
i
)
α
2
α
3
-
(
658.6
+
8
i
)
α
3
2
]
}
-
1
-
2
[
-
2179.867
α
1
375716.354
α
1
α
2
-
117.78
α
2
2
+
408.75
α
1
α
2
-
29.6
i
α
2
α
3
α
2
(
-
168.9853
α
2
+
0.821
i
α
3
)
+
(
61899.185
α
1
α
2
-
4510.5
i
α
2
α
3
)
×
(
(
-
375716.354
α
1
+
117.78
α
2
)
×
(
3295.26
α
1
α
2
-
15.562
i
α
1
α
3
)
)
-
1
-
2179.867
α
1
375716.354
α
1
a
2
-
117.78
α
2
2
]
=
l
3
(
α
1
,
α
2
,
α
3
)
+
m
3
(
α
1
,
α
2
,
α
3
)
,
where L
3
(
α
1
,
α
2
,
α
3
) is the Lyapunov exponent function about α
1, α
2, and α
3 and l
3
(
α
1
,
α
2
,
α
3
) is the real part. The following is the figure of l
3
(
α
1
,
α
2
,
α
3
).(1) Fix the value ofα
1 and α
2
∈
(
0,0.05
)
,
α
3
∈
(
0,0.05
).(i)
α
1
=
0.02, α
2
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 26).
(ii)
α
1
=
0.02, α
2
=
0.048, and α
3
∈
(
0,0.05
) (see Figure 27).
(iii)
α
1
=
0.02, α
2
=
0.02, and α
3
∈
(
0,0.05
) (see Figure 28).
(iv)
α
1
=
0.02, α
2
=
0.004, and α
3
∈
(
0,0.05
) (see Figure 29).Figure 26
The change ofl
3
(
α
1
,
α
2
,
α
3
) with α
1, α
2, and α
3.Figure 27
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.048.Figure 28
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.02.Figure 29
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.004.(2) Fix the value ofα
2 and α
1
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
2
=
0.004, α
1
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 30).
(ii)
α
2
=
0.004, α
3
=
0.01, and α
1
∈
(
0,0.05
) (see Figure 31).
(iii)
α
2
=
0.004, α
3
=
0.03, α
1
∈
(
0,0.05
) (see Figure 32).
(iv)
α
2
=
0.004, α
3
=
0.05, and α
1
∈
(
0,0.05
) (see Figure 33).Figure 30
The change ofl
3
(
α
1
,
α
2
,
α
3
) with α
1, α
2, and α
3.Figure 31
The change ofl
3
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.01.Figure 32
The change ofl
3
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.03.Figure 33
The change ofl
3
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.05.(3) Fix the value ofα
3 and α
1
∈
(
0,0.05
), α
2
∈
(
0,0.05
).(i)
α
3
=
0.03, α
1
∈
(
0,0.05
), and α
2
∈
(
0,0.05
) (see Figure 34).
(ii)
α
3
=
0.03, α
1
=
0.01, and α
2
∈
(
0,0.05
) (see Figure 35).
(iii)
α
3
=
0.03, α
1
=
0.02, and α
2
∈
(
0,0.05
) (see Figure 36).
(iv)
α
3
=
0.03, α
1
=
0.03, and α
2
∈
(
0,0.05
) (see Figure 37).Figure 34
The change ofl
3
(
α
1
,
α
2
,
α
3
) with α
1, α
2, and α
3.Figure 35
The change ofl
3
(
α
2
) with α
2 when α
1
=
0.01, α
3
=
0.03.Figure 36
The change ofl
3
(
α
2
) with α
2 when α
1
=
0.02, α
3
=
0.03.Figure 37
The change ofl
3
(
α
2
) with α
2 when α
1
=
0.03, α
3
=
0.03.Premisingα
3
=
0.01, the stable area in the α
1, α
2 plane is determined by inequality group (10). Figure 38 shows the system fixed point in the area of the asymptotic stability.Figure 38
α
3
=
0.01, the stable area in α
1, α
2 plane.


#### 3.1.2. Numerical Simulation

Using Matlab, parameters influence on the system (8) can be analyzed through numerical simulation.(1)
α
1,
α
2,and
α
3
Influence on the Collecting Market. The second collector on the assumption that the manufacturer and the first collector stem variables parameters is fixed with α
1
=
0.01, α
3
=
0.01, her system variables parameters α
2 is at [
0,0.05
]. Manufacturers’ initial collecting price p
m, collector 1 and 2 for manufacturers of the collecting price of depreciate rate for 20, 0.35, respectively, 0.3. And, as shown in Figure 39 changes Figures 39(a), 39(b), and 39(c) shows, the corresponding Lyapunov index as shown in Figure 39(d). Figure 39 showing, when 0
<
α
2
<
0.0097, the system is in stable state. after multi game, p
m, r
1 and r
2 is stable at the point of (17.4, 0.4904, 0.4825). When α
1
=
0.01
,
α
2
=
0.0097
,
α
3
=
0.0097, the system (8) occurs the first bifurcation, then after cycle 2, 4 cycle, and so forth, the system is gradually into the chaotic state. Figure 39(d) is Lyapunov index spectrum distribution which can also confirmed the phenomena. at α
2
=
0.0099, the max Lyapunov index is zero, the system is in critical condition. At α
2
=
0.0133, most of Lyapunov index greater than zero which explains system at chaos state (Figure 39 and Figure 44).The max Lyapunov index.
(a)
α
1
=
0.01, α
3
=
0.01, and α
2
=
[
0.005
], p
m,
(b)
α
1
=
0.01, α
3
=
0.01
,  and α
2
=
[
0.005
], r
1
,
(c)
α
1
=
0.01, α
3
=
0.01, and α
2
=
[
0.005
], r
2,
(d)
α
1
=
0.01, α
3
=
0.01,   and α
2
=
[
0.005
]When initial setup and other parameters are the same andα
1, α
3, and α
3 are adjusted in the interval [
0,0.05
] simultaneously, system variables of the manufacturers and two collectors change as shown in Figures 40(a), 40(b), and 40(c), and Lyapunov index is as in Figure 40(d). When 0
<
α
2
<
0.0098, the system is in stable state at the point of (17.4, 0.4904, 0.4825). When the value of α
1, α
2, and α
3 is 0.0098, the system (8) occurs the first bifurcation, then after cycle 2, 4 cycle, and so forth, the system is gradually into the chaotic state. Figure 40(d) is Lyapunov index spectrum distribution which can also confirm the phenomena. At α
2
=
0.0099, the max Lyapunov index is zero, and the system is in critical condition. Comparing with Figure 39, system steady equilibrium does not change, keeps for (17.4, 0.4904, 0.4825).(a)α
1
,
α
2
,
α
3 is [
0
,
0.05
], p
1, (b) α
1
,
α
2
,
α
3 is [
0
,
0.05
], r
1, (c) α
1
,
α
2
,
α
2 is [
0
,
0.05
], r
2, and (d) α
1
,
α
2
,
α
3 is [
0
,
0.05
], the max Lyapunov index.
(a)
(b)
(c)
(d)(2) System Variables Power Spectrum. According to the numerical simulation, not only system variables graph in the phase space along with time can be drawn, but also power spectrum graph of system variables can be estimated by period chart method. When α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, power spectrum graph of the manufacturer and two collectors is shown as in Figures 41(a), 41(b), and 41(c).(a)α
1
=
0.01,  α
2
=
0.014, and α
3
=
0.01, power spectrum of p
m; (b) α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, power spectrum of r
1; and (c)α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, power of spectrum r
2.
(a)
(b)
(c)According to the numerical results, no matter how large the collection price adjustment parameters are, the collection price traverses the whole value area over time, butp
m, r
1, and r
2 in the system (10) always limit in a certain range, which also verifies boundedness and ergodicity of chaos. Because α
2
>
α
3 and r
1 amplitude is larger than r
2, the adjustment of system variables must be limited in a certain range.(3) The Influence of Initial Setup. The butterfly effect is an important symbol of chaotic motions, which is sensitive to initial state, and small changes of the initial conditions may cause the adjacent orbital evolution to index form separate after multiple episodes. For the above conditions, when α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, sensitive dependence on initial value of p
m, r
1, and r
2 is as shown in Figures 42(a), 42(b), and 42(c). p
m value is, respectively, taken 20 and 20.001, after 96 cycles iteration, p
m difference value is 0.1335 which is 133.5 times than the initial difference value. r
1 and r
2 value is, respectively, taken 0.35 and 0.351, 0.3 and 0.301, respectively, experience in 44 and 45 iteration, r
1 and r
2 difference values is respective 0.2826 and 0.1364 which are 282 times and 136 times than the initial difference values. Furthermore, from Figure 42 we can conclude that system variables values are approximate in the first period, with increase of iterations the difference increase obviously.(a)α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, the initial value sensitivity of p
m, (b) α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, the initial value sensitivity of r
1; and (c) α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, the initial value sensitivity r
2.
(a)
(b)
(c)


### 3.2. Stackelberg Equilibrium

#### 3.2.1. Model and the Analysis

Suppose that the manufacturer and collectors have principal and subordinate relationship, the manufacturer is Stackelberg leader, and collectors are followers; then they process sequential dynamic game; the game equilibrium is Stackelberg equilibrium. In this game, the manufacturer makes the decision of sales price and collection price according to the market informa Lyapunov exponent function abouttion; then two collectors make decision according to the decision-making of the manufacture:(24)
max
⁡
∏
r
1
=
(
p
m
-
p
r
1
(
t
)
-
A
1
)
(
k
1
+
β
1
p
r
1
(
t
)
-
δ
1
p
r
2
(
t
)
)
max
⁡
∏
r
2
=
(
p
m
-
p
r
2
(
t
)
-
A
2
)
(
k
2
+
β
2
p
r
2
(
t
)
-
δ
2
p
r
1
(
t
)
)
s
.
t
.
i
p
m
∈
arg
max
⁡
∏
r
1
s
.
t
.
i
p
m
∈
arg
max
⁡
∏
r
2
.The result is(25)
r
1
*
*
=
(
p
m
(
2
β
1
β
2
-
2
β
2
δ
1
+
δ
1
β
2
-
δ
1
δ
2
)
+
2
β
2
k
1
+
2
β
1
β
2
A
1
+
k
2
δ
1
+
β
2
A
2
δ
1
)
×
(
p
m
(
4
β
1
β
2
-
δ
1
δ
2
)
)
-
1
,
r
2
*
*
=
(
+
2
β
1
k
2
+
2
β
1
β
2
A
2
+
k
1
δ
2
+
β
1
A
2
δ
2
)
p
m
(
2
β
1
β
2
-
2
β
1
δ
2
+
δ
2
β
1
-
δ
1
δ
2
)
+
2
β
1
k
2
+
2
β
1
β
2
A
2
+
k
1
δ
2
+
β
1
A
2
δ
2
)
×
(
p
m
(
4
β
1
β
2
-
δ
1
δ
2
)
)
-
1
.Formula (25) is reaction function of the collectors. Put formula (25) into formula (2) and the following can be obtained:
(26)
p
m
*
*
=
(
p
0
-
c
r
)
2
-
(
(
β
1
-
δ
2
)
(
2
β
2
k
1
+
2
β
2
β
1
A
1
+
δ
1
k
2
+
δ
1
β
2
A
2
)
+
(
β
2
-
δ
1
)
(
2
β
1
k
2
+
2
β
2
β
1
A
2
+
δ
2
k
1
-
δ
2
β
1
A
1
)
-
(
4
β
1
β
2
-
δ
1
δ
2
)
(
k
1
+
k
2
)
)
×
(
(
2
β
1
β
2
+
2
β
2
δ
1
-
δ
1
β
2
)
(
β
2
-
δ
1
)
+
(
2
β
1
β
2
+
2
β
1
δ
2
-
δ
2
β
1
)
(
β
1
-
δ
2
)
)
-
1
.With Nash equilibrium value the parameters values arep
0
=
70, c
r
=
30, A
1
=
1, A
2
=
0.8, k
1
=
1, k
2
=
1, β
1
=
0.6, β
2
=
0.55, δ
1
=
0.2, δ
2
=
0.2. Then we put the values into Stackelberg equilibrium (
p
m
*
*
,
r
1
*
*
,
r
2
*
*
); the equilibrium value is (
19.17
,
0.48
,
0.48
).The manufacturer makes decision based on limited rational expectations. She adjusts the game process on the basis of marginal gains. If the marginal profits oft period are positive, in the t
+
1 period this strategy action will be used. And if the t period marginal profits are negative, in the t
+
1 period the collection price will be lowered. The repeated game model of dynamic adjustment is
(27)
p
m
(
t
+
1
)
=
p
m
(
t
)
+
α
p
m
(
t
)
×
(
(
(
1
-
r
1
)
(
β
1
-
δ
2
)
+
(
1
-
r
2
)
(
β
2
-
δ
1
)
)
×
(
p
0
-
2
p
m
(
t
)
-
c
r
)
-
k
1
-
k
2
(
(
1
-
r
1
)
(
β
1
-
δ
2
)
+
(
1
-
r
2
)
(
β
2
-
δ
1
)
)
)
.In formula (27), α is adjustment coefficient of collection price and r
1
,
r
2 are as follows:
(28)
r
1
(
t
)
=
(
2
β
2
k
1
+
2
β
1
β
2
A
1
+
k
2
δ
1
+
β
2
A
2
δ
1
p
m
(
t
)
(
2
β
1
β
2
-
2
β
2
δ
1
+
δ
1
β
2
-
δ
1
δ
2
)
+
2
β
2
k
1
+
2
β
1
β
2
A
1
+
k
2
δ
1
+
β
2
A
2
δ
1
)
×
(
p
m
(
t
)
(
4
β
1
β
2
-
δ
1
δ
2
)
)
-
1
,
(29)
r
2
(
t
)
=
(
2
β
1
k
2
+
2
β
1
β
2
A
2
+
k
1
δ
2
+
β
1
A
2
δ
2
p
m
(
t
)
(
2
β
1
β
2
-
2
β
1
δ
2
+
δ
2
β
1
-
δ
1
δ
2
)
+
2
β
1
k
2
+
2
β
1
β
2
A
2
+
k
1
δ
2
+
β
1
A
2
δ
2
)
×
(
p
m
(
t
)
(
4
β
1
β
2
-
δ
1
δ
2
)
)
-
1
.
Substitution parameters value into system (17) can be concluded which is as follows:
(30)
p
m
(
t
+
1
)
=
p
m
(
t
)
+
α
p
m
(
t
)
(
17.403
-
0.908
p
m
(
t
)
)
,
r
1
(
t
)
=
0.398
+
1.617
p
m
(
t
)
,
r
2
(
t
)
=
0.391
+
1.703
p
m
(
t
)
.
From system (30) it can be concluded that the manufacture firstly makes decision and the decision variables directly relate to α. But collectors’ decision variables directly relate to p
m. The following simulates complex dynamics characteristics of system (30) through numerical simulation.


#### 3.2.2. Numerical Simulation

(1)
α
Influence on the Collecting Market. When α
=
[
0,0.05
], the initial values of p
m, r
1, and r
2 are 20, 0.4, and 0.5; p
m, r
1, and r
2 are as shown in Figures 43(a), 43(b), and 43(c). From the figures, we can conclude that the system is in stable state when 0
<
α
<
0.113, and the stable values of p
m, r
1, and r
2 are 19.17, 0.4824, and 0.4799. When α
=
0.113, the system occurs the first bifurcation, then after cycle 2, 4 cycle, and so forth, the system is gradually into the chaotic state.(a)α = [0, 0.05], p
m, (b) α = [0, 0.05], r
1; and (c) α =  [0, 0.05], r
2.
(a)
(b)
(c)Figure 44
α =  [0, 0.05], max Lyapunov spectrum.(2) System Variables Power Spectrum. According to the numerical simulation, not only system variables graph in the phase space along with time can be drew, but also power spectrum graph of system variables can be estimated by period chart method. When α
=
0.15, power spectrum graph of the manufacturer and two collectors is shown in Figures 45(a), 45(b), and 45(c). According to the numerical results, no matter how large the collection price adjustment parameters are, the collection price traverses the whole value area over time, but p
m, r
1, and r
2 in the system (19) always limit in a certain range, which also verifies boundedness and ergodicity of chaos.(a)α
=
0.15, the power spectrum of p
m, (b) α
=
0.15, the spectrum of r
1; and (c) α
=
0.15, the power spectrum of r
2.
(a)
(b)
(c)(3) The Influence of Initial Setup. When α
=
0.15, sensitive dependence on initial value of p
m, r
1, and r
2 is as shown in Figures 46(a), 46(b), and 46(c). p
m values are, respectively, taken to be 20 and 20.001, after 36 cycles iteration, p
m difference value is 3.35 which is 3350 times the initial difference value. r
1 and r
2 difference values are 69 times and 73 times the initial difference values. Furthermore, from Figure 46 we can conclude that system variables values are approximate in the first period; with increase of iterations the difference increases obviously. The influence to the manufacturer is far greater than collectors because of the collectors’ decision-making later than the manufacturer’s. After the collectors observe action of the manufacturer, they make decision. The collectors’ decision-making is superior to the manufacturer.(a)α
=
0.15, the initial value sensitivity of p
m, (b) α
=
0.15, the initial value sensitivity of r
1
; and (c) α
=
0.15, the initial value sensitivity of r
2.
(a)
(b)
(c)


### 3.3. Model Comparison and Analysis

(1) Profits Analysis of Nash Equilibrium. When α
=
[
0,0.05
], profits of the manufacturer and two collectors are as shown in Figures 47(a), 47(b), 47(c), and 47(d). Figure 47(a) shows that profits of the manufacturer in the stable state are not the best; it means Nash equilibrium is not optimal solution; the optimal solution of single cycle may not bring enterprise maximum returns. Furthermore, from Figure 47(a) it can be concluded that, with collecting price increasing, profits state instead into chaos state, the optimal price is 19.64 in the chaotic state, and the corresponding profits are 208.6. Compared to Figure 47(b), profits of the manufacture before chaotic state are 216.4, so the internal randomness of chaotic system cannot react to the actual situation of the market.(a)∏
m graph with p
m, (b) ∏
m graph with α
1, (c) ∏
r
1 graph with α
2, and (d) ∏
r
2 graph with α
3.
(a)
(b)
(c)
(d)(2) Profits Analysis of Stackelberg Equilibrium. When α
=
[
0,0.05
], profits of the manufacturer and two collectors are as shown in Figures 48(a), 48(b), 48(c), and 48(d). Figure 48 shows that profits of the manufacturer in the stable state are the best; it means Stackelberg equilibrium is the optimal solution, and the optimal solution of the closed-loop supply chain is optimal. Figures 48(a), 48(b), 48(c), and 48(d) are consistent.(a)∏
m graph with p
m, (b) ∏
m graph with α
1, (c) ∏
r
1 graph with α
2, and (d) ∏
r
2 graph with α
3.
(a)
(b)
(c)
(d)Through numerical simulation and comparison analysis, it can be concluded that Stackelberg game equilibrium is better than Nash equilibrium from the point of view of the profits and collecting price. Three conclusions can be summarized as follows. First, when the manufacturer and the two collectors make decision by static game, the collection price and profits of the three parties are lower than the optimal of dynamic game equilibrium. Second, as the manufacturer and the two collectors independently make decisions, the three parties will try to lower the collecting price in order to obtain the maximum profits. Third, as the manufacturer and the two collectors independently make decisions, the system more easily reaches chaos, and the cycle of stable state is shorter.


## 3.1. Nash Equilibrium

### 3.1.1. Model and Analysis

(
1) Model. The manufacturer and collectors are equal in closed-loop supply chain; that is, the two sides make decision at the same time: manufacturer’s decision is to choose the collection price p
i to maximize returns, and collectors’ decision is to choose p
r
i to maximize their returns.Nash equilibrium can be concluded by the first-order conditions of three reaction functions, which meet∂
Π
m
/
∂
p
m
=
0, ∂
Π
r
i
/
∂
r
i
=
0 simultaneously. The Nash equilibrium is the optimal decision of a game in all kinds of may guess, thus make it to gain maximum benefit. The Nash equilibrium is
(5)
(
p
m
*
,
r
1
*
,
r
2
*
)
=
(
-
B
+
B
2
-
4
A
C
2
A
,
t
1
+
t
2
p
m
*
,
t
3
+
t
4
p
m
*
)
,
where A, B, and C are
(6)
A
=
2
(
β
1
+
β
2
-
δ
1
-
δ
2
)
(
p
0
-
c
r
)
h
h
h
-
2
t
1
(
β
1
-
δ
2
)
-
2
t
3
(
β
2
-
δ
1
)
,
B
=
k
1
+
k
2
-
(
β
1
+
β
2
-
δ
1
-
δ
2
)
(
p
0
-
c
r
)
h
h
h
+
(
p
0
-
c
r
)
(
t
1
(
β
1
-
δ
2
)
+
t
3
(
β
2
-
δ
1
)
)
h
h
h
-
2
t
2
(
β
1
-
δ
2
)
-
2
t
4
(
β
2
-
δ
1
)
,
C
=
(
p
0
-
c
r
)
t
2
(
β
1
-
δ
2
)
+
t
4
(
β
2
-
δ
1
)
.In formula (6), t
1, t
2, t
3, t
4 can be expressed as
(7)
t
1
=
2
β
1
β
2
-
2
β
2
δ
1
+
δ
1
β
2
-
δ
1
δ
2
4
β
1
β
2
-
δ
1
δ
2
,
t
2
=
2
β
2
k
1
+
2
β
2
β
1
A
1
+
δ
1
k
2
-
δ
1
β
2
A
2
4
β
1
β
2
-
δ
1
δ
2
,
t
3
=
2
β
1
β
2
-
2
β
1
δ
2
+
δ
2
β
1
-
δ
1
δ
2
4
β
1
β
2
-
δ
1
δ
2
,
t
4
=
2
β
1
k
2
+
2
β
2
β
1
A
2
+
δ
2
k
1
-
δ
2
β
1
A
1
4
β
1
β
2
-
δ
1
δ
2
.In reality, the game between node enterprises in closed-loop supply chain is continuous, enterprises’ decision-making is a long-term repeated process, and its action has long-term memory. And each node enterprise does not completely control the market information and also cannot fully expect future market changes, so based on limited rational expectations decision we adjust process with marginal gains. They can make the next-period price decision on the basis of the local estimate to his marginal profits in current period. Their price adjustment processes are(8)
p
m
(
t
+
1
)
=
p
m
(
t
)
+
α
1
p
m
(
t
)
×
(
(
(
1
-
r
1
(
t
)
)
(
β
1
-
δ
2
)
+
(
1
-
r
2
(
t
)
)
(
β
2
-
δ
1
)
)
×
(
p
0
-
2
p
m
(
t
)
-
c
r
)
-
k
1
-
k
2
(
(
1
-
r
1
(
t
)
)
(
β
1
-
δ
2
)
+
(
1
-
r
2
(
t
)
)
(
β
2
-
δ
1
)
)
)
,
r
1
(
t
+
1
)
=
r
1
(
t
)
+
α
2
r
1
(
t
)
×
(
-
2
β
1
r
1
(
t
)
p
m
2
(
t
)
+
p
m
(
t
)
h
h
h
×
(
k
1
+
(
β
1
-
δ
1
+
δ
1
r
2
(
t
)
)
h
h
h
×
p
m
(
t
)
+
β
1
A
1
)
(
-
2
β
1
r
1
(
t
)
p
m
2
(
t
)
+
p
m
(
t
)
)
,
r
2
(
t
+
1
)
=
r
2
(
t
)
+
α
3
r
2
(
t
)
×
(
-
2
β
2
r
2
(
t
)
p
m
2
(
t
)
+
p
m
(
t
)
h
h
×
(
k
2
+
(
β
2
-
δ
2
+
δ
2
r
1
(
t
)
)
p
m
(
t
)
+
β
2
A
2
)
(
-
2
β
2
r
2
(
t
)
p
m
2
(
t
)
+
p
m
(
t
)
)
,
where p
1
(
t
+
1
),  r
1
(
t
+
1
), r
2
(
t
+
1
), and p
1
(
t
), r
1
(
t
), r
2
(
t
) are used to present the price in the (
t
+
1
)th period and in the (
t
)th period. Where α
i are positive parameters, i
=
1,2
,
3, which denotes adjustment speed, respectively, for the manufacturer and two collectors.From the system (8) we can be conclude that the optimal collecting price of the manufacturer is related to the collection price adjustment coefficient, the two collectors’ collecting price, sale price, and consumers’ environmental protection awareness. Similarly, the optimal collecting price of collectors is related to collecting price adjustment parameters, collecting price of the manufacturer, collecting function, consumers’s environmental protection awareness, and collecting costs.(2) Model Analysis. In closed-loop supply chain, any enterprise decision-making is all according to the maximum profits, the equilibrium less than zero is not practical significance, such as manufacturers to make power mainly because of benefits, collecting business is also benefit with collecting. So only research system (8) is the equilibrium; system (8) of the eight fixed points, respectively, is E
0
(
0,0
,
0
),  E
1
(
p
m
1
,
0,0
),  E
2
(
0
,
r
1
2
,
0
),  E
3
(
0
,
0
,
r
2
3
),  E
4
(
p
m
4
,
r
1
4
,
0
), E
5
(
p
m
5
,
0
,
r
2
5
), E
6
(
0
,
r
1
6
,
r
2
6
), E
7
(
p
m
*
,
r
1
*
,
r
2
*
). The analysis shows that by E
0, E
1, E
2, E
3, E
4, E
5 and E
6, is no real significance. Therefore, we discuss Nash equilibrium which is E
7
(
p
m
*
,
r
1
*
,
r
2
*
); E
7 is the reaction function intersection of the manufacturer and collectors, which means marginal profits of both sides are zero, but this does not mean that the result of the game will tend to equilibrium. Instead, a party rational behavior change may cause game process to occur very complex phenomenon. E
7 has local stability, this stability region by α
1, α
2, α
3 decided. The stability of E
7 will be studied. The first we calculate Jacobi matrix of system (8), means
(9)
J
=
(
J
11
J
12
J
13
J
21
J
22
J
23
J
31
J
32
J
33
)
.We putE
7
(
p
m
*
,
r
1
*
,
r
2
*
) into formula (9), according to the Jury stability criterion, the sufficient and necessary conditions of E
7 asymptotic stability are that the zeros of characteristic polynomial F
(
λ
)
=
λ
3
+
A
λ
2
+
B
λ
+
C are in a circular area. Moreover, the characteristic polynomial will satisfy the conditions as follows:
(10)
F
(
1
)
=
1
+
A
+
B
+
C
>
0
,
(
-
1
)
3
F
(
-
1
)
=
1
-
A
+
B
-
C
>
0
,
1
>
|
C
|
,
1
-
C
2
>
|
B
-
A
C
|
.At present, only a few simple dynamic systems are analyzed with analytical method, but complex dynamic system mainly uses the numerical analysis method. This paper processes numerical simulation on system (8) by Matlab to expresses the dynamic characteristics. Because remanufacturing is still initial stage in china, consumers environmental protection awareness is lower, and manufacturing cost (including the fixed and variable cost) is higher; for autoparts collecting, the parameters of system (8) can be defined as p
0
=
70, c
r
=
30, A
1
=
1, A
2
=
0.8, k
1
=
1, k
2
=
1, β
1
=
0.6, β
2
=
0.55, δ
1
=
0.2, and δ
2
=
0.2 to study the local stability of equilibrium point. Get the parameter value into E
7
(
p
m
*
,
r
1
*
,
r
2
*
) and formula (10) to conclude the value for (17.385, 0.491, 0.489).Change the discrete system (8) into a continuous system and the points of Jacobi matrix for
(11)
J
=
[
-
13.317
α
1
-
17.755
α
1
-
15.562
α
1
-
0.775
α
2
-
177.888
α
2
29.68
α
2
-
0.821
α
3
29.599
α
3
-
164.65
α
3
]
.Consider the continuous-time nonlinear dynamical system(12)
x
˙
=
A
x
+
N
(
x
)
,
A
=
(
a
i
j
)
n
×
n
,
x
∈
R
n
.Let the functionN
(
x
) be written as
(13)
N
(
x
)
=
1
2
B
(
x
,
x
)
+
1
6
C
(
x
,
x
,
x
)
+
O
(
∥
x
∥
4
)
,
where B
(
x
,
y
) and C
(
x
,
y
,
z
) are bilinear and trilinear functions, respectively. In coordinates, we have
(14)
B
i
(
x
,
y
)
=
∑
j
,
k
=
1
n
∂
2
N
i
(
ξ
)
∂
ξ
j
∂
ξ
k
|
ξ
=
0
x
j
y
k
,
C
i
(
x
,
y
,
z
)
=
∑
j
,
k
,
l
=
1
n
∂
3
N
i
(
ξ
)
∂
ξ
j
∂
ξ
k
∂
ξ
l
|
ξ
=
0
x
j
y
k
z
l
.Supposing thatA has a pair of complex eigenvalues on the imaginary axis, they are λ
1,2
=
±
i
ω
(
ω
>
0
), and these eigenvalues are the only eigenvalues with Re
λ
=
0. Let q
∈
C
n be a complex eigenvector to λ
1
=
i
ω:
(15)
A
q
=
i
ω
q
,
A
q
¯
=
-
i
ω
q
¯
.And the adjoint eigen vectorp
∈
C
n admits the properties
(16)
A
T
p
=
-
i
ω
p
,
A
T
p
¯
=
i
ω
p
¯and satisfies the normalization〈
q
,
p
〉
=
1.The first Lyapunov coefficient at the origin is defined by(17)
l
1
=
1
2
ω
Re
[
〈
p
,
C
(
q
,
q
,
q
¯
)
〉
-
2
〈
p
,
B
(
q
,
A
-
1
B
(
q
,
q
¯
)
)
〉
hhhhhhhh
+
〈
p
,
B
(
q
¯
,
(
2
i
ω
E
-
A
)
-
1
B
(
q
,
q
)
)
〉
]
.Next, we calculatei
w
q, J
q
-
=
-
i
w
q
-, J
T
p
=
-
i
w
p, and  J
T
p
-
=
i
w
p
-; If q
=
(
q
1
,
q
2
,
q
3
)
T
≠
0, calculate q
,
q
-
,
p
,
p
-:
(18)
[
-
13.317
α
1
-
i
w
-
17.755
α
1
-
15.562
α
1
-
0.775
α
2
-
177.888
α
2
-
i
w
29.68
α
2
-
0.821
α
3
29.599
α
3
-
164.65
α
3
-
i
w
]
×
[
q
1
q
2
q
3
]
=
[
0
0
0
]
having nonzero solutions.It means(19)
|
-
13.317
α
1
-
i
w
-
17.755
α
1
-
15.562
α
1
-
0.775
α
2
-
177.888
α
2
-
i
w
29.68
α
2
-
0.821
α
3
29.599
α
3
-
164.65
α
3
-
i
w
|
=
0
.It can be got from calculating(20)
f
1
(
w
)
=
-
375716.354
α
1
α
2
α
3
+
117.77755
α
2
2
α
3
+
13.317
α
1
w
2
+
177.888
α
2
w
2
+
164.65
α
3
w
2
=
0
,
f
2
(
w
)
=
-
2368.9345
α
1
α
2
w
+
0.6
α
2
2
w
-
28410.761
α
2
α
3
w
-
2179.87
α
1
α
3
w
+
w
3
=
0
.(1) When{
w
=
α
1
f
1
(
w
)
=
0
f
2
(
w
)
=
0, matrix J has pure imaginary eigenvalues ±
i
α
1. Calculate q, q
-, p, p
-.The result is as follows.(21)
L
2
(
α
1
,
α
2
,
α
3
)
=
1
2
ω
[
〈
p
,
C
(
q
,
q
,
q
-
)
〉
-
2
〈
p
,
B
(
q
,
A
-
1
B
(
q
,
q
-
)
)
〉
h
h
h
h
h
h
+
〈
p
,
B
(
q
-
,
(
2
i
ω
E
-
A
)
-
1
B
(
q
,
q
)
)
〉
]
=
(
408.75
-
29.6
i
)
α
1
α
2
(
0.82
i
α
1
-
168.9853
α
2
)
α
2
-
2
[
(
408.75
-
29.6
i
)
α
1
α
2
(
0.82
i
α
1
-
168.9853
α
2
)
α
2
-
2179.87
α
1
375716.354
α
1
α
2
-
117.7775
α
2
2
+
(
+
α
1
α
2
(
1.2
i
α
2
+
(
375716.354
+
56821.5
i
)
α
3
)
}
)
-
1
(
61899.185
-
4510.5
i
)
α
2
×
(
5846897.898
i
α
1
2
-
(
1.24
×
1
0
9
+
1832.854
i
)
α
1
α
2
+
388107.8
α
2
2
)
-
1
(
(
61899.185
-
4510.5
i
)
α
2
)
(
408.75
-
29.6
i
)
α
1
α
2
(
0.82
i
α
1
-
168.9853
α
2
)
α
2
]
+
(
+
α
1
α
2
(
1.2
i
α
2
+
(
375716.354
+
56821.5
i
)
α
3
)
}
)
-
1
{
α
1
[
(
-
4
+
26.634
i
)
α
1
h
h
h
h
h
h
h
h
+
(
2179.87
+
329.3
i
)
α
3
]
}
×
(
{
(
-
53
-
8
i
)
α
1
3
+
α
1
2
[
(
-
711.55
+
4737.8
i
)
α
2
-
(
658.6
-
4359.7
i
)
α
3
]
-
117.8
α
2
2
α
3
+
α
1
α
2
(
1.2
i
α
2
+
(
375716.354
+
56821.5
i
)
α
3
)
(
-
53
-
8
i
)
α
1
3
}
)
-
1
)
+
{
(
-
(
407.31
+
29.7
i
)
α
2
)
×
(
-
1.55
i
α
1
α
2
-
151.97
α
2
α
3
)
}
×
(
{
(
15.6
i
α
1
+
3295.3
α
2
)
(
-
53
-
8
i
)
α
1
3
×
{
(
-
53.3
-
8
i
)
α
1
3
(
-
53
-
8
i
)
α
1
3
+
α
1
2
[
(
-
711.552
+
4737.9
i
)
α
2
-
(
658.6
-
4359.74
i
)
α
3
]
-
117.77755
α
2
2
α
3
+
α
1
α
2
×
[
1.2
i
α
2
+
(
375716.354
+
56821.52
i
)
α
3
]
(
-
53
-
8
i
)
α
1
3
}
}
)
-
1
=
l
1
(
α
1
,
α
2
,
α
3
)
+
m
1
(
α
1
,
α
2
,
α
3
)
,
where L
1
(
α
1
,
α
2
,
α
3
) is the Lyapunov exponent function about α
1, α
2, and α
3 and l
1
(
α
1
,
α
2
,
α
3
) is the real part. The following is the figure of l
1
(
α
1
,
α
2
,
α
3
).(1) Fix the value ofα
1 and α
2
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
1
=
0.02, α
2
=
0.048, and α
3
∈
(
0,0.05
) (see Figure 3).
(ii)
α
1
=
0.02, α
2
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 2).Figure 2
The change ofl
1
(
α
1
,
α
2
,
α
3
) with α
1
,
α
2, and α
3.Figure 3
The change ofl
1
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.048.We can know from Figure3 that when α
1 =  0.02, α
2 = 0.048, and α
3 is between zero and 0.05, l
1
(
α
3
) is always equal to 1.375.(iii)
α
1
=
0.02, α
2
=
0.02, and α
3
∈
(
0,0.05
) (see Figure 4).
(iv)
α
1
=
0.02, α
2
=
0.004, and α
3
∈
(
0,0.05
) (see Figure 5).Figure 4
The change ofl
1
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.02.Figure 5
The change ofl
1
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.004.From Figures4 and 5 we can observe the same change of l
1
(
α
3
) as in Figures 2 and 3.(2) Fix the value ofα
2 and α
1
∈
(
0,0.05
), α
3
∈
(
0,0.05
). α
2
=
0.004, α
1
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 6). α
2
=
0.004, α
3
=
0.01, and α
1
∈
(
0,0.05
) (see Figure 7). α
2
=
0.004, α
3
=
0.03, and α
1
∈
(
0,0.05
) (see Figure 8). α
2
=
0.004, α
3
=
0.05, and α
1
∈
(
0,0.05
) (see Figure 9).Figure 6
The change ofl
1
(
α
1
,
α
2
,
α
3
) with α
2, α
1, and α
3.Figure 7
The change ofl
1
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.01.Figure 8
The change ofl
1
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.03.Figure 9
The change ofl
1
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.05.(3) Fix the value ofα
3 and α
1
∈
(
0,0.05
), α
2
∈
(
0,0.05
).(i)
α
3
=
0.03, α
1
∈
(
0,0.05
), and α
2
∈
(
0,0.05
) (see Figure 10).
(ii)
α
3
=
0.03, α
1
=
0.01, and α
2
∈
(
0,0.05
) (see Figure 11).
(iii)
α
3
=
0.03, α
1
=
0.02, and α
2
∈
(
0,0.05
) (see Figure 12).
(iv)
α
3
=
0.03, α
1
=
0.03, and α
2
∈
(
0,0.05
) (see Figure 13).Figure 10
The change ofl
1
(
α
1
,
α
2
,
α
3
) with α
3, α
1, and α
2.Figure 11
The change ofl
1
(
α
2
) with α
2 when α
1
=
0.01, α
3
=
0.03.Figure 12
The change ofl
1
(
α
2
) with α
2 when α
1
=
0.02, α
3
=
0.03.Figure 13
The change ofl
1
(
α
2
) with α
2 when α
1
=
0.03, α
3
=
0.03.(2) When{
w
=
α
2
f
1
(
α
2
)
=
0
f
2
(
α
2
)
=
0, matrix J has pure imaginary eigenvalues ±
i
α
2. Calculate q, q
-, p, p
-.The result is as follows:(22)
L
2
(
α
1
,
α
2
,
α
3
)
=
1
2
ω
[
〈
p
,
C
(
q
,
q
,
q
-
)
〉
-
2
〈
p
,
B
(
q
,
A
-
1
B
(
q
,
q
-
)
)
〉
+
〈
p
,
B
(
q
-
,
(
2
i
ω
E
-
A
)
-
1
B
(
q
,
q
)
)
〉
]
=
(
0.000851
-
0.175
i
)
+
37.57
+
0.18
i
375716.354
α
1
-
117.8
α
2
+
(
2.42
+
0.012
i
)
α
1
α
2
+
(
0.01293
-
2.74
i
)
α
2
375716.354
α
1
2
-
117.78
α
1
α
2
+
4359.74
α
1
375716.354
α
1
α
2
-
117.77755
α
2
2
+
(
{
α
2
2
[
(
-
0.014
+
0.000066
i
)
α
2
+
(
0.0065
+
1.37
i
)
α
3
]
+
α
1
α
2
[
(
-
4
+
0.19
i
)
α
2
+
(
18.784
+
329.2
i
)
α
3
]
+
α
1
2
(
26.6
i
α
2
+
2179.87
α
3
)
}
)
×
{
α
1
α
2
[
(
-
53.268
+
4737.87
i
)
α
1
α
2
-
(
711.552
+
9.2
i
)
α
2
2
+
(
375716.354
+
4359.74
i
)
α
1
α
3
-
(
776.38
-
56821.52
i
)
α
2
α
3
]
}
-
1
=
l
2
(
α
1
,
α
2
,
α
3
)
+
m
2
(
α
1
,
α
2
,
α
3
)
,
where L
2
(
α
1
,
α
2
,
α
3
) is the Lyapunov exponent function about α
1, α
2, and α
3 and l
2
(
α
1
,
α
2
,
α
3
) is the real part. The following is the figure of l
2
(
α
1
,
α
2
,
α
3
).(1) Fix the value ofα
1 and α
2
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
1
=
0.02, α
2
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 14).
(ii)
α
1
=
0.02, α
2
=
0.048, and α
3
∈
(
0,0.05
) (see Figure 15).
(iii)
α
1
=
0.02, α
2
=
0.02, and α
3
∈
(
0,0.05
) (see Figure 16).
(iv)
α
1
=
0.02, α
2
=
0.004, and α
3
∈
(
0,0.05
) (see Figure 17).Figure 14
The change ofl
2
(
α
1
,
α
2
,
α
3
) with α
3, α
1, and α
2.Figure 15
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.048.Figure 16
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.02.Figure 17
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.004.(2) Fix the value ofα
2 and α
1
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
2
=
0.004, α
1
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 18).
(ii)
α
2
=
0.004, α
3
=
0.01, and α
1
∈
(
0,0.05
) (see Figure 19).
(iii)
α
2
=
0.004, α
1
∈
(
0,0.05
), and α
3
=
0.03 (see Figure 20).
(iv)
α
2
=
0.004, α
3
=
0.05, and α
1
∈
(
0,0.05
) (see Figure 21).Figure 18
The change ofl
2
(
α
1
,
α
2
,
α
3
) with α
2, α
1, and α
3.Figure 19
The change ofl
2
(
α
1
) with α
1 when α
2
=
0.004
,
α
2
=
0.01.Figure 20
The change ofl
2
(
α
1
) with α
1 when α
2
=
0.004, α
2
=
0.03.Figure 21
The change ofl
2
(
α
1
) with α
1 when α
2
=
0.004, α
2
=
0.05.(3) Fix the value ofα
3 and α
1
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
3
=
0.03, α
1
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 22).
(ii)
α
3
=
0.03, α
1
=
0.01, and α
2
∈
(
0,0.05
) (see Figure 23).
(iii)
α
3
=
0.03, α
1
=
0.02, and α
2
∈
(
0,0.05
) (see Figure 24).
(iv)
α
3
=
0.03, α
1
=
0.03, and α
2
∈
(
0,0.05
) (see Figure 25).Figure 22
The change ofl
2
(
α
1
,
α
2
,
α
3
) with α
2, α
1, and α
3.Figure 23
The change ofl
2
(
α
2
) with α
2 when α
1
=
0.01, α
3
=
0.03.Figure 24
The change ofl
2
(
α
2
) with α
2 when α
1
=
0.02, α
3
=
0.03.Figure 25
The change ofl
2
(
α
2
) with α
2 when α
1
=
0.03, α
3
=
0.03.(
3
) When {
w
=
α
3
f
1
(
α
3
)
=
0
f
2
(
α
3
)
=
0, matrix J has pure imaginary eigenvalues ±
i
α
3. Calculate q, q
-, p, p
-.The result is as follows:(23)
L
3
(
α
1
,
α
2
,
α
3
)
=
1
2
ω
[
〈
p
,
C
(
q
,
q
,
q
-
)
〉
-
2
〈
p
,
B
(
q
,
A
-
1
B
(
q
,
q
-
)
)
〉
+
〈
p
,
B
(
q
-
,
(
2
i
ω
E
-
A
)
-
1
B
(
q
,
q
)
)
〉
]
=
408.75
α
1
α
2
-
29.6
i
α
2
α
3
α
2
(
-
168.9853
α
2
+
0.821
i
α
3
)
+
{
(
7183233.84
+
87766
i
)
α
1
2
α
2
+
(
61899.185
+
631.33
i
)
α
1
α
2
2
-
(
414.48
-
33923.1
i
)
α
1
2
α
3
-
(
13181
-
1085129.6
i
)
α
1
α
2
α
3
-
(
46
-
4510.5
i
)
α
2
2
α
3
-
(
5124.6
+
62.25
i
)
α
1
α
3
2
}
×
{
-
(
658.6
+
8
i
)
α
3
2
]
(
3295.26
α
1
α
2
+
15.562
i
α
1
α
3
)
×
[
-
(
658.6
+
8
i
)
α
3
2
(
375716.354
+
4737.87
i
)
α
1
α
2
-
(
117.78
+
1.2
i
)
α
2
2
-
(
53.3
-
4359.74
i
)
α
1
α
3
-
(
711.552
-
56821.5
i
)
α
2
α
3
-
(
658.6
+
8
i
)
α
3
2
]
}
-
1
-
2
[
-
2179.867
α
1
375716.354
α
1
α
2
-
117.78
α
2
2
+
408.75
α
1
α
2
-
29.6
i
α
2
α
3
α
2
(
-
168.9853
α
2
+
0.821
i
α
3
)
+
(
61899.185
α
1
α
2
-
4510.5
i
α
2
α
3
)
×
(
(
-
375716.354
α
1
+
117.78
α
2
)
×
(
3295.26
α
1
α
2
-
15.562
i
α
1
α
3
)
)
-
1
-
2179.867
α
1
375716.354
α
1
a
2
-
117.78
α
2
2
]
=
l
3
(
α
1
,
α
2
,
α
3
)
+
m
3
(
α
1
,
α
2
,
α
3
)
,
where L
3
(
α
1
,
α
2
,
α
3
) is the Lyapunov exponent function about α
1, α
2, and α
3 and l
3
(
α
1
,
α
2
,
α
3
) is the real part. The following is the figure of l
3
(
α
1
,
α
2
,
α
3
).(1) Fix the value ofα
1 and α
2
∈
(
0,0.05
)
,
α
3
∈
(
0,0.05
).(i)
α
1
=
0.02, α
2
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 26).
(ii)
α
1
=
0.02, α
2
=
0.048, and α
3
∈
(
0,0.05
) (see Figure 27).
(iii)
α
1
=
0.02, α
2
=
0.02, and α
3
∈
(
0,0.05
) (see Figure 28).
(iv)
α
1
=
0.02, α
2
=
0.004, and α
3
∈
(
0,0.05
) (see Figure 29).Figure 26
The change ofl
3
(
α
1
,
α
2
,
α
3
) with α
1, α
2, and α
3.Figure 27
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.048.Figure 28
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.02.Figure 29
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.004.(2) Fix the value ofα
2 and α
1
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
2
=
0.004, α
1
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 30).
(ii)
α
2
=
0.004, α
3
=
0.01, and α
1
∈
(
0,0.05
) (see Figure 31).
(iii)
α
2
=
0.004, α
3
=
0.03, α
1
∈
(
0,0.05
) (see Figure 32).
(iv)
α
2
=
0.004, α
3
=
0.05, and α
1
∈
(
0,0.05
) (see Figure 33).Figure 30
The change ofl
3
(
α
1
,
α
2
,
α
3
) with α
1, α
2, and α
3.Figure 31
The change ofl
3
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.01.Figure 32
The change ofl
3
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.03.Figure 33
The change ofl
3
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.05.(3) Fix the value ofα
3 and α
1
∈
(
0,0.05
), α
2
∈
(
0,0.05
).(i)
α
3
=
0.03, α
1
∈
(
0,0.05
), and α
2
∈
(
0,0.05
) (see Figure 34).
(ii)
α
3
=
0.03, α
1
=
0.01, and α
2
∈
(
0,0.05
) (see Figure 35).
(iii)
α
3
=
0.03, α
1
=
0.02, and α
2
∈
(
0,0.05
) (see Figure 36).
(iv)
α
3
=
0.03, α
1
=
0.03, and α
2
∈
(
0,0.05
) (see Figure 37).Figure 34
The change ofl
3
(
α
1
,
α
2
,
α
3
) with α
1, α
2, and α
3.Figure 35
The change ofl
3
(
α
2
) with α
2 when α
1
=
0.01, α
3
=
0.03.Figure 36
The change ofl
3
(
α
2
) with α
2 when α
1
=
0.02, α
3
=
0.03.Figure 37
The change ofl
3
(
α
2
) with α
2 when α
1
=
0.03, α
3
=
0.03.Premisingα
3
=
0.01, the stable area in the α
1, α
2 plane is determined by inequality group (10). Figure 38 shows the system fixed point in the area of the asymptotic stability.Figure 38
α
3
=
0.01, the stable area in α
1, α
2 plane.


### 3.1.2. Numerical Simulation

Using Matlab, parameters influence on the system (8) can be analyzed through numerical simulation.(1)
α
1,
α
2,and
α
3
Influence on the Collecting Market. The second collector on the assumption that the manufacturer and the first collector stem variables parameters is fixed with α
1
=
0.01, α
3
=
0.01, her system variables parameters α
2 is at [
0,0.05
]. Manufacturers’ initial collecting price p
m, collector 1 and 2 for manufacturers of the collecting price of depreciate rate for 20, 0.35, respectively, 0.3. And, as shown in Figure 39 changes Figures 39(a), 39(b), and 39(c) shows, the corresponding Lyapunov index as shown in Figure 39(d). Figure 39 showing, when 0
<
α
2
<
0.0097, the system is in stable state. after multi game, p
m, r
1 and r
2 is stable at the point of (17.4, 0.4904, 0.4825). When α
1
=
0.01
,
α
2
=
0.0097
,
α
3
=
0.0097, the system (8) occurs the first bifurcation, then after cycle 2, 4 cycle, and so forth, the system is gradually into the chaotic state. Figure 39(d) is Lyapunov index spectrum distribution which can also confirmed the phenomena. at α
2
=
0.0099, the max Lyapunov index is zero, the system is in critical condition. At α
2
=
0.0133, most of Lyapunov index greater than zero which explains system at chaos state (Figure 39 and Figure 44).The max Lyapunov index.
(a)
α
1
=
0.01, α
3
=
0.01, and α
2
=
[
0.005
], p
m,
(b)
α
1
=
0.01, α
3
=
0.01
,  and α
2
=
[
0.005
], r
1
,
(c)
α
1
=
0.01, α
3
=
0.01, and α
2
=
[
0.005
], r
2,
(d)
α
1
=
0.01, α
3
=
0.01,   and α
2
=
[
0.005
]When initial setup and other parameters are the same andα
1, α
3, and α
3 are adjusted in the interval [
0,0.05
] simultaneously, system variables of the manufacturers and two collectors change as shown in Figures 40(a), 40(b), and 40(c), and Lyapunov index is as in Figure 40(d). When 0
<
α
2
<
0.0098, the system is in stable state at the point of (17.4, 0.4904, 0.4825). When the value of α
1, α
2, and α
3 is 0.0098, the system (8) occurs the first bifurcation, then after cycle 2, 4 cycle, and so forth, the system is gradually into the chaotic state. Figure 40(d) is Lyapunov index spectrum distribution which can also confirm the phenomena. At α
2
=
0.0099, the max Lyapunov index is zero, and the system is in critical condition. Comparing with Figure 39, system steady equilibrium does not change, keeps for (17.4, 0.4904, 0.4825).(a)α
1
,
α
2
,
α
3 is [
0
,
0.05
], p
1, (b) α
1
,
α
2
,
α
3 is [
0
,
0.05
], r
1, (c) α
1
,
α
2
,
α
2 is [
0
,
0.05
], r
2, and (d) α
1
,
α
2
,
α
3 is [
0
,
0.05
], the max Lyapunov index.
(a)
(b)
(c)
(d)(2) System Variables Power Spectrum. According to the numerical simulation, not only system variables graph in the phase space along with time can be drawn, but also power spectrum graph of system variables can be estimated by period chart method. When α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, power spectrum graph of the manufacturer and two collectors is shown as in Figures 41(a), 41(b), and 41(c).(a)α
1
=
0.01,  α
2
=
0.014, and α
3
=
0.01, power spectrum of p
m; (b) α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, power spectrum of r
1; and (c)α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, power of spectrum r
2.
(a)
(b)
(c)According to the numerical results, no matter how large the collection price adjustment parameters are, the collection price traverses the whole value area over time, butp
m, r
1, and r
2 in the system (10) always limit in a certain range, which also verifies boundedness and ergodicity of chaos. Because α
2
>
α
3 and r
1 amplitude is larger than r
2, the adjustment of system variables must be limited in a certain range.(3) The Influence of Initial Setup. The butterfly effect is an important symbol of chaotic motions, which is sensitive to initial state, and small changes of the initial conditions may cause the adjacent orbital evolution to index form separate after multiple episodes. For the above conditions, when α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, sensitive dependence on initial value of p
m, r
1, and r
2 is as shown in Figures 42(a), 42(b), and 42(c). p
m value is, respectively, taken 20 and 20.001, after 96 cycles iteration, p
m difference value is 0.1335 which is 133.5 times than the initial difference value. r
1 and r
2 value is, respectively, taken 0.35 and 0.351, 0.3 and 0.301, respectively, experience in 44 and 45 iteration, r
1 and r
2 difference values is respective 0.2826 and 0.1364 which are 282 times and 136 times than the initial difference values. Furthermore, from Figure 42 we can conclude that system variables values are approximate in the first period, with increase of iterations the difference increase obviously.(a)α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, the initial value sensitivity of p
m, (b) α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, the initial value sensitivity of r
1; and (c) α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, the initial value sensitivity r
2.
(a)
(b)
(c)


## 3.1.1. Model and Analysis

(
1) Model. The manufacturer and collectors are equal in closed-loop supply chain; that is, the two sides make decision at the same time: manufacturer’s decision is to choose the collection price p
i to maximize returns, and collectors’ decision is to choose p
r
i to maximize their returns.Nash equilibrium can be concluded by the first-order conditions of three reaction functions, which meet∂
Π
m
/
∂
p
m
=
0, ∂
Π
r
i
/
∂
r
i
=
0 simultaneously. The Nash equilibrium is the optimal decision of a game in all kinds of may guess, thus make it to gain maximum benefit. The Nash equilibrium is
(5)
(
p
m
*
,
r
1
*
,
r
2
*
)
=
(
-
B
+
B
2
-
4
A
C
2
A
,
t
1
+
t
2
p
m
*
,
t
3
+
t
4
p
m
*
)
,
where A, B, and C are
(6)
A
=
2
(
β
1
+
β
2
-
δ
1
-
δ
2
)
(
p
0
-
c
r
)
h
h
h
-
2
t
1
(
β
1
-
δ
2
)
-
2
t
3
(
β
2
-
δ
1
)
,
B
=
k
1
+
k
2
-
(
β
1
+
β
2
-
δ
1
-
δ
2
)
(
p
0
-
c
r
)
h
h
h
+
(
p
0
-
c
r
)
(
t
1
(
β
1
-
δ
2
)
+
t
3
(
β
2
-
δ
1
)
)
h
h
h
-
2
t
2
(
β
1
-
δ
2
)
-
2
t
4
(
β
2
-
δ
1
)
,
C
=
(
p
0
-
c
r
)
t
2
(
β
1
-
δ
2
)
+
t
4
(
β
2
-
δ
1
)
.In formula (6), t
1, t
2, t
3, t
4 can be expressed as
(7)
t
1
=
2
β
1
β
2
-
2
β
2
δ
1
+
δ
1
β
2
-
δ
1
δ
2
4
β
1
β
2
-
δ
1
δ
2
,
t
2
=
2
β
2
k
1
+
2
β
2
β
1
A
1
+
δ
1
k
2
-
δ
1
β
2
A
2
4
β
1
β
2
-
δ
1
δ
2
,
t
3
=
2
β
1
β
2
-
2
β
1
δ
2
+
δ
2
β
1
-
δ
1
δ
2
4
β
1
β
2
-
δ
1
δ
2
,
t
4
=
2
β
1
k
2
+
2
β
2
β
1
A
2
+
δ
2
k
1
-
δ
2
β
1
A
1
4
β
1
β
2
-
δ
1
δ
2
.In reality, the game between node enterprises in closed-loop supply chain is continuous, enterprises’ decision-making is a long-term repeated process, and its action has long-term memory. And each node enterprise does not completely control the market information and also cannot fully expect future market changes, so based on limited rational expectations decision we adjust process with marginal gains. They can make the next-period price decision on the basis of the local estimate to his marginal profits in current period. Their price adjustment processes are(8)
p
m
(
t
+
1
)
=
p
m
(
t
)
+
α
1
p
m
(
t
)
×
(
(
(
1
-
r
1
(
t
)
)
(
β
1
-
δ
2
)
+
(
1
-
r
2
(
t
)
)
(
β
2
-
δ
1
)
)
×
(
p
0
-
2
p
m
(
t
)
-
c
r
)
-
k
1
-
k
2
(
(
1
-
r
1
(
t
)
)
(
β
1
-
δ
2
)
+
(
1
-
r
2
(
t
)
)
(
β
2
-
δ
1
)
)
)
,
r
1
(
t
+
1
)
=
r
1
(
t
)
+
α
2
r
1
(
t
)
×
(
-
2
β
1
r
1
(
t
)
p
m
2
(
t
)
+
p
m
(
t
)
h
h
h
×
(
k
1
+
(
β
1
-
δ
1
+
δ
1
r
2
(
t
)
)
h
h
h
×
p
m
(
t
)
+
β
1
A
1
)
(
-
2
β
1
r
1
(
t
)
p
m
2
(
t
)
+
p
m
(
t
)
)
,
r
2
(
t
+
1
)
=
r
2
(
t
)
+
α
3
r
2
(
t
)
×
(
-
2
β
2
r
2
(
t
)
p
m
2
(
t
)
+
p
m
(
t
)
h
h
×
(
k
2
+
(
β
2
-
δ
2
+
δ
2
r
1
(
t
)
)
p
m
(
t
)
+
β
2
A
2
)
(
-
2
β
2
r
2
(
t
)
p
m
2
(
t
)
+
p
m
(
t
)
)
,
where p
1
(
t
+
1
),  r
1
(
t
+
1
), r
2
(
t
+
1
), and p
1
(
t
), r
1
(
t
), r
2
(
t
) are used to present the price in the (
t
+
1
)th period and in the (
t
)th period. Where α
i are positive parameters, i
=
1,2
,
3, which denotes adjustment speed, respectively, for the manufacturer and two collectors.From the system (8) we can be conclude that the optimal collecting price of the manufacturer is related to the collection price adjustment coefficient, the two collectors’ collecting price, sale price, and consumers’ environmental protection awareness. Similarly, the optimal collecting price of collectors is related to collecting price adjustment parameters, collecting price of the manufacturer, collecting function, consumers’s environmental protection awareness, and collecting costs.(2) Model Analysis. In closed-loop supply chain, any enterprise decision-making is all according to the maximum profits, the equilibrium less than zero is not practical significance, such as manufacturers to make power mainly because of benefits, collecting business is also benefit with collecting. So only research system (8) is the equilibrium; system (8) of the eight fixed points, respectively, is E
0
(
0,0
,
0
),  E
1
(
p
m
1
,
0,0
),  E
2
(
0
,
r
1
2
,
0
),  E
3
(
0
,
0
,
r
2
3
),  E
4
(
p
m
4
,
r
1
4
,
0
), E
5
(
p
m
5
,
0
,
r
2
5
), E
6
(
0
,
r
1
6
,
r
2
6
), E
7
(
p
m
*
,
r
1
*
,
r
2
*
). The analysis shows that by E
0, E
1, E
2, E
3, E
4, E
5 and E
6, is no real significance. Therefore, we discuss Nash equilibrium which is E
7
(
p
m
*
,
r
1
*
,
r
2
*
); E
7 is the reaction function intersection of the manufacturer and collectors, which means marginal profits of both sides are zero, but this does not mean that the result of the game will tend to equilibrium. Instead, a party rational behavior change may cause game process to occur very complex phenomenon. E
7 has local stability, this stability region by α
1, α
2, α
3 decided. The stability of E
7 will be studied. The first we calculate Jacobi matrix of system (8), means
(9)
J
=
(
J
11
J
12
J
13
J
21
J
22
J
23
J
31
J
32
J
33
)
.We putE
7
(
p
m
*
,
r
1
*
,
r
2
*
) into formula (9), according to the Jury stability criterion, the sufficient and necessary conditions of E
7 asymptotic stability are that the zeros of characteristic polynomial F
(
λ
)
=
λ
3
+
A
λ
2
+
B
λ
+
C are in a circular area. Moreover, the characteristic polynomial will satisfy the conditions as follows:
(10)
F
(
1
)
=
1
+
A
+
B
+
C
>
0
,
(
-
1
)
3
F
(
-
1
)
=
1
-
A
+
B
-
C
>
0
,
1
>
|
C
|
,
1
-
C
2
>
|
B
-
A
C
|
.At present, only a few simple dynamic systems are analyzed with analytical method, but complex dynamic system mainly uses the numerical analysis method. This paper processes numerical simulation on system (8) by Matlab to expresses the dynamic characteristics. Because remanufacturing is still initial stage in china, consumers environmental protection awareness is lower, and manufacturing cost (including the fixed and variable cost) is higher; for autoparts collecting, the parameters of system (8) can be defined as p
0
=
70, c
r
=
30, A
1
=
1, A
2
=
0.8, k
1
=
1, k
2
=
1, β
1
=
0.6, β
2
=
0.55, δ
1
=
0.2, and δ
2
=
0.2 to study the local stability of equilibrium point. Get the parameter value into E
7
(
p
m
*
,
r
1
*
,
r
2
*
) and formula (10) to conclude the value for (17.385, 0.491, 0.489).Change the discrete system (8) into a continuous system and the points of Jacobi matrix for
(11)
J
=
[
-
13.317
α
1
-
17.755
α
1
-
15.562
α
1
-
0.775
α
2
-
177.888
α
2
29.68
α
2
-
0.821
α
3
29.599
α
3
-
164.65
α
3
]
.Consider the continuous-time nonlinear dynamical system(12)
x
˙
=
A
x
+
N
(
x
)
,
A
=
(
a
i
j
)
n
×
n
,
x
∈
R
n
.Let the functionN
(
x
) be written as
(13)
N
(
x
)
=
1
2
B
(
x
,
x
)
+
1
6
C
(
x
,
x
,
x
)
+
O
(
∥
x
∥
4
)
,
where B
(
x
,
y
) and C
(
x
,
y
,
z
) are bilinear and trilinear functions, respectively. In coordinates, we have
(14)
B
i
(
x
,
y
)
=
∑
j
,
k
=
1
n
∂
2
N
i
(
ξ
)
∂
ξ
j
∂
ξ
k
|
ξ
=
0
x
j
y
k
,
C
i
(
x
,
y
,
z
)
=
∑
j
,
k
,
l
=
1
n
∂
3
N
i
(
ξ
)
∂
ξ
j
∂
ξ
k
∂
ξ
l
|
ξ
=
0
x
j
y
k
z
l
.Supposing thatA has a pair of complex eigenvalues on the imaginary axis, they are λ
1,2
=
±
i
ω
(
ω
>
0
), and these eigenvalues are the only eigenvalues with Re
λ
=
0. Let q
∈
C
n be a complex eigenvector to λ
1
=
i
ω:
(15)
A
q
=
i
ω
q
,
A
q
¯
=
-
i
ω
q
¯
.And the adjoint eigen vectorp
∈
C
n admits the properties
(16)
A
T
p
=
-
i
ω
p
,
A
T
p
¯
=
i
ω
p
¯and satisfies the normalization〈
q
,
p
〉
=
1.The first Lyapunov coefficient at the origin is defined by(17)
l
1
=
1
2
ω
Re
[
〈
p
,
C
(
q
,
q
,
q
¯
)
〉
-
2
〈
p
,
B
(
q
,
A
-
1
B
(
q
,
q
¯
)
)
〉
hhhhhhhh
+
〈
p
,
B
(
q
¯
,
(
2
i
ω
E
-
A
)
-
1
B
(
q
,
q
)
)
〉
]
.Next, we calculatei
w
q, J
q
-
=
-
i
w
q
-, J
T
p
=
-
i
w
p, and  J
T
p
-
=
i
w
p
-; If q
=
(
q
1
,
q
2
,
q
3
)
T
≠
0, calculate q
,
q
-
,
p
,
p
-:
(18)
[
-
13.317
α
1
-
i
w
-
17.755
α
1
-
15.562
α
1
-
0.775
α
2
-
177.888
α
2
-
i
w
29.68
α
2
-
0.821
α
3
29.599
α
3
-
164.65
α
3
-
i
w
]
×
[
q
1
q
2
q
3
]
=
[
0
0
0
]
having nonzero solutions.It means(19)
|
-
13.317
α
1
-
i
w
-
17.755
α
1
-
15.562
α
1
-
0.775
α
2
-
177.888
α
2
-
i
w
29.68
α
2
-
0.821
α
3
29.599
α
3
-
164.65
α
3
-
i
w
|
=
0
.It can be got from calculating(20)
f
1
(
w
)
=
-
375716.354
α
1
α
2
α
3
+
117.77755
α
2
2
α
3
+
13.317
α
1
w
2
+
177.888
α
2
w
2
+
164.65
α
3
w
2
=
0
,
f
2
(
w
)
=
-
2368.9345
α
1
α
2
w
+
0.6
α
2
2
w
-
28410.761
α
2
α
3
w
-
2179.87
α
1
α
3
w
+
w
3
=
0
.(1) When{
w
=
α
1
f
1
(
w
)
=
0
f
2
(
w
)
=
0, matrix J has pure imaginary eigenvalues ±
i
α
1. Calculate q, q
-, p, p
-.The result is as follows.(21)
L
2
(
α
1
,
α
2
,
α
3
)
=
1
2
ω
[
〈
p
,
C
(
q
,
q
,
q
-
)
〉
-
2
〈
p
,
B
(
q
,
A
-
1
B
(
q
,
q
-
)
)
〉
h
h
h
h
h
h
+
〈
p
,
B
(
q
-
,
(
2
i
ω
E
-
A
)
-
1
B
(
q
,
q
)
)
〉
]
=
(
408.75
-
29.6
i
)
α
1
α
2
(
0.82
i
α
1
-
168.9853
α
2
)
α
2
-
2
[
(
408.75
-
29.6
i
)
α
1
α
2
(
0.82
i
α
1
-
168.9853
α
2
)
α
2
-
2179.87
α
1
375716.354
α
1
α
2
-
117.7775
α
2
2
+
(
+
α
1
α
2
(
1.2
i
α
2
+
(
375716.354
+
56821.5
i
)
α
3
)
}
)
-
1
(
61899.185
-
4510.5
i
)
α
2
×
(
5846897.898
i
α
1
2
-
(
1.24
×
1
0
9
+
1832.854
i
)
α
1
α
2
+
388107.8
α
2
2
)
-
1
(
(
61899.185
-
4510.5
i
)
α
2
)
(
408.75
-
29.6
i
)
α
1
α
2
(
0.82
i
α
1
-
168.9853
α
2
)
α
2
]
+
(
+
α
1
α
2
(
1.2
i
α
2
+
(
375716.354
+
56821.5
i
)
α
3
)
}
)
-
1
{
α
1
[
(
-
4
+
26.634
i
)
α
1
h
h
h
h
h
h
h
h
+
(
2179.87
+
329.3
i
)
α
3
]
}
×
(
{
(
-
53
-
8
i
)
α
1
3
+
α
1
2
[
(
-
711.55
+
4737.8
i
)
α
2
-
(
658.6
-
4359.7
i
)
α
3
]
-
117.8
α
2
2
α
3
+
α
1
α
2
(
1.2
i
α
2
+
(
375716.354
+
56821.5
i
)
α
3
)
(
-
53
-
8
i
)
α
1
3
}
)
-
1
)
+
{
(
-
(
407.31
+
29.7
i
)
α
2
)
×
(
-
1.55
i
α
1
α
2
-
151.97
α
2
α
3
)
}
×
(
{
(
15.6
i
α
1
+
3295.3
α
2
)
(
-
53
-
8
i
)
α
1
3
×
{
(
-
53.3
-
8
i
)
α
1
3
(
-
53
-
8
i
)
α
1
3
+
α
1
2
[
(
-
711.552
+
4737.9
i
)
α
2
-
(
658.6
-
4359.74
i
)
α
3
]
-
117.77755
α
2
2
α
3
+
α
1
α
2
×
[
1.2
i
α
2
+
(
375716.354
+
56821.52
i
)
α
3
]
(
-
53
-
8
i
)
α
1
3
}
}
)
-
1
=
l
1
(
α
1
,
α
2
,
α
3
)
+
m
1
(
α
1
,
α
2
,
α
3
)
,
where L
1
(
α
1
,
α
2
,
α
3
) is the Lyapunov exponent function about α
1, α
2, and α
3 and l
1
(
α
1
,
α
2
,
α
3
) is the real part. The following is the figure of l
1
(
α
1
,
α
2
,
α
3
).(1) Fix the value ofα
1 and α
2
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
1
=
0.02, α
2
=
0.048, and α
3
∈
(
0,0.05
) (see Figure 3).
(ii)
α
1
=
0.02, α
2
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 2).Figure 2
The change ofl
1
(
α
1
,
α
2
,
α
3
) with α
1
,
α
2, and α
3.Figure 3
The change ofl
1
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.048.We can know from Figure3 that when α
1 =  0.02, α
2 = 0.048, and α
3 is between zero and 0.05, l
1
(
α
3
) is always equal to 1.375.(iii)
α
1
=
0.02, α
2
=
0.02, and α
3
∈
(
0,0.05
) (see Figure 4).
(iv)
α
1
=
0.02, α
2
=
0.004, and α
3
∈
(
0,0.05
) (see Figure 5).Figure 4
The change ofl
1
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.02.Figure 5
The change ofl
1
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.004.From Figures4 and 5 we can observe the same change of l
1
(
α
3
) as in Figures 2 and 3.(2) Fix the value ofα
2 and α
1
∈
(
0,0.05
), α
3
∈
(
0,0.05
). α
2
=
0.004, α
1
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 6). α
2
=
0.004, α
3
=
0.01, and α
1
∈
(
0,0.05
) (see Figure 7). α
2
=
0.004, α
3
=
0.03, and α
1
∈
(
0,0.05
) (see Figure 8). α
2
=
0.004, α
3
=
0.05, and α
1
∈
(
0,0.05
) (see Figure 9).Figure 6
The change ofl
1
(
α
1
,
α
2
,
α
3
) with α
2, α
1, and α
3.Figure 7
The change ofl
1
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.01.Figure 8
The change ofl
1
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.03.Figure 9
The change ofl
1
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.05.(3) Fix the value ofα
3 and α
1
∈
(
0,0.05
), α
2
∈
(
0,0.05
).(i)
α
3
=
0.03, α
1
∈
(
0,0.05
), and α
2
∈
(
0,0.05
) (see Figure 10).
(ii)
α
3
=
0.03, α
1
=
0.01, and α
2
∈
(
0,0.05
) (see Figure 11).
(iii)
α
3
=
0.03, α
1
=
0.02, and α
2
∈
(
0,0.05
) (see Figure 12).
(iv)
α
3
=
0.03, α
1
=
0.03, and α
2
∈
(
0,0.05
) (see Figure 13).Figure 10
The change ofl
1
(
α
1
,
α
2
,
α
3
) with α
3, α
1, and α
2.Figure 11
The change ofl
1
(
α
2
) with α
2 when α
1
=
0.01, α
3
=
0.03.Figure 12
The change ofl
1
(
α
2
) with α
2 when α
1
=
0.02, α
3
=
0.03.Figure 13
The change ofl
1
(
α
2
) with α
2 when α
1
=
0.03, α
3
=
0.03.(2) When{
w
=
α
2
f
1
(
α
2
)
=
0
f
2
(
α
2
)
=
0, matrix J has pure imaginary eigenvalues ±
i
α
2. Calculate q, q
-, p, p
-.The result is as follows:(22)
L
2
(
α
1
,
α
2
,
α
3
)
=
1
2
ω
[
〈
p
,
C
(
q
,
q
,
q
-
)
〉
-
2
〈
p
,
B
(
q
,
A
-
1
B
(
q
,
q
-
)
)
〉
+
〈
p
,
B
(
q
-
,
(
2
i
ω
E
-
A
)
-
1
B
(
q
,
q
)
)
〉
]
=
(
0.000851
-
0.175
i
)
+
37.57
+
0.18
i
375716.354
α
1
-
117.8
α
2
+
(
2.42
+
0.012
i
)
α
1
α
2
+
(
0.01293
-
2.74
i
)
α
2
375716.354
α
1
2
-
117.78
α
1
α
2
+
4359.74
α
1
375716.354
α
1
α
2
-
117.77755
α
2
2
+
(
{
α
2
2
[
(
-
0.014
+
0.000066
i
)
α
2
+
(
0.0065
+
1.37
i
)
α
3
]
+
α
1
α
2
[
(
-
4
+
0.19
i
)
α
2
+
(
18.784
+
329.2
i
)
α
3
]
+
α
1
2
(
26.6
i
α
2
+
2179.87
α
3
)
}
)
×
{
α
1
α
2
[
(
-
53.268
+
4737.87
i
)
α
1
α
2
-
(
711.552
+
9.2
i
)
α
2
2
+
(
375716.354
+
4359.74
i
)
α
1
α
3
-
(
776.38
-
56821.52
i
)
α
2
α
3
]
}
-
1
=
l
2
(
α
1
,
α
2
,
α
3
)
+
m
2
(
α
1
,
α
2
,
α
3
)
,
where L
2
(
α
1
,
α
2
,
α
3
) is the Lyapunov exponent function about α
1, α
2, and α
3 and l
2
(
α
1
,
α
2
,
α
3
) is the real part. The following is the figure of l
2
(
α
1
,
α
2
,
α
3
).(1) Fix the value ofα
1 and α
2
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
1
=
0.02, α
2
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 14).
(ii)
α
1
=
0.02, α
2
=
0.048, and α
3
∈
(
0,0.05
) (see Figure 15).
(iii)
α
1
=
0.02, α
2
=
0.02, and α
3
∈
(
0,0.05
) (see Figure 16).
(iv)
α
1
=
0.02, α
2
=
0.004, and α
3
∈
(
0,0.05
) (see Figure 17).Figure 14
The change ofl
2
(
α
1
,
α
2
,
α
3
) with α
3, α
1, and α
2.Figure 15
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.048.Figure 16
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.02.Figure 17
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.004.(2) Fix the value ofα
2 and α
1
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
2
=
0.004, α
1
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 18).
(ii)
α
2
=
0.004, α
3
=
0.01, and α
1
∈
(
0,0.05
) (see Figure 19).
(iii)
α
2
=
0.004, α
1
∈
(
0,0.05
), and α
3
=
0.03 (see Figure 20).
(iv)
α
2
=
0.004, α
3
=
0.05, and α
1
∈
(
0,0.05
) (see Figure 21).Figure 18
The change ofl
2
(
α
1
,
α
2
,
α
3
) with α
2, α
1, and α
3.Figure 19
The change ofl
2
(
α
1
) with α
1 when α
2
=
0.004
,
α
2
=
0.01.Figure 20
The change ofl
2
(
α
1
) with α
1 when α
2
=
0.004, α
2
=
0.03.Figure 21
The change ofl
2
(
α
1
) with α
1 when α
2
=
0.004, α
2
=
0.05.(3) Fix the value ofα
3 and α
1
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
3
=
0.03, α
1
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 22).
(ii)
α
3
=
0.03, α
1
=
0.01, and α
2
∈
(
0,0.05
) (see Figure 23).
(iii)
α
3
=
0.03, α
1
=
0.02, and α
2
∈
(
0,0.05
) (see Figure 24).
(iv)
α
3
=
0.03, α
1
=
0.03, and α
2
∈
(
0,0.05
) (see Figure 25).Figure 22
The change ofl
2
(
α
1
,
α
2
,
α
3
) with α
2, α
1, and α
3.Figure 23
The change ofl
2
(
α
2
) with α
2 when α
1
=
0.01, α
3
=
0.03.Figure 24
The change ofl
2
(
α
2
) with α
2 when α
1
=
0.02, α
3
=
0.03.Figure 25
The change ofl
2
(
α
2
) with α
2 when α
1
=
0.03, α
3
=
0.03.(
3
) When {
w
=
α
3
f
1
(
α
3
)
=
0
f
2
(
α
3
)
=
0, matrix J has pure imaginary eigenvalues ±
i
α
3. Calculate q, q
-, p, p
-.The result is as follows:(23)
L
3
(
α
1
,
α
2
,
α
3
)
=
1
2
ω
[
〈
p
,
C
(
q
,
q
,
q
-
)
〉
-
2
〈
p
,
B
(
q
,
A
-
1
B
(
q
,
q
-
)
)
〉
+
〈
p
,
B
(
q
-
,
(
2
i
ω
E
-
A
)
-
1
B
(
q
,
q
)
)
〉
]
=
408.75
α
1
α
2
-
29.6
i
α
2
α
3
α
2
(
-
168.9853
α
2
+
0.821
i
α
3
)
+
{
(
7183233.84
+
87766
i
)
α
1
2
α
2
+
(
61899.185
+
631.33
i
)
α
1
α
2
2
-
(
414.48
-
33923.1
i
)
α
1
2
α
3
-
(
13181
-
1085129.6
i
)
α
1
α
2
α
3
-
(
46
-
4510.5
i
)
α
2
2
α
3
-
(
5124.6
+
62.25
i
)
α
1
α
3
2
}
×
{
-
(
658.6
+
8
i
)
α
3
2
]
(
3295.26
α
1
α
2
+
15.562
i
α
1
α
3
)
×
[
-
(
658.6
+
8
i
)
α
3
2
(
375716.354
+
4737.87
i
)
α
1
α
2
-
(
117.78
+
1.2
i
)
α
2
2
-
(
53.3
-
4359.74
i
)
α
1
α
3
-
(
711.552
-
56821.5
i
)
α
2
α
3
-
(
658.6
+
8
i
)
α
3
2
]
}
-
1
-
2
[
-
2179.867
α
1
375716.354
α
1
α
2
-
117.78
α
2
2
+
408.75
α
1
α
2
-
29.6
i
α
2
α
3
α
2
(
-
168.9853
α
2
+
0.821
i
α
3
)
+
(
61899.185
α
1
α
2
-
4510.5
i
α
2
α
3
)
×
(
(
-
375716.354
α
1
+
117.78
α
2
)
×
(
3295.26
α
1
α
2
-
15.562
i
α
1
α
3
)
)
-
1
-
2179.867
α
1
375716.354
α
1
a
2
-
117.78
α
2
2
]
=
l
3
(
α
1
,
α
2
,
α
3
)
+
m
3
(
α
1
,
α
2
,
α
3
)
,
where L
3
(
α
1
,
α
2
,
α
3
) is the Lyapunov exponent function about α
1, α
2, and α
3 and l
3
(
α
1
,
α
2
,
α
3
) is the real part. The following is the figure of l
3
(
α
1
,
α
2
,
α
3
).(1) Fix the value ofα
1 and α
2
∈
(
0,0.05
)
,
α
3
∈
(
0,0.05
).(i)
α
1
=
0.02, α
2
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 26).
(ii)
α
1
=
0.02, α
2
=
0.048, and α
3
∈
(
0,0.05
) (see Figure 27).
(iii)
α
1
=
0.02, α
2
=
0.02, and α
3
∈
(
0,0.05
) (see Figure 28).
(iv)
α
1
=
0.02, α
2
=
0.004, and α
3
∈
(
0,0.05
) (see Figure 29).Figure 26
The change ofl
3
(
α
1
,
α
2
,
α
3
) with α
1, α
2, and α
3.Figure 27
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.048.Figure 28
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.02.Figure 29
The change ofl
2
(
α
3
) with α
3 when α
1
=
0.02, α
2
=
0.004.(2) Fix the value ofα
2 and α
1
∈
(
0,0.05
), α
3
∈
(
0,0.05
).(i)
α
2
=
0.004, α
1
∈
(
0,0.05
), and α
3
∈
(
0,0.05
) (see Figure 30).
(ii)
α
2
=
0.004, α
3
=
0.01, and α
1
∈
(
0,0.05
) (see Figure 31).
(iii)
α
2
=
0.004, α
3
=
0.03, α
1
∈
(
0,0.05
) (see Figure 32).
(iv)
α
2
=
0.004, α
3
=
0.05, and α
1
∈
(
0,0.05
) (see Figure 33).Figure 30
The change ofl
3
(
α
1
,
α
2
,
α
3
) with α
1, α
2, and α
3.Figure 31
The change ofl
3
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.01.Figure 32
The change ofl
3
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.03.Figure 33
The change ofl
3
(
α
1
) with α
1 when α
2
=
0.004, α
3
=
0.05.(3) Fix the value ofα
3 and α
1
∈
(
0,0.05
), α
2
∈
(
0,0.05
).(i)
α
3
=
0.03, α
1
∈
(
0,0.05
), and α
2
∈
(
0,0.05
) (see Figure 34).
(ii)
α
3
=
0.03, α
1
=
0.01, and α
2
∈
(
0,0.05
) (see Figure 35).
(iii)
α
3
=
0.03, α
1
=
0.02, and α
2
∈
(
0,0.05
) (see Figure 36).
(iv)
α
3
=
0.03, α
1
=
0.03, and α
2
∈
(
0,0.05
) (see Figure 37).Figure 34
The change ofl
3
(
α
1
,
α
2
,
α
3
) with α
1, α
2, and α
3.Figure 35
The change ofl
3
(
α
2
) with α
2 when α
1
=
0.01, α
3
=
0.03.Figure 36
The change ofl
3
(
α
2
) with α
2 when α
1
=
0.02, α
3
=
0.03.Figure 37
The change ofl
3
(
α
2
) with α
2 when α
1
=
0.03, α
3
=
0.03.Premisingα
3
=
0.01, the stable area in the α
1, α
2 plane is determined by inequality group (10). Figure 38 shows the system fixed point in the area of the asymptotic stability.Figure 38
α
3
=
0.01, the stable area in α
1, α
2 plane.


## 3.1.2. Numerical Simulation

Using Matlab, parameters influence on the system (8) can be analyzed through numerical simulation.(1)
α
1,
α
2,and
α
3
Influence on the Collecting Market. The second collector on the assumption that the manufacturer and the first collector stem variables parameters is fixed with α
1
=
0.01, α
3
=
0.01, her system variables parameters α
2 is at [
0,0.05
]. Manufacturers’ initial collecting price p
m, collector 1 and 2 for manufacturers of the collecting price of depreciate rate for 20, 0.35, respectively, 0.3. And, as shown in Figure 39 changes Figures 39(a), 39(b), and 39(c) shows, the corresponding Lyapunov index as shown in Figure 39(d). Figure 39 showing, when 0
<
α
2
<
0.0097, the system is in stable state. after multi game, p
m, r
1 and r
2 is stable at the point of (17.4, 0.4904, 0.4825). When α
1
=
0.01
,
α
2
=
0.0097
,
α
3
=
0.0097, the system (8) occurs the first bifurcation, then after cycle 2, 4 cycle, and so forth, the system is gradually into the chaotic state. Figure 39(d) is Lyapunov index spectrum distribution which can also confirmed the phenomena. at α
2
=
0.0099, the max Lyapunov index is zero, the system is in critical condition. At α
2
=
0.0133, most of Lyapunov index greater than zero which explains system at chaos state (Figure 39 and Figure 44).The max Lyapunov index.
(a)
α
1
=
0.01, α
3
=
0.01, and α
2
=
[
0.005
], p
m,
(b)
α
1
=
0.01, α
3
=
0.01
,  and α
2
=
[
0.005
], r
1
,
(c)
α
1
=
0.01, α
3
=
0.01, and α
2
=
[
0.005
], r
2,
(d)
α
1
=
0.01, α
3
=
0.01,   and α
2
=
[
0.005
]When initial setup and other parameters are the same andα
1, α
3, and α
3 are adjusted in the interval [
0,0.05
] simultaneously, system variables of the manufacturers and two collectors change as shown in Figures 40(a), 40(b), and 40(c), and Lyapunov index is as in Figure 40(d). When 0
<
α
2
<
0.0098, the system is in stable state at the point of (17.4, 0.4904, 0.4825). When the value of α
1, α
2, and α
3 is 0.0098, the system (8) occurs the first bifurcation, then after cycle 2, 4 cycle, and so forth, the system is gradually into the chaotic state. Figure 40(d) is Lyapunov index spectrum distribution which can also confirm the phenomena. At α
2
=
0.0099, the max Lyapunov index is zero, and the system is in critical condition. Comparing with Figure 39, system steady equilibrium does not change, keeps for (17.4, 0.4904, 0.4825).(a)α
1
,
α
2
,
α
3 is [
0
,
0.05
], p
1, (b) α
1
,
α
2
,
α
3 is [
0
,
0.05
], r
1, (c) α
1
,
α
2
,
α
2 is [
0
,
0.05
], r
2, and (d) α
1
,
α
2
,
α
3 is [
0
,
0.05
], the max Lyapunov index.
(a)
(b)
(c)
(d)(2) System Variables Power Spectrum. According to the numerical simulation, not only system variables graph in the phase space along with time can be drawn, but also power spectrum graph of system variables can be estimated by period chart method. When α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, power spectrum graph of the manufacturer and two collectors is shown as in Figures 41(a), 41(b), and 41(c).(a)α
1
=
0.01,  α
2
=
0.014, and α
3
=
0.01, power spectrum of p
m; (b) α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, power spectrum of r
1; and (c)α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, power of spectrum r
2.
(a)
(b)
(c)According to the numerical results, no matter how large the collection price adjustment parameters are, the collection price traverses the whole value area over time, butp
m, r
1, and r
2 in the system (10) always limit in a certain range, which also verifies boundedness and ergodicity of chaos. Because α
2
>
α
3 and r
1 amplitude is larger than r
2, the adjustment of system variables must be limited in a certain range.(3) The Influence of Initial Setup. The butterfly effect is an important symbol of chaotic motions, which is sensitive to initial state, and small changes of the initial conditions may cause the adjacent orbital evolution to index form separate after multiple episodes. For the above conditions, when α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, sensitive dependence on initial value of p
m, r
1, and r
2 is as shown in Figures 42(a), 42(b), and 42(c). p
m value is, respectively, taken 20 and 20.001, after 96 cycles iteration, p
m difference value is 0.1335 which is 133.5 times than the initial difference value. r
1 and r
2 value is, respectively, taken 0.35 and 0.351, 0.3 and 0.301, respectively, experience in 44 and 45 iteration, r
1 and r
2 difference values is respective 0.2826 and 0.1364 which are 282 times and 136 times than the initial difference values. Furthermore, from Figure 42 we can conclude that system variables values are approximate in the first period, with increase of iterations the difference increase obviously.(a)α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, the initial value sensitivity of p
m, (b) α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, the initial value sensitivity of r
1; and (c) α
1
=
0.01, α
2
=
0.014, and α
3
=
0.01, the initial value sensitivity r
2.
(a)
(b)
(c)


## 3.2. Stackelberg Equilibrium

### 3.2.1. Model and the Analysis

Suppose that the manufacturer and collectors have principal and subordinate relationship, the manufacturer is Stackelberg leader, and collectors are followers; then they process sequential dynamic game; the game equilibrium is Stackelberg equilibrium. In this game, the manufacturer makes the decision of sales price and collection price according to the market informa Lyapunov exponent function abouttion; then two collectors make decision according to the decision-making of the manufacture:(24)
max
⁡
∏
r
1
=
(
p
m
-
p
r
1
(
t
)
-
A
1
)
(
k
1
+
β
1
p
r
1
(
t
)
-
δ
1
p
r
2
(
t
)
)
max
⁡
∏
r
2
=
(
p
m
-
p
r
2
(
t
)
-
A
2
)
(
k
2
+
β
2
p
r
2
(
t
)
-
δ
2
p
r
1
(
t
)
)
s
.
t
.
i
p
m
∈
arg
max
⁡
∏
r
1
s
.
t
.
i
p
m
∈
arg
max
⁡
∏
r
2
.The result is(25)
r
1
*
*
=
(
p
m
(
2
β
1
β
2
-
2
β
2
δ
1
+
δ
1
β
2
-
δ
1
δ
2
)
+
2
β
2
k
1
+
2
β
1
β
2
A
1
+
k
2
δ
1
+
β
2
A
2
δ
1
)
×
(
p
m
(
4
β
1
β
2
-
δ
1
δ
2
)
)
-
1
,
r
2
*
*
=
(
+
2
β
1
k
2
+
2
β
1
β
2
A
2
+
k
1
δ
2
+
β
1
A
2
δ
2
)
p
m
(
2
β
1
β
2
-
2
β
1
δ
2
+
δ
2
β
1
-
δ
1
δ
2
)
+
2
β
1
k
2
+
2
β
1
β
2
A
2
+
k
1
δ
2
+
β
1
A
2
δ
2
)
×
(
p
m
(
4
β
1
β
2
-
δ
1
δ
2
)
)
-
1
.Formula (25) is reaction function of the collectors. Put formula (25) into formula (2) and the following can be obtained:
(26)
p
m
*
*
=
(
p
0
-
c
r
)
2
-
(
(
β
1
-
δ
2
)
(
2
β
2
k
1
+
2
β
2
β
1
A
1
+
δ
1
k
2
+
δ
1
β
2
A
2
)
+
(
β
2
-
δ
1
)
(
2
β
1
k
2
+
2
β
2
β
1
A
2
+
δ
2
k
1
-
δ
2
β
1
A
1
)
-
(
4
β
1
β
2
-
δ
1
δ
2
)
(
k
1
+
k
2
)
)
×
(
(
2
β
1
β
2
+
2
β
2
δ
1
-
δ
1
β
2
)
(
β
2
-
δ
1
)
+
(
2
β
1
β
2
+
2
β
1
δ
2
-
δ
2
β
1
)
(
β
1
-
δ
2
)
)
-
1
.With Nash equilibrium value the parameters values arep
0
=
70, c
r
=
30, A
1
=
1, A
2
=
0.8, k
1
=
1, k
2
=
1, β
1
=
0.6, β
2
=
0.55, δ
1
=
0.2, δ
2
=
0.2. Then we put the values into Stackelberg equilibrium (
p
m
*
*
,
r
1
*
*
,
r
2
*
*
); the equilibrium value is (
19.17
,
0.48
,
0.48
).The manufacturer makes decision based on limited rational expectations. She adjusts the game process on the basis of marginal gains. If the marginal profits oft period are positive, in the t
+
1 period this strategy action will be used. And if the t period marginal profits are negative, in the t
+
1 period the collection price will be lowered. The repeated game model of dynamic adjustment is
(27)
p
m
(
t
+
1
)
=
p
m
(
t
)
+
α
p
m
(
t
)
×
(
(
(
1
-
r
1
)
(
β
1
-
δ
2
)
+
(
1
-
r
2
)
(
β
2
-
δ
1
)
)
×
(
p
0
-
2
p
m
(
t
)
-
c
r
)
-
k
1
-
k
2
(
(
1
-
r
1
)
(
β
1
-
δ
2
)
+
(
1
-
r
2
)
(
β
2
-
δ
1
)
)
)
.In formula (27), α is adjustment coefficient of collection price and r
1
,
r
2 are as follows:
(28)
r
1
(
t
)
=
(
2
β
2
k
1
+
2
β
1
β
2
A
1
+
k
2
δ
1
+
β
2
A
2
δ
1
p
m
(
t
)
(
2
β
1
β
2
-
2
β
2
δ
1
+
δ
1
β
2
-
δ
1
δ
2
)
+
2
β
2
k
1
+
2
β
1
β
2
A
1
+
k
2
δ
1
+
β
2
A
2
δ
1
)
×
(
p
m
(
t
)
(
4
β
1
β
2
-
δ
1
δ
2
)
)
-
1
,
(29)
r
2
(
t
)
=
(
2
β
1
k
2
+
2
β
1
β
2
A
2
+
k
1
δ
2
+
β
1
A
2
δ
2
p
m
(
t
)
(
2
β
1
β
2
-
2
β
1
δ
2
+
δ
2
β
1
-
δ
1
δ
2
)
+
2
β
1
k
2
+
2
β
1
β
2
A
2
+
k
1
δ
2
+
β
1
A
2
δ
2
)
×
(
p
m
(
t
)
(
4
β
1
β
2
-
δ
1
δ
2
)
)
-
1
.
Substitution parameters value into system (17) can be concluded which is as follows:
(30)
p
m
(
t
+
1
)
=
p
m
(
t
)
+
α
p
m
(
t
)
(
17.403
-
0.908
p
m
(
t
)
)
,
r
1
(
t
)
=
0.398
+
1.617
p
m
(
t
)
,
r
2
(
t
)
=
0.391
+
1.703
p
m
(
t
)
.
From system (30) it can be concluded that the manufacture firstly makes decision and the decision variables directly relate to α. But collectors’ decision variables directly relate to p
m. The following simulates complex dynamics characteristics of system (30) through numerical simulation.


### 3.2.2. Numerical Simulation

(1)
α
Influence on the Collecting Market. When α
=
[
0,0.05
], the initial values of p
m, r
1, and r
2 are 20, 0.4, and 0.5; p
m, r
1, and r
2 are as shown in Figures 43(a), 43(b), and 43(c). From the figures, we can conclude that the system is in stable state when 0
<
α
<
0.113, and the stable values of p
m, r
1, and r
2 are 19.17, 0.4824, and 0.4799. When α
=
0.113, the system occurs the first bifurcation, then after cycle 2, 4 cycle, and so forth, the system is gradually into the chaotic state.(a)α = [0, 0.05], p
m, (b) α = [0, 0.05], r
1; and (c) α =  [0, 0.05], r
2.
(a)
(b)
(c)Figure 44
α =  [0, 0.05], max Lyapunov spectrum.(2) System Variables Power Spectrum. According to the numerical simulation, not only system variables graph in the phase space along with time can be drew, but also power spectrum graph of system variables can be estimated by period chart method. When α
=
0.15, power spectrum graph of the manufacturer and two collectors is shown in Figures 45(a), 45(b), and 45(c). According to the numerical results, no matter how large the collection price adjustment parameters are, the collection price traverses the whole value area over time, but p
m, r
1, and r
2 in the system (19) always limit in a certain range, which also verifies boundedness and ergodicity of chaos.(a)α
=
0.15, the power spectrum of p
m, (b) α
=
0.15, the spectrum of r
1; and (c) α
=
0.15, the power spectrum of r
2.
(a)
(b)
(c)(3) The Influence of Initial Setup. When α
=
0.15, sensitive dependence on initial value of p
m, r
1, and r
2 is as shown in Figures 46(a), 46(b), and 46(c). p
m values are, respectively, taken to be 20 and 20.001, after 36 cycles iteration, p
m difference value is 3.35 which is 3350 times the initial difference value. r
1 and r
2 difference values are 69 times and 73 times the initial difference values. Furthermore, from Figure 46 we can conclude that system variables values are approximate in the first period; with increase of iterations the difference increases obviously. The influence to the manufacturer is far greater than collectors because of the collectors’ decision-making later than the manufacturer’s. After the collectors observe action of the manufacturer, they make decision. The collectors’ decision-making is superior to the manufacturer.(a)α
=
0.15, the initial value sensitivity of p
m, (b) α
=
0.15, the initial value sensitivity of r
1
; and (c) α
=
0.15, the initial value sensitivity of r
2.
(a)
(b)
(c)


## 3.2.1. Model and the Analysis

Suppose that the manufacturer and collectors have principal and subordinate relationship, the manufacturer is Stackelberg leader, and collectors are followers; then they process sequential dynamic game; the game equilibrium is Stackelberg equilibrium. In this game, the manufacturer makes the decision of sales price and collection price according to the market informa Lyapunov exponent function abouttion; then two collectors make decision according to the decision-making of the manufacture:(24)
max
⁡
∏
r
1
=
(
p
m
-
p
r
1
(
t
)
-
A
1
)
(
k
1
+
β
1
p
r
1
(
t
)
-
δ
1
p
r
2
(
t
)
)
max
⁡
∏
r
2
=
(
p
m
-
p
r
2
(
t
)
-
A
2
)
(
k
2
+
β
2
p
r
2
(
t
)
-
δ
2
p
r
1
(
t
)
)
s
.
t
.
i
p
m
∈
arg
max
⁡
∏
r
1
s
.
t
.
i
p
m
∈
arg
max
⁡
∏
r
2
.The result is(25)
r
1
*
*
=
(
p
m
(
2
β
1
β
2
-
2
β
2
δ
1
+
δ
1
β
2
-
δ
1
δ
2
)
+
2
β
2
k
1
+
2
β
1
β
2
A
1
+
k
2
δ
1
+
β
2
A
2
δ
1
)
×
(
p
m
(
4
β
1
β
2
-
δ
1
δ
2
)
)
-
1
,
r
2
*
*
=
(
+
2
β
1
k
2
+
2
β
1
β
2
A
2
+
k
1
δ
2
+
β
1
A
2
δ
2
)
p
m
(
2
β
1
β
2
-
2
β
1
δ
2
+
δ
2
β
1
-
δ
1
δ
2
)
+
2
β
1
k
2
+
2
β
1
β
2
A
2
+
k
1
δ
2
+
β
1
A
2
δ
2
)
×
(
p
m
(
4
β
1
β
2
-
δ
1
δ
2
)
)
-
1
.Formula (25) is reaction function of the collectors. Put formula (25) into formula (2) and the following can be obtained:
(26)
p
m
*
*
=
(
p
0
-
c
r
)
2
-
(
(
β
1
-
δ
2
)
(
2
β
2
k
1
+
2
β
2
β
1
A
1
+
δ
1
k
2
+
δ
1
β
2
A
2
)
+
(
β
2
-
δ
1
)
(
2
β
1
k
2
+
2
β
2
β
1
A
2
+
δ
2
k
1
-
δ
2
β
1
A
1
)
-
(
4
β
1
β
2
-
δ
1
δ
2
)
(
k
1
+
k
2
)
)
×
(
(
2
β
1
β
2
+
2
β
2
δ
1
-
δ
1
β
2
)
(
β
2
-
δ
1
)
+
(
2
β
1
β
2
+
2
β
1
δ
2
-
δ
2
β
1
)
(
β
1
-
δ
2
)
)
-
1
.With Nash equilibrium value the parameters values arep
0
=
70, c
r
=
30, A
1
=
1, A
2
=
0.8, k
1
=
1, k
2
=
1, β
1
=
0.6, β
2
=
0.55, δ
1
=
0.2, δ
2
=
0.2. Then we put the values into Stackelberg equilibrium (
p
m
*
*
,
r
1
*
*
,
r
2
*
*
); the equilibrium value is (
19.17
,
0.48
,
0.48
).The manufacturer makes decision based on limited rational expectations. She adjusts the game process on the basis of marginal gains. If the marginal profits oft period are positive, in the t
+
1 period this strategy action will be used. And if the t period marginal profits are negative, in the t
+
1 period the collection price will be lowered. The repeated game model of dynamic adjustment is
(27)
p
m
(
t
+
1
)
=
p
m
(
t
)
+
α
p
m
(
t
)
×
(
(
(
1
-
r
1
)
(
β
1
-
δ
2
)
+
(
1
-
r
2
)
(
β
2
-
δ
1
)
)
×
(
p
0
-
2
p
m
(
t
)
-
c
r
)
-
k
1
-
k
2
(
(
1
-
r
1
)
(
β
1
-
δ
2
)
+
(
1
-
r
2
)
(
β
2
-
δ
1
)
)
)
.In formula (27), α is adjustment coefficient of collection price and r
1
,
r
2 are as follows:
(28)
r
1
(
t
)
=
(
2
β
2
k
1
+
2
β
1
β
2
A
1
+
k
2
δ
1
+
β
2
A
2
δ
1
p
m
(
t
)
(
2
β
1
β
2
-
2
β
2
δ
1
+
δ
1
β
2
-
δ
1
δ
2
)
+
2
β
2
k
1
+
2
β
1
β
2
A
1
+
k
2
δ
1
+
β
2
A
2
δ
1
)
×
(
p
m
(
t
)
(
4
β
1
β
2
-
δ
1
δ
2
)
)
-
1
,
(29)
r
2
(
t
)
=
(
2
β
1
k
2
+
2
β
1
β
2
A
2
+
k
1
δ
2
+
β
1
A
2
δ
2
p
m
(
t
)
(
2
β
1
β
2
-
2
β
1
δ
2
+
δ
2
β
1
-
δ
1
δ
2
)
+
2
β
1
k
2
+
2
β
1
β
2
A
2
+
k
1
δ
2
+
β
1
A
2
δ
2
)
×
(
p
m
(
t
)
(
4
β
1
β
2
-
δ
1
δ
2
)
)
-
1
.
Substitution parameters value into system (17) can be concluded which is as follows:
(30)
p
m
(
t
+
1
)
=
p
m
(
t
)
+
α
p
m
(
t
)
(
17.403
-
0.908
p
m
(
t
)
)
,
r
1
(
t
)
=
0.398
+
1.617
p
m
(
t
)
,
r
2
(
t
)
=
0.391
+
1.703
p
m
(
t
)
.
From system (30) it can be concluded that the manufacture firstly makes decision and the decision variables directly relate to α. But collectors’ decision variables directly relate to p
m. The following simulates complex dynamics characteristics of system (30) through numerical simulation.


## 3.2.2. Numerical Simulation

(1)
α
Influence on the Collecting Market. When α
=
[
0,0.05
], the initial values of p
m, r
1, and r
2 are 20, 0.4, and 0.5; p
m, r
1, and r
2 are as shown in Figures 43(a), 43(b), and 43(c). From the figures, we can conclude that the system is in stable state when 0
<
α
<
0.113, and the stable values of p
m, r
1, and r
2 are 19.17, 0.4824, and 0.4799. When α
=
0.113, the system occurs the first bifurcation, then after cycle 2, 4 cycle, and so forth, the system is gradually into the chaotic state.(a)α = [0, 0.05], p
m, (b) α = [0, 0.05], r
1; and (c) α =  [0, 0.05], r
2.
(a)
(b)
(c)Figure 44
α =  [0, 0.05], max Lyapunov spectrum.(2) System Variables Power Spectrum. According to the numerical simulation, not only system variables graph in the phase space along with time can be drew, but also power spectrum graph of system variables can be estimated by period chart method. When α
=
0.15, power spectrum graph of the manufacturer and two collectors is shown in Figures 45(a), 45(b), and 45(c). According to the numerical results, no matter how large the collection price adjustment parameters are, the collection price traverses the whole value area over time, but p
m, r
1, and r
2 in the system (19) always limit in a certain range, which also verifies boundedness and ergodicity of chaos.(a)α
=
0.15, the power spectrum of p
m, (b) α
=
0.15, the spectrum of r
1; and (c) α
=
0.15, the power spectrum of r
2.
(a)
(b)
(c)(3) The Influence of Initial Setup. When α
=
0.15, sensitive dependence on initial value of p
m, r
1, and r
2 is as shown in Figures 46(a), 46(b), and 46(c). p
m values are, respectively, taken to be 20 and 20.001, after 36 cycles iteration, p
m difference value is 3.35 which is 3350 times the initial difference value. r
1 and r
2 difference values are 69 times and 73 times the initial difference values. Furthermore, from Figure 46 we can conclude that system variables values are approximate in the first period; with increase of iterations the difference increases obviously. The influence to the manufacturer is far greater than collectors because of the collectors’ decision-making later than the manufacturer’s. After the collectors observe action of the manufacturer, they make decision. The collectors’ decision-making is superior to the manufacturer.(a)α
=
0.15, the initial value sensitivity of p
m, (b) α
=
0.15, the initial value sensitivity of r
1
; and (c) α
=
0.15, the initial value sensitivity of r
2.
(a)
(b)
(c)


## 3.3. Model Comparison and Analysis

(1) Profits Analysis of Nash Equilibrium. When α
=
[
0,0.05
], profits of the manufacturer and two collectors are as shown in Figures 47(a), 47(b), 47(c), and 47(d). Figure 47(a) shows that profits of the manufacturer in the stable state are not the best; it means Nash equilibrium is not optimal solution; the optimal solution of single cycle may not bring enterprise maximum returns. Furthermore, from Figure 47(a) it can be concluded that, with collecting price increasing, profits state instead into chaos state, the optimal price is 19.64 in the chaotic state, and the corresponding profits are 208.6. Compared to Figure 47(b), profits of the manufacture before chaotic state are 216.4, so the internal randomness of chaotic system cannot react to the actual situation of the market.(a)∏
m graph with p
m, (b) ∏
m graph with α
1, (c) ∏
r
1 graph with α
2, and (d) ∏
r
2 graph with α
3.
(a)
(b)
(c)
(d)(2) Profits Analysis of Stackelberg Equilibrium. When α
=
[
0,0.05
], profits of the manufacturer and two collectors are as shown in Figures 48(a), 48(b), 48(c), and 48(d). Figure 48 shows that profits of the manufacturer in the stable state are the best; it means Stackelberg equilibrium is the optimal solution, and the optimal solution of the closed-loop supply chain is optimal. Figures 48(a), 48(b), 48(c), and 48(d) are consistent.(a)∏
m graph with p
m, (b) ∏
m graph with α
1, (c) ∏
r
1 graph with α
2, and (d) ∏
r
2 graph with α
3.
(a)
(b)
(c)
(d)Through numerical simulation and comparison analysis, it can be concluded that Stackelberg game equilibrium is better than Nash equilibrium from the point of view of the profits and collecting price. Three conclusions can be summarized as follows. First, when the manufacturer and the two collectors make decision by static game, the collection price and profits of the three parties are lower than the optimal of dynamic game equilibrium. Second, as the manufacturer and the two collectors independently make decisions, the three parties will try to lower the collecting price in order to obtain the maximum profits. Third, as the manufacturer and the two collectors independently make decisions, the system more easily reaches chaos, and the cycle of stable state is shorter.


## 4. Centralized Control Decision-Making Model

Centralized control is that the manufacturers and the collectors codetermine to realize profits maximization of the supply chain system. We get formulas (2) and (3) into the formula (4) and attain the profits function of supply chain as follows:
(31)
∏
=
∏
m
+
∏
r
1
+
∏
r
2
=
k
1
(
p
0
-
c
r
-
A
1
)
+
k
2
(
p
0
-
c
r
-
A
2
)
+
p
m
(
t
)
(
1
-
r
1
(
t
)
)
×
(
(
p
0
-
c
r
)
(
1
-
δ
2
)
-
k
1
-
A
1
-
A
2
δ
2
)
-
p
m
(
t
)
2
(
β
1
(
1
-
r
1
(
t
)
)
2
+
β
2
(
1
-
r
2
(
t
)
)
2
)
+
p
m
(
t
)
(
1
-
r
2
(
t
)
)
×
(
(
p
0
-
c
r
)
(
1
-
δ
1
)
-
k
2
-
A
2
-
A
1
δ
1
)
+
p
m
(
t
)
2
(
δ
1
+
δ
2
)
(
1
-
r
1
(
t
)
)
(
1
-
r
2
(
t
)
)
.
The second-order derivatives about p
m, r
1, and r
2 in formula (31) are ∂
2
Π
/
∂
p
m
(
t
)
2
=
-
2
(
β
1
(
1
-
r
1
(
t
)
)
2 + β
2
(
1
-
r
2
(
t
)
)
2
)
+
2
(
δ
1
+
δ
2
)
(
1
-
r
1
(
t
)
)
(
1
-
r
2
(
t
)
)
∂
2
Π
/
∂
r
1
(
t
)
2 = -
2
β
1
p
m
(
t
)
2
<
0, ∂
2
Π
/
∂
r
2
(
t
)
2
=
-
2
β
2
p
m
(
t
)
2
<
0. When meeting the condition of (
δ
1
+
δ
2
)
(
1
-
r
1
(
t
)
)
(
1
-
r
2
(
t
)
)
<
β
1
(
1
-
r
1
(
t
)
)
2
+
β
2
(
1
-
r
2
(
t
)
)
2, Π is strictly concave function about p
m, r
1, and r
2. Therefore, the optimal solutions of p
m, r
1, and r
2 can be obtained by ∂
Π
/
∂
p
m
(
t
)
=
0, ∂
Π
/
∂
r
1
=
0
∂, and Π
/
∂
r
2
=
0 The process is as follows:
(32)
∂
Π
∂
p
m
=
(
1
-
r
1
(
t
)
)
(
(
p
0
-
c
r
)
(
1
-
δ
2
)
-
k
1
-
A
1
-
A
2
δ
2
)
+
(
1
-
r
2
(
t
)
)
(
(
p
0
-
c
r
)
(
1
-
δ
1
)
-
k
2
-
A
2
-
A
1
δ
1
)
-
2
p
m
(
t
)
(
β
1
(
1
-
r
1
(
t
)
)
2
+
β
2
(
1
-
r
2
(
t
)
)
2
)
+
2
p
m
(
t
)
(
δ
1
+
δ
2
)
(
1
-
r
1
(
t
)
)
(
1
-
r
2
(
t
)
)
=
0
,
∂
Π
∂
r
1
=
-
p
m
(
t
)
(
(
p
0
-
c
r
)
(
1
-
δ
2
)
-
k
1
-
A
1
-
A
2
δ
2
)
+
2
p
m
(
t
)
2
β
1
(
1
-
r
1
(
t
)
)
-
p
m
(
t
)
2
(
δ
1
+
δ
2
)
(
1
-
r
2
(
t
)
)
=
0
,
∂
Π
∂
r
2
=
-
p
m
(
t
)
(
(
p
0
-
c
r
)
(
1
-
δ
1
)
-
k
2
-
A
2
-
A
1
δ
1
)
+
2
p
m
(
t
)
2
β
2
(
1
-
r
2
(
t
)
)
-
p
m
(
t
)
2
(
δ
1
+
δ
2
)
(
1
-
r
1
(
t
)
)
=
0
.Value with above the top:p
0
=
90, c
r
=
30,A
1
=
1, A
2
=
0.8, k
1
=
1, k
2
=
1, β
1
=
0.6, β
2
=
0.55, δ
1
=
0.2, δ
2
=
0.2
.
p
m
(
1
)
=
20, r
1
(
1
)
=
0.35, and r
2
(
1
)
=
0.3. The equilibrium solution of formula (22) is (
p
-
m
,
r
-
1
,
r
-
2
)
=
(
20,0.35,0.35
).Manufacturers and collecting business make decision based on rational expectations, and the decision-making basis is marginal profits. Figure16 indicates that p
m, r
1, and r
2 three-dimensional changes with time. X, Y, and Z, respectively, show r
1, r
2, and p
m. Figure 49 describes system profits Π; r
1 and r
2 change with time. Maximum system profits are 919.3. In this condition, r
1, r
2 are, respectively, 0.1411 and 0.1411. In the centralized control decision-making model, we can conclude that the collecting business profits are suffered, but the manufacturer profits are superior to decentralized control decision-making model. Therefore, not only the system profits maximization but also business profits are considered in reality.(a)p
m, r
1, and r
2 conditions with time and (b) ∏
,
r
1
,
r
2 condition with time.
(a)
(b)


## 5. Conclusions

This paper researches on that used products collecting pricing game model and complexity analysis in a closed-loop supply chain which consist of a manufacturer and two collectors. Corresponding to the discrete system, we have established and analyzed in detail the corresponding continuity system.Through quantitative analysis of collecting price, time sequence diagram, and Lyapunov index, the paper describes game evolution rule in the closed-loop supply chain. We describe a series of chaotic system characteristics, based on which, we have a detailed analysis of the Lyapunov index which is under the condition of different parameter combination, and draw the figure under different conditions. The analysis found that, first, when the collection price is to a critical value, the system into chaos state. Second, the sale price of remanufacturing products are more than a critical value, the system into chaos state. Last, the collection price system of the initial value of which is sensitive, small changes in the initial value may cause market price fluctuations.It also draws conclusions that business profits, supply chain system profits, collecting price in Stackelberg game model are superior to Nash game model. In Nash equilibrium, the manufacturer and the collectors all lower the collecting price to gain, respectively, maximum profits. The system is easier into chaos state than Stackelberg, and stable state period is shorter.Through comparing centralized decision-making with decentralized decision-making, we have concluded that the collecting price and system profits of the latter are lower. This conclusion for the supply chain guide enterprises sure each phase of the best provide decision basis for the collection price, in order to achieve maximize returns.



---

*Source: 750179-2014-06-19.xml*