w

r \
LIBS-am. 7
Michigan State
University

K J

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.

DATE DUE DATE DUE DATE DUE

 

 

"SEP 2 8 1994

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmative Action/Equal Opportunity Institution
cmmt

 

 

A MULTIPORT APPROACH TO MODELING AND
SOLVING LARGE-SCALE DYNAMIC SYSTEMS

By

Yanying Wang

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Mechanical Engineering Department

1992

ABSTRACT

A MULTIPORT APPROACH TO MODELING AND
SOLVING LARGE-SCALE DYNAMIC SYSTEMS

by

Yanying Wang

One of the major challenges in the simulation of Large-Scale Dynamic Systems (LSDSs)
is to increase the efficiency of computation while maintaining the desired accuracy of
solution. In particular, when repeated runs are made of the same model with varying
input conditions and parameter values, then the time required for simulation becomes a
very important factor. The simulation of a LSDS typically includes three steps:
generating a computer-based model, sorting the system equations for solution, and
solving the equations numerically. There is a great practical benefit to improve the

efficiency with which any of these steps is executed.

In this work the modeling of LSDSs by means of bond graphs is extended. Two new
system graph node types, the dynamic block and the dynamic multiport, are introduced.
Each node type is defined by a set of differential-algebraic equations. An implementation
of the new node types has been made in an existing software, namely, ENPORT. The
equations sorting algorithm has been extended to include the new node types. A path-

order matrix has been generated to assist in finding an efficient solution order and to

reduce the amount of ealculation required to evaluate the Jacobian matrix.

Models that include dynamic nodes can be partitioned into several linked submodels, each
of which is itself a dynamic system. A plan has been made to assign to each dynamic
subsystem its own integrator (i.e. , integration algorithm) and its own step size. Hence
it is possible to organize the system solution by assigning multiple integrators and
multiple step sizes to the complete model. The multiple step size feature has been
implemented in software. An example that illustrates the potential for increasing the

efficiency of simulations by using multiple step sizes is presented.

To my parents, Qihao Wang and lie Shao.

iv

ACKNOWLEDGMENTS

The author wishes to express her deep appreciation to Dr. Ronald Rosenberg for his
advice, encouragement and support throughout the course of this research work as well

as guidance through the entire graduate program.

Special thanks are also due to the other members of the guidance committee Davor
Hrovat, Ford Motor Company, Hassan Khalil, Department of Electrical Engineering,
Michigan State University and Philip Fitzsimons, Mechanical Engineering Department,
Michigan State University for suggestions and taking an active part in the process of
completion of this dissertation. Their invaluable advice and suggestions are fully

appreciated.

Grateful acknowledgement is extended to Mechanical Engineering Department of
Michigan State University and Rosencode Associates Inc. for the financial support they

provided.

Finally, the author is most grateful to her husband, Ye Tian, who assisted in editing, her
daughter, Iris, and her son Eric. Their generous love, unlimited patience and endless
support brought her to the completion of the final chapter.

V

TABLE OF CONTENTS

LIST OF FIGURES ...................................... ix
LIST OF TABLES ...................................... xi
LIST OF ABBREVIATIONS ................................ xii
1. INTRODUCTION ..................................... 1
1.1 The Problem Statement ........................... 1

1.2 Literature Review .............................. 2

1.2.1 Modeling of LSDSs ........................ 2

1.2.2 Computational Methods for Solving LSDSs .......... 5

1.3 Dissertation Organization ......................... 6

2. DYNAMIC NODES IN SYSTEM GRAPHS ..................... 8
2.1 Generalizing Dynamic Node Types ................... 8

2.2 Defining DBN and DMN as Modeling Tools ............. 10

2.2.1 Dynamic Block Nodes (DBNs) .................. 10

2.2.2 Dynamic Multiport Nodes (DMNs) ............... 12

2.2.3 Equations of DNs .......................... 13

2.2.4 Causal Considerations ....................... 16

2. 3 Software Implementation .......................... 16

2.4 An example .................................. 20

2.4.1 Model ................................. 20

2.4.2 Simulation Comparison ...................... 24

3. ORGANIZING SYSTEM EQUATIONS FOR SOLUTION ............ 26
3.1 The Computational Graph ......................... 26

vi

3.2 Depth-First Search for Sorting ...................... 30
3.2.1 Basic Algorithm ........................... 30
3.2.2 Algebraic Loops .......................... 32
3.2.3 Dynamic Nodes ........................... 37
3. 3 Modified Breadth-First Search for Path Identification ........ 39
3.3.1 Basic Algorithm ........................... 40
3.3.2 Solving Order ............................ 40
3.3.3 Path-Order Matrix ......................... 42
3.4 Generation of Jacobian Status Matrix .................. 48
4. IMPROVEMENT OF COMPUTATIONAL EFFICIENCY FOR LSDSs
DESCRIBED BY BOND GRAPHS ........................ 55
4.1 Solution of Sorted System Equations ................... 56
4.1.1 Submodels .............................. 56
4.1.2 Structure Decomposition ..................... 58
4.2 Multi-rate Solutions ............................. 63
4.2.1 Direct Solution Methods ...................... 63
4.2.2 Multiple Integrators ........................ 66
4.3 Software Implementation .......................... 72
4.4 An Example ................................. 73
5 . EFFICIENT COMPUTATION OF ALGEBRAIC LOOPS ............ 78
5 .1 Problem Description ............................ 78
5.2 Iterative Method ............................... 79
5.3 An Algorithm for Finding an Efficient Set of
Iteration Variable .............................. 81
6. CONCLUSION ....................................... 87
6.1 Summary and Discussion of Results ................... 87
6.2 Suggestions for Future Research ..................... 9O

vii

APPENDIX A. BASIC SYSTEM GRAPHS ....................... 92

A.1 Block Diagrams ............................... 92
A2 Bond Graphs ................................. 96
APPENDIX B. AN EXAMPLE OF CAUSALITY ASSIGNMENT ........ 99
APPENDIX C. LISTING OF SORTING SUBROUTINES ............. 102
APPENDIX D. LISTING OF MSS CODE ....................... 118

APPENDIX E. ALGORITHM FOR OBTAINING THE PATH-ORDER MATRIX 147

THE BIBLIOGRAPHY ................................... 151

viii

LIST OF FIGURES

Eism nae:
2.1 Symbol of a DN ..................................... 10
2.2 Symbol of a m-input, n—output DBN ......................... 11
2.3 Symbol of a n—port DMN ............................... 12
2.4 Hierarchical structure of system graphs ...................... 14
2.5 List of function types in modified ENPORT .................... 20
2.6 An example of system graph with DNs ....................... 21
2.7 System graph with standard nodes .......................... 23
2. 8a Model with standard node types ........................... 24
2. 8b Model with DNs .................................... 25
3.1 Bond graph model of a power transmission system ................ 27
3.2 Computational graph and its edge-node table .................... 30
3.3 An example of SCC ................................... 32
3.4 Two-damper spring system .............................. 33
3.5 Bond graph of two-damper spring system ...................... 34
3.6 Computational graph of two—damper spring system ................ 36
3.7 An example of standard computational graph ................... 37
3.8 Characteristics of dynamic nodes ........................... 38
3.9 Representation of DNs in CGs ............................ 39
3.10 SCG for power transmission system ......................... 42
3.11 An example of path identification .......................... 45
3.12 Mechanical system and its bond graph model ................... 50
3.13 CG of a mechanical system .............................. 52
4.1 Submodel structure of system model ......................... 57
4.2 Structure model ..................................... 60
4.3 Flow diagram of surfing system equations ..................... 64
4.4 Diagram of multiple integration approach ...................... 67
4.5a Program flow chart, part 1 .............................. 69
4.5b Program flow chart, part 2 .............................. 70
4.5c Program flow chart, part 3 .............................. 71
4.6 Submodel structure of system graph ......................... 74
4.7 Time response ...................................... 77

ix

5.1 An example of SCC ................................... 82

5 .2 An example for finding the near-minimum independent set ........... 85
A.l Diagram of actual physical elements ........................ 93
A.2 Block diagrams of input-output relation ....................... 94
A.3 Basic multiport fields .................................. 98
B.1 DMN used to explain causality assignment ..................... 99
El An acych digraph .................................. 148
E.2 Construction of solving order for an acyclic digraph .............. 149

LIST OF TABLES

Table was
2.1 Classification of equation structure for dynamic nodes .............. 15
2.4 Causality and input-output relationship ....................... 25
3.1 Equation data listing of bond graph ......................... 28
3.2 SCCs in a two-damper spring system ........................ 36
3.3 List of SCC ........................................ 41
3.4 List of system equations ................................ 51
3.5 Output—input paths .................................... 53
4.1 Integration methods for

the ordinary-differential equations .......................... 65
4.2 CPU times for solution ................................. 76
5.1 List of independent sets in SCC ........................... 86
A.1 Basic building blocks used for modeling systems ................. 95
A.2 Bond graph atom nodes ................................ 97
B.1 Causality and input-output relationship ...................... 100

BFS
CG
DFS
DBN
DMN
DN
GMM
LSDS

MSS
SCAP
SCC
SCG

LIST OF ABBREVIATIONS

Breadth-First Search

Computational Graph

Depth-First Search

Dynamic Block Node

Dynamic Multiport Node

Dynamic Node

Graph-defined Macro Multiports
Large-scale dynamic system

multiple integration algorithm method
multiple step sizes method

Sequential Causality Assignment Procedure
Strongly Connected Component
Standard Computational Graph

xii

l . INTRODUCTION

1.1 The Problem Statement

There is widespread interest in the analysis and simulation of large-scale dynamic
systems (LSDSs). For this reason, a large number of software programs have been
developed to perform simulation of such systems. Of particular interest are the so-called
"physical systems” which are discussed in this dissertation. Such systems have an energy
basis, e.g. , electrical networks, rigid-body mechanical and fluid power systems. Since
such physical systems often include multiple energy domains, the bond graph technique
has proven to be a great help in modeling them. The modeling process with bond graphs
is straightforward and, furthermore, the bond graph also implies a complete set of system

equations with physically meaningful state variables.

When methods of LSDSs are used in iterative design processes, such as parameter
optimization and/or controller design, the time required for the simulation becomes an
important factor. Improvements in simulation of LSDSs that increase efficiency while
maintaining desired accuracy of solution are valuable to engineering practice. Such

improvements are the subject of this work.

It is useful to divide the simulation of LSDSs into three main steps. These are:

l

2
(1) Construct a computer-based model of the system to be simulated. From the

model derive a set of system equations.

(2) Sort the system equations into a solution order. Details of equation ordering
typically affect the solution efficiency.

(3) Perform a simulation of the response of the system by solving the equations

numerically.

Improvement in executing any of these three steps benefits engineers and scientists in a
wide variety of industries. Even small improvements have a large multiplier in terms of

engineering productivity.

1.2 Literature Review

1.2.1 Modeling of LSDSs

There are important challenges in the analysis and simulation of LSDSs. Since the
amount of computational effort required to analyze a LSDS usually grows at a rate
greater than the size of a system, simulating LSDSs may become very time-consuming
and possibly impractical. Many books have been written, research papers published, and
computing algorithms developed concerning the analysis and simulation of LSDSs
(Laddle, 1975; Siljak, 1978; DeCarlo and Sacks, 1981; Kobayashi et al., 1978;
Constantinescu, 1982; Wang and Jamshidi, 1982; William, 1983; Martinez-Benet and
Puigianer, 1988). Problem solving by modeling and simulation is an iterative procedure

(Broenink, 1990). It can be described as follows:

(1). formulate a quantitative model;
(2). carry out a numerical simulation;
(3). check the results to see if they satisfy the desired requirements;

(4). if not, modify the model and repeat steps 1 through 3.

In order to deal with problems involving large size and complexity in a systematic and
efficient way, the manner of description of objects and processes is important. A model
to enhance our understanding of a problem can take several forms. Block diagrams often
are used to display mathematical models in a form that allows us to understand the
interactions occurring among system’s elements (William, 1983). A mathematical (i.e.,
equation) model, which is a description of mathematical relations among system
variables, is often found to be useful and is widely used, due to its generality.

Typically, there is a variety of mathematical descriptions that can be applied to a given
system, and engineers must be prepared to decide what form and level of complexity are
most consistent with the objectives of the study and the available solution resources. The
nature of system involved has a strong influence on the selection of modeling and
simulation methods. The process of using a mathematical model to determine certain
features of the cause-and-effect relationships of a system is referred to as ’solving the

model’ (Close and Frederick, 1978).

The bond graph is a modeling tool introduced by Paynter (1960). In bond graph

modeling information concerning the interconnection among components of a system is

4
given by the (power) bonds. The node symbols in the graph denote the action of

physical components and effects. A standard bond graph is composed of elements from
the basic set {C, I, R, Se, Sf, TF, GY, 0, l}, which includes multiport capacitance,
inertia, dissipation, sources of effort and flow, modulated transformers and gyrators, and
the ideal junction elements 0 and 1, respectively (Rosenberg, 1971). Several texts
describe the fundamentals of bond graph modeling (Rosenberg and Karnopp, 1983;

Karnopp and Rosenberg, 1975).

Many engineers and scientists have found it useful to apply the bond graph technique to
solve a variety of problems. For example, the application of bond graphs to mechanisms
shows two advantages. First, the analyst may consider a mechanism-dynamics problem
from a point of view that often provide new insights into the complex behavior of such
a device. Second, the kinematics and dynamics of a mechanism are represented in a
form that is currently being used in a wide range of both engineering and non-
engineering disciplines. A bond graph study can be made of the kinematics and
dynamics of a general mechanism treated as a component of a dynamic system. Once
the kinematic mechanism multiport is developed for a particular mechanism, a general
dynamic model can be constructed parametrically in the I—, R-, and C- matrices (Allen
and Dubowsky, 1977). Then the multiport can be reused. Bond graph modeling has
been used to simulate electro-hydraulic systems (Dransfield, 1979), large mechanieal
systems (Bus and Tiemego, 1985), automotive power trains (Hrovat et al. , 1985), and

electromagnetic actuators (Karnopp, 1985). A recent bibliography cites many

5
applications of bond graph modeling to engineering system problems (Filippo et al. ,

1991).

1.2.2 Computational Methods for Solving LSDSs

Due to the complexity of LSDSs, considerable effort has been devoted to the
development of numerical techniques to solve them. Several approaches have been used
to simulate the transient response of LSDSs, while retaining a reasonable computational
accuracy. Three important approaches are the component connection method, the
perturbation method, and the diakoptics method. These are methods for model—order

reduction. Another important approach is to use multiple-time scales in the integration.

The component connection model of a dynamic interconnected system (DeCarlo and
Sacks, 1981) is a set of equations describing the dynamics of independent components
and a set of algebraic equations describing the interconnection properties. There are two
principal integration algorithms used to simulate this model, namely, the Sparse Tableau
algorithm and the Relaxation algorithm. Unfortunately, many important engineering

problems cannot be represented by this model.

Perturbation methods are useful for dealing with a system that can be approximated
effectively by a system of simpler structure. Perturbations are divided into two classes:
regular and singular perturbations. In regular perturbation (Kokotovic et al. , 1969) the

system is connected by some weak connections. It can be decomposed into two (or

I“

 

6
more) completely independent sub-systems by ignoring the weak connections. In singular

perturbation there is a perturbation to the left-hand side (i.e. , derivative term) of a
differential equation. Ignoring the perturbation term leads to a reduced ”slow” sub-
system. The ”fast" sub-system can be obtained by stretching the time scale. Therefore,
one has to solve the "slow” and the "fast” (or boundary-layer) models (Kokotovic et al. ,
1986). The main difficulty is finding a form of the system equations from which the

standard form, with the "slow” and the ”fast" variables separated, can be derived.

The original idea of diakoptics was suggested to solve a problem in two steps: (1) at sub-
system levels: tear a system apart into logical groups and solve the sub-systems
independently; and (2) at interconnection levels: combine these results with the
connection matrix to obtain an overall solution (Wu, 1976). This approach is applied

frequently in circuit analysis, but is not so useful for other types of models.

1.3 Dissertation Organization
This dissertation is devoted to a multiport approach to modeling and solving large-scale

dynamic systems with emphasis on improving computational efir‘ciency in solution.

In Chapter 2 the modeling tool to be used is presented. Two new types of dynamic
nodes are defined, a block diagram element and a bond graph element. The forms of the
equation set that can be used to support the dynamic nodes are described. For bond

graph elements, the causal constraint at the system interconnection level are considered

7
and discussed. The standard method for formulating mathematical models from system

graphs is extended to include system graphs with dynamic nodes. Software

implementation is described and an example is given.

Chapter 3 presents the fundamentals of two algorithms in graph theory and their
application to organizing system equations for solution. The models can include
algebraic loops and dynamic nodes (DNs), which are identified from the system graph.
These effects need specific computational treatment. A path-order matrix is created to
establish the solution order. The use of a path-order matrix to find the Jacobian Status

Matrix is also described.

In Chapter 4 the definitions of structural submodels and structure decomposition are
presented. A graph-oriented decomposition method is applied. Based on the
decomposed model, a multiple step size method and a multiple integration algorithm
method are suggested to solve for system time response. A flow chart depicting the

algorithms is given in this chapter and an example is shown.
In Chapter 5 the digraph analysis method is extended to finding an efficient set of
iteration variables for algebraic loops. The algorithm is developed and its application is

explained by examples.

Chapter 6 concludes this dissertation with suggestions for future studies.

2. DYNAMC NODES IN SYSTENI GRAPHS

2.1 Generalizing Dynamic Node Typos

A multienergetic model described by eight types of basic block elements (listed in
Appendix A, Table A. 1) and nine types of bond graph atom nodes (listed in Table A.2)
is referred to as a standard system graph. Such a combination of elements is useful in
modeling systems that contains control elements such as transfer functions. It is also

useful in showing nonlinear modulating effects.

To obtain a standard system graph, certain interconnection rules between blocks and
multiports have to be followed:
(1) For {C, I, R, TF, GY, Se, Sf} nodes, only input signals are allowed.

(2) For {0, 1} nodes, only output signals are allowed.

A standard system graph can be transformed into a mathematical (i.e. , equation) model.
Equations of a mathematical model can be derived from the defining relations of the
blocks and nodes. Block diagrams are direct pictorial representations of equations, while
bond graphs represent the model equations but have no input/ output details on them. To
generate an input/ output set of system equations one can use the Sequential Causality
Assignment Procedure (SCAP, Rosenberg and Karnopp, 1983).

8

9
In standard bond graph and block diagram representations, the types of nodes are pre-

defined. The equation(s) corresponding to each node also has (have) specific form.
General approaches to simulation of physical systems are often effective, but when it
comes to model large-scale systems, difficulties arise. The process of aggregating parts
of a model into a single unit is one way to reduce complexity. To accomplish this, a

new atom type node and a new block element are developed in this section.

The Dynamic Node (DN) is a new modeling tool which is defined in the framework of
both block diagram and bond graph language. Similar to the graphic definition of
general nodes, a DN may have ports connected to it and may have signals in and signals
out. It is worth mentioning that the -C and -I nodes are the special examples of DNs.
A graphical representation of DN is shown in Figure 2.1. The general DN has a set of
bonds (B) with an associated set of inputs (0,) and outputs (YB). There is a set of
signals in (8,) with variables Us and a set of signals out (So) with variables Y3. A state

vector X is associated with DN.

10

 

(X)
(YB) : N: 30

U
< .) a.)

Si (Us)

 

 

 

 

 

Figure 2.1 Symbol of a DN

2.2 Defining DBN and DMN as Modeling Tools
Depending upon coupling specifications with other elements, two types of DN can be

defined: the Dynamic Block Node (DBN) and the Dynamic Multiport Node (DMN).

2.2.1 Dynamic Block Nodes (DBNs)
In some cases all the connections between a dynamic submodel and the rest of the system

carry no power information. Thus they can be represented by signals. The term DBN

1 1
is used for this type of submodel.

Definition
* A DBN is a block node with m input signals and n output signals. The input
and output variables, together with a state vector, are related by a set of ordinary

difl’erential equations and a set of algebraic output equations.

The graphical symbol for a DBN is a block, defined as type DBN, with one or more
signal inputs and one or more signal outputs. An example is shown in Figure 2.2.

The state vector is implicit. The block label is arbitrary.

 

(X)
”1 : DBN : V‘
U

 

 

 

 

 

 

 

 

 

 

 

Figure 2.2 Symbol of a m-input, n-output DBN

12
2.2.2 Dynamic Multiport Nodes (DMNs)

The node type DMN used here should be distinguished from the macro type of multiport,
which is defined by a set of bond graph atom nodes {C, I, R, SE, SF, TF, GY, 0, l}
and their connections. The properties of Graph-defined Macro Multiports (GMM) are

derived from the details of their node sets.

Definition

"' A DMN is a multiport node with n ports and m signals in. The efi'ort and flow
variables at the ports, together with a state vector X, are related by a set of

ordinary difl’erential equations and output equations.

 

9 GM

DMN 4——————

:3/ '\

Figure 2.3 Symbol of a n-port DMN

 

 

 

 

13
The graphical symbol for a DMN is a multiport, defined as type DMN, with as many .

ports and input signals as nwded. Figure 2.3 shows an example of a n-port DMN with

m signals.

2.2.3 Equations of DNs
Based on the discussion in 2.2.1 and 2.2.2, the concept of system graph with newly
developed DBNs and DMNs is established. The hierarchical structure of the system

graph is summarized in Figure 2.4.

The equation forms supporting DNs are independent of the type of elements, DBN or
DMN, provided the DMNs are causally oriented. Therefore, the equations can be
described in general format. Four classes of equations characterizing DN models are

listed in Table 2.1.

14

 

 

SystemGraph‘

 

