# Neimark-Sacker Bifurcation in Demand-Inventory Model with Stock-Level-Dependent Demand

Authors: Piotr Hachuła; Magdalena Nockowska-Rosiak; Ewa Schmeidel
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

An analysis of dynamics of demand-inventory model with stock-level-dependent demand formulated as a three-dimensional system of difference equations with four parameters is considered. By reducing the model to the planar system with five parameters, an analysis of one-parameter bifurcation of equilibrium points is presented. By the analytical method, we prove that nondegeneracy conditions for the existence of Neimark-Sacker bifurcation for the planar system are fulfilled. To check the sign of the first Lyapunov coefficient of Neimark-Sacker bifurcation, we use numerical simulations. We give phase portraits of the planar system to confirm the previous analytical results and show new interesting complex dynamical behaviours emerging in it. Finally, the economical interpretation of the system is given.


---

## Body

## 1. Introduction

Great economic development after the second world war has released a need of development of mathematical methods to support optimization of economic and business processes. Material flow, production, and inventory are aspects of a business, which, to make it profitable, need to be optimized. Therefore, many models of supply chain were created in the mid-20th century. To mention the most noticeable, we ought to list works of Wagner and Whitin [1], Brown [2], and Holt et al. (HMMS model) [3]. They have laid the foundation for supply chain modelling. Those models, despite being relatively simple, have become an inspiration for contemporary researchers: to redefine models in order to fit to the current challenges and to analyse them using available computation power, for example, [4–6].Prediction of future demand and inventory is an important aspect of running and managing manufacturing or trade company. Methods supporting those tasks have been developed by economists already in the mid-20th century; nonetheless, they are still being improved, as economy still changes and creates new challenges. One of those methods is modelling of economic phenomena using mathematical formulations. The topic of demand and inventory needs a contextual approach, since many factors can influence it and different views may be needed. Therefore, the models are created with the usage of different mathematical tools. We can mention here recent works related to demand and inventory that investigate and describe specific economic cases: Chen and Hu in [7] consider an inventory and pricing model and develop polynomial time algorithms to maximize the total profit; Mondal et al. in [8] consider generic algorithm for the case of inventory of a deteriorating item; Qi et al. in [4] consider supply chain that experiences a disruption in demand during the planning horizon.Ma and Feng in [6] proposed the dynamical model of demand and inventory with mechanism of demand stimulation and inventory limitation. Nowakowska in [9] showed application of such a model to electrical energy market. Hachuła et al. analysed the stability of equilibrium points in [10, 11].The model describes demand and inventory of a product at one echelon of supply chain at the retailer. The considered supply chain consists of three echelons: manufacturer, retailer, and customers. The following rules are applied to the model: customers buy a good from a retailer; a retailer orders a product in the forecasted amount and forecast is prepared using single exponential smoothing Brown model [12]; a manufacturer produces and delivers exactly the ordered amount and the production capacity is unlimited; customers’ demand depends on a retail price, which can be changed by a discount; price cannot be arbitrarily changed but the retailer can offer a discount depending on stock volume.The model takes a form of the following system of difference equations:(1)Dn+1=ATA+1T-SnkDn,Sn+1=Sn-Dn+Dˇn,Dˇn+1=αDn+1-αDˇn,where n=0,1,2,… indicates instance of time; Sn is a stock volume at n; Sn is a real number; interpretation of nonnegative values is obvious and the negative value of stock informs what amount is missing to fulfill demand and needs to be overcompensated in the next periods; Dn is a demand volume at n, which assumes constant demand elasticity (see [13, p. 280]), Dn≥0; Dˇn is a forecast of demand at n and order placed by a retailer at a manufacturer; moreover, by assumption of unlimited capacity, it is also the delivered quantity at n, Dˇn≥0; A is a parameter for discount steering, A>0; T is a parameter for defining the target stock, T>0; k is a price elasticity coefficient that regulates dependence between price, discount, and demand, k>0; α is a forecast smoothing coefficient of Brown model, α∈(0,1).The first equation describes the influence of relation between current stockSn and target stock T on future demand Dn+1, in a way that when stock is high, the retailer offers a discount to encourage customers to buy a product. The second equation lets us calculate future stock Sn+1 as a sum of current stock Sn and delivery of amount of forecasted demand Dˇn diminished by actual sales Dn. The third equation describes a method of demand forecasting. Future demand is forecasted using single exponential smoothing by Brown [2]. Smoothing coefficient α split importance between actual demand (sales) Dn and forecasted demand Dn.We analyse properties of a given model (1) on the ground of theory of difference equations. Three-dimensional difference systems were studied by many authors, for example, in [14–17]. Under some assumptions, such systems can be rewritten as third-order difference equations. Asymptotic properties of these types of equations were investigated in [18–23]. For the background of difference equations theory, see monographs [24, 25]. Firstly, we construct a positive invariant set of the three-dimensional system. We notice that the given system has an invariant plane and can be reduced to the family of the planar systems, which we construct positive invariant sets and provide the equilibrium points for. We recall the properties of equilibrium points obtained in [11]. The aim of this paper is to analyse the dynamics of planar systems on the ground of bifurcation theory (see [26]) and numerical simulations. It was shown that nondegeneracy conditions for Neimark-Sacker bifurcation are fulfilled and existence of transcritical bifurcation is suggested. The analytical results are confirmed by numerical simulations. Finally, the economical interpretation of the system is given.


