Revealing instabilities in a generalized triadic supply network: A bifurcation analysis

Supply networks have become increasingly vulnerable to risks in the past decades caused, for instance, by increased demand volatility, reduction of buffers, increased outsourcing, and globalization.^1–3 A number of supply networks have experienced major supply crises in the recent years. Ericsson suffered a loss of around $$400$ million in 2000 due to a fire at its chip supplier’s facilities,⁴ the 2002 strike at a California port resulted in widespread retail shortages,⁵ and the 2011 earthquake in Japan caused a total damage of $$235$ billion to the Japanese economy.⁶ Disruptions can be caused by a variety of external causes, ranging from singular localized events, such as a fire at a supplier’s facility, to systematic shortages of raw materials and large-scale disasters.^3,7,8

In addition to these external causes, supply problems can also develop internally due to the dynamics of supply networks.^9,10 The dynamics of material flows in supply networks have been investigated in the extant literature, focusing on two aspects: (a) the bullwhip effect,^11,12 that is, the amplification of the variation in the order quantities, and consequently in the inventory levels, as one moves up in the supply chain,^9,13,14 and (b) the stability of inventory dynamics.^15–18 Our focus here is on the stability of inventory dynamics, which characterizes whether considered steady-state inventory levels can be sustained over time, i.e., “the ability to ensure continuity.”¹⁹ The lack of stability leads to divergence of inventory levels from their planned steady-state levels and oscillations, which are costly due to uncontrolled stockovers and stockouts.²⁰ If instabilities are encountered, one can expect firms to adapt their policies over time, for instance, by altering the inventory replenishment periods and investing into capacity to counteract instability. Although this can bring the network over time to stable operational conditions, the instability would still impact the firms transiently, leading to stockouts and/or excess inventories before corrective actions are taken and they become effective. Furthermore, future changes in demand and supply might again drive the network to instability. Therefore, it is important to investigate the stability properties of supply networks.

Stability of small serial supply chains has been investigated in the extant literature using control theoretic models that employ mainly linear relationships between material flow rates and state variables, i.e., inventory levels and demand forecasts.^15,17 The small number of studies that have considered non-linearities demonstrate, for instance, that supply chains are subject to border-collision bifurcations, chaos, and amplification of chaotic dynamics.^21–24 Although the extant literature provides valuable insights into the dynamic behavior and stability of supply networks, obtaining a comprehensive account of the general instability potential of supply networks is still challenging. Control-theory models that concern supply network stability focus on specific ordering and production policies, which are represented by specific functional forms. This is restrictive due to several reasons. First, default management policies of different firms differ considerably and it is challenging to gather information about these policies across the supply network. Even the prime entities of a supply network typically lack the full information available, due to, for instance, conflicts of interests between different stakeholders and the difficulty of integration in supply networks.^25,26 Therefore, it is in general infeasible to precisely know the management policies and hence the associated functional forms related to the control policies of all firms in a supply network. Second, management policies may deviate in practice from the prescribed default policies under capacity and material availability limitations.^27–30 For instance, behavioral effects such as hoarding^31,32 lead to deviations when the network is under strain.

The limitations of reliance on specific functions have been encountered in various other fields when considering network stability, most notably in ecology. This problem has been addressed by generalized modeling (GM),^33,34 which has been applied to several other areas, including metabolism,^35,36 cell signaling,³⁷ endocrinology,³⁸ and history.³³ A generalized model employs unspecified instead of specified functions, and provides a systematic method for analyzing the stability of a large spectrum of specific models within a single model. The real appeal in a supply network context is that a generalized model can be set up solely using the topology of the network, which captures the inputs and outputs of the firms. With this basic information, a generalized model can be set up, which is based on a set of intuitive generalized parameters that characterize the time-scales of change in the system and the elasticities, i.e., the sensitivity of functions to state variables in the steady state.

Two common and complementary approaches exist for analyzing generalized models: (a) statistical ensemble approach³⁴ and (b) bifurcation analysis.³³ In the statistical ensemble approach, the generalized parameters are initialized randomly and the stability of the network for each initialization of these parameters is determined. This method constitutes a fast and efficient computational approach to understand the overall impact of generalized parameters on stability, which is especially useful for large and complex networks.³⁴ It can also be used to identify the most crucial components (e.g., firms and parts) of the network for stability.³⁹ Bifurcation analysis is a standard method of the dynamical systems theory⁴⁰ and it complements the statistical ensemble approach by explicitly calculating the conditions that partition the generalized parameter space into regions with different stability characteristics. A bifurcation is a qualitative change in the behavior of a system. For instance, Hopf bifurcations mark the onset of oscillations and saddle-node bifurcations typically mark the onset of bistability, hysteresis, and catastrophic shifts.⁴¹ 

The central research question addressed in this study is to investigate whether certain topological structures induce instabilities in supply networks and is to characterize the dynamical consequences. We set up a generalized supply network model, similar to Demirel et al.⁴² In contrast to Demirel et al.,⁴² which investigates large network structures using the statistical ensemble approach, we perform an in-depth bifurcation analysis of basic supply network motifs, namely dyads, serial supply chains, and a triadic motif. We consider only feed-forward flows for a single end product in these networks. Hence, the triadic motif forms a cycle in the underlying undirected graph, but not a directed cycle (Fig. 1).