Dynamic
5...... { ,
Atoms m TRFN
”W i
_ Bloch Encoded { DBN
Static
Means
Dynamic
ma {
Static
Am ‘[ Dynamic - DMN
" Mm Extended
M Static

 

Figure 2.4 Hierarchical structure of system graphs

 

15

Table 2.1 Classification of equation structure
for dynamic nodes

          

 

D.E. explicit D.E. implicit A

  
      
 

   

 

   
 

Output 2, = «X... r, U.) o = o(x_, r, U“, 22,)
explicit ‘1'(X,., t. U.) Y... = ‘l'(X,,,. t. U.)

Output 2,, = o(x_, t, U.) o = cog, t, U_, X.)
implicit 0 = ‘1'(X.., t. U... Y.) 0 = 70!... t. U... Y.)

 

 

 

In Table 2.1 the following definitions are used:
XIII is state vector of the node,
U, is input vector to the node,
Ym is output vector from the node,

and t is time (independent variable).

From the functional point of view, DN is different from other dynamic nodes (e.g. , C
and I elements) in that its equations are not predefined. The order of differential
equations and the class of equation structure may vary. In later chapters, we discuss how
the newly defined DNs can bring efficiency and flexibility to modeling and solving the

whole system.

16
2.2.4 Causal Considerations

For DBNs the input and output variables are defined explicitly by the graph, since they
are signal related. On the other hand, the inputs and outputs related to the ports of a

DMN can not be determined prior to causality assignment.

A mathematical model can be derived from standard system graphs by means of a
standard formulation method (Karnopp and Rosenberg, 1975). In this process the
causality of all bonds must be specified. During causality assignment in bond graph
modeling, according to the SCAP, causath is assigned to all sources and storage
elements, and extended as far as possible into junction structures. Note that among atom
nodes, some (such as SE, SF, 0, 1) are constrained with respect to possible causalities,
some (such as C and I) have preferred causalities, and some (such as R) are indifferent
to causal orientation. The causal orientation of a DMN must be consistent with its

equation definitions. An example is presented in Appendix B.

2.3 Software Implementation
The existing ENPORT software for system graph modeling and simulation has been
extended to include the new node types DBN and DMN. Modeling details are described

at the graph level and at the equation level.

(1) Graph implementation

DBN and DMN nodes are created in the system graph in either line code or

17
graphic form. Required data are: node name, node type, connector names,

connector types, and connector orientations. One more piece of information
needs to be known for DMNs, i.e. , the eausality requirements at the ports. The

causality is defined by asking the user to answer the following question for each

port:

Enter ’E' if input is effort, ’F’ if input is flow, ’N' if input is indifferent.
Bond -1 (N): F <ret>

Bond - 2 (N): E <ret>

Subsequently, the program will analyze the causality of a whole system to check
whether there exists a causality conflict. If there is, the program will give a clear

error message and stop.

(2) Equation specification
Each DMN is defined by a FORTRAN subroutine. An example is shown below.
The executable statements associated with DOF, NUML, and DESC are text that
can be displayed during execution to remind the user of the meaning of the
particular DN definition. The executable statements in the rest of the subroutine
define the state vector size, and they evaluate the state derivative vector and the
output vector. Note that the dimensions of the input and output vectors are

known from the graph.

18

CDN>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
SUBROUTINE DNlTlME,X,P,NSX,SX,DSX,Y,DOF,NUML,DESC)

C

C---- PROGRAMMING: Your name. The date.

C

C---- DESCRIPTION: Short description of subroutine.
C

C---- INPUTS: TIME, current time

C X, input values

P, parameter values

SX, state variables

DOF, = .TRUE. if description requested,
= .FALSE. if evaluation requested.

C

C

C

C

C

C---- OUTPUTS: Y, output values

C NSX, no. of state variables
C DSX, derivative of state variables

C NUML, number of lines in description
C DESC, description of function

C

C

---- DECLARATIONS:
CHARACTER DESCIZO) ' 72
DOUBLE PRECISION TIME, X(20), SXIZO), DSX(20), P(20), Y(20)

INTEGER NUML, NSX
LOGICAL DOF
EQGGDN‘I‘I‘I‘I‘QG‘I"!§§**§***§I**!§I*§*§*§*ifiifiiirl-‘I‘I‘liil-iiiiiii
C
IF (DOF) THEN
C ------ Description section (max 20 lines)
NUML= 5
DESC(1)=' DN: set as 1st order differential equation.’
DESClZ) = ' DSXI1) = -(5.0/0.1)*SX(1)-(1.0/0.25)*SX(2)’
DESCl3) = ’ DSXlZl =(1.0/O.1)”SX(1l-(O.3/O.25l*SX(2l-X(2)'
DESCl4) = ’ Y(1) = SXl1)/O.1’
DESCISl = ’ Y(2) = SX(2)/O.25’
C
ELSE
C
C ------ Set the number of state variables
NSX= 2
C

C ------ Evaluation section

19

DSX(1)= -(5.0/O.1)*SX(1 H1 .O/O.25)*SX(2)
DSXlZ) = (1 .O/O.1 )‘SXll l-(O.3/O.25) *SX(2)-X(2)

C
Y(1) = SX(1)/0.1
Y(2l = SX(2)/O.25
C
ENDIF
C
RETURN '
END

C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

These subroutines should be compiled and linked with the main program for the

complete definition of DNs in software.

Modifications have been made to ENPORT in graphical modeling and equation
definition. In the graphical modeling DBNs and DNle are defined with the same
manner as other nodes. Equations to support DNs are collected in a file of such
subroutines. Each distinct routine is named ZZDN,j (i = 0, 1, ..., 9;j = l, 2, ..., 9).
Up to 99 sets of ordinary-differential-output equations are allowed in this file. The use
of these subroutines is similar to the use of User_Defined_Subroutines, ZZSUfi, in
ENPORT (see Rosenberg, 1990, and Appendix B for details). As a result, the function

types in ENPORT are enhanced, as shown in Figure 2.5.

2O

 

ind-Lu] W W
ZZSU W

W

ll fl 0102...99 0102...”

 

 

 

Figure 2.5 List of function types in modified ENPORT

2.4 An Example

2.4.1 Model

A system graph containing DBNs and DMNs has been built in Modified ENPORT and
is shown in Figure 2.6. This model includes one DBN and one DMN. Equations for
DBNl model a feedback controller and equations for DMNl model a motor. These

equations are listed as Equations (2.1) and (2.2) respectively.

21

 

 

 

 

 

 

Figure 2.6 An example of system graph with DNs

The equations for the DBN are:

“r .
-2.0 x + u2

x'
y

where x denotes THETA, ul denotes W1, uz denotes SI, and y denotes S2.

(2:. 1)

22
The equations for the DMN are:

121 =ul -50x1-4.01t2

x,=rox, -1.21t2-u2
yl=10x1 (2.2)
3‘45

where x, denotes P.E2, x2 denotes P.M2, ul denotes E.1, u, denotes E.2, yl denotes F. 1,

and y, denotes F.2.

The load, composed of an inertia (I) and a friction effect (R), is driven by the motor
through a stiff shaft (C). The input voltage to the motor is generated by node SE. The
feedback controller takes in the desired position 81 and the load velocity W1 and outputs

an actuator signal 82.

An alternative system graph with more details is given in Figure 2.7. Each DN has been
expanded into a set of standard nodes. The DBNl is defined by the standard block atoms
{SUM, GAIN, INT}. The DMNl is defined by the standard multiport set {R1, 1A, 11,
GY, R2, 1B, 12}. The physical parameters of the standard nodes were used to derived
the DBNl and DMNl equation details. Parameters chosen for the physical components
are as follows:

bl = 5.0 (fl) resistance

L = 0.1 (Henry) inductance

b2 = 0.3 (N.s/m)
m2 = 0.25 (Kg)

s4 = l0.0*s3

s2 = sl-s4

sl = 1.0 (V)

k = 100.0 (N.m/rad)
m = 1.0 (Kg.m2)

b3 = 0.5 (N.m.s/rad)

E1 = l.0*sl

23
frictiOn coefficient
motor inertia
feedback gain
negative feedback
reference position
shaft stiffness
load inertia
rotational friction coefficient

actuator voltage gain

 

 

SRC

 

 

 

 

Figure 2.7 System graph with standard nodes

 

24
2.4.2 Simulation Comparison

A simulation run is made from initial time (0) to final time (8 seconds). Figure 2.8
shows the behavior of the load velocity W1 and the load position THETA. The load
velocity is adjusted by controller to have the load approach a constant position. Results

obtained using the DBNl and DMNl match those obtained by the model with standard

node types.

SCALING

wr ‘,-—-' ' “"' -----
1.312-01 “\I\\\N ’.r"
-s.072-03 ,.'"
mam . "
5.233-01 ,"
0.002+00 'x'

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

/ f.” \\
0.00 0.20 0.40 0.60 0.80 1.00
' TIME * l OE l

LEGEND: W1 —— THETA - -

Figure 2.8a Model with standard node types

25

 

 

SCALING

at . ..... - - - . . -
Lair—or \
-s.In-03
mm
5.23l-OI .'
0.00:900 ::\\\

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.00 0.20 0.40 0.60 0.80 1.00
TIME '10P. l

LEGEND: Hl— THE'I'A - -

 

Figure 2.8b Model with DNs

 

3. ORGANIZING SYSTEM EQUATIONS FOR SOLUTION

3.1 The Computational Graph

A directed graph can be constructed to represent the structure of the system equations
derived from a bond graph model. This computational graph (CG) has system input
variables and state variables as starting nodes and derivatives of state variables and
system outputs as ending nodes. The construction of the graph can be done as described

when the bond graph model contains integral causalities.

Suppose we have a power transmission system as depicted in Figure 3.1. This is a
model of an inertial load driven by a motor through a slipping clutch. Bond M is the
shaft connection from the motor to the clutch on the upstream side and bond L is the
shaft connection of the clutch to the load on the downstream side. The motor is modeled
by a torque source (SEM) and an internal equivalent resistance (RM). The clutch is
modeled by a viscous friction coupling (RC). The load is modeled by a rotational inertia

(IL) and friction effect (RL).

For a given system graph model, a complete set of system equations can be developed.
The list of such equations, specified as relations between input-output variables, is given

in Table 3.1 for the example of Figure 3.1. In the table, "E” refers to effort, "F" refers

26

27

to flow, and "P” refers to momentum. The suffixes represent the names of bonds.

 

SEM—NlMl—AOC—AIL

w 1%!”

IL RL

 

 

 

Figure 3.1 Bond graph model of a power transmission system

28

Table 3.1 Equation data listing of bond graph

 

 

 

 

# of equation Output vbl name Function type Input vbl name I
1 E.Ml CON ---
2 F.Ml ASGN F.M2
3 E.M2 SUM E.Ml

EC
4 RM ASGN F.M2
5 F.M2 GAIN E.M2
6 EM ASGN EC
7 RC SUM F.M2

F.Ll
8 EL ASGN EC
9 EC GA2N EC
10 FL ASGN F.Ll
11 E.Ll SUM E.L2

EC
12 F.L2 ASGN F.Ll
l3 F.Ll A'l'l‘ RM
14 PM INTEG E.Ll
15 E.L2 GAIN F.Ll

 

 

 

 

 

29
Figure 3.2 depicts the input-output relationships in a CG. Each directed edge represents

an entry in the input list. Each node represents an equation and its associated output.
The node numbers in Table 3.1 correspond to equation numbers in Figure 3.1. The

variable labels are shown only for ease of interpretation.

The problem of sorting the system equations into a suitable order for sequential solution
is transformed into the problem of finding a precedence order for associated CG. There

are two distinct situations to consider.

(1) The CG has no cycle (i.e. , directed loops). This case is relatively simple and can
be treated by any one of several search algorithms.

(2) The CG contains one or more cycles. Each such cycle must be identified (e.g. ,
as a strongly connected component) and a reduced CG generated. Then the case

is that of (1) above.

30

 

rm Fm EN EL

E.MI

an an,“

..QOU’UNH
p

.0..U‘.HU

..ddauouuu

u
p
p

.Ud

p
O

 

 

 

 

Figure 3.2 Computational graph and its edge-node table

3.2 Depth-First Search for Sorting

3.2.1 Basic Algorithm

The Depth-First Search (DFS) technique is a method of scanning a finite graph (Even,

1979). For the development of sorting algorithms, the basic idea and definitions of DFS

are briefly summarized in this section.

31
Definition 1:

A finite directed graph G(V,E) consists of a finite set of vertices V= {v,,v,, - - ~vJ
and a finite set of edges E = {e,,e,, - - -e,,},' each edge e is incident to the elements

of the ordered pair of vertices (u, v), fiom u to v.

(A directed graph is often referred to as a digraph). With DFS, let us select and visit
a vertex a, and then visit a vertex b adjacent to a, and then continue with a vertex c
adjacent to b (but different from a), and then an ”unvisited" (1 adjacent to c, and so forth.
As we go deeper into the graph, we will eventually visit a vertex y which has no
unvisited neighbors. When this happens, we return to the vertex immediately preceding
y in the search and continue the procedure. If the particular search terminates, then a
new starting vertex is sought and the procedure starts again. If vertices are labeled

sequentially as they are visited, then the labels can be used to derive a searching order.

One of the very useful application of DFS is for identifying strongly-connected

components in a digraph.
Definition 2:
Let G(V,E) be a finite digraph. G, is a strongly-connected component of G if

given (u,v e V (6)), there exists a directed path from u to v and from v to it.

An example is given in Figure 3.3. The vertices a, b, c belong to a SCC of the digraph.

32
The vertex d belongs to a SCC of the digraph.

 

an \/ ed

 

 

 

 

Figure 3.3 An example of SCC

3.2.2 Algebraic Loops
An algebraic loop arises in a system model when the value of a variable depends on
itself. A simple example of a system which contains an algebraic loop is given in Figure

3.4. The bond graph with causality is shown in Figure 3.5.

33

 

 

///////////

R1 K R2

 

F(t) x

 

Figure 3.4 Two-damper spring system

 

34

 

 

 

 

Figure 3.5 Bond graph of two-damper spring system

For this example the equations of the system are

1 =1; (3.1)
fl =fa
f2 =f3 (3.2)
f. =f3

‘3 = ‘1 " ‘2 ‘ ‘4 (3-3)

35

e, = F(t) (3.4)
c2 = ¢.,(f2) (3.5)
e. = 4.0:) (3.6)

f, = 0;: ca.) (3.7)

where "f" is velocity, "e" is force, and ”x” is displacement.

From equations (3.1) - (3.7), we can derive

f, = 0;: < Fe) - 0.1a.) - «w» (3.8)

Equation (3.8) is the mathematical representation of an algebraic loop. In the general

nonlinear case one has to use an iterative method to find f3, given x and to.

Figure 3.6 depicts the CG of the two damper spring system that represents the system
equations (3.1) - (3.7). It can be seen easily that path ®@®® is a directed loop.
Therefore, the node set S, = {5,6,3} is a strongly-connected component. By applying
the algorithm stated in section 3.2.1 we can discover all of the SCCs in Figure 3.6.

They are listed in Table 3.2. It should be noted that the simplest SCC is a single node

36
of the CG.

 

 

 

 

 

Figure 3.6 Computational graph of two-damper spring system

Table 3.2 SCCs in a two-damper spring system

SCC Name I Node Name

 

 

 

37
By taking each SCC as a node the CG can be reconstructed. Such a CG, which does not

include any algebraic loops, is referred to as a standard computational graph (SCG).
Figure 3.7 is an example of a SCG derived from the CG of Figure 3.6. The node labels
in the SCG are different from those of the CG because they are the labels of the SCCs,

as given in Table 3.2.

 

 

 

 

 

Figure 3.7 An example of standard computational graph

3.2.3 Dynamic Nodes
The surfing algorithm described previously must incorporate the system equations derived

from models that include DNs. Each DN contributes a node to the CG, since it relates

38
a set of inputs (3:) to a set of outputs (y). In addition, it must be classified as a node that

generates derivatives of the state variables. Figure 3.8 shows the algebraic structure of

the DNs.

 

 

0

 

IA I?

 

ll
~i
i

 

 

 

 

Figure 3.8 Characteristics of dynamic nodes

In general a DN has multiple inputs and multiple outputs. For a particular DN with m
inputs and n outputs a single node will be created in the CG, as shown in Figure 3.9.
The sorting algorithm recognizes nodes with multiple outputs and processes them

accordingly.

39

 

 

 

 

 

Figure 3.9 Representation of DNs in CGs

3.3 Modified Breadth-First Search for Path Identification

For a typical LSDS there is a great amount of calculation carried out in the simulation.
It is desirable to organize the equations in a way that reduces the computational work to
the minimum necessary to achieve the desired results. To increase solution efficiency
a modification to the existing sorting algorithm was made. The main idea of this

algorithm comes from a Breadth-First Search (BFS) approach.

40
3.3.1 Basic Algorithm

Consider the ease of a finite directed graph G, in which two vertices s and t are
specified. The goal is to find a path from s to t, if there is any, which uses the least
number of edges. The algorithm is as follows (Even, 1979).

1 Label vertex s with 0.

2 i ‘— 0.

3 Search all unlabeled vertices adjacent to at least one vertex labeled i. If

none are found, stop.
4 Label all the vertices found in (3) with H 1.
5 If vertex t is labeled, stop;

If not, i¢- i+1; go to (3).

The index i is referred to as the BFS number. If a vertex is labeled with )t(v) = k, then
there is a path of length k from s to v. Now, the BFS number has another meaning, the
length from s to v. On the path from s to t, if )t(v) < Mw), we know that in order to
reach w, vertex v has to be visited first. If Mq) = Mp) < MW), vertices q and p have
to be visited before reaching w. Thus, the BFS number also gives the order of vertices

to be visited on the path.

3.3.2 Solving Order
As an example, we inspect the power transmission system shown in Figure 3.1 again.

The CG is given in Figure 3.2. There exists one algebraic loop. After applying DFS

41
to the SCG, the algebraic loop is identified as a SCC and all of the simple nodes are also

identified as SCCs. Thus the SCG shown in Figure 3.10 is constructed based on
redefined SCCs. The SCG nodes are related to the original equation nodes as listed in
Table 3.3. The result in Table 3.3 was obtained by running the sorting program on the

example (see Appendix C for a listing of the sorting program).

Table 3.3 List of SCC

 

 

 

Name of SCC SCC pointer Equation nodes
1 l
2 2
3 3
4 4
5 5
6 6
7 10
8 1 1
9 12
10 13
1 l 14
12 15

 

 

 

 

Definition
Source vertices - vertices which have only outward edges

Sink vertices - vertices which have only inward edges.

The source vertices in Figure 3.10 are 7 and 12. The sink vertices are 1, 2, 3, 4, 5, 8,

and 9.

42

 

 

 

 

 

Figure 3.10 SCG for power transmission system

Each vertex in the SCG represents a variable or a group of variables because it
represents an equation (or possibly a set of equations). To find the solving order for a
sink vertex is to identify the directed paths to that sink vertex from all source vertices,

including the ordering of vertices to reach the sink vertex along these paths.

3.3.3 Path-Order Matrix
A path-order matrix can be developed from a SCG that is acyclic, i.e., an SCG which

by definition has no direct cycles. Recall that the basic concept of a path is a sequence

43

of vertices and directed edges from a source vertex 3 to a sink (or target) vertex t. The
source vertex is defined to be at layer 1; subsequent vertices are assigned to layers
according to when they are reached. The path-order matrix is a convenient form to use

for deriving solution order information by computer.

So far, we have defined the solving order of an output variable. The directed path from
the output t to all related inputs from a directed tree with root at t. On this tree, vertices
having same indices are defined to be at the same layer. A mathematical notation is

adopted to record sets of paths from all outputs to related inputs.

Definition
A Path-order matrix is a matrix whose entries are positive integers (includes 0).
It has n rows, corresponding to the n source vertices. It has m columns,
corresponding to the m sink vertices. The vertices are associated with a SCG. If

pa, = k, then vertex i is in layer k on the path to sink vertex j.

An algorithm to derive the path—order matrix for a SCG is given in Appendix E.

The FORTRAN program implementing the algorithm is listed in Appendix C. An
example is presented to give some insight into the discussion of path matrix. For
convenience, the power transmission system in Figure 3.1 is considered again. Its SCG

is shown in Figure 3.10. From Figure 3.10 we see that the sink vertex set S is {1, 2,

44
3, 4, 5, 8, 9}, the source vertex set T is {7, 12}, and the number of vertices in the SCG

is 12.

Applying the path-finding algorithm to this graph, we identify all paths from the sink
vertex set to the related source vertex set, as given in Figure 3.11. Column headings in

italic denote layers.

 

 

 

 

 

 

a a a a a
taste...
I. t A. test
e a e e a as

Figure 3.11 An ex

 

 

ample of path identification

46
According to the information in Figure 3.11, a 12 by 7 path-order matrix can be

developed as follows:

--- Sink vertices ---

Vertices 1 2 3 4 5 8 9
1 1 0 0 0 0 0 0
2 0 l 0 0 0 0 0
3 0 0 1 0 0 0 0
4 0 0 0 1 0 0 0
5 0 0 0 0 1 0 0
6 2 2 2 2 2 0 0
7 4 4 4 4 4 O 0
8 0 0 0 0 0 l 0
9 0 0 O 0 0 0 1
10 0 0 0 0 2 0 0
11 3 3 3 3 3 2 2
12 L 4 4 4 4 4 3 3

 

 

From the computational point of view for equation solving, the solving process starts
from input vertices. In other words, it starts from the highest layer of a path. To obtain
the desired results we reverse the layer numbers of each path to make inputs the lowest

layer and outputs the highest layer. The final form of the path-order matrix is:

47

--- Sink vertices ---

Vertices 1 2 3 4 5 8 9
1 4 0 0 0 0 0 0
2 0 4 0 0 0 0 0
3 0 0 4 0 0 0 0
4 0 0 0 4 0 0 0
5 0 0 0 0 4 0 0
6 3 3 3 3 3 0 0
7 1 1 1 1 1 0 0
8 0 0 0 0 0 3 0
9 0 0 0 0 0 0 3
10 0 0 0 0 3 0 0
11 2 2 2 2 2 2 2
12 1 l 1 1 1 l 1

 

 

For example to find the output variable associated with vertex 5 we first solve the
equations associated with the vertices 7 and 12, then solve the equations associated with

the vertices 11, 10 and 6, as well as 5.

As we discuss in more detail in the next section, an important aspect of the path-order
matrix is that it provides the basis for a method to sort equations and get the Jacobian

in an efficient way.

48
3.4 Generation of Jacobian Status Matrix

For an implied set of explicit differential equations of the form

