To Main Procee
igs Docume!
Implementing formal model analysis
Paulo Goncalves, Chalermmon Lertpattarapong’’, and Jim Hines”
Abstract:
There have been a number of important contributions to the field of system dynamics in the
recent past. One of the most exciting has been the development of eigenvalue elasticity theory
(Forrester 1982, Kampmann 1996), which promises to increase our capacity to understand our
models, enhancing our ability to get insight from modeling. In this paper, we critically review the
current state of knowledge and provide an improved approach to determining loop gain
elasticities.
1, Introduction
There have been a number of important contributions to the field of system dynamics in
the recent past. A significant one was the introduction of advanced software
environments that made use of modern graphical computers. IThink/STELLA, Vensim,
and Powersim link diagrams to equations, simplifying considerably the process of
building models. These software environments have more recently permitted new ways
of running models by making it possible host models on the world wide web and by
making it possible to create interactive, multimedia management games. The user
friendly software environments have made modeling more accessible, and the number of
new practitioners in system dynamics has consequently multiplied.
Although the improved software has lowered the barriers to entry to the field and
reduced the cost of creating and running models, analysis of models remains difficult for
the inexperienced and time-consuming for everyone. The creators of Vensim made a
valuable contribution by automating and illuminating the process of causal tracing. But,
beyond this there has been essentially no advance in analytical tools since 1960, when
non-overlapping data plots were introduced in DY NAMO.’ The analysis of models in a
rigorous, quick, reliable, and standard manner is still today one of the most important and
challenging problems in system dynamics.
A research area that promises to increase our capacity to analyze our models is the
idea of eigenvalue elasticity.’ (N. Forrester 1982, Kampmann 1996) The process
essentially is to linearize the model at a point in time, calculate the eigenvalues and then
consider how the eigenvalues change as link gains change for the linearized system. The
resulting understanding of how link-gains impact eigenvalues can be used directly in
deciding which links to strengthen or weaken, or indirectly to figure out which loops
contribute most heavily to the eigenvalues. Software capable of carrying out the process
was developed some fifteen years ago in an experimental version of DY NAMO.
Unfortunately, too few people used it, and eventually Pugh- Roberts discontinued the
experimental version without ever determining how useful the approach would be on
real-world problems.
* Ph.D. Student, System Dynamics Group, Massachusetts Institute of Technology, 30 Wadsworth St.,
E53-358A, MA 02139, USA.
“M.Sc. Student, Leaders for Manufacturing, Massachusetts Institute of Technology.
“ Senior Lecturer, System Dynamics Group, Massachusetts Institute of Technology.
‘ The new software environments, especially Vensim, have also increased the ease of incorporating data
into models and have provided new tools for model calibration and policy optimization.
* Another promising area for formal model analysis is pathway participation (Mojtahedzadeh 1996), but we
have not reviewed this technique here.
Forrester showed that a complete description of link elasticities allows one in
principle to calculate loop elasticities. This suggestion though never implemented in
software, promised to provide an answer to how model structure, that is a set of feedback
loops, determines model behavior.
The particular calculation that Forrester suggested is actually not feasible. As he
realized later, Forrester’ s suggested approach results in a system of equations that is over-
determined - an effect of the fact that the number of loops increases much faster than the
number links. Kampmann discovered that a small subset of loops is sufficient to
uniquely describe eigenvalues (i.e. the behavior) of a system dynamics model.
(Kampmann 1996) Using such an “independent loop set” produces a smaller system of
equations, a system that can be solved.
A problem remains even with Kampmann’s approach. Usually, there will be a
number of different independent loop sets. While any one of them can be used to
discover how eigenvalues change as members of the set change, not all loops will be
members of the set. The question remains which loopset to choose, or whether an
arbitrarily chosen independent loopset can provide intuitively- appealing insight into a
system.
In this paper we present an alternative way to compute loop-gain eigenvalue
elasticities. This approach provides information on every loop in the model. We sidestep
the problem of an over-determined set of equations by eliminating the intermediate step
of using the relatively small number of link elasticities to get the relatively large number
of loop gains. Instead we express the characteristic equation directly in terms of loop
gains; and then use this equation to go directly to loop gain elasticities. Our approach
calculates a unique gain elasticity for every loop in the system, not just for a subset of
loops.
2. Eigenvalue elasticity theory
The formal structure of a system dynamics model with a vector of state variables x(t),
where x(t) =(x,, X,, ..., X,)’, a vector of auxiliary variables y(t), where y(t) =(y,, y,, ..., Vp
)’, avector of nonlinear functions f(t), where f(t) =(f,, £,, ..., f,)’ and g(t), where g(t) =
(Gus Gyr vee g, )’, can be represented by the set of nonlinear differential equations:
(Kampmann 1996)
X(t) =£(x(t), y(0), t), (1)
y(t) =g(x(t), y(t), t).
When the system is not influenced by exogenous variables, f and g are not
functions of time. So, we can simplify the notation.
¥ =f(x,y),
Bea. (2)
y =g(x,y).
Given the complexity of the system of nonlinear differential equations, we gain
tractability by linearizing it around an operating point x and y where, Y= (X,Y).
Linearization allows every variable (v,) to be expressed as a linear combination of other
variables (v, where j = 1, 2,..., i,...,n) in the model, such that:
wy,
v, “2 a,v,, where, ay ve (3)
and, a, is the partial derivative of variable v, with respect to variable v,.
We will illustrate with a version of the familiar workforce inventory model. The
model captures a simple production system. The model attempts to maintain desired
inventory by adjusting production via hiring and firing workers. More precisely:
Inventory integrates the difference between production and shipments. Shipments are
determined by demand reduced by stock-outs, should inventory fall too low. Production
depends on the workforce. And the workforce is “anchored” to the level necessary to
meet expected demand. The workforce is increased above this anchor if inventory is too
low and conversely workforce is decreased below the anchor if inventory is too high.
Expected demand is a smooth of actual demand.
A stock and flow diagram of the model is shown below. The model is composed
of three state variables, four flows, three auxiliary variables, two exogenous variables,
and five constants.
Minimun Sales
Time (MST)
| Inventory 7;
1
Producing id Stes
(P) (Ss) . Demand
. Desired (D)
Correction Inventory
Time (CT) Di)
Inventory
connen
i E ‘ted
sara Demand (ED)
Change in
E ter
Workforce Desired Demand
w) Progusng (CED)
Hiring/Firin DP
Rate HFR) > a
Time to Change
D id A
workforce E pen s
Hire/FireTime (DW)
(HFT)
Figure 1 - Diagram of a linear system dynamics model.
I=P—-S=PDY-W —-D IC =(DI -1)/CT
W =HFR =(DW -W)/HFT DP =IC +ED
ED =CED =(D -ED)/TCE DW =DP/PDY
In terms of equation (2) f(x, y) =(f,, f,, f, )’, where, for example, f,(I,W,ED, IC,
DP, DW) = PDY.W- D andg&x, y) =(g,, 9,, 9, )’, where g,(I,W, ED, IC, DP, DW) = (DI
- D/CT.
Such a system of differential equations can be represented it in matrix form as:
Yy AB b
: x “ay
where the tilde represents the deviation from the operating point, the entries in matrices
A,B, C, and D are given by the partial derivatives, A,,, =on/ ox, Bi xp =a/ ay,
Cyyn =y/0K, Dy,» =0y/ dy, and b=f@,y).”
The system of differential equations in figure 1 in matrix form becomes:
I 0. PDY 0 I 00 0 4
W =0 4/HFT 0 W +0 0 1/HFT DP + 0 D
ED 0 0 -Y/TCE ED 00 O DW 1/TCE
Ic -A/cT 00 1 0 0 0 1 1/CT
D= 0 O01W+4#+1 0 O DP + 0 DI
DW 0 00 ED 0 1/PDY 0 DW 0
Or, simply:
i 0 PDY 0 0 0 0 1 4 0
w 0 A/HFT 0 0 0 1/HFT W 0 0
ED 0 0 <-/TCE 0 0 0 ED 1/TCE 0 D
Ic ct 0 0 0 0 0 i * 0 wer DI
DP 0 0 1 1 0 0 oP 0 0
DW 0 0 0 0 1/PDY 0 DW 0 0
0 PDY 0 00 0 4/CT 0 0 0 0 0
A=0 -/HFT 0 ,B=0 0 1/HFT,C= 0 01,D=1 0 0
0 0 -/TCE 00 oO 0 00 0 1/PDY 0
* Note that when we have models with exogenous constants, we can still have functions f and g
independent of time, nevertheless they can modify the differential and algebraic equations. In such cases,
the exogenous constants are gathered in a new vector (c and d) and a matrix (E and F) multiplies it to
capture the relationships they have with the system. Under this condition equation (4) becomes:
Y ABE b Ec
yc py*to*rFa
We can also represent the system in a more compact diagrammatic way using graph
theoretical notation. Thus, the system in figure one becomes:
Figure 2 - A graph theoretical representation of a system dynamics model in non-reduced
form.
Additionally, Any system can be represented in reduced form by eliminating y
from equation (4):
X=Ax +By +b,
y =Cx+Dy
(1-D)y =Cx
Assuming that (1-D) is invertible we can write
y =(1-D)"Cx
Substituting into the equation for X we have:
X= Ax +B[(1-D)'Cx]+b
And, rewriting the equation above, we obtain the result in reduced form:
¥ =Jx +b,
(5)
J =A +B(1-Dy'C.
where, the matrix J =0¥/ox is known as the Jacobian of the system.
Returning to our example, the Jacobian (J) leads to the following relation:
0 PDY 0
J = -/HFT-PDY-CT -1/HFT YHFT-PDY ,
0 0 A/TCE
The eigenvalues (A) of the matrix J describe the behavior modes inherent in the
model. The eigenvalues are the solutions of the characteristic polynomial P(A):
P(A) =|Al, —J| =0 (6)
The characteristic equation for the example provided is given by:
A -PDY 0
P(A) =|Al, -J|= YHFT-PDY -CT A+1/HFT A/HFT-PDY ,
0 0 A+/TCE
=(A +1/TCE)(¥ +A/HFT +1/(HFT -CT))
1 1
=F + + x +
ce Tarr” *
1 1 1
+ A+
TCE- HFT Ch-HFT TCE-HFT-CT
And, the eigenvalues, for the example, in terms of the link gains are:
A, =——
TCE
,-—! or
2 2-HFT 2VHFT? CT -HFT
jae a
3 2-HFT 2VHFT? CT-HFT
Since the Jacobian (J) is computed from A, B, C, and D and the entries of those
matrices are the partial derivatives or the link gains (a,) in the system dynamics model, it
is simple to obtain the characteristic polynomial P (A) in terms of link gains (a). Thus,
Forrester suggested measuring the sensitivity of an eigenvalue with respect to a specific
link by simply computing the partial derivative of the eigenvalue with respect to the link
gain. This would allow one to understand how the strength of a link could impact specific
modes of behavior.
Sy =", (7)
Additionally, one could normalize the sensitivity measure to isolate the effect of
the change from the sizes (values) of eigenvalues and link gains. This normalization
could be obtained multiplying the sensitivity by the ratio of the link gain to the
eigenvalue. He defined this measure eigenvalue elasticity with respect to link gain or link
gain elasticity.
OA, ay
4 =3F 8
Ey 2a, A, (8)
“This can be directly observed from the example above.
Although it is not common practice to write the characteristic polynomial P (A) in
terms of the loop gains (g,) and Forrester never explicitly wrote P (A) in terms of the loop
gains, he did suggest that he could obtain the loop gain elasticities from Mason’s rule.
Thus, he extended the results for link sensitivity and link elasticity to loops.
=. and Ey i (9)
® in Mk
Considering that loop gains g, are determined by the product of the gains of links
in the loop and that, according to Mason’s rule, the characteristic equation can be written
in terms of the loops gains, we obtain that the eigenvalues are a function of link gains:
A, =hG,@a,)) (10)
Given also that any link can be a part of several distinct loops, Forrester rewrote
Ey:
dA, oy
E, = se mA 4 (11)
4 over all n oY, 0a
The partial derivative of ihe wi gain (g,) with respect to the gain of a link (a,)
not in the loop is zero (dg,/0a, =0). And, the partial derivative of the loop gain with
respect to the gain of a link (a,) in is loop is the ratio of g, by a, (09,/da, =9,/a, ).
Then, we can obtain a final equation for the link gain elasticity:
oA, Io Bi a. 9,20
E,
Ata, By @) Ay Sta, Dy Ak
inn notinn
A, 9,
Rare Ha
where A, is the set of all loops containing the link from v, to v,.
‘quation (12) shows that the link gain elasticity i ig equal to the sum of the loop
gain elasticities of all the loops composed by link a,. (N. Forrester, 1983) At this point,
the theory of eigenvalue elasticity could tell how link elasticities could influence
behavior, by its impact on the eigenvalues. Ultimately, it would be interesting for
practitioners to be able to say how structure drives behavior. Thus, it would be a stronger
result to say tell how specific loops would influence behavior. While Forrester hinted at
the problem, he didn’t provide a clear the solution.
In 1996, Kampmann introduced a few developments to the theory that had the
potential to answer how specific loops could drive model behavior. First, Kampmann
formalized the existing notation in the field through the introduction of graph theory. This
allowed him to describe the topology of a system dynamics model in a simpler way using
directed graphs.
Then, he explored the interdependencies between the number of loops and links in a
maximally connected system. He observed that the number of loops increased much
faster than the number of links even in modestly connected systems. Additionally, he
knew that a single link could be part of many different loops. So, he suggested that
considering loop gain elasticities would only make sense when it was possible to change
one loop without influencing other ones. In his own words he asks:
“Under what circumstances can a set of loops have their strengths determined
or changed independently by an appropriate assignment or change in the
strengths of individual links in the system?” (Kampmann 1996)
He finds a solution to that question when he analyzes the system of algebraic
equations relating loop gains (g,) and link gains (a, =e).
q
9. =[] gle) (14)
N q
Ing, =) inlate } (15)
Inlg =Clnld (16)
where equation (14) shows that the gain of loop n is just the product of link gains
composing the loop; g is the vector of loops gains, where g = (4, g,,.... 9)’ with L
number of loops in the system; e is the vector of link gains, where e =(e,, €,, ..., y)’ with
N number of links in the system; and, C isa L x N matrix.
Since there are more loops than links (L > N), the system above is over-
determined. The rank p of the matrix is lower than L and often the system does not have
a solution. A solution only exists, when the left-hand side is a linear combination of the
rows of C. In this context, the loop gains were not independent from one another. Thus,
to achieve independence it was necessary to obtain a set of loops that were linearly
independent. Kampmann called this set the independent loop set. And, extending an
important result from graph theory, he showed that in a strongly connected graph (G),
with N arcs and n nodes, the number of loops in the independent loop set was N-n+1.
The impact of the independent loop set result could be potentially high, because it
permits us to calculate loop elasticities from link elasticities. Unfortunately anyone using
this approach must decide which independent loopset to use. The problem is that any one
independent loopset will exclude certain loops (the dependent loops) from the analysis.
At this point we do not know how (or whether it is possible) to use an independent
loopset to calculate the elasticities of the corresponding dependent loops. Without being
able to perform such a calculation, its uncertain whether practitioners can understand the
information coming from Kampmann’s approach in terms of the loops that are most
intuitively understandable.
A way out of these difficulties is to sidestep the problems inherent in going from
link elasticities. Kampmann (1996) provided a general formula for expressing the
characteristic polynomial P(A) in terms of the loop gains of a reduced or non-reduced
system on a theorem, that is a direct consequence of Mason’s rule. His theorem stated
that:
“Theorem: Let J be the n x n Jacobian matrix for a system with graph
representation G, which may be either reduced or non-reduced, and let P(A)
be the corresponding characteristic polynomial
P(A) =det(Al, —J) =A +> p, A"
i=}
Then p, => (MoD, )
where the summation is taken over all cycle compositios D, of order i.”
Before we provide an example illustrating the application of the theorem we
introduce some notation. The order of a cycle is the number of state variables it contains.
A cycle composition (D) is a collection of node-disjoint cycles that covers all nodes in G.
the gain of a cycle composition (g(.)) is given by the product of the link gains forming the
cycle composition. The order of a cycle composition (i) is the sum of the order of its
constituent cycles. Dj represents all cycle compositions of order i. The size (|Dj |) of a
cycle composition is the number of cycles in D.
Example 2 shows a graph theoretic diagram of a second order system in reduced
form. According to the theorem above, the characteristic polynomial will have second
order (n=2). The general formula for the characteristic equation of the system will be
P(A) =A +p,-A +p). The term in p, will be the sum of all cycles compositions (D) of
order one (i=1). This corresponds to the sum of the gains of all first order loops in the
diagram below. Additionally, the size of this cycle composition will be equal to one (|D3|
=1). Thus, the signs of the loop gains will be negative. The term in p, will be the sum of
all cycles compositions (D) of order two (i=2). This corresponds to the sum of the gains
of single second-order loops (loops involving both two state variables) and the product of
gains of all pairs of disjoint first-order loops. The size of the former type of cycles will be
equal to one (|D2| =1) and have a negative sign. The size of the latter type of cycles will
be equal to two (|D2| =2) and have a positive sign.
a2l
all CE 2 a22
al2
8” (aj; +Ay))*S +{ Ay Ay +ay; “Ayy)
s’ —(G, +9,)8 HG, +0; G2)
The theorem allows us obtain an implicit characterization of the loop gain elasticity.
Through the application to several system dynamics models we hope to gain further
insight about how the loop structure drives model behavior. Now, we apply the theorem
to example one. This will illustrate the results of the theorem to a system in non-reduced
form. Table 1 shows the cycle compositions, order, gains, sizes, and signs for
determining the coefficients of the characteristic polynomial for example 1.
Table 1 - Cycle compositions, cycle order, cycle gains, sizes, and signs for determining
the coefficients of the characteristic polynomial.
Cycle Composition Order Gain Size | Sign
@) @ Gg) (Di) | CD?!
__ 1 g,=-1/HFT 1 “1
Ww
-1/HFT
-1/TCE 1 g,=-1/TCE 1 1
ED
2 g,=-1/(CT*HFT) | 1 ol,
2 2 +1
9, T=
1/(TCE*HFT)
10
Table 1 - Cycle compositions, cycle order, cycle gains, sizes, and signs for determining
the coefficients of the characteristic polynomial.(continued)
Cycle Composition Order Gain Size | Sign
@) @ @() (Dip | GD
3 G) G= 2 +1
1/(CT*HFT*TCE)
-1/TCE
“p
From the entries in the table above and the general formula in the theorem we obtain the
terms p, , p, and p,: 5
p, => (=P g(D,),
D,
P, =- Gyr
P, =+ 9,9, 9;
P3 = + 9.93
Now, it is possible to obtain the same result from the Jacobian on the reduced
system. We simply compare the solution for the characteristic polynomial with the
general form of the polynomial, to obtain p,, p,, and p,. Where the results should be equal
to the one derived above.
P(A) =A’ +p, -¥ +p, -A+p,
P(A) =” +9, ~9,)4 +(Gige —Gs)A +92
where,
ee
‘HFT TCE
1 1
P) “HET-TCE HFT -CT
_ 1
Ps “IFT -CT -TCE
And, the eigenvalues, for the example, in terms of the loop gains are:
A, ie
A, =i tvh —4g,
1
A; =“ =) gy —4g,
11
Obtaining the characteristic polynomial P(A) in terms of the loop gains (g,) allows
us to implicitly differentiate the equation P(A; g,) =0 with respect to the loop gain (g,).
From which we obtain:
dP (A;g,) _OP(A;g,) dA , oP(A;g,)
=0 17
dg, OA dg, Dy, mm
m 7 a
GA __ Pig) PAG,) (a8)
dg, Con A)
4
AP(A;g,) OP(A;g,) Ax
Ei es ee 19
* , A 9)
3. Discussion
Forrester suggested deriving loop elasticities from link elasticities. He described
how to create a system of equations involving link elasticities that in principle could be
solved for loop elasticities. Unfortunately, the large number of loops relative to links
meant that in practice the system will usually be over-determined. Kampmann solved
this problem by introducing the idea of an independent loopset. Unfortunately there will
usually not be a single independent loopset and Kampmann’ s approach leaves us with the
problem of choosing one independent loopset from many.
We suggest sidestepping this problem by going directly to loop elasticities from
Mason’s rule and the characteristic equation, rather than taking the intermediate step
through link elasticities. Kampmann noted that it would be possible to express the
characteristic polynomial P(A) in terms of the loop gains (g,) and obtain the loop gain
elasticities without addressing the issue of an independent loop set. Kampmann, however,
found this approach to fall short of his own standards:
“,, the formalism [of obtaining the loop gain elasticities directly from the
characteristic polynomial] is only valid if it is in fact possible, by an
appropriate change to the gains of individual links in the system, to change
the gains of one loop independently of all others.”
We think that Kampmann set the bar a bit too high, putting forward a standard
that no approach can meet. Even an independent loopset can contain a loop whose gain
can not be changed independently of the gains of the other loops; because even an
independent loopset can contain a loop composed only of links shared with other loops in
the independent loopset. Fortunately, the standard that Kampmann suggests is higher
than necessary: In practice people routinely alter loop gains, even though doing so
affects the gains of other loops. In fact it is unlikely that we would ever find a loop in the
real world whose gain could vary without affecting the gains of some other loops.
The only relevant standard to judge the usefulness of an approach is whether it
does in fact help us to improve the system by making structural (including parametric)
changes. Such improvement requires an understanding of which structures in a system
are most responsible for creating the behavior of concern. There is no requirement that
we must find structural sources of behavior that are completely disconnected from the
12
rest of the system. The yardstick by which we ought to measure any analytical approach
is whether it helps us in the task of understanding the connection between structure and
behavior. Examining the elasticity of eigenvalues (i.e. behavior modes) with respect to
loops promises to help us achieve this understand more clearly, more rigorously, and
more quickly than is possible using traditional system dynamics approaches. Finding the
elasticities directly, rather than indirectly, seems to us to be the best approach to finding
the loop gain elasticities.
4. References
Eberlein, R.L., “Simplification and Understanding of Models.” System Dynamics Review,
1989. 5(1).
Forrester, N. A Dynamic Synthesis of Basic Macroeconomic Policy: Implications for
Stabilization Policy Analysis. Ph.D. Thesis, M.I.T.,1982. Cambridge, MA.
Forrester, N., “Eigenvalue Analysis of Dominant Feedback Analysis.” Proceeding of the
1983 International System Dynamics Conference, Plenary Session Papers. System
Dynamics Society: Albany, NY . pp178-202.
Kampmann, C.E., “Feedback Loop Gains and System Behavior.” In: Richardson GP and
Sterman JD (Eds.) Proceeding of the 1996 International System Dynamics Conference
Boston. System Dynamics Society: Albany, NY. pp. 260-263.
Mojtahedzadeh, M.T. A Path Taken: Computer-Assisted Heuristics for Understanding
Dynamic Systems. Unpublished Ph.D. Dissertation, Rockefeller College of Public A ffairs
and Policy, State University of New Y ork at Albany, 1996, Albany NY.
Reinschke, KJ. Multivariate Control: A Graph Theoretical Approach. Lecture Notes in
Control and Information Sciences. 1988. Berlin: Springer-Verlag.
Richardson, G.P., “Loop Polarity, Loop Dominance, and the Concept of Dominant
Polarity.” In Proceeding of the 1984 International System Dynamics Conference. Oslo,
Norway.
Richardson, G.P., “Dominant Structure.” System Dynamics Review, 1986. 2(1): 68-75.
5. Appendix
A graph (G) is a mathematical object consisting of points or nodes(R) and links or edges
(N). We can assign properties to the elements above to expand its applications. For
instance, it is possible to assign weights or gains and directions to links. A weighted
graph is formed using the former assignments and a directed graph or digraph is formed
using the latter. A directed edge is called an arc. A series of disjoint nodes (1;) connected
by arcs, > >...>1q_ is called a directed path. Two nodes 1 and ry are said to be
strongly connected when there is a path from 7 to rm and from rx to 4. A digraph is said
to be strongly connected if for every node 1; and r, there is a path from 7, to r,. When the
initial node (r,) and the terminal node (rq ) in a directed path are the same (1, =rq) we
13
obtain a closed path. A closed path is called a directed cycle when going along the path
one reaches the initial node only once. An arc connecting a node to itself (1, > 4; ) is
called a self-loop.
In a graph representation of a system dynamics model, nodes are state variables
x(t) or auxiliary variables y(t) and edges are the connecting links. Arcs determine the
causality assignments and link capacities provide link gains. A directed cycle maps a
loop, and a self-loop is also a minor loop.
14
