Modelling supply networks and business cycles as unstable transport phenomena

Physical concepts developed to describe instabilities in traffic flows can be generalized in a way that allows one to understand the well-known instability of supply chains (the so-called ‘bull-whip effect’). That is, small variations in the consumption rate can cause large variations in the production rate of companies generating the requested product. Interestingly, the resulting oscillations have characteristic frequencies which are considerably lower than the variations in the consumption rate. This suggests that instabilities of supply chains may be the reason for the existence of business cycles. At the same time, we establish some links to queueing theory and between micro- and macroeconomics.


they have transferred a continuous macroscopic traffic model and reinterpreted the single terms and variables.
We believe that this is a very promising approach to understand business cycles, but instead of simply transferring macroscopic traffic models we suggest derivation of the equations for business cycles from first principles, which means deriving the dynamics on the macroscopic level from microscopic interactions in production systems. This would also make some contribution to the goal of understanding macroeconomics based on microeconomics (or, in a wider sense, based on the 'elementary interactions' of individuals; here, production managers).
In order to make some progress in this direction, we will generalize [13] an idea suggested by Daganzo to describe the dynamics of supply chains [9]. Like the work by Armbruster et al [10], his approach is related to traffic models as well, but he focuses on models in discrete space in order to reflect the discreteness of successive production steps. In order to describe the non-linear dynamics of production processes or interrelated economic sectors, we will have to generalize these ideas to complex supply networks. In section 2, we will first discuss 'macroscopic' business cycles in a sectorally structured economy and compare them with stop-and-go traffic. Afterwards, in section 3, we will develop a more fine-grained, 'microscopic' description of the management of dynamically interacting production units and relate it to classical queueing theory (which mainly focuses on stochastic fluctuations of production processes in a stationary state). Later on, in section 4, we will construct a mathematical relation between the microscopic and the macroscopic level of description, while some further research directions and other potential applications are indicated in the outlook of section 6.

Economic production sectors
We will first investigate a simplified economic system with U production sectors B generating certain kinds of products I ∈ {1, . . . , P}. The market for products I as a function of time t will be represented by the stock level ('inventory') N I (t). Let us assume that, in each production cycle, the production sector B generates p I B products of kind I and requires c J B products of kind J ('educts') for production. If Q B (t) is the number of production cycles per unit time, i.e. a measure of the production rate, productivity or 'throughput' of sector B, the quantity of products generated per unit time is p I B Q B (t), while the quantity of educts consumed per unit time is c J B Q B (t) (see figure 1). The temporal change of the quantity of products I in the market is, therefore, given by the conservation equation The production rate Q B (t) will be specified later on in section 3. For the time being, we will assume that Q B (t) agrees with the actual feeding ('arrival') rate λ B discussed in the next subsection: 90.4 Figure 1. Schematic illustration of the virtual and actual product flows for the case of two kinds of products, I and J . The diagram displays the various model variables related to production sector B, which is represented by the oval structure in the centre. Rectangles correspond to markets, arrows to product flows.

The feeding rates
In the absence of capacity constraints, we would have the relation λ B = ρ 0 B µ B for the timedependent feeding rate λ B , where µ B is the processing rate of sector B, i.e. the potential production rate in the case of no inefficiencies. ρ 0 B reflects the desired utilization. A value of ρ 0 B = 0.7 is reasonable, as it guarantees a relatively high production rate at moderate and reliable waiting times for finishing the products (see section 3.1). However, the actual utilization ρ B = λ B /µ B does not necessarily agree with the desired one, ρ 0 B . In order to illustrate this, let us assume that T J B is the average time needed to deliver to production sector B products of kind J , which are required for production (they will be called 'educts' in the following). If N J is the quantity of actually available products of kind J , the resulting maximum quantity of products that can be delivered per unit time is N J /T J B . As c J B educts are required per production cycle, one cannot have a higher production rate than N J /(T J B c J B ) per unit time, which defines the maximum feeding rate λ J B of educt J into production unit B: Moreover, as the products can only be finished if all the required educts are available in a production cycle, the actual feeding rate λ b is given by the minimum of these values and of the desired feeding rate ρ 0 B µ B (see figure 1):

