# Complexity Analysis and Control in Time-Delay Vaccine Supply Chain considering Cold Chain Transportation

Authors: Daoming Dai; Xuanyu Wu; Fengshan Si
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Due to the important role of the vaccines in the prevention of global epidemics, this paper focuses on a vaccine transportation supply chain composed of one distributor and one retailer. Based on the assumption that the decision-maker does not make a decision instantaneously, we present a decision-making time-delay model. Firstly, we captured some sufficient conditions of delay-induced bifurcation for the model by regarding different combinations of the decision delay periods as the bifurcation parameters and analyzed how the speed of decision-making adjustment of distributor or retailer affects the critical point of system stability. Secondly, we made a numerical simulation on the model by using a two-dimensional bifurcation diagram, largest Lyapunov exponent, and entropy and chaotic attractor, respectively. Finally, we used two coordination methods to control chaos and compared them. The results show that when the decision delay exceeds a certain threshold, the system will lose stability or go into chaos. The precipitous speed of decision variable adjustment of the distributor or retailer will increase the entropy of the system and lead the system into a chaotic state. When the vaccine supply chain is in chaos, the effect of external control on chaos is better than that of internal control on chaos.


---

## Body

## 1. Introduction

Due to the important role of vaccines in the prevention of infectious diseases spread, researchers have been interested in the vaccine supply chain from many perspectives. It has been proven that vaccination is the most effective method to stop the transmission of all kinds of infectious diseases [1]. For example, regular measles vaccination and supplementary immunization activities greatly reduce the incidence rate and mortality of measles worldwide. It was estimated that measles vaccination had averted about 23.2 million deaths during 2000–2018 [2]. Unfortunately, there are still tens of millions of children who have been not covered by routine vaccines in low and middle-income countries and there are also considerable challenges in achieving a high level of vaccination coverage. One of the main obstacles to the delivery of vaccines is that strict cold chain condition requires vaccines to be stored between 2°C and 8°C at all times [3].The transport conditions of vaccines are noticeable problems at present. Lee and Haidari [4] explored how vaccine supply chains may make substantial influence on the decision-making of ten examples of different members in the vaccine community. According to many open questions, including how the devices should be sized, the interaction between their price and size, and how often they should be replenished, Chen et al. [5] examined a model that captured the various tradeoffs in order to better guide the development of passive cold storage devices. Lemmens et al. [6] examined whether decisions at the strategic, tactical, and operational levels can figure out key issues in the vaccine supply chain, such as cold chain delivery, limited shelf life, and access to remote region. Shittu et al. [7] established a simulation model to examine the effects of variance in vaccine supply or demand on storage capacity requirements, and the results showed that the storage of vaccine can be reduced to 30% by reorganizing the supply chain. Leidner et al. [8] assessed the ability of different types of vaccine storage devices to maintain an appropriate temperature for vaccine storage and identify deviations from the recommended temperature. In addition, unreliable power systems and poor road conditions often lead to cold chain failures in many developing countries [9].From operation management perspective, some researchers have explored coordination of a vaccine supply chain. Tebbens et al. [10] proposed a mathematical framework of a vaccine stockpile over time and discussed the problems of capacity constraints, production and filling delay in the development and use of the polio vaccine reserve. Brison and Letallec [11] thought that one of the common factors limiting full and fair access to effective immunization was the gap between cold chain and logistics system. Lin et al. [12] examined that the distributor decided whether a cold chain or noncold chain is used to transport vaccines, while the retailer checked when the vaccine in a vaccine supply chain is received. Dai et al. [13] presented the delivery-time-dependent quantity flexibility contract and the late-rebate contract when there were uncertainties in the design, delivery, and demand of influenza vaccine based on the background of the influenza vaccine industry. Based on the profit goal and social responsibility of foreign vaccine suppliers, Niu et al. [14] explored two typical channel structures (i.e., exclusive retailing and competitive retailing) with the overseas vaccine supplier’s profit and social responsibility objectives and acquired equilibrium on vaccine prices, supply quantities, and vaccine supplier’s utilities.Many scholars have conducted chaotic behavior analysis on various aspects of supply chain. Xie and Ma [15] investigated the game between two recyclers and one processor in a duopoly market of color TV recycling and analyzed the response of the system to government decision-making. Wu and Ma [16] considered the dynamics of the game model of the epiphytic supply chain with two players and horizontal diversification of products and analyzed the equilibrium point and stable regions. Li et al. [17] established Nash game model and Stackelberg game model in multichannel supply chain and discussed the entropy diagram, the largest Lyapunov exponent and the chaotic attractor. Tu et al. [18] considered a dynamic hybrid supply chain model composed of two manufacturers and one retailer and explored the Nash equilibrium point and its stability region. Zhou and Chen [19] established two dynamic Stackelberg game models without (with) unit profit allocation and analyze the influence of relevant parameters on dynamic system stability, complexity, and channel revenue.Some scholars have been interested in stability problems of system with time delay. Zou et al. [20] established a delayed susceptible-infectious-recovered model for the transmission of porcine reproductive respiratory syndrome virus and acquired properties of the Hopf bifurcation by the normal form theory and center manifold theorem. Zhang et al. [21] proposed a delayed tobacco smoking model containing users in the form of snuffing and presented global exponential stability results for smoking present equilibrium by LMI techniques. Rajchakit and Rajchakit [22] investigated mean square robust stability of stochastic switched discrete-time-delay systems with convex polytopic uncertainties. Li and Liu [23] examined the bullwhip effect problem of supply chain with vendor order placement lead time delays, and the results showed that vendor order placement lead time delays played an important role in supply chain dynamics and contributed to its turbulence and volatility. References [24–28] analyzed the stability of several complex neural networks with kinds of time delays by using a Lyapunov-Krasovskii. In addition, Pratap et al. [29] explored the global existence of Filippov solutions and the robust generalized Mittag-Leffler synchronization of fractional order neural networks with discontinuous activation and impulses. References [30, 31] are devoted to the results of dissipative control for kinds of switched neutral time-delayed systems whose sufficient conditions for boundedness are obtained by using Lyapunov stability theory. References [32, 33] discussed the local stability and Hopf bifurcation of kinds of viruses (heroin or worm) model with time-delay and investigated the properties of the Hopf bifurcation by using center manifold theorem and the normal form theory. References [34, 35] investigated the problems of robust dissipative control for kinds of discrete-time systems and analyzed a set of sufficient conditions by using Lyapunov technique and linear matrix inequality approach.As far as we know, most existing researches on vaccine supply chains are related to medical research and public environmental health. Many scholars have been interested in vaccine supply chains from operation management perspective, but few of them are involved to stability of the channel system with time-delay. Many existing researches assume that decision-makers make decisions instantaneously, but decisions with time-delay often occur for information transmission delay. Thus, it is more practical that time-delay is introduced into a vaccine supply chain. In this paper, we consider a time-delay supply chain composed of a distributor and a retailer, in which the for-profit distributor purchases temperature-sensitive vaccines from the manufacturer and then resells them to the retailer. Because the vaccines become inactive in storage at unsuitable temperature, the distributor chooses the transportation mode with cold chain and increases the order quantity of the vaccines. The retailer often checks the effect of cold chain transportation and makes the distributor choose an appropriate cold chain transportation mode. We established a time-delay model of dynamic decision-making and analyzed the influence of time-delay parameter on the stability of supply chain system. We not only illustrated the Hopf bifurcation diagram of the system but also derived the condition that the adjustment speeds of the decision variables led to the period-doubling bifurcation of the system. In addition, we analyzed the chaotic behavior of the system by abundant simulation diagrams, such as Lyapunov exponent, entropy, and chaotic attractor.The remainder of this paper is organized as follows. The model conceptualization and the assumptions of the problem are described in Section2. The problem is modeled in Section 3. Equilibrium points and Hopf bifurcation are analyzed in Section 4. Numerical analysis and managerial insights are presented in Section 5. Chaos control is discussed in Section 6. Section 7 concludes this paper with discussion of the possible future research direction.


## 2. Description and Assumption of the Problem

### 2.1. Model Construction and Assumptions