## 2. Main Results

For the sake of convenience, let us write the model as a system of three first-order difference equations withx, y, and z replacing D, S, and Dˇ, respectively, in system (1),(2)xn+1=ATA+1T-ynkxn,yn+1=yn-xn+zn,zn+1=αxn+1-αzn,with xn,zn⩾0, yn<(A+1)T, or as (xn+1,yn+1,zn+1)=f(xn,yn,zn) with a mapping:(3)fx,y,z=ATkA+1T-ykx,y-x+z,αx+1-αz.It is easy to see that (0,y,0) for y<(A+1)T and (x,T,x) for x≥0 are equilibrium points of the system (xn+1,yn+1,zn+1)=f(xn,yn,zn).Let us shortly elaborate on the meaning of those equilibrium points. The goal of the retailer is to reach target stockT. If additionally actual demand xn and forecasted demand zn are equal and greater than zero, the retailer’s stock level does not change in the next period, which means that system (2) is in the equilibrium point (x,T,x). On the other hand, when xn=zn=0, the retailer’s stock level y cannot be changed. This means that the retailer remains in equilibrium point (0,y,0).We recall thatBINVf is the positive invariant set of the system (xn+1,yn+1,zn+1)=f(xn,yn,zn) if BINVf⊂Df and f(BINVf)⊂BINVf.Proposition 1 (see [11]).
Here,(4)BINVf=x,y,z∈R3:x,z≥0,y+α-1z<A+1T∪x,T,x:x≥αAT is the positive invariant set of the system (xn+1,yn+1,zn+1)=f(xn,yn,zn), where f is given by (3).One can observe that the second and the third equation of system (2) are linear ones. Therefore, we treat them with some linear transformations in order to obtain a more convenient form for analysis: both sides of the second equation multiplied by α and the second and the third equation added. As an output, we get(5)αyn+1-yn=zn-zn+1,which enables us to state the following.Corollary 2.
System (2) has the invariant plane αy+z=αy1+z1 for any y1∈R,z1≥0,αy1+z1<α(A+1)T. Hence, system (2) can be reduced to a family of planar systems:(6)xn+1=ATA+1T-ynkxn,yn+1=1-αyn-xn+αs,with s<(A+1)T, where s≔y1+z1/α, y1∈R,z1≥0.Proof.
From (5), we get the first-order recursion of the form αyn+1+zn+1=αyn+zn which implies that αyn+zn=αy1+z1. Constraint of the expression αy1+z1 is a direct consequence of the positive invariant set construction.Let us later write our systems (6) in the form of a family of mappings:(7)gsx,y=ATA+1T-ykx,1-αy-x+αs,for s<(A+1)T. Taking into consideration economic interpretation, gs:Dgs→R2, where Dgs≔[0,∞)×(-∞,(A+1)T).Proposition 3.
LetA,T,k>0, α∈(0,1), and s<(A+1)T. The set BINVgs=Dgs is the positive invariant set of the system (xn+1,yn+1)=gs(xn,yn), where gs is given by (7).Proof.
Lets<(A+1)T. For any x≥0, y<(A+1)T, we have (gs)1(x,y)≥0 and (gs)2(x,y)=(1-α)y-x+αs<(1-α)(A+1)T+α(A+1)T=(A+1)T.Mapping (7) may have at most two equilibrium points:(i)
(0,s)(ii)
(α(s-T),T) iff s>TStability of the equilibrium points of (2) has been examined in [11]. Based on the findings of [11], one can formulate the conclusions about the existence and stability of equilibrium points of planar systems (6), which are presented in Table 1. Now, using these results, we analyse the dynamics of the planar systems using bifurcation theory. We are especially interested in proving the existence of Neimark-Sacker bifurcation.Table 1
Overview of conditions for existence and stability of the equilibrium points of (6).
Conditions on parameters ( 0 , s ) ( α ( s - T ) , T ) C1 T > s Stable Does not exist C2 T = s Unstable Does not exist C3 T < s < ( A + 1 ) T,  k≤1 Unstable Stable C4 T < s < ( 1 + A / k ) T,  k>1 Unstable Stable C5 ( 1 + A / k ) T < s < ( 1 + A ) T,  k>1 Unstable UnstableGraphs illustrating the cases in Table1 are presented in Figures 1–4.Figure 1
Examples illustrating conditions (C1) and (C2) from Table1.Figure 2
Example illustrating condition (C2) from Table1.Figure 3
Examples illustrating conditions (C3) and (C4) from Table1.Figure 4
Examples illustrating condition (C5) from Table1. The right figure is a “big picture” of the same case.The planar system is dependent on five parametersα, A, T, k, and s. The most interesting case in the bifurcation analysis for the planar system (6) is dependence on parameter s with fixed A0,T0,k0>0, α0∈(0,1), which now we will carry on.From the conditions in Table1, in order to examine bifurcations on parameter s, we have to consider two cases dependent on k0.Case 1 (k0∈(0,1]).
Fors<T0, the equilibrium point (0,s) is the only equilibrium point of the system (xn+1,yn+1)=gs(xn,yn) which is stable. For s>T0, the considered system possesses two equilibrium points: unstable (0,s) and stable (α0(s-T0),T0). This suggests that s=T0 is the bifurcation parameter for the transcritical bifurcation.Case 2 (k0>1).
Fors=T0, probably the transcritical bifurcation occurs like in the previous case. For s<(1+A0/k0)T0, the equilibrium point (α0(s-T0),T0) is stable and loses its stability for s>(1+A0/k0)T0. We prove that, for s=(1+A0/k0)T0, the Neimark-Sacker bifurcation occurs. We use classical analytical results from [26] and numerical simulations. For the sake of convenience, let us quote from [26] Theorem 4.5 as Lemma 4 and Theorem 4.6 as Lemma 5.Lemma 4 (Theorem 4.5, [26]).
Suppose that a two-dimensional discrete-time systemx→f(x,α), x∈R2, α∈R1, with smooth f, has, for all sufficiently small |α|, the fixed point x=0 with multipliers μ1,2(α)=r(α)e±iϕ(α), where r(0)=1, ϕ(0)=θ0. Let the following conditions be satisfied:
(C.1)
r ′ ( 0 ) ≠ 0. (C.2)
e i k θ 0 ≠ 1 for k=1,2,3,4.
Then, there are smooth invertible coordinate and parameter changes transforming the system into(8)y1y2⟶1+βcos⁡θβ-sin⁡θβsin⁡θβcos⁡θβy1y2+y12+y22cos⁡θβ-sin⁡θβsin⁡θβcos⁡θβaβ-bβbβaβy1y2+Oy4,with θ(0)=θ0 and a(0)=Re⁡(e-iθ0c1(0)), where c1(0) is given by formula (4.21) in [26].Lemma 5 ((Theorem 4.6, [26]) (generic Neimark-Sacker bifurcation)).
For any generic two-dimensional one-parameter systemx→f(x,α), having at α=0 the fixed point x0=0 with complex multipliers μ1,2(α)=e±iθ0, there is a neighbourhood of x0 in which a unique closed invariant curve bifurcates from x0 as α passes through zero.Remark 6 (see [26]).
The genericity conditions assumed in the theorem are the transversality condition (C.1) and the nondegeneracy condition (C.2) from Lemma4 and the additional nondegeneracy condition:(C.3)a(0)≠0.We start with translation mappinggs to the mapping which possesses fixed point in (0,0) for sufficiently small s. We consider the C2 mapping Gs:R×(-∞,A0T0)→R2 with parameter s∈R:(9)Gsx,y≔A0T0k0A0T0-yk0x+α0s+A0T0k0-α0s+A0T0k0,1-α0y-x.The Jacobi matrix of Gs at (0,0) has the form (10)JGs0,0=1α0k0A0T0s+α0-11-α0, with the complex eigenvalues for sufficiently small s of the form(11)λ1,2s=1-12α0±i24α0k0A0T0s+1-α02,and modulus (12)rs=1+α0k0A0T0s. Notice that r(0)=1, r′(0)=α0k0/2A0T0≠0, which means that condition (C.1) of Lemma 4 is satisfied. Now, we check that condition (C.2) (called no strong resonance) is satisfied. We have to prove that eilθ0≠1, where l=1,2,3,4 and θ0 is such that λ1,2(0)=r(0)e±iθ0. We get for α0∈(0,1)(13)1-12α0+i24α0-α02=eiθ0≠1,1+α012α0-2+i1-12α04α0-α02=e2iθ0≠1,1-12α0α0-32+i2α02-4α0+34α0-α02=e3iθ0≠1,1+12α0α03-8α02+20α0-16+i-12α03+3α02-5α0+24α0-α02=e4iθ0≠1, because 2,2±2 are solutions to equation -1/2α03+3α02-5α0+2=0. Fulfilling conditions (C.1) and (C.2) from Lemma 4 means that nondegeneracy conditions are satisfied. To get the existence of Neimark-Sacker bifurcation, we have to check the last condition (C.3) of the form a(0)=Re⁡[eiθ0c(0)]≠0, named the first Lyapunov coefficient, given in Lemma 5. If a(0)<0, we get the supercritical Neimark-Sacker bifurcation, and if a(0)>0 we get the subcritical Neimark-Sacker bifurcation. Due to complicated conditions for any parameters A0,T0>0, α0∈(0,1), and k0>1, we check the sign of the first Lyapunov coefficient only for specific A0, T0, α0, and k0 in Section 3.Remark 7.
In the analogous way, we can prove the existence of the Neimark-Sacker one-parameter bifurcation of the planar system for any other parametersA, T, k, and α with the bifurcation parameter which fulfills the equality s=(1+A/k)T.