90.5
If c J B = 0, we set λ J B → ∞, so that the corresponding term does not have any impact on the value of the minimum function (4).
Note that production has some analogy to chemical reactions, where one also requires a certain quantity of educts to produce other (chemical) products. For chemical reactions in three-dimensional space, however, the reaction rate is given by multiplication of (a power of) the chemical concentrations of the educts, i.e. the above minimum function is replaced by an algebraic product. This difference is comparable with the probabilistic AND and the fuzzy AND in fuzzy logic, which are applied when several conditions are to be met at the same time. The probabilistic AND is based on a multiplication of the logical values, while the fuzzy AND corresponds to their minimum.

Adaptation of production speeds and transport capacities
One important aspect is the adaptation of production to the time-dependent inventories N J (t). The adaptation of the desired utilization ρ 0 B is delayed (see section 3), and a change of the processing rate µ B , which requires an adaptation of production sectors or capacities, takes a considerable time as well. We will assume with a typical adaptation time τ B . W B (1/Z B , . . .) reflects the production rate of sector B at steady state as a function of the inventories N J . For the time being, we will assume i.e. the production speed is adapted only to the inventories N J of the products J generated by sector B, for which p J B > 0. X J B is an additional and constant prefactor. This formula is defined in a way which always gives positive values Z B (t), if not all coefficients p J B are zero and the N J (t) stay positive (see below). Normally, the production speed W B (1/Z B , . . .) will increase with decreasing inventories N J , but it saturates due to financial or technological limitations and inefficiencies. In this study, we will assume where Z B could be called 'scaled inventories' and A B , B B and D B are suitable parameters (see figure 2). Furthermore, we will assume here that the transport capacities are adapted in parallel with the production speed and represent the proportionality constants by V J B . This implies corresponding to the specification 1/T J B = ρ 0 B µ B V J B , which results in the following formula for the feeding rates: Herein, the minimum extends over all indices J .