Using this framework, this study makes three major contributions to the literature. First, we show that topology-induced instabilities do not emerge in dyads and serial supply chains, while the triadic motif is subject to these instabilities (Sec. III). This implies that the triadic motif is one of the smallest supply networks that is able to generate topology-induced instabilities. Second, using bifurcation analysis, we investigate the impact of production rate elasticities on the stability of the triadic motif (Sec. IV). We show that static (saddle-node type) and oscillatory (Hopf) bifurcations take place in meaningful generalized parameter regions. In addition, we identify the intersections of these bifurcations, which are codimension-two bifurcation points, since these can indicate the presence of global bifurcations. To identify global bifurcations, we conduct further analysis on an illustrative specific model. Third, the managerial implications of the different identified dynamical scenarios are discussed.

We start by describing the triadic supply network and reviewing the generalized model for this network, which was introduced in Demirel et al.⁴² We then briefly summarize the concepts of stability and bifurcations (Sec. II C).

The triadic network consists of three firms, one prime company and two suppliers, that conjointly manufacture a single final product (Fig. 1). Firms are represented by nodes, while the flows of products are represented by arcs, with arrows pointing in the direction of the material flows that correspond to the manufacturing of parts, shipping of parts between firms, and sale of the final product to the external market. Products $1$ and $2$ are considered as external parts, thus their manufacturing is not described by the model. Firm I manufactures products $3$ and $4$ that both use the external part $1$⁠. Firm II manufactures product $5$⁠, using part $3$⁠, which is supplied by firm I, and external part $2$⁠. The prime company (firm III) assembles the final product $6$ from products $4$ and $5$⁠, supplied by firms I and II, respectively. The final product $6$ is then sold to the external market.

We make the following assumptions in formulating the model for the dynamics of the triadic supply network. First, we consider a simple bill-of-materials (BOM) structure that involves only the products $1$ to $6$⁠, as described above. Although the same triadic network can be used to manufacture multiple end products with more complicated BOM structures, we focus on this simple BOM for a single end product since the resulting model captures the inter-relatedness between the first-tier suppliers of a focal firm, i.e., suppliers I and II of the firm III, which is a common property of real world supply networks,^43–45 but is at the same time sufficiently low dimensional that enables bifurcation analysis. Second, the external demand for the final product is assumed to be fixed, which is reasonable for products with stationary demand and low variability. Third, we assume continuous and instantaneous material flows. We consider total system inventory levels of products $1−6$ and presume that the parts are produced and transported instantaneously, i.e., time delays are not captured, and that there is no batching of purchase and production orders. Therefore, we do not aim to capture instabilities that might be caused by discrete delays^17,18 and batching.¹³ We instead develop a high-level phenomenological model that focuses on instabilities that are induced by the structure of the network.

The generalized model describes the dynamics of the inventory levels of products, i.e., $P1$ to $P6$⁠. The inventory level of product $i$⁠, $Pi$⁠, increases due to processes where product $i$ is an output and decreases due to processes where $i$ is an input. To describe the changes in inventory levels, we introduce material flows that result in a product $i$ and consider that their flow rate $Fi$ depends on the inventory levels of (a) the parts used for the production, (b) the product itself, and (c) other products that are manufactured by the same firm. Using these material flows and approximating the stock levels $Pi$ by continuous numbers, which is a valid assumption if the material flow rates $Fi$ are high enough, we can express the changes in the inventory levels of the products by the ordinary differential equations

which can be rewritten in the matrix form,

where $P=(P1,…,P6)$ and $F(P)=(F1,…,F7)$ are the vectors of inventory levels and material flows, respectively, and

is the stoichiometric matrix of the material flows.

We have set up a model describing material flows on the triadic supply network. However, we have not specified functional forms for the flows. In fact, a whole spectrum of linear and nonlinear policies for inventory replenishment and order satisfaction can be captured by the model. Instead of focusing on a single policy, we now apply generalized modeling and capture the whole spectrum of alternative functional forms, characterized by a set of generalized parameters.

We now consider the stability of a steady state with stock levels $Pi∗$ and flow rates $Fi(P∗)=Fi∗$ where $dP/dt|∗=0$⁠. Stoichiometric balance in the steady state $dP/dt|∗=0$ requires that flows have to balance in the steady state, e.g., if firm III uses one unit of product $4$ and one unit of product $5$⁠, it has to produce exactly one unit of product $6$⁠. Thus, the number of free flow parameters reduces to one, that is $F2∗=F3∗=F4∗=F5∗=F6∗=F7∗:=F∗$ and $F1∗=2F∗$⁠. Focusing on realistic non-vanishing steady states, i.e., $Pi∗>0$ for all $i$⁠, and by normalizing all quantities with respect to their respective steady-state values, the system is transformed into a set of generalized variables,

with $p=(p1,…,p6)$⁠. Using this transformation, all variables evaluate to one in any steady state, i.e., $pi∗=1$ and $fi∗(p)=1$ for all $i$⁠. With these generalized variables, Eq. (1) can be rewritten as

where $f=(f1,…,f7)$⁠. Note that $Λ$ denotes the turnover rates of products (⁠$∑jFj∗Nij/Pi∗$⁠) times the branching factor for each outflow that determines how much of the product goes in the corresponding direction, which constitute the first type of generalized parameters.

