Progress in Eigenvalue Elasticity Analysis as a Coherent Loop
Dominance Analysis Tool
Burak Giineralp
University of Illinois at Urbana-Champaign
Department of Natural Resources and Environmental Sciences
S. Goodwin Ave. Urbana IL 61801 USA
Telephone: ++1 217 333 93 49
Fax: ++1 217 244-3219
guneralp @uiuc.edu
Abstract
Formal model analysis tools are essential elements in understanding how structure drives behavior.
Conventional model analysis relies heavily on a time-consuming experimental iterative process. Current
formal tools are not mature enough for application to most models. This paper presents a loop dominance
analysis approach based on eigenvalue elasticity analysis. Eigenvalue elasticity analysis, although a
potentially strong formal model analysis tool, has drawn criticisms over the years for a number of reasons.
The approach proposed in this study attempts to bring proper solutions to the issues raised by those
criticisms. To this end, a ten-step procedure is proposed. Among the most prominent features of the
proposed procedure is the ability to track the influences of feedback loops on a specific variable of
interest. Others include the ability to track the loop dominance dynamics over time and an attempt to the
codification of the proposed features of the eigenvalue elasticity analysis. Several programming codes,
written for this purpose in mind, are presented in the paper. The application of the proposed approach is
demonstrated using a simple economic long wave model and two other models, all chosen from earlier
methodological studies on formal loop dominance analysis. The results of these applications also facilitate
the comparison of the proposed approach to other formal model analysis tools.
Keywords: Loop dominance, nonlinear dynamics, complex systems, eigenvalue elasticity analysis,
formal model analysis, model structure, economic long wave.
1. Background
1.1. Motivation
The concept of feedback in systems thinking essentially stems from the notion that endogenous sources
are responsible for the creation of system behavior. Akin to this view, the premise of the system dynamics
has been to bring understanding to how structure drives behavior. “Understanding model behavior” has
claimed the first rank in a list of eight problem areas put forward as currently deserving the attention of
system dynamics practitioners (Richardson 1996).
Although the importance of structure in unfolding of model behavior is acknowledged, the formal tools to
reveal how structure generates behavior have largely been missing in the system dynamics applications.
The established practice in model analysis in system dynamics field is to formulate relevant hypotheses
and to perform repeated simulations to test them. Furthermore, this experimental approach has to be
coupled with intuition and experience on the part of the analyst in order the approach to be fruitful. The
difficulty of carrying out this process becomes burdensome even for experienced modeler in dealing with
complex models. While this approach is a necessary part of model analysis and probably will never be
completely replaced it should be supported by formal analysis tools in order to increase the efficacy of the
existing practice. Over the years, there have been attempts towards a better understanding of structure-
behavior linkages and the formulation of formal tools for this purpose (Forrester 1982; Davidsen 1991;
Richardson 1995; Ford 1999; Saleh 2002; Mojtahedzadeh et al. 2004; Oliva 2004).
The earliest and mathematically the most rigorous of the formal model analysis tools is the eigenvalue
elasticity analysis (Forrester 1982). This study is an attempt to advance the eigenvalue elasticity analysis
as a methodology for tracking the dynamics of feedback loop dominance in high-order nonlinear models.
The overall goal is to develop a theoretical and methodological framework to uncover structure-behavior
relationships in highly complex nonlinear models. Specific objectives of the research are to determine the
dominant behavior modes in highly nonlinear models and the contribution of models’ feedback loops to
these dominant modes.
The research holds promise of opening ways for possible theoretical developments in how certain
feedback loops operate together to give rise to various behavior modes in nonlinear models (Giineralp
2004). Practical contributions will be identification of parameters that play a significant role in generating
the dynamics of a model. This improvement would lead to identification of efficient policy leverage
points in the particular system under study. Finally, successful implementation of the proposed
methodology on a model would increase the credibility of the model for use in policy analysis and
decision-making practices. Accomplishment of these objectives will contribute to the analysis of complex
nonlinear systems by providing a numerical tool to identify significant causal structures that underlie the
dominant behavior modes in highly nonlinear systems.
This paper presents the results of the first stage of the research on loop dominance dynamics. The later
stages of the research will concentrate on the coupling of the eigenvalue elasticity analysis with the error
analysis tools (Gertner 1987). Important objectives of this type of analysis include the identification and
quantification of different sources of uncertainty that propagate through the model affecting model
behavior patterns. The coupling will enable to identify the effects of parameter uncertainty on the loop
dominance dynamics in a model. The broader impact of the research when completed will be twofold: Its
contribution to the analysis of complex systems by bringing a new understanding to how structure drives
behavior in nonlinear models, and by identifying the effect of parameter uncertainty on these structure-
behavior relations.
1.2. Loop dominance analysis
Analyzing the structure-behavior relationship in a model essentially reveals what parts (which feedback
structures) are responsible for which behavior mode of the model and to what extent. The central theme in
uncovering structure-behavior relations is the “dominant loop” concept. The model behavior evolves as a
consequence of the dynamic interactions between various feedback loops of the model. As a result, the
influence of any feedback loop at any time on the generation of model behavior may change. The most
influential structure in determining some segment of the dynamics of a system is called dominant and the
analysis that aims to uncover these loops is called the loop dominance analysis.
Uncovering dominant feedback structures is important in a number of ways. For instance, they may assist
in formulating principles on how certain feedback loop combinations function in concert to generate
observed behavior patterns. Graham (1977) constitutes the first and only attempt in organizing a set of
such principles. Recently, Giineralp (2004) presented interesting results on the structure-behavior
relationships on a set of linear models that may also have potential implications for nonlinear systems.
Loop dominance analysis has attracted the attention of a handful of system dynamics researchers. These
studies approached the subject from different points of view. On the one hand, there are structure-oriented
tools such as eigenvalue elasticity analysis (EEA) and pathway participation metrics (PPM). On the other,
there is the behavior-oriented approach of Ford (1999).
Forrester (1982) is the first to lay out the foundations of the eigenvalue elasticity analysis in the system
dynamics field. The practice consists of linearizing the model under study at any point in time, calculating
its eigenvalues and then noting how the eigenvalues change as causal link gains change in the linearized
model. The eigenvalues can be thought of as different behavior modes the superposition of which gives
rise to the observed behavior of the system. Thus, EEA, by forming a connection between the model
structure and behavior, provides a means to figure out the dominant structure in the model.
The EEA approach requires identifying all feedback loops and all links they pass through in the model. If
the model is highly complex with tens of feedback loops this procedure becomes prohibitively difficult.
To alleviate the problem, first Kampmann (1996) then Oliva (2004) came up with techniques that allow
for using a relatively small number of loops to determine the most influential structure in the model.
The rigor in the mathematical foundations of the EEA also seems to be one of its disadvantages. The
approach is repeatedly criticized over its reliance on abstract mathematical constructs such as eigenvalues,
which make little sense to the less mathematically advanced. Another weakness of the approach is that it
is conventionally used to identify dominance at the level of the model. In other words, it fails to relate the
identified dominant structure to any selected variable of interest (Ford 1999).
The most recent contribution to the loop dominance analysis has come from Saleh (2002). Saleh, besides
introducing a new measure, called behavior pattern index to quantify the behavior pattern of a state
variable, refined many aspects of the eigenvalue analysis including its application to nonlinear models.
His study once more showed that the EEA is perfectly suited for application to nonlinear models.
A model at its outset may contain a number of redundant substructures. The dominant behavior mode(s)
of the system under study could also be captured without these substructures. Hence, eigenvalue analysis
can also be used to simplify linear models by retaining selected behavior modes (Eberlein 1984).
Mojtahedzadeh (1997) has proposed another mathematical tool that would aid in understanding structure-
behavior linkages. His method makes use of a pathway participation metric (PPM), a measure devised
based on ideas first elaborated in Richardson (1995). In contrast to the EEA, which traditionally
approaches the structure-behavior relations at the level of model, PPM adopts a minimalist approach in
that pathways between two state variables are considered as the primary building blocks of influential
structure. Eventually, some combinations of pathways define the most influential system structure, which
contains a single feedback loop, in determining the behavior of a state variable at a given time interval.
PPM stands as the only approach whose features are implemented in an experimental piece of software,
Digest (Mojtahedzadeh et al. 2004).
One criticism directed towards the PPM method is that it identifies only a single feedback loop at any
instant to explain the behavior of a selected variable (Saleh 2002). However, there are actually many
interacting feedback loops that are influential in generating the behavior of the selected variable at any
time interval.
Another problem about the PPM method is its somewhat myopic approach to structure-behavior
relationship. In other words, by confining itself to a single path of dominance, the method misses what is
going on in the rest of the model and does not capture the model-wide dynamics. It is also possible that as
will be discussed later in the paper, PPM approach, in its current form, fails to properly reveal the causes
of oscillatory behavior.
Approaching the loop dominance question from the behavioral perspective, the distinguishing feature of
Ford’s approach is its exclusive emphasis on behavior patterns rather than the structure or system
conditions in models (Ford 1999). The behavioral approach defines dominance with three basic behavior
modes: linear, exponential and logarithmic. Furthermore, it identifies simultaneous dominance by
bringing in the concepts of multiple dominance and shadow feedback structures. While suffering from the
similar shortcomings such as the lack of automation and rigor to be applied to large scale complex models
the behavioral approach, by adopting a different vantage point and by raising important questions
regarding the purpose and implementation of feedback loop analysis, offers new prospects to enhance the
potential of existing loop dominance analysis tools.
2. The methodology
Eigenvalue elasticity analysis was first proposed by Forrester (1982) as a method for relating the strengths
(gains) of individual feedback loops to the behavior modes of a linear system. Forrester concentrates on
the gains of the causal links that constitute the structure of the model. The causal structure of a linear
model can be represented as a gain matrix (Eq. 1). Each entry in the gain matrix, which is equivalent to
the Jacobian used in linear analysis near an equilibrium point, represents a compact net gain that
represents the slope of the relationship between the net rate of the state variable p and the state variable q,
ie. the change in the net rate of the state variable p in response to a change in the level of the state
variable q, dap [xq .
Each eigenvalue of the gain matrix G represents an elemental behavior mode the system is capable of
generating. These elemental behavior modes (or behavior mode, in short) that may be present in a model
are:
(a) Monotonic convergent behavior mode associated with a real negative eigenvalue;
(b) Monotonic divergent behavior mode associated with a real positive eigenvalue;
(c) Sustained oscillatory behavior mode associated with a complex conjugate pair of eigenvalues with
zero real parts;
(d) Convergent oscillatory behavior mode associated with a complex conjugate pair of eigenvalues with
negative real parts; and
(e) Divergent oscillatory behavior mode associated with a complex conjugate pair of eigenvalues with
positive real parts.
These five behavior modes are shown in Figure | in the order they are listed. For monotonic convergent
behavior mode, the eigenvalue represents its decay rate, the reciprocal of which is the characteristic time
constant t of the behavior mode. Similarly, a negative real part of a complex eigenvalue pair determines
the decay rate of an exponential envelope around the oscillations. If the real part of a complex eigenvalue
pair is positive the oscillations will grow with time. On the other hand, the imaginary part of a complex
eigenvalue pair determines the observed frequency of the oscillations (Franklin et al. 2002).
Wt,
Time Time Time Time Time
(a) (b) (c) (d) (e)
Figure 1. Five behavior modes of Eigenvalue Elasticity Analysi
The superposition of these behavior modes gives rise to the observed behavior patterns of the model
variables. The gains are combinations of the model variables. Therefore, a change in the value of a certain
variable changes values of certain gains and thus modifies the behavior modes of the model (i.e. changes
the eigenvalues of the model). This may result in an alteration of the overall model behavior. Thus, in the
EEA, eigenvalues and eigenvectors characterize the global link between the model structure and behavior.
ox, Cha
Ox, 7 ox,
G=| , (1)
ox, dx,
ox Ox,
The EEA can be used for formal model analysis even when the system under study is highly nonlinear
(Kampmann 1996, Saleh 2002). The gain matrix of a nonlinear model can be imagined as one the entries
of which change over the simulation time due to nonlinear interactions between its elements. However, at
sufficiently small time intervals such as the simulation time step, the entries of the gain matrix of a
nonlinear model can be assumed to be constants. In other words, the dynamics of the original model is
approximated through time by the behaviors of a series of linear models with varying entries in their gain
matrices. Then the EEA can be applied to these series of gain matrices over simulation time.
Eigenvalue elasticity is a convenient measure of transient-response sensitivity of the model to parameter
changes. The elasticity values are dimensionless, so they are comparable with each other. The elasticity of
a complex conjugate eigenvalue pair is also a complex conjugate pair. In such a case, the real part of the
elasticity gives the effect on the exponential envelope around oscillations, while the imaginary part gives
the effect on the empirical frequency of oscillations (Saleh 2002). The magnitude of the elasticity gives
the overall sensitivity of the cyclic mode. In particular, causal links with large elasticities are important. If
a small number of links have distinctly larger magnitudes than others, this indicates that they define a
dominant subset of the model structure. Most of the time, these distinct links happen to form feedback
loops in the model. A detailed explanation on the computation of eigenvalue elasticities is provided in
Appendix A.1.
A complete description of link elasticities allows one to calculate loop elasticities. This idea is based on a
property of elasticity that the sum of all link elasticities arriving at a node, that is a variable in the model,
equals the sum of all link elasticities departing that node. This allows one to set up a system of
simultaneous linear equations the solution of which identifies the loop elasticities of the model (Forrester
1983, Kampmann 1996). In general, the magnitude of elasticity measures the overall importance of a loop
to a behavior mode. The magnitude can be used to rank loops according to their dominance over each
mode.
The implementation of this procedure to calculate loop elasticities becomes unreasonably difficult if the
model consists of many feedback loops. Kampmann (1996) discovered that a small subset of loops is
sufficient to explain the most influential substructure of a system dynamics model. However, the
“independent loop set” obtained is not unique in most instances and making matters worse, the resulting
structure-behavior explanations depend on the particular independent loop set. Recently, Oliva (2004)
came up with an improvement over the “independent loop set” concept where he, drawing upon the graph
and network theory concepts, proposes an algorithm to put feedback loops of a model in a hierarchical
order. The algorithm produces a “shortest independent loop set (SILS)”. The number of loops identified
this way and their length is unique. Although the SILS is not unique the SILS algorithm does generate a
unique outcome from a given model structure. The SILS algorithm resolves the problem of arbitrarily
forming an independent loop set. It also eliminates the need for using somewhat complicated approaches
such as Mason’s rule and solving the characteristic equation to directly calculate loop elasticities as
proposed by Gongalves ef al. (2000). Although the system of equations will still be overdetermined for a
typical-sized system dynamics model, it will be consistent and hence will yield a unique solution in terms
of loop elasticities (Kampmann 1996).
Oliva and Mojtahedzadeh (2004) has shown on a simple economic long wave model that the most
influential loops determined by the PPM method are contained in the SILS of the model. Therefore, it is
possible to explain the core dynamics based on the loops selected using the SILS method. However, as
they mentioned in their paper, coupling SILS with the EEA would allow for an exploration of interacting
feedback loops based on the simplest, most granular and intuitive loops. Therefore, the resulting analysis
will be system-wide rather than local and based on an intuitive set of loops. This paper presents, among
other examples, the results of and a discussion on a joint implementation of SILS and the EEA on a
simple long wave model.
As with every other loop dominance analysis tool, the EEA is still in its formative phase and because of
this, criticisms over its implementation have not been missing in the literature (Ford 1999, Mojtahedzadeh
et al. 2004). The problems these criticisms raise and how they are tackled in this study are explained in
the following paragraphs. It is hoped that the proposed approach will advance the EEA as a methodology
for tracking the dynamics of feedback loop dominance in highly nonlinear models.
First, a typical application of the EEA amounts to a loop dominance analysis at the level of the model
under study. In other words, while behavior modes and the most influential parts of the model are
identified one could not be able to explicitly relate these findings to any state variable in the model. Saleh
showed that it is possible to break the relative changes in net rate of any state variable at any time down to
its constituents each of which is due to each behavior mode of the system. In this paper, it is shown that
by further elaboration, the contribution of each behavior mode on any selected state variable can be
quantified at any time over the simulation run. The following is a brief explanation of the approach; the
details are presented in Appendix A.
During simulation, the dominance of each behavior mode may change over time as a result of the
dynamics among various feedback loops of the system. These changes in turn are reflected in the overall
behavior pattern of the system. As mentioned above, the model behavior is a linear combination of all
behavior modes represented by the eigenvalues of the system (Eq. 2). Thus the contributions of behavior
modes on the overall behavior at any instant can be examined independently as follows.
S=Q0,+..4G¥,+.4+Q,¥, (2)
where s is the slope vector, r; is the right eigenvector associated with the i" behavior mode, and @, is the
coefficient of the i" behavior mode, given by
a, =ane*) i=lLon (3)
where f, is the initial time, a are the initial values of @; at 1, and A, is the eigenvalue representing the
i” behavior mode.
sy = aren, = arr, @)
sit = A ety (5)
At each time step interval, Eq. 4(5) gives the portion of initial (final) slope of the state variable of interest
p coming from the i" behavior mode. The change in the state variable p due to the i” behavior mode is
then
As, =s,° —S. (6)
The resulting changes in overall behavior reveal the contribution of each behavior mode on the state
variable of interest. The comparison between the contributions of behavior modes becomes easier when
the changes they induce are normalized by the sum of the absolute values of those changes. Then the ratio
of the change due to that behavior mode and the sum of absolute changes gives the relative contribution
of each behavior mode to the overall behavior of the variable of interest (Eq. 7). The rescaled changes
vary between —1 and 1. A negative rescaled change at any time step means that behavior mode decreases
the slope of the state variable of interest in that interval. A positive change, however, means that the
behavior mode increases the slope of the state variable of interest in that interval. The closer the rescaled
change is to —1 or 1, the greater the change in the slope due to that behavior mode. It is worth mentioning
that the rescaled contribution values tell nothing about the nature of the behavior mode it belongs or that
of the overall behavior pattern.
As, .
=? _ i=lon (7)
x [Asp
m=l,p
Second, there is the question of how dominant behavior modes play out throughout the simulation. It is
possible that not every state variable of a model exhibits the same behavior mode. As mentioned in the
previous paragraph, the specification of a variable of interest also addresses the issue of different
variables simultaneously exhibiting different behavior patterns.
Third, it is also possible that the behavior mode exhibited by any state variable be composed only of a
single dominant behavior mode at all times during the simulation. There may be brief periods where the
dominant behavior mode of a state variable changes from one behavior mode to another or even periods
of considerable length during which more than one behavior mode dictates the behavior mode of one or
more state variable. These situations certainly have implications on determining the dominant feedback
structure. The conventional approach in EEA is to pick up a single behavior mode as the dominant one,
look at which loop(s) it is most sensitive to and then claim them (those) as the dominant substructure.
This approach, however, simply fails when there is more than a single dominant behavior mode to explain
the behavior of the selected variable. Nevertheless, the quantification of the contribution of each behavior
mode to the behavior exhibited by the selected state variable offers a way out. The rescaled contribution
values of behavior mode may serve as weights for those behavior modes (Eq. 7). This way one can say,
for example, that 60 percent of change in the state variable of interest is due to behavior mode a while the
remaining 40 percent is due to behavior mode b at a particular time step. These values can be used as
weights for the loop elasticities of each behavior mode for each loop in the “shortest independent loop
set”. The resulting values would then reflect the overall influence of loops over the behavior of the state
variable of interest at any time step:
0e, =e k=1.K (8)
isl
where é; is the elasticity of i” behavior mode to loop k.
Since what we really care about is the relative importance of the feedback loops in creating the system
dynamics it would be better if the overall loop elasticities obtained through Eq. 8 are normalized by the
sum of their absolute values (Eq. 9). The rescaled overall loop elasticities vary between —1 and 1. A
negative loop elasticity at any time step means that that loop drives the behavior of the state variable of
interest overall in that interval in negative direction. A positive elasticity, however, means that the loop
drives the behavior overall in positive direction. The closer the overall loop elasticity is to —1 or 1, the
greater the influence of the corresponding loop.
=— k=1.K (9)
Fourth is the lack of automation, which remains as a stringent obstacle to widespread implementation not
only of EEA but of other loop dominance analysis methods as well (Ford 1999, Saleh 2002). While a
necessity, the automation issue is neither something to be resolved in any one study nor in a short period
of time. The new tool that is presented in this study for the analysis of the structure-behavior relationship
should therefore be seen as a humble attempt towards the full automation of EEA features. Nevertheless,
it serves as a starting point for more advanced implementations of the methodology. For this purpose,
several pieces of computer code have been written. A piece of code was written in C language whose
main function is to provide the communication between the modeling/simulation software Vensim® and
MATLAB? where the actual loop dominance analysis is carried out. Vensim® is selected because of its
advanced data handling capabilities. Its Monte Carlo sensitivity analysis module will also be used in the
later stages of this research. MATLAB®, on the other hand, is a tool for algorithm development, data
visualization, data analysis, and numerical computation with matrices and vectors. The pseudo-codes for
the codified portions of the analysis are provided along with brief descriptions of their functions in
Appendix B.
3. The feedback loop dominance analysis procedure
The proposed methodology is formulated as a ten-step procedure. Apart from reflecting where the EEA
currently stands in loop dominance analysis this sequential approach is a step forward towards the
establishment of the EEA as a formal loop dominance method. While tools to carry out some of these
steps are available in some of the commercial SD software others still require manual labor (e.g.
determining the pathways between state variables and the gains of those pathways).
The first four steps stand as exploratory structural analysis because they amount to identifying structural
characteristics of the system dynamics model. The remaining six steps perform on the structural layout
produced by the previous steps and are computationally intensive.
1. Identify and list all nodes (i.e. all elements except the constants and table functions) including the
state variables (i.e. stocks) in the model. At this stage, the state variable of interest for which the
feedback loop dominance dynamics is to be analyzed can be selected.
2. Identify and list all causal links (except the ones that involve constants and flow-to-stock links) and
pathways in the model. The causal links are the structural building blocks of feedback loops and a
link elasticity is the sum of elasticities of all loops which that link is a part of. They need to be
identified so that their elasticities will be computed and used in obtaining loop elasticities later on.
Pathways are causal structures linking one state variable to itself or another. There may be multiple
pathways on which a state variable interacts with itself or other state variables in a model. A pathway
10.
may consist of one or more causal links and a causal link may lie on several pathways. Identification
of these pathways allows for relating eigenvalue elasticities computed from the compact gain matrix
G to the elasticities of causal links.
Identify and list the feedback loops in the Shortest Independent Loop Set (SILS) for the model using
the structural analysis procedure proposed by Oliva (2004). This procedure is an improvement over
Kampmann’s independent loop set concept in the sense that it greatly reduces the subjectivity in the
selection of loops that form the most granular yet comprehensive representation of the model
structure. Interested readers are referred to Oliva (2004) for the details of the SILS procedure.
Form the gain matrix, which represents the links between state variables in their most compact form.
The pathways linking the state variable pairs can be aggregated into single links. The gains of these
composite links are the elements of the gain matrix. In linear models, the equations for these
composite gains are fairly simple and the values they attain are constant; in nonlinear models,
however, they may get quite complicated and their values change over time. This amounts to writing
down the model equations and taking the partial derivatives with respect to each state variable.
Simulate the model and read the gain matrix, pathway gains, net flows of state variables over time
from Vensim into MATLAB. Vensim has built-in functions that allow for exporting simulation
output out to be used by other software. A C code was written which reads data necessary for the
analysis from Vensim and transfer to MATLAB where most of the rest of the analysis is performed.
Determine the characteristics of all elemental behavior modes of the system and the contribution of
each behavior mode to the behavior of the state variable of interest. The identification of behavior
modes of the system amounts to the computation of the eigenvalues and eigenvectors of the model at
each analysis time step. The contribution of each behavior mode at any time step can be determined
by sequentially deactivating all but one of the behavior modes and noting the change in the slope
vector of the selected variable between consecutive time steps. After a suitable rescaling so that they
take on values between —1 and 1, they serve as the weights of the behavior modes that constitute the
overall behavior mode of the selected variable for the concerned time interval. As with the rest of the
analysis this is repeated at every time step. See Appendix A.1 for a detailed discussion and Appendix
B.2 for the pseudo-code of an algorithm that does the relevant computations for third-order models.
Compute the elasticities of behavior modes (eigenvalues) with respect to the compact links in the gain
matrix. For a detailed explanation and the corresponding equation, Eq. A.9, see Appendix A.1. The
pseudo-code of the algorithm that does these computations is given in Appendix B.3.
Compute the eigenvalue elasticities with respect to each causal link selected in Step 2. For a detailed
explanation and the corresponding equation, Eq. A.10, see Appendix A.1. The pseudo-code of the
relevant algorithm is given in Appendix B.4. The algorithm is modified for the analysis of a simple
long wave model.
Using the links identified in Step 2 and the loops in SILS, form the directed cycle matrix. Directed
cycle matrix is a notion borrowed from the graph theory in mathematics and is used to present the
information on which cycles (loops) are formed by which links in matrix form. The rows in a directed
cycle matrix are links, and the columns are loops where the element dj of the matrix is | if link 7 is an
element of loop j and 0 otherwise. The loop elasticities are related to the link elasticities through the
directed cycle matrix (Appendix A.2). The matrix will most surely be overdetermined for almost any
system dynamics model. Nevertheless, its rank will be equal to the number of columns, thus will
always have a unique solution (Kampmann 1996).
Compute and plot the overall loop elasticity values over time and evaluate the findings. Having
determined the loop elasticities in Step 9, the overall loop elasticity of any loop can be computed by
multiplying that loop’s elasticity for each eigenvalue by the corresponding eigenvalue contribution
weight (see Step 6) and rescaling between —1 and 1. The resulting loop elasticity values are plotted
over time, which enables visualizing the evolution of loop dominance dynamics acting upon the
selected model variable over time. The pseudo-code of the algorithm that does the computations of
this step for the simple long wave example is given in Appendix B.4.
4. Application of the proposed methodology
The application of the proposed methodology is illustrated using three models selected based on the
criteria that they should come from earlier methodological loop dominance studies and would give the
most insight (which also happened to come from the most recent studies) on the potential of the EEA
approach. A system story will be weaved around the findings of each model’s analysis. This will result in
coherent, dynamically correct explanations of how certain parts of the system structure create observed
patterns of the system behavior as advocated by Mojtahedzadeh er al. 2004. In addition, comparisons with
earlier approaches will be presented. The model variable names are in italics in the text. The equations of
the models used as examples are given in Appendix C.
Yeast Population model
This is a simple second-order nonlinear model of yeast population dynamics. It is previously used in the
loop dominance analysis context by Saleh (2002). The stock-flow diagram is presented in Figure 2. The
model provides a classical example of overshoot-then-collapse dynamics.
divisiontime lifetime
births deaths
effalc birth effale death
lalcoholgeneration
alcoholpere
ellgeneration
Akohol
Figure 2. Stock-flow diagram of Yeast model.
1. Identify and list all nodes and choose the state variable of interest.
There are eight nodes in the model: Cells, Alcohol, births, deaths, alcoholgeneration,
alcoholpercellgeneration, eff alc birth, and eff alc death. Cells is chosen as the state variable of interest.
2. Identify and list all causal links (except the ones that involve constants and flow-to-stock links) and
pathways in the model.
There are seven such causal links in the model and they are listed in Table 1.
10
Table 1. Causal links of Yeast model relevant to the analysis.
Link no. Variable sequence
1 Alcohol — eff alc birth
2 Alcohol — eff alc death
3 Cells — alcoholgeneration
4 Cells — births
5 Cells — deaths
6 eff alc birth — births
7 eff alc death — deaths
There are two pathways from Cells to itself, one pathway from Cells to Alcohol and two from Alcohol to
Cells (Table 2).
Table 2. Pathways in Yeast model.
Pathway no.
Variable sequence
Oye
ak
Cells, births, Cells
Cells, deaths, Cells
Cells, alcoholgeneration, Alcohol
Alcohol, eff alc birth, births, Cells
Alcohol, eff alc death, deaths, Cells
3. Identify and list the feedback loops in the Shortest Independent Loop Set (SILS).
The loops in the SILS of Yeast model are listed in Table 3.
Table 3. Feedback loops in the Shortest Independent Loop Set of Yeast model.
Loop no. Variable sequence
LI Cells, births
L2 Cells, deaths
L3 Alcohol, eff alc birth, births, Cells, alcoholgeneration
L4 Alcohol, eff alc death, deaths, Cells, alcoholgeneration
4. Form the gain matrix, which represents the links between state variables in their most compact form.
The model is a second-order nonlinear one. Therefore, the dimensions of its gain matrix are 2X2 and the
elements of the matrix change continuously over time. The matrix and the partial derivatives are in Eq. 10
and Eq. 11, respectively. The first and second terms in Eq. 11.a are the gains of the Pathways | and 2
respectively. Similarly, the first and second terms in Eq. 11.b are the gains of the Pathways 4 and 5
respectively. The partial derivative in Eq. 11.c equals the gain of Pathway 3.
11
acey/” acaity/,
G _ dCells dAlcohol
yeast = . . (10)
rc) Aleonol/. tc) Aleohol//
Cells dAlcohol.
a calty/ eff ale birth eff ale death/
Cells ( IY visiontime) *( fiferime) (11a)
a cay, =-0.1+(Cells -eff alc death*Cells
dAlcohol OL ( vcs) +( We cetime) (1b)
a Aleonol/° 7 i :
Cells = alcoholpercellgeneration (11.c)
a Aleonol/, =
@Alcohol =? (ld)
5. Simulate the model and read the gain matrix, pathway gains, net flows of state variables over time.
The behavior of the variable of interest over 90 minutes is shown in Figure 3. The yeast cells grow at a
fast rate during the initial phases of the simulation but then as more and more alcohol accumulates in the
environment eventually they are wiped out after an even faster decline; thus the overshoot-then-collapse
dynamics. The data required to further the analysis is read from Vensim and transferred into MATLAB
using the main function whose pseudo-code is in Appendix B.1.
Cells
40
30
20
LI L113 L2u L2
10
0
O 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90
Time (Minute)
Cells : yeast
Figure 3. Reference behavior of state variable of interest, Cells.
12
6. Determine the characteristics of all elemental behavior modes of the system and the contribution of
each behavior mode to the behavior of the state variable of interest.
The behavior modes that are present at any time step and their relative contributions (weights) to the
overall behavior of the variable of interest are provided on Figure 4. The evolution of these behavior
modes over time is shown in Figure C.1 of Appendix C.1.
The Contributions of Elemental Behavior Modes to the Behavior of Cells
Two
Divergent
Two exponential
‘onvergenfexponential
divergent
‘oscillatory
oscillatory
convergent
behavior mode —_foeh, modefpeh. modes
behavior modes
30%, 60. 70
Cevi
Figure 4. Rescaled contributions of behavior modes on behavior of state variable of interest, Cells.
7. Compute the eigenvalue elasticities with respect to compact links in the gain matrix.
The elasticity of each behavior mode is in matrix form as shown in Figure 5. For brevity, the graphs for
the elasticities are not shown.
evl evi v2 ev? eve eve
evl a1 ein |e en eve a en
yeast | evi ev | jev2 _,ev2 yeast eve eve
2&2 ey ean en ean
(a) (b) (c)
Figure 5. Elasticity matrices for two real eigenvalues (a and b) and complex pair (c).
8. Compute the eigenvalue elasticities with respect to each causal link.
The set of equations that gives the elasticities with respect to each causal link are listed below. For
brevity, the graphs for the link elasticities are not shown.
812, pathway
BAA, parhiway?
@,* 2 5 @s2
" [ "| a
13
x] 812, ag
82
en (5 vA }
u
where g,, is the gain of compact link from state variable ¢ to p, 8pqpanway iS the gain of pathway i from
state variable q to p, and e,, is the elasticity to compact link from state variable g to p. The superscripts
that identify eigenvalues are omitted for clarity.
9. Using the links identified in Step 3 and the loops in SILS form the directed cycle matrix.
The directed cycle matrix for this model is given in Figure 6. It is a 7X4 matrix and its rank equals 4, the
number of loops on the SILS.
[LI 12 13 LA
ci | 0 0 I 0
cl2] 0 0 0 1
cl3 |] 0 0 1 1
cl4.| 1 0 0 0
cl5 | 0 1 0 0
cl6 | 0 0 1 0
cl7| 0 0 0 1
Figure 6. Directed cycle matrix of Yeast model.
10. Compute and plot the overall loop elasticity values over time and evaluate the findings.
The overall loop elasticities are computed and resulting loop dominance dynamics are shown in Figure 7.
The Loop Dominance Dynamics
os
n
08
F Pre
ig 1 LT 3
Soa eee
302
g
0 Let
ieee ue 40 50
02
04
Sasa
Figure 7. Evolution of loop dominance dynamics for Yeast model over time.
The dynamics that give rise to the well-known overshoot-then-collapse behavior in this particular case is
beautifully revealed by the loop dominance analysis. One important distinction of this analysis from the
one presented in Saleh (2002) is that the analysis in the latter study was done for only four time steps.
Each of these time steps is picked arbitrarily from each period during which different behavior modes are
dominant (Figure 4). In contrast, the analysis presented above is truly continuous in time.
14
Nevertheless, the findings of the two study are in agreement except the third analysis time step at t = 70 of
Saleh (2002)'. First of all, both studies identify four phases of model behavior each characterized by
distinct behavior modes: Exponential growth, diverging oscillatory, converging oscillatory, and
exponential decay behaviors (See Figure 4 and Figure C.1).
As for the loop dominance dynamics the two studies identify the same loops as dominant except the
discrepancy mentioned in the previous paragraph. While loop L/ is dominant during the exponential
growth phase loops L/ and L3 are co-dominant during the divergent oscillation phase. The important
point is that from the rescaled dominance graph it is not possible to differentiate which loop supports the
divergent mode. Close examination of elasticities reveals that loop L/ supports the divergent mode while
loop L3 has a restraining effect. During the second phase, loop L/ is mainly responsible for the expansion
of the envelope (i.e. its elasticity is substantially high only for the real part) while loop L3 generates the
oscillations (i.e. its elasticity is nonzero only for the imaginary part) although not as strong as loop L/.
At minute 65.5, the real part becomes negative, oscillatory mode becomes converging and we pass onto
the third phase (Figure C.1). The corresponding elasticity values of the loops change sign accordingly.
These loops are L/ and L2 since loop L3 and loop L4 are only responsible for the generation of the cyclic
pattern. This makes sense because loop L3 and loop L4 are major loops, in other words pathways for the
propagation of disturbances from one state variable to another*. Until this point, loop L3 has lost its
influence giving way to loop L4; loops L/ and L2 are of equal importance though now the elasticity of the
first is negative, that of the latter is positive. At minute 68.87, gain of loop LJ becomes negative and
hence its elasticity becomes positive but by now loop L2 has become the most influential loop on the
convergence of the exponential envelope. The influence of loop L4 is also increasing but it is still smaller
than that of loop L2.
At minute 77.9, the complex conjugate eigenvalue pair bifurcates and two real and negative eigenvalues
emerge, hence the fourth phase: the exponential decay. The elasticities of loops L/ and L2 are positive for
the larger eigenvalue, negative for the smaller. It is the opposite for loops L3 and L4. However, in terms
of magnitude loops L2 and L4 have the two largest elasticities. Towards the end of the simulation loop L2
becomes completely dominant while the resistance of loop L4 gradually dies away.
The application of the proposed methodology in this simple example enabled showing how loop
dynamics evolve over time for a second-order nonlinear model. In addition, this allows tracking the
continuous change in relative influence of loops in SILS. Another improvement is the fact that a variable
of interest, that is Cells, is specified and the loop dominance analysis is done for that variable. This is
important since the EEA is often critiqued for being too general in its conclusions meaning that the
resulting dominant loops could not be associated with a particular variable of interest in the model. Last
but not the least, using the contributions of behavior modes on the overall behavior of the variable of
interest as weights for the loop elasticities for each of these behavior modes is something new and the
preliminary results are promising.
An example with nonlinear oscillations (A Predator-Prey model)
This example is another second-order nonlinear system known as Lotka-Volterra model. It is used by
Mojtahedzadeh (1997) to illustrate the application of PPM approach on a nonlinear oscillatory model;
thus the following analysis not only illustrates the application of the proposed methodology on a model
' The discrepancy is apparently caused by a slight computational error in Saleh (2002).
? This observation is in line with the one stated in Giineralp (2004): In a second-order oscillatory system with one
major and two minor loops, the negative-polarity major loop is only responsible for the imaginary part whereas the
minor loops are mainly responsible for the real part and reduce only slightly the frequency of oscillations.
15
with sustained nonlinear oscillatory dynamics but it also provides an opportunity to compare the two
approaches. The stock-flow diagram is presented in Figure 8.
prey birth
A
prey birth
rate
moO
prey death
\aAt S
prey interaction
L5 constant
Predator. =}
predator birth predator death
La “4 A
predator interaction predator death
constant
rate
Figure 8. Stock-flow diagram of Predator-Prey model.
1. Identify and list all nodes and choose the state variable of interest.
There are six nodes in the model: Predator, predator birth, predator death, Prey, prey birth, and prey
death. Prey is chosen as the state variable of interest.
2. Identify and list all causal links (except the ones that involve constants and flow-to-stock links) and
pathways in the model.
There are six such causal links in the model and they are given in Table 4.
Table 4. Causal links of Predator-Prey model relevant to the analysis.
Link no. Variable sequence
Prey — prey birth
Prey — prey death
Prey — predator birth
Predator — predator birth
Predator — predator death
AunkYNe
Predator — prey death
There are a total of six pathways: Two pathways from Prey to itself, two from Predator to itself, one from
Prey to Predator and one from Predator to Prey (Table5).
16
Table 5. Pathways in Predator-Prey model.
Pathway no. Variable sequence
Prey, prey birth, Prey
Prey, prey death, Prey
Predator, predator birth, Predator
Predator, predator death, Predator
Predator, prey death, Prey
Austwne
Prey, predator birth, Predator
3. Identify and list the feedback loops in the Shortest Independent Loop Set (SILS)*.
The loops in the SILS of Predator-Prey model are listed in Table 6.
Table 6. Feedback loops in the Shortest Independent Loop Set of Predator-Prey model.
Loop no. Variable sequence
LI Prey, prey birth
L2 Prey, prey death
L3 Predator, predator birth
L4 Predator, predator death
L5 Predator, prey death, Prey, predator birth
4. Form the gain matrix, which represents the links between state variables in their most compact form.
The model is a second-order nonlinear one. Therefore, the dimensions of its gain matrix are 2X2 and the
elements of the matrix change continuously over time. The matrix and the partial derivatives are shown
below:
Cc) Prey a Prey
_ oPrey Predator
G predator prey = . . (13)
0 Predator, a Predator/”
oPrey oPredator
tc) Prey = saapier Aedes vaxrit di ,
Yi = prey birth rate + (-prey interaction constant * Predator ) (14a)
* For an alternative stock-flow representation of this model, SILS procedure identified only the first four loops
leaving out the major loop connecting the two stocks. The alternative representation is given in Figure C.2 of
Appendix C.2. Without the major loop the two stocks would be decoupled and there would not be a predator-prey
system. This fact and the following analysis illustrate that it is actually an essential loop for the system to oscillate
not to mention to exist as it is. Although this is a low-order system it gives sufficient reason to further scrutinize the
SILS approach with high-order systems. This example should be regarded as a warning against the faith-based
application of the SILS approach —or of any method for that matter.
17
Co) Prey _ . : ,
Predator = PTY interaction constant * Prey (14b)
a Preiaty/ = predator interaction constant * Predator (14c)
Prey
a Predator dator = Predator interaction constant * Prey — predator death rate (14d)
5. Simulate the model and read the gain matrix, pathway gains, net flows of state variables over time.
The behavior of the variable of interest over 60 months is shown in Figure 9. Because the oscillations are
sustained, the following analysis is conducted over a period of 30 months, which is approximately the
period of oscillations.
Prey
3.25 V\
25
1.75
LW L3 12 | 4 LI, 13 12 “4 LI
1
b+ —
0.25
0 6 12 18 24 30 36 42 48 54 60
Time (Month)
Prey: pp individual/acre
Figure 9. Reference behavior of state variable of interest, Prey.
6. Determine the characteristics of all elemental behavior modes of the system and the contribution of
each behavior mode to the behavior of the state variable of interest.
The only behavior mode for this model is sustained oscillations, which is reflected in the complex
eigenvalue pair of the gain matrix. However, these oscillations are different in nature than the elemental
sustained oscillatory behavior mode in the sense that the first are caused by the nonlinear mechanism in
the system. Although there is only one complex eigenvalue pair the real and imaginary parts continuously
change over time (in fact they oscillate too) because of the nonlinearity in the model. These oscillatory
changes in the real and imaginary parts of the eigenvalue pair are the cause of the observed sustained
oscillations. Since there is only one complex eigenvalue pair of the model at any time there is no point in
calculating the contribution weight, which is obviously one. The figure that shows how the exponential
envelope and the frequency of the oscillations (associated with the real and imaginary parts of the
complex conjugate eigenvalue pair, respectively) evolve over time, is given in Figure C.3 of Appendix
€2.
18
7. Compute the eigenvalue elasticities with respect to compact links in the gain matrix.
The elasticity of each behavior mode is in matrix form as shown in Figure 10. For brevity, the graphs for
the elasticities are not shown.
evteal —_evreal evimag ev imag
Ew | eu e10 Bevis | en ein
predator-prey =| ev real ev real predator-prey =| evimag _,ev imag
e e Coe
21 &2 21
(a) (b)
Figure 10. Elasticity matrices for real (a) and imaginary (b) parts of complex eigenvalue pair.
8. Compute the eigenvalue elasticities with respect to each causal link.
The set of equations that gives the elasticities with respect to each causal link are listed below. For
brevity, the graphs for the link elasticities are not shown.
Ges »(e pathway 1 zi eur = Cy +{ Brinn A ey; =n
§ 22, pathway3 § 22, pathway 812, pathways
C1, =O * sP 3 Gyr = B5% MP « €.=65% ia %
ae [ “%.) eee [ “,] aoe [ “4
9. Using the links identified in Step 3 and the loops in SILS form the directed cycle matrix.
The directed cycle matrix for this model is given in Figure 11. It is a 6X5 matrix and its rank equals 5, the
number of loops on SILS.
[ul 12 2 4 1s
a; ilo o oO Oo
a2) 0 1 0 0 0
a) 0 0 0 0 1
4) 0 0 1 0 0
dw) 0 0 0 1 0
6 | 0 0 OO oO 1
Figure 11. Directed cycle matrix of Predator-Prey model.
10. Compute and plot the overall loop elasticity values over time and evaluate the findings.
The overall loop elasticities are computed and resulting loop dominance dynamics impacting the
exponential envelope and the observed frequency of oscillations are shown in Figure 12. The loop
dominance dynamics impacting the oscillatory behavior as a whole are shown in Figure 13.
19
The Loop Dominance Dynamics for The Loop Dominance Dynamics for
the Exponential Envelope the Frequency of Oscillations
a D i
3
(a) (b)
Figure 12. Evolution of loop dominance dynamics for exponential envelope (a) and for frequency of
oscillations (b) of Prey behavior over time.
The Loop Dominance Dynamics
os Let bs pri Td | ug LI
0.35
loop elatoty
Figure 13. Evolution of loop dominance dynamics for overall behavior of Prey over time.
The results of the analysis show how the nature of the sustained oscillations in this nonlinear system
differs from the nature of those in linear systems. The oscillations are indeed sustained but this is not
because the exponential envelope is non-existent (i.e. real part of the eigenvalue pair is zero) but because
the exponential envelope keeps changing due to the nonlinearity in the model (Figure C.3).4 The
exponential envelope is a diverging one during two segments of the period of oscillations. The first of the
two segments is between months 0 and approximately 6.80 and the second is between months 20.5 and
30. Between these two segments the exponential envelope is contracting (Figure C.3). Figure 12a shows
that both minor birth loops, that is, loops L/ and L3 support the diverging exponential envelope but resist
when it is in the contracting phase. In contrast, both minor death loops, that is, loops L2 and L4 support
the contracting envelope but resist when it is in the divergent phase. In addition, it is worth noting that the
major loop LS plays no role over the exponential envelope.
* The nonlinearity is a simple multiplication of the two state variables reflecting the interaction of the two species
(See Appendix C.2 for the model equations).
20
The frequency of oscillations too changes over time (Figure C.3). The major loop LS is virtually the only
responsible structural piece for the frequency for almost the whole period (Figure 12.b). The elasticities of
the prey birth loop L/ and the predator death loop L4 are positive (potentially leading to higher frequency
of oscillations) for almost half the period and negative (potentially leading to lower frequency) for the
rest. The picture is the opposite for the prey death loop L2 and the predator birth loop L3. The minor death
loop L2 is, however, at least as influential as the major loop around month 20 (Figure 12b).
A word of caution in interpreting loop dominance graphs is in order. In Figure 13, although the two minor
loops, that is the loops L/ and L2, associated with the state variable Prey seems to be the most dominant
pair, actually their influences are opposite to each other as revealed in Figure 12a°. Because of this and
the fact that they act upon the same state variable (i.e. stock), their opposite influences actually cancel out
each other when their magnitudes are approximately equal. This is, of course, valid for the minor loops of
the state variable Predator too. Thus, the actual influence of minor loops on the behavior of variable of
interest is in fact revealed by the sum of the influences of minor loops with a common state variable.
Accordingly, at any analysis time step, if the sum of the influences of the minor loops of Prey(Predator)
is larger then the sum of the influences of the minor loops of Predator(Prey) then the minor loop of
Prey(Predator) with the larger magnitude is the dominant loop at that time step. The time segments
during each of which a different minor loop of the model is dominant are shown on Figures 9 and 13.
Thus the minor prey birth loop L/, the minor predator birth loop L3, the minor prey death loop L2 and the
minor predator death loop L4 are each dominant during the fast growth, peaking, slowing decline, and
grounding phases, respectively. It is worth noting, though, that the shifts in dominance between these
loops are not instantaneous but gradual.
The above explanation is almost totally in agreement with the one based on the PPM analysis
(Mojtahedzadeh 1997). The only disagreement may stem from the interpretation of the role of the major
loop L5. According to the PPM analysis the major loop plays relatively less significant role. The
significance of its role is actually high, but it becomes evident only when the influences of each loop on
the frequency of oscillations (i.e. the imaginary part of the eigenvalue pair) are plotted (Figure 12b). This
fact, however, is not revealed with the PPM approach. It must be emphasized that the oscillations would
not occur at all if there were no such loop that connects the two stocks (Figure 12b). The major loop's
overall influence is not prominent compared to the minor loops (see Figure 13) but here the situation is
certainly not as simple as that by the fact that the minor loops would not function as they do now without
the presence of a major loop.
Another difference between the EEA and the PPM approaches is that the first gives a more complete
picture in the sense that one is able to trace the gradual shifts in loop dominance through simulation and
at the same time observe the relative influences of each loop on the behavior of variable of interest.
On the other hand, this example also shows that one needs to be careful in interpreting the results of the
EEA. The overall loop elasticity of the oscillatory behavior does not give much detail on how model
structure generates the behavior (Figure 13). One needs to consider the elasticities of the components of
the behavior (i.e. the exponential envelope and the frequency) separately in order to reach a better
understanding of the structure-behavior dynamics. Furthermore, the relations of loops with each other
also need to be considered. Once again, the formal methods should not be regarded as magical devises
that reveal all that needs to be known automatically but as tools that merely help to get a better
understanding of the dynamics of systems under study.
5 Recall that the magnitude of the complex elasticity value gives the overall elasticity of the oscillatory behavior;
hence, the non-negative loop elasticities in Figure 13.
21
A simple model of Economic Long Wave
The proposed methodology is also applied to a simple long wave model (Sterman 1985). The model is
chosen for a number of reasons. First it is widely known and therefore serves as a suitable test for the
proposed methodology. Second, although it has a relatively simple structure it is capable of generating
complicated dynamics because of the nonlinearities. Third, it was recently used in a concurrent
application of SILS and PPM approaches (Oliva and Mojtahedzadeh 2004). Therefore, it would be
interesting to see how the results of a joint application of SILS and the EEA compare to the findings of
that study. Finally, the model is also used by Kampmann (1996) and Ford (1999) to illustrate the
applications of their approaches. Thus, the application of the EEA on the long wave will also provide an
opportunity to compare the proposed approach with the Ford’s behavioral approach. It is assumed that the
reader is already familiar with the model. The details of the model and its dynamics can be found in
Sterman (1985). The stock-flow diagram is presented in Figure 14.
relative
SS orders
Supply io
Capital orders
12 LN
Capital
Acquisitions
Ly Depreciation [12
desired
orders
le L3 capacity
capital
adjustment
in
Ls
capacity
utilization supply
adjustment
Backlog desired \ E'S
Capital orders Production L4 capital 11
Backlog
d .
peeing desired
supply line
Figure 14. Stock-flow diagram of a simple Long Wave model.
1. Identify and list all nodes and choose the state variable of interest.
There are nineteen nodes in the model. The list of the nodes is in Table C.1 in Appendix C.3. Capital is
chosen as the state variable of interest.
2. Identify and list all causal links (except the ones that involve constants and flow-to-stock links) and
pathways in the model.
There are twenty-six such causal links in the model and they are listed in Figure C.4 in Appendix C.3.
There are thirty-six pathways between the state variables. The pathways originating from Capital, Supply
and Backlog are listed in Tables C.2-4, respectively, in Appendix C.3.
22
3. Identify and list the feedback loops in the Shortest Independent Loop Set (SILS).
There are sixteen loops in the set (Note that 34 {the total number of links} — 19 {total number of nodes} +
1 = 16. See Appendix A.2 for a brief discussion on this relationship). The loops are listed in Table C.5 in
Appendix C.3. The loops are also shown on the stock-flow diagram in Figure 14. The numbers of loops
that are most influential at one time or another during the simulation are printed with a larger font on the
figure. In addition, the corresponding pathways are depicted with thicker lines.
4. Form the gain matrix, which represents the links between state variables in their most compact form.
The model is a third-order nonlinear one. Therefore, the dimensions of its gain matrix are 3X3 and the
elements of the matrix change continuously over time. The form of the matrix is shown in Eq. 13. The
partial derivative equations are given in Eqs. C.1-9 in Appendix C.3.
2 Capital 2 Capital a Capital
dCapital OSupply OBacklog
G =| 9Supply od Supply od Supply (13)
lanigwave dCapital OSupply OBacklog ”
to) Backlog Cl) Backlog re) Backlog
dCapital OSupply OBacklog
5. Simulate the model and read the gain matrix, pathway gains, net flows of state variables over time.
The behavior of the variable of interest Capital is sustained oscillations with a period of about 45 years
(Figure 15). Therefore, the following analysis is conducted over a period of 46 years from year 136 to
182.
Capital
14
213] 4 & 1] 2 144 Ss 1} 2 43]4
ul
8
5
os =
2
100 110 120 130 140 150 160 170 180 190 200
Time ( Year )
Capital: LW unit
Figure 15. Reference behavior of state variable of interest, Capital.
6. Determine the characteristics of all elemental behavior modes of the system and the contribution of
each behavior mode to the behavior of the state variable of interest.
Like the Predator-Prey model, the observed behavior for this model is sustained oscillations. However,
the oscillations in the long wave are not caused by a single complex conjugate eigenvalue pair. Rather,
the oscillatory behavior is composed of various different behavior modes each of which becomes
23
gradually dominant for certain time segments of the simulation (Figure C.5 and Figure 18.a). There are
five distinct segments in a period; hence, the long-wave cycle will be divided into five phases similar to
the segmentation of the model behavior in Kampmann (1996).° These phases in sequential order are
1.“self-order initiation”, 2.“capital growth”, 3.“capital deceleration”, 4.“capital peak-then-decline”, and
5.“capital depreciation”. They are also shown on the behavior of Capital on Figure 15.
The Contributions of Elemental Behavior Modes to the Behavior of Capital
1
oof]
3a lS] 7 2
136 198 140 142 144 146 149 150 152 184 155 159160 162 164 166 168170 172 174 176178 180 182
Cea Ce
Figure 16. Rescaled contributions of behavior modes on behavior of state variable of interest, Capital.
The relative contributions (weights) of the behavior modes that are present at any time step to the overall
behavior of the variable of interest are provided on Figure 16. The behavior modes corresponding to each
phase are listed in Table 7. The evolution of these behavior modes over time is shown in Figure C.5. Note
that there are three convergent behavior modes in the beginning. Two of these modes disappear right after
year 138. At the same time, an oscillatory mode, which is initially almost instantaneously converging then
for the rest of the phase diverging, emerges (due to the coalescence of the two real positive eigenvalues to
give birth to a complex eigenvalue pair) and the “self-order initiation” phase begins. Then oscillatory
mode gives way to two behavior modes around year 142 (a consequence of the bifurcation of the complex
pair to two real negative eigenvalues) and the “capital growth” phase is initiated. One of these behavior
modes is divergent throughout the second phase while the other is convergent for most of the phase
except the beginning and ending of the phase. Towards year 150 these two behavior modes disappear.
They are replaced by another divergent oscillatory mode, which heralds the beginning of the “capital
deceleration” phase. Around year 152, the system enters the “capital peak-then-decline” phase when this
oscillatory mode too gives way to two convergent behavior modes. The fourth phase, lasts until about
year 157 and the dominant behavior mode is the one with a time constant that varies from about 6 years to
1 year during the phase. However, after a short transitional period the convergent behavior mode with the
time constant of about 20 years becomes dominant; hence, the final “capital depreciation” phase.
° Kampmann divided the long wave into four segments based on a visual inspection of the model behavior and the
changes in eigenvalues. However, calculating the contribution of each eigenvalue in this study reveals that there are
actually five distinct phases of the long wave.
24
Table 7. Behavior modes present in each phase of the long wave.
Phase no. Phase name Behavior modes
1 Self-order initiation One divergent oscillatory, one convergent mode
2 Capital growth One divergent, two convergent modes
3 Capital deceleration One divergent oscillatory, one convergent mode
Capital peak-then-decline | Three convergent modes
Ces
Capital depreciation Three convergent modes
7. Compute the eigenvalue elasticities with respect to compact links in the gain matrix.
The elasticity of each behavior mode is in matrix form as shown in Figure 17. For brevity, the graphs for
the elasticities are not shown.
evi jevl ev ev2jev2_ev2
a1 2 a3 a 42 As
el _| jevl jevl ev e2 _| jev2jevd ev
longwave =| ©21 22 a3 ‘longwave =| ©2122 a3,
evl jevl ev v2 jev2 ev?
31 32 33 63, G39 33
(a) (b)
evs jev3. ev eve neve eve
41 42 3 a1 42 3
evs | jev3. jev3 ev eve _| eve eve eve
longwave =| 21 22, 3 Eqmnewave =| 21 2 es
v3 4ev3 ev eve eve eve
es Gp 33 S323
(c) (d)
Figure 17. Elasticity matrices for three real eigenvalues (a, b, and c) and complex pair (d).
8. Compute the eigenvalue elasticities with respect to each causal link.
The set of equations that gives the elasticities with respect to each causal link are given in Eqs. C.10-35 in
Appendix C.3. For brevity, the graphs for the link elasticities are not shown.
9. Using the links identified in Step 3 and the loops in SILS form the directed cycle matrix.
The directed cycle matrix for this model is given in Figure C.4. It is a 16X26 matrix and its rank equals
16, the number of loops on the SILS.
10. Compute and plot the overall loop elasticity values over time and evaluate the findings.
The overall loop elasticities are computed and resulting loop dominance dynamics impacting the
exponential envelope and the observed frequency of oscillations over a whole cycle are shown in Figure
18.a. The loop dominance dynamics over each phase are given in Figures 18.b-f in sequential order.
As the last phase of the previous cycle gives way to the first phase of the next cycle the capacity falls
below orders. At this early stage of the self-order initiation phase, the Supply first-order control loop (L2)
becomes dominant briefly. The reason is that self-ordering mechanism begins to work and increases
Supply stock. However, acquisitions do not rise fast enough for capacity to catch up with orders and this
loop quickly loses its influence; the Capital self-ordering loop (L/4) becomes the most dominant. After a
brief period during which there is no markedly dominant loop the first phase ends (Figure 18.b). The
previous analyses identify the importance of the self-ordering loop in this phase but none of them
identifies the role played by the Supply first order control loop in the beginning (Kampmann 1996, Ford
1999, Oliva and Mojtahedzadeh 2004).
25
The Loop Dominance Dynamics of Long Wave Model
‘overal oop elasticity
us uM us
ie_uo
Lue
“averal op alasiclly
The Loop Dominance Dynamics of Long Wave Model - Phase 1
os
oa
#
7 18
us ug
=
uo
2 18 lal 1
oyu us ua
(a)
(b)
The Loop Dominance Dynamics of Long Wave Model - Phase 3
150 10.5 181 Isis
0-2 3
io uo nu
is 6
us us
78
usd
12
(d)
The Loop Dominance Dynamics of Long Wave Model - Phase 4
53
=o =
134 1551555
1325
1533
= ha Se —
ra
Ie uo Uist
U2 us ua
The Loop Dominance Dynamics of Long Wave Model - Phase 5
136158 10162 164186 1681701728178
=o <7 =a
ust
2G a ts He
tio uu ust
(e)
(p
Figure 18. Evolution of loop dominance dynamics for a whole cycle (a) and for phase | (b), 2 (c), 3 (d), 4
(e), and 5 (f) of Capital over time.
The capital growth phase is driven by the Capital expansion (L5) loop. In addition, the Economic growth
(L3) and the Order fulfillment (L8) loops have relatively large positive elasticities and hence they are also
important in creating the behavior of Capital. The influences of these loops are countered to a certain
extend by the Supply first order control (L2), the Capital decay (L/), and the Backlog expansion (L6)
26
loops. Nevertheless, their combined effect is not enough to match that of the loops with positive
elasticities and Capital continues to increase exponentially during this phase. The ital self-ordering
loop, however, becomes completely dormant. In addition, the capacity utilization is at its maximum
effectively keeping the first-order Production scheduling (L4) loop in check (Figure 18.c). This
explanation regarding the importance of the Capital expansion and the economic growth loops agrees
with that of Kampmann. However, the two do not agree in the details, the reason probably being that the
loop sets used in the two studies are not exactly the same. Interestingly, in his comparative analysis, Ford
confirms that the loop L3 is dominant but concludes that the loop L5 is not. He does not consider other
influential loops such as the overtime loop (the Production scheduling loop in this study) as mentioned in
Kampmann. This phase corresponds to the phases II, III, and IV in Oliva and Mojtahedzadeh (2004)
where they identify the Supply first-order control, the Capital expansion, and the Economic growth loops,
respectively, as the most influential structures. Their most influential structures agree with the ones
identified in this analysis. Nevertheless, different from the above analysis, they identify a single dominant
structure at each phase, a characteristic of the PPM approach. Another difference is that the loops
identified as the most influential structures in this study continue to be so throughout the second phase,
which is chopped into three more phases each dominated by another loop in Oliva and Mojtahedzadeh
(2004).
At the third, capital deceleration phase the net rate of increase in Capital slows down and even begins to
decrease towards the end of the phase. At the start of the phase, the relatively stronger influences of
Backlog expansion (L6) and Self-ordering (L/4) loops can be distinguished. Then Self-ordering loop
becomes the most influential one but looses its influence towards the end of the phase. In contrast, the
influences of Supply first-order, Order fulfillment, Economic growth and Production scheduling loops
show increase (Figure 18.d). This corresponds to the beginning of the decrease in the rate of increase in
Capital. This phase is part of phase III in Kampmann. The two studies agree in that the self-ordering is
the most influential on average over the whole interval but somewhat differs in the details. The same
conclusion holds regarding the Ford’s analysis. The capital deceleration phase coincides with the
beginning of the phase V of Oliva and Mojtahedzadeh (2004). They identify the Supply first order control
loop as the most influential during this phase.
The Production scheduling loop becomes dominant at the first half of the Capital peak-then-decline
phase. During this phase, the drop in Backlog that began in the previous phase continues; this relaxes the
capacity utilization and the Production scheduling loop becomes dominant. The loop exerts its influence
by drawing Capital stock downward thus it first slows down; then the Supply first-order loop takes the
place of the Production scheduling loop as the most dominant, stops the growth of the stock and then
causes it to drop precipitously right after it peaks. Nevertheless, this fast decline does not last long
because the dominance of the loop drops while that of capital decay loop increases enough to shift the
behavior of Capital to exponential decay. This happens around year 157. This phase coincides with the
second half of Kampmann’s phase III and most of the phase V of Oliva and Mojtahedzadeh (2004).
It is worth to note that Sterman also reports that when the capacity utilization is at its maximum the
Capital expansion (L5) is dominant. In time, as excess capacity develops, the capacity utilization drops
and dominance shifts to the Production scheduling loop which limits the amplitude of the cycle (footnote
26 on p. 45 of Sterman 1985). This line of events corresponds to the second phase through the first half of
the fourth phase of the EEA and except the third phase where the Self-ordering loop (L14) is dominant
the two explanations agree. They differ, however, in the corresponding dominant eigenvalues as discussed
in the next section.
The transition from the fourth phase to the fifth is not sharp. The loops exchange dominance gradually
from approximately year 156 to 158. During this phase, most of the feedback loops are shut off by the
nonlinearities in the system and Capital slowly falls under the influence of capital decay loop until the
cycle repeats itself. The findings for this phase agree with those from the previous studies.
27
5. Discussion and Conclusion
In spite of the guiding principle of the system dynamics that the behavior is generated as a result of the
interactions between various feedback mechanisms of the system, the accurate depiction of the
relationship between a model’s structure and its behavior still rests upon an experimental process of
hypotheses testing by repeated simulations. This approach is time consuming and prone to error on the
analyst’s part. Therefore, System Dynamics has much to gain from development of formal methods for
model analysis. The methodology proposed in this paper consists of ten steps and advances the
application of the EEA in loop dominance identification. The improvements proposed in this paper are
summarized in the following.
1. It is possible to relate the loop dominance dynamics directly to any selected state variable through its
net rate (i.e. slope). One criticism directed to the EEA has been its inability to relate the loop dominance
to a variable of interest. This improvement addresses this criticism.
2. An improved approach based on Oliva (2004) in identifying the independent loop set is incorporated
into the model analysis. The improved approach almost always finds a unique set and hence the problem
of arbitrary selection of the members of the set is alleviated.
3. A new loop dominance measure is proposed which takes all behavior modes of the model under study
into account.
4. For the selected models, the loop elasticities are actually plotted over time thus allowing for the
visualization of how loop dominance dynamics unfold over time.
5. Many parts of the proposed procedure are already codified. Tools for others are either available in
commercial software packages or it is the author’s opinion that they can be readily codified if sufficient
interest develops.
In addition, three models selected from the previous loop dominance studies were used to demonstrate the
implementation of the proposed methodology. This also facilitated the comparison of the proposed
version of the EEA with other formal model analysis approaches giving an opportunity to identify the
strengths and weaknesses of each approach. Such comparative studies should continue to be conducted in
the future in order to refine the existing formal analysis tools.
Of the three analyses, the one on the long wave deserves a closer look. It is generally believed that
complex eigenvalue pairs are responsible for the turning points in a cycle. In this respect, the EEA of the
long wave model leads to a rather counterintuitive observation on the generation of the long wave cycle.
For the minimum of the cycle the analysis confirms the mentioned belief: During the first phase, there is a
complex conjugate pair of eigenvalues and the pair is responsible for the transition of the Capital from
decay to exponential growth (see Figures 15-16).
The situation, however, is different for the maximum of the cycle. The Capital reaches its peak within
Phase four where there are three distinct negative eigenvalues. In order to understand the formation of the
peak we need to look at what happens earlier, i.e. in the phases two, three, and four. Figure 19 presents
the eigenvalues together with the Capital and its rate (slope), respectively, through these three phases. In
Phase two, there are three distinct real eigenvalues of which the largest positive one is the dominant (see
Figures C.5-16). The result is the ever-increasing rate, which causes the Capital to increase faster and
faster. Then two of the eigenvalues coalesce and form a complex conjugate pair of eigenvalues. The
complex pair is almost solely responsible for the behavior of the Capital in this third phase (see Figures
C.5-16). Their effect can be observed on the behavior of the rate of the Capital, the increase of which first
slows down; then, reverses and the slope begins to decrease (Figure 19.a). This, of course, exhibits itself
as an inflection point on the behavior of the Capital after which the increase in the Capital begins to slow
down (Figure 19.b). The third phase is, then in a sense where the complex pair prepares via its rate the
28
downturn of the Capital, which is to happen in the next, fourth phase. At around year 152, another
bifurcation occurs resulting in three distinct negative eigenvalues; hence the fourth phase. Now, under the
influence of the negative eigenvalues the rate exhibits decay until about year 157 which also marks the
transition from the fourth to the last phase of the cycle (Figure 19.a). The rate exhibits decay; however, as
long as it is positive the Capital continues to grow albeit slower and slower (Figure 19.b). Towards year
154, the rate becomes negative continuing its decay. At that instant, the Capital reaches its peak and then
begins to decrease. The behavior exhibited by the Capital until about year 157 is actually exponential
decline because its rate becomes more and more negative. In the last phase, the rate no longer drops but
this time exhibits decay in the opposite direction, i.e. towards zero. Consequently, the behavior of the
Capital changes from exponential decline to exponential decay.
What the EEA analysis tells us about the peak in the Capital is that the system prepares itself for the
peaking by entering and then leaving the oscillatory mode well in advance of the time the Capital peaks.
Rather the complex eigenvalue pair, with some influence from the negative eigenvalue, causes the peak in
its rate in the third phase. In the forth phase in which the peak in the Capital occurs, the system has
already left the oscillatory mode and now is in fact in a goal-seeking mode, which is evidenced by the
behavior of the rate of the Capital. However, the system is so charged during —to some extent- the second,
-but mainly- third phases that the crossing of the rate from positive to negative during its decay in the
forth phase causes the Capital to slow down its growth and then to decline, causing the peak in the
Capital.
Although the loops identified being dominant in this study and that of Sterman’s mostly agree the two
differ is their findings regarding the dominant behavior modes (i.e., eigenvalues), especially in the second
through forth phases of the cycle. According to the EEA, the dominant behavior modes for the first and
third phases are divergent oscillations (i.e., complex conjugate eigenvalue pairs with real parts). These
phases correspond to the beginning and end of the time interval during which the capacity utilization is at
its maximum. The rest of this interval constitutes the second phase and the dominant behavior mode is
monotonic divergent (There are three real distinct eigenvalues the most positive of which is the
dominant). The second and third phases correspond to what Sterman termed “the expansion phase” for
which he mentions only the divergent oscillatory behavior mode (i.e., a complex eigenvalue pair with real
parts) as dominant. In the forth phase as capacity utilization drops, there are three monotonic convergent
behavior modes (i.e., real and distinct negative eigenvalues). Therefore the dominant behavior mode
during the peak of the Capital is not a damped oscillation as suggested in Sterman but a combination of
exponential decay modes. In other words, the amplitude of the Capital is limited by monotonic
convergent behavior modes, not by a highly damped oscillation.
One of the problems of the earlier EEA applications was the lack of an explicitly specified variable of
interest in the analysis. Specifying a variable of interest makes the results of the analysis more rigorous
and relevant. In a typical EEA study, only one behavior mode is regarded as dominant at each analysis
time step. Hence, the resulting explanation on loop dominance would be based on that behavior mode. In
this study all behavior modes are considered to the extend of their contribution to the behavior of the
variable of interest at each time step and a conglomerate measure for loop dominance is devised based on
this approach. Finally, the previous applications of the EEA fell short of providing results that show the
continuous and gradual change in loop dominance dynamics over time. That gap is also closed in this
study. The differences between the proposed methodology and the earlier applications of the EEA are
demonstrated on the long wave model analysis.
29
The Elemental Behavior Modes and the Rate of the Capital
pos
4 15
2
&
unit/year
“0.5
: 3
-l -1
146 48 150
Time (year)
[—evt —- ev2 ----ev3 — evi —rate]
(a)
The Elemental Behavior Modes and the Capital
eS
a
x
ese
ROS
ieee)
es
VW
0
301 h—- 10
B02 \ =
g i 9 =
3 03 5
2 \
@ 0.4 os 8
0.5
\ 7
-0.6 ¥
-0.7 \. 6
-0.8
09 2 3 5
Ai i 4
146 148 150 152 154 156 158 160
Time (year)
[ev —— ev2 ----ev3 — evi ~- Capital]
(b)
Figure 19. The behavior modes (eigenvalues), the rate of change in Capital (a) and Capital (b) in the third,
fourth phases and parts of second and fifth phases.
There is a close relationship between the eigenvalues and the total participation metric in linear systems.
Mojtahedzadeh er al. 2004 has shown that for any linear system regardless of the order of the system the
total participation metric is equal to the largest eigenvalue of the system in the steady-state condition.
30
Still, the most important difference between the EEA and the PPM approaches is that the first approaches
structure-behavior relations from a system-wide perspective while the latter adopts a localized
perspective. The results one gets from the two methods may not always be in agreement because of the
different perspectives the two approaches take on in loop dominance analysis. For instance, although in
the Predator-Prey example the explanations from the two approaches are mostly in agreement it seems the
PPM approach do not appropriately relate relevant parts of the model structure to the model behavior in
case of oscillatory behavior. In another comparative study, a linear oscillatory model was used and there
was even more disagreement between the results of the two approaches over the structural source of
oscillations (results not shown). It may be the case that the level of agreement between the two is higher
when they are applied to nonlinear models. It is, however, worth mentioning that the application of both
EEA and PPM on a simple long wave model resulted in the same loops as the most influential ones in the
generation of the cycle. Moreover, these loops are the ones selected by the loop selection algorithm
(SILS) proposed by Oliva (2004). This suggests that it may be possible for other —more complicated—
system dynamics models to reach an intuitive understanding of model dynamics (including the loop
dominance dynamics) based on the loops selected in the SILS of those models.
Although there have been notable improvements in the EEA as a formal model analysis tool since it was
first advocated more than twenty years ago, there is still need for further refinement of the methodology.
For instance, the proposed overall loop dominance measure performed well in the examples but their real
utility will show itself only after many applications on various system dynamics models. SILS approach
is no doubt useful in loop dominance analysis but needs to be used with scrutiny and may require suitable
revisions in the future. More importantly, there is still much to be done for fully computerized
implementation of EEA routines. As an example, setting up equations for the causal link elasticities
seems to be the routine that is most time-consuming. The successful automation of this step would greatly
reduce the duration of analysis. In addition, most of the codes presented are either tailored for models of a
specific order or for a specific model. While useful the set of codes presented in this study only serves as
a starting point and more elaborate ones are required for the proper computerized implementation of EEA.
The work presented in this paper is part of a larger study. The next phase will focus on understanding how
uncertainty in parameter values in a model affect its loop dominance dynamics and will look into devising
ways to differentiate the contribution of each uncertainty source to the uncertainty of dominance
dynamics by coupling the EEA approach with the error analysis.
Finally, it should be emphasized that existing approaches to loop dominance analysis may have much to
learn from each other. As a result, it would not be surprising to witness the integration of best aspects of
the different formal tools in the future.
References
Davidsen P. 1991. The structure-behavior graph. The System Dynamics Group, MIT, Cambridge.
Eberlein R. 1984. Simplifying Dynamic Models by Retaining Selected Behavior Modes. PhD Dissertation,
MIT, Cambridge, MA.
Ford D. 1999. A behavioral approach to feedback loop dominance analysis. System Dynamics Review
15(1): 3-36.
Forrester N. 1982. A Dynamic Synthesis of Basic Macroeconomic Theory: Implications for Stabilization
Policy Analysis. PhD Dissertation, MIT, Cambridge, MA.
. 1983. Eigenvalue analysis of dominant feedback loops. The 1983 International System Dynamics
Conference, Plenary Session Papers, pp. 178-202. System Dynamics Society, Albany, NY.
Franklin GF, Powell JD, Emami-Naeini A. 2002. Feedback Control of Dynamic Systems. Fourth Edition.
Prentice Hall Inc., NJ.
31
Gertner GZ. 1987. Approximating precision in simulation projections: An efficient alternative to Monte
Carlo methods. Forest Science 33: 230-238.
Gongalves P, Lertpattarapong C, Hines J. 2000. Implementing formal model analysis. Proceedings of the
2000 International System Dynamics Conference, Bergen, Norway. System Dynamics Society,
Albany, NY.
Graham AK. 1977. Principles of the Relationship Between Structure and Behavior of Dynamic Systems.
PhD Dissertation, MIT, Cambridge, MA.
Giineralp B. 2004. A principle on structure-behavior relations in system dynamics models. Proceedings of
the 2004 International System Dynamics Conference, Oxford, UK. System Dynamics Society, Albany,
NY.
Kampmann C. 1996. Feedback loop gains ans system behavior. Proceedings of the 1996 International
System Dynamics Conference, Boston, MA. System Dynamics Society, Albany, NY.
Mojtahedzadeh MT. 1997. A Path Taken: Computer Assisted Heuristics for Understanding Dynamic
Systems. PhD Dissertation, University at Albany, SUNY, Albany, NY.
Mojtahedzadeh MT, Andersen DF, Richardson GP. 2004. Using Digest to implement the pathway
participation method for detecting influential system structure. System Dynamics Review 20(1): 1-20.
Oliva R. 2004. Model structure analysis through graph theory: partition heuristics and feedback structure
decomposition. System Dynamics Review 20(4): 313-336.
Oliva R, Mojtahedzadeh MT. 2004. Keep it simple: dominance assessment of short feedback loops.
Proceedings of the 2004 International System Dynamics Conference, Oxford, UK. System Dynamics
Society, Albany, NY.
Richardson GP. 1995. Loop polarity, loop dominance, and the concept of dominant polarity. System
Dynamics Review 11(1): 67-88.
—— . 1996. Problems for the future of system dynamics. System Dynamics Review 12(2): 141-157.
Saleh MM. 2002. The Characterization of Model Behavior and its Causal Foundation. PhD Dissertation,
University of Bergen, Bergen, Norway.
Sterman JD. 1985. A behavioral model of the economic long wave. Journal of Economic Behavior and
Organization 6: 17-53.
32
Appendix A.
A.1. Mathematical foundation of the eigenvalue elasticity analysis
The causal structure of a linear (or linearized) model can be represented as a gain matrix. Let G be the
gain matrix of ann" order system dynamics model. Then G will be an nXn square matrix. The entries of
the gain matrix are the net gains of compact links between the state variables (stocks) of the model to be
analyzed. Forrester (1982) proposed using the eigenvalues of the gain matrix as a bridge between model
structure and behavior. Then the model can be represented in matrix form:
x=Gx (Al)
where X is the vector of the net rate (flow) of state variables and x is the state variables vector.
While the eigenvalues of the gain matrix G stand for the elemental behavior modes of the system the
eigenvalue elasticities reflect the influence of different parts of model structure on these behavior modes.
Saleh (2002) showed how overall model behavior could be partitioned into its components on the
eigenspace. The departure point in this approach is differentiating both sides of Eq. A.1, which leads to
the observation that the behavior pattern of any state variable is determined by its time derivative (slope)
and its second time derivative (curvature).
c=Gs (A.2)
where ¢ is the curvature (the second derivative) of x and s is the net rate (slope) (the first derivative) of x.
The x in Eq. A.1 constitutes the slope vector of the state variables. In the eigenspace, the slope vector
can be expressed as a linear combination of the right eigenvectors of the gain matrix, G.
a
s=)'ar, (A3)
isl
In the eigenspace, @, are the new components of the slope vector, s.
Differentiating Eq. A.3 gives:
c (A.4)
Now, in the eigenspace the @, are the new components of the curvature vector.
Substituting Eq. A.3 into Eq. A.4 and utilizing the fact that Gr, = Ay, results in
(A.5)
Then along a particular coordinate (spanned by a right eigenvector) the dynamics that unfolds can be
described by,
on (A6)
A.l
The solution of Eq. A.6 is then,
a, =arer i=lin (A.7)
where f, is the initial time, @ is the initial value of @, at t,,
Eqs. A.2 through A.7 shows that the only factor determining the dynamics along a particular coordinate
on the eigenspace is the eigenvalue associated with that coordinate.
Finally, substituting Eq. A.7 into Eq. A.3 yields the time trajectory of the slope:
s= diate At (A.8)
Thus, the slope trajectory is decomposed into several behavior modes, each expressed by an eigenvalue
and its associated right eigenvector.
The elasticity of an eigenvalue with respect to a variable measures the percentage change in the
eigenvalue for a given percentage change in that variable. The sensitivity S,,,; of an eigenvalue 2, with
respect to a variable pq is given by the partial derivative of that eigenvalue with respect to the variable
(Eq. A.9). The eigenvalue elasticity is then defined as the sensitivity of the eigenvalue with respect to the
variable normalized for the size of the variable and the size of the eigenvalue (Eq. A.10). As shown in the
equation, eigenvalue elasticity can be computed using left and right eigenvectors and the partial derivative
of the linear system matrix, G, with respect to gain g,,. The partial derivative of G with respect to a
variable can be found by calculating the matrix before and after a small change in the parameter. Another
way is to derive the expressions for the partial derivative of G and using them to compute elasticity
values.
0A,
pai = ie (A.9)
8 pg
_ Sia OG, Soa
es Soa eM (A.10)
where €y9,i =the elasticity of eigenvalue i, A, to compact gain link pq, &pq
A, = eigenvalue i
8pq = gain of the link from p to g
1 ranspose of the i'” left eigenvector (1 Xn vector)
G_ = linear(ized) system matrix (nXn)
r= 7" right eigenvector (nx 1).
A more useful formulation of eigenvalue elasticity is, however, possible without having to calculate the
partial derivative of the gain matrix, G (Eq. A.11).
es =H, (p) 4m (q) #22
th
(AL)
tl
where 1, (p) the p" element of the i left eigenvector (1 Xn vector)
= the q" element of the i” right eigenvector (1 Xn vector)
wt
=
i]
A2
Eq. A.11 comes from the relation that the sensitivity matrix §; of the eigenvalue /, is equal to the product
of the i" left eigenvector and the i right eigenvector of the gain matrix, G (Eq. A.12). The theorem
behind this relation and its proof is given in Saleh (2002).
S,=1-1 (A.12)
The compact link elasticities are not much of use by themselves in revealing the dominant structure. They
are, by definition, compact giving no hint of the detailed model structure. One needs to make use of the
elasticities to individual causal links in the model to bring in the structural detail into the analysis. It turns
out that it is possible to obtain causal link elasticities using pathway gains and compact link elasticities.
The general expression for causal link elasticities is
(A.13)
where e,, = the elasticity of eigenvalue i, A, to causal link J;
19 = the elasticity of eigenvalue i, A, to compact gain link pq, &pq
Spans) = gain of pathway r(s) from state variable p to q
= the number of pathways causal link J; is a part of
R = the total number of pathways between state variables p and g
A.2. Directed Cycle Matrix
A directed cycle matrix, in system dynamics context, reflects the membership of causal links (except the
ones that involve constants and flow-to-stock links) to the feedback loops of a model (Kampmann 1996).
If only the loops in an “independent loop set” as in Kampmann or the “shortest independent loop set” as
in Oliva (2004) is used then the columns of the directed cycle matrix are linearly independent. Thus one
can show the relation of links with the loops as in Eq. A.14.
4 Ip,
1, |=D| Ip, (A.14)
I Px
. dye.
D= ij (A.15)
where dj, equals | if link j is an element of loop k and 0 otherwise.
A3
Finally, the set of equations that relates the link elasticities to loop elasticities is
Gs lip,
e, |=D] ey, (A.16)
4, op
fi
where ¢, is the elasticity of link Jj, ¢,, is the elasticity of loop /p; and D is the directed cycle matrix as
shown in Eq. A.15.
The number of links is almost always larger than the number of loops in system dynamics models. This
means that the directed cycle matrix will be overdetermined. In fact, it is shown by Kampmann that the
number of independent loops is equal to (total number of links — total number of nodes + 1) and in a
typical system dynamics model the number of links far exceeds the number of nodes. Refer to Kampmann
(1996) for a proof of this relationship. The linear independence of columns in a directed cycle matrix,
however, allows one to obtain a unique solution for loop elasticities even when the matrix is
overdetermined.
A4
Appendix B.
Pseudo-codes for eigenvalue elasticity analysis implementation
B.1. Pseudo-code for the main function
The main function serves as an intermediary between Vensim and MATLAB. The pseudo-code of the last
MATLAB function called in the main function, tofile.m is not given below. It simply writes the results of
the analysis to files.
ep < engOpen(NULL)
dt ~ TIMESTEP
ord — ORDER
for time = 0: FINAL_TIME*(1/TIMESTEP)
timein — time * TIMESTEP
for i = 1: total number of compact gains and pathway gains
vensim_get_sens_at_time("data.vdt", "gain(i)’, "Time", timein, gain(i), 1)
end
for i = 1: total number of state variables (stocks)
vensim_get_sens_at_time("data.vdt", "slope(i)”,
"Time", timein, slope(i), 1)
end
t © timein
ge gain
sip < slope
call Evcont.m
call Elast.m
call LoopElast.m
call tofile.m
end
engClose(ep)
B.2. Pseudo-code of the function Evcont.m
start MATLAB connection
read timestep into MATLAB
read order of model into MATLAB
for each timestep throughout simulation run
actual time in simulation
for each compact and pathway gain
read simulated gain data from Vensim;
populate gain matrix
end
for each stock
read simulated slope data from Vensim;
populate slope matrix
end
read current time into MATLAB
read gain matrix into MATLAB
read slope matrix into MATLAB
call MATLAB function Evcont.m
ns " Blast.m
LoopElast.m
tofile.m
end time loop
close MATLAB connection
The function calculates the eigenvalues, eigenvectors, initial alphas and the contribution of each
eigenvalue on the behavior of the selected state variable. Its pseudo-code tailored for a third-order system
is
given below.
function Evcont.m
[v,d] — eig(G)
viny < inv(y)
alp < vinv * slp
voi < {1,2,0r3}
if (imag(d) == 0)
s1(1) < alp(1) * v(voi,1) * exp(d(1,1) * (0-0)
s2(1) — alp(2) * v(voi.2) * exp(d(2.2) * (0-0))
$3(1) — alp(3) * v(voi,3) * exp(d(3.3) * (0-0))
s1(2) — alp(1) * v(voi.1) * exp(d(1,1) * (dt-0))
s2(2) — alp(2) * v(voi.2) * exp(d(2.2) * (dt-0))
$3(2) — alp(3) * v(voi.3) * exp(d(3.3) * (dt-0))
diff_s(1) — s1(2) - s1(1)
diff_s(2) — s2(2) - 21)
diff_s(3) — s3(2) - s3(1)
absdiff — abs(diff_s(1)) + abs(diff_s(2)) + abs(diff_s(3))
for i= 1:3
AS
calculate eigenvalues and (right) eigenvectors
take inverse of eigenvector matrix
calculate initial alphas
select variable of interest
if all eigenvalues are real
calculate initial value of slope component along
eigenvector associated with first eigenvalue
calculate initial value of slope component along
eigenvector associated with second eigenvalue
calculate initial value of slope component along
eigenvector associated with third eigenvalue
calculate final value of slope component along
eigenvector associated with first eigenvalue
calculate final value of slope component along
eigenvector associated with second eigenvalue
calculate final value of slope component along
eigenvector associated with third eigenvalue
calculate change along first slope component
i ” "second
gs " "third"
calculate sum of magnitudes of changes
for each eigenvalue
cont(i) — diff_s(i) / absdiff calculate its contribution
end end
else if there is complex (eigenvalue) pair
betal < real(alp(1))*real(v(voi,2))-abs(imag(alp(1)))*abstimag(v(voi,2))) calculate the coefficient of cosine component
gammal < real(alp(1))*abstimag(v(voi.2)))+ abs(imag(alp(1)))*real(v(voi,2)) “ @ "sine z
sr(1) — real(alp(3))*real(v(voi,3))*exp(real(d(3,3))*(0-0)) calculate initial value of slope component along
eigenvector associated with real eigenvalue
se(1) — 2*exp(real(d(1,1))*(0-0))*(betal *cos(abs(imag(d(1,1)))*(0-0)*((2*pi/360)) calculate initial value of slope component along
-gamma | *sin(abs(imag(d(1,1)))*(0-0)*((2*pi)/360))) eigenvector associated with complex pair
sr(2) real(alp(3))*real(v(voi,3))*exp(real(d(3,3))*(dt-0)) calculate final value of slope component along
eigenvector associated with real eigenvalue
se(2) — 2*exp(real(d(1,1))*(dt-0))*(betal *eos(abs(imag(d(1 ,1)))*(dt-0)*((2*pi)/360)) calculate final value of slope component along
-gamma | *sin(abs(imag(d(1,1)))*(dt-0)*((2*pi)/360))) eigenvector associated with complex pair
diff_real — sr(2) - se(1) calculate change due to real eigenvalue
diff_cmp — se(2) - se(1) bs ” "complex pair
cont_real < diff_sreal / (abs(diff_real) + abs(diff_cmp)) calculate contribution of real eigenvalue
cont_cmp € diff_scmp / (abs(diff_real) + abs(diff_cmp)) cE ” complex pair
end end
B.3. Pseudo-code of the function Elast.m
The function calculates the compact link elasticities of each eigenvalue regardless of the order of the
model under study.
function Elast.m
get left eigenvectors
for each eigenvalue
for each row
for each column
Elij.ev) — Mev) * vinv(jev) * Gij) / dlev,ev) calculate the compact link elasticity values
end end
end end
end end
B.4. Pseudo-code of the function LoopElast.m
The function calculates the elastic’ of each eigenvalue to the links and then to the loops in the Shortest
Independent Loop Set (SILS). It is modified for the analysis of the simple long wave model.
function LoopElast.n
The for loop below calculates, for each eigenvalue, the elasticities to causal links (except the ones that involve constants and flow-to-stock links)
based on the compact link elasticities
for i= 1:3
el(1.i) — E(1,3,i) * (g13p! / G(1,3)) + EQ,3,/) * (g23p1 / G2,3))
el(2,i) — E(3,3,/) * ((g33p1 + g33p3 + g33p4) / G(,3)) + E(2,3,i) * ((g23p2 + 223p3 + 223p4) / G(2,3)) + E(L.3,i) * (g13p2/ G(L,3))
el(3,i) — E(2,3,i) * (g23p5 / G(2,3)) + EG,3,/) * (233p2/ G3.3))
el(4.i) — ECL) * (gl 1p2/ G(1,1)) + EQ,1,i) * (g21p2 + g21p8) / G(2,1)) + EG,L.i) * ((g31p2 + g31p8) / GG,1))
el(5,i) — ECL * (gl Ip! / G(1)) + EQ,1,i) * (g21pl + g21p7) / G2,1)) + EG,Li) * ((g31p! + 231p7) / GG,1))
el(6,i) — E(1,Li) * (gl 1p2/ G(1,1)) + E(L,3,/) * (g13p2/ G(,3)) + EQ.1.i) * ((g21p2 + g21p8)/ GQ,1))
+E(2,3,i) * ((g23p2 + g23p3) / G(2,3)) + EL.) * ((g31p2 + 231p8) / G1) + EG.3.i) * ((g33p1 + 233p3) / GB,3))
el(7,i) — (1,1) * (gl ipl + gl 1p2)/ GU,1)) + EQ,1) * ((g2Lp1 + g21p2 + g21p7 + 221p8)/ G2,1))
+EQG,1,i) * (g31pl + g31p2 + g31p7 + g31p8)/GG,1))
el(8,i) — E(2,1,) * (g21p4 / G(2,1)) + EG,L,i) * (g31p4/ GG,1)
el(9,i) — E(1,1i) * (gl 1p3 / G(L,1)) + EQ,Li) * ((g21p3 + g21p5 + g21p6 + g21p9) / G2,1))
+E(3,1i) * ((g31p3 + g31p5 + g31p6 + 231p9) / GG,1))
el(10,i) — EQ,1i) * (g21p4/ G2,1)) + EQ,3,i) * (g23p4 / G(2,3)) + EG,L,) * (g31p4 /GG,1)) + EG,3,/ * (233p4/ G3,3))
el(1 1,1) — EG.1i) * (g31p3 + g31p4 + g31p5 + g31p6 + g31p7 + g31p8 + g31p9) / GG,1)) + EG.2./)
+ E(3,3./) * ((g33p2 + 233p3 + g33p4) / GG.3))
el(12,i) — EQ,1,i) * (g21p3 / GQ,1)) + EG,1,i) * (231p3/GG,1)
A6
el(13,) — EQ2,1,i) * (g21p5 / G(2,1)) + EG,1,i) * (231p5 /GG,1)
el(14,) — EQ2,1i) * (g21p6 / G(2,1)) + EG,1,i) * (g31p6 / G3,1))
el(15,) — EQ2,1,i) * (g21p9/ G(2,1)) + EG,L,i) * (g31p9/GG,1))
el(16,i) — E(2,3,i) * (g23p4/ G(2,3)) + EG.3,/) * (233p4/ GG,3))
el(17,i) — EQ2,1i) * (g21p4 + g21p5 + g21p6 + g21p7 + g21p8) / G(2,1)) + E(2,2.i) * (22p1 / G2.2))
+ E(2,3.i) * ((g23p3 + 223p4 + 223p5) / G(2,3)) + EG,Li) * (g31p4 + g31p5 + g31p6 + g31p7 + g31p8)/ G,1)) + EG.2./)
+ E(,3,/) * ((g33p2 + 233p3 + £33p4) / GG.3))
el(18,) — E(1,3,i) * (g13p2/ G(,3)) + EQ,3,i) * ((g23p2 + g23p3) / G(2,3)) + E3.3.i) * ((g33p1 + g33p3) / G3,3))
el(19,) — EQ2,3,i) * (g23p4/ G(2,3)) + EG,3,i) * (233p4/ G3,3))
€l(20,i) — E(2,1,i) * (g21p6 + g21p7 + g21p8) / G(2,1)) + EG, 1.) * ((g31p6 + 231p7 + 231p8) / GB,1))
+ E(2,3,i) * ((g23p3 + £23p5) / G(2,3)) + EG.3,i) * ((g33p2 + g33p3) / GG,3))
el(21,) — E(1,1i) * (gl ipl + gl1p2)/G(,1)) + E(2,1,0) * (g2Ip! + g21p2) / G2,1)) + E(1,3,) * (g13p2/ G13)
+ E(2,3.i) * (g23p2/ G2,3))
el(22,i) — EQ2.1,i) * ((g21p7 + g21p8)/ GQ,1)) + EG.) * ((g31p7 + 23 1p8) / GB.) + EQ.3,i) * (223p3 / G(2,3))
+EG,3,i) * (933p3/ GG,3))
€l(23,i) — EQ2,1.i) * ((g21p4 + g21p5 + g21p6 + g21p7 + g21p8 + g21p9) / G(2,1)) + E(2,2,i) * (g22p1 / G2.2))
+ E(2,3,/) * ((g23p3 + g23p4 + g23p5) / G(2.3)) + EG.1.i) * ((g31p4 + g31p5 + g31p6 + g31p7 + g31p8 + g31p9) / GG,1))
+ E(3.2.i) + EG.3.i) * ((g33p2 + 233p3 + 933p4) / GB.3))
el(24,i) — E(1,2,i) + EQ.2,i) * (g22p2 / G(2,2))
el(25,i) — EQ,2,i) * (g22p1 / G2,2)) + EG.2,i)
€l(26,i) — EQ2,1,i) * ((g21p6 + g21p7 + g21p8) / G(2,1)) + E2,2,i) * (g22p1 / G(2,2)) + (2.3.1) * ((g23p3 + 223p5) / G(2,3))
+EQG,1i) * (g31p6 + g31p7 + g31p8) / GG,1)) + EG.2.i) + EG,3,i) * ((g33p2 + 233p3) / GB,3))
end
dem — create directed cycle matrix
fo 0 0 60 0 1 0 1 0 0 0 0 0 0 0 0
o 0 0 1 0 60 0 0 60 0 0 0 0 1 0 0
o 0 0 0 0 60 0 60 60 0 1 0 0 0 0 0
o 0 0 0 0 0 0 60 1 0 060 0 0 0 0 0
oo 1 0 0 0 0 1 00 0 0 0 0 0 1
oo 0 1 60 0 0 60 1 0 0 0 0 0 0 0
oo 1 0 0 0 0 1 1 0 0 0 0 0 0 1
o 0 0 0 60 6 0 0 60 0 60 0 1 0 0 0
1 0 0 0 1 1 0 0 0 1 0 1 0 0 1:0
o 0 0 0 00 0 00 0 0 0 1 1 0 0
o 0 0 0 60 1 0 0 0 0 1 0 0 1 0 0
oo 0 0 1 1 0 0 0 00 0 0 0 0 0
o 0 0 0 0 60 0 0 60 0 0 1 0 0 0 0
o 0 0 0 0 060 0 60 60 0 60 0 0 0 1 0
o 0 0 0 0 0 0 00 1 60 0 0 0 0 0
o 0 0 0 60 0 0 060 60 0 60 0 0 1 0 0
o 0 0 0 0 0 1 0 0 0 1 1 1 1 1 1
o 0 0 1 60 6 0 0 60 060 60 0 60 0 0 0
o 0 0 0 60 0 0 00 0 0 0 0 1 0 0
o 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1
oo 1 0 00 0 0 1 0 0 0 0 0 0 0
o 0 0 0 00 0 0 0 00 0 0 0 0 1
oo 0 0 0 0 1 0 0 1 1 1 1 1 1 1
o 1 0 0 1 0 0 00 1 0 21 1 oo 1 1
o 0 0 0 0 0 1 60 60 60 60 60 0 0 0 0
o 0 0 0 0 0 1 06 0 0 1 0 0 0 1 yy
. for each eigenvalue
2) = dem \ el(:,i) calculate elasticities to loops in SILS
end
for i= 1:16 for each loop in SILS
overall_elps(i) — cont(1) * elps(i,1) + cont(2) * elps(i,2) + cont(3) * elps(i,3) calculate overall elasticity
abs_elps += abs(overall_elps (1)) + abs(overall_elps 2)) + abs(overall_elps(3)) calculate sum of magnitudes
end end
for i= 1:16 for each loop in SILS
res_overall_el(i) = overall_elps(i) / abs_elps calculate rescaled overall elasticity
end end
AT
Appendix C.
C.1. Yeast model
Equations of Yeast model (based on Vensim notation)
TIME STEP = 0.01; Integration method: Euler
Stocks
Cells = INTEG (births-deaths, 1)
Alcohol = INTEG (alcoholgeneration, 0)
Flows
births
(Cells/divisiontime)*eff
alc birth
deaths = (Cells/lifetime)*eff alc death
alcoholgeneration = Cells*alcoholpercell generation
Auxiliaries
eff alc birth = (-0.1*Alcohol)+1.1
eff alc death = EXP(Alcohol-11)
lifetime = 30
divisiontime = 15
alcoholpercellgeneration = 0.01
The change in the eigenvalues of Yeast model over time
The Eigenvalues of Yeast Model
Phase I
Phase 2
—aevl —— ev2----evi
Phase 3
Phase 4
Figure C.1. Eigenvalues of Yeast model.
AB
C.2. Predator-Prey model
Equations of Predator-Prey model (based on Vensim notation)
TIME STEP = 0.01 months; Integration method: Euler
Stocks
Predator = INTEG (predator birth-predator death, 1)
Prey = INTEG (prey birth-prey death, 2)
Flows
predator death = predator death rate*Predator
predator birth = predator interaction constant*Prey*Predator
prey birth = prey birth rate*Prey
prey death = prey interaction constant*Prey*Predator
Auxiliaries
predator interaction constant = 0.1
prey interaction constant = 0.2
predator death rate = 0.15
prey birth rate = 0.35
An alternative stock-flow representation of Predator-Prey model where the major loop L5 is
concealed in the structure
{i ul
Prey
wey bit prey death
prey birth 12 a4
™ L ———/ prey interaction
interaction constant
Predator
predator birth
predator death
if
‘ L4 x.
predator predator death
interaction constant ~
rate
Figure C.2. An alternative stock-flow representation of Predator-Prey model.
Ad
The change in real and imaginary parts of the complex eigenvalue pair of Predator-Prey Model
The Complex Eigenvalue Pair of Predator-Prey Model
Phase | Phase\2 Phase 3
jaginary
Figure C.3. Complex eigenvalue pair of Predator-Prey model.
AO
C.3. Simple Long Wave model
Equations of Long Wave model (based on Vensim notation)
TIME STEP = 0.25 years; Integration method: Euler
Stocks
Capital = INTEG (Acquisitions-Depreciation, (capital output ratio*avg lifetime of capital)
/(avg lifetime of capital-capital output ratio))
Backlog = INTEG (Capital orders Backlog+goods orders-Production, normal delivery delay)
Supply = INTEG (Capital orders-Acquisitions, (Backlog/Production)*Depreciation)
Flows
Capital orders = Depreciation*relative orders
Acquisitions = (Supply*Production)/Backlog
Depreciation = Capital/avg lifetime of capital
Capital orders Backlog = Capital orders
Production = capacity*capacity utilization
goods orders = 1
Auxiliaries
capacity = Capital/capital output ratio
capacity utilization = capacity utilization fnc(desired production/capacity)
capital adjustment = (desired capital-Capital)/capital adjust time
desired capital = desired production*capital output ratio
desired orders = Depreciation+capital adjustment+supply adjustment
desired production = Backlog/normal delivery delay
desired supply line = Depreciation*(Backlog/Production)
relative orders = relative orders fnc(desired orders/Depreciation)
supply adjustment = (desired supply line-Supply)/supply adjust time
capacity utilization fne(
[(0,0)-(2,1.))],
(0,0),(0.2,0.3),(0.4,0.6),(0.6,0.8),(0.8,0.9),(1,1),(1.2,1.03),(1.4,1.05),(1.6,1.07),(1.8,1.09),(2,1.1)
relative orders fnc(
[(-1,0)-(40,6)],
(-1,0),(-0.5,0),(0,0.2),(0.5,0.5),(1,1),(1.5,1.5),(2,2),(2.5,2.5),(3,3),(3.5,3.5),(4,4),(4.5,4.4),(5,4.8),(5.5,5.2),
(6,5.5),(6.5,5.65),(7,5.7),(7.5,5.75),(8,5.8),(40,6))
avg lifetime of capital = 20
capital adjust time = 1.5
capital output ratio = 3
normal delivery delay = 1.5
supply adjust time = 1.5
The list of nodes, causal links and loops of Long Wave model
Table C.1. Nodes of Long wave model.
Nodes
Acquisitions capital adjustment desired production
Backlog Capital orders desired supply line
capacity Capital orders backlog Production
capacity utilization Depreciation relative orders
capacity utilization fne desired capital relative orders fne
Capital desired orders Supply
supply adjustment
Table C.2. Pathways originating from Capital.
Pathway no. Variable sequence
cle Capital, capacity, Production, Acquisitions, Capital
c2c Capital, capacity, capacity utilization, Production, Acquisitions, Capital
c3c Capital, Depreciation, Capital
cls Capital, capacity, Production, Acquisitions, Supply
c2s Capital, capacity, capacity utilization, Production, Acquisitions, Supply
c3s Capital, Depreciation, Capital orders, Supply
c4s Capital, capital adjustment, desired orders, relative orders, Capital orders, Supply
cs Capital, Depreciation, desired orders, relative orders, Capital orders, Supply
c6s Capital, Depreciation, desired supply line, supply adjustment, desired orders, relative
orders, Capital orders, Supply
c7s Capital, capacity, Production, desired supply line, supply adjustment, desired orders,
relative orders, Capital orders, Supply
c8s Capital, capacity, capacity utilization, Production, desired supply line, supply
adjustment, desired orders, relative orders, Capital orders, Supply
c9s Capital, Depreciation, relative orders, Capital orders, Supply
clb Capital, capacity, Production, Backlog
c2b Capital, capacity, capacity utilization, Production, Backlog
c3b Capital, Depreciation, Capital orders, Capital orders backlog, Backlog
c4b Capital, capital adjustment, desired orders, relative orders, Capital orders, Capital
orders backlog, Backlog
c5b Capital, Depreciation, desired orders, relative orders, Capital orders, Capital orders
backlog, Backlog
c6b Capital, Depreciation, desired supply line, supply adjustment, desired orders, relative
orders, Capital orders, Capital orders backlog, Backlog
c7b Capital, capacity, Production, desired supply line, supply adjustment, desired orders,
relative orders, Capital orders, Capital orders backlog, Backlog
c8b Capital, capacity, capacity utilization, Production, desired supply line, supply
adjustment, desired orders, relative orders, Capital orders, Capital orders backlog,
Backlog
c9b Capital, Depreciation, relative orders, Capital orders, Capital orders backlog, Backlog
A.12
Table C.3. Pathways originating from Supply.
Pathway no. Variable sequence
sls Supply, Acquisitions, Supply
s2s Supply, supply adjustment, desired orders, relative orders, Capital orders, Supply
sle Supply, Acquisitions, Capital
slb Supply, supply adjustment, desired orders, relative orders, Capital orders, Capital
orders backlog, Backlog
Table C.4. Pathways originating from Backlog.
Pathway no. Variable sequence
blb Backlog, desired production, capacity utilization, Production, Backlog
b2b Backlog, desired supply line, supply adjustment, desired orders, relative orders,
Capital orders, Capital orders backlog, Backlog
b3b Backlog, desired production, capacity utilization, Production, desired supply line,
supply adjustment, desired orders, relative orders, Capital orders, Capital orders
backlog, Backlog
b4b Backlog, desired production, desired capital, capital adjustment, desired orders,
relative orders, Capital orders, Capital orders backlog, Backlog
ble Backlog, Acquisitions, Capital
b2c Backlog, desired production, capacity utilization, Production, Acquisitions, Capital
bls Backlog, Acquisitions, Supply
b2s Backlog, desired production, capacity utilization, Production, Acquisitions, Supply
b3s Backlog, desired production, capacity utilization, Production, desired supply line,
supply adjustment, desired orders, relative orders, Capital orders, Supply
b4s Backlog, desired production, desired capital, capital adjustment, desired orders,
relative orders, Capital orders, Supply
b5s Backlog, desired supply line, supply adjustment, desired orders, relative orders,
Capital orders, Supply
Table C.5. Feedback loops in the Shortest Independent Loop Set of a simple Long wave model.
Loop no. Loop name Variable sequence
Li Capital decay Capital, Depreciation
L2 Supply line-1* order control Supply, Acquisitions
L3 Economic growth Production, Acquisitions, Capital, capacity
L4 Production scheduling Backlog, desired production, capacity utilization, Production
LS Capital expansion Capital orders, Supply, Acquisitions, Capital, Depreciation
Lo Backlog expansion Capital orders, Capital orders backlog, Backlog, Acquisitions, Capital,
Depreciation,
L7 Supply line adjustment Capital orders, Supply, supply adjustment, desired orders, relative
orders
L8 Order fulfillment Backlog, Acquisitions, Capital, capacity, Production
Lo Demand balancing Acquisitions, Capital, capacity, capacity utilization, Production
L10 Steady-state Capital Depreciation, relative orders, Capital orders, Supply, Acquisitions,
Capital
Lil Hoarding Capital orders backlog, Backlog, desired supply line, supply
adjustment, desired orders, relative orders, Capital orders
L12 Capital replenishment Depreciation, desired orders, relative orders, Capital orders, Supply,
Acquisitions, Capital
L13 Capital adjustment Capital, capital adjustment, desired orders, relative orders, Capital
orders, Supply, Acquisitions
Li4 Capital Self-ordering Capital orders backlog, Backlog, desired production, desired capital,
capital adjustment, desired orders, relative orders, Capital orders
LIS Steady-state Supply line Depreciation, desired supply line, supply adjustment, desired orders,
relative orders, Capital orders, Supply, Acquisitions, Capital
L16 Supply line adjustment (via
Production)
Production, desired supply line, supply adjustment, desired orders,
relative orders, Capital orders, Supply, Acquisitions, Capital, capacity
A.l4
The list of causal links and the directed cycle matrix of simple Long wave model
Tink: Variable sequence Loop no.
I) LI L2 L3 14 LS L6 L7 L8 L9 LIOLII LI2L13 L14L15 L16
cll Backlog - Acquisitions 0 0 0 0 0 0 0 0 0
cl2 Backlog - desired production
cl3 Backlog - desired Supply line
cl4 Capacity - capacity utilization
lS capacity - Production
cl6 | capacity utilization - Production
cl7 Capital - capacity
cl8 Capital - capital adj
cl9 Capital - Depreciation
cll0 | capital adjustment-desired orders
stment
0
0
0
0
0
0
0
0
1
0
clI1 | Capital orders-Capital orders Backlog | 0
cll2 Depreciation - Capital orders 0
cll3 Depreciation - desired orders 0
0
0
0
0
0
0
0
0
0
0
0
0
cll4 | Depreciation-desired Supply line
cls Depreciation - relative orders
cll6 | desired capital-capital adjustment
cli7 | desired orders - relative orders
cl18 | desired production-capacity utilization
cll9 | Desired production-desired capital
cl20 | desired Supply line-Supply adjustment
cl21 Production - Acquisitions
cl22 | Production - desired Supply line
123 relative orders - Capital orders
cl24 Supply - Acquisitions
cl25 Supply - Supply adjustment
cl26 | Supply adjustment-desired orders | 0
scoscooo oH CoO OOO oO COCO OOOH OOOH
ecoroosooosooeosoocoosoeooscoeosoooeoscooeosos
ecocooscHrocoosceoeccoeceo coco oH oH scooce
ecorocooecooscooscoocooOrKFoOoOHF oOo oOoOooOSCS
ecosceoecoeoeceoseccocooooHH oH coo Seooeco oH
rreorococoocooreooooooeoooeoooeosco
cosocosceoeoosceoscoeco ooo ooOHF OH Soo
ececocooscoroosecosecooeoscoecooHHOoOHOSCS
corr ooeoscocoooroeoooOeoH ooo oSeSoS
roorooreooFrFocooo oF oOo oOoO OOOH OS
corr ococeocooreooeoHrcooHrooocoece
corer ooeoscooHroscoeocoooOH oH Sooo SoSoSO
cooHrooorFoOFH ooo OF HF OOO OOO OF
KFOoFHF OOH OOH OOH SCO OOH OOO OCS
rOorFRrFHF ORF SCOoOFSCOCOOSoC OOOO OF OH OCS
°
Figure C.4. List of causal links and directed cycle matrix of Long Wave model.
The partial derivative equations of simple Long wave model
ac_S df.(x) ax |_|
pia xt Cx Sus SA
aC cor*B {t (rem dx = altc
8C _C¥ fay)
as cor*B
ac | SC [fas , Fa 35
0B cor B dx 0B
OS FAK) C df,(Kax § {tatoo La |
ac alte alte dk OC cor*B dx OC
as qd Co)
SS als _© Fou 2)
as alte. dk OS cor*B
aoc Ha OK $C {fa 1 fal 94}
QB alic. dk OB cor B B dx oB
OB _ folk), C Hoi Oe {Jab col Hats) 3
ac alte alte. dk OC cor cor dx 0C
aB _ _C df,() Ox
as altc dk oS
aB_ C df, (ax _C df.,(x) ax
0B altc dk OB cor dx OB
Cc cor*B C B cor C 1 Ss y filte
= + —+ —-—
(< ndl*cat cat f.,(x) C alte sat S| Cc
(CL)
(C.2)
(C3)
(C4)
(C5)
(C.6)
(C.7)
(C.8)
(C9)
The rest of the symbols are summarized in Table C.6.
Table C.6. Symbols used in Eqs. C.1-9.
Symbol Corresponding variable/parameter in the model
B Backlog
Cc Capital
Ss Supply
Sew Nonlinear function for capacity utilization
Seo Nonlinear function for relative orders
alte average lifetime of capital
cor capital-output ratio
ndl normal delivery delay
cat capital adjustment time
sat supply adjustment time
The change in the eigenvalues of simple Long wave model over time
The Eigenvalues of Simple Long Wave model
0s
04
0.3 i
0.2
0.1 erry
0 SanEEB)
0.1 fy
9215 2 3) 4 5
03 q
04 ~ H
0.5 : :
0.6 ‘
7} \ a
08 + of if
09 at
" mes
126 188 140 142 144 146 148 160 182 154 185 160 160 162 166 166 168 170.172 178 176 178 180 162
Sei — evi]
Figure C.5. Eigenvalues of simple Long wave model.
The set of equations for the causal link elasticities of simple Long wave model
eu) nen St Jes o(s ) (C10)
813 823
§23,dis
+ +
bap Oa (8:0 833.n3n t 833 nav) te, +] 2 Lb {= } (Cl
833 83 813
bn =en [fe ) _ { fos (.f0)
823 833
2te2s + Baicas 3tcab + 83108
ena = Ey {22 Jesus (821, 2s F 821,08 ) $e, * (gs: 2» + 831, ,) (C.13)
Su 821 831
; atets + Briers ste + 831.7
Sain {4 Jens (82, 1s F 821,07 ) $e * (gs Ib * 831, ») (C.14)
gu Sa 81
2 2teas t 8 r1.¢8s pas + 823,53
evs = {22 )esne(S2 Jess P (2, 2s F 821.08 ) $ey4% (2392 823,03 )
Su 13 821 $23
(C.15)
2» + 831.081 big + b
$e, (85162 831, ws) +0 (83301 S302)
831 833
yy =O (Sines Pires) +e, (Bia + 8 o1.025 + 821.075 ae)
Su 821
(C.16)
(Ss1e1 + 8s1e20 + Sse + Sica»)
+, *
831
Lug Rth {fu Je (&) C17)
821 831
cy =e Sista Waar fae face ave), {lf ne (18)
mn 21 31
+ ey) + €y3 *
833
Srtc3s S31.03b
Carig = ln * +e, *
821 831
Bac5s 831.05
cng = 1 * +e, *
821 831
S 21.065 § 31.060
Cag =n * +e, *
81 831
§ 21,098 831.0%
Cas = 1 * +e, *
821 831
8 23.b4s 833,540
C116 = Crs * +3, *
8a 833
8 5
Y Sareis DY Sani
=o. x| Ve | 22205 +| M3
Cou = a1 Hey +53
821 8x 825
8 4
Y Ssrein Y Sasa
te, *| Ft __* |+e, +e, x] V2 4
831 833
813,02 (ae + 83,535 ) (g35 bist Sassi)
cane +e); * +e; *
813 83 833
823,b4s 833.040
Carig = C3 * +e, *
823 833
A.19
(C.19)
(C.20)
(C.21)
(C.22)
(C.23)
(C.24)
(C.25)
(C.26)
(C.27)
(C.28)
e120 = €n * +23 ‘
821
et anes + 821028 23,b2s
en) =e * (Sincere Sruzze) wey e[ a Jes ¥ (801, 1st 81 ) eye Ba Jes
gu §13 Sa 83
oe (821m + S2re0) tege[ fat reve (8:12 + 8s.cs0) sone") (C31)
821 823 831 833
823 831 833
‘
fated, . [Zen] be {ees tte) ea
4 3
123 = Ca *
821 822 803
(C32)
9 4
[Feu] [Zen]
te, *| “+ —~ |+e,, +e, #|
831 833
Cetra = C1n + Cnn {fe } (C.33)
&n
Corrs = €22 (2 Jes (C.34)
8
(Be)
Baicis
Crap =e, *| VE $e, »(f24 # (Cae at fan.)
821 8x 823
(C35)
i
[Zam
4
es beg ten + Hee Sel
833
A.20
CC BY-NC-SA 4.0