MODELING AND ANALYZING THE CHAOTIC BEHAVIOR IN SUPPLY CHAIN NETWORKS: A CONTROL THEORETIC APPROACH

. Supply chain network (SCN) is a complex nonlinear system and may have a chaotic behavior. This network involves multiple entities that cooperate to satisfy customers demand and control network inventory. The policy of each entity in demand forecast and inventory control, and constraints and uncertainties of demand and supply (or production) signiﬁcantly aﬀects the complexity of its behavior. In this paper, a supply chain network is investigated that has two ordering policies: smooth ordering policy and a new policy that is designed based on proportional-derivative controller. Two forecast methods are used in the network: moving average (MA) forecast and exponential smoothing (ES) forecast. The supply capacity of each entity is constrained. The eﬀect of demand elasticity, which is the result of marketing activities, is involved in the SCN. The inventory adjustment parameter and demand elasticity are the most important decision parameters in the SCN. Overall, four scenarios are designed for modeling and analyzing the chaotic behavior of the network and in each scenario the maximum Lyapunov exponent is calculated and drawn. Finally, the best scenario for decision-making is obtained.


1.
Introduction. Supply chain network (SCN) consists of multiple entities. The operations of each entity are receiving goods, satisfying demands, controlling the inventory, and placing an order. These operations need interaction between the entities of the SCN. There are one or two inventory feedback control loops in each entity. There is also a time delay between order placing and goods delivery. Various types of uncertainties exist in the network, such as demand uncertainty, production uncertainty, and delivery uncertainty. These characteristics turn the supply chain into a complex and nonlinear dynamic system that may have a chaotic behavior.
One type of dynamic behavior is caused by marketing and competition activities that create interaction between suppliers and customers. The interaction may generate a chaotic behavior in the supply chain [17,22]. A price change has a fundamental effect on customers demand. Usually the demand goes down as price increases and vice versa. Wu and Zhang [35] showed the chaotic behavior of the supply chain by simulating the interaction between customers and suppliers where customers respond to the price discount offer made by the supplier and the supplier adjusts the price according to stock held.
Another type of dynamic behavior is a phenomenon in which the demand fluctuation is amplified along the upstream of the supply chain. The phenomenon is referred to as the bullwhip effect or Forrester effect [8]. The bullwhip effect is the culprit behind chaos in inventory levels and production strategies [1]. Four major factors of this effect are: demand forecast updating, reorder batching, price fluctuating, and shortage gaming [20]. The negative impacts of the bullwhip effect can be summed up as follows [13]: (1) inaccurate forecasting which may cause periods of low capacity utilization or periods of having no enough capacity; (2) inadequate customer service; (3) high inventory costs. Dejoncheere et al. [6,7] changed the supply chain to several linear subsystems to study the bullwhip effect and, by obtaining transfer function and drawing frequency response curves, measured this effect and investigated the way of avoiding it. Similar studies have been carried out by Daganzo [4,5]. In recent years, a variety of strategies have been studied to reduce or avoid the bullwhip effect [18,3,26].
Hwarng and Xie [14] introduced five main factors that influence the supply chain, namely, demand pattern, ordering policy, demand-information sharing, lead time, and supply chain level. They showed the chaotic behavior of the supply chain and its sensitivity to small changes of inventory control parameters using the beer distribution model [16] and Sterman dynamic equations [30]. They also quantified the degree of chaos using the maximum Lyapunov exponent (LE) across all level of the supply chain.
Previous studies have generally focused on a serial supply chain with four successive levels-factory, distributor, wholesaler, and retailer-where each level has only one entity. But the fact is that supply chain is a vast network with many entities in each level. Recently, new studies have been conducted on the chaotic behavior of the SCN. Ouyang and Li [24] analyzed the bullwhip effect of supply chains with a general network topology and general linear ordering policy that would allow for modeling complex interactions among suppliers at different levels. Tarokh et al. [31] investigated the role of supply network configuration in creating chaotic behavior. Their proposed model consisted of more than one entity at each network level with two agents between each two successive level for receiving and dispatching orders. Depending on its functionality, agents were called Orders Receiving-Dispatching (ORD) agent and Products Receiving-Dispatching (PRD) agent.
The control or reduction of chaotic behavior in the SCN is an important issue in supply chain management (SCM). This behavior can lead to inventory instability in the entities and disrupt decision-making. In the chaotic SCN, when the market demand increases, suppliers are often unable to support manufacturers; while market demand is slowdown, the suppliers often continue to overproduction, which resulting in overstock. Spiegler et al. [28] developed more accurate simplified linear representations of complex nonlinear supply chain models and reduced the complexity. The synchronization of two identical chaotic supply chain management systems based on their mathematical model is presented in [11]. Wu et al. [36] proposes a method for analyzing the operational complexity in supply chains by using an entropic measure based on information theory.
In this paper, a supply chain network is provided that has four successive levels based on the beer distribution model. Each level has various entities. Each entity must satisfy demand, control inventory and place an order through interactions with adjacent entities. The production constraint and supply constraint are facts that have been overlooked in previous studies and are included in the model. Therefore in this model, production capacity is limited and each entity has a supply capacity at any given time based on its inventory and can satisfy demand to the same quantity. The effect of marketing activities on the SCN is modeled with a nonlinear function commonly used in economics and the interaction between the supply network and customers and the effect of demand elasticity are examined. Two methods are used to forecast demand: moving average (MA) forecast method and exponential smoothing (ES) forecast method. In this model, each entity has a feedback control loop of inventory and two ordering policies are used: smooth ordering policy and a new policy that has been designed based on the proportional-derivative (PD) controller [2]. In short, four scenarios are designed for examining the chaotic behavior of the SCN: smooth ordering policy and ES forecast method (Scenario 1), smooth ordering policy and MA forecast method (Scenario 2), PD controller-based ordering policy and ES forecast method (Scenario 3), and PD controller-based ordering policy and MA forecast method (Scenario 4). Graphical methods such as time series and phase plots are used to show the existence of chaos. In each scenario, by calculating the maximum LE, the effects of the inventory adjustment parameter and the elasticity of demand on the SCN behavior are investigated and the best scenario (Scenario 4) is selected.
The superiority of this study to similar studies is more comprehensive modeling and improving behavior of the SCN in chaotic conditions. The supply chain is modeled as a network with many entities in each level. The production constraint and supply constraint are included in the model. The relationship between market activities and SCN using the elasticity of demand is a new topic. It is shown that the new ordering policy based on PD controller has good performance in reducing the chaotic behavior in the SCN. This is very important for decision makers to minimize inventory costs. This paper is organized as follows: section 2 briefly introduces chaotic systems and the feedback control theory. In section 3, a model of the SCN is described and its dynamic equations are obtained. The block diagram of entity operations is provided. Demand forecasting methods and ordering policies and the effect of demand elasticity are investigated. In section 4, a ten-entity SCN is simulated with four different scenarios and their results are compared. The paper ends with a conclusion. cooperation.