This paper mainly considers a vaccine supply chain composed of a distributor and a retailer in which the for-profit distributor purchases temperature-sensitive vaccines from the manufacturer and then resells them to the retailer. Because the vaccines become inactive in storage at unsuitable temperature, the distributor chooses the cold chain mode of transportation and increases the order quantity of the vaccine. The retailer often checks the effects of cold chain transportation that the distributors chose to make them choose a cold chain transportation mode that is as appropriate as possible. The structure of vaccine supply chain system is shown in Figure1.Figure 1
The vaccine supply chain system.The main assumptions of this paper are as follows.(i)
The distributor and the retailer are bounded rational.(ii)
This paper only considers one-time investment. The distributor uses a cold chain to transport vaccines, the cold chain transportation level that needs to be invested isy, and distributor incurs the investment cost k1y2/2. The retailer checks whether vaccines are transported by cold chain with the inspection level s and incurs the inspection cost k2s2/2, where k1 and k2 are investment parameter and inspection parameter, respectively.(iii)
In order to ensure that distributor and the retailer can make normal profits, letp>w>c.(iv)
In order to ensure the vaccines’ activity, the distributor must choose the cold chain to transport the vaccines, and retailer will check the cold chain transportation of vaccines to ensure the vaccines’ activity.(v)
The influence ofy on demand function is greater than the influence of s, so η1>η2.


### 2.2. Symbolic Description

The meanings ofα,β,p,w, and c are described concisely in Table 1.Table 1
The parameters description for the system.
αThe initial demandpRetail price of unit vaccineβSensitivity coefficient of consumers to the retail pricewWholesale price of unit vaccinecThe initial demandThe demandD is considered as a function on transportation level, inspection level, and selling price of the vaccines. The functional form of market demand on the vaccines can be written as follows [36]:(1)D=α1−βpη1y+η2s,where η1 and η2 represent the influence coefficients of the distributor’s cold chain transportation level y and the retailer’s inspection level s on demand, respectively.


## 2.1. Model Construction and Assumptions

This paper mainly considers a vaccine supply chain composed of a distributor and a retailer in which the for-profit distributor purchases temperature-sensitive vaccines from the manufacturer and then resells them to the retailer. Because the vaccines become inactive in storage at unsuitable temperature, the distributor chooses the cold chain mode of transportation and increases the order quantity of the vaccine. The retailer often checks the effects of cold chain transportation that the distributors chose to make them choose a cold chain transportation mode that is as appropriate as possible. The structure of vaccine supply chain system is shown in Figure1.Figure 1
The vaccine supply chain system.The main assumptions of this paper are as follows.(i)
The distributor and the retailer are bounded rational.(ii)
This paper only considers one-time investment. The distributor uses a cold chain to transport vaccines, the cold chain transportation level that needs to be invested isy, and distributor incurs the investment cost k1y2/2. The retailer checks whether vaccines are transported by cold chain with the inspection level s and incurs the inspection cost k2s2/2, where k1 and k2 are investment parameter and inspection parameter, respectively.(iii)
In order to ensure that distributor and the retailer can make normal profits, letp>w>c.(iv)
In order to ensure the vaccines’ activity, the distributor must choose the cold chain to transport the vaccines, and retailer will check the cold chain transportation of vaccines to ensure the vaccines’ activity.(v)
The influence ofy on demand function is greater than the influence of s, so η1>η2.


## 2.2. Symbolic Description

The meanings ofα,β,p,w, and c are described concisely in Table 1.Table 1
The parameters description for the system.
αThe initial demandpRetail price of unit vaccineβSensitivity coefficient of consumers to the retail pricewWholesale price of unit vaccinecThe initial demandThe demandD is considered as a function on transportation level, inspection level, and selling price of the vaccines. The functional form of market demand on the vaccines can be written as follows [36]:(1)D=α1−βpη1y+η2s,where η1 and η2 represent the influence coefficients of the distributor’s cold chain transportation level y and the retailer’s inspection level s on demand, respectively.


## 3. Multiperiod Decision-Making Game Model with Delay

The profit functions of the distributor and retailer are expressed as follows:(2)πd=w−cα1−βpη1y+η2s−k1y22,(3)πr=p−wα1−βpη1y+η2s−k2s22.From (2) and (3), the decision variable of the distributor is y, and the retailer’s decision variables are p and s. Thus, the marginal profit function can be written as(4)∂πd∂y=η1w−cα1−βp−k1y,∂πr∂p=α1−2βp+βwη1y+η2s,∂πr∂s=η2p−wα1−βp−k2s.Because of the complexity of the market, the distributor and retailer do not fully know the market demand and the decisions of the competitors in the most cases. Therefore, they adjust their own decisions according to their marginal profit. The distributor and retailer make decisions based on bounded rationality. The values of decision variables in periodt + 1 under bounded rational decision are the values of decision variables at period t plus the changes of decision variables at period t.(5)yt+1=yt+v1yt∂πd∂y,pt+1=pt+v2pt∂πr∂p,st+1=st+v3st∂πr∂s,where v1,v2,v3 are the decision variables adjusting speed of the distributor and the retailer, respectively.In classical research on supply chains by using dynamics, scholars always assumed that decision-makers make decisions instantaneously; however, this assumption seems not to be realistic, because decision-makers are often bounded rational and risk averse and need to be considered before making a decision.τ is the time of the decision-maker's recognition of the need to take a decision to making a decision. So, we introduce the time-delay parameter into the differential equations.(6)yt+1=yt+v1ytη1w−cα1−βpt−k1yt−τ,pt+1=pt+v2ptα1−2βpt+βwη1yt−τ+η2st−τ,st+1=st+v3stη2pt−wα1−βpt−k2st−τ.The flow diagram of model (6) is as shown in Figure 2.Figure 2
The flow diagram of model (6).


## 4. Equilibrium Points and Hopf Bifurcation Analysis

### 4.1. Positive Equilibrium Characteristic and Characteristic Equation of Model (6)

We studied the equilibrium points of model (6). According to the calculations, we obtained nine equilibrium points of model (6):(7)E10,0,0,E20,0,−αη2wk2,E30,1β,0,E4αη1w−ck1,0,0,E50,1+βw2β,αη21−βp24βk2,E6αη1w−ck1,0,−αη2wk2,E7αη11−βpw−c2k1,1+βw2β,0,E8y1∗,p1∗,s1∗,E9y2∗,p2∗,s2∗,where(8)y1∗=αc−wβw−1η12k1,p1∗=βw+12β,s1∗=αβw−12η24βk2,y2∗=αc−wη1βk2η12c−w+k1η22βw−1k12η22,p2∗=k1η12c−w+wk1η22k1η22,s2∗=−αc−wη1βk2η12c−w+k1η22βw−1k12η23.In order to analyze the stability ofE1∼E8, the Jacobian matrix of model (6) is given as follows:(9)J=A1B10B2A2B30B4A3,where(10)A1=1+v1C1,A2=1+v2C2,A3=1+v3C3,B1=v1yη1αβc−w,B2=v2pαη11−2βp+βw,B3=v2pαη21−2βp+βw,B4=v3sαη21−βp−η2αβp−w,C1=η1w−cα1−βp−2k1y,C2=α1−4βp+βwη1y+η2s,C3=η2p−wα1−βp−2k2s.According to the eigenvalues of the Jacobian matrix evaluated at the corresponding equilibrium points, if all the nonzero eigenvalues of the Jacobian matrix are less than 1, this equilibrium point is stable.TakeE1 as an example. The values of E1 are taken into Jacobian matrix J, and then the Jacobian matrix J1 at the equilibrium point E1 is acquired. Thus, J1 can be written as(11)J1=1+v1η1w−cα00010001−v3η2wα.Obviously, the eigenvalues ofJ1 are not all less than 1, so the equilibrium point E10,0,0 is an unstable equilibrium point. Similarly, it can be proved that E2∼E7 are unstable points.Because of the high cost of vaccine production,w≫1 and βw≫1, so we can see from the calculation that y1∗<0 and the values of decision variables by the decision-maker is impossible to be nonpositive in economics, which means that the equilibrium point E8 is unstable.Next, we analyzed the stability of the equilibrium pointE9y2∗,p2∗,s2∗, and the characteristic polynomial of the Jacobian matrix in E9 can be written as(12)Fλ=λ3−T1λ2+T2λ+T3,where(13)T1=C1+C2+C3,T2=C1C2+C1C3+C2C3−B3B4−B1B2,T3=C1B3B4+C3B1B2−C1C2C3.According to the Jury criterion, all eigenvalues of the characteristic equation are in the unit circle of the complex plane, which is a necessary and sufficient condition for the stability of the discrete system. The conditions thatE9y2∗,p2∗,s2∗ is stable must be met as follows:(14)F1>0,F−1<0,1>T3,T32−1>−T3T1−T2.For the sake of convenience, letu1=yt−y∗,u2=pt−p∗,u3=st−s∗, and we can transform the stability of model (6) at the equilibrium point E9y2∗,p2∗,s2∗ into the stability at the point (0, 0, 0) to study. Let u1=yt,u2=pt,u3=st, and model (6) can be linearized at equilibrium point E9y2∗,p2∗,s2∗.(15)yt+1=a1yt+a2pt+b1yt−τ,pt+1=a3pt+b2yt−τ+b3st−τ,st+1=a4pt+a5st+b4st−τ,where(16)a1=1+v1η1w−cα1−βp∗−k1y∗,a2=−v1y∗η1βw−cα,a3=1+v2α1−4βp∗+βwη1y∗+η2s∗,a4=v3s∗αη21−2βp+βw,a5=1+v3η2p∗−wα1−βp∗−k2s∗,b1=−v1y∗k1,b2=v2p∗η1α1−2βp∗+βw,b3=v2p∗η2α1−2βp∗+βw,b4=−v3s∗k2.Next,(17)λ−a1−b1e−λτ−a20−b2e−λτλ−a3−b3e−λτ0−a4λ−a5−b4e−λτ.Then, we obtain the characteristic polynomial of model (6):(18)fλ=λ3−λ2t1+λt2−a1a3a5+−λ2b1+b4+λt3+t4e−λτ+λb1b4+t5e−2λ^τ,where(19)t1=a1+a3+a5,t2=a1a3+a1a5+a3a5,t3=a1b4+a5b1+a3b1+a3b4−a4b3−a2b2,t4=a1a4b3+a2a5b2−a1a3b4−a3a5b1,t5=a4b1b3+a2b2b4−a3b1b4.