*1 = f: (1» u) (3.9)
the Jacobian is defined as follows:
J = [Jul (3.10)
where
_ 3f: 1 = 1, 2, n
“a .-12.,.. am

In the solution of nonlinear differential equations most numerical integration algorithms
make use of the local Jacobian repeatedly. One way the Jacobian can be estimated is to
use a difference approximation to the derivatives. This is relatively easy to implement
and quite general, but it is computationally costly. An increase in solution efficiency will

result if the cost of evaluating the Jacobian can be reduced.

A FORTRAN computer program I AC has been developed to generate the analytical state
equations of a system along with its system Jacobian matrix using the bond graph
representation of the system model (Hamilton, 1984; Sobhi, 1985). A symbol

manipulation technique was used in the program I AC. Availability of the Jacobian in

49

symbolic form increases the efficiency of the system analysis.

Here we introduce a status matrix associated with the Jacobian, SI = [SJ ii], whose
elements are either 0 or 1. In the calculation of a Jacobian, there are three possible types
for each entry: always zero, always constant, or a state-dependent or time-varying

function. Interpretations are made below.

If It is zero, then the j-th input does not affect the i-th output.

If I, is constant, the i-th output changes proportionally to the j-th input; those entries
have to be calculated only once.

If 11' is a function of the state or time, then the i-th output varies in response to the j-th
input according to the test state so that these entries have to be updated

continuously (at every time step).

Each entry of the status matrix of a Jacobian, 813, is defined as follows:

SJ-- = 0 : when 1,, is zero or constant,

81- = l : when 15 is a function of state or time.

First we consider models whose CGs are acyclic (SCGs). The idea is that if the paths
from each output to each input are identified, then we can trace the path to identify the

vertices along this path. Since each vertex refers to a function with local input and

50

output, by testing the type of functions we can obtain the status matrix of the Jacobian.
This idea has been implemented in the ENPORT package. In ENPORT the function for
each vertex in the CG derives from one of two sources: the standard function library or

user-defined subroutines.

 

 

 

////

 

 

 

 

 

 

 

Figure 3.12 Mechanical system and its bond graph model

Fact 1: For the i-th output, if there is no path to the j-th input, then J, is = 0 and
SJ, = 0.

Fact 2: For the i-th output, there exists a path to the j-th input. If the fitnction

51
types of the path vertices are all proportional, then J, is constant and SI,

= 0.
Fact 3: For the i-th output, the path exists to the j-th input. If any of the path
vertices has a non-linear fitnction type, then J, is a junction of state or

timeandSl,=1.

Table 3.4 List of system equations

 

 

 

# of Eqn Output Input Function type
1 El time SIN
E4 E1 SUM
E.2
E3
3 Q.2 F.2 INTEG
4 E.2 Q.2 ATT
5 R4 E.4 INTEG
6 E4 P.4 ATT
7 F.3 F.4 ASGN
8 F.2 F.4 ASGN
9 E3 F.3 DIODE

 

Let us consider a mechanical system and its bond graph model shown in Figure 3.12.
The system equations of this model are listed in Table 3.4 in abstracted forms. The
function types are all members of the standard library. A CG is constructed based on
these system equations (Figure 3.13). With the application of the modified BFS

algorithm to this SCG, the paths from P.4 and Q.2 to R4 and Q.2 can be identified

52
(Table 3.5).

 

 

 

 

 

 

Figure 3.13 CG of a mechanical system

Since only one node, 9 (representing E.3), contains a nonlinear function, the status

matrix associated with the Jacobian

 

 

 

am am‘
6R4 2
J = . 60 (3.5)
60.2 60.2
1 am 602 ,

53

is

s.r=[l 0] (3.6)

Table 3.5 Output-input paths

 

 

 

input/output 2 8
5 2+9+7+6-u5 8-.6-t5

 

 

 

Now we turn our attention to bond graph models which contain algebraic loops and/or
DNs. It has been pointed out previously that algebraic loops and DNs can be treated as
simple nodes (SCCs) for sorting purposes. To calculate the status matrix, we assume
that CG vertices corresponding to algebraic loops and DNs are nonlinear, since algebraic
loops have a group of equations which can not be decoupled explicitly in general and

DNs have an arbitrary set of differential-algebraic equations.

The method to simulate this type of bond graph model is similar to that used to simulate
standard bond graph models. The actual procedure includes the following steps:

- construct a CG for the bond graph model;

54
identify algebraic loops and DNs;

condense the CG to form a SCG;
find SCG paths from each output to each input;
check function types; and

generate SJ by the standard method.

4. MOVEMENT OF COMPUTATIONAL EFFICIENCY
FOR LSDSs DESCRIBED BY BOND GRAPHS

Much of the literature on simulation of LSDSs is concerned with numerically solving
large sets of differential-algebraic equations. A common purpose of much of the
research reported in the literature is to reduce the amount of computational work. The
engineering design problems with which we are concerned usually are not one-time
simulation runs. Typically they require repeated runs of the same model with varying
input and parameter conditions. For such problems there is great practical benefit to
reducing the computation time. In addition to the progress made on computer hardware
and operating system technology, two major innovations have had profound impact on
the study of LSDSs by computational methods. These are methods for model-order
reduction and methods that use multiple-timescales. We will not discuss model-order

reduction methods in this work. We will discuss multiple-timescales methods.

Thus far we have assumed implicitly that, given a set of differential-algebraic equations,
a particular numerical integration algorithm will be selected and used during the entire
simulation process. As was pointed out by Chua and Lin (1975), "under this assumption,
the step size for each time step may be optimized by choosing the largest possible value

of h for which the local truncation error remains bounded below the user-specified

55

56

maximum allowable error, and for which the algorithm remains numerically stable. For
large systems of equations, the amount of computation does not increase substantially
when the order of the algorithm is increased. Consequently, it often turns out to be more
efficient to vary both the order and the step size during each time step. " Here order

might refer to order of a Runge-Kutta algorithm.

Consider the system graph models we have discussed previously. It is desirable to
develop an efficient computational method to simulate them in a production mode (i.e. ,
for multiple solution runs). In this chapter we consider the increases in computational
efficiency achievable by selecting the integration algorithm and by selecting the step size
for each dynamic submodel. We shall refer to this as a multiple-integration-algorithm

method. We shall also include multiple-step-sizes as a tool.

4.1 Solution of Sorted System Equations

4.1.1 Submodels

Structural decomposition of a system is based on the graphical description and the
detailed model equations. Graph theory has been applied to CGs derived from the model
to sort the system equations and arrange them in a suitable calculation order. According
to the equation structure of each vertex in the CG we consider three types of nodes in

a system model. For illustration see Figure 4.1.

57

For given inputs and outputs submodels provide specific information about the status of
submodel equations. For instance, standard_node submodels have a set of explicit
algebraic equations, algebraic_loop submodels are represented by sets of implicit
algebraic equations, and DN submodels are described by a set of differential-algebraic
equations. The differences among the equation structure of these submodels are

significant. They will lead to different handling strategies at the solution stage.

 

 

 

 

 

 

 

 

Computational
Algebraic Differential-algebraic ebraic
m A18 loops

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4.1 Submodel structure of system model

58
It helps us to understand the characteristic of submodels by reealling how they are

constructed in the Modified ENPORT program. The construction of submodels occurs

in two ways:

( 1) User_constructed: Sets of differential-algebraic or algebraic equations are
group-coded in a FORTRAN file which is linked with the main program.
Each subroutine in this file determines a dynamic or algebraic submodel.
(2) Program_constructed: After the whole model is declared in Graph option
and equations relate to the model are defined in Equation option, the
algorithms embedded in the program are activated to identify simple nodes

and algebraic loops.

It is possible for a DN node to be included as part of an algebraic-loop set, based on its

algebraic input output structure. This is not the usual case.

4.1.2 Structure Decomposition

The above definitions of submodels can be extended to obtain a structure decomposition,
where a complete model is decomposed into an interconnection of submodels. In
general, each submodel interacts with the rest of a model in the same way as external

input, and output variables.

For numerical efficiency reasons it is often profitable to handle DNs and the rest of a

59
model separately. Recall the definition of DNs given in subsection 2.2.3. Detailed types

of equations were listed in Table 2.1. At this point we discuss DNs as part of a
complete system. A system may be decomposed into an interconnected set of submodels.

Each submodel contains exactly one DN, or it contains all of the remaining nodes.

Without loss of generality, let us consider a system with three dynamic nodes, as
depicted in Figure 4.2. For consistency in notation DNi (i > 0) is used to represent the

explicit dynamic nodes, while DNO is used to represent the rest of system.

Each submodel i has an input vector made up of two subvectors, namely, U, and U,,.
U, represents inputs from other submodels. U,, represents external (system) inputs. Each
submodel i has an output vector mode up of two subvectors, namely, Y, and Y,,. Y,
represents outputs to other submodels. Y,, represents external (system) outputs. Figure
4.2 shows that a submodel output component may not be feedback to its input. From

Figure 4.2 we can write the following equation sets:

DN,, has the equations

X0 = ¢1(xor ”or U”)
Y, = i.e.. U... U...) (4.1)
Yga = ¢w(xor an U”)

 

 

 

 

 

 

—— DNo

Yso
. ) 1 l
UO (x° YOI

DN1 [——

(X1)
Y1 I

.1: 0.02 7;,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4.2 Structure model

 

61
DN, has the equations

XI = ¢1(x19 ”1, Us!)

Y1 = *1(X1s U‘, U”) 4 2

Y8] = *;1(x1’ U1, U”) ( ' )
DN, has the equations

x2 = ¢2(x29 ”2, U“)

Y2 = *2(x2s ”2’ U“) (4.3)

Yd = *a(x2, U2, U“)

Furthermore, each U,, is composed of elements from Y,,, j #5 i.

The three structural submodels are connected with each other by following relationship,

 

 

”0 -W00 W01 W021 rYoI
”1 = W10 W11 W12 Y1
U 2 (W20 W21 W22, _Y2,

 

 

 

 

where W,- is a matrix of 0 and 1 elements. A "1" appears if an element of Y, appears

as a member of U,- else W is 0.

62

The numerieal solution of such a decomposed but interconnected system consists of
integrating sets of equations over a desired time interval. This is most commonly done
by discretizing the time interval into intervals marked by time-points t,, k=l,2,..., k.
Without loss of generality let us consider a time-point at the time t,=t,,. At this time

point, equations (4.1), (4.2) and (4.3) can be expressed as

xo (’0) = 4’0 (x000): U000‘): U,0(t0))
Y0 (t0) = *0 (X0(t0)9 U000): Uw(t0)) (4'4)
Y” (‘0) = 'w (X0(t0)r Uo(‘0)r Um(t0))

X1 (‘0) = ¢1 (X100). U1(t0)9 U,1(t0))
Y, (to) = 11,0000. 0,0,), U.,<:,)) <45)
Y“ (to) = *3] (X100). ”1(t0)9 U,1(t0))

X, (to) = 0,090.). 0,0,). U,,(t,))
190,) = macro). U,(:,). U, 0.)) (46)
ya (‘0) = *32(x2(t0)r U200): Ug(10))

According to the definition of a SCG as given in section 3.1, the starting vertices of the

CG consist of X, and U.,, and the ending vertices consist of X, and Y,,. Therefore, the

63

CG is a composition of static structure submodels at a certain time-point. The term
"static” is used here because the computations involved in reaching ending vertices are
algebraic. Numerical integration methods have to be used to obtain Xo(t,+,), X,(t,,,), and

x,(1,,,) from the values of x00), x,(1,), X2(t,), and U(t,) (i <i+1).

4.2 Multi-rate Solutions
4.2.1 Direct Solution Methods
In the case of a system graph with explicit integration implied everywhere, the state-

space and output equations of the model can be written as

x = F(X, Us, 1)

4.
Y, = 60!, 0,, t) ( 7)

Figure 4.3 indicates a procedure for organizing the system equations for solution. The

properties of differential-algebraic equations typically encountered in LSDSs are:

large: high order of X exists.

00'". G)
50‘. U)

sparse: most entries of the matrix are zero.

 

 

 

Input node equations

1

Construct the first level

computational graph
(06)

l

 

 

 

 

 

 

 

 

Identify DNs and algebraic loops

 

 

 

1

Construct the stande
computational graph
(see)

I

 

 

 

 

 

Generate the path-order matrix

 

 

 

 

Figure 4.3 Flow diagram of sorting system equations

 

65

- stiff: there are greatly differing time constants present.

One complication that can occur in some problems is that same components X, of the
derivative vector X appear implicitly in Equation (4.7). If we restrict each DN’s
equations to be explicit, however, this complication does not occur in system graph
analysis since the derivatives are explicit as well. We note that DN,, must also be explicit

in X. Thus, derivative causality in the bond graph part of DN,, is not permitted.

A wide variety of algorithms are available for direct integration of Equation (4.7). The

most commonly used ones are listed in Table 4.1.

Table 4.1 Integration methods for ordinary-differential equations

 

 

 

 

 

 

 

‘ Linear Explicit One-step Euler
; Explicit Runge-Kutta
. . . Adam Bashforth
Non-linear 39$: (“gimp Runge—Kutta—4
p tep Adams-Moulton
Stiff Implicit Multi-step Backward-Difference
Formulas

 

When all parts of a model have a similar time scale, it is usually possible to find an

efficient time step for solution that meets accuracy requirements.

4.2.2 Multiple Integrators

The main motivation to apply multiple integrators to solve the mathematical equations
listed in Equations (4.7) is the reduction of solution computational effort. In particular,
multiple integrators can be used to benefit models that have DNs, by taking advantage

of their distinct features.

In the prior discussion the differential equations for system graphs with DNs were
restricted to be purely explicit. To better show the advantage of the multiple integrator
method, the explicit restriction on DN equations is relaxed in this section. That is, any

DN may have implicit differential equations.

The basic idea of the multiple integrator approach is described in Figure 4.4. From the
computational point of view, there are two decisions to be made for each DN in the
system. They are (1) what step size to use for computation and synchronization, and (2)

what integrator to use.

Decisions about step size we refer to as MSS (Multiple Step Sizes).

Decisions about integration we refer to as MIA (Multiple Integration Algorithms).

For instance, the M88 method can be used in stiff but partitionable linear systems, which
contain some DNs with rapid dynamics compared to others. It is typically necessary to

use the smallest integration time step globally to obtain the solution of this type of model.

67
On the other hand, the MIA method can be applied when some DNs of a system have

large nonlinearities while others do not.

 

 

DNo

 

 

 

 

MSS

 

 

 

 

 

Main
controller

DN1

 

 

 

 

 

 

 

/

MIA

 

 

 

 

 

DN2

 

 

 

 

 

 

Figure 4.4 Diagram of multiple integration approach

Some solution properties of the DN equations (e.g. , linear or nonlinear, explicit or
implicit) can be derived from modeling phase. Some properties derived from values of
parameters (stiff or non-stiff) are often not known in advance of solution. Nevertheless,
a direct solution method can be used to make a test run, so as to get some insight into

the solution properties of various parts of the system.

68

Figure 4.5 shows a flow diagram for solving LSDSs using a multiple-integration
approach. In this diagram, the subscript index ”0" refers to equations not associated with
explicit DNs. In general, subscript "i” refers to DN,. Part 1 of Figure 4.5 shows the
initial setup. In block 3 ”to” denotes initial time for solution, "T" denotes final time,
”X,(t,,)" is the initial condition vector, and ”at,” is the step size for reporting for DN,.
Part 2 of Figure 4.5 shows the detailed procedure if a single integration (INTO) is used
for all DNs, but solution step size varies for each DN,, 120. Part 3 of Figure 4.5 shows

details if multiple integrators are used.

 

 

Read in nation solution
ordeICIrom SCG

 

 

 

l

 

 

Choose integrator (INTO)
for DN 0

 

 

l

 

 

Read in to, T, XS (t0),
and Ato

 

 

 

l

Select M88 or MIA
for each DN

.433 VIA
B

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4.5a Program flow chart, part 1

 

70

 

 

 

 

X, (‘91)
" X, (‘1)

 

 

 

 

 

A
l

 

 

 

 

 

 

 

 

 

 

 

 

 

 

I Evaluate U, (toI—.——
I

 

“0!. ti + Ate

 

Bum Y, (t,)

 

 

 

In

r—L—I 1- hi (Lu-«H

 

 

Evaluate x, (1,)

 

 

l

 

 

saverux,(t,.,)

 

 

l

 

 

’ W Y. (II): xlal)

 

 

 

 

 

 

 

 

 

 

 

 

Given X j (to), at,
J U
at, zato? I
.22..
at;

 

 

 

i-i+ll
moan”)

 

 

 

 

Solve for
xc (tit!)

 

 

 

 

 

 

 

_—-—IBvflmU,g
I

 

‘M' tl " “I
my,”

 

 

 

 

 

 

 

 

Evaluate X, (t,) XI 0 1.1)
I '1‘! (trot)

 

   

 

 

 

 

 

 

 

 

“01- I; + AI,
“thl (1,.)

l

m V. 0.1. x, (to

 

 

 

 

 

‘1' T7
yu

Stop

 

 

 

 

 

Figure 4.5b Program flow chart, part 2

 

71

 

 

B
l

Givm X j (to), at J - etc

I

(hemeintegtatcrfcr DN,
(INTj)

l

i-O

l

Evaluate U, (h) "‘

l

hor' t1 "’ At:
Us:rll~l'l‘,t11unlve)(j (rm)

I i-i+1

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

13mm: Y, (1.), x, (1,) MI". “Mi
Evaluate Y. 0.). X. 0.1

 

 

 

 

 

 

.! . °° Use 1mm Solve
f“ X. (‘9‘)

 

 

 

 

 

 

 

1,1?
I,-
Swp

 

 

 

 

Figure 4.5c Program flow chart, part 3

 

72
4.3 Software Implementation

A software design was made for implementing both the MIA and the M88 methods. The
host simulation software chosen to imbed the implementation was ENPORT. ENPORT
provides a well-tested base that supports bond graph and block diagram modeling and has
a suitable equation surfing procedure. We described earlier the implementafiun of DNs

within ENPORT.

In this dissertafiun the author did not implement the enfire mulfiple integrator approach
operafiunally. However, the M88 method is implemented. As an example, integrator
RK-4 was chosen to illustrate the M88 method and several subruufines for DNs were
coded in FORTRAN. These subroufines were then incorporated into the ENPORT

simulafion package.

The MSS program contains four subroufines. These subroufines are discussed briefly

below. Lisfings of these subroufines is provided in Appendix D.

1. Subruufine MCDDRV -- the main driver for MSS solufion phase.
The purpose of this subroufine is to control the selecfiun of cumputafion step sizes

for each DN.

2. Subruufine MCDINT -- control the integrafion process.

MCDINT uses the Runge-Kutta method to compute the state variables of DN

73

from a given inifial fime, t,,, to a final fime t, at DTMJ intervals. It stores results

at DTSTR intervals.

3. Subroufine RKMCD --- perfurrn integrafion with Runge-Kutta method.
The inputs to the subroufine include the DN equafions, computafion fime step,
current fime and the state variables at current fime. The subroufine provides the

values of state variables at the new fime.

4. Subruufine INIMCD --- get the inifial cundifions from user for a typical DN.

4.4 An Example

Next, an example is shown which uses the M88 method to obtain simulafion results. The
pusifion control model presented in Chapter 2 is considered again. The system graph,
which contains a set of standard nodes and connectors, is depicted in Figure 2.7. In this
example, the dynamic block nude DBNl is defined by the standard block atoms {SUM,
GAIN, INT}. The dynamic mulfiport node DMNl is defined by the standard mulfipurt
set {R1, 1A, 11, GY}. The MACRO is defined by {1B, R2, 12, 0A, C, 1C, I, R}. The
rest of the system (DN,) is defined by {MACRO, SE, SRC}. The model consists of
three submodels, according to the previous discussion of structure decomposifion. They
are DBN 1 , DMNl, and DN,. The submodel structure of this system is illustrated in

Figure 4.6.

74

 

 

 

SRC

 

 

 

 

 

 

DNO

 

 

 

 

.._-@

 

 

 

 

 

 

 

 

 

Figure 4.6 Submodel structure of system graph

More detailed infurmafion on the generafion of this parficular graph model with modified
ENPORT is contained in Chapter 2, in which all the required parameters and equafions
are specified. Given the detailed data, modified ENPORT proceeds to analyze and solve

the problem to obtain the fime response.

Equafions for each submodel can be derived according to the infurmafiun available to the

system. For DBNl, with x, = THETA, we get equafion (4.8)

75

dxlldt = wl
s2 = -2.0'lrxl + S] (4.3)
8 = x1

For DMNl, where x,= P.E2, we get equafion (4.9)

dx1/dt = -50.0=0:x1 + 5.1 - 1.0*F.M1
F.1 10.0*x, (49)
EM] 10.0 *x,

For the rest of the system (DN,), where x,= P.M2, x2: Q3, and x3= P.5, we get

equafiun (4.10)

dxlldt = -1.2:lrx,-100.01rx2 + E.MI