Preliminaries.
2.1. Chaotic systems. Chaotic systems are deterministic systems with high complexity that show irregular behavior and belong to the category of nonlinear dynamic systems. Chaotic waveforms consist of patterns that are repeated at nonperiodic intervals and are not asymptotically stable. Responses of these systems are extremely sensitive to changes of initial conditions and parameters. This characteristic is one of the main features of chaotic systems and is used as an important indicator for identifying them. Chaotic systems are often divided into two categories: continuous-time and discrete-time. Discrete-time systems must be at least first-order and continuous-time systems must be at least third-order.
The features of chaotic systems can be summed up as follows [34,9]: (1) nonlinearity and non-randomness; (2) high sensitivity to initial conditions; (3) nonperiodic or non-quasi-periodic; (4) strange attractors; (5) impossibility of the prediction of long-term behavior of the system; (6) chaos as inherent feature of these systems not because of foreign factors such as noise; (7) non invertible: impossibility of accessing the last of chaotic systems.
Usually there are two ways to identify and measure chaos: graphical methods [33] and quantitative methods [33,29]. Using these methods, it can be revealed whether the systems behavior is stable, periodic, quasi-periodic or chaotic.
Graphical methods are visible but less accurate, while in quantitative methods the degree of chaos can be computed. Various types of plots such as time series and phase plots are used to study the behavior of dynamic systems. Time series of a stable system moves toward a fixed point (equilibrium point) which is its basin of attraction. Stable limit cycle is a behavior of the dynamic system that its time series is periodic and basin of attraction is a circle. Quasi-periodic systems have quasiperiodic time series and their basin of attraction is a torus. In chaotic systems, the time series is very sensitive to changes in initial conditions and parameters. The basin of attraction is a fractal structure and state space trajectories become divergent and convergent continuously in a finite state space. These trajectories are not repeated and are thus unpredictable.
The Lyapunov exponent (LE) is the most important quantitative method that measures the sensitivity of initial conditions and is a standard quantifier for determining and classifying the behavior of nonlinear systems. It is the average rate of convergence and divergence in near state routes and usually shown by the symbol λ. In practice, a wide range of LEs can be obtained, but what matters is the largest LE which is calculated as follows [29]: ∆R n and ∆R n+1 are the distances between two nearby trajectories at times n and n + 1 respectively. The number of LEs also increases with the systems dimension. In chaotic systems, at least one LE or the largest LE is positive. .If all LEs are negative, the system will be stable. If one LE is zero and the others are negative, the limit cycle will be stable and if two LEs are zero and the others are negative, there will be a stable torus.

