Stability of the Supply Chain Using System Dynamics Simulation and the Accumulated Deviations from Equilibrium

. We propose and demonstrate a new methodology to stabilize systems with complex dynamics like the supply chain. This method is based on the accumulated deviations from equilibrium (ADE). It is most beneﬁcial for controlling system dynamic models characterized by multiple types of delays, many interacting variables, and feedback processes. We employ the classical version of particle swarmoptimizationas the optimizationapproach due to its performance in multidimensionalspace, stochastic properties, and global reach. We demonstrate the e ﬀ ectiveness of our method based on ADE using a manufacturing-supply-chain case study.


Introduction
The current world economic climate has led to highly volatile markets and fluctuating demands for manufactured goods.Consequently, demand forecasting and inventory control have become major concerns for supply chain managers.In recent papers [1][2][3], the authors demonstrated that the classical supply-chain approach to forecasting future inventory demand [4] is not enough to avoid huge swings in finished goods and in-process inventories.The authors implied that a new approach to modeling supply chain behavior based on a deep understanding of the structure of the problem, the causalities among the system variables, and the different types of time delays is needed.This echoes results from [5] who said "behavioral functions that underlie market actions can be more successfully interpreted through nonlinear decision functions" and "that theory should be drawn from actual causality relationships".Modeling such nonlinear decisions and causal relationships is the focus of systems dynamics (SD).We believe that system dynamics combined with elements from nonlinear control theory, calculus, and optimization can be used to understand, monitor, and stabilize supply-chain behavior.
System dynamics models of supply chains provide a mathematical interpretation of the intricate relationships between numerous variables.Sets of those variables can form feedback loops based on the causality structure of the problem.These intermingled feedback loops can act simultaneously, but at different times, they may have diverse effects.Therefore, transitions from a regime to another are widespread.Most system dynamics models used for supply chains are complex and nonlinear [6].Nonlinearities due to the multiplying effect of multiple causal variables upon another variable are "quite common in SD and difficult to tackle" [7].The number of variables and feedback loops, the discontinuities, and time delays increases the level of complexity [8].Time delays being added to "negative feedback loops increases the tendency for the system to oscillate" as stated by Sterman [9].Therefore, this level of complexity makes more difficult the task to design stable supply chains.
Instability is a major cause for poor performance in supply chains [4,10].Traditional methods in structural analysis can be used to investigate and stabilize supply chains using SD modeling.One of these methods is loop knockout [10].Loop knockout is a brute force method that relies on testing various feedback loops.This analysis of the different feedback loops is used to visualize the long-term effects of changes that could be implemented in order to stabilize the supply chain.This method can be effective for very simple models.However, it is not enough for medium to higher levels of complexity.Another method is eigenvalue analysis.Eingenvalue analysis is used to characterize the modes of Modelling and Simulation in Engineering behavior [11][12][13].The analysis is based on the linearization of the mathematical description at infinitesimal points and then study the evolution of the eigenvalues and their elasticities [14].The utilization of this method cannot be easily generalized to nonlinear models.Furthermore, methods in structural analysis rely on sensitivity analysis to determine the parameter values of the stabilization policies.Therefore, we propose a new alternative to stabilize the supply chain modeled as a dynamic system based on the ADE.
According to [15], the concept of stability is linked to the notion of equilibrium points-"an equilibrium point (EP) is stable if all solutions starting at nearby points stay nearby; otherwise it is unstable".He also described the notion of asymptotic stability which occurs when all solutions tend to the EP as time goes to infinity.In nonlinear control theory, stability of equilibrium points can be determined by (1) linearization around the EP and the values of the associated eigenvalues or (2) Lyapunov functions [15] when linearization is inadequate to approximate the global behavior of the system [16].
Note, however, that it is not always easy to find or construct a Lyapunov candidate function for a specific system.Because of this, we propose the use accumulated deviations from equilibrium.It is very well known from optimization-based control theory that minimizing the deviations of controlled variables from some desired level can have the effect of generating stable solutions (see the Appendix).We have developed a methodology based on this idea-stability analysis based on the accumulated deviations from the equilibrium (SADE)-in order to optimize and stabilize nonlinear systems.The optimization engine for this methodology is particle swarm optimization (PSO).
We will use the SADE methodology together with the concept of asymptotic stability to minimize oscillatory behaviors of specific control (state) variables-such as inprocess or finished goods inventory.If necessary, stability can be extended to the whole system by using a weighted average function that includes all state variables.This also allows higher weights to be assigned to those variables considered more important.Since our methodology does not require linearization of the system or eigenvalue calculations, it can be applied as a general procedure to linear or nonlinear dynamic systems.This method is novel from the perspective of SD modeling.
The paper is organized as follows.In Section 2, we describe the SADE methodology including the optimization problem and PSO implementation.In Section 3, we give details about the supply-chain case study and relevant systems dynamics literature.Those details include the definition of the optimization problem, procedures for testing policy robustness, and the respective comparison with a local search algorithm.We conclude this paper with further discussion of some key points and proposed future work.