dledt = 4.01ml - 1t3
dxydz = 100.0*x2 - x, (4'10)
F.MI = 4.0 *x,

It is known from a frequency analysis of the entire system that the state variables in DN,,
show a high frequency component of same importance. The MSS methods is employed
so that a smaller step size is used to integrate this part of the model. The CPU fime is
recorded for integrafiun runs from 0 to 10 seconds when various step sizes for the

submodels are used. The fime responses are shown in Figure 4.7. They show strong

76
overall similarity of behavior. Comparing the results obtained with different step sizes,

we can see that the accuracy is within the bound of 0.45 96. The CPU fimes for sulufion

are listed in Table 4.2.

Table 4.2 CPU fimes for sulufion

        

Submodels Calculafion fime CPU fime
step size (second)

DNO 0.0125

 

 

 

 

 

Table 4.2 shows that the computafion cost for comparable accuracy is reduced by more
than 50% when at for DBN 1 and DMNl is increased by a factor of four. It should be
noted that for producfion runs, i.e. , repeated simulafions, the aggregate computafion fime
can be greatly reduced. This is a small size problem, but even so it shows the

possibilifies for a gain in solufiun efficiency with the M88 approach.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

SCALING
m7!
5.50r-01
0.00900
101
mos-01
owes-0)
0.00 0.20 0.40 0.60 0.80 1.00
TIME '10: 1
1.36880: mara— III --
TIN! STEP DT-0.0125
SCALING
nuns
mos-01
0.00900
111
Loot-01
oraor-or
0.00 0.20 0.40 0.60 0.80 1.00
TIME '10! I
LEGEND: THETA-—— HI...
TIRE STEP DT-0.025
SCALING
”(TA
,.,“‘°I ”107-... fi":
0.00000 :-‘ /
A. .‘ /
111 l .
Liar-01 .1 /
-s.m-01 5' ..
F [I I
: //f I l '*..
l/ .... " l
0.00 0.20 0.40 0.60 . 0.80 1.00

TIME '108 I
LEGEND: THE?A-—— H1 ..

TIHE STEP 3T-0.05

Figure 4.7 Time response

5. EFFICIENT CONIPUTATION OF ALGEBRAIC LOOPS

5 .1 Problem Description

Systems containing algebraic loops arise quite naturally in many applieafions. The
existence of algebraic loops in equafions of a physical system may not be detected unfil
the equafion surfing process starts in most convenfional simulafiun approaches. The

definifion of an algebraic loop is a set of algebraic equafions of the form

Y = G(X, U, Y) (5.1)

where it is not possible to reorder or solve the equafions into the modified explicit form

Y = H(X, U) (5.2)

After algebraic loops are idenfified during the process of surfing into SCCs in the CG,
the next important step is to obtain solufions to each SCC. Unfortunately, the task for
algebraic loops can not always be accomplished easily. The computafion usually

demands a large amount of computer-fime. There are two commonly used approaches

78

to reduce cumputer-fime.

(1)

(2)

5.2

Modify the system model to avoid implicit algebraic equafions. Burreto and
Leferre (1985) introduced two methods to handle this situafion: (a). imposing
restricfions in the set of admissible solufions; and (b). preparing a model to
simulate the system with explicit methods. These two methods require a lot of
mathemafical handling and model parfifioning. Granda (1984) proposed a method
that an algebraic loop can be broken by introducing a parasific physical element
into bond graph model. This method, however, may introduce sfiffness to

system differenfial equafions.

Use iterafive numerical methods to solve the algebraic loops directly. Meanwhile

try to improve the computafional efficiency within implicit solufions.

Iterative Methods

Consider the mathemafical representafion of an algebraic loop, as stated previously, and

repeated bellow:

(5.3)

At a certain fime t,, X(t,), U(t,) are known, and Y(t,) is to be fund. The principal steps

of using iterafive methods to solve Y (t,) are as follows:

0

80
(1) make an inifial guess of Y“’(t,) (predictor);

(2) use a pre-chosen formula for calculafing Ya’(t,) (corrector);
(3) 031‘:th 1' = G(Y“’(t,)) ' (KYOTO);
(4) check | r | s r,;

if I r | > r,, let Y“’(t,) = Ym(t,), go to (2);

if I r I S 1'11: Y“’(ta) = Ya’m). stop.

In step (4) r, is the desired accuracy.

When this iterafive method is applied to solving algebraic loops, the vector Y is an
iterafion vector F. Generally speaking, for nonlinear systems the higher the dimension
of Y, the higher the order of F. The higher the F dimension, the more computafion will
be needed. In the case of a linear system (for which an analyfical solufion can be
generated), however, this is not true. In general, can we find a minimum set of iterafion
variables for F to solve the iterafion problem for Y? If the answer is 'yes" , then how

can this be achieved? If the answer is "no” , then what are the altemafive methods?

Some research has been conducted on this subject. Rules have been developed to assign
causality to R-fields in bond graph models (Zhou, 1988). With these rules, the minimum
set of iterafion variables of an algebraic loop can be found. This approach, which has
been applied to bond graph models directly, has been limited by the structure of bond

graphs (e.g. , one-port or mulfiport, internal or external bond, juncfion structure). Such

81

approaches also do not cover loops involving blocks and signals.

In this chapter, the above quesfions are discussed from the perspecfive of directed

graphs. An algorithm has been developed by author to implement the new method.

5 .3 An Algorithm for Finding an Efficient Set of Iteration Variables

From the previous discussion, it is known that a SCC is a digraph representafion of an
algebraic loop (including the limit case of one equafion, a simple SCC). In order to
discuss the algorithm we introduce definifions of some elementary concepts and

terminology which are commonly used in digraph theory.

Definitiom:

Reachability: If a directed path leading from vertex x, to vertex x, exists, we say
that x, is reachable from x,.

Acyclic digraph: If a digraph has no cycle (i.e. , it does not contain two mutually
reachable verfices), then it is referred to as an acyclic digraph.

cyclic edge: A edge is cyclic if and only if it lies on a cycle.

A minimum independent set: A set of cyclic edges in a SCC is an independent set
if removing these edges leads the SCC to become
an acyclic digraph. A minimum independent set is

an independent set which contains the minimum

82
numbers of edges.

An example is given in Figure 5.1 to illustrate how a set of iterafion variables can be

obtained with the aid of digraph theory,

 

 

 

 

 

 

Figure 5 .1 An example of SCC

Assume that the SCC in Figure 5.1 is idenfified from the CG of a system model. It
represents an algebraic loop. The simplest way to solve this problem is to take all nodes
as starfing variables and begin the iterafion process. This is not an efficient method, but

it is organizafionally simple.

Assume that an inifial iterafion set F contains only node 6). An iterafion process cannot

83
be carried on. This is because informafion from both nodes 6) and ® is required to make

6) knowable. Therefore, it is natural to add node G into the iterafion set F so that the
algebraic loop is solvable. It will be discussed later, however, that the size of this

iterafive set is not minimum.

Now let us invesfigate this problem from a graph theory point of view. In this example,
thaeexisttwocycles. OneisO 0 0 Mauritania 0 0 0 0. Lotus

consider the first choice of iterafion set which contains node 6). Removing the edge
incident on node 6), say E1 or 132, the cycle ©+®+©+®+® is broken but the cycle
®+®+®+©intheSCC isleft unaffected. To maketheSCC solvable, onemethodisto
add another node to the iterafiun set, for example, node Q). Hence, after the edges
incident on G) and Q), i.e. , E4 and E1 or E2, are removed, both cycles are broken. The
digraph remaining becomes acyclic. The size of the iterafion set, however, is sfill not
minimum. Let us consider other possible choices of iterafive sets, namely, {1}, {3}, or
{4}. Starfing from any one of these three iterafion sets, variables on each node can be
solved. It is also expected that removing edges incident on the iterafion set, E3 or E4,
the SCC becomes acyclic. From the discussion in this example, we can reach a

conclusion as follows:

Idenfifying a minimum set of iterafion variables in an algebraic loop is equivalent

to finding a minimum independent set in the SCC associated with the algebraic

loop.

84

Finding a minimum independent set for SCCs is NP-complete (Even, 1979). However,

algorithms to find a near-minimum independent set can be developed.

The DFS (depth-first search) technique, a very useful algorithm for scanning a finite
graph, was introduced in Chapter 3. In the example in Chapter 3, DFS was performed
on a SCC to idenfify back edges (Even, 1979). An arbitrary vertex as a starfing node
in SCC is chosen and DFS is applied to obtain a set of back edges. Such a set of edges

is an independent set.

The following procedure describes a proposed algorithm for finding a near-minimum

independent set.

1) Start from each vertex of a SCC and apply the DFS algorithm repeatedly to
idenfify back edges as well as the number of back edges. For each scanning,
store the back edges and the number of back edges in V,- and N, respecfively.

2) Compare the N, to get the smallest one, say N,. The V, is the near-minimum

independent set.

An example is given in Figure 5 .2 to show the applicafion of the algorithm.

85

 

 

 

 

 

 

 

 

 

Figure 5 .2 An example for finding the near-minimum independent set

The results described above are listed in Table 5.1.

86

Table 5 .1 List of independent sets in the SCC

   
 
 
  
  

Independent set No. of Near-minimum ,
edges independent set '

 

 

{e4, e8, e5, e10}
{e1, e5, e8, e10}
{e3, e7, e8, e10}
{e2, e10}

{e3, e6, e8, e10}
{62. e9}

{e2, e10}

 

 

 

 

.N-vak-bh

 

Comparing all the independent sets, a near-minimum independent set can be obtained.
From the results listed in Table 5 .1 we know that the near-minimum independent set in

a SCC is not unique.

The algorithm to determine a near-minimum set of iterafion variables of an algebraic loop
has been developed based on the use of computafional graphs and SCC concepts. This
algorithm can be applied to the solufion of the algebraic loops. Reducing the number of
iterafion variables will contribute to improving computafional efficiency for large scale

non-linear dynamic systems.

A computer program to realize the algorithm is coded in FORTRAN; this file is listed
in Appendix C. This subroufine and some other subroufines which are used to sort

equafions are saved in one file and named SORTCG.FOR.

6. CONCLUSION

6.1 Summary andDiscusionoftheRcsults

A new tool to improve flexibility and generality in modeling LSDSs was defined and
developed. Two new system graph element types, the DB (Dynamic Block) and the DM
(Dynamic Mulfiport), were introduced. Each new node is defined by a set of
differenfial-algebraic equafions. A system graph simulafion environment containing the

new modeling tools has several advantages over one without it. They are:

(1) a complex subsystem can be represented in a compact graphical fashion, due to

the assignment of mulfiple equafions to a single node;

(2) the graphical complexity of a system model can be reduced further, since the new

node type can be incorporated in macroelements; and

(3) models of subsystems which have been developed only in equafion form can be
included as single nodes, thus avoiding the potenfially difficult task of expressing

the system as a set of standard nodes and their funcfions.

The CG (cumputafional graph) is constructed from the node equafions of system graphs,

87

88
which may include DB3 and DM3. A systemafic approach to the surfing of model

equafions for solufion was developed. For the SCG (standard CG), each algebraic loop
(SCC) and DN was idenfified as a specific vertex. Finding a SCC in a digraph is
equivalent to finding an algebraic loop in the system equafions. For DNs the relafionship

of connecfions with other verfices is algebraic and gets incorporated into the CG.

A path-order matrix was introduced, associated with the CG. The entries of the matrix
show the order to reach sink verfices from source verfices in the CG. An algorithm to
generate the path matrix was developed and explained. Two types of CGs were
considered, CGs containing no acyclic sub-digraph and C63 containing acyclic sub-
digraphs. With the help of the path matrix, the state equafions can be organized for

numerical solufion at setup in a highly efficient manner.

The path-order matrix also provided a way to evaluate the Jacobian status matrix. The
status indicates what entries of the Jacobian have to be calculated once and which have
to be calculated every fime step during the solufion. For large-scale non—linear systems,
reducing the calculafions of Jacobian in the whole simulafion process can make a valuable

contribufion to the increasing solufion efficiency.

Based on the new modeling tools developed in this work, some system computafional
aspects, such as graph-oriented decomposifion and mulfiple methods of integrafion, were

addressed. A construcfion of structure submodels was suggested. This model

89
decomposifion provided a method to deal with DNs and the rest of system in parallel.

An improved strategy for solving large sets of equafions involved in system graphs with
DNs was presented. The proposed strategy is based on mulfiple independent but
synchronized integrators. Two types of methods were proposed. They were the mulfiple
step size (MSS) method and the mulfiple integrafion algorithm (MIA) method. The

advantages of the computafional method are

(1) a complex system with different local dynamic properfies can be decomposed into
several dynamic submodels;

(2) improved solufion efficiency can be obtained by applying suitable integrafion
methods to different submodels; and

(3) improved solufion efficiency can be achieved by selecfing suitable solufiun step-

sizes to different submodels.

As another cuntribufion to improving the solufion of LSDSs, an algorithm was developed
to idenfify a near-minimum set of iterafion variables for algebraic loops. This algorithm
uses concepts of graph theory to search for a set of edges to break all the cycles in the
loop. A program implementafion of the algorithm was presented. The combinafion of
implemented new methods was shown to produce significant reducfion in solufion fime

for comparable accuracy in an example.

90
6.2 Suggestions for future research

(1) In the discussion of structuring issues in modeling, we made the assumpfion that
equafions of each DN are a fixed set which we cannot modify or perhaps even access in
detail. In that sense, a given DM may require a specific causal orientafion, like a Se or
Sf node does. But if the port of a DM has a R-type causality and the local environment
does not assign a specific causality under the SCAP, then the DM belongs to an algebraic
loop. Since the DM is described mathemafieally by a set of differenfial—algebraic
equafions, the integrafion has to be operated with each iterafion step. A computafional

iterafion method is needed to handle this situafion.

(2) The mathemafical models to describe the DNs are limited to explicit state
equafions in this work. Since some theorefical work has been done on the development
of Lagrangian bond graphs, from which Lagrange’s equafion can be derived, the author
recommends that the representafion of mathemafical models of DN s be extended to allow
Lagrangian form. Thus system graphs with DNs can be used for the formulafion

equafions of mofion for dynamic systems in a more general way.

(3) The Jacobian status matrix provides a reference rule for whether or not to
calculate each term of the Jacobian matrix at the several fime steps. Now, a complete
computafional scheme for calculafing the Jacobian should be implemented. A study
should be done to invesfigate how much the solufion efficiency can be improved by using

the status matrix.

91

(4) The software implementafion of mulfiple step size methods is sfill in its infancy
and therefore needs fime and work to mature into a reliable piece of software. A
program implementafion for the mulfiple integrafion algorithm method can be developed
according to the flow diagram shown in Figure 4.5 . These two pieces of software can

be combined and tested with many models to make it more user friendly.

(5) It is suggested that the path method and the solving order be used to derive a set
of symbolic system equafions. The descripfion could be output to symbolic manipulafion
programs. The advantage of symbolic manipulafion systems is that they allow engineers

to analyze systems both parametrically and numerically.

(6) It is suggested that effort be made to integrate the modeling, surfing and solving
algorithms for general equafiuns of system graphs, and the method of finding a set of
near-minimum iterafion variables for algebraic loops, into a simulafion framework. The
ENPORT software could be modified accordingly, to improve the overall efficiency and

make it a more powerful tool for engineers.

APPENDICES

APPENDIX A. BASIC SYSTEM GRAPHS

A.1 Block Diagrams

A block diagram is a graphical presentafion of equafion informafion. It is often used to
display a system model in a form that allows us to understand interacfions occurring

between the system’s elements.

A physical system is consist of a number of elements and input—output relafionships, each
of them can be represented by a funcfional block. The transfer funcfions of these
elements are usually entered in corresponding blocks, which are connected by arrows to
indicate the direcfion of the flow of signals. Note that signals can only pass in the

direcfion of arrows.

The diagrams in Figure A.1 represent some actual physical elements. They are an
electrical resistor (a), a mechanical spring (b) and a moving mass (c) driven by an -
external force. The Figure A.2 depicts block diagram presentafions of some physical
elements and mathemafical processes. Figure A.2 (a) represents a resistor with an input
v (voltage) and an output i (current). They are related by a constant l/R. Figure A.2

(b) represents a spring whose resisfing tensile force f is proporfional to its extension x

92

93

 

 

_._.a
—V k 15x

_i- f_.. m
o—Wfl f 7///////

(a) (b) (c)

 

 

 

 

 

 

Figure A.1 Diagram of actual physical elements

so that f = kx. Figure A.2 (c) shows the relafiunship between the input force to an
object and the output accelerafion. The governing funcfion, in this case, is Newton’s
law, f = ma. Not only can block diagrams be used to describe actual physical elements,
but also can they be used to display mathemafical processes. Two examples of this are

shown in Figure A.2 (d) and 2.2 (e). If we integrate an accelerafion, a, over fime, the

velocity v is obtained as v = fad! . Similarly, integrafion of velocity over fime

94

produces displacement x, x = [M . Thus, in a block diagram presentafion symbols

in boxes represent operafions that must be performed on inputs to obtain outputs.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.1/k 4.3;. k _Lf—-1/m —a~
(a) (b) (c)
3., INT —"» L INT x
(d) (e)

 

 

 

Figure A.2 Block diagrams of input-output relafion

To model a mulfiport system by considering each element which has force and velocity

as input and output, the number of directed lines connecfing each block will be two.

Table A.1 Basic building blocks used for modeling systems

       
     

Functlon

IDistributor

 

 

3’13“
33:“
Y3=u

 

 

IPuncfion u IE) Y Y a H”)

 

KU

+31“ LEI—q Y

 

tegrator "N Y Y = fUdt

 

 

no output equafion

 

U
ISignal Sink @
ISignal Source -Y

 

u,
IllSummer EM”

 

U
Transfer Funcfion Y

 

 

 

L.__.___._a._@__=

96

Therefore, block diagram provides us an explicit way to show the power flow paths.
There are eight basic building blocks are commonly used in block diagrams for modeling

systems. They are listed in Table A. 1.

A.2 Bond Graphs

A bond graph is composed of a set of basic mulfiport nodes. They are the atoms of bond
graphs. In this dissertafion, the term atom will be used due to its un-splitted property.
Table A.2 gives a list of the nodes used in bond graph. The first two atom types are
called dynamic nodes because an integral or a derivafive equafion describes these nodes.
The third atom type models the energy dissipafion, the fourth and fifth atom types model
external inputs. The last four atom types model juncfion structures which enforce a
power—conserving constraint. Consider a dynamic system in bond graph form. We can
parfifion a graph into these four major groups menfioned above. This idea is represented
in Figure A.3, where dots represent set of all bonds that join juncfion structure to a given
field. Bonds that connect fields to juncfion structures are referred to as external bonds,
and bonds that joint one element of a juncfion structure to another are referred to as

internal bonds.

Table A.2 Bond graph atom nodes

 

 

 

 

 

 

 

 

 

 

 

 

 

 

_
Name Node
Capacitance 0 i C
I

Inertance O : I

f
Resistance 0 : R

f
Effort Source 9 2 SE

f
Flow Source 2 SF

1‘
Transformer 9 t 3 III: 0.. 1A

f 1 f 3
Gyrator 9, 5 av 9; 5
f . fl
Junction 1 fl 3 f2 7 f3
I 121 + e2 + c3 = 0
Juncfion 0 e, = e2 2 e,
I ft t f: tfs = 0

 

 

 

98

 

 

 

Source

Se, Sf

 

 

 

 

 

 

 

 

 

 

C,I

O,1,TF,GY

R

 

 

 

 

 

Storage

 

 

 

I 1mcfion, Structure

 

 

 

Dissipafion

 

Figure A.3 Basic mulfiport fields

 

APPENDIX B. AN EXAMPLE OF CAUSALITY ASSIGNIWENT

A simple DMN is considered and its equafions are written in explicit algebraic-
differenfial format to illustrate the idea of assigning causality to DNs. Figure BI is used
as an example to explain the assignment of causality of a DMN. Figure B.1 depicts a
bond graph presentafion of a DMN and its associated equafions are presented in Equafion

B.l.

 

:é:DMN—"—'

 

 

 

Figure B.l DMN used to explain causality assignment

X=f(x.r.U)
i=3 U) (3'1)

Any input and output variable chosen for this DMN has to be one of the four definifions

listed in Table B.l.

100

Table B.1 Causality and input-output relafionship

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Cl f
s U- c, , r- ;I
1 , 2
——>|§| DMN .,
f1 91
, . U. . .1
1L. DMN "— s:
3‘ A]
Us f2 , Y- e
S! 2
ft ‘1
U= e, , r= I
31 f2

 

 

101

If there is no requirement for assigning inputs and outputs, the DMN becomes indifferent

to causality. In summary, the rules to assign causalifies to DMNs are

(1) use equafions as constraint to assign causalifies to each bond. Therefore,

causalifies are fixed like SE and SF nodes.

(2) If there is no constrain, treat assigned causalifies as R nodes.

Now, let us consider a system graph with DMNs. When the process of assigning
causality results in causal conflicts at bonds of DMNs, the model is regarded as being
ill posed. When the second rule is applied to assign causalifies to DMNs, that is, a
causality can be chosen arbitrarily, it results in an algebraic-differenfial loop. Simulafion
of this type of models needs special skill, which is beyond the coverage of this

dissertafion.

APPENDIX C LISTING OF SORTING SUBROUTINES

 

CFILE:SORTCG

C

C---- PURPOSE: Sorting procedures for solution module.
C

C--- CONTENTS:
REDCMG Redefine the computational graph based on SCC
listing
REDATA Called by REDCMG and ALYSCC
ltempararily define SCC datal
MBFSLB Find the path from given output to
related inputs
lDENUY Called by MBFSLB. Identify the input nodes
and output nodes in new C6 in such an order
that [Xlil. UI' and [Xmil/dt, Yl'.
Identify the type of each node.
GENJAC Generate the Jacobian Status Matrix
ALYCMG Generate the sub-computational graph of SCC
NSTRIN Count the length of a string
WRTSTR Write a string on screen
ALYSCC Seach a near-minimun set of variables as the
initial guesses of the solution for each
algebraic loop
FINSRT Final sort for solution efficiency

OOOOOOOOOOOOOOOOOOO

C --- Last Modification: Aug. 12, 1991. YyW
CEOFH:SORTC€}

 

C>>>>>

C
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C

C Wriiten by: Yanying Wang, 04/10/91

C Last change: Yanying Wang, 04/12/91

 

C
C
SUBROUTINE REDCMG
C
C-—- PURPOSE: Redefine the computational graph including Strongly
C Connected Components. Replace each SCC by a single node.
C The modified computational graph should be ready to be
C applied the modified Breadth-First-Search algorithm.
C

102

103

C--- INPUTS: Computational graph data base,

C Strongly Connected Components data base.

C