Stability is in general associated with the response of the system to perturbations. Here, we are concerned mainly with the stability of steady states against small perturbations, which is captured by the concept of local asymptotic stability. If a small perturbation triggers a departure from a steady state, the state is said to be unstable. If the system instead returns to the original steady state, the state is said to be locally asymptotically stable. Local asymptotic stability of a steady-state solution is characterized by the eigenvalues of the Jacobian, $J$⁠, of the dynamical system. If the real parts of all eigenvalues are negative, the steady state is stable; otherwise, it is unstable.⁴¹ In what follows we refer to the eigenvalue with the largest real part as the leading eigenvalue.

The Jacobian $J$ of the normalized system is defined as

where

with $Θfp$ containing the logarithmic derivatives of the original material flows evaluated at the steady state:

The parameter $fipj$ measures the elasticity of the flow $Fi$ to the stock level $Pj$ at the steady state.^46,47 Elasticities constitute the second type of generalized parameters and are key to the stability analysis of generalized models.³³ Elasticities can be interpreted as the exponent of a power-law fit to the flow rate close to the steady state, i.e., $Fi∝Pj(fipj)$⁠. For instance, $fipj=1$ describes a linear relationship between $Fi$ and $Pj$⁠, while $fipj=0$ states that the material flow is independent of the stock level. In our generalized supply network model, there are three different types of elasticities: elasticity to supply, elasticity to inventory level, and elasticity to co-production.

Elasticity to supply $fi,Spj$ characterizes the dependency of the production rate of a product on the availability of parts $j$ that are required for its production. We can distinguish between two cases where it is (a) positive and (b) zero. If the inventory level of parts used for the manufacturing of the product is insufficient to meet a desired production rate, the production rate needs to be adjusted according to the inventory level of parts used, i.e., $fi,Spj>0$⁠. For instance, a situation where all supply parts that become available in stock are used directly for production is represented by $fi,Spj=1$⁠. Furthermore, production rates are determined by the inventory levels of ingredients if they are deteriorating, which also leads to a positive elasticity to supply. On the other hand, if supplied parts are abundant and non-deteriorating, the production is not generally influenced by material availability, i.e., $fi,Spj=0$⁠. Elasticity to supply parameters are $f3p1$⁠, $f4p1$⁠, $f5p2$⁠, $f5p3$⁠, $f6p4$⁠, $f6p5$⁠, and $f7p6$⁠.

Elasticity to inventory level $fi,Ipi$ captures the sensitivity of the production rate to the inventory level of the output product. Since the inventory level refers to the total amount of products held at both the buyer and the supplier facilities (due to the no-delay assumption), the elasticity to inventory level captures the impact of both types of inventories. A high inventory level at the buyer’s site will lead to a lower demand for the part and hence orders. Similarly if the inventory level at the producer’s end is high, its production rate will be decreased in order to achieve a target inventory level over some replenishment period. We note that although high levels of inventories can be built on purpose due to economic reasons, this is not against the self-inhibiting role of inventory control, which is a common property of inventory control models.^48,49 Instead, this corresponds to higher steady-state levels of inventories, which appear in the first type of generalized parameters, i.e., turnover rates. An exception, where the order amount can lead to an increase in the production rate despite an increase in the inventory level, is the impact of batching around the discontinuity points, which is not considered here due to the continuous flow assumption. Therefore, we consider the elasticity to inventory level to be always non-positive, i.e., $fi,Ipi≤0$⁠. A limiting case is the absence of elasticity to inventory level, $fi,Ipi=0$⁠, which takes place when the production rate is fully dictated by a constraint on production such as part availability or production capacity. Elasticity to inventory level parameters are $f1p1$⁠, $f2p2$⁠, $f3p3$⁠, $f4p4$⁠, $f5p5$⁠, and $f6p6$⁠.

Elasticity to co-production $fi,Cpj$ describes how the production rate of one part produced by a given supplier is impacted by the inventory level of another part produced by the same supplier. In contrast to the other elasticities, it can be positive or negative. On the one hand, an organization may use the same limited resources, such as available machinery and personnel, to produce different parts and hence may have capacity constraints. Thus, an increase in the production rate of one part may need to be compensated by a decrease in the production rate of another part. This is captured by a positive co-production elasticity, $fi,Cpj>0$⁠. On the other hand, the production of one part can induce the production of the other, which is described by negative co-production elasticity, $fi,Cpj<0$⁠. For instance, one product can be a by-product of the other and, therefore, while the main product is manufactured, the by-product is accumulated. In contrast to negative and positive elasticities to co-production, if the production of different parts does not share any common resources or the production rate is not set by constraints, the two production rates are independent, hence $fi,Cpj=0$⁠. Co-production elasticity parameters are $f3p4$ and $f4p3$⁠.

Annotating the type of elasticity to the corresponding entries in Eq. (7), we arrive at the final notation for $Θfp$⁠,

Before analyzing the dynamics of the generalized triadic supply network, we provide a brief introduction to the concepts of bifurcation theory, while mathematical definitions are provided in  Appendix.

Bifurcations are qualitative changes in dynamical behavior, such as the emergence of new steady states and changes in the stability of steady states, as model parameters are changed. Bifurcations can be broadly classified into local and global bifurcations. Local bifurcations involve changes in the dynamical behavior around a steady state. In contrast, global bifurcations are not confined to the local neighborhood of steady states and influence the dynamics over distant regions of state space.