### 4.2. Case 1,τ=0, Sufficient Conditions for Local Stability at Equilibrium Point E9y2∗,p2∗,s2∗

Whenτ=0, fλ can be simplified as follows:(20)λ3−λ2t1+b1+b4+λt2+t3+b1b4+t4+t5−a1a3a5.Letfλ=0; according to Routh-Hurwitz criterion, if H1:−t1+b1+b4>0,t2+t3+b1b4>0,−t1+b1+b4t2+t3+b1b4>t4+t5−a1a3a5 and t4+t5−a1a3a5>0, hold, then model (6) is asymptotically stable at the equilibrium point.


### 4.3. Case 2,τ>0, Sufficient Conditions for Local Stability at Equilibrium Point E9y2∗,p2∗,s2∗

Letfλ=0 and multiply both sides by eλτ, and then(21)−λ2b1+b4+λt3+t4+λ3−λ2t1+λt2−a1a3a5eλτ+λb1b4+t5e−λτ=0.Assume thatλ=iωω>0 is the root of (21); then,(22)Δ1cosωτ+Δ2sinωτ=−ωt3,Δ3cosωτ+Δ4sinωτ=−ω2b1+b4−t4,where(23)Δ1=−ω3+ωb1b4+ωt2,Δ2=ω2t1−a1a3a5−t5,Δ3=ω2t1−a1a3a5+t5,Δ4=ω3+ωb1b4−ωt2.By the aid of (22), we can obtain(24)sinωτ=−ω2b1+b4−t4Δ1+ωt3Δ3Δ1Δ4−Δ2Δ3,cosωτ=−ωt3Δ4+ω2b1+b4+t4Δ2Δ1Δ4−Δ2Δ3.By the aid of (24), we have(25)ω12+ω10r1+ω8r2+ω6r3+ω5r4+ω4r5+ω2r6+r7=0,where(26)r1=b1+b42+2t12−4t2,r2=6t22−2b1b42+t14−4t1a1a3a5−t12b1+b42−t32+2t1t3−2t4b1+b4+2b1b4+2t2b1+b42,r3=b1b4+t222b1b4−2t2−b1b4−t222b1b4+2t2−4t13a1a3a5−2a1a3a52+2t52−t322b1b4−2t2−2t12t4b1+b4+2t1b1+b42,r4=2t1t3b1+b4,r5=b1b4+t22b1b4−t22+2t12a1a3a52+t52+4t12a1a3a52−t52−t32b1b4−t22+t12t42+a1a3a5+t52b1+b42+t422b1b4+2t2,r6=−2t1a1a3a5+t52a1a3a5−t5−2t1a1a3a5−t52a1a3a5+t5+2t1t42a1a3a5+t52−t42b1+b42,r7=a1a3a5+t52a1a3a5−t52.Definegω=ω12+ω10r1+ω8r2+ω6r3+ω5r4+ω4r5+ω2r6+r7. Without loss of generality, we assume that H2:gω has at least one positive real root. In order to establish the main results of this article, we assume gω has n positive roots, denoted by g1,g2,…,gn,0<n≤12. From (24), we can get(27)τkj=1ωkarccos−ωt3Δ4+ω2b1+b4+t4Δ2Δ1Δ4−Δ2Δ3+2jπωk,k=1,2,3,…,n;j=0,1,….Let(28)τ0=minτkjk=1,2,3,…,n;j=0.1,…=minτk0k=1,2,3,…,n=τk00.Differentiating both sides of (21) with regard to τ, we can get(29)dλdτ−1=−2λb1+b4+t3+3λ2−2λt1+t2eλτ+b1b4e−λτλeλτλ3−λ2t1+λt2−a1a3a5−λe−λτλb1b4+t5−τλ.Whenτ=τ0 and substituting λ=iω0 into (31), it can be achieved that(30)Redλdττ=τ0−1=X1X2+Y1Y2X12+Y12,where(31)X1=ω04−ω02t2+ω02b1b4cosω0τ0+−ω03t1+ω0a1a3a5−ω0t5sinω0τ0,Y1=ω03t1−ω0a1a3a5−ω0t5cosω0τ0+ω04−ω02t2−ω02b1b4sinω0τ0,X2=−3ω02+t2+b1b4cosω0τ0+2ω0t1sinω0τ0+t3,Y2=−2ω0t1cosω0τ0+−3ω02+t2−b1b4sinω0τ0−2ω0b1+b4.To ensure the condition of the occurrence for Hopf bifurcation, we have the following hypothesisH3:X1X2+Y1Y2≠0, and then we have the following results.Theorem 1.
For model (6), if the conditions H1,H2, and H3 hold, the equilibrium point E9y2∗,p2∗,s2∗ is asymptotically stable when τ∈0,τ0; when τ=τ0, model (6) undergoes a Hopf bifurcation at equilibrium point E9y2∗,p2∗,s2∗; and model (6) is unstable at equilibrium point E9y2∗,p2∗,s2∗ whenτ>τ0.


## 4.1. Positive Equilibrium Characteristic and Characteristic Equation of Model (6)