## 3. Numerical Simulations

In this section, by using numeral simulation, we give the bifurcation diagrams and phase portraits of system (6) to confirm the previous analytical results and show new interesting complex dynamical behaviours emerging for system (6).Now, let us give illustration to the cases provided in Table1. Graphs in Figures 1–4 depict behaviour of (6) depending on the parameter s. The other parameters remain unchanged at A=1, T=1000, and α=0.5, and, for C1, C2, C4, and C5, k=4; for (C3), k=0.5. The graphs clearly confirm the analytical results.Interesting behaviour can be observed at border points=T depicted in Figures 1-2 in Graphs (C2a-b). At the first glance, in Graph (C2a), one can assume stability of the equilibrium point (0,s), which would be a contradiction of the analytical result. However, looking from a closer perspective in Graph (C2b) one can see that the orbits leave the surroundings of the equilibrium point, herewith confirming the instability proved with analytical methods.Let us have again a look at the border case ofs=(1+A/k)T depicted in Figure 5. In each of the graphs, the same parameters are applied in order to observe the behaviour of the system for varying initial conditions of x and y (A=1, T=1000, α=0.5, s=1250, and k=4). Each time, we set the initial conditions closer and closer to the equilibrium point (x∗∗,y∗∗). Appearance of this kind of circle suggests supercritical bifurcation of Neimark-Sacker type. Indeed, negative sign of the first Lyapunov coefficient a(0) obtained numerically confirms that suggestion.Figure 5
Dynamics of the cases=(1+A/k)T with increasing zoom on the equilibrium point (x∗∗,y∗∗).
(a) (b) (c) (d)


