A Comprehensive Eigenvalue A nalysis of System Dynamics Models
Mohamed Saleh Pal Davidsen Khaled Bayoumi
Decision Support Department SystemDynamicsGroup Information Technology Institute
Cairo University University of Bergen Cairo, Egypt
Email: m.saleh@ fci-cu.edu.eg Email: davidsen@ifi.uib.no Email: kbayoumi@ iti.net.eg
Abstract
In this paper, we develop a comprehensive eigenvalue analysis for linear models, in order to
identify the leverage points in models. The analysis is comprehensive as we develop a closed-
form analytic function relating the behavior of any state variable to all parameters in the model.
Moreover, by decomposing the behavior into several modes of behavior - each characterized by
an eigenvalue and an eigenvector - it is possible to develop a closed-form analytic function
relating a certain mode of behavior to all parameters in the model. In the first section of this
paper, we explain the mathematical foundation of eigenvalue analysis. In the second section we
identify the structural origin of the modes of behavior. This enables us to pinpoint the leverage
points of the model. Finally, in the last section, for illustration purpose, we apply the method to a
linearized version of the market growth model. The analysis of this linearized model enables us
to explain the model behavior as a superimposition of a number of behavior modes, and set the
stage for analyzing the original, non-linear version of the market growth model.
Keywords: linear model analysis, eigenvalue analysis, leverage points
I. Mathematical Background
The first step in the analysis is to express the behavior as a summation of several modes of
behavior, i.e. to decompose the behavior into several modes of behavior. In this paper, we will
decompose the time trajectory of a state variable into several modes of behavior. The time
trajectory of the state variable is a mathematical function that specifies the value of the state
variable at any time instant.
The point of departure is the structure of the model. The structure of any linear model can be
represented by the following compact matrix equation:
X¥=Gx+b (1)
Where:
X is the vector of state variables.
X is the vector of first time derivates of state variables, which is also called the slope vector, 5.
D is a constant vector.
G is the gain matrix, which is defined as follows:
ee
|
OX, ‘ ‘ Xx,
Each element, gij (Yon ), in the above matrix constitutes a compact gain; i.e. the change in the
J
net rate (slope) of state variable i, in response to an infinitesimal change in the level (value) of
state variable j, in the model.
Note that, in this paper, we denote scalars using variables in small letters; vectors in over-scored,
small letters; and matrices in capital letters.
Differentiating equation (1) with respect to time yields:
xoCr
That is:
t=Gs (2)
Where
TC is the curvature vector, which is the vector of second time derivates of state variables; i.e. ¥.
5 is the slope vector, which is the vector of first time derivates of state variables; i.e. x.
That is the gain matrix, G, relates the slope vector to the curvature vector. In other words, the
gain matrix, G, transforms the § vector into the ¢ vector in an n-dimensional standard space.
The solution to the system of differential equations, specified by equation (2), provides us with
the time trajectory of the slope vector of the model. In general, each curvature, c, is a linear
combination of all the slopes (i.e. function of the entire slope vector, s). The implication is that
the curvature, c, associated with a single coordinate (axis) in the standard space is determined by
the slopes associated with all the coordinates (component mixing). Therefore, to be able to solve
for the time trajectory of the slope, we change the coordinate system using the eigenvalue
method.
From the gain matrix, G, we can derive the eigenvalues and the right eigenvectors, as, per
definition;
GR=aqKE
If G is an n x n matrix (i.e. we have n state variables), we will have n eigenvalues, Ax, each
associated with a right eigenvector, %. The default case is to have n distinct eigenvalues, and
consequently the right eigenvectors will be linearly independent (Luenberger, 1979), and span an
n-dimensional space, R®; i.e. the right eigenvectors will form a new coordinate system, which is
known as the “eigen-space”. Note that, in general, the coordinates in this system are not orthogo-
nal, yet they are of unity length. The eigenvectors only specify a variety of directions in space
along which the dynamics of the model unfold (as explained below).
As the right eigenvectors span the n-dimensional space, then the slope vector, s, (at any point of
time) can be expressed as a linear combination of the right eigenvectors; i.e.:
Sait +027, +... +0n F, (3)
In this new coordinate system (eigen-space) the alphas (cx) will be the new components of the
slope vector.
Differentiating equation (3) with respect to time yields:
C=, +6,5, +...+0,h,
In this new coordinate system, , will be the new components of the curvature vector.
Substituting 5 -- as expressed in equation (3) -- into equation (2) yields:
t=Gla., toit+..tor |
By rearranging, we obtain:
C=aG RF +mGF, +... tonG fF,
By utilizing the fact thatG % =A, ] , we obtain:
C=aA1R +o2A27, +... tOndn I,
Andas:
C=, +655, +..+0,f,
Then, along a particular coordinate (specified by a right eigenvector), the dynamics that takes
place, can be described by the following differential equation:
Gy, = Ay Oy
The solution of the above differential equation is:
0 Ak (tz)
Ok =O €
Where t is the initial time (i.e. the starting time of the analysis period); and ox? are the initial
values of ax (at time t).
It is clear that the only factor determining the dynamics along a particular coordinate (i.e. a right
eigenvector) is the eigenvalue associated with that coordinate itself. Substituting the solution of
the dynamic behavior of each alpha (c,), into equation (3) yields the time trajectory of the slope,
ie:
tote? On + tae OF (4)
Integrating the above slope trajectory equation, with respect to time, yields:
X= (afar) (et -1} Ft. + (an an) {2 © -1} B+, (5)
Where x, is a constant vector that contains the initial values of the state variables.
The above equation is our target mathematical function that decomposes the state trajectory into
several modes of behavior, each expressed by an eigenvalue with its associated right eigenvector.
In the next section, we will interpret this equation.
II. Structural Foundation of behavior
In the first section of this paper, we developed a mathematical equation for the state trajectory
behavior (equation 5). In this section, we will identify the origin of each component of the state
trajectory. As portrayed in the figure below, the basic components are:
* The eigenvalues (1)
* The right eigenvectors ( j )
| e
+ The initial values of alphas (o,°)
>
Fig. 1: The origin of the components of the state trajectory
G X
|
x,
Note that we used delayed links, in the above figure, to indicate that the values of eigenvalues,
right eigenvectors, and alphas control the future trajectory of state variables (as indicated by
equation 5).
Note that, in our previous, work we have been focusing our efforts solely on identifying the
origin of eigenvalues, and ignoring the origin of the eigenvectors (and consequently also
ignoring the origin of the alphas). For example in Saleh & Davidsen (2001b) we identified the
marginal impact of each parameter, link, and loop, in the business cycles model, to the dominant
eigenvalue that was generating the oscillatory behavior observed, in the model. Nevertheless,
this partial analysis (i.e. by solely focusing on the dominant eigenvalue) gave us valuable
insights about the marginal impacts of each parameter, link, and loop, in the model, on the
frequency of oscillation, and on the damping ratio of the envelope within which the oscillatory
behavior is seen to unfold. However, in general, eigenvectors can play a significant role in
shaping behavior. For example, Saleh & Davidsen (2001a) demonstrated that a simple linear
model that has a single positive feedback loop can generate exponential decay behavior, rather
than exponential growth; if the initial slope vector was orthogonal to the right eigenvector
associated with the positive eigenvalue (that is responsible of generating the exponential growth
mode of behavior). For this reason, in this paper, we included the eigenvectors in our analysis.
In the rest of this section, we will explain the origin of each of the basic components of the state
trajectory.
Regarding eigenvalues, they originate solely from the structure of the model; or more
specifically, from the gain matrix, G. Recall that the gain matrix is used as a condensed
representation of the model structure, and was used, in section 1, as the point of departure to
drive the state trajectory. Also note that (as portrayed in figure 1), in linear models, the gain
matrix, G, depends on the model’s parameters (constants in the model). Hence, for each
eigenvalue, it is possible to formulate a mathematical function relating the eigenvalue (dependent
variable) to the model’s parameters (independent variables). Moreover, for simplification,
instead of formulating a single complicated function relating an eigenvalue to all parameters in
the model, it is possible to develop many mathematical functions; where each function relates the
eigenvalue to a single parameter (this is done by substituting the values of the rest of the
parameters, in the gain matrix, G). In our research, we utilized the symbolic toolbox in Matlab to
automate the above process.
Similarly, as eigenvectors also originate solely from the gain matrix, it is possible to formulate a
mathematical function relating any eigenvector (dependent variable) to any parameter
(independent variable).
The initial value of each of the alphas (i.e. ox” ) represents the projection of the initial slope
vector along a particular right eigenvector (i.e. a particular coordinate in the eigenspace); i.e. the
initial values of the alphas are dependent on the initial values of net rates (the values at the start
of the simulation), and the right eigenvectors. The initial values of net rates, in turn, are functions
of the initial state of the model (i.e.x, ), and the model parameters. The right eigenvectors, in
tum, are functions of the model parameters (as indicated above). Hence, it is possible to
formulate a mathematical function, relating the initial value of each alpha (dependent variable) to
any parameter (independent variable). Note that the initial values of sates variable are substituted
into this mathematical function.
To summarize, for each component in the state trajectory (ie. equation 5), it is possible to
formulate a mathematical function relating the component to any parameter. This enables us to
develop a compound mathematical function relating the future value, at any point of time, of a
state variable (dependent variable) to any parameter (independent variable). The partial
differentiation, at any point of time, of this compound function with respect to the parameter,
yields the sensitivity of the future value (of the state variable) to that parameter. To automate the
above process we developed a full-fledged program using the symbolic toolbox in “Matlab”
software.
III. Analyzing the Market Growth Model
IIL.I. Model Overview
The market growth model represents the nature of a single firm competing in a potentially
boundless market. The firm consists of three sectors, sales force, order fulfillment, and capacity
acquisition. Figure 2 shows the relationship between the firm three sectors and the market sector.
The initial version of the market growth model, analyzed in this paper, is taken from the CD of
the Business Dynamics book (Sterman 2000); the name of the file is “mktgrow1.sim” (Powersim
model). The model is described in details in chapter 15 in the book. We recommend the reader to
review the model, and the associated explanation, in the Business Dynamics book. However, we
made some modifications to this model (mktgrow1.sim) in order to linearize it. These
modifications are described below.
( >) TF
i
;}————__Customer
Sales Contacts
WY
Revenue Ant)
| Delivery |
(YY
Order Delay Customer
Fulfillment 483) Demand
Capacity w_|
B2 }—Orders—
BL Delivery
Utilization Delay
Capacity
Acquisition
ne:
The Firm The Market
Fig. 2: The sectors of the market growth model
IILII. Model Simplification
Here is the list of modifications which was introduced in order to linearize the model.
1. We assumed that "sales effectiveness" is a constant variable in the model, rather than an
endogenous variable whose value is determined by the current state of the model.
2. We assumed that "capacity utilization" is a constant variable in the model, rather than an
endogenous variable.
3. We assumed that the "pressure to expand capacity" is based on the "perceived backlog" (by
the company) rather than the "perceived delivery delay".
4. We linearized the table function of the "effect of pressure to expand capacity on the desired
capacity".
5. We used an additive formulation rather than a multiplicative one for "desired capacity".
For more information about this linearized model the reader may refer to the supplementary
powersim model, which is attached to this paper.
IILIII. Analysis Results
In this study, we will focus on the Backlog state variable, whose behavior is plotted in the figure
below.
25,000
20,000-4
15,0004
Backlog
10,0004
5,000+
0 20 40 60 80 100
Time
Fig. 3: The behavior of the backlog.
Using the eigenvalue approach, the backlog behavior, can be decomposed into the following
modes of behavior:
1. An exponential growth monotonic mode of behavior, shown in figure 4 (characterized by
a single eigenvalue).
2. An expanding oscillation mode of behavior, shown in figure 5 (characterized by a pair of
complex conjugate eigenvalues).
3. A decaying oscillation mode of behavior, shown in figure 6 (characterized by a pair of
complex conjugate eigenvalues).
4. Another decaying oscillation mode of behavior, shown in figure 7 (characterized by a
pair of complex conjugate eigenvalues).
Note that superimposing these four modes of behavior yields the total backlog behavior.
_ a
a 0 a
Fig. 4: The exponential growth mode associated with behavior of the backlog
/
i
\
Time
Fig. 5: The expanding oscillation mode associated with behavior of the backlog
Fist Decaying_Oscllation, Mode
Time
Fig. 6: The first decaying oscillation mode associated with behavior of the backlog
}ecaying_Oscilation_Mode
Second_D
0 20 a0 60 80 10
Time
Fig. 7: The second decaying oscillation mode associated with behavior of the backlog
Now, one can compute the sensitivities of the total backlog behavior to the model parameters; or
one can compute the sensitivities of a certain mode of behavior (e.g. the exponential growth
mode) to the model parameters. In this study, we focused on the sensitivities of the total backlog
behavior to the model parameters. However, in principle, one can easily choose a certain mode
of behavior to study rather than the total behavior. Such a deep analysis can only be performed
by using the eigenvalue approach.
As we mentioned in the previous section, the Matlab algorithm that we developed compute the
sensitivity - across time - of the backlog behavior (i.e. the sensitivity of all future values of
backlog) to a change in any parameter at the start time of the simulation. In this paper, for
demonstration purpose, we will plot the sensitivities of the total backlog behavior to the
following seven parameters in this model:
CapacityA cquisitionDelay
CompanyG oalForBacklog
SalesEffectivness
CostPerSalesRepresentative
Price
RevenueReportingDelay
SalesForceA djustmentTime
SSP Quip oh
CapacityAcquisitionDelay
a
fh 7 13 19 25 31 37 43 49 55 61 67 73 79 a5 WY 97
Sensitivity
°
Time
Fig. 8: The sensitivity - across time - of the backlog behavior to C apacityA cquisitionDelay
CompanyGoalForBacklog
Sensitivity
8
Time
Fig. 9: The sensitivity - across time - of the backlog behavior to C ompanyG oalF orBacklog
10
SalesE ffectivness
1 #7 #13 19 25 31 37 43 49 55 61 67 73 79 85 91 97
Time
Fig. 10: The sensitivity - across time - of the backlog behavior to SalesE ffectivness
CostPerSalesRepresentative
0
25 7_13 19 25 31 37 43 49 S5~6l 67 73 79 85 91 97
5
8 10
Bas
-20
Time
Fig. 11: The sensitivity - across time - of the backlog behavior to
CostPerSalesR epresentative
Price
15
10
Sensitivity
1 7 #13 19 25 31 37 43 49 55 61 67 73 79 85 91 97
Time
Fig. 12: The sensitivity - across time - of the backlog behavior to Price
11
RevenueReportingDelay
a 32 1 2 2 Z 73 79 85 91 97
Time
Fig. 13: The sensitivity - across time - of the backlog behavior to RevenueReportingDelay
SalesForceAdjustmentTime
7 #13 19 25 31 37 43 49 5) 67 73 79 85 91 97
-1000
Time
Fig. 14: The sensitivity - across time - of the backlog behavior to
SalesF orceA djustmentTime
12
References
Abdel-Aleem, B. (2004): An automated system to analyze system dynamics models -- Case
Study: Commodity Production Cycles. M. Phil. Thesis. University of Bergen. Norway.
Abdel-Gawad, A. (2004): An automated system to analyze system dynamics models. M. Phil.
Thesis. University of Bergen. Norway.
Barlas, Y.; Kanar K. (2000): Structure-oriented Behvior Tests in Model Validation. Proceedings
of the 2000 Intemational System Dynamics Conference. Bergen.
Davidsen, P. (1991): The Structure-Behavior Graph. The System Dynamics Group, MIT.
Cambridge.
Eberlein, R. (1984): Simplifying Dynamic Models by Retaining Selected Behavior Modes. Ph.D.
Thesis, M.I.T., Cambridge, MA.
Forrester. N. (1982): A Dynamic Synthesis of Basic Macroeconomic Policy: Implications for
Stabilization Policy Analysis. Ph.D. Thesis, M.I.T., Cambridge, MA.
Forrester. N. (1983): Eigenvalue Analysis of Dominant Feedback Loops. The 1983 International
System Dynamics Conference, Plenary Session Papers, pp. 178-202.
Ford, D. (1999): A Behavioral Approach to Feedback Loop Dominance Analysis. System
Dynamics Review. Volume 15, Issue 1, pp. 3-36.
Goncalves, P.; Lertpattaraong, C.; Hines, J. (2000): Implementing Formal Model Analysis.
Proceedings of the 2000 International System Dynamics Conference. Bergen.
Kampmann, C. (1996): Feedback Loop Gains and System Behavior. Proceedings of the 1996
Intemational System Dynamics Conference. Boston.
Kreyszig, E. (1979): Advanced Engineering Mathematics. Fourth Edition. John Wiley & Sons.
Luenberger, D. (1979): Introduction to Dynamic Systems: Theory, Models and Applications.
John Wiley & Sons.
Mojtahedzadeh, M. (1996): A Path Taken: Computer-Assisted Heuristics for Understanding
Dynamic Systems. Ph.D. Thesis. Rockefeller College of Public A ffairs and Policy. Albany NY.
Myrtveit, M.; Saleh, M. (2000): Superimposing Dynamic Behavior on Causal Loop Diagram of
System Dynamics Models. Proceedings of the 2000 Intemational System Dynamics Conference.
Bergen.
Ogata, K. (1997): Modem Control Engineering. Third Edition. Prentice-Hall.
13
Press, W.; Teukolsky, S.; Vetterling, W.; Flannery, B. (1992): Numerical Recipes in C. Second
Edition. Cambridge Univ. Press.
Reinschke, K. (1988): Multivariable Control: A Graph-theoretical Approach. Lecture Notes in
Control and Information Sciences. Springer-V erlag.
Richardson, G. (1984): Loop polarity, loop dominance, and the concept of dominant polarity.
System Dynamics Review Vol. 11; pp. 67-88.
Saleh, M.; Davidsen P. (2001a): The origins of behavior pattems. Proceedings of the 2001
Intemational System Dynamics Conference. Atlanta, USA.
Saleh, M.; Davidsen P. (2001b): The origins of business cycles. Proceedings of the 2001
Intemational System Dynamics Conference. Atlanta, USA.
Sterman, J. (2000): Business Dynamics: Systems Thinking and Modeling for a Complex World.
McGraw-Hill.
14