We studied the equilibrium points of model (6). According to the calculations, we obtained nine equilibrium points of model (6):(7)E10,0,0,E20,0,−αη2wk2,E30,1β,0,E4αη1w−ck1,0,0,E50,1+βw2β,αη21−βp24βk2,E6αη1w−ck1,0,−αη2wk2,E7αη11−βpw−c2k1,1+βw2β,0,E8y1∗,p1∗,s1∗,E9y2∗,p2∗,s2∗,where(8)y1∗=αc−wβw−1η12k1,p1∗=βw+12β,s1∗=αβw−12η24βk2,y2∗=αc−wη1βk2η12c−w+k1η22βw−1k12η22,p2∗=k1η12c−w+wk1η22k1η22,s2∗=−αc−wη1βk2η12c−w+k1η22βw−1k12η23.In order to analyze the stability ofE1∼E8, the Jacobian matrix of model (6) is given as follows:(9)J=A1B10B2A2B30B4A3,where(10)A1=1+v1C1,A2=1+v2C2,A3=1+v3C3,B1=v1yη1αβc−w,B2=v2pαη11−2βp+βw,B3=v2pαη21−2βp+βw,B4=v3sαη21−βp−η2αβp−w,C1=η1w−cα1−βp−2k1y,C2=α1−4βp+βwη1y+η2s,C3=η2p−wα1−βp−2k2s.According to the eigenvalues of the Jacobian matrix evaluated at the corresponding equilibrium points, if all the nonzero eigenvalues of the Jacobian matrix are less than 1, this equilibrium point is stable.TakeE1 as an example. The values of E1 are taken into Jacobian matrix J, and then the Jacobian matrix J1 at the equilibrium point E1 is acquired. Thus, J1 can be written as(11)J1=1+v1η1w−cα00010001−v3η2wα.Obviously, the eigenvalues ofJ1 are not all less than 1, so the equilibrium point E10,0,0 is an unstable equilibrium point. Similarly, it can be proved that E2∼E7 are unstable points.Because of the high cost of vaccine production,w≫1 and βw≫1, so we can see from the calculation that y1∗<0 and the values of decision variables by the decision-maker is impossible to be nonpositive in economics, which means that the equilibrium point E8 is unstable.Next, we analyzed the stability of the equilibrium pointE9y2∗,p2∗,s2∗, and the characteristic polynomial of the Jacobian matrix in E9 can be written as(12)Fλ=λ3−T1λ2+T2λ+T3,where(13)T1=C1+C2+C3,T2=C1C2+C1C3+C2C3−B3B4−B1B2,T3=C1B3B4+C3B1B2−C1C2C3.According to the Jury criterion, all eigenvalues of the characteristic equation are in the unit circle of the complex plane, which is a necessary and sufficient condition for the stability of the discrete system. The conditions thatE9y2∗,p2∗,s2∗ is stable must be met as follows:(14)F1>0,F−1<0,1>T3,T32−1>−T3T1−T2.For the sake of convenience, letu1=yt−y∗,u2=pt−p∗,u3=st−s∗, and we can transform the stability of model (6) at the equilibrium point E9y2∗,p2∗,s2∗ into the stability at the point (0, 0, 0) to study. Let u1=yt,u2=pt,u3=st, and model (6) can be linearized at equilibrium point E9y2∗,p2∗,s2∗.(15)yt+1=a1yt+a2pt+b1yt−τ,pt+1=a3pt+b2yt−τ+b3st−τ,st+1=a4pt+a5st+b4st−τ,where(16)a1=1+v1η1w−cα1−βp∗−k1y∗,a2=−v1y∗η1βw−cα,a3=1+v2α1−4βp∗+βwη1y∗+η2s∗,a4=v3s∗αη21−2βp+βw,a5=1+v3η2p∗−wα1−βp∗−k2s∗,b1=−v1y∗k1,b2=v2p∗η1α1−2βp∗+βw,b3=v2p∗η2α1−2βp∗+βw,b4=−v3s∗k2.Next,(17)λ−a1−b1e−λτ−a20−b2e−λτλ−a3−b3e−λτ0−a4λ−a5−b4e−λτ.Then, we obtain the characteristic polynomial of model (6):(18)fλ=λ3−λ2t1+λt2−a1a3a5+−λ2b1+b4+λt3+t4e−λτ+λb1b4+t5e−2λ^τ,where(19)t1=a1+a3+a5,t2=a1a3+a1a5+a3a5,t3=a1b4+a5b1+a3b1+a3b4−a4b3−a2b2,t4=a1a4b3+a2a5b2−a1a3b4−a3a5b1,t5=a4b1b3+a2b2b4−a3b1b4.


## 4.2. Case 1,τ=0, Sufficient Conditions for Local Stability at Equilibrium Point E9y2∗,p2∗,s2∗

Whenτ=0, fλ can be simplified as follows:(20)λ3−λ2t1+b1+b4+λt2+t3+b1b4+t4+t5−a1a3a5.Letfλ=0; according to Routh-Hurwitz criterion, if H1:−t1+b1+b4>0,t2+t3+b1b4>0,−t1+b1+b4t2+t3+b1b4>t4+t5−a1a3a5 and t4+t5−a1a3a5>0, hold, then model (6) is asymptotically stable at the equilibrium point.


## 4.3. Case 2,τ>0, Sufficient Conditions for Local Stability at Equilibrium Point E9y2∗,p2∗,s2∗

Letfλ=0 and multiply both sides by eλτ, and then(21)−λ2b1+b4+λt3+t4+λ3−λ2t1+λt2−a1a3a5eλτ+λb1b4+t5e−λτ=0.Assume thatλ=iωω>0 is the root of (21); then,(22)Δ1cosωτ+Δ2sinωτ=−ωt3,Δ3cosωτ+Δ4sinωτ=−ω2b1+b4−t4,where(23)Δ1=−ω3+ωb1b4+ωt2,Δ2=ω2t1−a1a3a5−t5,Δ3=ω2t1−a1a3a5+t5,Δ4=ω3+ωb1b4−ωt2.By the aid of (22), we can obtain(24)sinωτ=−ω2b1+b4−t4Δ1+ωt3Δ3Δ1Δ4−Δ2Δ3,cosωτ=−ωt3Δ4+ω2b1+b4+t4Δ2Δ1Δ4−Δ2Δ3.By the aid of (24), we have(25)ω12+ω10r1+ω8r2+ω6r3+ω5r4+ω4r5+ω2r6+r7=0,where(26)r1=b1+b42+2t12−4t2,r2=6t22−2b1b42+t14−4t1a1a3a5−t12b1+b42−t32+2t1t3−2t4b1+b4+2b1b4+2t2b1+b42,r3=b1b4+t222b1b4−2t2−b1b4−t222b1b4+2t2−4t13a1a3a5−2a1a3a52+2t52−t322b1b4−2t2−2t12t4b1+b4+2t1b1+b42,r4=2t1t3b1+b4,r5=b1b4+t22b1b4−t22+2t12a1a3a52+t52+4t12a1a3a52−t52−t32b1b4−t22+t12t42+a1a3a5+t52b1+b42+t422b1b4+2t2,r6=−2t1a1a3a5+t52a1a3a5−t5−2t1a1a3a5−t52a1a3a5+t5+2t1t42a1a3a5+t52−t42b1+b42,r7=a1a3a5+t52a1a3a5−t52.Definegω=ω12+ω10r1+ω8r2+ω6r3+ω5r4+ω4r5+ω2r6+r7. Without loss of generality, we assume that H2:gω has at least one positive real root. In order to establish the main results of this article, we assume gω has n positive roots, denoted by g1,g2,…,gn,0<n≤12. From (24), we can get(27)τkj=1ωkarccos−ωt3Δ4+ω2b1+b4+t4Δ2Δ1Δ4−Δ2Δ3+2jπωk,k=1,2,3,…,n;j=0,1,….Let(28)τ0=minτkjk=1,2,3,…,n;j=0.1,…=minτk0k=1,2,3,…,n=τk00.Differentiating both sides of (21) with regard to τ, we can get(29)dλdτ−1=−2λb1+b4+t3+3λ2−2λt1+t2eλτ+b1b4e−λτλeλτλ3−λ2t1+λt2−a1a3a5−λe−λτλb1b4+t5−τλ.Whenτ=τ0 and substituting λ=iω0 into (31), it can be achieved that(30)Redλdττ=τ0−1=X1X2+Y1Y2X12+Y12,where(31)X1=ω04−ω02t2+ω02b1b4cosω0τ0+−ω03t1+ω0a1a3a5−ω0t5sinω0τ0,Y1=ω03t1−ω0a1a3a5−ω0t5cosω0τ0+ω04−ω02t2−ω02b1b4sinω0τ0,X2=−3ω02+t2+b1b4cosω0τ0+2ω0t1sinω0τ0+t3,Y2=−2ω0t1cosω0τ0+−3ω02+t2−b1b4sinω0τ0−2ω0b1+b4.To ensure the condition of the occurrence for Hopf bifurcation, we have the following hypothesisH3:X1X2+Y1Y2≠0, and then we have the following results.Theorem 1.
For model (6), if the conditions H1,H2, and H3 hold, the equilibrium point E9y2∗,p2∗,s2∗ is asymptotically stable when τ∈0,τ0; when τ=τ0, model (6) undergoes a Hopf bifurcation at equilibrium point E9y2∗,p2∗,s2∗; and model (6) is unstable at equilibrium point E9y2∗,p2∗,s2∗ whenτ>τ0.