C-- OUTPUTS: Modified computational graph (NE, EOTlil, ElN(il )
C

INCLUDE 'SIZEBK.CBK'
INCLUDE 'BFSMBK.CBK'
INCLUDE 'CPGFBK.CBK'

INTEGER l0, I1, J. LWK
CHARACTER DFILE’B
C
EXTERNAL REDATA
CREDCMG............................................................
C
WRITEI6,'I3X,"Enter your input data file name please:"l'l
READIS,'(A6)'I DFILE
OPENIUNIT = 1 ,NAME = DFILE/I'.CMG',FORM = 'FORMATTED',STATUS = 'UNKNOWN’I
OPENIUNIT = 3,NAME = DFILE/I'.RCM',FORM = 'FORMATTED',STATUS = 'UNKNOWN'I
C
l = 0
110 l= l+1
READI1,100,END=99) NE, EOTlIl, ElNlll
C TYPE',EOT(ll,EINlll
GOTO 110
99 CONTINUE
100 FORMATI3I4)
C
TYPE‘, 'NE=',NE
CALL REDATA
TYPE','NSCC=', NSCC
C
C--- Rename the EOTlil with the new node index in SCC data base.
C
DO 10 l0= 1, NE
DO 20 ll =1, NSCC
L1 = SCCPllI)
L2: SCCPlll +1)- 1
DO 25 J = L1, L2
lFlEOTllOl.EO.SCCL(Jll THEN
NEWEOll0l= l1
GOTO 10
ENDIF
25 CONTINUE
20 CONTINUE
10 CONTINUE
C
C--- Rename the ElNlil with the new node index in SCC data base.
DO 30 l0= 1, NE
DO 35 ll =1, NSCC
L1 = SCCPll1l

104

L2= SCCPlll +1l-1
DO 40 J = L1, L2
lFlElNllOl.EO.SCCL(J)l THEN
NEWElll0l= l1
GOTO 30
ENDIF
40 CONTINUE
35 CONTINUE
30 CONTINUE
C
C-- Exclude the edges having identical input and output nodes.
C
LWK == NE
DO 45 l= 1.NE
50 lFlNEWEOlIl.EO.NEWEIllll THEN
NEWEOlll= NEWEOILWKI
NEWEIlll= NEWEIILWKI
LWK: LWK-1
GOTO 50
ENDIF
lFll.EO.LWKl GOTO 55
45 CONTINUE
C
55 CONTINUE
NOE= LWK
NON = NSCC
C
C-—- I did it!
C
C--- Write the new CMG data into a file.
WRITEl3,'(” The new computational graph:"l'l
WRITEl3,'l3Xl'l
WRITEl3,1010l NON, NOE
WRITEl3.'l3Xl'l
WRITEl3/l" NE EO El”l'l
WRITEI3.'I3XI'I
DO 65 l0== 1.NOE
WRITEl3.300l l0. NEWEOIIOI,NEWEIIIOI
65 CONTINUE
300 FORMATl3x,3l6l
1010 FORMATlZX,’ NON =', l4,'NOE =’,l4l
C

WRITEl6,’(" New computational graph generated."l'l
C RETURN

END
CEND:REDCMG<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
C Wriiten by: Yanying Wang, 04/12/91

C Last change: Yanying Wang, 04/12/91
C ..........................................................................

 

C
SUBROUTINE REDATA
C
Cu- PURPOSE: Prepare the Strongly Connected Component data for the use
C of redefining new CMG.

C Called by REDCMG.FOR
C

C-- OUTPUTS: NSCC, SCCP. SCCL
C

INCLUDE 'SIZEBK.CBK'
INCLUDE 'CPGFBK.CBK'
INCLUDE 'BFSMBK.CBI('

C
CREDATAooeeooeeeeeeoeoeeoeoeoonnoooneonecooeooooooeoeoooeoeeoecocoon
C
NSCC= 12
SCCPlll= 1
SCCPl2l= 2
SCCPl3l= 3
SCCPI4I= 4
SCCPISI= 5
SCCPl61= 6
SCCPI7I= 10
SCCPl8l= 11
SCCPl9l= 12
SCCP110l= 13
SCCPl11l= 14
SCCPI12l= 15
SCCPl13l= 16
C
SCCLIII = 2
SCCLlZl = 4
SCCLl3l = 6
SCCLl4l = 8
SCCLISI = 11
SCCLl6l = 3
SCCLl7l = 5
SCCLl8l = 7
SCCLl9l = 9
SCCLI10I= 1
SCCLl11l= 10
SCCLl12l= 12
SCCLI13I= 15
SCCLll4l= 13
SCCLl151= 14
C

RETURN
END

106

CEND:REDATA<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C

C Wriiten by: Yanying Wang, 02/18/91

C Last change: Yanying Wang, 04/15/91

 

C
C
SUBROUTINE MBFSLB
C
C»- PURPOSE: Find the path from given outputs to related inputs, label
C the nodes with level number so that the output can be
C obtained by tracing the level index starting at 1.
C The idea of the algorithm comes from Breadth-First-Soarch.
C
C--- INPUTS:
c--- OUTPUTS:
C
INCLUDE 'BFSMBK.CBK'
INCLUDE 'SlZEBK.CBl<'
INCLUDE 'SOLNBK.CBK’
C
PARAMETER (NOE =12,NON =12)
CHARACTER DFILE'B
INTEGER lAlNONl. IB(NON), lClNOE), OUEINOE)
C
EXTERNAL GENJAC, lDENUY
C
CMBFaBOOOOOOOO00......OOOOOOOOOOOOOOOOOOOQOOOOGOOOOOOO0.06.000000...
C
0... Road in node number from computatinal graph
C
WRITEl6,'l3X,”Enter your input data file name please:")')
READ(5,'(A6l'l DFILE
OPENIUNIT =1,NAME = DFILE/I'.CMG'.FORM = 'FORMA'I‘I’ED’.STATUS = 'UNKNOWN'I
OPENlUNIT = 2,NAME = DFILE/l' .OUT',FORM = 'FORMATI'ED',STATUS = 'UNKNOWN')
C
I = 0
10 I = l+ 1
READI1,100,END =99) NOE, NEWEOll), NEWEIII)
GOTO 10
99 CONTINUE
100 FORMATl3l4l
C
DO 8 l= 1.NON
DO 7 J = 1.NON
LEVlI,J) = 0
7 CONTINUE
8 CONTINUE
C

CALL lDENUY

107

C
DO 500 IT= 1,lO
KO= IOVRIITI
DO 20 l= 1,NON
lAIll= 0
lBlll= 0
20 CONTINUE
DO 30 I= 1,NOE
OUElll= 0
lClll= 0
30 CONTINUE
C
C--- Fill IA, ID, lC arrays
C
IPT= 1
lAlKOl= lPT
DO 35 l=1,NOE
lFlNEWEllll.EO.l(Ol THEN
lCllPT)= NEWEOIII
lPT = lPT+1
ENDIF
35 CONTINUE
lBlKOl= lPT
C
lOUE= 1
lLEV= 1
DO 40 J = 1.NOE
JJ= ICIJI
IFllJJ.NE.0l.AND.llA(JJ).E0.0)l THEN
KO= JJ
lAlJJl= lPT
DO 45 I: 1,NOE
lFlNEWEllll.EO.KOl THEN
ICIIPTl= NEWEOIII
IPT= lPT+1
ENDIF
45 CONTINUE
IBIKOl= lPT
lFIlAlKO).EO.lB(KOll THEN
IBlKOl= 0
OUEIIOUEI= KO
IOUE= IOUE+1
LEVIKO,IT)= 1
ENDIF
ENDIF
40 CONTINUE
C
l2= IOUE
DO 50l1= 1,NOE
lFl|2.GT.l1l THEN
KO= OUEllll

108

CC II = II +1
DO 55 J =1, NON
lFllBlJl.GT.lAlJll THEN
J2= lBlJl-l
DO 60 J1 =lAlJl,lB(J)-1
lFllClJl l.EO.KOl THEN
lClJ1) = ICIJZI
lBlJl= J2
ENDIF
60 CONTINUE
ENDIF
lFlllAlJl.EO.lB(Jll.AND.llBlJl.NE.0)l THEN
IBIJl= 0
LEVlJ,ITl = LEVlKO,lTl + 1
OUElI2)= J
l2 = l2 + 1
ENDIF
55 CONTINUE
ENDIF
50 CONTINUE
C
500 CONTINUE
C
C--- Finish all of the paths searching...
C
WRITEIZ,’(" The output nodes are:")')
WRITE(2,200l llOVRlll.l=1, lOl
WRITEl2.'(3Xl’l
WRITE(2,'(" The input nodes are:")')
WRITEl2,200l (llVRlll,l=1, llNl
WRITEl2.'l3Xl'l
WRITEl2,'(" The output-input path index matrix:")')
IOP= IO + 1
DO 65 l0= 1,NON
WRITEl2,200l l0. (LEVll0.Jl. J =1,IOl
65 CONTINUE
200 FORMATI20l4l
C
CC CALL GENJAC
RETURN
END
CEND:MBFSLB<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C
C
CGENJAC> > > > > > > > > > > > > > > > > > > Last changed: 03/25/91 wa
C
SUBROUTINE GENJAC
C
C-~- Purpose: Find Jacobian Status Matrix. The difinition of
C Jacobian matrix here is JAClij) = DDXlillDXIII.
C

109

C-- INPUTS: LEVIIJI. output-input path index matrix

C NO, number of output variables

C NlN, number of input variables

C

C--- OUTPUTS: JACSTAli,i) =0 if JACli,il = 0

C JACSTAli,j) =1 if JACli.il -- something else
C

INCLUDE 'SIZEBK.CBK'
INCLUDE 'SOLNBK.CBK'
INCLUDE 'BFSMBK.CBK'

INTEGER NODEP, JACSTAI3,8l
LOGICAL LINFLG
C
COOOOOOOOO9......GOO...OO.9.0.0.900...0.0000000000000000000000000000000
C
TYPE', 'IO=',IO, 'llN =', "N
OPENlUNIT = 9,NAME = 'JACSTA.OUT',FORM = 'FORMATTED',STATUS = 'UNKNOWN'I
DO 20 l= 1,IO
DO 15 J = 1.IlN
JACSTAII,J) = 1
15 CONTINUE
20 CONTINUE
C
C--- Find the fixed zeros (only take first NXI terms).
C
DO 30 l= 1,NXl
DO 25 J = 1,NXI
NODEP = llVRlJl
lFlLEVlNODEP,ll.E0.0l JACSTAll,Jl = 0
25 CONTINUE
30 CONTINUE
C
DO 35 l=1, NXI
WRITE(9.10001 (JACSTAII,JI. J =1,NX|)
35 CONTINUE
1000 FORMATll X,8(1X,I2ll
C
RETURN
END
CENDzGENJAC<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C
C
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
C Wriiten by: Yanying Wang, 04/10/91
C Last change: Yanying Wang, 06/24/91
C

 

C
SUBROUTINE lDENUY
c .

110

C-- PURPOSE: Identify the starting vertices from computational graph

C in such an order IIVRIiI = [Xli), U]'.

C Identify the ending vertices from computational graph
C in such an order IOVRli) = [DXliI/DT, Yl'.

C

C»- INPUTS:

C--- OUTPUTS:

C

INCLUDE ’BFSMBK.CBK'

INCLUDE 'SIZEBK.CBK'

INCLUDE 'SOLNBK.CBK'
C

O'DENUYOQO..OOOOOOOOOOOOGOO0.0000000000.09.09...OOOOOOOOOOOOOOGOOOOO

C
NON = NSCC
DO 6 l= 1.NSCC
IIVRIII= 0
6 IOVRIII= 0
IO = 0
IIN = 0
C
On Identify the inputs from Computational Graph.
C--- Count the state variables first.
DO 50 1= 1. NXI
MP= XIXlII
DO 55 J=1, NSCC
JJ = SCCLISCCPIJII
IFIVOTPIJJI.EO.MPI THEN
C--- State variable found here

VTYPIJI= 'IN'
C
IIN = IIN + 1
IIVRIIINI= J
GOTO 50
ENDIF
55 CONTINUE
50 CONTINUE
C

C--- Then count the other starting vertices.
D019I=1, NSCC
NEOS = SCCLlSCCPlIII
DO 21 J =1,NEWNE
lFII.EO.NEWEIIJII GOTO 19
21 ‘ CONTINUE
IFIFNIZCIFNTPINEOSII.NEO.’INTEG'I
(3... System input found here
VTYPIII= 'U'
C
IIN = IIN + 1
IIVRllINI= I

111

19 CONTINUE
C
C--- Identify the outputs from the computational graph.
C-- Put the derivative of state variables at the first.
DO 18 I= 1, NXI
MP= DXIXIII
DO 19 J= 1, NSCC
JJ = SCCLISCCPIJII
IFIVOTPIJJI.EO.MPI THEN
C-- Derivative of state variable here

VTYPIJI= 'OU'
C
IO= IO+ 1
IOVRIIOI= J
GOTO 18
ENDIF
19 CONTINUE
18 CONTINUE
C

C--- Append the other outputs.
DO 5 I=1, NSCC
DO 2 J =1,NEWNE
IFII.EO.NEWEOIJII GOTO 5
2 CONTINUE
DO 7 K= 1, NXI
IFII.EO.IOVRIK)I GOTO 5
7 CONTINUE
C--- System output found
VTYPIII = 'Y'
IO = IO +1 .
IOVRIIOI = I
5 CONTINUE
C
C--- Identify the algebraic loops
DO 60 I=1, NSCC
J = SCCPII+1I- SCCPIII
IFIJ.GT.1I VTYPIII = 'AL'
60 CONTINUE
C--- Rest of vertices must be simple
DO 65 I=1, NSCC
IFIVTYPIII.EO.'###’I VTYPIII ='SIM'
65 CONTINUE
C
C--- MACRO node identification is not ready yet II
C
RETURN
END
CEND:IDENUY<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C
C
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

112

C
C Wriiten by: Yanying Wang, 06/04/91
C Last change: Yanying Wang, 06/23/91

 

C
C
program ALYCMG
C
C"- PURPOSE: Identify each SCC subgraph from Computational Graph. The
C SCC sub-graph is different with SCC in that it contains
C not only vertices also edges.
C Called by ANAALG.FOR
C

C--- INPUTS: Old computational graph data base and SCC data base.
C
c..- OUTPUTS: NOTLIJ with index of algebraic loop defined

C ALNEINALPSI number of edges in each algebraic loop
C ALEOII,JI leaving vertices in Jth Al-Ioop
C ALEIII,JI ariving vertices in Jth AI-Ioop
C NALPS number of algebraic loops
C
INCLUDE 'SIZEBK.CBK'
INCLUDE ’CPGFBK.C8I('
INCLUDE 'BFSMBK.CBK'
C

INTEGER I. J .K . NOTLIMAXVOTI

D INTEGER ALEOI10,10I, ALEIIIOJOI, ALNEIIOI. NALPS
CHARACTER MESSAGE'40. DFILE‘B
LOGICAL OKAY

C
EXTERNAL REDATA, WRTSTR
C
CALYCMGOOOOOOOOOOOOOOOOOOOOOOOOO0.00.0.0.9...OOOOOOOOOOOOOOOOOOOOOO
C
C--— Initialize
DO 12 I=1, NE
NOTLIII = 0
DO 14 J=1, MAXALP
ALEOII.JI= 0
ALEIII,JI= 0
ALNEIJI = 0
14 CONTINUE
12 CONTINUE
C
WRITEl6.’I3X,"Enter your input data file name please:"l’I
name/(Aer) DFILE
OPENIUNIT =1,NAME = DFILE/I'.CMG'.FORM = 'FORMA'ITED',STATUS = 'UNKNOWN'I
C
l = 0
5 l= I+1

READI1,100,END=99I NE, EOTIII. EINIII
GOTO. 5

113

99 CONTINUE
100 FORMATI3l4I
C
CALL REDATA
C--- Identify the algebraic loops
NALPS = 0
DO 10l =1. NSCC
II = SCCPII+1I - SCCPIII
IFll1.GT.1I THEN
IFINALPS.EO.MAMAII THEN
CALL WRTSTRI' "' Too many algebraic Ioops'l
OKAY: .FALSE.
CC RETURN
ENDIF
NALPS= NALPS+ 1
DO 20 II a SCCPIII. SCCPII+1)-1
N = SCCLIIII
NOTLINI= NALPS
20 CONTINUE

ENDIF
1O CONTINUE
C TYPE’, 'NALPS=', NAPLS

C
C--- Identify the sub-graph which including the algebraic loop
DO 25I=1, NSCC
ALNEIII= 0
25 CONTINUE
C
DO 40 K= 1, NALPS
DO 30 I= 1, NE
N1 = EOTIII
N2= EINIII
IFINOTLINI I.EO.I<.AND.NOTLIN2I.EO.KI THEN
ALNEIKI= ALNEIKI+ 1
ALEOIALNEIKI,KI= N1
ALEIIALNEIKI,KI= N2
ENDIF
30 CONTINUE
40 CONTINUE
C
C--- Print the results on screen
DO 50 I=1, NALPS
DO 60 J = 1, ALNEIII
WRITEl6,1020) J, ALEOIJ,II, ALEIIJJ)
60 CONTINUE
50 CONTINUE
1020 FORMATl2x,3l5I
C
STOP
END
C

114

CENDIALYCMG<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C
C
CNSTR'N e90990eeeeoeooeeeeeeeeeeeeeeeea... LOST 0.18009110/02/90 Y.Wang
C
FUNCTION NSTRINISTRING)
C
C--- NSTRING FINDS AND RETURNS THE NUMBER OF CHARACTERS IN
C STRING WHICH PRECEDE ANY TRAILING BLANKS. INTERNAL BLANKS
C ARE COUNTED AS CHARACTERS.
C
CHARACTER’I'I STRING
LOGICAL BLANK
INTEGER NSTRIN, LENGTH

INTRINSIC LEN
...NSTRIN ..............................................................

GET ACTUAL LENGTH OF STRING AND INITIALIZE BLANK

00000 0

LENGTH = LENISTRINGI

BLANK= .TRUE.
C
C--- SCAN BACK FROM END OF STRING FOR FIRST NON-BLANK CHARACTER
C

10 CONTINUE
IFISTRINGILENGTH:LENGTHI.NE.' 'I THEN
BLANK = .FALSE.
ELSE
LENGTH = LENGTH-1
EN DIF
IFIIBLANKI.AND.(LENGTl-I.GT.OII GO TO 10
C
NSTRIN = LENGTH
C
RETURN
END

CEND:NSTRIN<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C
C
CWRTSTReeeeeeeeeeeeeeeee99090009909099.e. Last change: 10,04/90 YWANG
C

SUBROUTINE WRTSTRISTRINGI
C

CHARACTER'I’I STRING

INTEGER LSTRNG. NSTRIN

EXTERNAL NSTRIN

C. . .WRTSTR ................................................................

115

C
LSTRNG = NSTRINISTRINGI
IF ILSTRNG.GT.OI THEN
WRITEI6.'IAI'ISTRINGII :LSTRNGI
ELSE
WRITE16.'I3XI'I
ENDIF
C
RETURN
END

CEND:WRTSTR<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

C

C _
C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C

C Wriiten by: Yanying Wang, 06/19/91

C Last change: Yanying Wang, 06/24/91

 

C
C
SUBROUTINE ALYSCC
C
C--- PURPOSE: For each complex SCC, apply Depth-First Search on it and
C find the minimum number of back edges. Those variables on
C the back edges are starting nodes for solving the
C algebraic loops.
C
C-- INPUTS: Computational graph structure of SCC.
C
C--- OUTPUTS: NOB number of back edges
C BIXII.JI name of outword vertices related to the
C back edges for the Jth AI-Ioop
C
C-- NOTATION: KIV), DFS index of vertex V
C FIV), father of vertex V
C LIV), lowerpoint of vertex V
C VONSIVI, location of V on stack S (=0 if not on S)
C ETIEI, edge type of edge E
C
INCLUDE 'SIZEBK.CBK'
INCLUDE 'CPGFBK.C8I(’
INCLUDE 'BFSMBK.CBK'
C
INTEGER LS. VX, EX, V, INI, U, BX, BEGXIMAXALP)
INTEGER KIMAXEONI, FIMAXEONI. LIMAXEONI. SIMAXEONI
INTEGER ET IMAXVOTI, VONSIMAXEONI
LOGICAL OKAY
C
CALYSCCO9.9.9.0...09.9909999009999990.0.0.9..OOOOOOOOOOOOOOOOOOOOOOO
C

INI= 0
DO 5 l= 1, NALPS

116

NOBIII = 0
DO 6 J =1. NE
6 BIXIJJI = 0
5 CONTINUE
C .
c-.. Repeatly perform DFS by starting from each vertex
DO 999 IO=1. NSCC
BX = 0
11 = SCCPIIO+1I - SCCPIIO)
IFII1.GT.1I THEN
lNl= INI +1
DO 995 VX= SCCPIIOI. SCCPIIO+1I-1
C--— Initialize
DO 10 I=1, I1
Klll = 0
Fill = 0
Llll = 0
VONSIII = 0
10 CONTINUE
DO 15 EX=1, ALNEIINII
15 ETIEXI= 0

LL= 0
LS = 0
V = SCCLIVXI

100 LL= LL+ 1
KW) = LL
LlVI = LL
LS= LS+ 1
SILSI= V
VONSIVI= LS
1 10 CONTINUE
C
C--- Check unused incident edges
DO 30 EX= 1, ALNEIINII
ETX= ET IEXI
IFIALEOIEXJNII.EO.V.AND.ETX.E0.0I THEN
U= ALEIIEXJNII

GOTO 35
ENDIF
30 CONTINUE
GOTO 200
35 CONTINUE

IFIKIUI.E0.0I THEN
C--- Tree edge here
ETIEXI = 1
FIU) = V
V= U
GOTO 100
ENDIF

117

IFIKIUI.GE.KIVII THEN
C--- Forward edge here
ETIEXI= 2
ELSE
IFIVONSIUI.GT.0I THEN
C--- Back edge here
BX = BX + 1
LIVI= MINILIVI.KIUII
ETIEXI= 3
C--- Relate the edge to the outward vertex
BEGXIBXI= ALEOIEXJNII
ELSE
C-—- Must be a cross edge
ETIEXI = 4
ENDIF
ENDIF
GOTO 110

200 IFIFIVI.GT.OI THEN
LIFIVII = MINILIFIVII.LIVII
V= FIVI
GOTO 1 10
ELSE
DO 40 NIX= SCCPIIOI. SCCPIIO+1I- 1
IFIKINIXI.E0.0I THEN
V = NIX
GOTO 100
ENDIF
4O CONTINUE
ENDIF

IFIBX.GT.NOBIINIII THEN
DO 50 I = 1,BX
50 BIXIIJNI) = BEGXIII
NOBIINII = BX
ENDIF
995 CONTINUE
ENDIF
999 CONTINUE
RETURN
END
CEND:ALYSCC<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C
C

APPENDIX D. LISTING OF MSS CODE

 

CFILE:MCDSOL ##CDSOL

 

C
Cm PURPOSE: Direct the system to obtain a solution of Macro_dynamic_
C node.

C

Cm CONTENTS: MCDDRV main driver for solution phase
C MCDINT controls integration

C - RKMCD oversees Runge-Kutta integration
C INIMCD sets initial conditions

C

Cm INDEX:

C INIMCD

C MCDDRC

C MCDINT

C RKMCD

C .

Cm Last revision: January 27,1992. Y-Y.WANG

C

 

CEOFHzMCDSOI:
CMCDDRV>>>>>>>>>>>>>>>>>>>>>>>> Last Change: 01/23/92 wa
C SUBROUTINE MCDDRVIDTCALC,NCALCI

Cm PURPOSE: Controls the selection of computation steps for M_C_D

Cm INPUT: DTSTR, storage time interval for whole system

C NSAV, number of storage intervals

C

Cm OUTPUT: DTM, computation time interval for M_C_D
C MCALC, no. of steps per system interval

C SOLMCD. flags to see the results on screen

C OPT4M, option of selection

C

INCLUDE 'SIZEBK.CBK'
INCLUDE 'SOLNBK.CBK'
INCLUDE 'MCDCBK.CBK’

CHARACTER STRING'70. FNAM'B. CH1'1, FNIZC'B