Sade Methodology
2.1.Description.The SADE methodology and its general functioning are shown in Figure 1.The supply chain environment represents the actual participants, structure, strategies, policies, objectives, variables, constraints, and parameters of the real system.This environment captures all of different scenarios of the supply chain (SC) over time.These scenarios, together with the associated decisions, produce changes in the behavior of the supply chain.By behavior, we mean the observed patterns in the state variables over time.
The SD model replicates the dynamic behavior of the supply chain environment.We chose the SD modeling approach based on the advice in [5]-it can capture the causal relationships, feedback processes, and multiple time delays necessary to track accurately behavioral evolution of the system.As exogenous events occur in the SC environment, their impacts on the behavior of the system are predicted using the SD simulation.If those predicted impacts do not show any instability patterns, no actions are taken.Otherwise, a new management policy must be found to remove the instability or minimize its impact.
This new policy is found using a PSO algorithm.That algorithm modifies the set of parameters that constitute the current policy until the ADE is minimized.In every iteration of the algorithm, the parameter set is sent to the SD model in order to calculate, through simulation, the value of the ADE (objective function).Simulation is used due to the difficulty of solving the complex dynamic equations by analytical methods.The optimization problem and the PSO algorithm are described in Sections 2.2 and 2.3.
Once the best setting of parameters (stabilization policy) is obtained, then it is implemented in the actual supply chain.