Bull-whip effect and stop-and-go traffic
For a linear supply chain with c I B = δ B,I +1 and p I B = δ B,I (where δ K ,L = 1 if K = L, and δ K ,L = 0 otherwise) we obtain the particular set of equations In the capacity-constrained case (characterized by small values of V J B ), this leads to with and d(ρ 0 Interestingly, equations (11) and (14) basically agree with the macroscopic traffic flow model by Hilliges and Weidlich [14], where (11) is analogous to the equation for the vehicle density and (14) corresponds to the equation for the average vehicle speed V I in street segment I . These equations are linearly unstable with respect to perturbations of N I (t), if τ I exceeds a certain threshold (see figure 3), which depends on the maximum slope W (X I I /N I ) of W (X I I /N I ) [14]. As drivers (over-)react with a time delay to a changing traffic situation in front, stop-and-go traffic can emerge. The frequently observed instability of supply chains, called the 'bull-whip effect', occurs for similar reasons (e.g. in the 'beer distribution game' [15,16]) (see figure 5).  figure 5, but the perturbation amplitude has been chosen five times smaller. The difference between both curves is the amplitude of the bull-whip effect, i.e. the maximum time-dependent variation in the inventories. Note that there is a critical adaptation time, below which perturbations are not amplified. In this case, the investigated linear supply chain is stable.
The mechanism behind this instability is the delay τ I in the adaptation of the production speeds and transport capacities, which implies over-or underproduction. The repeatedly or periodically resulting high inventories are due to temporary bottlenecks in the supply chain and could be avoided by appropriate control functions [17,18]. Note that, instead of unstable production with relatively high inventories and low production speeds, one may reach the same average, but stable throughput at low inventories and high production speeds (see figure 4). However, this is normally related to higher energy and maintenance costs, so that production tends to operate in the linearly unstable regime.
If the transport capacities are not a limiting factor (i.e. the parameters V J B are large), instead of (11) we have the set of equations Together with (14), this basically corresponds to a particular microscopic traffic model, the socalled optimal velocity model [19], if we restrict our comparison to the linear regime around the steady state and identify ρ 0 I µ I with the actual velocity of vehicle I , but X I I /N I with the inverse distance to the next car ahead (apart from proportionality constants). It is known [17,19] that this model is linearly unstable if i.e. if the adaptation time τ I exceeds a certain threshold which depends on the slope W I (X I I /N I ) of W I (X I I /N I ).  For large enough values of D B , the curve has a maximum at finite inventories. In this case, high steady-state flows Q B appear twice: for low (and linearly stable) inventories and higher (but potentially unstable) inventories (see also figure 6).

Business cycles
The most interesting point is the reaction of the system to a perturbation in the throughputs Q B (t), e.g. a periodic perturbation of the so-called consumption rate Q U +1 (t) (see sections 4.3 and 4.4). The resulting oscillations in the inventories can be much slower and are usually synchronized among different production sectors (see figures 5, 11 and 12(a)). Therefore, they may explain business cycles as a self-organized phenomenon with slow dynamics on the time scale of several years.
The reduction in the resulting oscillation frequency compared with the perturbation frequency is also known from stop-and-go traffic. It has been explained by the non-linearity in the model equations. However, there is a significant difference between the dynamics of traffic flows and supply chains: while stop-and-go waves have a characteristic amplitude independently of the average vehicle density, the amplitude of oscillations in the inventories changes continuously in the capacity-constrained case (see figure 6). In other words, the phase transition from stable to unstable traffic flow is hysteretic (i.e. of first order) [20,21], while the phase transition from stable to unstable supply chains in the capacity-constrained case appears to be continuous (i.e. of second order). This seems to be a particular property of the Hilliges-Weidlich model, while the optimal-velocity model mentioned in section 2.4 is known to display a hysteretic transition [19]. Therefore, a hysteretic transition is found for supply chains which are not constrained by their transport capacities. The transition between these two different regimes would be interesting to investigate.
Note that, in order to find the emergence of slow oscillations, i.e. of business cycles, we do not need to have a linear supply chain. Supply networks can display similar features (see section 5). The only requirement is that no stationary state exists or the stationary state is linearly unstable with respect to perturbations. This is being intensively studied for particular network types in a recent scientific collaboration [18].   figure 5, but the perturbation amplitude has been chosen five times smaller (of the order of 2%). The difference between both curves reflects the amplitude of the maximum time-dependent variation in the inventories. Note that, in the capacity-constrained Hilliges-Weidlich case investigated here, the change of the amplitude (and the associated kind of non-equilibrium phase transition) is continuous. The supply chain is stable with respect to perturbations, where both curves agree. This is the case at small and large steady-state inventories.

90.10
Summarizing our present insights, preconditions for the emergence of business cycles are • large values of the adaptation times τ B , • significant changes of the production speeds W B with changing inventories, • the presence of perturbations, if a steady system state exists, • no consideration of forecasts or poor forecasts of the future time development of the inventories N I (t), • a highly non-linear interaction among the production sectors.
We should finally emphasize that the above considerations can be applied to economic systems with any number of sectors and any network structure. We do not even need to assume a sectorally structured economy, as the same kinds of equations apply on the 'microscopic' level of production networks, only with a significantly higher number of equations (see section 4.1). The special cases discussed above have been chosen only for illustrative reasons.

'Microscopic' model of production processes
In this section, we will formulate a generalized model for the dynamical interaction of production processes. Compared with the description of interacting economic sectors in the previous section, this approach may be called microscopic. Note that the concept developed in this section is required for two reasons: first in order to be able to describe real production processes, which involves buffers and other variables and calls for a more complex model; second in order to allow the derivation of the macroeconomic dynamics from microeconomic assumptions regarding the production management of single companies (see section 4). Readers not interested in these aspects may skip this section and continue with section 5.

Production units in terms of queueing theoretical quantities
We will now investigate a system with u production units (machines or factories) b ∈{1, 2, . . . , u} producing p products i, j ∈ {1, 2, . . . , p}. The respective production process is characterized by parameters c j b and p i b : in each production cycle, production unit b requires c j b products ('educts') j ∈ {1, . . . , p} and produces p i b products i ∈ {1, . . . , p}. The number of production cycles of production unit b per unit time is a measure of the throughput and will be represented by Q b (t). For production in the steady state, we can often approximate this quantity in terms of variables from queueing theory [22]. For this, let λ b be the feeding ('arrival') rate, C b the number of parallel production channels, µ b the overall processing ('departure') rate (i.e. C b times the processing rate of a single production channel), the utilization and S b the storage capacity of production unit b. Then the following relationships apply (see figure 7): The functions p λ b and p µ b reflect a reduction in the efficiency of the feeding (arrival) and the processing (departure) in production unit b. They depend on ρ b , C b , S b and possibly other Figure 7. Schematic illustration of a production unit b as a queueing system with a limited storage capacity S b and C b parallel production channels. The arrival rate λ b , the departure rate µ b and the effects of inefficiencies are indicated. quantities as well. In some cases, p λ b is the probability of rejecting arrivals when the storage capacity S b is reached, i.e.
where P b (l) is the probability of having a queue of length l b in production unit b. On the other hand, the average number of active production channels is Consequently, the relative reduction p µ b in the processing rate due to In more complex production systems, the factors p λ b and p µ b would be generalized measures of inefficiencies in the production process, which would be functions of all relevant process parameters.
According to Little's law for queueing systems, is the average arrival rate, we can express the throughput also in terms of the average queue length For example, for a system with infinite storage capacity S b → ∞, we have For an M/M/1:(S b /FIFO) process (one channel with first-in-first-out serving, storage capacity S b , Poisson-distributed arrival times and exponentially distributed service intervals), one finds for Note that both the expected value and the standard deviation of the queue length and the waiting time diverge for ρ b → 1. Therefore, efficient production is normally related to ρ b ≤ 0.7 [22].

The feeding rates of production units
In the absence of capacity constraints, we just have the relation λ b = ρ 0 b µ b for the feeding rate, where the actual utilization ρ b agrees with the desired utilization ρ 0 b , e.g. ρ 0 b = 0.7. However, a lack of required educts may lead to a reduction of λ b . In that case, the feeding rate is limited by the minimum arrival rate A j b of required educts j , divided by the quantity c j b of educts needed for one production cycle. We will assume that the maximum arrival rates A j b are given by is the maximum transport rate due to capacity contraints V j b for getting educt j from the input buffer into the production unit b, and the prefactor ρ 0 b µ b suggests that the transport capacity is adapted to the production speed (generalizations are possible). Therefore, if λ denotes the maximum feeding rate for educt j into production unit b, the actually resulting feeding rate is where the minimum extends over all indices j . The product flows are illustrated in figure 8. When we assume that the flows S j ab of products j to b from the delivering production units a can be fed into the production process in parallel (instead of having to go through the input 90.13 buffer first), we have the generalized relationship with ε = 1. The previously discussed case (requiring delivery through the input buffer) corresponds to ε = 0. In any case, the fraction of the supply a S j ab , which is not needed to satisfy the production requirements c j b Q b , is delivered to the input buffer.

Input and output buffers
Let us assume each production unit b has input buffers for required products i ('educts') and output buffers a ('warehouse') for the products. We will assume the input buffers are filled with I i b (t) educts i ∈ {1, . . . , p} and the output buffers with O i b (t) products waiting to be delivered. If S i ab (t) denotes the delivery flow ('supply') of products i from production unit a to b, the change of an input buffer stock with time is given by the conservation equation as a S i ab (t) is the quantity of products i delivered from various sources (production units) a, and c i b Q b is the quantity of educts i used up for production per unit time (see figure 8). Analogously, the dynamics of an output buffer stock is determined by the equation as p j b Q b is the quantity of newly generated products j , and c S j bc (t) are deliveries to other production units c (see figure 8). The following specifications in this paper ensure the nonnegativity conditions I i b (t) ≥ 0 and O i b (t) ≥ 0 (see also sections 3.5 and 4.1). We will, for example, assume that the delivery flows S i ab are basically given by the order flows ('demands') D i ab , but limited by the maximum transport rates V i ab ρ 0 b µ b for delivering the available products from output buffer O i a to production unit b: If required, suitable specifications can also guarantee ) κ ], which stops the delivery flow of product i when the respective input buffer of production unit b is full. Analogously, we may multiply ) κ ], which stops the production when one of the output buffers of b is full. High values of κ produce a hard cutoff, while small values of κ ≥ 1 can describe cases where the production efficiency goes down even before the buffer size is fully used up.