## 5. Numerical Simulation

In order to further verify the obtained theoretical results, two-dimensional bifurcation diagram is used to analyze the influence of the time-delay and adjustment speed of decision variables on the stability of the system. We use the largest Lyapunov exponent and entropy to measure system complexity, respectively. The principle of the largest Lyapunov exponent is that when exponent value is less than zero the system is stable; in contrast, when the exponent value is greater than zero, the system is unstable. The rules of entropy to judge system complexity is as follows: the system is in a stable state when the entropy is low; in contrast, the system is in chaos when the entropy is large.


### 5.1. Hopf Bifurcation Diagram Caused by Delay

We set the parameters as follows:α=1,β=0.29,p=200,w=100,c=50,η1=0.6,η2=0.3,k1=0.1,k2=0.4 and by some computations we obtained τ0=0.304, ω0=1.68. First, the two-dimensional bifurcation diagram of the system stability domain changes with the time-delay as shown in Figure 3(a).Figure 3
The impact ofτ on the stability and complexity of model (6). (a) Bifurcation diagram of model (6) with (b) Lyapunov exponents with respect to ττ of Figure 3(a). (c) The entropy diagram of Figure 3(a).
(a)(b)(c)As can be seen from Figure3(a), model (6) gradually produces Hopf bifurcation from stability with the increase of τ and then changes to unstable state, and the critical point of bifurcation is τ0=0.304. When τ<0.304, equilibrium point E9y2∗,p2∗,s2∗ is asymptotically stable; and when τ>0.304, model (6) goes from Hopf bifurcation to the unstable state at equilibrium point E9y2∗,p2∗,s2∗.From Figure3(a), it has been shown that the vaccine supply chain equilibrium is locally asymptotically stable if the time-delay parameter below the critical value. The advantage of a stable state is that the decision-makers in the system generally know the development law of the decision variables; at this time, the decision-maker can control the vaccine transportation more easily. However, once the time-delay parameter exceeds a certain threshold, the system undergoes a Hopf bifurcation and goes into chaos. The emergence of Hopf bifurcation means that model (6) changes from the vaccines transportation equilibrium to a limit cycle. At this time, the vaccine transportation is out of control. Figure 3(b) shows the largest Lyapunov exponent when τ increases. Comparing Figure 3(a) with Figure 3(b), when bifurcation occurs and the system gets in chaos, it can be seen that the largest Lyapunov exponents are consistent with the bifurcation diagram shown in Figure 3(a). Similarly, comparing Figure 3(a) with Figure 3(c), it can be seen that the entropy of the system is increasing rapidly when bifurcation occurs and gets in chaos.


### 5.2. Bifurcation Diagram Caused by Adjusting Speed of Decision Variable

In this paper, we show the dynamic characteristics of the system in four ways: bifurcation diagram, maximum Lyapunov exponent, entropy, and chaos attractor.


#### 5.2.1. Chaos and Bifurcation Phenomenon

Whenτ=0.1 and other parameters remain unchanged, the influence of time-delay on the stability of the system can be eliminated. Assume v2=0.7, and v3=0.6, and the value of v1 ranges from 0 to 0.5. We examine the effect of the adjusting speed v1 of the cold chain transportation level on the system stability. Figure 4(a) demonstrates that when v1 increases from 0 to 0.251, the cold chain transportation level y, retail price of unit vaccine p, and inspection level s are stable at 5.2, 2.472, and 1.837, respectively; by contrast, when v1 increases to 0.251, the system turns into stable cycles of period 2, which means that each decision’s variables have two possible values; when v1>0.38, the cold chain transportation level y, retail price of unit vaccine p, and inspection level s have a lot of possible values.Figure 4
Dynamic evolution of decision variables withv1. (a) v2=0.7,v3=0.6. (b) Lyapunov exponents for v1 varying from 0 to 0.5 of Figure 4(a). (c) The entropy diagram of Figure 4(a). (d) Chaotic attractor of the system.
(a)(b)(c)(d)Ceteris paribus, assuming thatv1=0.45, and v3=0.32, Figure 5 presents the effect of the adjustment speed v2 of retail price of vaccine on system stability. When v2 changes at 0,0.5, the changes of decision variables are shown in Figure 5(a). When v2∈0,0.166, the cold chain transportation level y, retail price of unit vaccine p, and inspection level s are stable at 1.148, 1.627, and 2.809, respectively; when v2 increases to 0.166, the first bifurcation occurs, and the system turns into stable cycles of period 2, which means that each decision variables have two possible values; with the further increase of v2, the cold chain transportation levely, retail price of unit vaccine p. and inspection level s become chaotic; that is, each decision variables have a lot of possible values.Figure 5
Dynamic evolution of decision variables withv2. (a) v1=0.45,v3=0.32. (b) Lyapunov exponents for v2 varying from 0 to 0.5 of Figure 5(a). (c) The entropy diagram of Figure 5(a). (d) Chaotic attractor of the system.
(a)(b)(c)(d)Furthermore, assuming thatv1=0.19 and v2=0.36, Figure 6 analyzes the effect of v3 on system stability. When v3 changes at [0, 0.5], the changes of decision variables are shown in Figure 6(a). When v3∈0,0.151, the cold chain transportation level y, retail price of unit vaccine p, and inspection level s are stable at 2.492, 1.814, and 1.844, respectively; when v3 increases to 0.151, the system turns into stable cycles of period 2; with the increase of v3, eventually the cold chain transportation level y, retail price of unit vaccine p, and inspection level s become chaotic; that is, each decision’s variables have a lot of possible values.Figure 6
Dynamic evolution of decision variables withv3. (a) v1=0.19,v3=0.36. (b) Lyapunov exponents for v3 varying from 0 to 0.5 of Figure 6(a). (c) The entropy diagram of Figure 6(a). (d) Chaotic attractor of the system.
(a)(b)(c)(d)From Figures4(a)–6(a), it can be seen that as the adjustment speed of decision variables increases, the number of system cycles increases, and finally the chaotic state appears. The increasing of adjustment speed of decision variables means that there are irrational decisions in the vaccine supply chain. It is difficult to obtain the sustainable development of the vaccine supply chain structure.


#### 5.2.2. Largest Lyapunov Exponent

The largest Lyapunov exponent can be used to determine whether the system is in a chaotic state. If the largest Lyapunov exponent is less than 0, the system is in a stable state; if the largest Lyapunov exponent is close to 0, the system is in a periodic motion; if the largest Lyapunov exponent is greater than 0, the system is in a chaotic state. Take Figure4(b) for example; when the largest Lyapunov exponent first returns to the 0 axis, the system reaches two times periodic bifurcation point. When the system returns for the second time and passes through the 0 axis, the system enters a chaotic state. Comparing the largest Lyapunov exponent with Figure 4(a), when bifurcation occurs and the system becomes in chaos, it can be seen that the largest Lyapunov exponents are consistent with the bifurcation diagram shown in Figure 4(a). Figures 5(b)–6(b) show similar results, which are not repeated here.


#### 5.2.3. Entropy

When many scholars analyze the influence of time-delay parameter in the system, they only prove and give the conditions for the stability of the system. However, there are few people discussing what the chaotic behavior will bring to the system. The essence of entropy is the internal chaos of a system. According to information theory, the higher the entropy is, the higher the chaos of the system’s internal information is, which also results in the decision-maker making more effort to acquire useful information. From Figure4(c), we can know that when v1∈0,0.251, the entropy of model (6) is equal to 0, and the cold chain transportation levely, retail price of unit vaccine p, or inspection level s has a unique value. When v1>0.251, model (6) undergoes a doubly periodic bifurcation state and falls into chaos finally, the entropy of model (6) continues to increase, and the cold chain transportation level y, retail price of unit vaccine p, or inspection level s becomes unstable and has lots of values. From the information theory, the large entropy can increase uncertainty of information on decision variables, which may cause the distributor and retailer to miss the best decision-making opportunity to maintain vaccines’ activity in the process of vaccine transportation. Figures 5(c)–6(c) show similar results.


