4 MODELING, THE MECHANICS

So far the discussion has been largely qualitative, stressing orientation in problems and systems. This chapter develops basic quantitative tools to model systems. The System Dynamics paradigm is developed, and the DYNAMO programming language is introduced to quantify the models.

4.1 The Time Step and its Notation

Let us first clear up an important administrative matter in simulation, that of the time step. We have discussed the concept of an "iteration" to obtain values for a model's variables from one instant (T--T) to another (T) and how a series of iterations provides model behavior over time. This notation can be simplified to avoid subscripts or parentheses. The following convention, used in system dynamics modeling, does so.

The value of a variable VAR at time T, previously specified VAR(T), will be designated VAR.K. Similarly, VAR(T--T) will be denoted VAR.J, VAR(T+-T) will be VAR.L. A rate variable FLOW between the time points T--T and T will be denoted FLOW.JK, and FLOW between T and T+-T will be FLOW.KL.

The following tables show this symbolically.

Variable VAR at a time point:

Old notation: VAR(T--T) VAR(T) VAR(T+-T)

New notation: VAR.J VAR.K VAR.L

Variable FLOW between time points:

Old notation: FLOW(T--T,T) FLOW(T,T+-T)

New notation: FLOW.JK FLOW.KL

This "JKL" notation, which is easier to type for computer entry, will be used henceforth.

The intervals JK and KL each equal one time step, that is, a -T and not a time unit. A model is simulated by calculating the model variables each time step. Furthermore the time point K represents the current time, while J represents the last calculation time, and L the next. In other words, K is "now", J and JK are the past, and KL and L the future. An iteration causes the JKL time sequence to move with each time step:

a selected iteration J K L

one iteration later J K L

two iterations later J K L

--------------> TIME

A point concerning rates may avoid confusion. A rate of flow, in a simulation, is held constant for one -T, that is during the time step -T between model iterations. Yet the definition of a rate of flow -- that is the term "rate" -- is the rate per unit of time. If the time units used in the model are months, say, then the value for a variable's rate of flow is the flow per month, not the flow per -T. Yet the value changes at each iteration. On June 10th, the indicated monthly rate of flow of sales may, that day, be 230, on June 11th it may be 310, and so forth.

Since the rate of flow in the real process being modeled may change much more often than monthly, yet the basic time unit used in the model remains a month, it follows that the simulation must be iterated more often than monthly. This is done by making -T smaller than 1.0. A -T of .25 would calculate the model four times each month.

Now, since the rate of flow is defined as flow per time unit (month), but the calculation is iterated four times per time unit, then a rate of flow FLOW.JK would need to be multiplied by the value of the -T to determine the actual flow occurring during the -T between J and K.

We therefore multiply rates of flow by the value -T when they are accumulated in a level, as we saw in the previous section.

Finally, the value of the rate of flow held throughout the interval JK is the rate calculated at time point J, and the value of the rate of flow held throughout the interval KL is the rate calculated at time point K. This is relevant mathematically, as it means that a rate of flow, instead of being defined dynamically as a derivative of other variables, is instead an analytic function of those variables. On the other hand levels are defined as integral equations, incorporating the dynamic aspects of the model.

4.2 Levels and Rates, the Conservative Components

Levels and rates are the system components which contain the system's material -- as opposed to informational -- flows.

The value of a level at time point K will always be defined as the same level at time point J plus the inflows and less the outflows between J and K. The following sample DYNAMO equation applies:

L LEVEL.K=LEVEL.J+DT*(INFLOW.JK-OUTFLOW.JK) (4.2.1)

The initial L alerts the DYNAMO compiler to the fact that this is a level equation, causing the compiler to check the equation's format. This check enhances program debugging. Note that to be consistent with DYNAMO notation, -T has been replaced by DT. This is more compatible with standard keyboards, and we will use DT instead of -T henceforth.

A level can only be a function of its own previous value, and the rates flowing into and out of the level. It cannot be a direct function of another level. This is consistent with the fact that a level is only directly affected by rates flowing into or out of it, not by other levels.

There may be more than one inflow or outflow, or there may be no inflow, or no outflow. Note that a degenerative case of a level equation would be one which had neither inflows or outflows. In this event the level would be constant, equaling its last value plus the net zero flows.

Rates are the flows into or out of a level. Rate equations take many functional forms, depending on the system being modeled, but they always begin the same way. The equation for a variable named FLOW would begin