Adaptation of production speeds and transport capacities
One important aspect is the adaptation of production to changing demand. On the one hand, a change of the processing rate µ b is expensive and time consuming. On the other hand, the adaptation of the desired utilization ρ 0 b is delayed by the average queueing time T b and other factors. We will assume with a typical adaptation time τ b significantly greater than the adaptation times of the other variables. (That is why we do not consider possible time delays in the other mathematical relations, here.) W b (1/Z b , . . .) is some control function reflecting the strategy by the production manager in adapting the utilization and/or processing rate to the output buffer stocks O j a or other variables. In the following, we will assume i.e. the management strategy is only sensitive to the perceived stock levels N j b of the products j generated by unit b, for which p j b > 0. X j b are proportionality constants. The perceived stock levels will be specified by For δ = 0, the management takes into account only its own output buffer stock O j b , while for δ = 1, it also tracks the output buffer stocks O j a of competing production units a. Normally, the control function W b will increase with decreasing perceived stock levels N j b , but the function W b (1/Z b , . . .) saturates due to financial, spatial, or technological limitations and inefficiencies in the processing of high order flows. Again, we will use a function of the form with suitable parameters A b , B b and D b . A b corresponds to the maximum production speed.

Order flows, delivery networks and price dynamics
The flow of orders is basically given by the flow c i b Q b of educts required for the production by unit b. Deviations can be reflected by a correction function denotes the desired input buffer stock. If the proportions q i ab with a q i ab = 1 reflect orders from different producers a, we have and, therefore, (see equation (29)). When the first term is the smaller one, this implies with (27) dI