The simplest and most common bifurcations are codimension-one local bifurcations that require the change of a single model parameter only. Depending on the nature of leading eigenvalues, codimension-one local bifurcations can be divided into two different types: saddle-node type and Hopf. The leading eigenvalue is zero at a saddle-node type bifurcation, while a complex-conjugated pair of purely imaginary leading eigenvalues exists at a Hopf bifurcation (see  Appendix Fig. 7). In a saddle-node type bifurcation, different steady states collide, exchange stability, or annihilate each other. If a steady state inventory level is destabilized at a saddle-node bifurcation, even small deviations from this level will lead to exponentially diverging trajectories. Therefore, this will be characterized by quickly depleting stocks or building up inventories for different products in the network. For instance, instability can lead to the uncontrolled build up of part 3 at the expense of part 4 by the firm I [as shown in Fig. 4(a)]. In a Hopf bifurcation, a limit cycle emerges while the steady state changes its stability. Hopf bifurcations are therefore related to oscillatory behavior. The resulting dynamics are similar to the oscillations observed together with the bullwhip effect.⁵⁰ 

Codimension-two local bifurcations, such as Takens-Bogdanov, Gavrilov-Guckenheimer (also known as fold-Hopf), and cusp bifurcations, are the points where codimension-one bifurcations change their qualitative behavior, i.e., they are bifurcations of bifurcations. For instance, cusp bifurcations occur at the intersection of two saddle-node bifurcations, while a Takens-Bogdanov bifurcation happens at the intersection of a Hopf and a saddle-node bifurcation. In general, higher codimension bifurcations have implications for complex phenomena such as chaotic dynamics,⁴⁰ which has also been observed in models of supply chains.^21–24 

Before investigating the triadic motif, we consider a simple dyad, which can be easily extended to a serial supply chain. The triadic motif shown in Fig. 1 can be reduced into a dyad by removing supplier II and products $2,3,$and$5$⁠. The generalized model of the dyad is given by

and

Saddle-node bifurcations occur when the leading eigenvalue of $J$ is zero. Eigenvalue $λ$ satisfies $det(J−λI)=0$⁠. Therefore, the saddle-node bifurcation occurs when $det(J)=0$⁠, which for the dyad yields

The sum of all terms is always non-positive and $det(J)=0$ is fulfilled if and only if all elasticities in at least one of these sets {$f1,Ip1,f6,Sp6$}, {$f1,Ip1,f7,Sp6$}, {$f1,Ip1,f7,Sp6$}, {$f4,Ip4,f6,Sp6$}, {$f4,Ip4,f7,Sp6$}, and {$f6,Ip6,f7,Sp6$} are equal to zero. Therefore, saddle-node bifurcations occur only at the boundaries of the considered parameter space.

Hopf-bifurcations occur if the resultant $R$ of the Jacobian evaluates to zero, under the condition that the Hopf-number $χ$ is negative. The Hopf-number $χ$⁠, which can be found by the resultant method, is negative when the eigenvalue is complex (see  Appendix).

The resultant $R$ is given by

All terms are non-negative for all possible meaningful values of elasticities. Therefore, the resultant $R$ is also non-negative and reaches zero, if and only if, specific triplets of elasticities are zero and thus in certain parts of the boundary between the meaningful considered region of parameter space and the non-considered one. We note that the leading eigenvalue is zero on the boundary, as shown above. Both the resultant $R$ and the determinant $det(J)$ vanishing at these certain regions of the boundary indicate that Takens-Bogdanov bifurcations occur in these regions (see  Appendix). This is further confirmed by the fact that the Hopf number $χ$⁠, which is given by

evaluates to zero at these points.

All the bifurcations, saddle-node or Takens-Bogdanov, occur at the boundary between the considered and the non-considered parameter spaces. Therefore, the stability, hence the dynamic behavior, does not change within the considered parameter space. This indicates that we can check the stability of an arbitrary interior point of the considered parameter space and all other interior points must have the same stability property. Accordingly, we evaluated the leading eigenvalue at an interior point and found that it is negative. Thus the dyad is asymptotically stable in the interior of the considered parameter space. The leading eigenvalue reaches zero at the boundary, the system becomes neutrally stable, and a bifurcation point occurs, as discussed above. The bifurcation leads to instability in the non-considered parameter space, which corresponds to unrealistic or erroneous policies. In summary, the dyad does not change its dynamical behavior within the considered meaningful parameter space and is asymptotically stable except on the boundary.

These calculations can be extended to serial supply chains of different sizes, yielding similar results. In the following, we show that if elasticities of the same kind have identical values (i.e., $fi,Ipj=fI$⁠, $fi,Spj=fS$⁠), a supply chain of arbitrary size cannot become unstable for any given meaningful set of parameter values. The Jacobian $J$ of a serial supply chain of length $N$ is given by

which is a tridiagonal Toeplitz matrix.⁵¹ The eigenvalues $λk$ (⁠$k=1…N$⁠) are given by

Using meaningful parameters, the first term of Eq. (16) is always non-positive. Therefore, the leading eigenvalue, $λl$⁠, is obtained when the cosine function reaches its maximal value of $1$⁠:

