87 | Modelling Nonlinear Supply Chains Using Methods of System Dynamics Modeliranje nelinearnih oskrbovalnih verig z metodami dinamike sistemov Rudolf Pušenjak 1, *, Maks Oblak 1 1 Faculty of Industrial Engineering Novo mesto / Šegova ulica 112, SI-8000 Novo mesto E-Mails: rudolf.pusenjak@fini-unm.si; maks.oblak@fini-unm.si * The corresponding author: Full Prof. Dr. Rudolf Pušenjak, Kidričeva 5 Kidričeva 5, 9240 Ljutomer, Slovenia; Tel.: +386-2-582-14-06 Abstract: This paper treats a three level nonlinear supply chain, which consists from a manufacturer, a distributor and a retailer as a dynamical system. The dynamic behavior of supply chain is modelled by using the continuous Lorenz-like model with disturbances. The Lorenz model is known from physics in the study of dissipative hydrodynamical systems with excitation and represents a paradigmatic example of the deterministic system, which is able to exhibit chaotic motions. Very complex supply chain relationships are modelled due to the information flow distortion at the distributor, caused by disturbances, which appear at the retailer. By using the proposed nonlinear supply chain model, the steady state, the saturation and the chaotic state of the supply chain are analyzed, respectively. The range of parameter values, which are characteristic for the appearance of such states are determined. In the paper, the meaning of the regular state in which the supply chain cannot keep due to the disturbances and the risk of the saturation state are explained. In the analysis of such states, the multistage homotopic perturbation method is used. The gained results are checked analyticaly or by the help of standard numerical integration, such as the Runge-Kutta method . Key words: Nonlinear supply chains, dynamics of Lorenz-like model, characteristic states, multistage homotopic perturbation method. Povzetek: Članek obravnava trinivojsko nelinearno oskrbovalno verigo, ki sestoji iz proizvajalca, distributerja in trgovca na drobno kot dinamični sistem. Dinamično obnašanje oskrbovalne verige je modelirano z zveznim modelom z motnjami, ki je podoben Lorenzovemu sistemu. Lorenzov system je znan iz fizike v študiju disipativnih hidrodinamičnih sistemov z vzbujanjem in velja za enega izmed tipičnih determinističnih sistemov, ki lahko izvajajo kaotična nihanja. V oskrbovalni verigi modeliramo kompleksne razmere zaradi popačenja informacijskega toka pri distributerju, ki jih izzovejo motnje, s katerimi se sooča trgovec. S pomočjo modela je izvedena analiza ustaljenega stanja, stanja nasičenja in kaotičnega stanja oskrbovalne verige. Določena so območja vrednosti parametrov oskrbovalne verige za nastanek posameznih stanj. V članku je pojasnjen pomen regularnega stanja oskrbovalne verige, v katerem oskrbovalna veriga zaradi motenj v splošnem ne more vztrajati in stanje nasičenja kot enega najnevarnejših stanj oskrbovalne verige. V analizi teh stanj je v članku uporabljena večstopenjska homotopsko perturbacijska metoda. Dobljeni rezultati so preverjeni analitično oziroma s pomočjo standardne numerične metode Runge-Kutta. Ključne besede: Nelinearne oskrbovalne verige, dinamika sistema, karakteristična stanja, večstopenjska homotopsko perturbacijska metoda. Izvirni znanstveni članek TEHNIKA ANALI PAZU 3/ 2013/ 2: 87-93 www.anali-pazu.si MODELLING NONLINEAR SUPPLY CHAINS USING METHODS OF SYSTEM DYNAMICS | 88 1. Introduction Supply chains and networks consist of many players such as suppliers, manufacturers, distributors and retailers. Each individual player in supply chain is trying to achieve the maximum profit and to ensure the maximum customer satisfaction. Besides this individual goals, all players also have global objectives in order to create a resilient and competitive supply chain. According to Christopher [1], competition in the future belongs to competing supply chains, not to competing individual organizations. Supply chains are dynamic systems, which are characterized by the downstream material flow in the direct path and an information flow, running in both downstream direct path and upstream feedback path. Both flows are subjected to the various kind of temporaly varying uncertainties, such as demand uncertainty, supply uncertainty, delivery uncertainty and forecasting uncertainty, because supply chains are working in an increasingly dynamic and volatile global market-place. Supply chain structure, which consists from the direct path and feedback path is typical for the control systems, known from technical disciplines. Application of the nonlinear theory of dynamical systems together with a great computational power of today's computers makes modelling and simulation of supply chains possible and offers a deep understanding of their behavior. 2. Nomenclature k = the safety stock coefficient at the manufacturer r = the information distortion rate at the distributor m = the customer order satisfaction rate at the retailer i = the discrete time period x i = the quantity of products demanded by retailer in i th time period y i = the quantity of products, which can be supplied by distributor in i th time period z i = the quantity of products, which can be produced by manufacturer in i th time period in dependence on the order 3. Three level model of nonlinear supply chain Modelling of supply chains is faced in general with setting, optimization and control of the system model. Application of principles of dynamics in business systems was introduced by Forrester [2] in the sixties of 20th century. First attempts were based on the linear theory of dynamical systems with the most important principle of superposition. The application of the linear theory in supply chain modelling has led to important findings regarding the supply chain stability. Furthermore, linear theory was able to explain the bullwhip phenomenon. However, due to the superposition principle, it was unable to predict a wide variety of states, which can supply chain passes. The main reason why linear theory of supply chains is not sufficient is the fact, that linear systems cannot exhibit chaos. Although chaotic systems are not random but deterministic, they are sensitive on small changes of initially conditions, which cause unpredictable system responses. This phenomenon of unpredictability is the characteristic property of chaotic systems, which cannot happen in linear supply chains. In order to understand supply chain phenomena in a very complex system, the nonlinear supply chain model must be envisaged (see [3] and references therein). Supply chains have in general a multilevel structure or even are forming networks, however they are organized to connect three key players: manufacturers, distributors and retailers. In order to apply analytical tools in a very complex nonlinear supply chain, a three level supply chain, which consists of one manufacturer, one distributor and one retailer, is modelled in this paper. Two kinds of flow characterize the supply chain: the product flow, which flows downstream from the manufacturer via distributor to the retailer and information flow, which flows downstream and upstream, respectively through the supply chain. Both downstream flow running in a so called direct path, while the upstream information flow is running in the so called feedback. The presence of feedback makes the supply chain controllable. Beside of direct path and feedback, respectively, there exist several interconnections in each stage, which make the supply chain nonlinear. The structure of three level nonlinear supply chain treated in the paper is shown in Fig.1. The mathematical model of the supply chain consists from three equations, one equation for each player. Thus, if desired, the mathematical model can be expanded with additional equations, corresponding to the multilevel supply chain structure, but on the cost of the inevitable increase in the mathematical complexity. To simplify the structure of nonlinear supply chain on three level, the dynamics between the supplier and manufacturer is incorporated into the manufacturer. Due to the simplified structure, uncertainties, which can be generated by supplier, are not considered in the model. Instead of this, it is assumed, that the manufacturer holds the safety stock in order to reduce the small production batches and eventually to satisfy the increased demand. The height of the safety stock is regulated by means of the safety stock coefficient k at the manufacturer. Despite the above simplification, the model can exhibit the great complexity at different levels of the supply chain and contains the inventory concept at all stages. The emphasis in this model is on the information distortion rate r at the distributor, which plays the role of a factor that is influenced by uncertainties from the retailer. Before the demand information is presented to the manufacturer, it is distorted in order to avoid the uncertainties generated by retailer. The retailer self tries to avoid the uncertainties that pass from the customer and distributor, using his own intuition and experience. The customer order satisfaction rate m is introduced at the retailer as a means to avoid the uncertainties that are generated by the customer. It is also assumed, that the customer order satisfaction rate can be used for mitigating the uncertainties, which come from the distributor. ANALI PAZU, 3/ 2013/ 2, str. 87-93 Rudolf PUŠENJAK, Maks OBLAK 89 | In order to form a mathematical model of the supply chain, depicted on Fig.1, it is assumed that the demand information, which flows upstream the chain from the retailer to the distributor and from the distributor to the manufacturer, respectively, needs a unit time to be transmitted along links of the supply chain. In other words, the distributor processing of demand information within the i th time period is what retailer requested in the i-1 th time period. Similarly, the manufacturer processing of demand information within the i th time period is what distributor requested in the i-1 th time period. The delay in processing of demand information causes the information distortion. The important goal of managing the supply chain is not only to mitigate the influence of uncertainties on the ordering policy, but also the fulfillment of the market demand. The quantity of demanded products in the current i th time period is always dependent on the demand satisfaction in the previous i-1 th time period. This is the key principle in the managing of the supply chain because it determines the ordering pattern upstream the chain stages. Decision of the individual player in the supply chain on the order and his inventory in the current i th time period is made after considering the information distortion and his present inventory level. As is mentioned previously, the demand information is transmitted from the stage to the subsequent upstream stage of the supply chain with a delay of one unit time. As can be seen in Fig. 1, the outcoming order quantity is not the same as the requested incoming order quantity at any level. The order quantity x i of the retailer in the i th time period is linearly coupled with the distributor and depends on the demand my i-1, which is satisfied by distributor in the previous period of time and on the corresponding order quantity of the retailer mx i-1 in the previous time period, where m is the ratio at which the customer demand is satisfied. Thus, the equation of the retailer, which reflects the quantity demanded by customers, reads:   11 i i i x m y x   (1) The downstream as well as the upstream information flow between distributor, the producer and the retailer is nonlinear due to the strong coupling. The quantity of products y i, which can be supplied by distributor in i th time period is affected by the combined effect x i-1z i-1 of the retailers order quantity x i-1 and the manufacturers quantity of produced goods z i-1 in the previous time period, which give rise to the quadratic nonlinearity. The distributor takes into account that the order information received from the retailer might be distorted in the amount rx i-1 and reduces the inventory level for the amount x i-1z i-1. Thus the distributor equation of the quantity of products y i, which can be supplied in i th time period, is given by the relation:   11 i i i y x r z   , (2) where r is the information distortion coefficient. The third equation is the manufacturer equation, which gives the production quantity in the i th time period. The production quantity under given assumptions depends on the distributor's order and the safety stock. Of course, distributor's order depends again on the reailer's order x i-1 in the previous time period, so that the manufacturer must take into account the combined effect x i-1y i-1 of the retailer and distributor output. Before the manufacturer makes the production decision z i for the i th time period, the quadratic nonlinear contribution x i-1y i-1 is increased for the share of safety stock kz i-1: 1 1 1 i i i i z x y kz     , (3) where k is the safety stock coefficient at the manufacturer. System of Eqs. (1), (2) and (3) describes a discrete model of the nonlinear supply chain. However, the discrete time unit of one period is quit small in respect to the long term behavior. Thus, instead of the discrete system of equations, it is resonable, to study the continuous system, which can be derived from the discrete equations by approximating the difference quotients through time derivatives, defined in the following way: 1 1 1 1 1 1 1 1 1 1 1 1 ii T i i i ii T i i i ii T i i i xx x x x T yy y y y T zz y y z T                                  (4) By using substitutions (4), the discrete equations (1), (2) and (3) are rewritten on the continuous form:       1 1 x my m x y x r z y z xy k z               , (5) Fig. 1. Structure of three level nonlinear supply chain. MODELLING NONLINEAR SUPPLY CHAINS USING METHODS OF SYSTEM DYNAMICS | 90 where index i-1 is omitted on both sides of equations (5), because it is meaningless in the continuous system. The system of equations (5) is the Lorenz-like system, which can exhibit chaotic motions in the proper range of system parameters. In fact, the system (5) differs from original Lorenz equations:   x y x y rx xz y z xy bz       , (6) only in the first equation. Second and third equation in both systems are in fact the same, because the parameter b in third equation of the system (6) can be interpreted as the substitution for the term 1- k in the third equation of the system (5). However, the first equation of the system (5) can be rewriten on the form   1 x my m x      d yx   , where 1 1 1 1, 1 mm d m m               (see Ref. [3]). Thus, first equation in both systems differ in fraction 1/σ, which appears in the quotient m/σ = 1 - 1/σ. Parameters σ, r, b of the original Lorenz system can take positive values only. Due to the small differences between Lorenz system and the system (5), the values of parameters m, r and of the term 1-k take positive values, too. In the study, the values of parameters m, r, k will be fixed giving the so called "reference supply chain model". Until now, the time dependent perturbations of uncertainties are not considered at all. The external disturbances have the additive character and can affect any of the three levels of the supply chain. When uncertainties are added at one stage, they can propagate both upstream and downstream, which is the common phenomenon exhibited in a real supply chain. Due to the additivity, the uncertainties at each level of the supply chain can be simply analyzed by means of the generalized Lorenz-like system (5) in the following form:           1 2 2 1 1 x my m x d t y rx xz y d t z xy k z d t                         , (7) where primes denote the change of variables in respect to the variables of the reference supply chain model. In the following study, we restrict the form of disturbances d 1(t), d 2(t) and d 3(t) on the periodic (harmonic) case only and thus allow the periodic seasonal variations. Influence of internal disturbances on the behavior of the presented nonlinear supply chain model is not studied in this paper due to the limited space, although this study is quit possible, when the temporal dependencies of parameters m, r, k are allowed. 4. States of nonlinear supply chain in the reference model Equations (5) describe the evolution of the nonlinear supply chain as the time increases. Due to the nonlinear character of equations, the supply chain undergoes different states through the evolution. Generally speaking, the evolution of a (nonlinear) dynamic system is divided into two phases. The first phase belongs to the so called transient phenomenon, which appears after the system is initially disturbed and exponentially dies after a short period of time. After passing the transient phenomenon, a dynamic system enters the second phase, which may be a saturation state, a steady state, or chaotic state. In the sequel, these states are discussed only in respect to the Lorenz-like system (5). Saturation state corresponds to the equilibrium or fixed point of the system, where derivatives on the left hand side of Eqs. (5) vanishes. This means, that all dynamic forces of the system are mutually canceled and the system is at rest. While this state is a desired state in physics and engineering, it is a dangerous state in economics and supply chain managing. Existence and stability of fixed points depend on parameter r. For 1 1 m r , there exists only one fixed point at the origin x=0, y=0, z=0, which can be proved stable. For 1 1 m r , the origin loses its stability and a pair of new fixed points appear, which are stable in the range 1 1 c m rr    , where r c denotes the critical value of the parameter r. By solving the nonlinear algebraic system of equations:   1 0, my m x    0, rx xz y      10 xy k z    , where m>0, 1 1 m r , 1> k >0, the following nontrivial pair of fixed points:             1 1 1 1 1 1 , 1 1 , m m m m e e e m m m m x k r y k r z r                       (8) is obtained, where the subscript e denotes the equilibrium. Stability of solutions of equations (5) is examined through linearization of the system near the equilibrium point. By using Jacobian matrix, the linearized equations are written successfuly in the following matrix form:  10 1 1 ee ee x m m x y r z x y z y x k z                                   , (9) where δx, δy, δz are small perturbations around the equilibrium point, giving a nonequlibrium state ,, e e e x x x y y y z z z          and where ,, x y z    are corresponding perturbations of derivatives. From (9) we derive the characteristic equation of eigenvalues: (10) When the equilibrium point is the origin, the characteristic equation simplifies:       2 2 1 1 1 0 m m r k              , (11) giving three eigenvalues: ANALI PAZU, 3/ 2013/ 2, str. 87-93 Rudolf PUŠENJAK, Maks OBLAK 91 |       2 2 2 4 1 1 1 2,3 2 1, m m m r k              . (12) From (12) it is clear, that all 1 2 3 ,,    have negative values when 1 1 m r , (m>0,1> k >0), however 13 0, 0   and 2 0   when 1 1 m r . This proves, that the origin is stable when 1 1 m r and becomes unstable when 1 1 m r . The stability around the pair of equlibrium points     1 11 m e m x k r        ,       11 11 mm e mm y k r        ,   1 m e m zr   can be determined with investigation of the properties of eigenvalues of the characteristic equation (13) All three roots can be easily computed in Mathematica ® by using function Solve, however the displayed result is comprehensive. In the case of Lorenz system (6), Sparrow has been found the critical value     3 1 ,1 b c b rb         , where the subcritical Hopf bifurcation occurs [4]. All three roots of characteristic equation have negative real part if rr c, then two complex eigenvalues have positive real part and equilibria become unstable. Such kind of stability analysis in the case of system (5) is too complicate. Instead of this, roots, obtained by solving Eq. (13) with Mathematica ® can be evaluated numerically and the sign of real part of each root can be checked separately. If this test finds two complex roots with positive real parts, the instability occurs. Thus, the computer test is proposed, which consists from increasing the value of parameter r in the range r>0, computing corresponding eigenvalues by solving Eq. (13) and looking for the smallest value of r, where two complex roots have positive real parts. In the region of instability, unstable periodic and chaotic motions, respectively, are possible. Chaotic motions are significantly distinguished from unstable periodic motions because they are disordered and completely unpredictable. Until yet, we have discussed the saturation state and the chaotic state, but not the steady state. According to the definition, the steady state of the nonlinear supply chain occurs if the behavior of the supply chain stays same for the long enough time. In other words, the system is in the steady state if the same pattern is repeated in the long enough period of time. Having known solutions of Eqs. (5), the trajectories of the system can be shown in the phase space portrait. Phase space portrait of the system in the steady state takes the form of closed loops. The system is said to be in a regular state, if the trajectory of the system takes a cyclic pattern, which facilitates the future prediction. Phase space portraits of the chaotic state are very sensitive on the change of initial conditions and show disordered patterns so that a future prediction is impossible. Parameters m, r, k of the Lorenz-like model (5) in reference model of nonlinear supply chain must have fixed, but realistic values, because only in such a case the supply chain operates in normal operation conditions. Parameter m in Eq. (5), may be interpreted as percentage value. If m is fixed to have percentage value m = 15 %, then the retailer places orders assuming that least 15 % of his orders are satisfied by the distributor. This percentage value of course is relatively low, but it is chosen to allow the study of the steady state behavior in the supply chain model without disturbances. The parameter r means the rate of the information distortion at the distributor and can be interpreted as the percentage value in a natural way. When r = 30, this means that the distributor received the information about the order quantity from the retailer with 30 % distortion rate. Finally, when the value of parameter k is chosen to be k = ⅓, then the safety stock coefficient determines, that the manufacturer produces one third more (approximately 33 % more) than the order received to avoid uncertain situations. It is worth to note, that chosen parameter values m = 15, r = 30, k = ⅓ in the nonlinear supply chain reference model differ from the values σ = 10, r =28, b = 8/3 in the classic Lorenz atmospheric convection model [6]. Considering chosen values of parameter, we first compute the critical value r c. To do this, we compute eigenvalues of Eq. (13) and check them, when their positive real parts appear. By performing the computer test, we obtain r c = 26.9105. Because the chosen value of parameter r is higher than the critical value, r=30 > r c=26.9105, we cannot find any stable pair of equilibrium points. Therefore, the saturation state in the system (5) without disturbances does not exist. As is mentioned previously, in the range r>r c both unstable periodic as well as chaotic solutions appear. With chosen values of parameters, we obtain the steady state periodic solution after the very short transient phase dies. The time history of the phenomenon is depicted in Figure 2.a and the corresponding phase portrait in the Figure 2.b. Both the time responses as well as the phase diagrams show that three players of the supply chain fulfilled their objectives without hindering the reaching of objectives of other participants. Both time responses as well as phase portraits are computed by using the multistage homotopic perturbation method and drawn in figures with continuous line. The method is still presented in [6] and is compared in this paper with results of numerical integration of the system (5) using the Runge-Kutta method. Results of the numerical integration are shown in the Figures 2. a and 2. b with dashed curves and compare very well with results of the multistage homotopic perturbation method. MODELLING NONLINEAR SUPPLY CHAINS USING METHODS OF SYSTEM DYNAMICS | 92 Time histories on Fig. 2.a show the periodic nature of the steady state and from corresponding plots in Fig. 2.b it can be seen, that all phase portraits are cyclic. Thus, the steady state of the supply chain is a regular state. A direct proof of this conclusion can be made by construction of the so called Poincare map [7], which in this case is reduced in a point. Because the claim is evident from figures shown, the proof is omitted, however its sketch can be found in [8]. 5. States of nonlinear supply chain in the presence of external disturbances Supply chain cannot stays in the regular state, studied until now, for a long time. Such behavior is the consequence of the inevitable disturbances, which can affect any level of the supply chain. The time dependence of disturbances can be random, but to simplify the discussion a bit, we assume the presence of external seasonal (periodic) disturbances only, which can affect the each level of the supply chain additively in accordance with Eqs. (7). By using disturbances in the form             1 2 3 10cos 5 , 5cos 10 , 10cos 10 d t t d t t d t t    , but retaining the values of parameters 15, m  30, r  1 3 k  to be equal as in the reference model, we are surprised, that there exists a saturation state, which cannot be found in the reference model. With vanishing derivatives on the left hand side of Eqs. (7) and specified harmonic disturbances, we can find time courses of variables x(t), y(t) and z(t), which can annulate disturbances d 1(t), d 2(t) and d 3(t) on the right hand side of Eqs. (7). The desired functions are computed in Mathematica ® by applying the command Solve[{m*y- (m+1)*x+d 1==0, r*x-y-x*z+d 2==0, x*y-(1- k)*z+d 3==0},{x,y,z}]. The displays of computed functions are comprehensive, however plotting of time histories is a simple task, which is shown in Figs. 3.a,b. From these plots it follows, that all functions x(t), y(t) and z(t) have nonzero DC components (their derivatives of course are zero), which are coordinates of the saturation point in the Fig. 3.c. A simple computation in Mathematica ® , performed on obtained solutions x(t), y(t) and z(t), reveals, that coordinates of the saturation point in the Fig. 3.c are x e=4.21029, y e=4.49435, z e=29.0757. As mentioned in previous section, the saturation state of the supply chain is dangerous, because further changes cannot be effectively computed, when the supply chain enters in this state. Fig. 2. Steady state of the nonlinear supply chain without external disturbances. a) Time histories of retailer order's, quantity of products, which can be supplied by distributor and quantity of products made by manufacturer. b) phase portraits of the steady state. ________ results, computed by multistage homotopic perturbation method, - - - - results, computed by Runge-Kutta method. a) b) c) Figure 3. Saturation state of the nonlinear supply chain in the presence of disturbances. a) time histories of x( t) and y( t). _______ time history of variable x( t), - - - - time history of variable y( t). b) time history of variable z( t). c) position of fixed (equilibrium) point in the phase space. ANALI PAZU, 3/ 2013/ 2, str. 87-93 Rudolf PUŠENJAK, Maks OBLAK 93 | Besides the passage into the saturation state, nonlinear supply chain can exhibit chaotic motions in the presence of external disturbations. Chaotic motions in the case of Lorenz attractor are possible, if the system parameters take the characteristic values. To show the phenomenon of chaos in nonlinear supply chain, the values of parameters m=30, r=28.5 and k=⅓ are chosen and external disturbances are changed to follow functions             1 2 3 0.56cos 3 , 20cos 5 , 50cos 10 d t t d t t d t t    in dependence on time. The chaotic responses of the supply chain are again computed by the multistage homotopic perturbation method and compared in time as well as in the phase space with results of numerical integration of Eqs. (7) by using Runge-Kutta method. Time histories are shown in Fig. 4.a and the corresponding phase portraits are depicted in the Fig. 4.b, respectively, where results of multistage homotopic perturbation method are plotted by continuous line and results, obtained by Runge-Kutta method are shown with dashed lines. The comparison of computed results by application of both methods reveals an excellent agreement also in this chaotic case. 6. Conclusion The paper treats the modelling of the nonlinear supply chain, which consists from a manufacturerer, a distributor and a retailer by derivation of the Lorenz-like system of equations, which can exhibit chaos. The computer test to find the critical value of the parameter r, which determines stability behavior of the supply chain, is envisaged. The steady state of the supply chain without disturbances, which correspond to the regular state, is analysed and explained. Conditions for occurence of the saturation state in the supply chain, which is subjected to the periodic external disturbances, are found and the corresponding saturation state is effectively computed. In the supply chain with periodic disturbances, the chaotic state is also analysed. The time histories and phase portraits of investigated states of supply chain are computed by multistage homotopic perturbation method and compared with computations, performed by using Runge-Kutta method of numerical integration of governing differential equations. It is found, that both method are in an excellent agreement. References 1. Cristopher, M. Logistic and supply chain management: strategies for reducing costs and improving services. Financial Times, Pitman Publishing, London, 1992. 2. Forrester, J. W. Industrial dynamics. Portland: Productivity press, 1961. 3. Pušenjak, R., Oblak, M., Fošner, M. Osnove modeliranja dinamičnih procesov v logistiki II. Celje: Fakulteta za logistiko, 2013. http://blend.fl.uni-mb.si. 4. Sparrow, C. The Lorenz Equations: Bifurcations, Chaos, and Strange Attractors, Springer-Verlag New York Inc., New York, 1982. 5. Lorenz, E. N. Deterministic non-periodic flow, Journal of the Athmospheric Sciences 29, no. 2, 131—140, 1963. 6. Pušenjak, R., Oblak, M., Lipičnik, M. Competitive predator-prey systems with time-dependent coefficients: a multistage homotopy perturbation analysys = Tekmovalni sistemi plenilec-plen s časovno odvisnimi koeficienti: večstopenjska homotopsko perturbacijska analiza, Anali PAZU, vol. 2, no. 1, 22-31, 2012. 7. Wiggins, S. Global Bifurcations and Chaos: analytical methods, Springer-Verlag, New York, 1988. 8. Pušenjak, R. R., Oblak, M. M. Incremental harmonic balance method with multiple time variables for dynamical systems with cubic non-linearities, Int. j. numer. Methods eng., vol. 59, iss. 2, 255-292, 2004. a) b) Figure 4. Chaotic state of nonlinear supply chain. a) time histories of supply chain variables x(t), y(t) and z(t). b) phase portraits of the chaotic state. ________ results, computed by multistage homotopic perturbation method, - - - - results, computed by Runge-Kutta method.