90.15
The function U b (I i b , I i0 b ) should be chosen greater than one if I i b is smaller than the desired input buffer stock I i0 b , smaller than one for I i b > I i0 b , and otherwise one, e.g.
Without this correction function, the input buffer tends to be emptied in the course of time.
Together with c i b and p i b , the fractions q i ab characterize the delivery or supply network. One can imagine various scenarios. For example, the fractions could be assumed to be constant or modelled by an evolutionary selection equation with a selection rate ν [23]: The specification of the fitness F i ab depends on the relevant parameters. For example, it could be treated as constant. However, in some cases, it makes sense to relate the fitness F i ab to the inverse of the real or virtual costs ('prices') p i ab of product i , when delivered from production unit a to b: Assuming a law of supply and demand, one conceivable specification would be where p 0 i are the costs when the supply S i ab agrees with the demand D i ab . Reduced prices would require a slight generalization of equation (29).
Other specifications are, of course, possible as well, as the above formulae partly depend on the strategies of the human decision makers involved.

Calculation of cycle times
Apart from the productivity or throughput Q b of a production unit, production managers are also highly interested in the cycle time, i.e. the time interval between the beginning of the generation of a product and its completion. Let us first discuss the process cycle time t b between entering the queue of production unit b and leaving it, assuming that all c i b required educts i for one production cycle are transported together and located at the same place in the queue. The problem is similar to determining the travel times of vehicles entering a traffic jam [29].
According to equation (22), the average waiting time is determined as the quotient T b = L b /Q b of the average queue length L b and the throughput Q b , if production operates in the steady state. Let us now generalize this formula to situations in which the inflow into production unit b is time dependent and possibly differs from the time-dependent outflow (see figure 7). In reality, this time dependence results from fluctuations in the production process, breakdowns of machines, etc. In some cases, one can use the length-dependent formulae 90. 16 and (see section 3.1). We will now derive a delay-differential equation for the cycle time under varying production conditions. For this, let l b (t) be the actual length of the queue in production unit b at time t. The change of this length in time is given by the difference between the inflow and the outflow at time t: If the production unit b has C b channels, an educt entering the production queue at time t must move forward l b (t) − C b steps before it is finally served by one of the C b channels. If, after entering the queue at time t, t * b (t) denotes the waiting time until one of the channels is reached, and if the average processing rate of a single channel is µ b /C b , the average treatment time is given by The overall time t b required for the processing of the product is, therefore, given by the sum of the waiting time t * b and the treatment time by one of the channels: (which replaces the average value T b ). On the other hand, the waiting educts move forward Q out b steps per unit time. For this reason, the waiting time t * b (t) is given by the implicit equation Identifying the time derivative of this equation with equation (45) results in which leads to the delay-differential equation As the production initially starts with a waiting time of t * b (0) = 0 (when the factory or production unit b is opened), this equation can be solved numerically as a function of the outflow Q out b (t ). In this way, it is possible to determine the waiting time t * b (t) and process cycle time t b . In the future, approximation methods will be developed for cases where t * b (0) is not known or Q out b (t ) cannot be controlled or anticipated. A rough approximation would be to replace these values by mean values (cf equation (22)).
In a similar way, one can calculate the waiting time t i b in the input buffer The waiting time t j b in the output buffer O j b is given by dt