R FLOW.KL=

with the right hand side being a function of other model variables.

The R alerts the DYNAMO compiler to a rate equation, so the compiler can conduct consistency checks.

Rates provide values to be used in the future (the next DT), while levels are calculated based on past values. In fact, a rate determined for the next time step, FLOW.KL, will, in the next iteration, become the past rate FLOW.JK for use in calculating the level it effects. This convention prevents simultaneity in in the equations defining the system variables.

Before providing an example of a system of levels and rates, it will be useful to introduce the graphical symbology of system dynamics. Levels are symbolized as rectangles and rates as butterfly shaped "valves". Figure 4.2.1 applies.

The cloud shaped symbols represent the environment outside the system. The solid lines with directional arrows indicate conservative flows of material through the valves, and into or out of the levels. LEV1 has two inflows (RATE1 and RATE3), and one outflow (RATE2), while LEV2 has one inflow and two outflows.

The level equations for figure 4.2.1 would be written as follows:

L LEV1.K=LEV1.J+DT*(RATE1.JK+RATE3.JK-RATE2.JK)

L LEV2.K=LEV2.J+DT*(RATE2.JK-RATE3.JK-RATE4.JK)

Level equations can always be defined once the inflow and outflow rates are known. The rate equations for figure 4.2.1 cannot be written, however, because rates depend on informational elements, auxiliaries, and constants not shown.

4.3 Auxiliaries and Constants

Flows within a system are analytic (static) functions of other variables and parameters. They can be functions of the system's levels, or other rates, or information auxiliaries. Auxiliaries can, in turn, depend on variables in the system, (levels, rates, or other auxiliaries) and/or on parameters from outside the system. We need symbols for auxiliaries, and parameters. Figure 4.3.1 expands 4.2.1 to include both.

Auxiliaries that are determined endogenously are drawn as circles, like AUX2. If exogenously provided as functions only of time, they are drawn as AUX1 is. Parameters, such as K2 and K4, are shown as symbols representing the Greek uppercase letter theta.

The dotted lines indicate information flows, in contrast to the solid lines showing material flows. Figure 4.3.1 shows most of the basic symbols used in system dynamics modeling.

In this figure, the levels are still defined as in figure 4.2.1. Level equations always can be defined simply by knowing the system configuration. Rate and auxiliary equations, however, are not specifically known from merely seeing the system diagram. This is consistent with reality. Looking at a plumbing system does not tell the observer what is flowing through the valves or how they are controlled, yet any accumulation -- such as a bathtub filling -- is known to depend on its previous value, plus inflows less outflows.

We can, however, derive functional dependency from the diagram. For example, we know that RATE2 is a function of LEV1 and the constant K2, as it has information links (the dotted lines) to those elements. RATE3 is dependent on RATE2, AUX2 and K3, while AUX2 is dependent on LEV2 and RATE4. To demonstrate some specific functional relationships for a rate and an auxiliary, let us assume in the system modeled that RATE3 equals one-half of the sum of AUX2 and RATE2. Let us also assume that AUX2 is equal to RATE4 plus RATE4 divided by LEV2. Equations used in the DYNAMO program would then read,

R RATE3.KL=K3*(AUX2.K+RATE2.KL)

C K3=0.5

A AUX2.K=RATE4.KL+(RATE4.KL/LEV2.K)

Exactly why these equations hold is not important here, we are only interested in format so that we can begin to model. Note that the first single letters alert the DYNAMO compiler to the type of equation: R for Rate, C for Constant, and A for Auxiliary. Both rate and auxiliary equations depend only on values known for the time point K (recall that rates for the interval KL are determined at time K, and then hold until recalculated one iteration later, when the previous L has become the new K.)

4.4 Initial Values, Model Specifications, and Labels

To simulate a system model, the model equations must be initialized. That is, the initial values for the levels of the system must be given. This merely means that each level equation must be augmented by an initial value equation, called an "N" equation:

L LEV1.K=LEV1.J+DT*(RATE1.JK-RATE2.JK)

N LEV1=100

Additionally, a DYNAMO program requires several specifications, including the value for a DT, the length of time the simulation is to cover, and the variables to be saved and printed out as tables or graphs. A typical specification line, placed at the end of a program, is the following:

SPEC DT=.25/LENGTH=100/SAVPER=5

