Algorithmic styles for nonlinear function generation in lumped-parameter systems

TECHNICAL NOTE

Algorithmic styles for nonlinear function generation in lumped-parameter systems W. E. Cundiff School of Administration, Griffith University, Nathan, Brisbane, Queensland, 4111, Australia (Received Februao' 1987: revised October 1987)

Within the simulation community, it is often believed that integration is the limiting factor to program performance and run-time efficiency when simulating large-scale continuous feedback systems. Much attention is given to selection of the best approximation technique (e.g., Euler, RungeKutta, Adams-Bashforth). However, performance is shown to be up to six times more sensitive to impact on nonlinear function generation due to computational complexities associated with one- and two-dimensional table-look-up algorithms. Such algorithms utilize linear interpolation and sequential search that when combined, produce a number of approaches to algorithm construction. Algorithms and associated data structures have been formulated in contrasting general-purpose languages to illustrate stylistic differences and for direct application in nonlinear simulation of lumped-parameter systems.

Keywords: lumped-parameter systems, algorithm construction, programming style, system dynamics, table-look-up

Introduction Earlier work has focused on the embodiment of dynamic feedback models in the notation of APL and their further simulation directly via an APL interpreter, t These investigations prompted extended examination of other forms of simulation 2 and interest in the innerworkings of dynamic feedback models. Of special relevance was the aspect of finite differences and of integration, central to all continuous systems simulation, in particular that adhering to the system dynamics methodology? The application of simulation techniques in problem solving or, as a means toward better understanding of phenomena, is usually brought about in the absence of analytical solutions, relying on an established, elegant mathematical foundation. Nonlinearity generally characterizes candidate systems for computer simulation. As such various means for introducing nonlinear behavior into models via multiplier coefficients have been developed, most notably the table function. Function specification by way of a two-dimensional data structure allows great flexibility in setting model

340 Appl. Math. Modelling, 1988, Vol. 12, June

coefficients. These act as multipliers of rates affecting state variables dependent on the value of an independent variable at an advancing time slice or fixed time increment during the course of a simulation. Assemblages characterized by ordinary differential or difference equations and discrete state variables have been described as lumped-parameter systems.4 It is this class of system that is central to the majority of applications in environmental and engineering simulation, especially where feedback is present. These systems contrast with so-called distributed-parameter systems embodied in partial differential equations. Numerous alternative approaches, dictated by either availability or personal choice, have been applied to simulation studies and the documentation of methodology for modeling feedback systems. Subroutine repertoires for popular third generation languages, including FORTRAN, -~BAS1C, 6 and Pascal, 7 have been developed in the absence of, say, DYNAMO or CSMP. However, the aspect of nonlinear function generation, involving sequential search and interpolation for table look-up, is often either not discussed or is not sup-

© 1988 BunerworthPublishers

Technical note

ported by the software package or compiler. "9 This includes not only separately compiled libraries for general-purpose languages, but special-purpose languages and compilers for continuous systems simulation.

Model integrity Recently, increased attention has been extended to formal procedures for model specification. With origins in software engineering, the integrity of large-scale models should be subject to rigorous evaluation of program reliability, correctness, opportunities for systematic enhancement, and overall disciplined approaches to the management of complexity. ~° Frequently, discussions of model verification and validation appear in the literature, but such considerations are somewhat limited in comparison to analogues applied to more prosaic production level applications, e.g., management information systems for planning and control in organizations. Because of the advent of high-level, enduser simulation software, the adherence to any formalism is often bypassed in the pursuit of results. Nowhere is this more evident than in the use of powerful microcomputer-based software, in particular spreadsheets and logic-oriented financial modeling packages. So, while many simulation studies are based around languages that are outdated in terms of internal execution efficiencies, syntax, and conventions for structured program design, algorithms, along with associated comments and program documentation, can be obscure to even the program's creator.

Languages and style Given the variations in program design that are possible, noting the fundamental constructs of sequence, decision, and repetition, consistent style is of prime importance in building large-scale models and the subsequent ongoing maintenance effort in support of modelbased systems. The discussion that follows examines varied algorithmic styles via their programming in five diverse languages, the aim of which is to provide an addition to the programming toolbox of the modeler for deriving resultant multipliers through computation of nonlinear functions in table form. The elaboration of the five modules is not intended to be in any way comparative. Rather the focus is on style as expressed, albeit subject to the limitations or robustness of particular languages, in the individual algorithms. Before beginning the discussion, we first direct the reader to Table I, which documents the central importance of table look-up construction in continuous feedback models of the lumped-parameter class. The performance figures are from a simulation carried out on a Motorola 68000 CPU running CSSL IV. The model involves six state equations, and the researchers note: "As in many continuous system simulation programs, non-linear function generation dominates the execution time of the entire program. ''~ Second, the reader acquainted with the aforementioned work by the author may find avenues for in-

Computation performance

Table I *

Operation

Machine cycles

Addition Subtraction Multiplication Division Square Square root Integration (third-order Adams-Bashforth) Real pole (integration) Function of one variable table look-up Function of t w o variables table look-up

6 6 72 160 70 330 250 450 450 1500

*Adapted from Ref. 11.

corporating the algorithms so described in models of their design and construction. Depending on the locale, the availability of an APL interpreter may prevent or limit building models in this language. As such, the discussion is very much concerned with the practical aspect of expanding the work to date in APL to other general-purpose languages, especially those available in full-scale release for installation on microcomputers, including Apple II 8-bit processors and IBM-XT or upwardly compatible "work-alikes" incorporating 16bit chip technology. We now describe the five cases individually, drawing attention to language-specific characteristics and noting facets of their construction within the general context of programming style and practice.

Module 1 The discussion begins with APL, as noted, the original focus of work leading to the present topic. All examples, except for Module 5 are based on a 2 x N array supplied as a parameter or argument to the particular module. The fundamental mathematics underpinning linear interpolation is well known and found in many basic texts on statistics. For a very good description in the context of APL and programming, see Ref. 12. Computation begins with a sequential search along the X dimension (visualize the 2 rows of the data structure as individual X and Y axes as in Figure I).

/ / / ¢vLtJ 3-

F--

-~ D 5--

/ %.

d

EL

/

\

./

>{

/

\

/

\

/

\ 2-

/

\ \

/

1

2sbo saba zsbo 10doo FOOD

o

Figure I

Appl. Math. Modelling, 1988, Vol. 12, June 341

Technical note

MODULE 1 - A P L

[o] [1] [2] [3] [4] [5] [6]

.DL R V M W;CK;S;XY;MI;Z .GO ((.RO V[2;]) < C K V[2;].IO W)/LO R V[1;CK];.GO0 LO:S (Z Z[.GD Z V[2;],W]).IOW XY V[;V[2;].IO Z [ ( S - 1),S+ 1]] MI (XY[ 1;2] - XY[ 1;1 ])%XY[2 ;2] - XY[2; 1] R ( ( M I # W ) - M I # X Y [ 2 ; I ] ) + XY[I;I] .DL

In the APL function shown in module 1, two arguments are supplied to the function M. Internally V and W are local to M, where V is the matrix of coordinates as already described. W is locally the value of the independent variable as passed to M following computation elsewhere in the body of the model. In line 1, reading from right to left, the independent variable W is tested for presence in the second row of the array V, which has as its elements the X coordinates. If not found, the built-in function .IO (index of) will return a value equal to 1 + the number of elements contained along the first dimension of the array. Line 1 seeks equivalence between the independent variable and one of the values constituting the X vector. The built-in function " s h a p e " or .RO returns the number of elements in the second row of the matrix, again the X vector. The parenthetical expression will yield a Boolean resultant that, if 1, will evoke a conditional branch to line 3, labeled LO:. If a precise oneto-one mapping is found between W and an element in the X vector, the computation proceeds to line 2, where row/column correspondence with the first row sets the explicit result R to the Y value. The function ends computation at this point, and an unconditional branch to zero, i.e., .GO 0, effects an exit.

If the branch to line 3 is effected, the local variable Z is assigned the values of the X vector along with the concatenation of W, the independent variable at time t. The lengthened X vector is sorted in descending order by evoking the grade-down function (.GD). The sequential search along the X vector proceeds once more to determine the upper and lower ranges of values within which W is located. These are computed from S, as assigned in line 3, where S is the position of W in the sorted vector. Thus, S - 1 and S + I are the index values on either side of the position occupied by W. Once these two values are known, the standard formula for linear interpolation is applied in lines 5 and 6, yielding as an explicit result, R, the Y value (multiplier) corresponding to the X value as supplied to the function locally as W. Module 2

The A P L function just described is structurally based on sequential statement execution combining both conditional and unconditional branching for program exit, as well as normal exit following computation of the final line of code. The approach involves no iteration in contrast to a formulation in a language that embraces

MODULE 2 - PASCAL

[1] [2] [3]

[4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

342

F U N C T I O N M(V:MATRIX;W:REAL) REAL; TYPE MATRIX = ARRAY [1..2,1..6] OF REAL; VAR MI:REAL; NOTFOUND:BOOLEAN; BEGIN NOTFOUND: = TRUE; CK: = I; W H I L E N O T F O U N D AND (CK<6) DO BEGIN IF (W>MIN(V)) AND (W
Appl. Math. Modelling, 1988, Vol. 12, June

Technical note

the philosophy of structured design and programming such as Modula-2, Ada or, as portrayed here, Pascal. Strongly typed languages have gained considerable currency among the simulation community, and the aspect of single entry/exit promotes modularity. Distinguishing features include the integrative processing initiated in line 10 after declarations in lines 2-9. The leading-test loop involving the WHILE construct is controlled by Boolean tests involving both a counter and a sentinel, CK and NOTFOUND, respectively. Where the previous module relied on functions as part of the set of primitives built into the language, the Pascal module requires coding two functions in the application program embodying the model or calling them from a precompiled library. Within the loop, statement 12 invokes the two functions for progressively narrowing the range; the MIN/MAX functions apply only to the second row of the martix of XY coordinates. When the range is found, by satisfying both conditions of the I F . . . T H E N . . . ELSE construct, control is passed to the compound statement sandwiched by BEGIN and END at lines 13 and 17. The standard formula for linear interpolation is again used to compute the result, which is passed to the name defining the function, in this case M, after which a normal exit to END will occur. If the logical conditions are unmet through progressive iterations and the consequent incrementing of the counter, the value of the independent variable will not fall within the range sought. Under the unmet conditions, control is always transferred to the ELSE in line 18. However the Boolean sentinel set to true in line 8 will remain so through the final pass when CK will be valued at 5. Both conditions necessary for the assignment of the value of the input parameter to the function

M will be met and the function will halt execution via a normal exit from the loop to END. Module 3

The module written in BASIC combines both conditional and unconditional branching, as seen in the initial module in APL, with the structured convention of the I F . . . T H E N . . . ELSE central to the decision logic of the Pascal function. The approach also contrasts with the previous two modules in that the subroutine for computing the nonlinear function is integral to the main program body, where calls to the module are effected by branching statements (GOSUB). The construction of a model would include the statements in contrast to APL and, to a lesser extent, Pascal. An array is introduced in the BASIC example for storing initial values in contrast to defining global and local variables independently. Another aspect of style involves direct tests for maximum and minimum values across the sequential search of the X vector of independent variables. Although the overall logic was written as closely as possible to that of the Pascal function, a comparison shows considerably more complexity in the algorithm itself when written in BASIC. Because of often unacknowledged ancestral linkages to earlier versions of FORTRAN, the programming style is dictated as somewhat monolithic in comparison to other languages, past and present: hence, the grouping within the algorithm of relational comparisons along the second row of the coordinate array prior to tests for the range within which the independent variable resides. However, like APL, BASIC, in its interpretive form, provides an interactive programming environment conducive to algorithm construction and model building that is more responsive than compiled versions.

MODULE 3 - BASIC

[32Ol NOTFOUND = 1:CK = 1 [330] WHILE NOTFOUND AND C K < 6 [340] IF V[R2,CK]V[R2,CK + 1] THEN MAX = V[R2,CK] ELSE MAX = V[R2,CK + I] IF NOT (Z[4]>MIN AND Z[4]
MODULE 4 - PL/PC

[1] [2] [3] [4] [5] [6] [7]

OPERATOR M(V M) VAR CK MI BEGIN CK=(1 INDEX(W< =V[1;]))- 1; MI = V[0;CK + 1] - V[0;CK]%(V[1 ;CK + 1] - V[1 ;CK]); RETURN(MI*W - (MI*V[1 ;CK])) + V[0;CK1 END

Appl. Math. Modelling, 1988, Vol. 12, June

343

Technical note Module 4

PL/PC is a new language and programming environment with origins in two languages: APL and Modula2. There are over 130 built-in primitive functions and subroutines that perform tasks found in enhanced APL interpreters. The control constructs are derived from Modula-2, adhering closely to the conventions of topdown design and structured programming. By the pairing of syntactical characteristics from the two languages, procedures and operators can be defined by the user that are concise and expressive. The built-in functions are named keywords, rather than relying on a specialized symbol set as in APL or, alternatively, mnemonics such as are supported by DEC's APLSF as portrayed in Module I. In contrast to APL, PL/PC executes in the conventional left to right mode. Operators are either binary or unary in terms of input parameters, much like Pascal, with procedures permitting multiple parameters beyond the one or two limited to a defined operator. Line 1 indicates two input parameters as seen in the other examples so far. Local variables are again declared similar to that of popular structured languages; note that PL/PC supports only one conceptual numeric type. Of particular interest here is the use of the primitive INDEX in line 4. Within the inner parenthetical expression, the independent variable is tested for existence in the second row of a matrix, again identical in rank and dimension of previously described data structures. The module as described uses the primitive INDEX to find the first I in the Boolean vector passed by the logical test for W being greater than or equal to each element in the X vector (note that in PL/PC the index origin is set to zero). Once found, the local variable CK is assigned the value of the element occupying the lower end of the range within which W falls by subtracting I. From here on, CK is used in the standard

interpolation formula as previously shown returning the resultant multiplier at line 6. Module 5

The final module for consideration is written in a highlevel language with numerous built-in features for simulating discrete-event systems. Extensions to the SIMSCRIPT compiler have been implemented to also simulate continuous or combined discrete/continuous systems. ~3It is with these extensions that an algorithm for nonlinear function generation would be necessary to model systems requiring either the DYNAMO TABLE/TABHL functions or, in the case of CSMP, the arbitrary function generator (AFGEN). SIMSCRIPT can be described at several levels ranging from a BASIClike language through to scientifically oriented capabilities found in FORTRAN. Beyond this, SIMSCRIPT supports ALGOL-like features, list processing, and simulation-specific routines. So far, our discussion has focused on a single, archetypal data structure in the form of a matrix of X and Y coordinates. The SIMSCRIPT module departs from this orientation, whereby the input data is in the form of a single vector which will be called XY. The vector of independent values for X is constructed as follows: XY[0] = (XY[2] - X Y [ I ] ) / XY[3] where X Y [ I ] = lower end of the range of X values XY[2] = upper end of the range of X values XY[3] = the size of the interval in the range The remaining elements are the dependent Y or fiX) values, which must be of equal length to the number of X values as calculated above. They then form pairs of coordinates for the table describing the function. In the header statement the independent variable X and the entire vector described above are passed to the

M O D U L E 5 - S I M S C R I P T 11.5

[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19]

344

ROUTINE TO MODIFY.RATE GIVEN X AND XY YIELDING M DEFINE I,J AND K AS INTEGER VARIABLES DEFINE XY AS A REAL I-DIMENSIONAL ARRAY DEFINE X,M AND D AS REAL VARIABLES IF X LE XY(1) LET M = XY(4) ELSE IF X GE XY(2) LET M = XY(TRUNC.F(XY(0) + 4)) ELSE LET I = TRUNC. F((X - XY(I))/XY(3)) LET J = I + 4 LET K = J + i LET D= X - X Y ( I ) - XY(3)*I LET M = XY(J) + D*(XY(K) - XY(J))/XY(3) ALWAYS ALWAYS RETURN END

Appl. Math. Modelling, 1988, Vol. 12, June

Technical note

routine as parameters, noting that the interpolated resultant will be returned under the variable name M. In CSMP, the syntax would appear as FUNCTION XY = (Xt,Y2),(X2,Y2) . . . . . (X,,,Y,,) M = AFGEN(XY,X) or, in DYNAMO; A

M.K = TABLE(Y, XI,X,,,Xint . . . . .

T

Y=YI/Y2/.-.Y,,

I)

Returning to the module, variables are defined by type in statements 2-4. The SIMSCRIPT decision logic of I F . . . ELSE . . . ALWAYS determines whether the independent variable is at the low or high end of the X range; if it falls outside the range, the Y value is set to one of the extremes in the range. Control is passed to succeeding statements if neither condition is met, starting with line ! 1. First, with the knowledge that the independent variable (X) is within the range of X,_,, values, the difference between X and the lower end of the range of X values is calculated, then divided by the interval of the range. The subscript value J is incremented by adding 4 to the result of the previous computation. This indexes into the Y values of the XY vector. Next, subtract the number of intervals multiplied by the size of the intervals from the result of the value of the low extreme in the set of X values minus the independent variable. We substitute actual values (refer to Figure 1) for the above and conclude with statement 15, where Xparameters: XY[0] = 5, XY[I] = 0, XY[2] = 15000, XY[3] = 3000

AND X = 10000

AND

case those necessary for accommodating nonlinearity in the continuous simulation of lumped-parameter systems. Various control constructs were discussed in light of syntactical rules pertaining to five languages for general-purpose system development and model specification. Topics addressed via the five modules and the accompanying discussion include contrasting modes of variable definition, conditional/unconditional transfer of program control, decision logic, range of built-in primitive functions, varying index origin and subscripting, notions of structured design, strong data typing versus single conceptual types, repetition by leading test loops versus noniterative one-pass approaches to algorithm construction, and alternative multistep algorithm styles in place of the incorporation of standard formulas.14 The algorithms as presented can be applied directly in simulation and effectively expand the avenues beyond earlier work in APL to those researchers familiar with the languages as described. Where a model builder is skilled in other languages, the modules can provide a template for translation depending on the shared characteristics of languages and characteristics of interpreters and compilers. Finally, it is hoped that the foregoing is useful when considering the stylistic factor in algorithm design. Perhaps some of the modules will find their way into the application toolbox of researchers, programmers, and model builders.

References I 2 3 4 5

Y parameters: XY[4] = 0.0, XY[5] = 0. I, 6

XY[6] = 0.75, XY[7] = 1.0, XY[8] = 3.0, XY[9] = 5.0

AND

7 8

I = 3, J = 7, K = 8andD = 1000

THEN M = I + 1000.(3-1)/3000

9 10 II

YIELDING M = 1.67 (the multiplier coefficient) 12

Conclusion

13

The stated aim of this work is to focus attention on stylistic differences in specifying algorithms, in this

14

Cundiff, W. E. Simulating system dynamics in A P E Appl. Math. Modelling 1979, 3, 312 Cundiff, W. E. On the specification of diverse models in A P E Simulation 1985, 45, 138 Cundiff, W. E. Model integration and algorithmic transparency in APL. Appl. Math. Modelling 1984, 8, 445 Karplus, W. J. The spectrum of mathematical modelling and systems simulation. Math. Comput. Simulation 1977, 19, 3-10 Roberts, E. B. Managerial Applications of System Dynamics. MIT Press, 1978, p. 72 Roberts, N., Anderson, D., Choate, J., Deal, R., Garet, M., and Shaffer, W. Instructor's Manual for Introduction to Computer Simulation. Addison-Wesley, 1982, pp. 34-54 Cundiff, W. E. The approach to system dynamics at Griffith University, Australia. DYNAMICA, 1983, 9-11, 91-94 Peterson, N. D. MIMIC, an alternative programming language for industrial dynamics. Socio-Economic Planning Sciences, 1972, 6, 319-327 Gordon, G. System Simulation. Prentice-Hall, 1978, pp. 67-79 Kreutzer, W. System Simulation: Programming Styles and Languages. Addison-Wesley, 1986, p. v Karplus, W. J. and Cheung, S. Performance evaluation tools for simulators consisting of networks of microcomputers. Procceedings of the Summer Computer Simulation Conference, Vol. I, North-Holland, 1984, p. 322 Hellerman, H. and Smith, 1. APL/360 Programming and Applications. McGraw-Hill, 1976, pp. 90-94 Delfosse, C. M. Continuous Simulation and Combined Simulation in SIMSCRIPT H.5. CACI Inc.-Federal, 1976 Makoui, A. Software interface for multiprocessor simulation. Engineering Report UCLA-C5D-860014, Los Angeles, 1986

Appl. Math. Modelling, 1988, Vol. 12, June

345