#### 5.2.4. Chaotic Attractor

The chaotic attractors are shown in Figures4(d)–6(d) and are analyzed. When the system goes into a chaotic state, the structure of the chaotic attractor can be more complicated. Take Figure 4(d) for example; when v1>0.251, the system is in a chaotic state, and the chaotic attractor is complicated. Figures 5(d)–6(d) show similar results.


### 5.3. Influence of Adjustment Speed of Decision Variables on Stability Critical Point

Compared with Figure4(a), Figure 7(a) shows the change of system stability with v1 when v3 is constant and v2 is reduced to 0.3. Obviously, with the decrease of retailer’s adjustment speed v2, the critical point of stability moves to the right. Figure 7(b) shows the change of system stability with v1 when v2 is constant and v3 is reduced to 0.27. In the same way, with the decrease of retailer’s adjustment speed v3, the critical point of the stability moves to the right. It is not hard to see that, in terms of the goal of affecting the critical point of the stable state of the vaccine supply chain, there exist reverse changes in the relationship between the adjustment speed of the distributor’s decision variable v1 and the adjustment speed of the retailer’s decision variable v2 andv3. This interesting conclusion can imply powerful suggestions for distributor and retailer to find a balance in the process of vaccines transportation.Figure 7
Dynamic evolution of decision variables withv1. (a) v2 is reduced to 0.3. (b) v3 is reduced to 0.27.
(a)(b)


## 5.1. Hopf Bifurcation Diagram Caused by Delay

We set the parameters as follows:α=1,β=0.29,p=200,w=100,c=50,η1=0.6,η2=0.3,k1=0.1,k2=0.4 and by some computations we obtained τ0=0.304, ω0=1.68. First, the two-dimensional bifurcation diagram of the system stability domain changes with the time-delay as shown in Figure 3(a).Figure 3
The impact ofτ on the stability and complexity of model (6). (a) Bifurcation diagram of model (6) with (b) Lyapunov exponents with respect to ττ of Figure 3(a). (c) The entropy diagram of Figure 3(a).
(a)(b)(c)As can be seen from Figure3(a), model (6) gradually produces Hopf bifurcation from stability with the increase of τ and then changes to unstable state, and the critical point of bifurcation is τ0=0.304. When τ<0.304, equilibrium point E9y2∗,p2∗,s2∗ is asymptotically stable; and when τ>0.304, model (6) goes from Hopf bifurcation to the unstable state at equilibrium point E9y2∗,p2∗,s2∗.From Figure3(a), it has been shown that the vaccine supply chain equilibrium is locally asymptotically stable if the time-delay parameter below the critical value. The advantage of a stable state is that the decision-makers in the system generally know the development law of the decision variables; at this time, the decision-maker can control the vaccine transportation more easily. However, once the time-delay parameter exceeds a certain threshold, the system undergoes a Hopf bifurcation and goes into chaos. The emergence of Hopf bifurcation means that model (6) changes from the vaccines transportation equilibrium to a limit cycle. At this time, the vaccine transportation is out of control. Figure 3(b) shows the largest Lyapunov exponent when τ increases. Comparing Figure 3(a) with Figure 3(b), when bifurcation occurs and the system gets in chaos, it can be seen that the largest Lyapunov exponents are consistent with the bifurcation diagram shown in Figure 3(a). Similarly, comparing Figure 3(a) with Figure 3(c), it can be seen that the entropy of the system is increasing rapidly when bifurcation occurs and gets in chaos.


## 5.2. Bifurcation Diagram Caused by Adjusting Speed of Decision Variable

In this paper, we show the dynamic characteristics of the system in four ways: bifurcation diagram, maximum Lyapunov exponent, entropy, and chaos attractor.


### 5.2.1. Chaos and Bifurcation Phenomenon

Whenτ=0.1 and other parameters remain unchanged, the influence of time-delay on the stability of the system can be eliminated. Assume v2=0.7, and v3=0.6, and the value of v1 ranges from 0 to 0.5. We examine the effect of the adjusting speed v1 of the cold chain transportation level on the system stability. Figure 4(a) demonstrates that when v1 increases from 0 to 0.251, the cold chain transportation level y, retail price of unit vaccine p, and inspection level s are stable at 5.2, 2.472, and 1.837, respectively; by contrast, when v1 increases to 0.251, the system turns into stable cycles of period 2, which means that each decision’s variables have two possible values; when v1>0.38, the cold chain transportation level y, retail price of unit vaccine p, and inspection level s have a lot of possible values.Figure 4
Dynamic evolution of decision variables withv1. (a) v2=0.7,v3=0.6. (b) Lyapunov exponents for v1 varying from 0 to 0.5 of Figure 4(a). (c) The entropy diagram of Figure 4(a). (d) Chaotic attractor of the system.
(a)(b)(c)(d)Ceteris paribus, assuming thatv1=0.45, and v3=0.32, Figure 5 presents the effect of the adjustment speed v2 of retail price of vaccine on system stability. When v2 changes at 0,0.5, the changes of decision variables are shown in Figure 5(a). When v2∈0,0.166, the cold chain transportation level y, retail price of unit vaccine p, and inspection level s are stable at 1.148, 1.627, and 2.809, respectively; when v2 increases to 0.166, the first bifurcation occurs, and the system turns into stable cycles of period 2, which means that each decision variables have two possible values; with the further increase of v2, the cold chain transportation levely, retail price of unit vaccine p. and inspection level s become chaotic; that is, each decision variables have a lot of possible values.Figure 5
Dynamic evolution of decision variables withv2. (a) v1=0.45,v3=0.32. (b) Lyapunov exponents for v2 varying from 0 to 0.5 of Figure 5(a). (c) The entropy diagram of Figure 5(a). (d) Chaotic attractor of the system.
(a)(b)(c)(d)Furthermore, assuming thatv1=0.19 and v2=0.36, Figure 6 analyzes the effect of v3 on system stability. When v3 changes at [0, 0.5], the changes of decision variables are shown in Figure 6(a). When v3∈0,0.151, the cold chain transportation level y, retail price of unit vaccine p, and inspection level s are stable at 2.492, 1.814, and 1.844, respectively; when v3 increases to 0.151, the system turns into stable cycles of period 2; with the increase of v3, eventually the cold chain transportation level y, retail price of unit vaccine p, and inspection level s become chaotic; that is, each decision’s variables have a lot of possible values.Figure 6
Dynamic evolution of decision variables withv3. (a) v1=0.19,v3=0.36. (b) Lyapunov exponents for v3 varying from 0 to 0.5 of Figure 6(a). (c) The entropy diagram of Figure 6(a). (d) Chaotic attractor of the system.
(a)(b)(c)(d)From Figures4(a)–6(a), it can be seen that as the adjustment speed of decision variables increases, the number of system cycles increases, and finally the chaotic state appears. The increasing of adjustment speed of decision variables means that there are irrational decisions in the vaccine supply chain. It is difficult to obtain the sustainable development of the vaccine supply chain structure.


### 5.2.2. Largest Lyapunov Exponent

The largest Lyapunov exponent can be used to determine whether the system is in a chaotic state. If the largest Lyapunov exponent is less than 0, the system is in a stable state; if the largest Lyapunov exponent is close to 0, the system is in a periodic motion; if the largest Lyapunov exponent is greater than 0, the system is in a chaotic state. Take Figure4(b) for example; when the largest Lyapunov exponent first returns to the 0 axis, the system reaches two times periodic bifurcation point. When the system returns for the second time and passes through the 0 axis, the system enters a chaotic state. Comparing the largest Lyapunov exponent with Figure 4(a), when bifurcation occurs and the system becomes in chaos, it can be seen that the largest Lyapunov exponents are consistent with the bifurcation diagram shown in Figure 4(a). Figures 5(b)–6(b) show similar results, which are not repeated here.


### 5.2.3. Entropy