The leading eigenvalue is negative or zero for all possible elasticities. Therefore, a supply chain of arbitrary length, with elasticities of the same type having identical values, can never become unstable in this model.

We now present the bifurcation analysis of the generalized triadic network. The model contains in total $17$ free parameters. To reduce the dimensionality of the analysis, we first assume that all elasticities of the same kind have identical values and analyze the bifurcations in this homogeneous system (Sec. IV A). Afterwards, to study the impact of the triadic interactions in more detail, we vary the two elasticity to co-production parameters independently, while keeping the other elasticities fixed (Sec. IV B).

For simplicity, we start by assuming that all parameters of the same kind have identical values, i.e., $fi,Ipj=fI$⁠, $fi,Spj=fS$⁠, $fi,Cpj=fC$⁠, and $Pi∗=P∗$⁠, for all $i$ and $j$⁠.

As above, saddle-node bifurcations are given by the solutions of $det(J)=0$⁠, which leads to

Avoiding the trivial case $fS=fI=0$⁠, this expression has three solutions,

each defining a saddle-node bifurcation surface, which are plotted in Fig. 2(a) (red and pink). Parameter values within the region that is bounded by the three surfaces and the sets $fS=0$ and $fC=0$ represent stable steady-state solutions.

The first saddle-node bifurcation surface SN1 corresponds to the insensitivity of the production rate to inventory levels. This also implies that positive elasticities to inventory level (that is the increase of production rates with the increase of inventory level) destabilize the system, as would be expected. The surfaces SN2 and SN3 involve, in addition to the elasticity to inventory level $fI$ the elasticity to co-production $fC$⁠. Both saddle-node bifurcations, SN2 and SN3, are felt the most in supplier I, which acts as both a first and a second tier supplier in the triadic network. This can be seen in the eigenvector corresponding to the leading eigenvalue (not shown), since products $3$ and $4$⁠, which are both products of supplier I, have the largest absolute value. The instability introduced by passing the saddle-node bifurcation surface SN2 concerns negative elasticity to co-production, $fC<0$⁠, while the instability due to passing SN3 occurs for positive elasticity to co-production, $fC>0$⁠.

Negative elasticity to co-production occurs when the production of part $3$ induces the production of part $4$⁠, or vice versa. In this case, the stabilizing influence of the elasticity to inventory level, or in general the stabilizing effect of the self-regulated stock control, is annihilated by the destabilizing influence of the elasticity to co-production, thus the network is unstable for $|fC|>|fI|$⁠.

Positive elasticity to co-production occurs when the production rates for parts $3$ and $4$ at supplier I are limited by the firm’s finite production capacity. As before, this can counteract the stabilizing influence of the other elasticities and destabilize the system. If the production is not limited by part availability, i.e., $fS=0$⁠, the system is unstable for $|fC|>|fI|$⁠. However, for larger $fS$ and $|fI|$⁠, the elasticity to co-production that is necessary to destabilize the network increases [see surface SN3 in Fig. 2(a)].

The Hopf bifurcation surfaces are given by a resultant $R6(fI,fS,fC)$ of zero and a negative Hopf-number $χ$⁠. We solve this numerically by using a triangulation method⁵² [see Figs. 2(a) and 2(b), blue surface HB]. Hopf-bifurcations occur in a very narrow parameter region, namely for $fI≈0$ and $fC<0$⁠. It is strongly felt in supplier I, since product $1$ has the highest contribution to the leading eigenvector (not shown). Hopf bifurcations are related to, at least temporary, oscillations of inventory levels as discussed in more detail below.

As stated in Sec. II, higher codimension bifurcations take place at the intersections of separate lower codimension bifurcations. We observe two types of codimension-2 bifurcations in the generalized supply network model. First, two cusp bifurcations occur at the intersection of the saddle-node bifurcation surfaces SN2 and SN3 with SN1 [Fig. 2(a)]. Second, Takens-Bogdanov bifurcations emerge at the intersection of the Hopf bifurcation surface HB with the saddle-node bifurcation surface SN1 [Fig. 2(b)]. Takens-Bogdanov bifurcation has implications for global bifurcations as discussed in Sec. IV B.

The analysis above indicates that the elasticity to co-production, and thus the inter-dependencies between the production rates of the two parts that firm I produces (parts 3 and 4), has a strong effect on stability. These inter-dependencies are captured by the parameters $f3,Cp4$ and $f4,Cp3$⁠. Due to their importance, we allow these two parameters to take different values and compute the bifurcation points. For simplicity, we keep, as before, uniform values for the remaining parameters, i.e., $fi,Ipj=fI$⁠, $fi,Spj=fS$⁠, and $Pi∗=1$⁠.

For this detailed analysis of the influence of the co-production elasticities, we are interested in global bifurcations as well as local bifurcations. Global bifurcations can emerge from high-codimension local bifurcations. Since the Jacobian of the generalized system only allows to compute local bifurcations, we complement this approach by using a specific model to illustrate typical dynamical scenarios with global bifurcations. Specifically we take,

This specific model has the advantage that all functions employed are parameterized directly in terms of the elasticity parameters (⁠$f3,Cp4$⁠, $f4,Cp3$⁠, $fI$⁠, and $fS$⁠). In addition, it has a steady state at $P∗=1$⁠, and the Jacobian of this model at $P∗=1$ is identical to the Jacobian of the generalized model. Therefore, all local bifurcations found in the specific model occur for exactly the same parameter values in the generalized counterpart.