2.2.
Feedback control theory. A typical control system has the feature that some output quantity is measured and then compared with a desired value, and the resulting error is used to correct the system's output. This concept is called feedback or closed-loop control [10]. A schematic of a control system is shown in Fig. 1. The controller compares the measured process variable y(t) with a desired variable r(t) that is often called the set point. The difference or error signal, e(t), is then processed to calculate a control signal u(t) as a new process input. This input will try to adjust the measured process variable back to the desired set point. Usually controllers can be classified as follows: (1) (4) proportional-integral-derivative (PID) controller. These controllers are often combined with logic, sequential functions, selectors, and simple function blocks to build the complicated automation systems used for energy production, transportation, and manufacturing [2]. A mathematical description of the PID controller in the continuous time domain is: The controller parameters are proportional gain K p , integral or reset time T i , and derivative time T d . The control signa u(t) is a sum of three terms: The proportional term (P), the integral term (I), and the derivative term (D). The integral, proportional and derivative part can be interpreted as control actions based on the past, the present, and the future.
The proportional term (P) gives a system control input proportional with the error. There is always a steady state error in proportional (P) control. The error will decrease with increasing gain, but the tendency towards oscillation will also increase. Using a too large P term gives an unstable system.
The proportional term (P) gives a system control input proportional with the error. There is always a steady state error in proportional (P) control. The error will decrease with increasing gain, but the tendency towards oscillation will also increase. Using a too large P term gives an unstable system.
The integral term (I) gives an addition from the sum of the previous errors to the system control input. The main reason for integral (I) control is to reduce or eliminate steady-state errors. The strength of integral action increases with decreasing integral time T i . The tendency for oscillation also increases with decreasing T i .The most common use of the I term is normally together with the P term, called a PI controller.
The derivative term (D) gives an addition from the rate of change in the error to the system control input. A rapid change in the error will give an addition to the system control input. This improves the response to a sudden change in the system state or desired variable. The D term is typically used in conjunction with the P or PI as a PD or PID controller. Using these controllers, damping increases and the stability of a system improves.
To derive the discrete form of (2), integral and derivative terms are approximated as where T s is a sampling time.
As a result, the PID controller in the discrete time domain is derived as follows: 3. Model.
3.1. Model specification. In this paper, a supply chain network is investigated which, like the beer distribution model, has four successive levels: factories, distributors, wholesalers, and retailers (Fig.2). Each level has several entities. There is information sharing between the entities of a level and between two successive levels. In this system, orders propagate from customers to factories and products flow from factories to customers. The products are distributed among the entities of a downstream level according to their actual supply lines (orders placed but not received). In each level, received orders from the downstream level and incoming products from the upstream level are distributed among its entities. Each entity in the SCN has a feedback control loop of inventory and can make decisions. Decisions are made locally based on information, such as the actual and desired inventory levels, backlog, expected orders, and incoming products. The entire ordering process is expressed as follows: at each time t , the entity receives products, satisfies demand, observes the new inventory level, and finally places an order. There is a time delay between order placement and delivery. Customers demand is randomly distributed among retailers. In upper levels, the demand of each entity is equal to a fraction of lower levels orders.