Optimization Problem.
The SD model can be described by an equation of the form ẋ(t) = f(x(t), p), where x(t) is the vector of state variables (dimension n) and p is a vector of adjustable parameters (dimension q) with lower and upper bounds p L and p U , respectively.
We can formulate an optimization problem that will find the parameter vector p * that causes the state variable x s to become asymptotically stable (see the Appendix) around the equilibrium point x eq s (p * ).We will find this optimal parameter vector by minimizing the ADE (see the Appendix for the mathematical definition of ADE) for predetermined time horizon T and making use of Theorem 1 (see the Appendix).That is, we will find the vector that makes ADE converge.The optimization problem is then stated as where w s ≥ 0, ( No actions needed

Supply chain environment
Capture The objective function J(p) is defined as the weighted average value of the ADE, and T is the time horizon.The use of weights, w s , means that J(p) will support the simultaneous stabilization of any subset of m state variables (m ≤ n).This allows higher weights to be assigned to the variables that are considered more important.
If we do not know the equilibrium point x eq s in advance, we can modify J(p) to include it as a variable (a s ) in the parameter vector.This step is supported by the results of Theorem 2 (see the Appendix).
The objective function defined in (1) can be incorporated very easily into any SD formulation by adding a "stock and flow" piece to the model that is linked to the state variables of interest as illustrated in Figure 2.Then, we define the variables DE and ADE as (2) 2.3.PSO Algorithm.Particle swarm optimization (PSO) was invented in the mid-1990s by Kennedy and Eberhart [17].PSO is conceptually simple and can be implemented in a few lines of code.In comparison with other stochastic optimization techniques like genetic algorithms (GAs) or simulated annealing, PSO has fewer complicated operations and fewer defining parameters [18].PSO has also been shown to be more computationally efficient than GAs when applied to unconstrained nonlinear problems with continuous variables Hassan et al. [19].
The specific algorithm we use is called "local best PSO" [20].This algorithm is based on a social network composed of neighborhoods related to each particle.The algorithm maintains a swarm of particles, where each particle represents a candidate solution to the optimization problem.These particles move across the search space communicating good positions to each other within the neighborhood and adjusting their own position and velocity based on these good positions.For this purpose, each particle keeps a memory of its own best position found so far and the neighborhood best position among all the neighbor particles.The goodness of a position is determined by using a fitness function.A typical stopping condition of the algorithm is when the maximum number of iterations has been exceeded.The basic elements of the algorithm are defined as follows.
(i) Particle.A particle i is represented by a n p -dimensional real-valued vector p i .This vector is composed of particle positions p i j ; that is, p i = [p i1 , p i2 , . . ., p inp ].Each particle position corresponds to one of the parameters of the parameter vector defined in the optimization Section 2.2.
(ii) Swarm Size.It is the number of particles in the swarm, and it is denoted by N .
(iii) Fitness Function.It is a mathematical function used to quantify how good the solution represented by a particle is.For a particle i, the fitness function is the objective function J(p i ) as defined in Section 2.2.
(iv) Personal Best Position.As a particle moves through the search space, it compares its fitness value at the current position to the fitness value it has ever attained so far, which is called the personal best position.For each particle i, the personal best position can be expressed as the real-valued vector y i = [y i1 , y i2 , . . ., y inp ], and it is determined so that (v) Neighborhood Size.Defines the extent of the social iteration within the swarm [20].Selection of neighborhoods was done based on particle indexes.Each particle has a neighborhood associated to, where B i defines the set of indexes for the neighbors of particle i.
(vi) Neighborhood Best Position.It is the best position among all the personal best positions in the neighborhood.It is denoted by the real-valued vector y i = [ y i1 , y i2 , . . ., y inp ], and it is determined so that J( y i ) ≤ J(y j ), j ∈ B i .
(vii) Global Best Position.It is the best position among all the personal best positions achieved so far in the entire swarm.It is denoted by the real-valued vector g = [g 1 , g 2 , . . ., g np ], and it is determined so that J(g) ≤ J(y i ), i = 1, . . ., N .
(viii) Particle Velocity.It is the velocity of the moving particle i represented by the real-valued vector This vector reflects both the experiential knowledge of the particle and socially exchanged information from the particle's neighborhood [20].The experiential knowledge of a particle is generally referred as the cognitive component, which quantifies the performance of particle i relative to past performances.It is represented by the term c 1 r 1 (y i − p i ).The socially exchanged information is referred as the social component of the velocity equation.It is represented by the term c 2 r 2 ( y i − p i ).
(ix) Acceleration Coefficients.The acceleration coefficients, c 1 and c 2 , together with the random vectors r 1 and r 2 , control the stochastic influence of the cognitive and social components on the overall velocity of a particle [20].The constants c 1 and c 2 are also referred to as trust parameters, where c 1 expresses how much confidence a particle has in itself, while c 2 expresses how much a particle has in its neighbors.The random vectors are defined as r 1 = [r 11 , r 12 , . . ., r 1np ] and r 2 = [r 21 , r 22 , . . ., r 2np ], where r 1 j and r 2 j are uniformly distributed random numbers in [0, 1].
(x) Inertia Weight.It is a parameter "w" that is used to control the influence in the new velocity of a particle by its previous velocity (flight direction).Thus, it influences the tradeoff between the global and local exploration abilities of the particles [21].For initial stages of the search process, where global exploration is required, it is recommended to set a large inertia weight, while for the last stages, the inertia weight should be reduced for better local exploration.A decrement function for decreasing the inertia weight at the iteration k can be given by w(k) = αw(k ), where α = 0.98 and k is the last iteration when the inertia changed.A parameter "iteration lag" is defined to set the number of iterations since the last change in the fitness function that are required to change the inertia weight.
The steps of the algorithm are described in the following lines.
(v) Set the initial value of the personal best position vector as y i (0) = p i (0), i = 1, . . ., N .
(viii) Set the initial value of the inertia weight w(0).
Step 3. Weight updating: if the fitness function has not decreased in a number of iterations equal to the "iteration lag," then update the inertia weight using w(k) = αw(k ).
Step 4. Velocity updating: calculate the velocity of particle i by using Step 5. Position updating: based on the updated velocities, each particle changes its position according to the following equation: Step 6.Personal best updating: determine the personal best position visited so far by each particle.
(ii) Set Step 7. Neighborhood best updating: determine the neighborhood best position y i (k) visited so far by the whole swarm by using the formula Step 8. Global best updating: determine the global best position g(k) visited so far by the whole swarm by using the formula Step 9. Stopping criteria: if the maximum number of iterations is achieved, then stop, g * = g(k) is the optimal solution; otherwise, go to Step 2.

Differences between Our Approach and Current PSO
Applications in Supply Chain Management.The contribution of our approach is the development of the ADE as objective function in order to solve the problem of instabilities in the supply chain.PSO has been selected as the mechanism to find a solution.The current PSO applications in supply chain management are very different from our approach.The current PSO applications emphasized straight forward solutions to static supply chain problems without regards to stability, robustness, and system dynamics [22][23][24][25][26].

Comparisons with Other Optimization
Algorithms.Initial comparisons were performed against Genetic Algorithms, PSO, PHC, and methods based on eigenvalue analysis.However, the best scheme was based on the combination of PSO and PHC.Our major contribution is the utilization of ADE.Therefore, PSO and PHC solve the problem very efficiently (timeliness and performance-based) [27].Our further research work will make emphasis on the utilization of other algorithms.

Supply Chain Case Study
Our case study focuses on inventory and labor concerns in a manufacturing supply chain.The nonlinear SD model of this supply chain is a simplified version of Sterman's original model [9].It has two submodels: inventory (Figure 3) and labor (Figure 4).The inventory management sector (Figure 3) is represented by two state variables: inventory and work in process inventory.The variable work in process inventory represents all the stages of the production process, where intermediate inventory is created.The variable Inventory represents the finished goods inventory.This model assumes that orders are filled as they arrive and the ones that cannot be filled immediately are lost as customers seek other sources of supply.The labor sector (Figure 4) is represented by two state variables: vacancies and labor.The stock of vacancies is the supply line or order of workers that have been placed but not yet filled.This states that workers cannot be instantly hired.Hiring takes time: positions must be authorized, and vacancies must be created.The labor force is a stock of people, which is increased by the hiring rate and decreased by the quit rate.This last rate includes voluntary quits and retirements, excluding the possibility of layoffs.
The supply-chain behavior is impacted by interactions between inventory-management policies and laboradjustment policies.To capture the impact of these policies, the model uses four state variables and several parameters.Two of those parameters-productivity and standard workweek-have constant values of 0.25 widgets/person and 40 hours/week, respectively.The remaining nine parameters (see Table 1) are considered variables that managers can set to design stabilization policies.The current values for these parameters are shown in Table 1.The goal is to find a policy that maintains the inventory and labor state variables at equilibrium and avoids large oscillations in the inventory.
Customer orders are arriving at the rate of 10000 widgets/week.After the system remains in equilibrium for the first five weeks, customer orders experience a linear increment for the next fifty weeks until reaching 20% of its original value, where they remain constant.The resulting behavior of the variables inventory and labor (Figure 5) shows several oscillatory fluctuations.These fluctuations are caused by delays in production.

Optimization Problem.
In order to determine a stabilization policy for these two state variables, we solved the following optimization problem using our proposed global search algorithm (PSO).The algorithm was run at the fifth week using the setting mentioned below: (1) number of iterations = 150,  These parameter values were chosen after doing some initial experiments with the empirical rules selected to guide the choice (see Table 2).
The new parameter set associated to the stabilization policy is shown in Table 3.The objective function (ADE) was improved by 82%.It took 189 seconds to calculate this policy after 150 iterations of the algorithm (the algorithm was executed on a 1.86 GHz Pentium PC with 1 GB of memory).
The minimization problem considered the first two state variables as the variables of interest.The following weights were assigned w 1 = 0.6, w 2 = 0.4 to represent the concern of management in the inventory Figure 6 shows the behavior of the state variables when the stabilization policy is applied at the fifth week.It clearly shows the convergence of the ADE with minimal oscillation.Two main actions have been taken to respond to the change in customer orders.First, the increment in production has been achieved by decreasing the manufacturing cycle time.Second, by decreasing the time to adjust, labor will help to  approach production more closely to the desired production rates.As a result, the EP of the inventory variable has increased from 40000 units to 93359 units.The EP for the labor variable remains not far from its original value of 1000 people.Stabilization of both variables has been achieved in approximately 80 weeks.

Testing for Policy Robustness.
The stabilization policy was tested by generating a sudden change in weeks 80 and 100 in the customer orders and showing the system's response to this change.The different percentage changes in customer orders and responses are shown in Table 4. Figure 7 depicts the robust behavior of the inventory and labor variables to the changes.These variables showed a sharp increase or decrease in their levels necessary to adapt to the changes before reaching new equilibrium points.Stability returns approximately in week 130. optimal solution is an important factor in selecting a search algorithm, the quality of such solution in terms of oscillation reduction has to be analyzed.For that reason, we decided to compare the results obtained by solving the optimization problem using our global search algorithm PSO (Figure 6) and the ones obtained by using the local search algorithm Powell hill climbing (PHC) [29].We tested two starting EPs for the PHC algorithm.In the first test (PHC1) it was used the original EP of the model (for the inventory and labor variables) as the starting point which was (a 1 , a 2 ) = (40000, 1000).For the second test (PHC2), we used the EP obtained by the PSO algorithm, which was (a 1 , a 2 ) = (93359, 1200).In both tests it took around 15 seconds to find the solution.We can see from Figure 8 that after some fluctuations policy PHC1 achieves stabilization of the inventory and labor variables approximately in 70 weeks.Policy PCH2 stabilizes the system in 120 weeks at a lower EP that the one it started.

Comparing
Comparing the policies obtained by PCH and PSO, we can see that the second one (PSO) shows less fluctuations before reaching the EP due to the higher inventory level.Both algorithms generate the same EP for the labor variable.However, the policy with PSO gets that by increasing the manpower smoothly after a small trough, while the policies with PCH obtains the stability after several decreasing fluctuations.

Conclusions
This paper proposes a methodology, based on the ADE, to eliminate or minimize oscillatory behaviors of the supply chain.Our approach utilizes the modeling flexibilities of system dynamics to model the supply chain and the potential of the PSO algorithm to scan the search space to solve the stabilization problem.
We tested our methodology by increasing and decreasing the customer orders in a manufacturing supply chain.We showed that our approach can stabilize the behavior of the state variables despite these fluctuations.We concluded that after a sudden perturbation of the system the stabilization  policy remains stable requiring a period of adaptation to the changes before reaching new EPs.We compared the results of our proposed approach that uses the global search algorithm PSO with the ones obtained by the local search algorithm PHC.Because PSO does not depend on the starting point of the state variables, it provides a more expanded and deeper search of the space to find the EP.We concluded that the PSO algorithm provided a better solution than the PHC algorithm in terms of fewer oscillations.PHC achieved faster stability when started from the original EP of the model.

Future Work
We propose to test the performance of this local best PSO algorithm using a real-size model of the supply chain.It will be interesting to compare the results of this algorithm with the ones obtained with other evolutionary algorithms such as GAs and PSO variations and hybrids.
Currently, this methodology is focused in generating stabilization policies at the strategic and tactical levels of the supply chain where SD modeling is more suitable.We propose to extend this research to the operational level by developing a methodology that will propagate stability from the higher to the lower level of the supply chain.We plan to meet this challenge by using hybrid simulation (SD and discrete event simulation) to model the different levels of the supply chain and a hierarchical approach to coordinate between the policies obtained at these levels.

Appendix Related Definitions and Theorems
Here, we present some important definitions and theorems that provide the support to understand our methodology.Definition 1.The point x eq ∈ R n is said to be an equilibrium point of the differential equation ẋ(t) = f(x(t))( ẋ(t) = ∂x(t)/∂t) if it has the property that once the corresponding system reaches x eq at time t eq , it will remain at x eq for all future time; in other words, f(x(t)) = 0 for all t ≥ t eq .Definition 2. Consider the system defined by ẋ(t) = f(x(t)); x(0) = x 0 , where x(t) ∈ R n ; f : R n → R n ; x(t) = [x s (t)] = [x 1 (t), x 2 (t), . . ., x n (t)] T , s = 1, . . ., n.The state variable x s is defined to be stable (around the EP x eq s ) if it is bounded; that is, there is a finite number M s such that |x s (t) − x eq s | ≤ M s (the symbol |c| represents the absolute value of c).If this condition holds for all state variables, then the system is said to be stable.Definition 3. Consider the system defined by ẋ(t) = f(x(t)); x(0) = x 0 , where x(t) ∈ R n ; f : R n → R n ; x(t) = [x s (t)] = [x 1 (t), x 2 (t), . . ., x n (t)] T , s = 1, . . ., n.The state variable x s is defined to be asymptotically stable (around the EP x eq s ) if it is both stable (satisfies Definition 2), and additionally, we have Lim t → ∞ (x s (t) − x eq s ) → 0. If these two conditions hold for all state variables, then the system is said to be asymptotically stable.Definition 4. Consider the system defined by ẋ(t) = f(x(t)); x(0) = x 0 , where x(t) ∈ R n ; f : R n → R n ; x(t) = [x s (t)] = [x 1 (t), x 2 (t), . . ., x n (t)] T , s = 1, . . ., n.For the state variable x s , the accumulated deviations from its EP x eq s is defined as Theorem 1.Consider the system defined by ẋ(t) = f(x(t)); x(0) = x 0 , where x(t) ∈ R n ; f : R n → R n ; x(t) = [x s (t)] = [x 1 (t), x 2 (t), . . ., x n (t)] T , s = 1, . . ., n.The state variable x s is asymptotically stable (around the EP x eq s ) if

Figure 1 :Figure 2 :
Figure 1: General procedure of the SADE methodology.

Figure 5 :
Figure 5: Behavior of state variables for the current policy.

Figure 7 :
Figure 7: Behavior of state variables after a sudden change in customer orders.

∞0
|x s (t) − a s |dt converges, then a s = x eq s .

Table 3 :
Parameter values for policy using PSO algorithm.
Polices with a Local Search Algorithm.This methodology does not require to find the global optimum to obtain satisfactory reduction in instability.Using a local search algorithm, we can obtain a quick convergence of the ADE in just few seconds.Although the time to find the

Table 4 :
New EPs reached after the system was perturbed.