We use this specific model and compute higher-order derivatives to determine the nature of Hopf bifurcations (i.e., subcritical or supercritical) and to identify global bifurcations. The results in this section are obtained using the numerical continuation software Auto.⁵³ 

We note that while local bifurcations of the specific model are also found in the generalized model, global bifurcation curves are particular to the specific model. However, qualitatively equivalent global bifurcation curves must occur for a large family of specific functional forms, which is implied by the existence of high codimension bifurcations.⁴⁰ 

We first consider a scenario with $fS=1$ and $fI=−0.01$ (Fig. 3). This represents, for instance, a strained steady state where part inventories are scarce so that the desired production rate of the downstream firm cannot be achieved and the available parts are immediately used for production. This scenario involves several bifurcations (Fig. 3), namely two Hopf bifurcation curves (HB1 and HB2, blue lines), two saddle-node bifurcation curves (SN1 and SN2, red lines), and two global bifurcations (HC1 and HC2, purple lines). These bifurcation curves partition the $f4,Cp3$-$f3,Cp4$ parameter space into six regions with different dynamical properties. The steady-state solution is unstable in the regions I-IV (white background), but locally stable in regions V and VI (gray background). Typical trajectories of inventory levels $P3$ (red) and $P4$ (black) for each region are plotted in Fig. 4, with initial conditions close to (solid) and further away from (dashed) the steady state.

If both elasticities $f4,Cp3$ and $f3,Cp4$ are negative (region I), the system is unstable. The same holds true if both elasticities have high positive values (region II). However, if one elasticity has a high enough positive value and the other is negative, the system becomes stable (regions V and VI). In unstable systems (regions I-IV) the stock levels escape from the steady state, which leads to either an overproduction or scarcity of products. In contrast, in stable systems (region V and VI) trajectories that start close to the steady state approach it (see Fig. 4 solid lines). Trajectories that start further away from the steady state (Fig. 4 dashed lines) are affected by global bifurcations, as in Region V. Regions V and VI in Fig. 3 are separated from one another by the global homoclinic bifurcation curve HC1 that emerges from the Takens-Bogdanov point TB1. This point occurs at the intersection of the Hopf bifurcation curve HB1 with the saddle-node bifurcation curve SN1. Since the Hopf bifurcation HB1 is subcritical, it indicates that an unstable limit cycle emerges at HB1, which surrounds the stable steady state in Region V (Fig. 5). The unstable limit cycle defines a basin of attraction⁵⁴ for the focal stable steady state $P∗=1$⁠. In a sufficiently small neighborhood of the steady state $P∗=1$⁠, inventories are dynamically pushed towards $P∗=1$ by the unstable limit cycle [Fig. 4(e) solid lines]. Otherwise, they escape towards the system boundary [Fig. 4(e) dashed lines]. As the parameter $f3,Cp4$ increases, the limit cycle grows in amplitude [Fig. 5] and hence the basin of attraction of the stable fixed point $P∗=1$⁠. The limit cycle eventually collides with the homoclinic orbit at the global homoclinic bifurcation point HC1 that emerges from the Takens-Bogdanov point TB1. Therefore, the two types of dynamics, i.e., convergence to vs. divergence from the steady state, coexist in region V. Above the HC1 curve (region VI), the fixed point $P∗=1$ is globally attracting so that the system eventually recovers back to this state. Analogous to the scenario regarding the HC1 curve, the homoclinic bifurcation curve HC2 emerges from the Takens-Bogdanov point TB2. We note that, since the HC2 and HB2 curves are very close, the region between them is not labeled in Fig. 3.

We now analyze the impact of the parameters $fS$ and $fI$ on bifurcations in the $f4,Cp3$-$f3,Cp4$ parameter space. For small elasticity to inventory $fI$⁠, the elasticity to supply $fS$ has two major effects [compare Figs. 3, 6(a), and 6(b) where $fS$ is 1, 0.5, and 0, respectively]. First, with decreasing $fS$ both Hopf bifurcation curves get closer to the corresponding axis where one of the parameters $f4,Cp3$ and $f3,Cp4$ is zero [$fS=1$ in Fig. 3 and $fS=0.5$ in Fig. 6(a)] until the curves vanish for small $fS$ [Fig. 6(b)]. Therefore, the Takens-Bogdanov points approach the origin. Second, with decreasing $fS$ the saddle-node bifurcations align with the axis where either $f4,Cp3$ or $f3,Cp4$ is zero [Fig. 6(b)]. As a consequence of these two mechanisms, as elasticity to supply $fS$ is reduced, (i) the regions of stability in the $f4,Cp3<0$ — $f3,Cp4>0$ and $f4,Cp3>0$ — $f3,Cp4<0$ quadrants grow, (ii) the region of stability in the $f4,Cp3>0$ — $f3,Cp4>0$ quadrant shrinks, and (iii) the region of stability in the $f4,Cp3<0$ — $f3,Cp4<0$ quadrant is not much affected. We note that for $fI>0$ the two saddle-node bifurcation curves never intersect at the origin (⁠$f4,Cp3=f3,Cp4=0$⁠), which always remains stable. The elasticity to inventory level $fI$ shifts the saddle-node bifurcation curves. With increasing $fI$ the curves move away from the origin, while the Hopf bifurcation ceases to exist for sufficiently high $fI$ [compare Fig. 6(d) with Fig. 6(b), and Fig. 6(c) with Fig. 3].