Factories
Wholesalers Retailers Customers Distributers In this model, the production capacity is constrained; as well as each entity has a supply capacity at each time t based on its inventory and can satisfy demand to the same quantity. Basic hypotheses of the model are: 1) An entity with sufficient inventory must fulfill orders.
2) Unfilled orders are kept in a backlog and will be filled when the inventory is sufficient. 3) Placed orders cannot be cancelled and incoming products cannot be returned.

Dynamic equations.
Each entity in the SCN receives some products after a time delay from the time of placing an order. At the same time, a new demand is received. Based on their supply capacity, entities fulfill all or part of the backlog and current demand. Operations of each entity are as follows: where x i (t) is the effective inventory (inventory level after fulfilling the backlog), is the incoming product , a part of order was made in the time t − τ in the past, and b i (t − 1) is the backlog at time t − 1.
x i (t) may be positive or negative (due to backlog), but It should be noted that an actual SCN is a nonlinear discrete time-delay system. Its capacity is limited and may be an inventory shortage. These features are included in this model.
There are N r retailers, N w wholesalers, N d distributors, and N f factories, where the number of entities at the upper level is less than that of the lower level: In a real SCN, there is no information sharing between customers and entities at the retailer level. Hence, customers demand is randomly distributed among retailers. If d(t) is the total customers demand, the demand of each entity at the retailer level is: where w i is a random number and: At the other levels, the demand of each entity is as follows: where N d is the number of entities and o d j (t) is the order of each entity at a downstream level. At upstream levels, each entity can accept any fraction of downstream orders and the equation (9) is established.
The supply capacity c i (t) and the backlog b i (t) for each entity are calculated as follows: Since there is information sharing between two successive levels, the incoming product f d j (t) at the entities of a downstream level is: where N u is the number of entities at the upstream level and k j is the product distribution coefficient which is obtained according to the actual supply line z j (t)(orders placed but not received) that can be defined as follows: Therefore, k j in a downstream level is given: Fig. 3 shows the block diagram of entity operations. Each entity receives new demand from lower level and previous order from upper level. At the same time, it determines the supply capacity (11) and backlog (12) and then updates their inventory and output according to the equations (5) and (6).o l j (t) (the most important decision variable in the SCN ) is the order quantity of j th entity in the lower level that is obtained with different ordering policies. This paper examines two ordering policies. A well-known one is the smooth ordering policy whose decision equation is as follows: where α i is the inventory adjustment parameter,d i (t) is the demand forecast, and e i (t) is the error between the actual inventory x i (t) and the desired inventoryx i (t): It should be noted that if x i (t) is greater thanx i (t), then α i e i (t) is negative. It may alsod i (t) is less than |α i e i (t)| and thusd i (t) + α i e i (t) will be negative and the order will not be done. The demand forecastd i (t) usually obtained from two forecast methods: 1) Moving average (MA) forecast method: whereT m is the number of periods used to compute the forecast; 2) Exponential smoothing (ES) forecast method: θ i is a parameter which determines how fast expectation are updated. Usually the following equation is used to make a relative comparison of these forecasting methods [7]: A new ordering policy based on the PD controller is used in the model: τ D i is called the derivative time constant which is adjustable. In this policy, the difference between two successive errors is used to accelerate decision-making, It should be noted that the production of a factory is constrained by c p (production capacity): In the retailer level, customers demand d(t) is placed instead of N l j=1 o l j (t). Nonlinear demand function is one of the most important functions in economics that states the relationship between demand and price: where p is the price, d max is the potential market demand (the maximum possible demand), and β is the price elasticity of demand which is defined as: If customers demand decreases, there will be an overstock at the retailers level and its effect will transfer toward upstream in the SCN. Consequently, there will be an overstock at the factories level and there will be an inevitable need to reduce the price. If the retailers inventory is: (25) and x T denotes the threshold stock in the retailers level, the overstock rate is defined as: Assuming that the constant c is to convert the overstock rate into the discount rate r(t) and considering a maximum discount rate r max , the discount rate is calculated as follows: If p o and d o are the price and the basic demand respectively and p(t) and d(t) equal the price and the new demand at time t, we have: By substituting (28) and (29) in (30) and through a series of computations, the relationship between the new demand and the basic demand is obtained as follows: which expresses the relationship between customers demand , the inventory of SCN and the elasticity of demand.
In short, the SCN activities can be summed up as follows: 1) Receiving a demand in each level: customers demand d(t) in the retailer level and total order of low-level entities N l j=1 o l j (t) in the other levels. 2) Determine the demand of each entity according to (8) and (10).
3) Receiving a part of order was made in the time t − τ in the past (13) in each entity. 4) Determine the supply capacity c(t) (11), backlog b(t) (12), and product output y(t) (6) in each entity. 5) Updating the inventory of each entity (5).