LOGICAL PROCFG, FULL, NEWLIN, ENDLIN, E7YORN
INTEGER INDEX. I. NC2I. NSMCD

118

119
REAL DTCALC. LO. HI, DTMT, DTMTEM, ABSDTM, DTHI, DTLO

EXTERNAL BLNKLN, WRTSTR, PROMPT, GET ANS, INVOPT
EXTERNAL GET IN, CONTUE, E7YORN, INIMCD. FNI2C
INTRINSIC INDEX
C
DATA OPT4M/‘T'l
C
COOOMCDDRVOOOOOOOOOOOOO00.0.9.0...0....OQOOOOOOOOOOOOOOOOOOO0......
C
OVFLCH = .FALSE.
FULL = .TRUE.
TIMADV = NCALC'DTCALC
C
C---- Select the integration method
50 CONTINUE
DO 20 N = 1.NONS
FNAM = FNI2CIFNTPINII
IF IFNAMI1:4).EO.'ZZDS'I THEN
C
C---- Check if valid M_C_D-type name IZZDSOI ZZDS99I
NF= NCHARSIFNAMI
IF INF.EO.5I THEN
CH1 = FNAMI5:5I
FNAMI5:6I = '0'//CH1
NF= 6
ENDIF
IF INF.EO.6I THEN
N1 = INDEXI'OI23456789',FNAMI5:5II
N2 = INDEXI'OI 23456789'.FNAMI6:6II
IF IN1.E0.0.0R.N2.E0.0I THEN

NC2I= 0
ELSEIF IN1.EO.1.AND.N2.EO.II THEN
NC2I= 0
ELSE
NC2I = 10‘IN1-1l +IN2-1l
ENDIF
ENDIF
CALL BLNKLN
WRITEISTRING,1024I FNAM
CALL WRTSTRISTRINGI
1024 FORMATI' "",A6.""’I
C
C -------- Ask user for the number of state vbl
CALL BLNKLN

NSMCD = NXMCDINCZII
WRITEISTRING.1025I NSMCD
1025 FORMATI' Enter the number of state variables I'.
1 l2,’):')
CALL PROMPTISTRINGI
NEWLIN = .TRUE.

120

CALL GETININSMCD.O.20,NEWLIN,ENDLINI
NXMCDINCZII = NSMCD

C
C-------- Ask user for the initial conditions
C
CALL INIMCDINCZII
C

123 DTMINC2II= DTCALC
CALL BLNKLN
WRITEISTRING,1030I DTMINCZII
1030 FORMATI' Enter the calculation interval I’,
1 1PE12.4,'I:'I
CALL PROMPTISTRINGI
LO= DTSTR/10000.0
Hl= DTSTR
NEWLIN = .TRUE.
DTMT= DTMINCZII
CALL GETRLIDTMT,LO,HI,NEWLIN,ENDLINI
DTMINC2II= DTMT

MCALCINCZI) = NINTIDTSTR/DTMTI

DTMTEM = DTSTR/ABSIMCALCINCZIII
DTHI = ABSIDTMTEMI + 0.1E-06
DTLO = ABSIDTMTEMI-0.1506
ABSDTM = ABSIDTMTI
IF IIABSDTM.GT.DTHIl.OR.IABSDTM.LT.DTLOII THEN
STRING = ' The MCD Dt must be the value of storage '
CALL WRTSTRISTRINGI
WRITEISTRING,2100I DTSTR
2100 FORMATI' time I'.
1 1PE12.4,') devided by an interger. Please try again')
CALL WRTSTRISTRINGI
GOTO 123
ENDIF
ENDIF
20 CONTINUE
C
CALL BLNKLN
DO 60 I= 1,3
60 SOLMCDIII = .FALSE.
SOLMCDI1 I = E7YORNI' Do you want to watch results on the '
1 Il'screen7',.FALSE.I
IF ISOLMCDII II THEN
SOLMCDI21= E7YORNI' The state variables?’,.TRUE.I
SOLMCDI3I= E7YORNI' The output variables?’,.TRUE.)
ENDIF
OVFMCD = E7YORNI' Do you want solution range checking?’,.FALSE.l
C
C---- Last chance for user to change mind
CALL CONTUEIPROCFGI

121

IF I.NOT.PROCFGI THEN
GOTO 50
ELSE
CALL BLNKLN
CALL WRTSTRI' integration for M_C_D will commence’)
ENDIF
C
RETURN
END
C>>>>>
C
CMCDINT >>>>>>>>>>>>>>>>>>>>>>>>>> LastChange:01/27/92wa
C
SUBROUTINE MCDINTIIFI'P,TIME,X,P,Y,NUMOT)

C

C--- PURPOSE: Computes Xlt) from TIN to TZSOLV at DTCALC intervals.
C Stores results at DTSTR intervals.

C Uses Runge-Kutta method.

C

Cm INPUTS: IFT P, function type index
TIME, current time

X, inputs

P, parameters

NUMOT, number of outputs
DTM, caculation interval

--- OUTPUTS: Y, outputs of M_C_D

00000000

INCLUDE 'SIZEBK.CBK'
INCLUDE 'SOLNBK.CBK'
INCLUDE 'UTILBK.CBK'
INCLUDE 'MCDCBK.CBK'

INTEGER IFTP, NBUFR, NSX
CHARACTER STRING'BO, 0E8I2OI'72
DOUBLE PRECISION TMCD, TINM
DOUBLE PRECISION X120). P120). YI2OI
DOUBLE PRECISION SXIZOI. 08XI2OI
LOGICAL DOF

EXTERNAL RKMCD, WRTSTR, BLNKLN

EXTERNAL 220801, 220802, 220803, 220804, 220805
EXTERNAL 220806, 220807, 220808, 220809, 220810
EXTERNAL 220811, 220812, 220813, 220814, 220815
EXTERNAL 220816, 220817, 220818, 220819, 220820
EXTERNAL 220821. 220822, 220823, 220824, 220825
EXTERNAL 220826, 220827, 220828, 220829. 220830
EXTERNAL 220831, 220832, 220833, 220834, 220835
EXTERNAL 220836, 220837, 220838, 220839, 220840
EXTERNAL 220841. 220842, 220843, 220844, 220845
EXTERNAL 220846, 220847, 220848, 220849, 220850

122

EXTERNAL 220851 . 220852, 220853. 220854. 220855
EXTERNAL 220856, 220857. 220858. 220859. 220860
EXTERNAL 220861 , 220862, 220863, 220864, 220865
EXTERNAL 220866, 220867, 220868, 220869, 220870
EXTERNAL 220871 , 220872, 220873, 220874, 220875
EXTERNAL 220876, 220877, 220878, 220879, 220880
EXTERNAL 220881 , 220882, 220883, 220884, 220885
EXTERNAL 220886, 220887, 220888, 220889, 220890
EXTERNAL 220891 , 220892. 220893. 220894, 220895
EXTERNAL 220896, 220897, 220898. 20899
C
COOCMCD'NTOOO0.00.0.9...OO0....0.00.90.00.00...00.0.00...0.0.0.0....
C
NC2I = -IIFTP + 99)
IF'IT = -IIFI'P + 99)
TINM = TIME- TIMADV
NSX = NXMCDINCZII
IF ITIME.EO.TINI THEN
DO 2 I = 1 , NSX
SXIII = XOMCDINCZIJI
2 CONTINUE
TINM = TIME
GOTO (8001 ,8002,8003,8004,8005.8006.8007.8008.8009.8010,
8011,8012,8013,8014,8015,8016,8017,8018,8019,8020.
8021 ,8022,8023,8024,8025,8026,8027,8028,8029.8030.
8031 ,8032,8033,8034,8035,8036,8037,8038,8039,8040.
8041 ,8042,8043,8044,8045,8046,8047,8048,8049,8050.
8051 ,8052,8053,8054,8055,8056,8057,8058,8059,8060.
8061 ,8062,8063, 8064,8065,8066, 8067,8068,8069,8070.
8071 .8072,8073,8074,8075,8076,8077,8078,8079,8080.
8081 ,8082,8083,8084,8085.8086.8087.8088,8089,8090.
8091 .8092.8093.8094.8095.8096.8097.8098.8099I.IFTT

(DQNOSU'I-DUN-I

C

8001 CALL 220801ITINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8002 CALL 220802ITINM,X,P,N8X,8X,08X,Y,DOF,NBUFR.0E8I
GOTO 8100

8003 CALL 220803ITINM,X,P,NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8004 CALL 220804ITINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,DES)
GOTO 8100

8005 CALL 2208051TINM,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0E8I
GOTO 8100

8006 CALL 220806ITINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,DE8I
GOTO 8100

8007 CALL 220807ITINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8008 CALL 2208081TINM,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0E8)
GOTO 8100

8009 CALL 220809ITINM,X,P,NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

123

8010 CALL 22081OITINM,X,P,N8X.SX.08X.Y.DOF.NBUFR,0E8)
GOTO 81 00

801 1 CALL 22081 1(TINM.X,P,NSX,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

801 2 CALL 22081 2ITINM,X,P,N8X.SX.08X,Y.DOF.NBUFR,DESI
GOTO 8100

8013 CALL 220813(TINM,X,P.N8X.SX.08X.Y.DOF,NBUFR.DES)
GOTO 81 00

8014 CALL 220814ITINM,X,P,N8X.8X,08X,Y,DOF,NBUFR,0E8)
GOTO 81 00

801 5 CALL 22081 5(TINM,X,P,N8X,8X.08X.Y.DOF.NBUFR.0ESI
GOTO 8100

801 6 CALL 220816ITINM.X.P.NSX,8X,08X.Y.DOF,NBUFR.DESI
GOTO 81 00

801 7 CALL 220817ITINM,X,P,N8X,SX,DSX,Y,DOF.NBUFR,DESI
GOTO 81 00

801 8 CALL 220818ITINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

801 9 CALL 22081 9ITINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8020 CALL 220820ITINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,DE8)
GOTO 8100

8021 CALL 220821 ITINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 81 00

8022 CALL 220822ITINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,DE8I
GOTO 8100

8023 CALL 220823ITINM,X,P,NSX,SX,08X,Y,DOF,NBUFR,0E8I
GOTO 81 00

8024 CALL 220824ITINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 81 00

8025 CALL 220825ITINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 81 00

8026 CALL 220826ITINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8027 CALL 220827ITINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 81 00

8028 CALL 2208281TINM,X,P,NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8029 CALL 220829IT|NM,X,P,N8X,SX,08X,Y,DOF.NBUFR,DE8I
GOTO 8100

8030 CALL 220830ITINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,DESI
GOTO 8100

8031 CALL 220831ITINM,X,P,N8X,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 8100

8032 CALL 220832ITINM,X,P,N8X,SX,08X,Y,DOF.NBUFR,0E8I
GOTO 8100

8033 CALL 220833ITINM,X,P,N8X.SX,08X.Y.DOF.NBUFR,0E8)
GOTO 8100

8034 CALL 220834ITINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,DE8)
GOTO 81 00

8035 CALL 2208351TINM,X,P,NSX,SX.08X,Y,DOF.NBUFR.0E8I

124

GOTO 8100

8036 CALL 220836ITINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8037 CALL 220837ITINM.X,P.NSX,8X,08X.Y.DOF.NBUFR.0ESI
GOTO 8100

8038 CALL 220838ITINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100

8039 CALL 220839ITINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ES)
GOTO 8100

8040 CALL 220840(TINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8041 CALL 220841ITINM,X,P,N8X,SX,DSX,Y,DOF,NBUFR,0ES)
GOTO 8100

8042 CALL 220842ITINM,X,P,NSX,8X,08X.Y,DOF,NBUFR,0ESI
GOTO 8100 -

8043 CALL 220843IT|NM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0ESI
GOTO 8100

8044 CALL 220844ITINM,X,P,N8X,SX,DSX,Y,DOF,NBUFR.0ESI
GOTO 8100

8045 CALL 220845(TINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8046 CALL 2208461TINM,X,P.N8X.SX,08X,Y,DOF,NBUFR,DE8)
GOTO 8100

8047 CALL 220847IT|NM,X,P,NSX,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100

8048 CALL 220848ITINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8049 CALL 220849ITINM.X,P.N8X,8X.08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8050 CALL 2208501TINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,DE8I
GOTO 8100

8051 CALL 220851ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 8100

8052 CALL 2208521TINM,X,P,N8X,SX,DSX,Y,DOF,NBUFR,0ESI
GOTO 8100

8053 CALL 220853ITINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8054 CALL 220854ITINM,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8I
GOTO 8100

8055 CALL 220855ITINM,X,P,NSX,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8056 CALL 220856ITINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8057 CALL 220857ITINM,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8058 CALL 2208581TINM,X,P,NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8059 CALL 220859ITINM,X,P,NSX,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8060 CALL 2208601TINM,X,P,N8X.SX.08X,Y,DOF.NBUFR.0ESI
GOTO 8100

125

8061 CALL ZZDS61ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8062 CALL ZZDS62ITINM,X,P,NSX,SX,DSX,Y.DOF,NBUFR,DES)
GOTO 8100
8063 CALL ZZDS63ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8064 CALL ZZDS64ITINM.X.P.NSX,SX.DSX,Y,DOF,N8UFR,DES)
GOTO 8100
8065 CALL ZZDS65ITINM,X,P,NSX,SX,DSX,Y,DOF,N8UFR,DES)
GOTO 8100
8066 CALL ZZDS66ITINM,X,P.NSX.SX,DSX.Y.DOF.NBUFR,DESI
. GOTO 8100
8067 CALL ZZDS67ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8068 CALL ZZDSGBITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 8100
8069 CALL ZZDS69ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 8100
8070 CALL ZZDS70ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8071 CALL ZZDS71(TINM,X,P,NSX,SX,DSX,Y,DOF.NBUFR,DES)
GOTO 8100
8072 CALL ZZDS72ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8073 CALL ZZDS73ITINM,X,P,NSX,SX,DSX.Y,DOF,NBUFR,DES)
GOTO 8100
8074 CALL ZZDS74ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 8100
8075 CALL ZZDS75ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8076 CALL ZZDS76ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 8100
8077 CALL ZZDS77ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 8100
8078 CALL ZZDS78ITINM,X.P.NSX,SX.DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8079 CALL ZZDS79ITINM,X.P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 8100
8080 CALL ZZDSBOITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8081 CALL ZZDS81lTINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8082 CALL ZZDSBZITINM,X,P,NSX,SX,DSX,Y,DOF,N8UFR,DES)
GOTO 8100
8083 CALL ZZDS83lTINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8084 CALL ZZDS84ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8085 CALL ZZDS85ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 8100
8086 CALL ZZDSB6ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)

126

, GOTO 8100
8087 CALL ZZDSB7ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8088 CALL ZZDS88ITINM.X,P,NSX.SX.DSX.Y,DOF,NBUFR,DES)
GOTO 8100
8089 CALL ZZDS89ITINM.X,P.NSX.SX.DSX.Y,DOF.NBUFR,DES)
GOTO 8100
8090 CALL ZZDSBOITINM,X,P,NSX,SX,DSX.Y,DOF,NBUFR,DES)
GOTO 8100
8091 CALL ZZDS91 IT INM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8092 CALL ZZDS921TINM.X,P,NSX.SX.DSX.Y.DOF,NBUFR.DES)
GOTO 8100
8093 CALL ZZD593ITINM,X,P.NSX,SX,DSX.Y,DOF,N8UFR,DESI
GOTO 8100
8094 CALL ZZD594ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 8100
8095 CALL ZZDS951TINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8096 CALL ZZDS96ITINM,X,P,NSX.SX,DSX,Y,DOF,NBUFR,DESI
GOTO 8100
8097 CALL ZZDS97ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 8100
8098 CALL ZZDS98ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 8100
8099 CALL ZZDS99ITINM,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
C
8100 DO 8 I=1, NSX
XTMCDINCZIJ) = SXII)
XDTMCDINC2I,II = DSXIII
8 CONTINUE
TMCD = TINM
GOTO 999
ELSE
DO 9 I=1, NSX
SXII) = XTMCDINC2l,I)
DSXIII= XDTMCDINCZIJ)
9 CONTINUE
ENDIF
C
C---- This loop for computing SX and DSX
TMCD= TINM
DO 150 NC= 1.MCALCINC2II
TMCD= TMCD-l- DTMINC2II
CALL RKMCDITMCD,IFTP,X,P,NSX,SX,DSX,Y)
150 CONTINUE
C
DO 151 I=1, NSX
XTMCDINC2l,I) = SXIl)
XDTMCDINCZI.II= DSXIII
151 CONTINUE

127

C
C-m-Dump to the screen upon request
999 CONTINUE
IF (SOLMCDI1 I) THEN
CALL BLNKLN
Cm- Write the values
IF (SOLMCDIZII THEN
WRITEISTRING,9810)NC2I,TMCD.IXTMCDINCZIJIJ =1,NSXI
CALL WRTSTRISTRINGI
9810 PORMATI' ZZDS’,I2,':IX)’,6(2X,1PE12.4))
9820 FORMATI' ZZDS',I2,':IY)',612XJPE12.4II
CALL WRTSTRISTRINGI
ENDIF
IF ISOLMCDI3II THEN
WRITEISTRING,9820) NC2I,TMCD.IYIII. l= 1,NUMOT)
CALL WRTSTRISTRINGI
ENDIF
ENDIF
C
C
C-—-- Goodbye and all that
C
C---- Error section
RETURN
END
C>>>>>
C
C
CRKMCD >>>>>>>>>>>>>>>>>>>>>>>>> Last Change: 01/20/92wa
C
SUBROUTINE RKMCDITMCD,IFTP,X,P,NSX,SX,DSX,Y)

C
Cm PURPOSE: Perform integration by Runge-Kutta method.
C
Cm INPUTS: DTM, computation time step
C TMCD, current time ITCALCI
C SX, SX at current time
C
Cm OUTPUTS: SX, at the new time
C IERRF, error flag (=0 if no errors)
C
INCLUDE 'SIZEBK.CBK'
INCLUDE 'SOLNBK.CBK'
INCLUDE 'UTILBK.CBK'
INCLUDE 'MCDCBK.CBK’
C

CHARACTER DESI20l’72

INTEGER N, NBUFR

DOUBLE PRECISION TMCD, TCALC, TIME, TXX

DOUBLE PRECISION K3120).KOI20I.K1I20).K2I20I, DSXI20I
DOUBLE PRECISION SXI20). XI20).PI20),YI20I, XSVI20I

128
LOGICAL DOF

EXTERNAL 220801, 220802, 220803, 220804. 220805
EXTERNAL 220806, 220807, 220808, 220809, 220810
EXTERNAL 220811, 220812, 220813, 220814, 220815
EXTERNAL 220816, 220817. 220818. 220819, 220820
EXTERNAL 220821, 220822, 220823, 220824, 220825
EXTERNAL 220826, 220827. 220828, 220829, 220830
EXTERNAL 220831, 220832, 220833, 220834, 220835
EXTERNAL 220836, 220837, 220838, 220839, 220840
EXTERNAL 220841, 220842, 220843, 220844, 220845
EXTERNAL 220846, 220847, 220848, 220849, 220850
EXTERNAL 220851, 220852, 220853, 220854. 220855
EXTERNAL 220856, 220857, 220858, 220859, 220860
EXTERNAL 220861, 220862, 220863, 220864, 220865
EXTERNAL 220866, 220867, 220868, 220869, 220870
EXTERNAL 220871, 220872, 220873, 220874, 220875
EXTERNAL 220876, 220877, 220878, 220879, 220880
EXTERNAL 220881, 220882, 220883, 220884, 220885
EXTERNAL 220886. 220887, 220888, 220889, 220890
EXTERNAL 220891, 220892, 220893, 220894, 220895
EXTERNAL 220896, 220897, 220898, 220899
C
COOORKMCDOOOO....00....OOOOOOOOOOOOOOOOOQOOOOOOOOOOOOOOOOOOOOOOOOOOO
C
C ----- Save the current 8X values
IFTT= -(IFTP+ 99)
TCALC = TMCD
DO 105 N= 1.N8X
XSVINI = SXINI
105 CONTINUE
C
C ------ Make initial estimate for XI
00110 N=1.NSX
KOIN)= DTMIIF'ITI’ DSXIN)
110 SXINI = XSVINI +0.5‘KOIN)

C ------ Make midpoint correction based on new DXI
.TXX= TCALC +0.5‘DTMIIFTI’)

TIME = TXX

GOTO (9001 ,9002,9003.9004.9005.9006,9007,9008,9009,9010,
9011,9012,9013,9014,9015,9016,9017,9018,9019,9020.
9021 ,9022,9023,9024,9025,9026,9027,9028,9029,9030.
9031 ,9032,9033,9034,9035,9036,9037,9038,9039,9040.
9041 ,9042,9043,9044,9045,9046,9047,9048,9049,9050.
9051 ,9052,9053,9054,9055,9056,9057,9058,9059,9060.
9061 .9062,9063,9064,9065,9066,9067,9068,9069,9070.
9071 ,9072,9073,9074,9075,9076,9077,9078,9079,9080,
9081 .9082,9083,9084,9085,9086,9087.9088,9089.9090.
9091 .9092.9093,9094.9095.9096.9097.9098.9099).IFI'T

CDQNODUI1th-i

129

9001 CALL 220801(TIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100

9002 CALL 220802ITIME,X,P.N8X,8X.DSX.Y.DOF.N8UFR,0E8)
GOTO 9100

9003 CALL 220803lTIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 9100

9004 CALL 220804(TIME,X,P,NSX,SX,08X,Y,DOF,NBUFR.DES)
GOTO 9100

9005 CALL ZZDSOSITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100

9006 CALL ZZDSO6ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100

9007 CALL 220807ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0ESI
GOTO 9100

9008 CALL 220808ITIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DES)
GOTO 9100

9009 CALL 220809ITIME,X,P,NSX,SX.DSX,Y,DOF,NBUFR.0E8)
GOTO 9100

