- Email: [email protected]
Computational techniques for analysis of system dynamics models of social systems
COMPUTATIONAL SYSTEM
TECHNIQUES FOR ANALYSIS DYNAMICS MODELS OF SOCIAL SYSTEMS
OF
JAMESR. BURNS* Department of Engineering Analysis and Design, Texas Tech. University, Lubbock, Texas 79409 and
DAVIDW. MALONE* Battelle Columbus Laboratories, Columbus. Ohio 43201
This paper presents a methodology for analysis of system dynamics models of social systems. The methodology consists ofcomputational algorithms. Specifically, the computational techniques of sensitivity analysis as extended to system dynamics models are shown to enhance an understanding of model behavior and to contribute to considerations of model validity, whereas optimization analyses are shown to contribute to model utility and policy formulation.
INTRODUCTION Recently, many renown statesmen and scientists are becoming increasingly alarmed by the rapidity with which our nation and the world is using up natural resources, polluting its environment, destroying the natural ecology, and overpopulating its landmasses. At the same time there is increasing recognition that the systems approach and its associated methodologies provide a technological resource capable of facilitating man’s understanding of these complex situations with the intent of ameliorating or alleviating such problems. Thusoneofthegreatestchallengesfacingmodernsystem theory today is the modeling and analysis of these socalled societal systems. Two basic system methodologies are required-a methodology for modeling the societal system, and a methodology for validating, verifying, and utilizing all such models. This paper presents a methodology for model validation and utilization via the extension of control theoretic notions which previously have enjoyed application only to mathematical models employed in the physical sciences. The recent works of Forrester and Meadows [l-3] have demonstrated the usefulness of continuous, nonlinear, dynamic, feedback models to represent societal systems. This modeling technique, which has been aptly named ‘system dynamics’, has achieved transdisciplinary appeal and comprehensibility as attested by the very broad range of critical review. There have recently appeared a number of publications critiquing these models or otherwise debating their contents and conclusions. References [4-S] all represent attempts to * The authors were formerly with the Center for Large Scale Systems, Purdue University. Lafayette, Indiana 47907.
invalidate or otherwise demonstrate the tenuous and perhaps dubious assumptions and approximations under which the models described in [2, 31 were developed. Reference [6] discusses the adequacy of system dynamics in general to represent world-wide interactions as a consideration independent of the accuracy and validity of any particular model. Reference [7] is a collection of edited papers which discuss the relative validity of the World 2 (Forrester) and the World 3 (Meadows) models. The structure of major subsystems within each of the models is analyzed and pertinent discussions of Malthusian approaches to understanding the state of the world system are provided. Reference [9] contains a collection of papers which analyze and critique the system dynamics model described in [ 11. None of these references, however, make use of analytical techniques for validating and controlling the dynamical behavior of a system dynamics model. In recognition of this fact, this paper presents algorithms for validation and utilization of system dynamics models of societal systems using the model described in [2] as the vehicle of demonstration. The purpose of the paper is to present a collection of algorithms which could comprise a methodology for validation and utilization of system dynamics models and to demonstrate their usefulness in the context of a controversial social system model. The objective of the paper is not, however, to support or invalidate the criticisms of the model found in the references mentioned above since these criticisms are qualitative and subjective, whereas the approach being suggested here is both quantitative and objective. Moreover, most of the criticisms cited are directed towards the structure of the model, whereas the methodology being advocated by this paper is more
215
JAMES R. BURNS
216
and
directly concerned with model validity and utility as related to the settings of parameters and initial conditions. The paper first discusses briefly Forrester’s world model 121. Then, algorithms for sensitivity analysis are developed for nonlinear dynamic models of social systems. The algorithms are used to compute sensitivity coefficients which, when plotted against time, provide time histories of changes in state (level) variables due to small but constant parameter changes introduced at a specified time. These sensitivity coefficients are the quantitative representation of what Forrester refers to as ‘counterintuitive’ system behavior. Then, the notion ofoptimality is introduced via a mathematical construct ofa value system.The techniques ofparameter optimization and variational calculus are then applied to extremize the construct. Parameter optimization offers considerable advantage over conventional trial and error schemes because the latter are dependent upon the analysts’ intuition. By treating the parameters as functions of time it is shown that, through deployment of the variational calculus, even better system behavior may be obtained than is possible when the parameters are treated as mere constants only. An open-loop policy of near-optimal control is determined which demonstrates that. by treating the parameters as time-varying. it i\ not ncccz5at-! to constrain all of the paramctcr~ in order to achicvc the desired behavior. These algorithms, which are derived from technology in everyday use by system theorists, are shown in this paper to produce useful insights into the behavior of a model whose form admits ready communication with a broad segment of the population. Further, use of these techniques provides more objective input to the intense controversy which has been generated by the publication of World Dynamics and The Limits to Growth, FORRESTER’S
WORLD
DAVID
W.
MALONE
on the six-dimension function space (x, t). Such a model can be described by a system of nonlinear, first-order, coupled differential equations. In order to facilitate the discussion which follows, the model shall be delineated in the following way: i = fCx(t), u,
tl,
.x(q)) = xg.
(2)
This reformulation did not change the operation or behavior of the model, as Forrester’s version is simply a difference equation formulation of equation (2). The complete model has slightly over 100 additional auxiliary variables and parameters, including some associated with the integration and printing logic. In Chapters V and VI of [2], Forrester presents an extensive discussion of the implications for system behavior of postulated changes in critical parameters ofthe model. The following analysis considers the same set of (analogous) control parameters as were examined by Forrester: Ul
birth rate normal, BRN
u2
natural resource usage normal, NRUN
u=
uj
= pollution
normal, POLN
capital investment in agriculture fraction normal, CONAF I.
capital investment CIGN.
a5
For purposesofthis parameter vector
generation
study the components of the control have been normalized about their Sea le
3
\~ \
\
-\
normal,
2
z
NR
2.5~10”
--
POLB CIAF I 5
.-.-. ---Cl
\
factors
P
x109 x IO’O x IO-’ xl09
MODEL
The model of the world developed by Forrester [Z] is a computerized abstraction of the real world which employs system dynamics as the simulation context within which the model is embedded. The model consists of a set of equations which describe the interaction between five major ‘level variables’ identified as population, natural resources, pollution, fraction of capital devoted to agriculture, and capital investment. Using state-space notation, the control theorist denotes these variables as state variables, or as components of a state vector defined as follows: population,
P
natural resources, = pollution,
NR
POL
(1)
capital investment CIAF
in agriculture,
capital investment,
CI
1973
2000
2025
2050
2075
2100
Time
Fig. I. Baseline model response assuming current policies prevail.
System dynamics models of social systems nominal 1970 values. Such normalization causes each component to assume a value of unity as its nominal value and facilitates the deployment of certain computational techniques to be subsequently described. By using the nominal (see [Z]) parameter settings and integrating the mode1 represented by (2) from the year to = 1973 to the year tf = 2100, the results shown in Fig. 1 were obtained, where each of the five states xi(t) are plotted for the time interval 1973 to the year 2100. Figure 1 is to be compared with Fig. 4-1, p. 70 123 where the same states (excepting x4, CIAF) are plotted for the time interval 1900-2100. In this analysis an eighth order, variable-step, Runge-Kutta integration algorithm was used to perform the integration of equation (2) on a CDC 6500. Since the algorithms to be discussed in the sequel are applicable to all models which can be represented by (2), a justification of equation (2) as sufficiently general to encompass all or nearly all system dynamics models is required. System dynamics models consist of two basic equations: the level equation,
ri(t + At) = ii(t) + At. /?T. r(t),
=
Wjft),
MO,
. .
. , L(t);
.
Uj,Ub..., u,, t], i = 1, . . . ,2n, where the rates may be functions of more than one level, and of a set of parameters denoted by uj, ulr, . . . , u,. An additional equation of system dynamics, called the auxiliary equation, can be incorporated in the rate equation without loss of generality. Thus the system is described by just two equation types which are rewritten here in vector notation: Z(t + At) = l(t) + At. P . r(t), r(t) = Wt), u, Cl, where P is an n x 2n matrix of zeros, ones, and minus ones, r(t) is a vector of length 2n, and I(t) is a vector of length n. These equations may be combined in the following way: z(t) 1
Z(t+ At) - I(t) At = P.r(t) = P * h[l(t), u, t] = g[l(th u, t].
By changing notation, one arrives at the representation given by (2): f = f[x(t), u, t]. There are two situations that arise in system dynamics models which require special attention. They are delays and clip functions. Delays can be accommodated by the following transformation. Consider the following general model with delays: l(t) = g[l(th l(t - ax u, t] S-E.P.S. R/‘&-o
where r(t - cr)isintended to denote the vector ofdelayed states. By defining XTA [l(t)=, l(t - c#], j-Lx(t), u, t-JTA (g[l(t), I(t - a), II, cl*, ict - a)? it is clear that f = fCx(0, % tl as before. Therefore the representation given by equation (2) is sufficiently general to encompass models containing delays. Clip and switch functions, on the other hand, cannot be accommodated in every instance. They are permissible as long as they do not produce states Xi with discontinuities or with discontinuous derivatives at specific points in time. The succeeding analysis is contingent upon the requirement that x(t) be continuous and differentiable over the whole time interval to I t 5 tf . Clip functions which violate this requirement must be removed in order to proceed with the analysis suggested in this paper.
i = 1,. . . , n,
where /&“.r(t) denotes the sum or difference of all rates affecting level Ii (pi is a vector of length 2n which consists of zeros, ones, or minus ones); and the rate equation, rift)
217
SENSITIVITY
ANALYSIS
Sensitivity analysis implies a study of the sensitivity of the system’s response to changes (or perturbations) in the values of its parameters and initial conditions. Reference [lo] is a genera1 theoretical treatment of the subject. By studying the sensitivity of the mode1 to changes in its parameters and initial conditions, it is possible to gain insight into the behavior and validity of the model. The fact that sensitivity analysis leads to simple initial value problems makes it ideal for mechanization on a digital computer. The investigation enables the analyst to determine which parameters and initial conditions most strongly affect a given level or state variable within the model. This section develops a generalized algorithm for sensitivity analysis of system dynamics models represented by (2). The algorithm is used to compute sensitivity coefficients which provide time histories of changes in state (level) variables due to a small but constant parameter or initial state change. Since the inputs to any system dynamics model are the parameter and initial state vectors, the objective is to measure the sensitivity of any given state variable to a perturbation in any given input which occurs at time to. Let sT = (uT, XT) be the vector of inputs. In order to determine the sensitivity of x, the solution of (2) to a perturbation in the values of the input vector s. (2) is differentiated with respect to s to obtain: di
z=
af a~ af ax'as+as'
Since x is differentiable and continuous,
a2 d ax as=iitds. 0
(3)
218
JAMES R. BAIRNSand DAVIII W. MALONF
The following definitions are required:
where
is the nominal state trajectory due to nominal initial condition and parameter settings. By inserting these definitions into (3), the following expression is obtained: g(B)
= B = F(+, t) . B + G($, t) (4) B(& to) = Bc,.
Given that x is a vector of length n, u a vector of length m, it should be apparent that in general B is a matrix of dimension n x /, F is a Jacobian matrix of dimension II x II. and G is a Jacobian matrix of dimension II x 1.where I = II + m. For the world model being considered, B has dimensions 5 x 10, F has dimensions 5 x 5, and G has dimensions 5 x 10. The matrices G($, t) and B,, require further explanation. As previously defined. . G($. I) 4 ;;; howcvcr.
This normalization permits ready comparison of the relative influence which any component of the input vector sT = (uT, xt) has on any given state xi(t). The complete sensitivity analysis algorithm involves the simultaneous integration of equations (2) and (4) via the digital computer and numerical integration methods. The algorithm requires that first partial derivatives dfJaxjl and afi/auj be computed using analytically determined expressions for each derivative. Alternatively, each sensitivity coefficient could be computed using numerical approximation methods. i.e. ? xAs + As, c) - xds, t) bij = z z As; ’ / I
where As“ = (O,O, , Asj, . ,O). The latter methodwasn’t used in thisanalysis because of the requirement to plot all sensitivity coefficients relative to a given state at any one time. The accomplishment of this requirement would entail large amounts of additional computer storage of sensitivity coefficient trajectories for later plotting. Secondly, use of algorithm (6) requires an additional five integrations of the system equations (2) for the model being considered. These disadvantages, coupled with the fact that (6) is an approximate technique, led to the implementation of algorithm (4) for the model being considered.
x(s + As. t) = .Y(s,t) + B(& t) As. Since
-af = O(nx n).;=[$,O(nx ax,
n)]=G@,t),
where O(n x n) is the n x n zero matrix, and af/au is a matrix of dimension n x tn. The matrix B, = B(4, to) is determined from a consideration of the original definition of B:
L
0
X0”
1
(7)
This equation exposes the interpretation of sensitivity coefficients as the time-dependent coefficients of a Taylor series expansion about the nominal state trajectory &s, t) due to nominal parameter and initial condition vector settings.
RESULTS
Thus
(6)
OF THE SENSITIVITY
ANALYSIS
Integration of the generalized algorithm represented by (4) together with the model equations represented by (2) leads to the results shown in Fig. 2. Figure 2(a) contains time histories of all parameter sensitivity coefficients associated with population for the time interval 1973-2100, whereas, Fig. 2(b) contains time histories of all initial condition sensitivity coefficients associated with population. The interested reader is referred to [ll] for plots of the remaining sensitivity coefficients. Figure 2, together with implications of equation (7), provides additional quantitative information which can be introduced into substantive discussionsof the validity of the model which produced the results shown in Fig. 1. For example, estimates of the changes in amplitude of population which would accompany specified small changes in the components of the parameter vector can be obtained almost by inspection of Fig. 2(a).
System dynamics models of social systems
Time
Fig. 2(a). Parameter sensitivity coefficients associated with population. P.
The following remarks briefly summarize the main observations obtained from the parameter sensitivity
analysis of the Forrester model: 1. The strongest influence ofa parameter change occurs at a considerable length of time after the change, sometimes as much as 100 yr later. 2. The most influential parameter is usually, but not always, that directly associated with the state: an exceptionis theeffect ofachangeincapitalinvestment generation normal (CIGN) on natural resources and pollution. This high sensitivity to CIGN has been noted elsewhere [ 151. 3. Changingaparameter to influence its associated state variable may have considerable influence ..on other states. Note that these observations have previously been made by Forrester in a qualitative sense. The sensitivity coefficients provide a quantitative assessment of the delay times, relative influence and second-order effects. The following observations summarize the main results of the initial condition sensitivity analysis for the Forrester world model: 1. The results demonstrate that all of the states are most sensitive to the initial setting of natural resources. The influence which this initial state has is seen to grow and become increasingly important as time progresses, except in the case of natural resources itself, where the influence. is observed to decrease with time. 2. Most state variables indicate a very high degree of sensitivity to the setting of the corresponding initial
219
state for the first few years of the time period. In every instance that sensitivity is seen to decrease as time progr&ses. Sensitivity analyses contribute both to model validity and to model utility. The sensitivity analysis of the initial states contributes to model validity because of the insight obtained regarding how small variations in the initial states influence the states at future points in time. Thus, a high degree of accuracy in the initial setting of natural resources is desirable because of the dominating influence which thisinitial condition is observed to have. On the other hand, the parameter sensitivity results provide insighj as to which parameters are the most influential for controlling a given state. Thus CIGN would represent an appropriate parameter for controlling pollution and for controlling the rate of natural resource usage. In summary, sensitivity analysis represents a mechanism for measuring the sensitivity of any given input with any given output, and describes how changes in the input parameter would influence future behavior. Relative to model validity, sensitivity analysis suggests data collection priorities to improve the accuracies of the more sensitive parameters and initial conditions in the model. Relative to model utility, sensitivity analysis provides insight as to how changesina particular control parameter would influence the model state at some future time. In Fig. 2(a), increases in POLN and NRUN are seen to have a negative effect on population. The fact that an increase in NRUN decreases population is counter to one’s intuition regarding the actual cfiect an
Scale
I unit=
factor
---.-.-. -.-
1
IQ‘73
--.F’\.
1,147~ IO9 PO NR, POL, CIAF, Cl0
I 2000
I 2025
I 2050
I 2075
I 2Kx)
Time
Fig. 2(b). Initial state sensitivity c&fficients associated with population. P.
JAMESR. BURNS and DAVID W. MALONE
220 increase
in NRUN would have on population. This and all such similar observations provides additional information regarding the dynamical behavior of the model andexplainswhysensitivitycoefficientsare the quantitative representation of what Forrester refers to as “counterintuitive” system behavior.
parameter vector. Precisely stated, the problem is to determine I*@*), where ‘I I*@*) = min rl(QL, P, u, t) df WM I It, subject to: x==(x,u,t),
OPTIMIZATION
ANALYSES
In this section two methods for optimizing system dynamics models of social systems are presented, using the Forrester world model as the vehicle of demonstration. Specifically, the extensions necessary to apply (1) parameter optimization techniques and (2) the variational calculus are discussed together wirh resultant implications for behavior of the world model. Suchanalytical techniquesgenerateadditional objective informationand insights regarding model behavior, and hence contribute to substantative deliberations regarding model validity. Before any optimization analysis can be performed, a mathematical construct of the value system used to measure the performance of the model must be developed. In this study the following index- was developed to measure the performance of the model:
where r/(QL. P. u, t) = -C,
.QL2 + Cz (P - P,,&* +
c3
.IIu
-
111;.
(8)
The index is intended to mathematically characterize the valuesystemused by Forrester to evaluate the model behavior.Thefirst term to appear in the index represents a desire to maximize the quality-of-life, QJ!,, for the time interval [to, ts]. The second term represents an intent to stabilize the population, while the third term measures the amount of control effort required to accomplish the objectives represented by the first two terms. In words this index stipulates that quality of life be maximized and population be stabilized while minimizing the extent of intervention into the system required to achieve these goals. The constants Cr, C2, and C3 are used to weight each term appropriately and Q is a five-dimensional identity matrix for the model being considered. It should be emphasized that other value systems could be used to measure model performance and that other different mathematical constructs exist for monitoring model performance. The parameter optimization problem requires choosing a parameter vector u to minimize the index of performance defined by (8) subject to the system differential constraints given by (2) and certain physical limitations on the choice of values available for the * The
values for BRN
and CIGN
are normalized
values.
x(t,)=x,,
M=
:u(ui20,i=
1,5).
The Powell method [12] was selected for use in this study because it was found to provide faster convergence than other function minimization schemes, yet it does not require analytic expressions for the first partial derivatives of the index of performance. In order to prevent the method from choosing parameter values which are physically unrealizable, penalty functions were employed. The penalty function defined below was chosen from several other alternative configurations which seemed less well suited to this problem. I+l+Cgi and 0, ui 2 0 gi = 1 c. l&l. ui < 0.
C = 1 x 105.
The parameter optimization program was started with an initial estimate for the parameter vector of ur = (1, 1, 1, 1, 1). This initial estimate corresponds to the situation in which the world is permitted to run its course in the absence of restraints being imposed upon its behavior. This initial estimate produced a value for the index of performance of 12,800. In the first 120 function evaluations, the index of performance decreased from 12,800 to -40. In the last 140 function evaluations, the index decreased from - 40 to - 46. The parameter specification at the end of the 260 function evaluations was determined to be u*r = (0.60, 0.52, 1.16, 0.89, 0.71). These values are to be compared with the set of values which were hand-picked by Forrester for the explicit purpose of producing an improvement in the system behavior relative to his implicit value system. These values (see [Z], p. 120) are: uT = (0.70, 0.25, 0.50, 0.80, 0.60).* Some interesting contrasts appear between the handpicked values of Forrester and the values picked by the computer. The computer-determined value for birth rate normal is more severe than that chosen by Forrester and therefore harder to achieve. However, the restraints applied to natural resource usage would be less severe under the computer-determined policy as opposed to the Forrester policy. Interestingly enough, the computer-determined policy suggests allowing pollution generation to increase rather than restraining it as the Forrester policy recommends. The resultant plot did not indicate excessive pollution. The values chosen for agricultural production CONAF and for CIGN are less restrictive than those corresponding values determined by Forrester, and therefore easier to achieve. It is interesting to note that the parameter values picked by the parameter optimization algorithm
System dynamics models of social systems
are in the interior of the feasible region of the parameter space M, and not on the boundary. In contrast to the previous approach, the strategy of the open-loop control approach is to treat the parameters as functions of time and to select the optimal set of time functions which will minimize the index of performance.Thustheproblemistoselectanadmissable u*(r) such that ‘I I*@*) = mm
u&f s to
rr specified
and i = 1,. . . (5. to 4 t 5 cl}.
Due to the high degree of nonlinearity and coupling which resides in the system equations, considerable computational irregularities enter into the ‘problem posed above. The table functions mf(x, u, r) were fitted with polynomials so as to provide continuous derivatives. The index of performance was converted from its Lagrangian format shown above into the form of Mayer. Then the optimization problem can be reformulated as follows: Z*(u*) = min x,(t,) uelu subject to : 2 = f’(z, u, t), .z(t,,)= zo, tf specified where: zT = (XT,X6) .fT = (.fr. rl(QL, P, u, 0). Weak necessary conditions for a stationary point within the control subspace M for the problem posed above require that the first variation of the augmented index of performance vanish. The augmented index is defined as Al = x,(t,) +
” IT{& I (0,
u, r) - i) dt.
ax&,) - (;.T.l%):~,= 0. AT $z,
LT
A = ‘;(z,
af * y-&
u)
(2, u)
=
u). (13)
These adjoint (Euler-Lagrange) equations are solved simultaneously with the system equations: z(t0) = zo.
(14)
The resultant problem is generally referred to as a twopoint boundary-value problem (TPBVP). A great many methods have been proposed for solution of the TPBVP, many of which were applied without success to the problem posed above. Any method which required the simultaneous integration of both the system and adjoint equations, whether forward integration, backward integration, or both is required, was generally met with failure. These disappointing results eventually led to the selectionofthefirst-ordergradientmethod developed by Denham and Bryson, and discussed in Citron [I33 andBrysonandHo[14].Thissteepestdescent technique relaxes condition (12) while requiring that the system and adjoint equations be satisfied at all times. The relaxation removes the coupling of the control vector which would otherwise exist between the system and adjoint equations. In addition the technique permits the system equations to be. integrated forward in time while the adjoint equations are integrated backward in time, thus taking advantage of the fact that the system equations are numerically stable during forward integration while the adjoint equations are stable during backward integration. The first-order gradient method consists of the following steps: 1. Estimate a control vector u(r) and integrate the systemequations(14) to the terminal time t/, storing the state trajectory x(t). 2. Integrate the adjoint equations (13) backward and store the resultant adjoint trajectory. 3. Computeanewcontrol trajectory using the following formula :
(9)
Taking the first variation of (9), noting that the continuity of z(r) permits interchanging the operations of differentiation and variation, and taking advantage of the functional independence of r(t) and u(r), it can be shown that the weak necessary conditions for the problem considered become:
AT +
1 = -AT . L, n(QT = (0, 0, 0, 0, 0, l),
z = s’(z, u, tb
subject to
M = [ulu&) 2 0,
From (10) it is determined that n(tJ = (0, 0,0,0,0,1). Thus the adjoint system is defined by:
tl(Qk P, u, 0 dt
i = f(x, u, t), tit,,) = x0,
221
(10)
0,
(11)
=I 0.
(12)
U new =
%ld
-
wi:p
where
and Q is the (m x m) unit matrix, dS is the step size. 4. Repeat steps 1, 2 and 3 until no further improvement in the index of performance Zcan be obtained. OPTIMIZATION ANALYSIS RESULTS
The open loop optimization program implementing the method just described was started from several
JAMESR. BURNS and
222 Sea le factors
I‘-__
--
-P -----.-.-.-
-.
/
. . .’ x ./
./’ 2cI
_ ror-
.’
.’ -
2x109 NR 25x IO” POL8 x IO’O QL 5 x IQ’ Cl 5 x 109 . .._‘N. ‘. ‘. ‘\ \ ‘\ ‘\ ‘\ ‘\\ ‘\ ‘\ ‘\ . ‘\ ‘1. .-. _-‘-‘\_
.A
C‘
I’
0’
1.
./
/
--
‘\
,/
-- --173
l
2000
I
2026
I
2050
I
2075
I
2100
Time
Fig. 3(a). Model response to control trajectories pictured in Fig. 3(b).
different points in the control subspace M and allowed
to continue its iteration until very little additional improvement in the performance index was realized. The terminal results were then plotted. Figures 3 and 4 show the terminal results obtained by starting from two different points, ml and m2 in M. In Fig. 3(a) is shown the model response induced by the control setting shown in Fig. 3(b) (such control setting determined by starting the iteration at point ml), and similarly for Figs. 4(a) and 4(b). The starting point m, EM, from which Fig. 3 was obtained is uT = (0.70,0.25,0.50,0.80, 0.60) where u is assumed to be constant over the whole time interval. The attentive reader will recognize these values as the optimal parameter values determined by the parameter optimization. The starting point m2 EM from which Fig. 4 was obtained is uT = (1, 1, 1, 1, l), or the nominal, unconstrained parameter values associated with the current policies of the early 1970%. In both runs the iterative steepest descent algorithm was allowed to continue until no further reduction in the index of performance could be obtained. The performance value associatedwithFig.3(b)is - 127whereas that associated with Fig. 4(b) is -102. In both figures the control trajectories were constrained to the nominal values corresponding to the present policies for the first 3 yr of the simulation interval by weighting heavily the third term of the index of performance for this period. The multimodal response surface defined by (8) is responsible for the obvious observation that the results in Figs. 3 and 4 are different. The consequence of starting from two different points in the control subspace is that convergence to two different local minima defined by the response surface has been achieved. As
DAVII)
W.
MALONI:
a result, the analyst is faced with the problem of choosing between two different control policies, each about equally ‘good’ insofar as goodness is defined by the performance measure (8). Apparently, the system trajectories must be evaluated for desirability, the control trajectories for feasibility. The system trajectories shown in Fig. 3 are more desirable than those shown in Fig. 4 because quality of life is higher and natural resource usage is lower. However, the control trajectories shown in Fig. 3 are less desirable than those shown in Fig. 4 because of the increased amount of control effort required to achieve the higher quality of life. Since the performance indexconsiders both of these factors (quality of life and control effort), it ‘measures’ the performance observed in Figs. 3 and 4 as about equal-a minus 127 for Fig. 3 compared with a minus 102 for Fig. 4. These numbers, as compared with the performance value of 12,800 associated with Fig. 1. are roughly the same. The continuous, time-varying control strategies suggested here appear to be much more realistic than discrete interventions into societal systems because discrete changes, especially if they are dramatic, are psychologically resisted. The implication is that continuous time-varying control of social systems is more reasonable and less restrictive than one is led to believe from the results of Forrester’s analyses. In particular, Fig. 4(b) suggests a time-varying control strategy which manipulates or constrains only two of the five control variables. The remaining three control parameters retain their nominal values of 1.
LcOLN______
/ -___
,CONAF
\y..-._.-._._L._.-.-_.--.-.-.-.__
0.90
0.30
,-
i
C ‘19 73
I
2000
I
2026
I
2050
I
2075
I
2100
Time
Fig. 3(b). Optimal control trajectories obtained by starting from point I in the control suhspncc M.
System dynamics models of social svstems
Scale A
_--0
1973
factors
I.20
P 2 x 10s NR 2.5~ IO”
----
3 A\
t
_-----__ I 2000
I 2025
223
0.30 I 2030
I 2073
I 2100
Time
Fig. 4(a). Model response induced by control shown in Fig. 4(b).
t
histories
0~0
Time
Fig. 4(b). Optimal control trajectories determined by starting from point 2 in the control subspace M.
CONCLUSION
This paper presents several algorithms for analyzing and optimizing system dynamics models of social systems. The usefulness of these generalized methods is demonstrated via their applicatian to Forrester’s moclel of the world. The need for such a demonstration was
discussed in the introduction. There attention was directed to the recognition amongst statesmen and scientists alike of the facility which system analysis has to offer to the study and analysis of social systems. In this framework, many such problems can be viewed asarisingfrom the inability to synthesize and implement policieswhichwouldsteerthesesystemsfromthepresent undesirable states to more ideal future states. The methodology suggested in this paper addresses this problem by demonstrating the use of automated techniques for policy determination. In addition, a review of previous related studies showed that computational algorithms have never been applied to the controversial system dynamics models Cl-33 of today for the purpose of contributing to deliberations regarding model validity. Again, the methodology suggested here is intended to address the urgent requirement for mechanized techniques which would contribute to considerations of model validity and enhance model understanding.
REFERENCES 1.
J. W. Forrester, ,._,,.I
(IYOYJ.
UrbanDynamics.
MIT Press, Cambridge
2. J. W. Forrester, World Dytuxmics. Wright-Allen Press. Cambridge (1971). 3. D. H. Meadows. D. L. Meadows (‘I ~1.. T/K* Limits to Growrh. Universe Books New York (1972). 4. R. Boyd, World dynamics: a note. Sci~ncu. 177, 516-519 (1972). 5. C. Starr and R. Rudman, Parameters of technological growth, Astron. Aeron. 20-27 (1973). 6. E. S: Boylan, The system dynamics approach to modeling world-wide interactions:8 critical analysis. Unpublished Report, Department ofMathematics, Rutgers University (1972). 7. H. S. D. Cole. C. Freeman et al.. Models o/ Doom Universe Books, New York (1972). 8. M. H. Rothkopf. World models won’t work. Simtrltrtio~~ -61 9.
ilY73).
A. P. Sage, Ed., SMC-2
10.
IEEE Trans. (2), 122-220 (1972).
Systems, Man,
R. Tomovic and M. Vukobratovic,
General
Cybernetics Sensitivity
Elsevier, New York (1972). 11. J. Burns, Applications of control theory to system dynamics models of social systems, Ph.D. Dissertation, Purdue University (1973). 12. M. J. D. Powell, An efficient method for finding the tiinimum of a function of several variables without calculating derivatives, Computer J. 7, 155-162 (1964). 13. S. J. Citron, Elements of Optimal Control. Holt, Rinehart and Winston, New York (1969). 14. A. E: Bryson, Jr. and Y. C. Ho. Applied Optimd Cortfrrd Blaisdell, Waltham (1969). 15. R.A. Burnett and P. J. Dionne, GLOBE6: a multircgion interactive world simulation, Simuhtion 192. 197 ( 1973). Theory.
TECHNIQUES FOR ANALYSIS DYNAMICS MODELS OF SOCIAL SYSTEMS
OF
JAMESR. BURNS* Department of Engineering Analysis and Design, Texas Tech. University, Lubbock, Texas 79409 and
DAVIDW. MALONE* Battelle Columbus Laboratories, Columbus. Ohio 43201
This paper presents a methodology for analysis of system dynamics models of social systems. The methodology consists ofcomputational algorithms. Specifically, the computational techniques of sensitivity analysis as extended to system dynamics models are shown to enhance an understanding of model behavior and to contribute to considerations of model validity, whereas optimization analyses are shown to contribute to model utility and policy formulation.
INTRODUCTION Recently, many renown statesmen and scientists are becoming increasingly alarmed by the rapidity with which our nation and the world is using up natural resources, polluting its environment, destroying the natural ecology, and overpopulating its landmasses. At the same time there is increasing recognition that the systems approach and its associated methodologies provide a technological resource capable of facilitating man’s understanding of these complex situations with the intent of ameliorating or alleviating such problems. Thusoneofthegreatestchallengesfacingmodernsystem theory today is the modeling and analysis of these socalled societal systems. Two basic system methodologies are required-a methodology for modeling the societal system, and a methodology for validating, verifying, and utilizing all such models. This paper presents a methodology for model validation and utilization via the extension of control theoretic notions which previously have enjoyed application only to mathematical models employed in the physical sciences. The recent works of Forrester and Meadows [l-3] have demonstrated the usefulness of continuous, nonlinear, dynamic, feedback models to represent societal systems. This modeling technique, which has been aptly named ‘system dynamics’, has achieved transdisciplinary appeal and comprehensibility as attested by the very broad range of critical review. There have recently appeared a number of publications critiquing these models or otherwise debating their contents and conclusions. References [4-S] all represent attempts to * The authors were formerly with the Center for Large Scale Systems, Purdue University. Lafayette, Indiana 47907.
invalidate or otherwise demonstrate the tenuous and perhaps dubious assumptions and approximations under which the models described in [2, 31 were developed. Reference [6] discusses the adequacy of system dynamics in general to represent world-wide interactions as a consideration independent of the accuracy and validity of any particular model. Reference [7] is a collection of edited papers which discuss the relative validity of the World 2 (Forrester) and the World 3 (Meadows) models. The structure of major subsystems within each of the models is analyzed and pertinent discussions of Malthusian approaches to understanding the state of the world system are provided. Reference [9] contains a collection of papers which analyze and critique the system dynamics model described in [ 11. None of these references, however, make use of analytical techniques for validating and controlling the dynamical behavior of a system dynamics model. In recognition of this fact, this paper presents algorithms for validation and utilization of system dynamics models of societal systems using the model described in [2] as the vehicle of demonstration. The purpose of the paper is to present a collection of algorithms which could comprise a methodology for validation and utilization of system dynamics models and to demonstrate their usefulness in the context of a controversial social system model. The objective of the paper is not, however, to support or invalidate the criticisms of the model found in the references mentioned above since these criticisms are qualitative and subjective, whereas the approach being suggested here is both quantitative and objective. Moreover, most of the criticisms cited are directed towards the structure of the model, whereas the methodology being advocated by this paper is more
215
JAMES R. BURNS
216
and
directly concerned with model validity and utility as related to the settings of parameters and initial conditions. The paper first discusses briefly Forrester’s world model 121. Then, algorithms for sensitivity analysis are developed for nonlinear dynamic models of social systems. The algorithms are used to compute sensitivity coefficients which, when plotted against time, provide time histories of changes in state (level) variables due to small but constant parameter changes introduced at a specified time. These sensitivity coefficients are the quantitative representation of what Forrester refers to as ‘counterintuitive’ system behavior. Then, the notion ofoptimality is introduced via a mathematical construct ofa value system.The techniques ofparameter optimization and variational calculus are then applied to extremize the construct. Parameter optimization offers considerable advantage over conventional trial and error schemes because the latter are dependent upon the analysts’ intuition. By treating the parameters as functions of time it is shown that, through deployment of the variational calculus, even better system behavior may be obtained than is possible when the parameters are treated as mere constants only. An open-loop policy of near-optimal control is determined which demonstrates that. by treating the parameters as time-varying. it i\ not ncccz5at-! to constrain all of the paramctcr~ in order to achicvc the desired behavior. These algorithms, which are derived from technology in everyday use by system theorists, are shown in this paper to produce useful insights into the behavior of a model whose form admits ready communication with a broad segment of the population. Further, use of these techniques provides more objective input to the intense controversy which has been generated by the publication of World Dynamics and The Limits to Growth, FORRESTER’S
WORLD
DAVID
W.
MALONE
on the six-dimension function space (x, t). Such a model can be described by a system of nonlinear, first-order, coupled differential equations. In order to facilitate the discussion which follows, the model shall be delineated in the following way: i = fCx(t), u,
tl,
.x(q)) = xg.
(2)
This reformulation did not change the operation or behavior of the model, as Forrester’s version is simply a difference equation formulation of equation (2). The complete model has slightly over 100 additional auxiliary variables and parameters, including some associated with the integration and printing logic. In Chapters V and VI of [2], Forrester presents an extensive discussion of the implications for system behavior of postulated changes in critical parameters ofthe model. The following analysis considers the same set of (analogous) control parameters as were examined by Forrester: Ul
birth rate normal, BRN
u2
natural resource usage normal, NRUN
u=
uj
= pollution
normal, POLN
capital investment in agriculture fraction normal, CONAF I.
capital investment CIGN.
a5
For purposesofthis parameter vector
generation
study the components of the control have been normalized about their Sea le
3
\~ \
\
-\
normal,
2
z
NR
2.5~10”
--
POLB CIAF I 5
.-.-. ---Cl
\
factors
P
x109 x IO’O x IO-’ xl09
MODEL
The model of the world developed by Forrester [Z] is a computerized abstraction of the real world which employs system dynamics as the simulation context within which the model is embedded. The model consists of a set of equations which describe the interaction between five major ‘level variables’ identified as population, natural resources, pollution, fraction of capital devoted to agriculture, and capital investment. Using state-space notation, the control theorist denotes these variables as state variables, or as components of a state vector defined as follows: population,
P
natural resources, = pollution,
NR
POL
(1)
capital investment CIAF
in agriculture,
capital investment,
CI
1973
2000
2025
2050
2075
2100
Time
Fig. I. Baseline model response assuming current policies prevail.
System dynamics models of social systems nominal 1970 values. Such normalization causes each component to assume a value of unity as its nominal value and facilitates the deployment of certain computational techniques to be subsequently described. By using the nominal (see [Z]) parameter settings and integrating the mode1 represented by (2) from the year to = 1973 to the year tf = 2100, the results shown in Fig. 1 were obtained, where each of the five states xi(t) are plotted for the time interval 1973 to the year 2100. Figure 1 is to be compared with Fig. 4-1, p. 70 123 where the same states (excepting x4, CIAF) are plotted for the time interval 1900-2100. In this analysis an eighth order, variable-step, Runge-Kutta integration algorithm was used to perform the integration of equation (2) on a CDC 6500. Since the algorithms to be discussed in the sequel are applicable to all models which can be represented by (2), a justification of equation (2) as sufficiently general to encompass all or nearly all system dynamics models is required. System dynamics models consist of two basic equations: the level equation,
ri(t + At) = ii(t) + At. /?T. r(t),
=
Wjft),
MO,
. .
. , L(t);
.
Uj,Ub..., u,, t], i = 1, . . . ,2n, where the rates may be functions of more than one level, and of a set of parameters denoted by uj, ulr, . . . , u,. An additional equation of system dynamics, called the auxiliary equation, can be incorporated in the rate equation without loss of generality. Thus the system is described by just two equation types which are rewritten here in vector notation: Z(t + At) = l(t) + At. P . r(t), r(t) = Wt), u, Cl, where P is an n x 2n matrix of zeros, ones, and minus ones, r(t) is a vector of length 2n, and I(t) is a vector of length n. These equations may be combined in the following way: z(t) 1
Z(t+ At) - I(t) At = P.r(t) = P * h[l(t), u, t] = g[l(th u, t].
By changing notation, one arrives at the representation given by (2): f = f[x(t), u, t]. There are two situations that arise in system dynamics models which require special attention. They are delays and clip functions. Delays can be accommodated by the following transformation. Consider the following general model with delays: l(t) = g[l(th l(t - ax u, t] S-E.P.S. R/‘&-o
where r(t - cr)isintended to denote the vector ofdelayed states. By defining XTA [l(t)=, l(t - c#], j-Lx(t), u, t-JTA (g[l(t), I(t - a), II, cl*, ict - a)? it is clear that f = fCx(0, % tl as before. Therefore the representation given by equation (2) is sufficiently general to encompass models containing delays. Clip and switch functions, on the other hand, cannot be accommodated in every instance. They are permissible as long as they do not produce states Xi with discontinuities or with discontinuous derivatives at specific points in time. The succeeding analysis is contingent upon the requirement that x(t) be continuous and differentiable over the whole time interval to I t 5 tf . Clip functions which violate this requirement must be removed in order to proceed with the analysis suggested in this paper.
i = 1,. . . , n,
where /&“.r(t) denotes the sum or difference of all rates affecting level Ii (pi is a vector of length 2n which consists of zeros, ones, or minus ones); and the rate equation, rift)
217
SENSITIVITY
ANALYSIS
Sensitivity analysis implies a study of the sensitivity of the system’s response to changes (or perturbations) in the values of its parameters and initial conditions. Reference [lo] is a genera1 theoretical treatment of the subject. By studying the sensitivity of the mode1 to changes in its parameters and initial conditions, it is possible to gain insight into the behavior and validity of the model. The fact that sensitivity analysis leads to simple initial value problems makes it ideal for mechanization on a digital computer. The investigation enables the analyst to determine which parameters and initial conditions most strongly affect a given level or state variable within the model. This section develops a generalized algorithm for sensitivity analysis of system dynamics models represented by (2). The algorithm is used to compute sensitivity coefficients which provide time histories of changes in state (level) variables due to a small but constant parameter or initial state change. Since the inputs to any system dynamics model are the parameter and initial state vectors, the objective is to measure the sensitivity of any given state variable to a perturbation in any given input which occurs at time to. Let sT = (uT, XT) be the vector of inputs. In order to determine the sensitivity of x, the solution of (2) to a perturbation in the values of the input vector s. (2) is differentiated with respect to s to obtain: di
z=
af a~ af ax'as+as'
Since x is differentiable and continuous,
a2 d ax as=iitds. 0
(3)
218
JAMES R. BAIRNSand DAVIII W. MALONF
The following definitions are required:
where
is the nominal state trajectory due to nominal initial condition and parameter settings. By inserting these definitions into (3), the following expression is obtained: g(B)
= B = F(+, t) . B + G($, t) (4) B(& to) = Bc,.
Given that x is a vector of length n, u a vector of length m, it should be apparent that in general B is a matrix of dimension n x /, F is a Jacobian matrix of dimension II x II. and G is a Jacobian matrix of dimension II x 1.where I = II + m. For the world model being considered, B has dimensions 5 x 10, F has dimensions 5 x 5, and G has dimensions 5 x 10. The matrices G($, t) and B,, require further explanation. As previously defined. . G($. I) 4 ;;; howcvcr.
This normalization permits ready comparison of the relative influence which any component of the input vector sT = (uT, xt) has on any given state xi(t). The complete sensitivity analysis algorithm involves the simultaneous integration of equations (2) and (4) via the digital computer and numerical integration methods. The algorithm requires that first partial derivatives dfJaxjl and afi/auj be computed using analytically determined expressions for each derivative. Alternatively, each sensitivity coefficient could be computed using numerical approximation methods. i.e. ? xAs + As, c) - xds, t) bij = z z As; ’ / I
where As“ = (O,O, , Asj, . ,O). The latter methodwasn’t used in thisanalysis because of the requirement to plot all sensitivity coefficients relative to a given state at any one time. The accomplishment of this requirement would entail large amounts of additional computer storage of sensitivity coefficient trajectories for later plotting. Secondly, use of algorithm (6) requires an additional five integrations of the system equations (2) for the model being considered. These disadvantages, coupled with the fact that (6) is an approximate technique, led to the implementation of algorithm (4) for the model being considered.
x(s + As. t) = .Y(s,t) + B(& t) As. Since
-af = O(nx n).;=[$,O(nx ax,
n)]=G@,t),
where O(n x n) is the n x n zero matrix, and af/au is a matrix of dimension n x tn. The matrix B, = B(4, to) is determined from a consideration of the original definition of B:
L
0
X0”
1
(7)
This equation exposes the interpretation of sensitivity coefficients as the time-dependent coefficients of a Taylor series expansion about the nominal state trajectory &s, t) due to nominal parameter and initial condition vector settings.
RESULTS
Thus
(6)
OF THE SENSITIVITY
ANALYSIS
Integration of the generalized algorithm represented by (4) together with the model equations represented by (2) leads to the results shown in Fig. 2. Figure 2(a) contains time histories of all parameter sensitivity coefficients associated with population for the time interval 1973-2100, whereas, Fig. 2(b) contains time histories of all initial condition sensitivity coefficients associated with population. The interested reader is referred to [ll] for plots of the remaining sensitivity coefficients. Figure 2, together with implications of equation (7), provides additional quantitative information which can be introduced into substantive discussionsof the validity of the model which produced the results shown in Fig. 1. For example, estimates of the changes in amplitude of population which would accompany specified small changes in the components of the parameter vector can be obtained almost by inspection of Fig. 2(a).
System dynamics models of social systems
Time
Fig. 2(a). Parameter sensitivity coefficients associated with population. P.
The following remarks briefly summarize the main observations obtained from the parameter sensitivity
analysis of the Forrester model: 1. The strongest influence ofa parameter change occurs at a considerable length of time after the change, sometimes as much as 100 yr later. 2. The most influential parameter is usually, but not always, that directly associated with the state: an exceptionis theeffect ofachangeincapitalinvestment generation normal (CIGN) on natural resources and pollution. This high sensitivity to CIGN has been noted elsewhere [ 151. 3. Changingaparameter to influence its associated state variable may have considerable influence ..on other states. Note that these observations have previously been made by Forrester in a qualitative sense. The sensitivity coefficients provide a quantitative assessment of the delay times, relative influence and second-order effects. The following observations summarize the main results of the initial condition sensitivity analysis for the Forrester world model: 1. The results demonstrate that all of the states are most sensitive to the initial setting of natural resources. The influence which this initial state has is seen to grow and become increasingly important as time progresses, except in the case of natural resources itself, where the influence. is observed to decrease with time. 2. Most state variables indicate a very high degree of sensitivity to the setting of the corresponding initial
219
state for the first few years of the time period. In every instance that sensitivity is seen to decrease as time progr&ses. Sensitivity analyses contribute both to model validity and to model utility. The sensitivity analysis of the initial states contributes to model validity because of the insight obtained regarding how small variations in the initial states influence the states at future points in time. Thus, a high degree of accuracy in the initial setting of natural resources is desirable because of the dominating influence which thisinitial condition is observed to have. On the other hand, the parameter sensitivity results provide insighj as to which parameters are the most influential for controlling a given state. Thus CIGN would represent an appropriate parameter for controlling pollution and for controlling the rate of natural resource usage. In summary, sensitivity analysis represents a mechanism for measuring the sensitivity of any given input with any given output, and describes how changes in the input parameter would influence future behavior. Relative to model validity, sensitivity analysis suggests data collection priorities to improve the accuracies of the more sensitive parameters and initial conditions in the model. Relative to model utility, sensitivity analysis provides insight as to how changesina particular control parameter would influence the model state at some future time. In Fig. 2(a), increases in POLN and NRUN are seen to have a negative effect on population. The fact that an increase in NRUN decreases population is counter to one’s intuition regarding the actual cfiect an
Scale
I unit=
factor
---.-.-. -.-
1
IQ‘73
--.F’\.
1,147~ IO9 PO NR, POL, CIAF, Cl0
I 2000
I 2025
I 2050
I 2075
I 2Kx)
Time
Fig. 2(b). Initial state sensitivity c&fficients associated with population. P.
JAMESR. BURNS and DAVID W. MALONE
220 increase
in NRUN would have on population. This and all such similar observations provides additional information regarding the dynamical behavior of the model andexplainswhysensitivitycoefficientsare the quantitative representation of what Forrester refers to as “counterintuitive” system behavior.
parameter vector. Precisely stated, the problem is to determine I*@*), where ‘I I*@*) = min rl(QL, P, u, t) df WM I It, subject to: x==(x,u,t),
OPTIMIZATION
ANALYSES
In this section two methods for optimizing system dynamics models of social systems are presented, using the Forrester world model as the vehicle of demonstration. Specifically, the extensions necessary to apply (1) parameter optimization techniques and (2) the variational calculus are discussed together wirh resultant implications for behavior of the world model. Suchanalytical techniquesgenerateadditional objective informationand insights regarding model behavior, and hence contribute to substantative deliberations regarding model validity. Before any optimization analysis can be performed, a mathematical construct of the value system used to measure the performance of the model must be developed. In this study the following index- was developed to measure the performance of the model:
where r/(QL. P. u, t) = -C,
.QL2 + Cz (P - P,,&* +
c3
.IIu
-
111;.
(8)
The index is intended to mathematically characterize the valuesystemused by Forrester to evaluate the model behavior.Thefirst term to appear in the index represents a desire to maximize the quality-of-life, QJ!,, for the time interval [to, ts]. The second term represents an intent to stabilize the population, while the third term measures the amount of control effort required to accomplish the objectives represented by the first two terms. In words this index stipulates that quality of life be maximized and population be stabilized while minimizing the extent of intervention into the system required to achieve these goals. The constants Cr, C2, and C3 are used to weight each term appropriately and Q is a five-dimensional identity matrix for the model being considered. It should be emphasized that other value systems could be used to measure model performance and that other different mathematical constructs exist for monitoring model performance. The parameter optimization problem requires choosing a parameter vector u to minimize the index of performance defined by (8) subject to the system differential constraints given by (2) and certain physical limitations on the choice of values available for the * The
values for BRN
and CIGN
are normalized
values.
x(t,)=x,,
M=
:u(ui20,i=
1,5).
The Powell method [12] was selected for use in this study because it was found to provide faster convergence than other function minimization schemes, yet it does not require analytic expressions for the first partial derivatives of the index of performance. In order to prevent the method from choosing parameter values which are physically unrealizable, penalty functions were employed. The penalty function defined below was chosen from several other alternative configurations which seemed less well suited to this problem. I+l+Cgi and 0, ui 2 0 gi = 1 c. l&l. ui < 0.
C = 1 x 105.
The parameter optimization program was started with an initial estimate for the parameter vector of ur = (1, 1, 1, 1, 1). This initial estimate corresponds to the situation in which the world is permitted to run its course in the absence of restraints being imposed upon its behavior. This initial estimate produced a value for the index of performance of 12,800. In the first 120 function evaluations, the index of performance decreased from 12,800 to -40. In the last 140 function evaluations, the index decreased from - 40 to - 46. The parameter specification at the end of the 260 function evaluations was determined to be u*r = (0.60, 0.52, 1.16, 0.89, 0.71). These values are to be compared with the set of values which were hand-picked by Forrester for the explicit purpose of producing an improvement in the system behavior relative to his implicit value system. These values (see [Z], p. 120) are: uT = (0.70, 0.25, 0.50, 0.80, 0.60).* Some interesting contrasts appear between the handpicked values of Forrester and the values picked by the computer. The computer-determined value for birth rate normal is more severe than that chosen by Forrester and therefore harder to achieve. However, the restraints applied to natural resource usage would be less severe under the computer-determined policy as opposed to the Forrester policy. Interestingly enough, the computer-determined policy suggests allowing pollution generation to increase rather than restraining it as the Forrester policy recommends. The resultant plot did not indicate excessive pollution. The values chosen for agricultural production CONAF and for CIGN are less restrictive than those corresponding values determined by Forrester, and therefore easier to achieve. It is interesting to note that the parameter values picked by the parameter optimization algorithm
System dynamics models of social systems
are in the interior of the feasible region of the parameter space M, and not on the boundary. In contrast to the previous approach, the strategy of the open-loop control approach is to treat the parameters as functions of time and to select the optimal set of time functions which will minimize the index of performance.Thustheproblemistoselectanadmissable u*(r) such that ‘I I*@*) = mm
u&f s to
rr specified
and i = 1,. . . (5. to 4 t 5 cl}.
Due to the high degree of nonlinearity and coupling which resides in the system equations, considerable computational irregularities enter into the ‘problem posed above. The table functions mf(x, u, r) were fitted with polynomials so as to provide continuous derivatives. The index of performance was converted from its Lagrangian format shown above into the form of Mayer. Then the optimization problem can be reformulated as follows: Z*(u*) = min x,(t,) uelu subject to : 2 = f’(z, u, t), .z(t,,)= zo, tf specified where: zT = (XT,X6) .fT = (.fr. rl(QL, P, u, 0). Weak necessary conditions for a stationary point within the control subspace M for the problem posed above require that the first variation of the augmented index of performance vanish. The augmented index is defined as Al = x,(t,) +
” IT{& I (0,
u, r) - i) dt.
ax&,) - (;.T.l%):~,= 0. AT $z,
LT
A = ‘;(z,
af * y-&
u)
(2, u)
=
u). (13)
These adjoint (Euler-Lagrange) equations are solved simultaneously with the system equations: z(t0) = zo.
(14)
The resultant problem is generally referred to as a twopoint boundary-value problem (TPBVP). A great many methods have been proposed for solution of the TPBVP, many of which were applied without success to the problem posed above. Any method which required the simultaneous integration of both the system and adjoint equations, whether forward integration, backward integration, or both is required, was generally met with failure. These disappointing results eventually led to the selectionofthefirst-ordergradientmethod developed by Denham and Bryson, and discussed in Citron [I33 andBrysonandHo[14].Thissteepestdescent technique relaxes condition (12) while requiring that the system and adjoint equations be satisfied at all times. The relaxation removes the coupling of the control vector which would otherwise exist between the system and adjoint equations. In addition the technique permits the system equations to be. integrated forward in time while the adjoint equations are integrated backward in time, thus taking advantage of the fact that the system equations are numerically stable during forward integration while the adjoint equations are stable during backward integration. The first-order gradient method consists of the following steps: 1. Estimate a control vector u(r) and integrate the systemequations(14) to the terminal time t/, storing the state trajectory x(t). 2. Integrate the adjoint equations (13) backward and store the resultant adjoint trajectory. 3. Computeanewcontrol trajectory using the following formula :
(9)
Taking the first variation of (9), noting that the continuity of z(r) permits interchanging the operations of differentiation and variation, and taking advantage of the functional independence of r(t) and u(r), it can be shown that the weak necessary conditions for the problem considered become:
AT +
1 = -AT . L, n(QT = (0, 0, 0, 0, 0, l),
z = s’(z, u, tb
subject to
M = [ulu&) 2 0,
From (10) it is determined that n(tJ = (0, 0,0,0,0,1). Thus the adjoint system is defined by:
tl(Qk P, u, 0 dt
i = f(x, u, t), tit,,) = x0,
221
(10)
0,
(11)
=I 0.
(12)
U new =
%ld
-
wi:p
where
and Q is the (m x m) unit matrix, dS is the step size. 4. Repeat steps 1, 2 and 3 until no further improvement in the index of performance Zcan be obtained. OPTIMIZATION ANALYSIS RESULTS
The open loop optimization program implementing the method just described was started from several
JAMESR. BURNS and
222 Sea le factors
I‘-__
--
-P -----.-.-.-
-.
/
. . .’ x ./
./’ 2cI
_ ror-
.’
.’ -
2x109 NR 25x IO” POL8 x IO’O QL 5 x IQ’ Cl 5 x 109 . .._‘N. ‘. ‘. ‘\ \ ‘\ ‘\ ‘\ ‘\\ ‘\ ‘\ ‘\ . ‘\ ‘1. .-. _-‘-‘\_
.A
C‘
I’
0’
1.
./
/
--
‘\
,/
-- --173
l
2000
I
2026
I
2050
I
2075
I
2100
Time
Fig. 3(a). Model response to control trajectories pictured in Fig. 3(b).
different points in the control subspace M and allowed
to continue its iteration until very little additional improvement in the performance index was realized. The terminal results were then plotted. Figures 3 and 4 show the terminal results obtained by starting from two different points, ml and m2 in M. In Fig. 3(a) is shown the model response induced by the control setting shown in Fig. 3(b) (such control setting determined by starting the iteration at point ml), and similarly for Figs. 4(a) and 4(b). The starting point m, EM, from which Fig. 3 was obtained is uT = (0.70,0.25,0.50,0.80, 0.60) where u is assumed to be constant over the whole time interval. The attentive reader will recognize these values as the optimal parameter values determined by the parameter optimization. The starting point m2 EM from which Fig. 4 was obtained is uT = (1, 1, 1, 1, l), or the nominal, unconstrained parameter values associated with the current policies of the early 1970%. In both runs the iterative steepest descent algorithm was allowed to continue until no further reduction in the index of performance could be obtained. The performance value associatedwithFig.3(b)is - 127whereas that associated with Fig. 4(b) is -102. In both figures the control trajectories were constrained to the nominal values corresponding to the present policies for the first 3 yr of the simulation interval by weighting heavily the third term of the index of performance for this period. The multimodal response surface defined by (8) is responsible for the obvious observation that the results in Figs. 3 and 4 are different. The consequence of starting from two different points in the control subspace is that convergence to two different local minima defined by the response surface has been achieved. As
DAVII)
W.
MALONI:
a result, the analyst is faced with the problem of choosing between two different control policies, each about equally ‘good’ insofar as goodness is defined by the performance measure (8). Apparently, the system trajectories must be evaluated for desirability, the control trajectories for feasibility. The system trajectories shown in Fig. 3 are more desirable than those shown in Fig. 4 because quality of life is higher and natural resource usage is lower. However, the control trajectories shown in Fig. 3 are less desirable than those shown in Fig. 4 because of the increased amount of control effort required to achieve the higher quality of life. Since the performance indexconsiders both of these factors (quality of life and control effort), it ‘measures’ the performance observed in Figs. 3 and 4 as about equal-a minus 127 for Fig. 3 compared with a minus 102 for Fig. 4. These numbers, as compared with the performance value of 12,800 associated with Fig. 1. are roughly the same. The continuous, time-varying control strategies suggested here appear to be much more realistic than discrete interventions into societal systems because discrete changes, especially if they are dramatic, are psychologically resisted. The implication is that continuous time-varying control of social systems is more reasonable and less restrictive than one is led to believe from the results of Forrester’s analyses. In particular, Fig. 4(b) suggests a time-varying control strategy which manipulates or constrains only two of the five control variables. The remaining three control parameters retain their nominal values of 1.
LcOLN______
/ -___
,CONAF
\y..-._.-._._L._.-.-_.--.-.-.-.__
0.90
0.30
,-
i
C ‘19 73
I
2000
I
2026
I
2050
I
2075
I
2100
Time
Fig. 3(b). Optimal control trajectories obtained by starting from point I in the control suhspncc M.
System dynamics models of social svstems
Scale A
_--0
1973
factors
I.20
P 2 x 10s NR 2.5~ IO”
----
3 A\
t
_-----__ I 2000
I 2025
223
0.30 I 2030
I 2073
I 2100
Time
Fig. 4(a). Model response induced by control shown in Fig. 4(b).
t
histories
0~0
Time
Fig. 4(b). Optimal control trajectories determined by starting from point 2 in the control subspace M.
CONCLUSION
This paper presents several algorithms for analyzing and optimizing system dynamics models of social systems. The usefulness of these generalized methods is demonstrated via their applicatian to Forrester’s moclel of the world. The need for such a demonstration was
discussed in the introduction. There attention was directed to the recognition amongst statesmen and scientists alike of the facility which system analysis has to offer to the study and analysis of social systems. In this framework, many such problems can be viewed asarisingfrom the inability to synthesize and implement policieswhichwouldsteerthesesystemsfromthepresent undesirable states to more ideal future states. The methodology suggested in this paper addresses this problem by demonstrating the use of automated techniques for policy determination. In addition, a review of previous related studies showed that computational algorithms have never been applied to the controversial system dynamics models Cl-33 of today for the purpose of contributing to deliberations regarding model validity. Again, the methodology suggested here is intended to address the urgent requirement for mechanized techniques which would contribute to considerations of model validity and enhance model understanding.
REFERENCES 1.
J. W. Forrester, ,._,,.I
(IYOYJ.
UrbanDynamics.
MIT Press, Cambridge
2. J. W. Forrester, World Dytuxmics. Wright-Allen Press. Cambridge (1971). 3. D. H. Meadows. D. L. Meadows (‘I ~1.. T/K* Limits to Growrh. Universe Books New York (1972). 4. R. Boyd, World dynamics: a note. Sci~ncu. 177, 516-519 (1972). 5. C. Starr and R. Rudman, Parameters of technological growth, Astron. Aeron. 20-27 (1973). 6. E. S: Boylan, The system dynamics approach to modeling world-wide interactions:8 critical analysis. Unpublished Report, Department ofMathematics, Rutgers University (1972). 7. H. S. D. Cole. C. Freeman et al.. Models o/ Doom Universe Books, New York (1972). 8. M. H. Rothkopf. World models won’t work. Simtrltrtio~~ -61 9.
ilY73).
A. P. Sage, Ed., SMC-2
10.
IEEE Trans. (2), 122-220 (1972).
Systems, Man,
R. Tomovic and M. Vukobratovic,
General
Cybernetics Sensitivity
Elsevier, New York (1972). 11. J. Burns, Applications of control theory to system dynamics models of social systems, Ph.D. Dissertation, Purdue University (1973). 12. M. J. D. Powell, An efficient method for finding the tiinimum of a function of several variables without calculating derivatives, Computer J. 7, 155-162 (1964). 13. S. J. Citron, Elements of Optimal Control. Holt, Rinehart and Winston, New York (1969). 14. A. E: Bryson, Jr. and Y. C. Ho. Applied Optimd Cortfrrd Blaisdell, Waltham (1969). 15. R.A. Burnett and P. J. Dionne, GLOBE6: a multircgion interactive world simulation, Simuhtion 192. 197 ( 1973). Theory.