3.3.
Characteristics of the model. . In this paper, modeling of the SCN is almost real. The production capacity of factories and supply capacity of other entities are limited. The interaction between the SCN and customers is modeled with a nonlinear function that included the elasticity of demand as an important exogenous parameter. A decentralized approach is used for decision-making. In this approach, entities make local, independent decisions regarding demand forecast, inventory control, and order quantity. A new ordering policy based on the PD controller is used in the model that can improve decision-making, especially in the critical condition such as chaos.

4.
Simulation. Consider a SCN with ten entities: one factory, two distributors, three wholesalers and four retailers. Each entity has two ordering policies and two demand forecast methods; so there are four scenarios for decision-making (Table 1). It is assumed that all entities simultaneously use one scenario and that their parameters are the same. The total initial inventory, the desired inventory, and the demand of each level are distributed equally among its entities. Initial values and parameters are set according to Table 2. The model is simulated with the MATLAB software and in each scenario, 2000 data points are used to calculate the maximum LE. First, the presence of chaos in the SCN is proven. Suppose all entities make decisions using the Scenario 1. In other words, each entity applies smooth ordering policy and ES forecast method. With α = 0.21 and β = 0.65, the maximum LE is negative ( −0.0104) and thus the behavior of the SCN is stable. Fig. 4 shows the time series plot of the distributers total inventory (DTI) and the factories' total inventory (FTI) in this case. Fig. 5 shows the phase plot of DTI-FTI that attracted to a fixed point.
If α is slightly decreased, i.e., by 0.01, the maximum LE becomes positive (λ max = 0.0195 ) and the behavior of the system becomes chaotic (Fig. 6). The phase plot of DTI-FTI is shown in Fig. 7 and the attractor becomes complex. These figures suggest that the supply chain network is very sensitive to small changes of parameters and may exhibit a chaotic behavior. Now with different scenarios, effects of inventory adjustment parameter and demand elasticity are investigated on the behavior of the supply chain network. First, assume that an ES forecast method is used in all entities and the inventory adjustment parameter is constant at 0.21. Change the elasticity of demand from 0 to 2 and use two ordering policies. By calculating the maximum LE, the chaotic behavior of the SCN is studied. Fig. 8 shows that the ordering policy based on the PD controller (Scenario 3) is suitable and the behavior of the SCN is stable in a greater range of β. Now the elasticity of demand is kept constant (β = 0.65 ) and α is changed from 0 to 1. Fig. 9 shows the effect of the inventory adjustment In addition, the stability range is larger with this scenario.
Continuing simulation, the MA forecast method is used in all entities. First, the inventory adjustment parameter is kept constant (α = 0.21) and β is changed from 0 to 2, and the maximum LE is calculated and drown (Fig. 10). Scenario 4 (ordering policy based on the PD controller and MA forecast method) leads to chaotic behavior in a small range and is suitable for decision-making. Then, β is kept constant (β = 0.65) and α is changed from 0 to 1 , and the maximum LE is once again calculated and drown (Fig. 11). In this case, both ordering policies create chaotic behaviors in a large range, but there is less chaos with the ordering Comparison of forecast methods in face of the changes in the elasticity of demand is simulated with two appropriate scenarios, i.e. scenarios 3 and 4 (Fig. 12). In this manner, the ordering policy based on the PD controller and the MA forecast method (Scenario 4) is suitable. In other words, a decision is made better in face of marketing activities and demand fluctuations. Comparison of forecast methods in face of changes in the inventory adjustment parameter is simulated with scenarios 3 and 4 (Fig. 13). In this case, again the Scenario 4 is suitable.
For greater certainty, maximum LE is recalculated by changing α and β with an increment of 0.05 from 0 to 1 and 0 to 2 in all the scenarios. The number of Maximum LEs in different ranges is shown in Table 3. In the stable state (λ max ≺ 0),   all scenarios are compared in Fig. 14. The chaotic behaviors of the SCN with less intense (0 ≤ λ max ≺ 0.01) are compared in Fig. 15. Simulation results indicate that once again Scenario 4 is the most suitable one.  Figure 11. Effect of the inventory adjustment parameter by the MA forecast method Finally, the results of simulations show that the Scenario 4 is the most suitable scenario among four suggested scenarios. In other words, the ordering policy based on the PD controller and the MA forecast method is effective in reducing the chaotic behavior of the SCN. With this scenario, better decisions can be made in face of marketing activities and demand fluctuations.