90.17
Finally, the transport time t i ab between the output buffer of production unit a and the input buffer of production unit b can be estimated by The total cycle time is, then, calculated as the sum of the respective waiting, treatment and transport times.

Relation between production processes and macroeconomics
We will now give some support for the 'macroscopic' model of supply networks applied in section 2 by relating it to the 'microscopic' dynamical production model developed in section 3.

Just-in-time production
Let us start by deriving a simplified model of production processes. For this, we will summarize the input and output buffers, defining the stock levels ('inventories') of markets for the products i by According to the balance equations (27) and (28), we find Moreover, typical for just-in-time production, we will assume negligible input buffers Consequently, for ε = 1 and V i ab = V i b , the feeding (arrival) rates become Furthermore, for δ = 1 we obtain N j b = a O j a = N j and the related adaptation equations We will now assume Q b = λ b , which is usually given for small enough channel utilizations ρ b or sufficient storage capacities S b . In that case, one would have to solve equation (57) together Note that formula (58) maintains non-negative inventories N i (t), as required. To show this, let us assume that t would be the first point in time where some inventory N j vanishes, say N i (t) = 0. The inventory N i (t) could only become negative for dN i /dt < 0, which would require c i b > 0. But then, the minimum in formula (58) would be V i b N i (t)/c i b = 0, which contradicts our assumption and proves non-negativity (at least for V j b > 0).

Micro-macro link
In the following, we will try to reduce the number of equations by an aggregation procedure, which keeps the structure of the above model equations. Let us first define production sectors B by summarizing those production units b which are characterized by the same adaptation times τ b and a proportional throughput Q b (t). With suitable constants τ B and k b , we can then assume The total throughput, processing and feeding rates of production sector B are defined by where b ∈ B indicates that b belongs to production sector B. The equations for the inventories keep their structure if we define and use the relation b = B b∈B . Moreover, in order to get the equations and we need to set for all j (requiring that the resulting values on the right-hand sides depend only on B, but not on b ∈ B), and We may also summarize all those markets i which show a proportional dynamics of N i (t) in time. This assumes if we denote the proportionality factors by X i and define

90.19
where i ∈ I indicates that i is part of the market sector I . Introducing and summing equation (61) over i ∈ I yields the further reduced set of equations In order to obtain the equations and we have to define and Based on assumptions (59), (65) and (67), one could achieve a considerable reduction in the number of equations, leading from a microscopic model of production processes involving management decisions to a macroscopic model of interacting economic sectors. Conditions (59), (65) and (67) presuppose similar production processes of the summarized production units b and proportional coefficients regarding the summarized markets i . While our aggregation method can be generalized to time-dependent parameters k b (t) and X i (t), we still require τ b = τ B , and the resulting values of V should depend only on B, but not b ∈ B. These are the main criteria for defining production and market sectors based on empirical data.

Dynamic input-output model
In equation (70), the sum over B extends from 1 to U + 1, with B = U + 1 representing the final consumer sector. Splitting it up, the resulting equations read as p I U +1 = 0. Similarly, J in formula (71) runs from 0 to P with J = 0 representing a market sector for basic resources. Therefore,  The circle summarizes the production sector (oval) and its market (rectangle). The product flow into its market is represented by a thick arrow, the product flows from markets to production sectors by thin arrows. Bottom: the particular case of a linear supply chain, in which each production unit receives goods only from the previous one.
Let us now concentrate again on the case Q B (t) = λ B (t) and define the consumption rate from market sector I by When we define the U production sectors B through the respective kind of products they produce, we can set p I B = k B δ B,I (see figure 9). Without loss of generality, we can choose k B = 1, which just defines the unit quantity in which we measure products of market sector I . In the stationary case with dN I /dt = 0 for all I , we then obtain the equation which corresponds to Leontief's input-output model of macroeconomics [24]. Therefore, the above model of supply networks can be considered as a dynamical generalization of this classical economic model. We should note that, in view of the instability of supply networks, i.e. the existence of the bull-whip effect and of business cycles, the applicability of steady-state concepts in economics is questionable. People have, therefore, tried to formulate dynamical input-output models for a long time in order to take into account investment strategies and other aspects. However, many of these approaches have turned out to be inconsistent or useless (cf [25]). In contrast, the above dynamical input-output model results naturally from a much more general model of supply and production networks. Investment strategies could, for example, be taken into account by a suitable specification of W B (· · ·), which may be chosen as a function not only of the inventories N I , but also of time derivatives such as dQ I /dt, dY I /dt or dN I /dt.