We have used generalized modeling to perform an in-depth bifurcation analysis of a dyad, a serial supply chain, and a simple triadic supply network. Our focus has been to investigate the range of dynamic behaviors these systems can exhibit and to identify conditions under which these occur. Our study has revealed three main findings.

First, we have shown that dyads and serial supply chains are remarkably stable. In particular, we have shown that dyads are not subject to topology-induced instabilities and are stable for all considered elasticities and thus managerial policies. In addition, we have proved that this is also true for serial supply chains of arbitrary length, if all firms have identical policies. These findings might be surprising in times of accumulating reports of supply chain disruptions and an increasing awareness of supply chain risks.^55,56 However, we note that we only investigate topology-induced instabilities and excluded instabilities caused by other factors including time-delays, batching, forecasting, and external disruptions.

Second, we demonstrated that the triadic supply network exhibits a number of different topology-induced bifurcations, such as saddle-node type, Hopf, cusp and Takens-Bogdanov bifurcations, and also global homoclinic bifurcations. Thus, combined with our first finding, the triadic supply network is one of the smallest motifs that can produce topology-induced instabilities. The importance of such bifurcations was previously identified in the supply chain management literature.^21–24 We identified conditions under which the different kinds of instabilities occur, i.e., which managerial policies may facilitate them. Thus, it complements previous studies, which focused on replenishment periods and lead times.^21–24 In addition, our analysis also enables us to provide interpretations of these instabilities.

For instance, saddle-node bifurcations are divergent instabilities and denote the overproduction or the absence of products. This would eventually lead to the collapse of the supply network, in the absence of management intervention to revise and modify inventory ordering policies accordingly. If this is done gradually, the network would still stay close to an instability and thus be more vulnerable to small disturbances, such as malfunctioning equipment. Hopf bifurcations lead to, at least temporary, oscillations of inventory levels, leading to alternation between shortage and surplus of supply. We note that the observed oscillations involve similar dynamical behavior as in the context of the bullwhip effect.^{12,17,57–59} However, the underlying cause is different, since we exclude the necessary drivers of the bullwhip effect in our model, e.g., time-delays, batching, and forecasting errors, and focus on topology-induced instabilities. Furthermore, we show that a global homoclinic bifurcation emerges from the Takens-Bogdanov point. This leads to two separate dynamic regimes where the steady state inventory levels can be sustained only against small perturbations in one and additionally against large perturbations in the other. Therefore, under certain conditions minor problems can be fixed with existing policies but dealing with major problems, when encountered, would necessitate the modification of the policies.

Third, we have shown that inter-dependencies of production rates within the same supplier, described by elasticities to co-production, strongly influence the stability of the triadic supply network. In particular, we have shown that the triadic motif is stable only in a specific parameter range. Therefore, it is of high importance to be aware of such underlying inter-dependencies between the rates of production carried out by firms in the network. Lack of awareness increases the inherent risk in the network^60,61. This risk is especially high, if co-dependencies have the identical sign in the same supplier, since it might cause the divergence of stock levels. Such dependencies result, for instance, if both products are competing for a finite resource, where both co-production elasticities are positive. As shown in Sec. IV, these instabilities emerge in the co-producing suppliers and, therefore, it is of high importance for the prime entities to monitor these suppliers closely. As the concern is capacity limitations at co-producing suppliers, prime entities can be advised to provide them incentives for capacity investments.^62,63 Furthermore, considering that co-producing suppliers might have limited resources, prime entities can be recommended to invest into such suppliers, which constitutes a criterion for choosing candidates for strategic supplier development activities.^64,65 Furthermore, the elasticity to inventory has a stabilizing impact, which indicates that the firms can be advised to more quickly adjust the orders /production rates to compensate for deviations from target inventory levels. However, we note that the effectiveness of this policy can be inhibited by delays.

All of our main findings have direct implications for the management of supply networks, since triadic relationships are a common feature of many supply networks.^43–45 Together with the finding that serial supply chains are remarkably stable, while triadic motifs can promote instability, our study implies that supply chain management should closely monitor and manage these triadic motifs, as discussed above. This supports recent studies that also emphasize the influence of triadic interactions.^66–68 However, in real supply networks firms can adapt their policies over time to overcome instabilities. Nevertheless, we expect the instabilities caused by triadic relations to impact performance adversely at least temporarily, as discussed in Sec. I. Hence, the empirical investigation of the impact of cycles on delivery performance and inventory costs is a promising area for further research.

When transferring our results to real supply networks, one has to be aware of the limitations of the underlying model. We neglect time delays due to transport or manufacturing processes, forecasting, and order batching, which are known to contribute to the bullwhip effect.⁵⁸ Therefore, the model captures only topology-induced instabilities and there may be other sources of instability. More detailed models that incorporate these aspects can be developed by incorporating discrete material delays and additional state variables for forecasts and work-in-process and transport inventories. In addition, we use continuous flows in our model and thus for small stock levels, the model will only be valid for large timescales. For very small stock levels and smaller timescales, stochastic models have to be applied, e.g., agent-based models,^69–71 and stochastic fluctuations can be included explicitly. However, note that our results are, due to the definition of stability itself, robust towards small fluctuations in stock-levels.