## 4. Economical Interpretation

System (6) depends on five parameters: A, α, T, k, and s. Each of the parameters carries the economical meaning described in Section 1. In a business practice, the parameters T and s are set in advance. Nonetheless, as it is shown in Table 1, changes of parameter s cause the change of behaviour of the equilibrium points. Therefore, we present bifurcation diagrams of x and y as s changes in Figures 6 and 7.Figure 6
Diagram of bifurcation of the equilibrium points with respect tox, for the parameter s varying from 800 to 2000, x1=0.1, y1=s for A=1, T=1000, α=0.5, and k=4.Figure 7
Diagram of bifurcation of the equilibrium points with respect toy, for the parameter s varying from 800 to 2000, x1=0.1, y1=s for A=1, T=1000, α=0.5, and k=4.Parameterk, which represents price elasticity, is defined depending on the types of a product and a market and often is determined experimentally and must be adjusted frequently. Therefore, it is interesting to analyse behaviour of system (6) at varying values of this parameter. For that purpose, let us illustrate changes of the dynamics of the equilibrium points in the bifurcation diagrams (see Figures 8 and 9). Model (6) includes steering mechanism, which allows the user to influence the behaviour of the system. With the parameter α, one can decide on the share of influence between the forecast and actuals of sales. The less α is, the greater the impact of actuals is, and vice versa. Figure 11 presents bifurcation diagram for changing parameter α. Similarly, in Figure 10, we present bifurcation for varying A. A is a parameter to control the influence of the discrepancy between the actual stock and the target one T.Figure 8
Diagrams of bifurcation of the equilibrium points with respect tox for the parameter k varying from 4.95 to 6.15, x1=100.1, y1=1000 for A=1, T=1000, α=0.5, and s=1200.Figure 9
Diagrams of bifurcation of the equilibrium points with respect toy for the parameter k varying from 4.95 to 5.85, x1=100.1, y1=1000 for A=1, T=1000, α=0.5, and s=1200.Figure 10
Diagrams of bifurcation of the equilibrium points with respect tox (a) and y (b), for the parameter A varying from 0 to 1.2, x1=0.1, y1=s for s=1250, T=1000, α=0.5, and k=4.
(a) (b)Figure 11
Diagram of bifurcation of the equilibrium points with respect tox (a) and y (b), for the parameter α varying from 0 to 0.3, x1=100.1, y1=s=1230 for A=1, T=1000, and k=4.
(a) (b)Business Example. Let us consider a product HR that follows demand(D)-price(p) dependence D2/D1=(p2/p1)-k with price elasticity k=4. The target T=1000 is set on stock volume and assumed demand at least z1=50. Actual demand is insufficient x1=1, and hence the actual stock exceeds the target y1=1100>T. Smoothing coefficient α=0.5, which means calculation of the amount to be delivered takes into account actuals and forecast with equal wages. Parameter A=1. One can easily calculate s=y1+z1/α=1200. From Table 1, it is observed that the product satisfies condition (C4), and hence the system has two equilibrium points (0,1200) and (100,1000), unstable and stable, respectively.



---

*Source: 8162865-2017-01-01.xml*