Specification of the boundary conditions
In the following, we will discuss some examples. If we assume something like a conservation of materials or value, this implies certain constraints, which reduce the potential complexity of the related supply networks. First of all, the quantity of products I consumed by the production sectors B should be generated somewhere, i.e.
Second, the quantity of educts consumed by some production sector B corresponds to the quantity of its generated products, i.e.
For closed systems (with c 0 B = 0 and c I U +1 = 0), this implies that the sum over all inventories N I is constant, i.e.
In the following, we will again discuss the case p I B = δ B,I . In this particular case, relations (80) and (79) As in section 4.3, we will extend the sums over I from 0 to U and the sums over B from 1 to U + 1. In this way, we define Values c 0 B > 0 allow us to describe the inflow of basic resources I = 0, while c I U +1 > 0 allows us to describe the depletion of products by a consumer sector B = U + 1, the ageing of products, or the generation of products of minor quality. The boundary conditions are completely defined by specifying N 0 (t) and Q U +1 (t), see equations (75) and (76). This treatment is fully consistent with the intuitive one of linear supply chains in section 2.4.

Impact of the supply network's topology
We will now discuss simulations based on the equations specified in sections 4.1-4.4 for the three different supply networks sketched in figures 10(a)-(c), each with five levels: (a) a linear supply chain with 5 production units, (b) a 'supply ladder' with 10 production units and (c) a hierarchical supply tree with 31 production units.
By introducing random variables ξ L K , which were assumed to be equally distributed in the interval [−η, η], we have taken into account a heterogeneity η in the individual parameters characterizing the different production units. Here, we have chosen N 0 (t) = N 0 = 20 = X J B , c 0 B and c I U +1 are defined in accordance with equation (83). These specifications guarantee that, for η = 0, i.e. if the production units are characterized by identical parameters, the dynamics of the inventories is the same for all three network topologies discussed (see figure 11). However, the topology matters a great deal if we have a heterogeneity η > 0 in the model parameters (see figures 12(b), (c)). As our dynamical model of supply networks assumes non-linear interactions, changes of η can have large effects. The same applies to small changes in N 0 (see figure 6) or in the relaxation times τ B (see figure 3).

Summary and outlook
In this paper, we have sketched a dynamical theory of supply networks. Here, we could only present first results and indicate possible future research directions, which may contribute to the interdisciplinary field of econophysics [26]. The proposed theory is developed to help understand the complex non-linear phenomena in production and supply networks, in particular their breakdowns, instabilities and inefficiencies. We have also sketched how to derive 'macroscopic' equations for the dynamics of a sectorally structured economics from 'microscopic' equations describing single production processes, involving strategic decisions of production managers reflected by the control functions W b . The resulting model is a dynamical generalization of the classical input-output model of macroeconomics. For the particular case of a linear supply chain,  figure 11 one can conclude that heterogeneity in supply networks can considerably decrease the undesired oscillation amplitudes in the inventories. The strongest effect by far is found for supply ladders, which is relevant for the design of robust supply networks.