SAVE LEV1/RATE/AUX2

DT=.25 causes an iteration to occur four times per time unit (time units are inherent in the model, they may be minutes, days, years, decades, etc), and LENGTH=100 will cause the simulation to run for 100 time units. Note this means there will be 400 iterations. SAVPER=5 means the value of each variable will be saved each fifth time unit (at time 0, 5, 10, etc), and the variables to be saved for output are LEV1, RATE1, and AUX2. Any number of variables can be saved, and during debugging runs it is wise to save many variables to provide logic checks.

Finally, we can partially document a model by placing comments after one or more spaces at the end of an equation. For example, if a rate of flow equaling the cash outflow for a proprietor's weekly labor costs is called LBUD, the equation with comment might be written:

A LBUD.K=WKFRCE.K*WAGE.K LABOR BUDGET ($/WK)

where the comment reminds the modeler or user what LBUD represents, and what its units of measure are.

It is not the intent here to be complete in explaining the DYNAMO language. Program manuals will provide details. But with the above introduction, we can begin to quantify some sample system models in chapter 5. Then, through the examples and exercises in chapter 6, the reader should become reasonably familiar with modeling system dynamics models using DYNAMO.

4.5 Logic and Test Functions

Modeling often requires observing constraints and/or discontinuities in functional relationships. For example, inventory levels cannot be negative, so a sales rate may be constrained to be no larger than remaining inventory. Or a 20 percent price reduction may be imposed suddenly when inventory reaches a level considered excessive. Or a variable may be constrained to lie between two limits. Such functional relationships can be called logic functions, in that they imply an "if-then" logic: "If inventory exceeds 100, then decrease price by 20 percent".

There are various forms of logic functions in simulation languages. The DYNAMO documentation provides details, but a few sample functions should prove useful.

One useful function is known as a "clip" function. It has four arguments, and takes the following form:

A VAR.K=CLIP(A,B,C,D). (4.5.1)

The arguments are related in the following way. If C is greater than or equal to D, then the variable VAR equals A, otherwise VAR equals B. The arguments A, B, C, and D can be variables or parameters.

For example, actual sales AS can be made the lesser of sales demand SD or inventory INV by the function

A AS.K=CLIP(INV.K,SD.K,SD.K,INV.K). (4.5.2)

This is read "if sales demand (SD) exceeds inventory (INV), then actual sales (AS) equals inventory, otherwise actual sales equals sales demand (SD).

To reduce price to 80 percent of its original value PRICE when inventory INV exceeds 100, one could define actual price APRICE as follows:

A APRICE.K=CLIP(PRICE.K,.8*PRICE.K,100,INV.K). (4.5.3)

CLIP's can be embedded within CLIPs. If price is reduced to 80 percent of its original value when inventory reaches 100, then remain fixed until inventory reaches 130, when price is reduced to 50 percent of original value, one could define an embedded clip EMCLIP as follows:

A APRICE.K=CLIP(PRICE.K,EMCLIP.K,100,INV.K) (4.5.4)

A EMCLIP.K=CLIP(.8*PRICE.K,.5*PRICE.K,130,INV.K) (4.5.5)

Two other useful logic functions are "MAX" and "MIN". MAX(A,B) takes the value of A or B that is greater, while MIN(A,B) takes the value that is lesser. To ensure sales are positive, on could define actual sales ASALES as:

A ASALES=MAX(0,CLIP(INV.K,SD.K,SD.K,INV.K)), (4.5.7)

where the embedded CLIP is expression 4.5.2 above.

Test functions are of several forms, including pulses, steps, sine waves, and random noise generators. Again, the DYNAMO documentation should be consulted. However two functions are particularly useful in testing the response of a system in steady state to either a sudden step increase in a variable, or a pulsed change to a variable.

The step increase is written:

A NVAR.K=VAR.K+STEP(A,B) (4.5.8)

Here, the variable VAR is changed by adding a step increase at time B, with the height of the step being A. Assume a system has reached steady state by time 20, with VAR equal to 12. The expression

A NVAR.K=VAR.K+STEP(1,23) (4.5.9)

would result in a variable NVAR which had value VAR up to time 23, but NVAR would equal VAR+1 starting at time 23. This would presumably cause dynamic reactions which effected all variables thereafter, (including VAR), and the system could be observed to determine the transient effects of such a step change in VAR.

