Extending eigenvalue analysis to nonlinear models via
incorporating higher order terms of Taylor series expansion
Mohamed Saleh Pal Davidsen
Decision Support Department System Dynamics Group
Cairo University University of Bergen
Saleh@SalehSite.info davidsen@ifi.uib.no
Abstract
This paper is a concept paper about a suggestion proposed by Nathan Forrester, in the last
year conference to extend eigenvalue analysis to nonlinear models. His idea was to consider
higher order terms of the Taylor series expansion when approximating nonlinear models. In
this paper, we demonstrate the feasibility of Nathan's idea. The main contribution of this
paper is to devise a pragmatic approach to solve the resulting equations of Taylor series
expansion. This pragmatic approach is based on our novel concept of 'smoothed Jacobian'
matrix, which is computed from both the ordinary Jacobian matrix and the set of Hessian
matrices. Recall that the elements of the ordinary Jacobian matrix represent slopes of
relationships, while the elements of the Hessian matrices represent curvatures of
relationships. So by integrating the elements the ordinary Jacobian with the elements of the
Hessian matrices, we are actually smoothing the slopes given the knowledge about
curvatures. Consequently we are smoothing the time trajectories of eigenvalues and
eigenvectors in nonlinear models.
Keywords: system dynamics, nonlinear model analysis, eigenvalue analysis, Taylor series
expansion, behavior modes.
1. Mathematical Background
Any nonlinear system dynamics model consists of a set of nonlinear functions between state
variables and net rates. In approximating any single nonlinear function in the model: If we
only use the first order term of Taylor series expansion — which is the current practice — we
have to compute the gradient vector associated with this function. However, if we also used
the second order term of Taylor series, we have to compute the Hessian matrix associated
with this function. In general, in a model with 'n' state variables, one has to compute 'n'
Hessian matrices — in addition to computing the Jacobian matrix (where each row of the
Jacobian is, in fact, a gradient vector).
Note that in many nonlinear functions, the higher order terms (i.e. the third and higher) will
equal zero. This assumption is based on the following two reasons:
1. Many table functions can be approximated by a second order polynomial (so
differentiating twice is enough).
2. A major source of nonlinearity is the multiplication of two state variables in a rate
equation (so also in this case differentiating twice is enough).
So, based on the above assumption, one may conclude that in a certain category of nonlinear
system dynamics models, all structural information of the model can be captured in the
Jacobian and the 'n' Hessian matrices.
h
The following equation demonstrates Taylor series expansion associated with the i" net rate.
X(i) =X, (i) + grad,"(x —x,)+0.5(x—x,)" H,(x-x,)
where
X(i) is the value of net rate 'i' at time t (t represents any point of time during the simulation)
xX is the state variables vector at time t.
X, is the state variables vector at a specified anchor point, in the n-dimensional space.
X,(i) is the value of net rate 'i' computed at Xo
grad," is the transpose of the gradient associated with net rate 'i' computed at Xo
H, is the Hessian matrix associated with net rate 'i' computed at Xy
Note that in this paper, we denote matrices by bold and capital letters, vectors by bold and
small letters, and scalar values by small letters.
In the above equation, taking X, as the equilibrium point, and denoting x-X, by vector y
yields:
YG) =grad,"y+0.5y" Hy
Note that net rates equal zero at the equilibrium point, and also note that y(i) = X(i)
Basically, vector y represents the deviation vector from the equilibrium point. It is important
here to stress that the above equation is valid at any point in the n-dimensional space — i.e. not
only in the neighborhood of the equilibrium point. Moreover, it does not matter if the model
actually reaches equilibrium, or never reaches it.
The main contribution of the paper is to devise a pragmatic approach to solve the above
equation. This is done via putting the above equation in the following compact form:
7S
y-J’y a
where
y is the net rate vector at any point of time
y is the deviation vector at any point of time
J$ is a nxn matrix. The i" row of this matrix, J‘(i,:), is defined as follows
Ji.) = grad," + 0.5y" H,
The beauty of equation | is that it is not an approximation, but rather an exact equation at any
point. Consequently, it is valid at any point in time.
Each element of matrix J‘ is either a zero, or a constant, or a linear function of the state
variables. In this paper, we call matrix J* the 'smoothed J acobian'.
Using a Symbolic Toolbox (e.g. Matlab Symbolic Toolbox) one can symbolically compute
the smoothed eigenvalues’ and smoothed eigenvectors associated with matrix J. Based on the
eigenvalues and eigenvectors, one can solve the time trajectory of vector y, as demonstrated
below:
Y = Wi exp(Ait) + Wo exp(Agt) +...+ Wn exp(Ant)
where:
Wi), W2, ... Wn are the weights vectors
Au, Aa, «+++ An are the smoothed eigenvalues
' The eigenvalues are called smoothed as they are derived from the smoothed Jacobian.
For more information about how to compute the weight vectors, the reader may refer to Saleh
et al. 2006.
Based on the above equation, the state variables vector trajectory can be expressed as:
X= Wy exp(Ait) + Wo exp(Agt) +...+ Wn exp(Ant) + Xo (2)
Recall that x, is the state variables vector at the equilibrium point.
Equation 2 enables us to decompose the behavior of state variables into several modes of
behavior. By plotting the smoothed eigenvalues across time, it is possible to detect
bifurcation points in behavior modes, and hence to chop the simulation time into distinct
phases. This process will be demonstrated on a simple model in the next section.
2. Case Study
In this paper, we will apply the above concept on a simple epidemic model described in
chapter 9, in Sterman 's book (Sterman, 2000). The model is called the "SIR" model as it
consists of the following three stocks: Susceptible Population (S), Infectious Population (I)
and Recovered Population (R). The following figure shows the stock and flow diagram of the
model, and the model itself is attached as a supplementary material.
Initial
Infectious
Population
Susceptible | ________sy_____gy Infectious |__| Recovered
Population | Population R
Population §
GS Infection 4r) GC} Recovery
Rate Rate
Depletion | |. Contagion Recovery
yy, A
Contact Rate c Kivétiigs
Total . Duration of
Population N intectiviry 4 liness d
Figure 1: Stock and flow diagram of the SIR model
As shown, in the above figure, there are the following 3 loops in the model:
1. The Balancing "Depletion" Loop
2. The Balancing "Recovery" Loop
3. The Reinforcing "Contagion" Loop.
By removing auxiliary variables, the model can be put in the following condensed form (the
reader may check attached model):
dS/dt = -0.00015*I*S;
di/dt = (0.00015*I*S)-(0.5*1);
dR/dt = 0.5*1;
From the above equation, it is clear that the source of nonlinearity in rate equations is the
multiplication of two state variables; hence as stated before, all structural information in the
model can be captured by the Jacobian matrix and the '3' Hessian matrices, which are as
follows:
—0.00015*I —0.00015*S 0
J =| 0.00015*I 0.00015*S-0.5 0
0 0.5 0
0 —0.00015 0
H, =| — 0.00015 0 0
0 0 0
0 0.00015 0
H, =| 0.00015 0 0
0 0 0
Note that all elements of the '3' Hessian matrices are either constants or zeros. This is another
proof that higher order terms in Taylor series expansion (i.e. third order and above) are zeros.
In our analysis, we selected the point S=0, I=0 and R =0 as our equilibrium point X,. Note
that it does not matter that the model will never reach this equilibrium point. In fact, the
model never approaches the neighborhood around the origin point; this is because the
summation of the '3' stocks is always constant, and equals the total population. In general,
one can select any other equilibrium point and carry out the same process.
In this case vector X equals vector y; and equation | can take the following form:
x=J*x
To compute matrix J* one has to first compute matrix J at x,. Recall that the gradient vectors
correspond to rows in the J matrix computed at Xo. Matrix J computed at X, is as follows:
0 0 0
J),,=|9 -0.5 0
0 05 0
And hence matrix J* matrix is as follows (see equation 1):
—0.000075*I —0.000075*S 0
J* =| 0.000075*I 0.000075*S-0.5 0
0 0.5 0
The first row of J* is computed from the following formula: grads" + 0.5 x’ Hy
The second row of J* is computed from the following formula: grad)" +0.5x' Hy
And the third row of J* is computed from the following formula: grada! +0.5x' Hp
Note that the "smoothed Jacobian" matrix J‘ is different than the "ordinary Jacobian" matrix
J. For example, the first element of the first row in matrix J* is equal -0.000075*I; while the
corresponding element in matrix J is equal -0.00015*I. This is an important issue. Recall that
the normal practice is to focus only on the ordinary matrix, and to ignore the 2" order term of
the Taylor series expansion. This can be considered an error term. And as the width of the
analysis interval increases, this error term can be significant. For this reason, we postulate
that in nonlinear models, it is preferred to use the smoothed Jacobian instead of the ordinary
Jacobian.
Now, there are '3' smoothed eigenvalues associated with matrix J*. By using the Matlab
symbolic toolbox we derived the following formulas:
dy =-1/4+3/80000*S-3/80000*I+1/80000*(4e+008-120000*S-120000*1+9*S”2-
18*1*$+9*192)(1/2)
do = -1/4+3/80000*S-3/80000*1- 1/80000*(4e+008-120000*S-120000*1+9*8/2-
18*1*S+9*1A2)\(1/2)
43 =0
Matlab computations showed that the third smoothed eigenvalue is always zero. Now, by
substituting the values of S, I, and R at ay time instant (obtained via simulation) in the above
formulas, we can obtain the corresponding values for the smoothed eigenvalues. . The figures
below plot the real and imaginary parts of the first and second smoothed eigenvalues across
time.
9 5 10 20 25 30 |—Re[A1]
—Re[A2]
Time
Figure 3: Real parts of the smoothed eigenvalues across time
imp}
— Imp]
Figure 4: Imaginary parts of the smoothed eigenvalues across time
The above figures show that there are two bifurcations points.
¢ Time 6: 2 real positive modes merge to one complex mode
e Time 12.5: 1 complex mode split to 2 real negative modes
So basically one can chop time into the following 3 phases:
1. Exponential growth [0,6]
2. Part of oscillation cycle [6,12.5]
3. Exponential decay [12.5, oo ]
3. Conclusion and F uture Work
In this paper, we adopted the idea of Nathan Forrester to consider the second order term of
Taylor series expansion, when approximating nonlinear models. Moreover, we invented a
new concept, which is the smoothed Jacobian matrix. This matrix is valid in any point in the
n-dimensional space -- in contrast to the ordinary Jacobian which is valid only in the
neighborhood of the current operating point; as in dealing with the ordinary matrix, we ignore
the 2" order term of the Taylor series expansion. This can be considered an error term, and as
the width of the analysis interval increases, this error term can be significant. In this paper we
postulate that in nonlinear models, it is preferred to use the smoothed Jacobian instead of the
ordinary Jacobian. Recall that the elements of the smoothed matrix can be considered
smoothed slopes of relationships. The smoothing process is actualized via incorporating the
curvatures specified by the Hessian matrices. Finally, by tracking the values of this smoothed
matrix across time; and symbolically computing the corresponding smoothed eigenvalues and
eigenvectors, one can identify the modes of behavior at any time instant.
In our future work, we will pursuit the following 4 directions:
1. Apply the approach to larger and more complex models.
2. Assess the impacts of changing the gains of loops on the smoothed eigenvalues. In
this point, we will follow the procedure outlined by Kampmann & Oliva 2006.
3. Assess the impacts of changing the values of parameters on the smoothed
eigenvalues and weights. In this point, we will follow the procedure outlined by
Saleh et al. 2006.
4. Extend the analysis to include higher orders terms of Taylor expansion (third and
above) which will be relevant in highly nonlinear models.
Bibliography
Abdel-Gawad A; Abdel-Aleem B; Saleh M; Davidsen P. (2005) "Identifying dominant
behavior patterns, links and loops: Automated eigenvalue analysis of system dynamics
models" Proceedings of the 2005 International System Dynamics Conference, Boston, Mass.
Chen, C.-T. (1970) Introduction to Linear Systems Theory. Holt, Rinehart and Winston, Inc.:
New York, NY.
Edwards C, Penny D. (2005) Differential equations and linear algebra. 2nd ed. Pearson
Education, Upper Saddle River, NJ.
Forrester N. (1982) "A Dynamic Synthesis of Basic Macroeconomic Policy: Implications for
Stabilization Policy Analysis" PhD Thesis, Sloan School of Management, Mass. Inst. of
Technology, Cambridge, MA.
Giineralp B. (2005) "Towards Coherent Loop Dominance Analysis: Progress in Eigenvalue
Elasticity Analysis" Proceedings of the 2005 Int. System Dynamics Conference. Boston.
Kampmann C.E. (1996) "Feedback Loop Gains and System Behavior" Paper presented at the
1996 International System Dynamics Conference, Cambridge, MA.
Kampmann C.E.; Oliva R. (2006) "Loop Eigenvalue Elasticity Analysis: Three Case Studies"
System Dynamics Review 22(2): 146-162.
Luenberg D. (1979) Introduction to dynamic systems: Theory, models and applications.
Wiley, New York, NY.
Ogata K. (1990) Modern Control Engineering. 2nd ed. Prentice Hall: Englewood Cliffs, NJ.
Oliva R, Mojtahedzadeh MT. (2004) "Keep it simple: Dominance assessment of short
feedback loops" Proceedings of the 2004 International System Dynamics Conference.
Oxford, UK.
Oliva R. 2006. Rogelio Oliva's System Dynamics Resource Page. August 1, 2006, 2006.
Available from http://iops.tamu.edu/faculty/roliva/research/sd/.
Saleh M, Davidsen P. (2001a) "The origins of behavior patterns" Proceedings of the 2001
International System Dynamics Conference. Atlanta, GA.
Saleh M, Davidsen P. (2001b) "The origins of business cycles" Proceedings of the 2001
International System Dynamics Conference. Atlanta, GA.
Saleh M, Oliva R., Davidsen P., Kampmann C. (2006) " Eigenvalue Analysis of System
Dynamics Models: Another Perspective" Proceedings of the 2006 International System
Dynamics Conference.Nijmegen, The Netherlands.
Sterman J.D. (2000) Business dynamics: Systems thinking and modeling for a complex world.
Irwin McGraw-Hill, Boston, MA
10