When many scholars analyze the influence of time-delay parameter in the system, they only prove and give the conditions for the stability of the system. However, there are few people discussing what the chaotic behavior will bring to the system. The essence of entropy is the internal chaos of a system. According to information theory, the higher the entropy is, the higher the chaos of the system’s internal information is, which also results in the decision-maker making more effort to acquire useful information. From Figure4(c), we can know that when v1∈0,0.251, the entropy of model (6) is equal to 0, and the cold chain transportation levely, retail price of unit vaccine p, or inspection level s has a unique value. When v1>0.251, model (6) undergoes a doubly periodic bifurcation state and falls into chaos finally, the entropy of model (6) continues to increase, and the cold chain transportation level y, retail price of unit vaccine p, or inspection level s becomes unstable and has lots of values. From the information theory, the large entropy can increase uncertainty of information on decision variables, which may cause the distributor and retailer to miss the best decision-making opportunity to maintain vaccines’ activity in the process of vaccine transportation. Figures 5(c)–6(c) show similar results.


### 5.2.4. Chaotic Attractor

The chaotic attractors are shown in Figures4(d)–6(d) and are analyzed. When the system goes into a chaotic state, the structure of the chaotic attractor can be more complicated. Take Figure 4(d) for example; when v1>0.251, the system is in a chaotic state, and the chaotic attractor is complicated. Figures 5(d)–6(d) show similar results.


## 5.2.1. Chaos and Bifurcation Phenomenon

Whenτ=0.1 and other parameters remain unchanged, the influence of time-delay on the stability of the system can be eliminated. Assume v2=0.7, and v3=0.6, and the value of v1 ranges from 0 to 0.5. We examine the effect of the adjusting speed v1 of the cold chain transportation level on the system stability. Figure 4(a) demonstrates that when v1 increases from 0 to 0.251, the cold chain transportation level y, retail price of unit vaccine p, and inspection level s are stable at 5.2, 2.472, and 1.837, respectively; by contrast, when v1 increases to 0.251, the system turns into stable cycles of period 2, which means that each decision’s variables have two possible values; when v1>0.38, the cold chain transportation level y, retail price of unit vaccine p, and inspection level s have a lot of possible values.Figure 4
Dynamic evolution of decision variables withv1. (a) v2=0.7,v3=0.6. (b) Lyapunov exponents for v1 varying from 0 to 0.5 of Figure 4(a). (c) The entropy diagram of Figure 4(a). (d) Chaotic attractor of the system.
(a)(b)(c)(d)Ceteris paribus, assuming thatv1=0.45, and v3=0.32, Figure 5 presents the effect of the adjustment speed v2 of retail price of vaccine on system stability. When v2 changes at 0,0.5, the changes of decision variables are shown in Figure 5(a). When v2∈0,0.166, the cold chain transportation level y, retail price of unit vaccine p, and inspection level s are stable at 1.148, 1.627, and 2.809, respectively; when v2 increases to 0.166, the first bifurcation occurs, and the system turns into stable cycles of period 2, which means that each decision variables have two possible values; with the further increase of v2, the cold chain transportation levely, retail price of unit vaccine p. and inspection level s become chaotic; that is, each decision variables have a lot of possible values.Figure 5
Dynamic evolution of decision variables withv2. (a) v1=0.45,v3=0.32. (b) Lyapunov exponents for v2 varying from 0 to 0.5 of Figure 5(a). (c) The entropy diagram of Figure 5(a). (d) Chaotic attractor of the system.
(a)(b)(c)(d)Furthermore, assuming thatv1=0.19 and v2=0.36, Figure 6 analyzes the effect of v3 on system stability. When v3 changes at [0, 0.5], the changes of decision variables are shown in Figure 6(a). When v3∈0,0.151, the cold chain transportation level y, retail price of unit vaccine p, and inspection level s are stable at 2.492, 1.814, and 1.844, respectively; when v3 increases to 0.151, the system turns into stable cycles of period 2; with the increase of v3, eventually the cold chain transportation level y, retail price of unit vaccine p, and inspection level s become chaotic; that is, each decision’s variables have a lot of possible values.Figure 6
Dynamic evolution of decision variables withv3. (a) v1=0.19,v3=0.36. (b) Lyapunov exponents for v3 varying from 0 to 0.5 of Figure 6(a). (c) The entropy diagram of Figure 6(a). (d) Chaotic attractor of the system.
(a)(b)(c)(d)From Figures4(a)–6(a), it can be seen that as the adjustment speed of decision variables increases, the number of system cycles increases, and finally the chaotic state appears. The increasing of adjustment speed of decision variables means that there are irrational decisions in the vaccine supply chain. It is difficult to obtain the sustainable development of the vaccine supply chain structure.


## 5.2.2. Largest Lyapunov Exponent

The largest Lyapunov exponent can be used to determine whether the system is in a chaotic state. If the largest Lyapunov exponent is less than 0, the system is in a stable state; if the largest Lyapunov exponent is close to 0, the system is in a periodic motion; if the largest Lyapunov exponent is greater than 0, the system is in a chaotic state. Take Figure4(b) for example; when the largest Lyapunov exponent first returns to the 0 axis, the system reaches two times periodic bifurcation point. When the system returns for the second time and passes through the 0 axis, the system enters a chaotic state. Comparing the largest Lyapunov exponent with Figure 4(a), when bifurcation occurs and the system becomes in chaos, it can be seen that the largest Lyapunov exponents are consistent with the bifurcation diagram shown in Figure 4(a). Figures 5(b)–6(b) show similar results, which are not repeated here.


## 5.2.3. Entropy

When many scholars analyze the influence of time-delay parameter in the system, they only prove and give the conditions for the stability of the system. However, there are few people discussing what the chaotic behavior will bring to the system. The essence of entropy is the internal chaos of a system. According to information theory, the higher the entropy is, the higher the chaos of the system’s internal information is, which also results in the decision-maker making more effort to acquire useful information. From Figure4(c), we can know that when v1∈0,0.251, the entropy of model (6) is equal to 0, and the cold chain transportation levely, retail price of unit vaccine p, or inspection level s has a unique value. When v1>0.251, model (6) undergoes a doubly periodic bifurcation state and falls into chaos finally, the entropy of model (6) continues to increase, and the cold chain transportation level y, retail price of unit vaccine p, or inspection level s becomes unstable and has lots of values. From the information theory, the large entropy can increase uncertainty of information on decision variables, which may cause the distributor and retailer to miss the best decision-making opportunity to maintain vaccines’ activity in the process of vaccine transportation. Figures 5(c)–6(c) show similar results.


## 5.2.4. Chaotic Attractor

The chaotic attractors are shown in Figures4(d)–6(d) and are analyzed. When the system goes into a chaotic state, the structure of the chaotic attractor can be more complicated. Take Figure 4(d) for example; when v1>0.251, the system is in a chaotic state, and the chaotic attractor is complicated. Figures 5(d)–6(d) show similar results.


## 5.3. Influence of Adjustment Speed of Decision Variables on Stability Critical Point

Compared with Figure4(a), Figure 7(a) shows the change of system stability with v1 when v3 is constant and v2 is reduced to 0.3. Obviously, with the decrease of retailer’s adjustment speed v2, the critical point of stability moves to the right. Figure 7(b) shows the change of system stability with v1 when v2 is constant and v3 is reduced to 0.27. In the same way, with the decrease of retailer’s adjustment speed v3, the critical point of the stability moves to the right. It is not hard to see that, in terms of the goal of affecting the critical point of the stable state of the vaccine supply chain, there exist reverse changes in the relationship between the adjustment speed of the distributor’s decision variable v1 and the adjustment speed of the retailer’s decision variable v2 andv3. This interesting conclusion can imply powerful suggestions for distributor and retailer to find a balance in the process of vaccines transportation.Figure 7
Dynamic evolution of decision variables withv1. (a) v2 is reduced to 0.3. (b) v3 is reduced to 0.27.
(a)(b)


## 6. Chaos Control

In the vaccine supply chain especially, chaos can make the vaccines transportation difficult to be carried out smoothly, which inactivates the vaccines. In order to avoid chaos, we try to control the chaos by adjusting the parameter control method and variable feedback control method, respectively. The chaotic control effects of adjusting parameters on model (6) are demonstrated by numerical simulation. As can be seen from the foregoing, let τ=0.5,v1=0.4,v2=v3=0.5. When other parameters keep fixed, the system is in a chaotic state.