PULSE has the form:

PULSE(A,B,C) (4.5.10)

where a momentary jolt of height A occurs at time B, and the jolt reoccurs every C time units thereafter. Two things should be noted about PULSE. First, the pulse height remains at the value A for only one DT, not for one time unit. So if DT is .25, and the total cumulative input jolt desired is one unit, then the pulse height A should be 4, as 4*.25 = 1.0. Second, if the modeler wants only one pulse to occur during the simulation, a large value of C should be chosen, one greater than the length of the simulation run.

This cursory introduction to logic and test functions is meant to show that such functions are available. The user will need to consult program documentation for a more complete listing of such functions, and practice their use to gain proficiency. The examples of chapters 5 and 6 provide exercises incorporating selected functions.

4.6 Delays in DYNAMO

DYNAMO provides several functions that cause variables to be delayed. Delays were discussed qualitatively in Section 3.6. Here the concern will be mainly with symbols and equation forms used to model delays.

Delays are used to delay the effect of either a rate or an auxiliary. Delays in DYNAMO are programmed "macros", lines of code that allow the modeler to use one expression to replace a series of longer operations. In this case, a delay function replaces one or more sets of level/rate pairs. For example, a "DELAY3" -- which is a third order delay of a rate -- takes the place of a string of three sets of levels and rates. In the following figure, the function

A DFLOW.K=DELAY3(FLOW.KL,9) (4.6.1)

shows the macro behind DELAY3, which delays the variable FLOW an average of 9 time units (not 9 DT's).

There are two basic types of delay macros. The first type, like DELAY3, is a material delay, used to delay rates of flow. "DELAY1" is a first order delay, and DELAYP, a "pipeline" third order delay which allows access to the sum of the materials held in the three levels of the delay. The pipeline delay is handy if, for example, one is calculating the fraction of a total population in each of several states -- such as in an epidemic model where the fractions of the population incubating may be needed. The incubating population may be held in a third order pipeline delay, and DELAYP would allow accessing the value of that population.

The second type of delay is an information delay. These delays are used to delay auxiliaries instead of rates. In DYNAMO, the first order information delay is called "SMOOTH", and the third order is called "DLINF3".

The form of each delay equation except DELAYP is as in 4.6.1, for example,

A DD.K=DLINF3(AUX.K,DTIME) (4.6.2)

makes DD a third order information delay of the auxiliary AUX, with the overall average delay time equal to DTIME, and

A SS.K=SMOOTH(AUX.K,DTIME) (4.6.3)

makes SS a first order delay, or "smoothing" of AUX, with the average overall delay time again being DTIME.

The pipeline delay has three arguments:

A PP.K=DELAYP(FLOW.KL,DTIME.K,SUM,K) (4.6.4)

makes PP a third order delay of the rate FLOW, delayed an average of DTIME.K (the ".K" on DTIME shows that the delay time need not be a constant), and the cumulative material in the pipeline is given the name SUM, and can be used elsewhere in the model.

Note that material delays apply to rates, and information delays to auxiliaries. What of levels? Levels cannot, strictly speaking, be delayed as material, they can only be delayed as information. For example, management may have a lagged perception of inventory size (a level) and may then act to change inventories only after a delay. But it is not the inventory itself being delayed for use later, but rather information about the inventory. The modeler, then, should set an information auxiliary equal to the level, and then delay the auxiliary.

For programming purposes, we will need new symbols for delayed variables, and also recognize a slight change in the equation format allowed for levels. For our purposes, delayed variables can be symbolized as follows.

The following figures show a 3rd order delay of a rate and of an auxiliary, each with time constant 5, and each with sample equations of how they might be used in a DYNAMO program.

Note particularly that the equation for the level "OLDER" has the term DFLOW.J as an argument. Previously we said all level equations have a form using only rates within the parentheses. The use of a macro has changed that. In principle, however, the delay is a rate that ultimately changes the level OLDER.

The lower part of the figure shows how the order rate (ORDER) may depend on a lagged inventory measure. Here, ORDER is the difference between desired inventories (DESINV) and management perception of inventory (PINV), which is a third order delay of actual inventory, the delay being 5 time units.

Finally, section 3.6 discussed the boxcar delay. No standard macro exists for this delay. The modeler could, however, build one by simply making each "boxcar" a level, the flow rates out of each level being equal to the level itself, and making DT=1.