9010 CALL 220810(TlME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DE8I
GOTO 9100

901 1 CALL 2081 1ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100

9012 CALL 22081 2(TIME,X,P,NSX,8X,DSX,Y,DOF,N8UFR,DES)
GOTO 9100

901 3 CALL 2081 3(TIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 9100

9014 CALL 220814ITIME,X,P,NSX.SX,DSX.Y.DOF,NBUFR.DESI
GOTO 9100

9015 CALL 22031 5(TIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR.0E8I
GOTO 9100

9016 CALL 22081 6(TIME,X,P,N8X,8X,DSX,Y,DOF,N8UFR,DES)
GOTO 9100

9017 CALL 220s1 7(TIME,X.P,N8X,SX,DSX,Y,DOF,NBUFR,0ESI
GOTO 9100

901 8 CALL 22081 8(TIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,0ESI
GOTO 9100

901 9 CALL 22081 9(TIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0ESI
GOTO 9100

9020 CALL 220820ITIME,X.P,NSX,SX,08X,Y,DOF,NBUFR,DES)
GOTO 9100

9021 CALL 220821(TIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100

9022 CALL 220822ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100

9023 CALL 220823ITIME,X.P,NSX.SX,DSX.Y.DOF,N8UFR.DES)
GOTO 9100

9024 CALL 220824ITIME,X,P,NSX,SX,08X,Y,DOF,N8UFR,DES)
GOTO 9100

9025 CALL 220825ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0ES)
GOTO 9100

9026 CALL 220826ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)

130

GOTO 9100
9027 CALL 220827ITIME,X,P,NSX,SX,DSX,Y,DOF,N8UFR,DES)
GOTO 9100
9028 CALL 220828ITIME,X,P,N8X,SX,DSX,Y,DOF.NBUFR,DES)
GOTO 9100
9029 CALL 220829ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,DES)
GOTO 9100
9030 CALL ZZDS3OITIME,X,P,NSX,SX,08X,Y,DOF,NBUFR,DES)
GOTO 9100
9031 CALL 220531 (TIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100
9032 CALL ZZDS32ITIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,0ES)
GOTO 9100
9033 CALL 220833lTIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 9100
9034 CALL 220834ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DE8)
GOTO 9100
9035 - CALL 220835ITIME,X,P.N8X,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100
9036 CALL 220836ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100
9037 CALL 220837ITIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,0ES)
GOTO 9100
9038 CALL 220838lTIME,X,P,NSX,SX,08X,Y,DOF,NBUFR,0ES)
GOTO 9100
9039 CALL 220839lTIME.X.P.NSX,8X.DSX.Y.DOF,NBUFR,DES)
GOTO 9100
9040 CALL 220840ITIME,X,P,N8X.SX.08X.Y.DOF,NBUFR,DES)
, GOTO 9100
9041 CALL 220541 (TIME,X,P.NSX.SX.DSX.Y,DOF,NBUFR,DE8)
GOTO 9100
9042 CALL 220842(TIME.X.P.NSX.SX.DSX,Y,DOF,NBUFR,DE8)
GOTO 9100
9043 CALL 220843ITIME.X,P,NSX,SX.DSX,Y.DOF,NBUFR,DES)
GOTO 9100
9044 CALL 220844ITIME.X.P,NSX,8X.DSX.Y,DOF,NBUFR,0ES)
GOTO 9100
9045 CALL 220845ITIME,X.P.NSX.SX.DSX.Y,DOF.NBUFR,DE8)
GOTO 9100
9046 CALL 220846ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DES)
GOTO 9100
9047 CALL 220847ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100
9048 CALL 220848ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100
9049 CALL 220849ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0ES)
GOTO 9100
9050 CALL ZZDSSOITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100
9051 CALL 220551 (TIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 9100

131

9052 CALL 220852ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8I
GOTO 9100

9053 CALL 220853ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 9100

9054 CALL 220854ITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 91 00

9055 CALL 2208551TIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 9100

9056 CALL 220856ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DES)
GOTO 9100

9057 CALL 220857ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DE8)
GOTO 9100

9058 CALL 220858ITIME,X,P,N8X,8X,DSX.Y.DOF,NBUFR,DE8I
GOTO 9100

9059 CALL 220859ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DESI
GOTO 9100

9060 CALL 220860ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 91 00

9061 CALL 220861(TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 9100

9062 CALL 220862ITIME.X,P,N8X.SX.08X,Y.DOF,NBUFR,0E8)
GOTO 9100

9063 CALL 220863ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 91 00

9064 CALL 220864ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DESI
GOTO 9100

9065 CALL 220865ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 9100

9066 CALL 220866ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DE8I
GOTO 9100

9067 CALL 220867ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8I
GOTO 91 00

9068 CALL 220868ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 9100

9069 CALL 220869ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DESI
GOTO 9100

9070 CALL 220870ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 9100

9071 CALL 220871(TIME,X,P,NSX,SX,08X,Y,DOF,NBUFR.0E8)
GOTO 9100

9072 CALL 220872ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8I
GOTO 9100

9073 CALL 220873ITIME.X.P.N8X.8X.08X.Y.DOF,NBUFR.DESI
GOTO 91 00

9074 CALL 220874ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 9100

9075 CALL 2208751TIME,X,P,NSX,SX,08X,Y,DOF,NBUFR,0ES)
GOTO 9100

9076 CALL 220876ITIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,0E8)
GOTO 91 00

9077 CALL 220877ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)

132

GOTO 9100

9078 CALL 220878ITIME,X,P,NSX,SX,DSX.Y.DOF,NBUFR.0ESI
GOTO 9100

9079 CALL 220879ITIME,X,P,NSX.SX,DSX,Y.DOF,N8UFR.DES)
GOTO 9100

9080 CALL ZZDSBOITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100

9081 CALL 220881(TIME,X,P,NSX,SX,DSX.Y,DOF,NBUFR,DESI
GOTO 9100

9082 CALL 220882lTIME,X,P,NSX,SX,DSX,Y,DOF,NBUPR,DE8)
GOTO 9100

9083 CALL 220883ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 9100

9084 CALL 220884lTIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100

9085 CALL 220885(TIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,DES)
GOTO 9100

9086 CALL 220886ITIME.X,P,NSX.SX.DSX,Y.DOF.NBUFR,DES)
GOTO 9100

9087 CALL 220887ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0E8)
GOTO 9100

9088 CALL 220888ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 9100

9089 CALL 220889ITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,DE8I
GOTO 9100

9090 CALL 220890ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 9100

9091 CALL ZZDS91(TIME,X,P,NSX,SX,08X,Y,DOF,NBUFR,DE8)
GOTO 9100

9092 CALL 220892ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 9100

9093 CALL 220893ITIME,X.P.N8X.SX.08X.Y.DOF,NBUFR,DE8)
GOTO 9100

9094 CALL 220894ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0E8I
GOTO 9100

9095 CALL 220895ITIME,X,P,NSX,SX,08X,Y,DOF,NBUFR,DES)
GOTO 9100

9096 CALL 220896ITIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 9100

9097 CALL ZZDS97ITIME,X,P,N8X,SX,DSX,Y,DOF,N8UFR,DES)
GOTO 9100

9098 CALL 2208981TIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,0E8I
GOTO 9100

9099 CALL 220899ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI

C

9100 00120 N=1,NSX

K1INI= DTMIIFI'T)‘ DSXINI

120 SXINI= XSVINI+ 0.5'K1INI

C

Cmm Make final estimate based on new DSX
TXX= TCALC +0.5’DTMIIFIT I

133

TIME = TXX
GOTO (7001 .7002.7003.7004.7005.7006.7007,7008.7009,7010,

1 7011,7012,7013,7014,7015,7016,7017,7018,7019,7020.
2 7021 ,7022.7023,7024.7025.7026,7027.7028,7029.7030.
3 7031 ,7032,7033,7034,7035,7036,7037,7038,7039,7040.
4 7041 ,7042,7043,7044,7045,7046,7047,7048,7049,7050.
5 7051 ,7052,7053,7054,7055,7056.7057,7058,7059,7060.
6 7061 ,7062,7063,7064,7065,7066,7067,7068,7069,7070,
7 7071 ,7072,7073,7074,7075,7076,7077,7078,7079,7080.
8 7081 ,7082,7083,7084.7085,7086,7087.7088,7089,7090.
9 7091 .7092,7093.7094,7095.7096.7097.7098.7099I,IFTT
C
7001 CALL 220801ITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,0ESI
GOTO 71 00
7002 CALL 220802ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0ESI
GOTO 71 00
7003 CALL 220803IT|ME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DESI
GOTO 71 00
7004 CALL 220804ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 71 00
7005 CALL ZZDSOSITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 71 OO '
7006 CALL 220806ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 71 00
7007 CALL 220807ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 71 00
7008 CALL 220808ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 71 00
7009 CALL ZZDSO9IT|ME,X,P,NSX,8X,08X,Y,0OF,NBUFR,DESI
GOTO 7100
7010 CALL 22081OITIME,X,P,N8X.8X,08X,Y,DOF,NBUFR,0E8)
GOTO 71 00
701 1 CALL 220811(TIME,X,P,N8X,8X,08X,Y,DOF.NBUFR,DESI
GOTO 71 00
7012 CALL 2081 2(TIME.X,P,NSX.8X.08X.Y.DOF,NBUFR,DE8)
GOTO 7100 A
7013 CALL 220813(TIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0E8)
GOTO 7100
7014 CALL 220814(TIME,X,P.NSX.SX.08X.Y,OOF.NBUFR.DESI
GOTO 71 00
7015 CALL 22081 5ITIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,0E8)
GOTO 7100
7016 CALL 220816ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0ESI
GOTO 71 00
7017 CALL 22081 7(TIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 71 00
701 8 CALL 22081 8(TIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0E8)
GOTO 71 OO ’
7019 CALL 220819(TIME.X.P,N8X.8X.08X.Y.DOF,NBUFR,0E8)
GOTO 71 00

7020 CALL ZZDSZOITIME.X,P,N8X,8X,08X,Y,DOF,NBUFR,DE8)

134

GOTO 7100
7021 CALL 220821(TIME,X,P,N8X,8X,08X.Y.DOF,NBUFR,DESI
GOTO 7100
7022 CALL 2208221TIME,X,P,N8X.SX,08X,Y,DOF,NBUFR,DES)
GOTO 7100
7023 CALL 220823ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 7100
7024 CALL 220824ITIME,X,P,NSX,SX,08X,Y,DOF.NBUFR,0ESI
GOTO 7100 .
7025 CALL 220825(TIME,X,P,NSX,SX,DSX.Y.DOF,NBUFR,0E8)
GOTO 7100
7026 CALL 220826ITIME,X,P,N8X,8X,08X,Y.DOF,NBUFR,DESI
GOTO 7100
7027 CALL 220827ITIME.X.P.N8X.8X.08X.Y.DOF.NBUFR.0ES)
GOTO 7100
7028 CALL 220828(TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 7100
7029 CALL 220829lTIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,0E8)
GOTO 7100
7030 CALL 22083OITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DESI
GOTO 7100
7031 CALL 220831(TIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DESI
GOTO 7100
7032 CALL 220832ITIME.X,P,NSX,8X,08X,Y,DOF,NBUFR,DES)
GOTO 7100
7033 CALL 220833ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0E8I
GOTO 7100
7034 CALL 220834ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 7100
7035 CALL 2208351TIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 7100
7036 CALL 220836ITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,0E8)
- GOTO 7100
7037 CALL 220837ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DE8I
GOTO 7100
7038 CALL 220838ITIME,X,P,N8X,8X,08X.Y,DOF,NBUFR,0E8)
GOTO 7100
7039 CALL 220839ITIME,X,P,N8X.8X,DSX,Y.OOF,NBUFR,DESI
GOTO 7100
7040 CALL 22084OITIME,X,P.NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 7100
7041 CALL 220841(TIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 7100
7042 CALL 220842ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0E8)
GOTO 7100
7043 CALL 220843ITIME.X.P.N8X.8X.08X.Y,DOF.NBUFR.DESI
GOTO 7100
7044 CALL 220844ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 7100
7045 CALL 2208451TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 7100

135

7046 CALL 220846ITIME,X,P,NSX,SX,DSX.Y.DOF,NBUFR,DES)
GOTO 7100

7047 CALL 220847ITIME.X.P,NSX.SX,08X.Y.DOF,NBUFR,DES)
GOTO 7100

7048 CALL 220848lTIME,X,P,NSX,SX.DSX,Y,DOF.NBUFR,DES)
GOTO 7100

7049 CALL 220849ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0E8)
GOTO 7100

7050 CALL 220850lTlME,X,P,NSX,SX,DSX.Y.DOF,NBUFR,DESI
GOTO 7100

7051 CALL 220851lTlME,X,P,NSX.SX,DSX,Y,DOF,NBUFR,DESI
GOTO 7100

7052 CALL ZZDS52ITIME,X,P.NSX.SX,DSX,Y.DOF.NBUFR,DES)
GOTO 7100

7053 CALL 220853ITIME,X,P,NSX.SX,DSX,Y,DOF,N8UFR,DESI
GOTO 7100

7054 CALL 220854ITIME, x P, NSX, 5x sz v DOF, NBUER, DESI
GOTO 7100

7055 CALL 220855lTIME, x, P, NSX, 5x, sz, v DOF, NBUER, DESI
GOTO 7100

7056 CALL 220856ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 7100

7057 CALL 220857ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 7100

7058 CALL 220858ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0ES)
GOTO 7100

7059 CALL 220859ITIME,X,P,NSX.8X,DSX,Y.DOF.NBUFR.DES)
GOTO 7100

7060 CALL 220860ITIME,X,P,N8X,8X,DSX,Y,DOF,N8UFR,DES)
GOTO 7100

7061 CALL 220561 (TIME.X.P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 7100

7062 CALL 220862ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DES)
GOTO 7100

7063 CALL 220863ITIME,X,P,NSX.SX,DSX,Y,DOF.NBUFR,DESI
GOTO 7100

7064 CALL 220864lTlME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,0ESI
GOTO 7100 ,

7065 CALL 220865(TIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 7100

7066 CALL 220866(TIME,X,P,N8X,SX,08X,Y,DOF,N8UFR,DES)
GOTO 7100

7067 CALL ZZDS67ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 7100

7068 CALL 220868ITIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DES)
GOTO 7100

7069 CALL 220869ITIME,X,P,N8X,SX,08X,Y,DOF,N8UFR,DESI
GOTO 7100

7070 CALL 220870lTIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,0E8)
GOTO 7100

7071 CALL 220571 (TIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0ES)

136

GOTO 7100

7072 CALL 220872ITIME.X.P.N8X,8X.08X.Y,DOF,NBUFR,DESI
GOTO 7100

7073 CALL 220873ITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR.0ESI
GOTO 7100

7074 CALL 220874ITIME,X.P,N8X.8X.08X.Y.DOF.NBUFR.OESI
GOTO 7100

7075 CALL 220875ITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 7100

7076 CALL 2208761TIME,X,P,N8X,SX,08X.Y.DOF,NBUFR.0E8I
GOTO 7100

7077 CALL 220877ITIME,X,P,NSX,8X,08X.Y,DOF,NBUFR,DESI
GOTO 7100

7078 CALL 220878(TIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DE8)
GOTO 7100 _

7079 CALL 220879ITIME,X.P.N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 7100 .

7080 CALL ZZDSBOITIME,X,P,N8X.8X,08X,Y,DOF,NBUFR,0E8)
GOTO 7100

7081 CALL 220881(TIME,X,P,N8X.8X.08X,Y.DOF,NBUFR,DESI
GOTO 7100

7082 CALL 2208821TIME,X,P,NSX,SX,08X,Y,DOF,NBUFR,DESI
GOTO 7100

7083 CALL 220883ITIME,X.P,N8X,8X,08X,Y,DOF,NBUFR,DES)
GOTO 7100

7084 CALL 2208841TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR.0E8I
GOTO 7100

7085 CALL 2208851T IME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 7100

7086 CALL 220886ITIME,X,P,N8X,8X,DSX,Y.DOF.NBUFR,DESI
GOTO 7100

7087 CALL 220887ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DESI
GOTO 7100

7088 CALL 220888ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DESI
GOTO 7100

7089 CALL 220889ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 7100

7090 CALL 22087OITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,DE8I
GOTO 7100

7091 CALL 220891(TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DES)
GOTO 7100

7092 CALL 2208921TIME.X.P.N8X.SX,08X.Y.DOF,NBUFR.0ESI
GOTO 7100

7093 CALL 220893ITIME,X,P,N8X,SX,08X,Y.DOF,NBUFR,DESI
GOTO 7100

7094 CALL 220894ITIME,X.P.N8X.8X.08X.Y.DOF.NBUFR,0E8)
GOTO 7100

7095 CALL 220895ITIME,X,P,NSX,SX.DSX,Y,DOF,NBUFR,DESI
GOTO 7100

7096 CALL 2208961TIME.X.P.N8X.SX.08X.Y.DOF.NBUFR.DESI
GOTO 7100

137

7097 CALL 220897ITIME,X,P,N8X.8X,08X.Y,DOF,NBUFR,0E8)
GOTO 7100
7098 CALL 2208981TIME.X.P.NSX.8X.08X.Y,DOF,NBUFR,0E8)
GOTO 7100
7099 CALL 220899ITIME,X.P,N8X,SX,08X,Y,DOF,NBUFR,0E8I
7100 DO 130 N=1,N8X
K2INI= DTMIIFTTI" DSXINI
130 SXINI = XSVINI +K2INI
TXX= TCALC +0TMI|FTTI
TIME = TXX
GOTO (8001 ,8002,8003.8004.8005,8006,8007,8008,8009,8010.

1 8011,8012,8013,8014,8015,8016,8017,8018,8019,8020.
2 8021 ,8022,8023,8024.8025,8026,8027.8028,8029,8030.
3 8031 ,8032.8033,8034,8035.8036.8037,8038.8039,8040.
4 8041 ,8042,8043,8044,8045,8046,8047,8048,8049,8050.
5 8051 ,8052,8053,8054,8055,8056,8057,8058,8059,8060.
6 8061 ,8062.8063.8064.8065,8066,8067,8068,8069,8070.
7 8071 ,8072,8073,8074,8075,8076,8077,8078,8079,8080.
8 8081 .8082,8083,8084,8085, 8086, 8087,8088,8089,8090.
9 8091 .8092.8093,8094.8095,8096,8097,8098,8099I.IFI'T
C
8001 CALL 220801ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 8100
8002 CALL 2208021TIME,X,P,N8X.8X,08X.Y.DOF,N8UFR.DES)
GOTO 8100
8003 CALL 220803ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 81 00
8004 CALL 220804(TIME.X,P,N8X.8X,08X,Y,DOF,NBUFR,0E8I
GOTO 81 00
8005 CALL 220805ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DE8)
GOTO 8100
8006 CALL 2208061TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DE8)
GOTO 8100
8007 CALL 220807ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100
8008 CALL 220808ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100
8009 CALL 220809ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,DE8I
GOTO 8100
8010 CALL 22081OITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100
801 1 CALL 220811(TIME,X,P,N8X,SX,08X,Y,DOF,NBUFR.0ESI
GOTO 8100 -
801 2 CALL 220812(TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100
801 3 CALL 220813(TIME,X.P.N8X.8X.DSX.Y.DOF,NBUFR.DESI
GOTO 8100
8014 CALL 220814(TIME,X,P,N8X,SX,08X,Y,DOF,N8UFR,DESI
GOTO 8100

8015 CALL 220815ITIME.X.P.N8X.8X.08X.Y.DOF,NBUFR,DESI
GOTO 8100

138

801 6 CALL 22081 6ITIME,X,P,N8X,8X,08X.Y.DOF,NBUFR,DESI
GOTO 8100

8017 CALL 22081 7ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

801 8 CALL 22081 8(TIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8019 CALL 22081 9ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100

8020 CALL 2208201TIME,X,P,NSX,8X,08X.Y,DOF.NBUFR.0ESI
GOTO 8100

8021 CALL 220821(TIME,X,P,N8X,SX,08X,Y,DOF.NBUFR.0E8I
GOTO 8100

8022 CALL 2208221TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100

8023 CALL 220823ITIME,X,P.N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100

8024 CALL 220824ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8025 CALL 2208251TIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0E8I
GOTO 8100

8026 CALL 2208261TIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8027 CALL 220827ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DE8)
GOTO 8100

8028 CALL 220828(TIME,X,P,N8X,SX,08X,Y,DOF,NBUFR.0E8I
GOTO 8100

8029 CALL 220829ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100

8030 CALL 22083OITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8I
GOTO 8100 -

8031 CALL 220831ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0E8)
GOTO 8100

8032 CALL 2208321TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8033 CALL 220833ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0ES)
GOTO 8100

8034 CALL 220834ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8035 CALL 220835ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100

8036 CALL 2208361T|ME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8037 CALL 220837ITIME,X.P,NSX,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100

8038 CALL 2208381TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8039 CALL 220839ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DE8)
GOTO 8100

8040 CALL 22084OITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100

8041 CALL 220841ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI

139

GOTO 8100

8042 CALL 220842ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8043 CALL 220843ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DE8)
GOTO 8100

8044 CALL 220844ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8045 CALL 220845ITIME.X.P.NSX.SX.08X.Y.DOF.NBUFR.DESI
GOTO 8100

8046 CALL 2208461TIME,X,P,NSX,8X,DSX.Y,DOF,NBUFR,DESI
GOTO 8100

8047 CALL 220847ITIME.X.P.NSX.SX.08X,Y,DOF,NBUFR.DESI
GOTO 8100

8048 CALL ZZDS48ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100 .

8049 CALL 220849(TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8050 CALL 220850ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8051 CALL 220851ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100

8052 CALL 220852ITIME,X,P,N8X,SX,08X,Y,0OF,NBUFR,0ES)
GOTO 8100

8053 CALL 220853ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DESI
GOTO 81 00

8054 CALL 220854(TIME,X,P,NSX,SX,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8055 CALL 220855IT|ME,X,P,NSX,SX,08X,Y,DOF,NBUFR,DESI
GOTO 8100

8056 CALL 220856ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8I
GOTO 8100

8057 CALL 220857ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DES)
GOTO 8100 '

8058 CALL 2208581TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DE8)
GOTO 8100

8059 CALL 220859IT|ME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ES)
GOTO 8100

8060 CALL 22086OITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8061 CALL 220861(TIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DE8)
GOTO 8100

8062 CALL 2208621T|ME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 8100

8063 CALL 220863ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100

8064 CALL ZZDS64ITIME,X,P,N8X,SX,DSX.Y.DOF.NBUFR,DESI
GOTO 8100

8065 CALL 220865ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100

8066 CALL 2208661TIME,X,P.N8X.SX.DSX.Y.DOF,NBUFR,DESI
GOTO 8100

140