Since we gained a basic understanding of the topological instabilities that can emerge in a small triadic motif, we can either proceed and study additional important motifs (building blocks) of supply networks (see Ref. 72 for examples) or investigate larger networks. In particular, it is important to understand, how the instabilities caused by triadic motifs impact the supply network. Furthermore, multiple end products and complex bill of materials can be considered. The model can be extended to such complex settings in a straightforward way, where the rate of change of inventory level of any part is expressed as the sum of the rates of processes where it is produced minus the rates of processes where it is consumed. Therefore, structurally similar but larger models can be developed for more realistic networks with large number of firms and complicated bill of materials structures. For instance, Demirel et al.⁴² considered networks of up to $106$ firms. However, with increasing network size, it becomes infeasible to perform a bifurcation analysis and the statistical ensemble approach has to be applied, as discussed in the introduction.

In summary, we used generalized modeling to investigate the dynamical behavior of a small supply network motif. We showed that, while serial supply chains are remarkably stable, the triadic motif is one of the simplest supply networks causing topology-induced instabilities. We have analyzed these instabilities in the network and identified their causes. We showed that the inter-dependency between the production rates of parts manufactured by a single supplier can be highly destabilizing. Therefore, we conclude that it is of high importance for supply network stability to monitor co-producing suppliers closely and to be aware of the network structure, i.e., adjacent tiers have to be aware of the underlying relationships.

This work was supported by the Engineering and Physical Sciences Research Council (EP/K031686/1).

The aim of this study is to analyze the stability of supply networks and changes in the dynamic behavior. To do so, we use basic concepts from bifurcation theory, which we briefly introduce in the following. For a more detailed discussion, we refer the reader to Refs. 40 and 41.

In general, an $N$-dimensional dynamical system can be described by a set of $N$ ordinary differential equations,

where $f=(f1,…,fN)$ represents the change of the state vector $x=(x1,…,xN)$⁠. Since our analysis focuses on the long time behavior of supply networks, we are interested in states that do not change their behavior over time, so called steady states. These fulfill the condition

The stability of a steady state is given by the response of the system to small perturbations. If a perturbation drives the system away from the steady state, the state is considered unstable, while it is considered locally asymptotically stable if it returns to the original state. The local asymptotical stability of a steady state can be mathematically determined by the eigenvalues $λ$ of the Jacobian $J$⁠. The elements of the Jacobian $Jij$ are defined as

The eigenvalues $λ$ are given by the roots of the characteristic polynomial, $ρn(λ)=det(J−λI)=∑k=0nckλk$⁠, with coefficients $ck$⁠. If the real part of all eigenvalues is negative, the steady state is locally asymptotically stable, otherwise it is unstable.

The eigenvalues, and thus the stability, often depend on certain parameters. In many systems, a threshold for parameter values can be found for which a steady state changes its dynamic stability and where, in addition, the number of steady states might change. These thresholds are called bifurcations points⁴⁰ and they occur, if the real parts of a subset of eigenvalues become positive. This can happen in two fundamental ways. Either the real part of a single eigenvalue becomes larger than zero [Fig. 7(c)], or a pair of complex conjugated eigenvalues crosses the imaginary axis. The first case corresponds to saddle-node type bifurcations (including saddle-node, transcritical, and pitchfork bifurcations) at which two steady states collide and either exchange stability (transcritical bifurcation, see Fig. 7(a) or cancel each other out (saddle-node and pitchfork bifurcation). This class of bifurcations occurs if the determinant of $J$ is equal to zero (⁠$det(J)=0$⁠). The second case is referred to as Hopf bifurcation, where a stationary steady state becomes unstable (stable) and a stable (unstable) limit cycle emerges, which can be observed as a contained oscillation [see Fig. 7(b)]. To determine if a Hopf bifurcation occurs, we use the resultant $RN$ of the system, which vanishes if two symmetric eigenvalues cross the imaginary axis, and the Hopf number $χ$⁠, which is negative if the eigenvalues have a non-vanishing imaginary part⁷³. The Hopf number $χ$ is thereby associated to the Hopf frequency, given by the square of the complex part of the eigenvalues. The conditions for a Hopf bifurcation are thus $RN=0$ and $χ<0$⁠. For a vanishing resultant $RN$ and a Hopf number $χ$ equal to zero, a Takens-Bogdanov point occurs (see below). To calculate the resultant $RN$ and the Hopf number $χ$⁠, we use the equations derived by⁷⁴, which result for a system with $N=6$ in

with $ck$ being the coefficient of the characteristic polynomial.

Both saddle-node type and Hopf bifurcations are called codimension-one bifurcations as they only require the modification of a single parameter. Bifurcations that necessitate the change of several parameters are referred to as higher codimension bifurcations, such as Takens-Bogdanov and cusp bifurcations. These occur at intersections of codimension-one bifurcations. Takens-Bogdanov bifurcations occur at the intersection of a Hopf bifurcation and a saddle-node type bifurcation and homoclinic orbits emerge in their vicinity. Cusp bifurcations emerge if two saddle-node type bifurcations collide and imply the presence of hysteresis phenomenon.