### 6.1. Adjustment Parameter Control Method

The control methods are mainly divided into feedback control method and nonfeedback control method. The adjusting parameter control method is one of the feedback control methods. This method is to change the internal structure of the original system to achieve control effect. The original system is as follows:(32)yt+1=yt+v1ytη1w−cα1−βpt−k1yt−τ,pt+1=pt+v2ptα1−2βpt+βwη1yt−τ+η2st−τ,st+1=st+v3stη2pt−wα1−βpt−k2st−τ.The system is transformed into the following one by parameter adjustment control [37]:(33)yt+1=1−uyt+v1ytη1w−cα1−βpt−k1yt−τ+uyt,pt+1=1−upt+v2ptα1−2βpt+βwη1yt−τ+η2st−τ+upt,st+1=1−ust+v3stη2pt−wα1−βpt−k2st−τ+ust.With the change of the added adjustment parameter, the system changes as shown in Figure8. When u=0, the system is in a chaotic state, which indicates that the distributor and the retailer cannot make joint decisions to control chaos. With the increase of u, the system reaches a stable state from chaotic state, which indicates that the distributor and retailer can effectively control the chaos by taking joint measures, such as signing contracts and so on.Figure 8
Systematic bifurcations with variation of adjustment parameteru.


### 6.2. State Variable Feedback Control Method

The state variable feedback control is a characteristic of modern control theory. State variable feedback control means that each state variable of the system is multiplied by the corresponding feedback coefficient and fed back to the input end and added to the reference input, and the sum is used as the control signal of the controlled system. The state variables of a system can show the internal characteristics of the whole system without knowing the internal structure of the system, and the control signals can be collected from the whole space or from the local. Therefore, compared with the traditional feedback control, the state variable feedback control can not only be a more excellent and effective control method, but also the form of its controller is more abundant [38].(34)yt+1=yt+v1ytη1w−cα1−βpt−k1yt−τ−uyt,pt+1=pt+v2ptα1−2βpt+βwη1yt−τ+η2st−τ−upt,st+1=st+v3stη2pt−wα1−βpt−k2st−τ−ust.Figure9 shows that the chaotic system gradually returns to a stable state through variable feedback control method. When u=0, the system is in a chaotic state, which indicates that the government does not take effective control measures for chaotic systems at this time. When u>0.253, the system returns to a stable state, which indicates that the government takes external intervention measures to accelerate the system to a stable state and ensure stable economic development.Figure 9
Systematic bifurcations with variation of adjustment parameteru.From the comparison between Figures8 and 9, it is obvious that the control system in Figure 9 enters the stable state earlier than in Figure 8, which means that the control effect of the variable feedback control method is better than that of adjusting parameter control method. The reason is that the structure of cold chain transportation in the vaccine supply chain is complex, and chaos control can bring about extra cost. It is difficult for the distributor and retailer to control the chaotic vaccine supply chain. At this time, it is imperative that external forces are used to control the chaos in the vaccine supply chain.


## 6.1. Adjustment Parameter Control Method

The control methods are mainly divided into feedback control method and nonfeedback control method. The adjusting parameter control method is one of the feedback control methods. This method is to change the internal structure of the original system to achieve control effect. The original system is as follows:(32)yt+1=yt+v1ytη1w−cα1−βpt−k1yt−τ,pt+1=pt+v2ptα1−2βpt+βwη1yt−τ+η2st−τ,st+1=st+v3stη2pt−wα1−βpt−k2st−τ.The system is transformed into the following one by parameter adjustment control [37]:(33)yt+1=1−uyt+v1ytη1w−cα1−βpt−k1yt−τ+uyt,pt+1=1−upt+v2ptα1−2βpt+βwη1yt−τ+η2st−τ+upt,st+1=1−ust+v3stη2pt−wα1−βpt−k2st−τ+ust.With the change of the added adjustment parameter, the system changes as shown in Figure8. When u=0, the system is in a chaotic state, which indicates that the distributor and the retailer cannot make joint decisions to control chaos. With the increase of u, the system reaches a stable state from chaotic state, which indicates that the distributor and retailer can effectively control the chaos by taking joint measures, such as signing contracts and so on.Figure 8
Systematic bifurcations with variation of adjustment parameteru.


## 6.2. State Variable Feedback Control Method

The state variable feedback control is a characteristic of modern control theory. State variable feedback control means that each state variable of the system is multiplied by the corresponding feedback coefficient and fed back to the input end and added to the reference input, and the sum is used as the control signal of the controlled system. The state variables of a system can show the internal characteristics of the whole system without knowing the internal structure of the system, and the control signals can be collected from the whole space or from the local. Therefore, compared with the traditional feedback control, the state variable feedback control can not only be a more excellent and effective control method, but also the form of its controller is more abundant [38].(34)yt+1=yt+v1ytη1w−cα1−βpt−k1yt−τ−uyt,pt+1=pt+v2ptα1−2βpt+βwη1yt−τ+η2st−τ−upt,st+1=st+v3stη2pt−wα1−βpt−k2st−τ−ust.Figure9 shows that the chaotic system gradually returns to a stable state through variable feedback control method. When u=0, the system is in a chaotic state, which indicates that the government does not take effective control measures for chaotic systems at this time. When u>0.253, the system returns to a stable state, which indicates that the government takes external intervention measures to accelerate the system to a stable state and ensure stable economic development.Figure 9
Systematic bifurcations with variation of adjustment parameteru.From the comparison between Figures8 and 9, it is obvious that the control system in Figure 9 enters the stable state earlier than in Figure 8, which means that the control effect of the variable feedback control method is better than that of adjusting parameter control method. The reason is that the structure of cold chain transportation in the vaccine supply chain is complex, and chaos control can bring about extra cost. It is difficult for the distributor and retailer to control the chaotic vaccine supply chain. At this time, it is imperative that external forces are used to control the chaos in the vaccine supply chain.


## 7. Conclusion

In this study, we focus on the vaccine supply chain composed of a distributor and retailer. We assume that the distributor chooses cold chain transportation to ensure vaccines activity, and the retailer checks whether the cold chain transportation effect that the distributor chose is qualified each time. We thought that the decisions of the distributor and retailer are not instantaneous and established a decision-making time-delay model. By using two-dimensional bifurcation diagram, Largest Lyapunov exponent, and entropy and chaotic attractors, we not only acquired the condition that the time-delay causes the system to produce a Hopf bifurcation, but also the condition that the adjustment speed of the decision variable causes the system to produce a period-doubling bifurcation. In addition, we used two different methods to control the chaotic system and compared the control effect. We drew the following conclusions.(1)
When time-delay parameterτ∈0,τ0, the vaccines transportation equilibrium is locally asymptotically stable. At this time, the distributor and retailer can cooperate well to ensure vaccines’ activity. When time-delay parameter τ≥τ0, the vaccine transportation system produces Hpof bifurcation and loses stability, which means that the vaccine supply chain changes from the vaccines transportation equilibrium to a limit cycle.(2)
If the distributor and retailer adjust the decision variables too quickly, the vaccine supply chain will bifurcate and fall into chaos and the entropy of the system will increase, which will cause the decision-maker to collect more additional information to find out the appropriate results, and it may make the decision-maker miss the best time to keep the vaccines active.(3)
In terms of the goal of affecting the critical point of the stable state of the vaccine supply chain, there exist reverse changes in the relationship between the adjustment speed of the distributor’s decision variablev1 and the adjustment speed of the retailer's decision variable v2 (or v3).(4)
Compared with the internal joint control of the distributor and retailer, external control, such as government intervention, will have a better control effect of the chaos of the system.In order to simplify the problem, we established a simplified vaccine supply chain model, while ignoring many important features of the vaccine supply chain. For example, there may be some loss of vaccines during transportation. A supply chain considering the wastage rate of vaccines will make the model closer to reality. Also, vaccines are different from ordinary profitable products, and a supply chain considering social responsibility is also a research hotspot. We will study these issues in the future.



---

*Source: 4392708-2020-10-06.xml*