8067 CALL 220867ITIME.X.P,N8X.SX.08X.Y.DOF,NBUFR,0ES)
GOTO 8100
8068 CALL 2208681TIME.X,P,NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100
8069 CALL 220869ITIME,X,P,N8X,SX,08X.Y,DOF,NBUFR,0ESI
GOTO 8100
8070 CALL 220870ITIME.X,P.N8X.8X.08X.Y,DOF,NBUFR,0E8)
-GOTO 8100
8071 CALL 220871(TlME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100
8072 CALL 2208721TIME.X.P.N8X.SX.08X.Y.DOF.NBUFR.0ESI
GOTO 8100
8073 CALL 220873ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0ESI
GOTO 8100
8074 CALL 220874(T|ME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
- GOTO 8100
8075 CALL 2208751TIME.X.P.N8X,8X,08X.Y,DOF,NBUFR,0ESI
GOTO 8100
8076 CALL 220876ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 8100
8077 CALL 220877ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100
8078 CALL 220878ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100
8079 CALL 2208791TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 81 00
8080 CALL 2208801TIME.X.P.N8X.8X,08X.Y.DOF,NBUFR.0ESI
GOTO 8100
8081 CALL 220881ITIME,X,P.N8X,SX.08X.Y.DOF,NBUFR,0E8)
GOTO 8100
8082 CALL 220882ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DE8)
GOTO 8100
8083 CALL 220883ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 8100
8084 CALL 220884ITIME.X.P.NSX.SX.08X.Y,DOF.NBUFR,0ESI
GOTO 8100
8085 CALL 2208851T|ME,X.P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 81 00
8086 CALL 2208861TIME.X,P,N8X,8X,08X,Y,DOF,NBUFR,0ES)
GOTO 8100
8087 CALL 220887ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 8100
8088 CALL 220888(TIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100
8089 CALL 2208891TIME,X.P.N8X,8X,08X,Y.DOF,NBUFR,0ES)
GOTO 8100
8090 CALL 220880ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 8100
8091 CALL 220891(TIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0E8I
GOTO 8100
8092 CALL 220892ITIME,X,P,N8X,SX,DSX.Y,DOF.NBUFR,0ESI

141

GOTO 8100

8093 CALL 220893ITIME,X,P,NSX,8X.DSX.Y,DOF,NBUFR,DES)

' GOTO 8100

8094 CALL ZZDS94ITIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DES)
GOTO 8100

8095 CALL ZZDS95ITIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,0ES)
GOTO 8100

8096 CALL 220896ITIME,X,P,N8X,SX,DSX.Y.DOF,N8UFR,DES)
GOTO 8100

8097 CALL 220897ITIME,X,P.NSX,8X.08X.Y,DOF,NBUFR,DES)
GOTO 8100

8098 CALL 220898ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0E8)
GOTO 8100

8099 CALL 220899lTIME.X.P.N8X.SX.DSX.Y.DOF.NBUFR,DES)

C . I

8100 00140 N=1,N8X

K3INI= DTMIIFIT)‘ DSXINI

140 SXIN) = XSVINI + (KOIN) + 2.0°(K1(N)+ K2lN)) + K3IN))/6.0

C

C ------ Update the 08X vector

GOTO (6001 .6002,6003.6004.6005.6006.6007,6008.6009,6010.
6011,6012,6013,6014,6015,6016,6017,6018,6019,6020.
6021 ,6022,6023,6024,6025,6026,6027,6028,6029,6030.

' 6031,6032,6033,6034,6035,6036,6037,6038,6039,6040.
6041 ,6042,6043,6044,6045,6046,6047,6048,6049,6050,
6051 ,6052,6053,6054,6055,6056,6057,6058,6059,6060.
6061 ,6062,6063,6064,6065,6066,6067,6068,6069,6070.
6071 ,6072,6073,6074,6075,6076,6077,6078,6079,6080.
6081 ,6082,6083,6084,6085,6086,6087,6088,6089,6090.
6091 ,6092,6093,6094,6095,6096,6097.6098,6099),IFIT

CDQNOIUIhUN-fi

C

6001 CALL 220801(TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8I
GOTO 6100

6002 CALL 2208021TIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 6100

6003 CALL 220803ITIME,X,P,NSX,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 6100 ‘

6004 CALL 220804(TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 6100 '

6005 CALL ZZDSOSITIME,X,P,N8X,SX,08X.Y,DOF,NBUFR,0E8I
GOTO 6100

6006 CALL 2208061TIME.X,P,N8X,8X,08X,Y,DOF,NBUFR,DES)
GOTO 6100

6007 CALL 220807ITIME,X,P,NSX,SX,08X,Y,DOF,NBUFR,0E8I
GOTO 6100

6008 CALL 220808ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8)
GOTO 6100

6009 CALL 220809ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0E8)
GOTO 6100

6010 CALL 22081OITIME,X,P,N8X,SX.DSX.Y,DOF.NBUFR,DE8)
GOTO 6100

142

6011 CALL 22081 1ITIME,X,P,N8X,8X.08X,Y.DOF.NBUFR,0E8)
GOTO 6100

601 2 CALL 22081 2(TIME,X,P,N8X,SX,08X,Y.DOF,NBUFR,DESI
GOTO 6100

601 3 CALL 220813ITIME,X, P,N8X,SX,08X,Y,DOF,NBUFR,DESI
GOTO 6100

601 4 CALL ZZDS1 4ITIME,X,P,N8X,SX,08X,Y.DOF.NBUFR,DESI
GOTO 6100

6015 CALL 22081 5ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0E8)
GOTO 6100

601 6 CALL 220816ITIME.X,P,NSX,SX,08X,Y,DOF,NBUFR,0E8)
GOTO 6100

601 7 CALL 22081 7ITIME.X,P.N8X.SX.DSX.Y.DOF,NBUFR.0ES)
GOTO 6100

6018 CALL 22081 8(TIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DESI
GOTO 6100

601 9 CALL 22081 9ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DESI
GOTO 6100 '

6020 CALL ZZDSZOITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DESI
GOTO 6100

6021 CALL 220821 (TIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0ESI
GOTO 61 00

6022 CALL 220822ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 6100 -

6023 CALL 220823ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0ES)
GOTO 6100

6024 CALL 220824ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ES)

' GOTO 6100

6025 CALL 220825ITIME.X,P,N8X,SX.08X.Y.DOF,NBUFR,0E8)
GOTO 6100

6026 CALL 2208261TIME,X,P,N8X,8X,08X,Y.DOF,NBUFR,0ESI
GOTO 61 00

6027 CALL 220827ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 6100

6028 CALL 220828IT|ME,X,P,N8X,SX.08X,Y,DOF,NBUFR,DESI
GOTO 61 00

6029 CALL 220829ITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 6100

6030 CALL 220830ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DESI
GOTO 6100

6031 CALL 220831 (TIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,DES)
GOTO 61 OO -

6032 CALL 220832ITIME,X,P,NSX.SX,08X,Y,DOF,NBUFR,0E8)
GOTO 6100

6033 CALL 220833ITIME,X,P,NSX.SX,08X,Y,DOF.NBUFR.DESI
GOTO 6100

6034 CALL 220834ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DESI
GOTO 6100

6035 CALL 220835ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0E8I
GOTO 6100

6036 CALL 220836ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0E8)

143

GOTO 6100

6037 CALL 220837ITIME,X,P,NSX.SX.DSX.Y.DOF,NBUFR.DESI
GOTO 6100

6038 CALL 220838ITIME.X.P,NSX,SX,DSX,Y,DOF,NBUFR,DE8)
GOTO 6100

6039 CALL 220839ITIME,X,P.NSX,SX,DSX,Y,DOF,NBUFR,0ESI
GOTO 6100

6040 CALL 220840(TIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DES)
GOTO 6100

6041 CALL 220841ITIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DES)
GOTO 6100

6042 CALL 220842ITIME,X,P,NSX,SX,DSX.Y.DOF,NBUFR,DES)
GOTO 6100

6043 CALL 220843lTIME,X.P.NSX.SX.DSX.Y.DOF,NBUFR,DESI
GOTO 6100

6044 CALL ZZDS44ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 6100

6045 CALL 220845ITIME,X.P,NSX,SX.08X,Y.DOF,NBUFR,DES)
GOTO 6100

6046 CALL 220846ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR.DES)
GOTO 6100 ,

6047 CALL 220847ITIME,X,P,NSX,SX,DSX,Y,DOF.NBUFR.DES)
GOTO 6100

6048 CALL ZZDS48ITIME,X,P,NSX,SX,08X,Y,DOF,NBUFR,DE8)
GOTO 6100

6049 CALL 220849ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0E8)
GOTO 6100

6050 CALL 220850lTIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 6100 A

6051 CALL 220851ITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 6100

6052 CALL 220852ITIME.X.P,NSX,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 6100

6053 CALL 220853lTIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 6100 ‘

6054 CALL 220854ITIME,X,P,N8X,SX,08X,Y,DOF.NBUFR.DES)
GOTO 6100

6055 CALL 220855(TIME,X,P,NSX,SX,08X,Y,DOF,NBUFR,DES)
GOTO 6100

6056 CALL 220856ITIME,X,P,N8X,SX,08X,Y,DOF,NBUFR,DES)
GOTO 6100

6057 CALL 220857ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DES)
GOTO 6100

6058 CALL 220858(TIME,X,P,NSX.SX,DSX,Y,DOF,NBUFR,DESI
GOTO 6100

6059 CALL 220859ITIME,X,P,N8X,SX,DSX.Y,DOF,NBUFR,DES)
GOTO 6100

6060 CALL 220860lTIME,X,P,NSX,SX,08X,Y,DOF,NBUFR,DES)
GOTO 6100

6061 CALL 220861(TIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DES)
GOTO 6100

144

6062 CALL 2208621TIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,0E8)
GOTO 6100

6063 CALL 220863ITIME.X,P.NSX,8X.08X,Y,DOF,NBUFR,DESI
GOTO 6100

6064 CALL 220864ITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,0ESI
GOTO 6100

6065 CALL 2208651TIME.X,P,N8X.SX.08X.Y.DOF.NBUFR.0ESI

‘ GOTO 6100

6066 CALL ZZDS66ITIME,X,P,N8X.SX.08X,Y,DOF,NBUFR,0E8)
GOTO 6100

6067 CALL 220867ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 6100

6068 CALL 220868ITIME,X,P,N8X,8X,08X.Y.DOF,NBUFR,DE8I
GOTO 6100

6069 CALL 220869ITIME,X,P,NSX,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 6100

6070 CALL 22087OITIME,X,P,N8X,SX,08X,Y,DOF.NBUFR.0ESI
GOTO 6100

6071 CALL 220871(TIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 6100

6072 CALL 220872ITIME,X,P,N8X,8X,08X.Y.DOF,NBUFR,0ESI
GOTO 6100 .

6073 CALL 220873ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 6100

6074 CALL 220874ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DE8)
GOTO 6100 A

6075 CALL 220875ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR.DESI
GOTO 6100

6076 CALL 220876ITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,DE8I
GOTO 6100

6077 CALL 220877ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,DESI
GOTO 6100

6078 CALL 220878ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 6100

6079 CALL 220879ITIME.X,P,N8X,8X,DSX,Y,DOF,NBUFR,DESI
GOTO 6100 ‘

6080 CALL 22088OIT|ME,X,P,N8X,SX,08X,Y,DOF,NBUFR,0ESI
GOTO 61 00

6081 CALL 220881(TIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0E8)
GOTO 6100

6082 CALL 220882ITIME,X.P,NSX,SX,DSX,Y,DOF,NBUFR,0E8)
GOTO 6100

6083 CALL 220883ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0E8)
GOTO 6100

6084 CALL 220884ITIME,X,P,NSX,8X,08X,Y,DOF,NBUFR,0ESI
GOTO 6100

6085 CALL 2208851TIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,0E8)
GOTO 6100

6086 CALL 220886ITIME,X,P,N8X,8X,08X,Y,DOF,NBUFR,0E8I
GOTO 6100

6087 CALL 220887ITIME,X,P,N8X,8X,DSX,Y,DOF,NBUFR,0ESI

145

GOTO 6100

6088 CALL 220888ITIME.X.P,NSX.SX,DSX,Y,DOF.N8UFR.DE8)
GOTO 6100

6089 CALL 220889ITIME,X,P,N8X,SX,DSX,Y,DOF.NBUFR,DESI
GOTO 6100

6090 CALL 220890ITIME,X,P,N8X.SX,DSX.Y,DOF.NBUFR,DESI
GOTO 6100

6091 CALL ZZDS91(TIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0ESI
GOTO 6100

6092 CALL ZZDS92(TIME,X,P,NSX,8X,DSX,Y,DOF,N8UFR,DESI
GOTO 6100 ‘

6093 CALL 220893ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DE8)
GOTO 6100

6094 CALL 220894ITIME,X.P,NSX.SX,08X.Y.DOF.NBUFR,DES)
GOTO 6100

6095 CALL 220895lTIME.X,P,N8X,SX,DSX,Y,DOF,NBUFR,DES)
GOTO 6100

6096 CALL 220896ITIME.X,P,NSX,SX.08X,Y,DOF.NBUFR,DES)
GOTO 6100

6097 CALL 220897ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,DESI
GOTO 6100

6098 CALL 220898ITIME,X,P,N8X,SX,DSX,Y,DOF,NBUFR,DES)

' GOTO 6100

6099 CALL 220899ITIME,X,P,NSX,SX,DSX,Y,DOF,NBUFR,0ES)

C

C ------ Every thing is alright here

6100 IERRF= 0

RETURN
C

END
C>>>>>
C
C

CINIMCD>>>>>>>>>>>>>>>>>>>>>>>> Last Change: 01/08/92 wa
C
SUBROUTINE INIMCDINCZII

C
Cm PURPOSE: Get the initial conditions from user for a typical
C MCD node.
C
Cm INPUTS: NC2I, index of 2208i]
C
Cm OUTPUTS: XOMCDINCZIJI, initial conditions
C
INCLUDE 'SIZEBK.CBK'
INCLUDE 'SOLNBK.CBK'
INCLUDE 'MCDCBK.CBK'
C

REAL RLO, RHI, RVAL
CHARACTER STRING'72
INTEGER N. NC2|

146

LOGICAL NEWLIN.ENDLIN
C
EXTERNAL BLNKLN, WRTSTR, VI2C8, NCHARS, PROMPT, GET RL
g0.0'NIMCDOOCO...O...OO‘OOOOOOOOOOOOOOOOOOOOOOOOOO0.00.00.99.00.0.90....
C
RLO =.-1.E25
RHI=1.E25
CALL BLNKLN
STRING = ' Enter the initial conditions:'
CALL WRTSTRISTRINGI
1200 FORMATI' XI',I2,'I 7',T15,'(',1PE12.4,'I:'I
C
00 10 I=1, NXMCDINCZI)
WRITEISTRING,1200I I,XOMCD(NC2I,I)
CALL PROMPTISTRINGI
RVAL=XOMCDINC2I,II
NEWLIN = .TRUE.
CALL GETRLIRVAL,RLO,RHI,NEWLIN,ENDLIN)
XOMCDINCZI,I)= RVAL
10 CONTINUE
C
RETURN
END
C>>>>>
C
C

APPENDIX E. ALGORITHIVI FOR OBTAINING
THE PATH-ORDER MATRIX

Path-Order Matrix
The path-order matrix is built only when all the paths have been identified. From a

graphical point of view, the problem can be stated as follows:

Case 1: SCG does not contain acyclic sub—digraph.

Definition:
If two directed paths have the same start and end vertices, this digraph is referred

to as acyclic.

An example of an acyclic digraph is given in Figure E.1

Here, one path is Vl 3% V3 3 V, and another is VI 53 V4 93 V5. Both paths share the

same start and end vertices.

To find path matrix, select a vertex and put it on an initially empty queue of vertices to

be visited. We repeatedly remove the vertex t at the head of the queue. Check incident

147

148
edges and then place onto the queue all the vertices adjacent to t.

 

 

 

 

 

Figure E.1 An acyclic digraph

Case 2: SCG contains acych sub-digraphs.

Let us continue our discussion with finding paths and labeling layers from outputs to

related inputs in SCG containing acyclic sub-digraphs.

149

As shown in Figure E.2, we take V. as an output, V1 and V2 as inputs.

 

 

0 %% Removed

«rt—Mg

 

 

 

®—=®—*

 

 

 

Figure E.2 Construction of solving order for an acyclic digraph

It is noticed that V3 is on both layer 4 and 5. It is resolved by removing the vertex
having smaller layer number, and the vertices following it. The method can be explained
by looking at the solving order of the output variable V3. Given VI and V2, we can
obtain V3. From V3, V6 is obtained, then V7 and V4. Knowing V7 and V4, V5 is reached

and further V3 is reached.

150

The algorithm designed to obtain the path matrix is as follows:

1)
2)

3)

4)

5)

6)

8)

Reverse the direction of each edge inEG to obtain 56.

Put all source vertices into a set S.

If there is no unlabeled sink vertex, STOP; otherwise, choose a unlabeled sink
vertext,puttintoasetV, setiequaltoO.

Let i = i-l-l and give i as index of v, veV.

If there is a v (eV), so that veS, delete this v from V; if V is empty, goto 3).
Through directed edge, search all vertices u incident to v that veV.

If there are some u’s which are labeled, erase the Old labels.

Put all u’s into V; goto 4).

It must be emphasized that the above algorithm works only on digraph containing no

circuit and self loops. Recall that SCG has every SCC as a vertex (discussed in section

3.2.2.). The SCG satisfies the restriction and, thus, the algorithm can be applied on to

it.

BIBLIOGRAPHY

THE BIBLIOGRAPHY

Allen, R.R., Dubowsky, S., 1977, Mechanisms as Components of Dynamic Systems: A
Bond Graph Approach, Trans. ASME J. Engineering Industry, Vol 99, No.1, pp104-
111.

Bos, A.M. , Tiemego, M.J.L. , 1985, Formula Manipulation in the Bond Graph Modeling
and Simulations of Large Mechanical Systems, J. Franklin Institute, Vol 319, No.1/2,
pp51-66.

Broenink, J.F. , 1990, Computer-Aided Physical Systems Modeling and Simulation: A
Bond Graph Approach, Febodruk, Enschede.

Burreto, J. , and Lefevre, I. , 1985, R-fields in The Solution of Implicit Equations, J .
Franklin Inst., Vol.3l9, No.1/2, pp227-237.

Chua L. O. and P. Lin, 1975, Computer-Aid Analysis of Electronic Circuits: Algorithm
and Computational Techniques, Prentice-Hall Inc.

Close, C.M., Frederick, D.K., 1978, Modeling and Analysis of Dynamic Systems,
Houghton Mifflin Company.

Constantinescu, J. , 1982, Study of the Transient Processes in Large-Scale Power
Systems, Rev. Roum. Sci. Techn-Electrotechn et Energ., Vol.27, No.2, pp211-227.

DeCarlo, R.A., Sacks, R., 1981, Interconnected Dynamical Systems, New York, Mrcel
Dekker.

Dransfield, P. , 1979, Using Bond Graphs in Simulating an Electra-Hydraulic System,
J. Franklin Institute, Vol.308, No.3, pp175-182.

Even, S. , Graph Algorithms, 1979, Computer Science Press.

Filippo, J.M., Delgado, M., Brie, C. and Paynter, H., 1991, Survey of Bond Graph
Theory, Application and Programs, J. Franklin Institute, Vol.328, No.5/6, pp 565-606.

Granda, J.J. , 1984, Bond Graph Modeling Solutions of Algebraic Loops and Differential
Causality in Mechanical and Electrical Systems, Proc. Applied Simulation and Modeling,

151

152
pp188-193, IASTED Conference, San Francisco, CA.

Hamilton, P.S. , 1984, Derivation of the algebraic system Jacobian matrix from bond
graph using a symbol manipulation technique, Ph. D. dissertation, The University of
Texas at Austin.

Hrovat, D., Tobler, W., Tsangarides, M., 1985, Bond Graph Modeling of Dominant
Dynamics of Automotive Power Trains, Dynamic Systems: Modeling and Control,
ASME DSC-Vol. 1,

Karnopp, D.C., Rosenberg, R.C., 1975, System Dynamics: A Unified Approach, Wiley,
New York.

Karnopp, D.C. , 1985 , Bond Graph Models for Electromagnetic Actuators, J . Franklin
Institute, Vol.319, No.1/2, ppl73-182.

Kobayashi, H., Muto, 8., Tamura, Y., Narita, 8., 1978, Decomposition Algorithm for
Determining Sensitivity Constants of Large-Scale Power Systems, Electrical Engineering
in Japan, Vol.98, No.1, pp45-51.

Kokotovic, P.V., Perkins, W.R., Cruz, J.B., and D’Ans, G., 1969, e—coupling Method
for Near-optimum Design of Large—scale Linear Systems, Proc. IEE 116, pp887-892.

Kokotovic, P.V., Khalil, H., O’Reilly, 1., 1986, Singular Perturbation Method in
Control: Analysis and Design, Academic Press Inc. , Orlando, Florida.

Laddle, G.S. , 1975 , Variational Comparison Theorem and Perturbations of Nonlinear
Systems, Proceedings of the American Mathematical Society, 52, pp181-187.

Martinez-Benet, J.M. , Puigianer, L. , 1988, A Powerful Improvement on The
Methodology for Solving large-Scale Pipeline Networks, Computational Chem. Engng. ,
Vol.12, No.2/3, pp261-265.

Paynter, N.M., 1960, Analysis and Design of Engineering Systems, M.I.T. Press,
Cambridge, Massechuset.

Rosenberg, R. C. , 1971, State-Space Formulation for Bond Graph Models of Multiport
Systems, ASME J. Dynamic Systems, Measurement, and Control, Vol.93, No.1, pp35-
40.

Rosenberg, R.C., Karnopp, D.C., 1983, Introduction to Physical System Dynamics,
MCGraw-Hill, New York.

153

Rosenberg, R.C. , 1990, The ENPORT Reference Manual, Rosencode Associates Inc. ,
Lansing, Michigan.

Siljak, 0.0., 1978, Large-Scale Dynamic Systems, North-Holland, New York.

Sobhi, A. , 1985, Symbolic derivation of the state equations and system J acobian using
bond graph, M.S. thesis, The University of Texas at Austin.

Wang, C.M., Jamshidi, M., 1982, A Computational Algorithm for a Class of Large
Scale Nonlinear Time Delay Systems, IEEE, 1982 large Scale Systems Symposium.

William, LR, 1983, Modeling, Analysis, and Control of Dynamic Systems, John Wiley
& Sons, Inc.

Wu, F.F. , 1976, Solution of Large-Scale Networks by Tearing, IEEE Transactions on
Circuits and Systems, Vol.CAS—23, No.12, Dec.

Zhou, T., 1988, Ph. D. dissertation, Michigan State University.

 

NICHIGRN STRTE UNIV. LIBRARIES
IIIIIIIIIIIIII IIIIIIIIIIIIIHIIIIIIIIIIIIIIIIHHI
31293008992798