90.25
one can relate it to the Hilliges-Weidlich model, which was originally developed for traffic flow. These equations can describe the 'bull-whip effect' due to their linear instability in a certain regime of operation. The underlying mechanism is a slow adaptation of the processing rate to changes in the order flows or stock levels in the market. Interestingly, the resulting oscillations in the inventories of the different products i have a characteristic frequency which can be much lower than the underlying fast variations in the consumption rate. Depending on the network structure, these oscillations are synchronized among different economic sectors and may explain business cycles as a self-organized phenomenon with slow dynamics. These conclusions are not restricted to linear supply chains, but can be generalized to many other supply networks which are linearly unstable with respect to perturbations. In reality, business cycles are, of course, less regular and of smaller amplitudes than in figures 5 or 11. However, the above model allows us to reflect these aspects in a natural way by inclusion of fluctuations, heterogeneity, additional capacity constraints and realistic network structures. Investigations with empirical data and for particular kinds of networks are on the way. It should also be noted that already deterministic models of supply networks can show irregular, non-periodic behaviour such as chaotic dynamics [15,27], which calls for suitable control concepts [28]. Compared with traffic dynamics, economic and production systems have some interesting new features: instead of a continuous space, we have discrete production units or sectors, and the production speed as a function of the inventories is different from the empirical velocity-density relation in traffic flow. Due to the minimum condition (56), production systems may operate in different regimes, and small changes of parameters may have tremendous effects. For example, we may have a transition from small oscillations of relatively high frequency to large oscillations of low frequency. Apart from this, the management strategies can vary to a large extent, and with this the control functions W b (· · ·). With suitable strategies, the oscillations can be mitigated or even suppressed [9,17,18]. Moreover, production systems are frequently supply networks with complex topologies rather than linear supply chains, i.e. they have additional features compared with (more or less) one-dimensional freeway traffic. They are more comparable with street networks of cities [29].
Apart from some equations which were not further applied in this study, most of the proposed model equations were conservation equations, equations given by the product flows, or relations derived from stochastic concepts used in queueing theory. They reflect the transport and interaction of products, so that the physics of driven many-particle systems and of complex systems can make some significant contributions to the new multidisciplinary field of selforganization phenomena in production and supply networks.
Future work will have to address questions such as the relevance of the network structure to the resulting dynamics, possible control strategies, the role of the market and pricing mechanism, etc. This particularly concerns the specification of the function W b , and the equations suggested in the subsection on delivery networks and price dynamics. In this paper, we have already seen that the supply network's topology and the level of heterogeneity in production systems have a significant impact on the resulting dynamics: parameter changes can have tremendous effects. For example, a heterogeneity in the parameters characterizing the different production units can stabilize the production and the market considerably, while in traffic flow, heterogeneity normally has a destabilization effect. The stabilization of supply networks through heterogeneity probably could explain why the variations in economic systems appear to be less dramatic than in our simulations, but additional inefficiencies and capacity restrictions (corresponding to finite values of S b , I i,max b and O i,max b ) are probably another reason. This and the role of heterogeneity for the micro-macro link will be studied in more detail in the future. Should it turn out that the micro-macro link requires a high degree of homogeneity in the parameters, it would be preferable to simulate economic dynamics based on a microscopic model of production and supply networks in the future. The tendency for globalization certainly increases the degree of homogeneity, but this also tends to generate larger oscillations in the inventories, i.e. more serious over-and underproduction. At least in some markets, there are definite signs of a development in this direction. Forthcoming publications will, therefore, investigate alternative control strategies (including forecasts). First results on how to decrease the instability of supply chains can be found in [9,17,18].

Advantages, extensions and potential applications
As the variables in our dynamical model of supply and production networks are operational and measurable, the model can be tested and calibrated with empirical data. Moreover, it is flexible and easy to generalize. For this reason, it can be adapted to various applications. Our approach can be related to microscopic considerations such as queueing theory or event-driven (Monte Carlo) simulations of production processes, but, since it focuses on the average dynamics, it is numerically much more efficient and, therefore, suitable for online control. Nevertheless, the formulae can be extended by noise terms to reflect stochastic effects. Our system of coupled differential equations would then become a coupled system of stochastic differential equations (Langevin equations), where the noise amplitudes would be determined via relationships from queueing theory.
The dynamical theory of supply chains appears to be a promising field with many research opportunities. It is not just useful for a deeper understanding of the origin and dynamics of business cycles or for the optimization of production processes and supply networks. It could also contribute to the further improvement of existing traffic control strategies in urban street networks or to the development of more robust routing algorithms for internet traffic. Apart from this, the model can be viewed as a dynamic multi-player game, where the individual control functions W b reflect the strategies of the players b in terms of quantities N i , which represent the information feedback available to them. Generalizing this idea, the above model of supply networks may serve as a basis for particular kinds of neural networks. We are now setting up different projects in these directions, and cooperation is very welcome.