5.
Conclusion. Modeling the supply chain as a network is closer to reality and the results are more valid for improving or controlling chaotic behavior in a real supply chain. Supply chain network is a complex system with nonlinear dynamics and may have a chaotic behavior. Each entity has the power to make decisions independently. The supply capacity of an entity is a significant variable in modeling and analyzing the chaotic behavior of the SCN. Time delay between order and delivery, interaction  Figure 13. Effect of the inventory adjustment parameter by the ordering policy based on the PD controller between entities and between the network and customers, the uncertainty of demand and production, and the impact of market activities make its behavior complex.
The ordering policy and the demand forecast method have the most important role in the behavior of the SCN. The ordering policy based on the PD controller plays a crucial role in stabilizing its behavior. This policy speeds up the decisionmaking process, especially demand fluctuations are larger. The MA forecast method is appropriate with this policy, with less intense chaotic behavior in the SCN.
The inventory adjustment parameter is an important internal decision variable and has a major role in controlling the inventory. The ordering policy based on the PD controller and the MA forecast method (Scenario 4) makes a better behavior of the SCN in face of changes in the inventory adjustment parameter. The elasticity of demand as an external parameter can indicate market activities and customers behavior. This parameter is not adjustable and the SCN should be robust against its changes. Again, the ordering policy based on PD controller and the MA forecast method (Scenario 4) is more appropriate.
Controlling chaotic behavior in the SCN by control methods such as robust control, adaptive control, and sliding mode control would be an interested area for future investigation.