Continuous System Simulation
Continuous System Simulation
SIMULATION
Continuous System
Simulation
D.J. Murray-Smith
Department of Electronics and Electrical Engineering
University of Glasgow
Glasgow
UK
© 1995 D.J. M u r r a y - S m i t h
O r i g i n a l l y published b y C h a p m a n & H a l l i n 1995
Typeset by E X P O H o l d i n g s , M a l a y s i a
Preface xi
Index 251
PREFACE
1.1 OBJECTIVES IN THE MODELING OF able information and can expose gaps and
SYSTEMS inconsistencies in the existing knowledge
base. Problems of this kind can only be
Much scientific and engineering work resolved by the design of new experiments. It
involves the formalization of knowledge and is clear, therefore, that mathematical model-
the building of models. From experiments ing should not be regarded simply as an alter-
and observations one can form abstract repre- native to experimentation and a solid base of
sentations and laws which can then provide a experimental data is needed in order to
basis for analysis or engineering design. develop a quantitative model. Further experi-
Mathematical models are one particular type mental data are needed in order to carry out
of abstract representation which has been the processes of model validation to establish
found to provide a very effective means of the range of conditions over which the model
encoding information about a real system. may be applied with confidence.
Such representations can also serve to com- The acceptability of mathematical model-
press large amounts of observational data ing, and the extent to which modeling tech-
and allow quantitative relationships to be niques are used, varies considerably from one
established which may not be directly evident field to another. Mathematical models are
from experimental data. used so extensively in the physical sciences
A mathematical model which embraces the and engineering, and are so central to the
essential features of a real-world system per- introductory teaching of many topics within
mits useful analysis to be carried out for the those fields, that, for some at least, there is a
range of conditions for which the model is danger that models can become more import-
believed to be accurate and useful. In this ant than the systems which they represent.
respect a model may be likened to a map: the Often it is only when students are exposed to
scale of the map determines both the accuracy more advanced experimental work in these
and the range of applications that are appro- fields that they begin to appreciate some of
priate. As with a map, some cautious extra- the limitations of the models which they are
polation may be possible with a model. In this using. However, in the biological sciences
case such extrapolation may allow investi- mathematical modeling is still viewed with
gation of situations which cannot be studied some suspicion, although there are many
experimentally because of safety constraints examples in which models have been used
or measurement problems, or because the very effectively to help to solve important
system in question has not yet been biological problems.
constructed. Objectives in the development and applica-
The development of a formal mathematical tion of mathematical Q10dels are very de-
description for a given system normally pendent upon the subject area. Generally,
involves a critical examination of all the avail- within pure science, mathematical models are
2 The principles of modeling
developed for one or more of the following tors, capacitors and inductors, equations
reasons: may be written down from Ohm's law and
Kirchhoff's laws describing currents and volt-
1. hypothesis testing;
ages within a given electrical network as
2. the development of new or improved
continuous functions of time. They are con-
experiments;
tinuous-variable models and involve mathe-
3. the provision of a concise quantitative
matical descriptions based on differential
description which extends the capabilities
equations. This is the general class of model
of the human brain to handle the com-
which is of primary importance in this book.
plexities of the real system under consid-
A second class of mathematical description,
eration.
known as discrete-event models, is also im-
Within engineering, mathematical modeling portant in certain types of computer simula-
is used as a tool for decision making and as a tion application. Changes in discrete systems
basis for design calculations. In this field of take place either periodically or in a random
application, questions concerning model fashion. Between the events which mark these
accuracy and credibility are of central import- changes all the variables of the system remain
ance, since decisions may be made on the constant. One example of a discrete system in
basis of mathematical model predictions which changes take place periodically arises
which have a direct bearing on the safety, in the modeling of a digital computer used
reliability, efficiency or effectiveness of some for real-time signal processing or direct digi-
new product or system. tal control. In this case a continuous variable
In general, whatever the field of appli- is sampled periodically using some form of
cation, mathematical modeling should be analog to digital converter, and calculations
seen as just one aspect of a systems approach are performed on the discrete numbers which
to problem solving. In most cases the steps have been input to the computer from the
involved in the modeling of a process are converter. The system cannot be modeled
linked closely to questions of measurement entirely using differential equations owing to
and experimentation. The integration of mod- the discrete nature of the processing actions
eling, experimentation and data analysis can within the computer and a difference equa-
often provide new insight and allow progress tion type of formulation is necessary for the
to be made which might not be achieved so model of the computer and the associated
readily using one of these three elements in analog to digital converters. In section 3.10
isolation. further consideration is given to the treatment
of discrete models of this kind, in which
changes take place on a regular periodic
1.2 CONTINUOUS-VARIABLE AND
basis. Problems of discrete systems in which
DISCRETE-EVENT MODELS
events occur randomly have led to the devel-
Models commonly used in the physical sci- opment of a separate form of simulation
ences, engineering and biology often involve methodology which is generally known as
variables which are continuous functions of discrete-event simulation. Examples of areas
time. One example is in applied mechanics, of application for discrete-event simulation
where equations of motion may be derived include traffic flow analysis of all kinds. For
from the application of Newton's laws of example, the times of arrival of vehicles at a
motion and involve quantities such as posi- set of road traffic signals and requests for con-
tion, velocity and acceleration which may be nections in a telephone communication
defined as continuous functions. Similarly, in system are both cases of discrete events
modeling electrical circuits involving resis- which can be regarded as occurring ran-
Types of continuous-variable model 3
domly. Event-oriented simulation of this kind variable (the 'output') and another (the
is outside the scope of this book. Further 'input'), there are also internal variables
details of such simulation methods may be which have counterparts in the real system.
found in any of the many available texts Unlike data models, which are based simply
dealing with this specialized area (e.g. refs 1 on observations, system models should incor-
and 2). porate all relevant available knowledge con-
cerning the structure and parameters of the
real system.
1.3 TYPES OF CONTINUOUS-VARIABLE
This book is concerned mainly with models
MODEL
of systems. Models of data may, in some
Mathematical descriptions involving continu- cases, provide a basis for submodels within a
ous variables may be divided, in general larger system description, but for most pur-
terms, into two broad classes. These are poses simulation models of the type being
models of data and models of systems. considered here belong to the class defined
The first of these classes involves math- above as system models. Many subgroups of
ematical descriptions which have been fitted such models can be defined. Some of these
to measured response data. The resulting categories relate to mathematical aspects of
models consist of mathematical functions, the description, while others are concerned
such as polynomials or exponentials, and the more with the field of application.
coefficients or constants appearing in these One of the most fundamental divisions, in
descriptions often have no direct link to rec- terms of mathematics, is between linear and
ognizable elements of the real system. nonlinear system models. In a linear system
Statistical models and other forms of 'black model responses are additive in their effects
box' description come within this class. Such and satisfy the principle of superposition. In
models simply express an observed relation- other words, the output is directly propor-
ship between two or more variables of the tional to the input, and doubling of the input
real system. They have an important role in leads to doubling of the output. In a non-
fields such as control systems engineering linear model the principle of superposition
where input-output descriptions, for example is not satisfied. It is important that questions
transfer functions, can be derived from meas- of linearity should be considered carefully in
urements and can provide a starting point for the initial stages of model development and
design calculations. They tend, however, to that assumptions of linearity should not be
be of less value for the study of problems at made without justification. The range of
the research level since they provide little linear operation of the system, if any such
information about the internal processes of range exists, has to be established by careful
the system. evaluation of what is known about the real
Models of systems, the second of the system or by carrying out additional tests on
classes being considered, are usually devel- the system.
oped by applying established laws and prin- Time invariance is a second important
ciples together with simplifying assumptions property which can provide a basis for classi-
and hypotheses about the structure and func- fication of models of systems. A time-invari-
tion of the system being considered. Such ant system is one in which the observed
models are essentially explanatory in nature performance of the system is independent of
and reflect complex causal relationships. the times at which the observations are made.
Although models of. systems may involve As with linearity, time invariance must be
descriptions in which particular attention is verified in some way before being adopted as
given to relationships between one measured a basis for model formulation.
4 The principles of modeling
Models which are both linear and time careful comparisons with the real system's
invariant are particularly important and behavior, can be of enormous value when
receive considerable attention in most intro- used appropriately.
ductory textbooks in fields such as applied A great deal has been written about mathe-
mathematics, systems engineering, automatic matical modeling methods and principles.
control, electrical circuit analysis and applied Useful texts concerned with general aspects
mechanics. Many practical systems do have of modeling dynamic systems include those
properties which may be approximated in a by Jacoby and Kowalik [3], Jeffers [4] and
satisfactory way by linear time-invariant Spriet and Vansteenkiste [5]. Since this text is
descriptions, even if only for specified and concerned with simulation techniques, rather
restricted sets of conditions. Linear time- than with the fundamentals of modeling, the
invariant models are also very attractive treatment of the model development process
because they can be analyzed using standard is approached in terms of practical issues
mathematical techniques. Nonlinear and which can arise when simulation tools are
time-varying dynamic models present used for developing and evaluating mathe-
far greater difficulty in terms of classical matical models.
mathematical solutions based on paper and The development of a mathematical model
pencil methods. It is for this reason, and for may be described conveniently using the ter-
this very important class of problems, that minology adopted in structured program-
computer simulation techniques have become ming methods in computing science. On one
so important. Computer simulation methods hand, a 'bottom-up' modeling approach
and numerical modeling can provide insight involves starting with a highly detailed
into and solutions to problems for which no description of the system under considera-
analytical solutions are currently available, tion. This incorporates all of the information
except possibly in very special and restricted available about the elements of the system
cases. and allows one to build up a comprehensive
set of equations. On the other hand, a 'top-
down' type of modeling method involves
1.4 MODELING AND SIMULATION
working from a statement of the modeling
PROCEDURES FOR CONTINUOUS-
objectives and leads to the formulation of a
VARIABLE MODELS
model which may involve a significantly
A critical assessment of well-established smaller set of equations. The resulting model
models, which have become accepted and has to be judged in the context of the in-
widely used, shows that they were usually tended application. Although the form of a
developed for specific applications. Such top-down model may be relatively simple at
models have clear boundaries and the under- the initial stages of development, additional
lying simplifying assumptions are well docu- features may be added if the first assessment
mented and fully understood by those using suggests that the model is deficient in some
the models. On one hand, if there are no state- significant way. Modifications are made, in
ments available about the modeling assump- an iterative fashion, until the model is
tions and about the model boundaries, any deemed adequate for the intended applica-
results obtained from a model are of limited tion and can be used with some confidence.
practical value and may possibly be very As already pointed out, useful analogies
misleading. On the other hand, a well- can be drawn between models and maps. A
documented model, which has been devel- political map which shows cities and national
oped for a particular type of application, boundaries may be very useful for many pur-
using an interactive approach involving poses but is of no value if one is interested in
Procedures for continuous-variable models 5
the physical geography of the same area. In unlikely to be adequate for AC studies in the
the latter case one needs a map which shows radio-frequency range. Eventually a point is
all the physical features of the region in terms reached in the modeling of a system of this
of mountain ranges, rivers and lakes. Political kind at which the whole concept of a lumped
boundaries are of secondary importance in parameter model becomes inappropriate and
this case. A combined map which includes it becomes necessary to use a distributed
both physical features and all the information parameter model based on partial differential
on the political map may well be compre- equations. One could then ensure that the
hensive but may not always be the most satis- description of the circuit is sufficiently
factory representation for a given application. general to cover the complete range of fre-
Similarly, the intended application of a model quencies that could possibly be of interest.
has a major bearing on the form which the However, the effort involved in this process
model should take and the aspects of the of developing a truly general model can be
system which are to be included in the mathe- considerable, and a compromise is often pre-
matical description. A fully comprehensive ferred in which the model has a well-defined
bottom-up type of model is not always best. range of validity and applicability. Models
In engineering, models are often developed used outside their valid range give inaccurate
using bottom-up design principles since and sometimes highly misleading predictions
detailed information is available about the of system performance.
characteristics of each element of the system One obvious problem with the develop-
and well-established physical laws and prin- ment of models using the bottom-up approach
ciples can be applied to the formulation of the is that it is difficult sometimes to obtain reli-
model equations. A good example of this is able information about the elements or sub-
the development of a mathematical model for systems. Major simplifying assumptions may
an electrical circuit involving lumped passive have to be made to obtain any kind of quanti-
elements such as resistors and capacitors. tative description, and this is particularly true
Well-known principles of circuit analysis in fields such as biology. In such cases there
(such as Kirchhoff's laws and Ohm's law) may also be problems in assessing the range
provide a simple means of establishing a set of applicability of the model and in carrying
of ordinary differential equations for a given out detailed comparisons with the real system,
network of circuit elements. The resulting since the number of variables in a typical
model may be regarded as valid provided it model may be much greater than the number
is used within the range of frequency and of quantities directly measurable. The physio-
applied voltage for which the description for logical mechanisms of oxygen and carbon
each circuit element is valid. In this respect it dioxide exchange in the human body provide
is important that the person developing the an interesting application which is discussed
model should be aware that the use of Ohm's in Chapter 12. In models of this highly
law in describing a lumped circuit element, in complex system the lungs and the body
terms of the single property known as resis- tissues are conventionally represented by
tance, involves implicit assumptions of linear- uniform compartments. Variables such as the
ity and frequency independence. The modeler concentration of each gas in the lung compart-
must remember that a simple circuit element ments or the blood gas concentrations within
such as a wirewound resistor has other the tissue compartments may not be directly
important properties as well as resistance; the measurable in the real system, and one may
effect of inductance will become significant as therefore have only limited scope for compar-
the frequency is increased and a model devel- ing these aspects of the model with the real
oped for use in DC operating conditions is system. Top-down modeling may, however,
6 The principles of modeling
be quite appropriate for dealing with situa- could be regarded simply as an isolated elec-
tions such as this, with simple models being trical load. Similarly, for a thermal system,
postulated initially and additional features boundaries may be defined at which tempera-
being added as new experimental evidence is tures are constant or independent of heat
accumulated. flow. Such a representation clearly requires
The 'principle of parsimony' is often careful consideration to be given to inputs
quoted in texts on mathematical modeling. that originate outside these predefined
This states, essentially, that models should be boundaries and may have significant effects
kept as simple as possible. The concepts of upon the system within. This may even
top-down modeling are generally consistent require the use of measured quantities from
with this principle. Certainly the inclusion of the real system being converted to digital
unnecessary detail in models is undesirable form and used as variables within the simu-
and leads to an excessive number of para- lation model. This is a topic of particular sig-
meters. Such overparameterization can cause nificance in real-time simulation applications
difficulties in the validation stages of model involving a combination of real-system hard-
development. ware and a simulation model. Examples of
Whatever the approach adopted, it is the use of auxiliary inputs of this kind can be
important to distinguish between simplifi- found in aerospace engineering.
cations involving the elimination of elements Validation is an important and integral part
which have a negligible effect on the overall of the modeling process. In some respects
behavior of the system and those introduced 'model calibration' is a better description of
for mathematical convenience or computa- this aspect of model development. No model
tional speed and efficiency. Although the can ever be completely validated and shown to
choice between complex comprehensive be a perfect representation of the system under
models and simpler models developed for a investigation. However, the performance of a
specific application is often controversial, it is simulation model can be evaluated in terms of
now generally accepted that both the top- equilibrium conditions, stability limits, tran-
down and bottom-up approaches are valu- sient response or frequency response. Findings
able. They complement each other and a of this kind may then be compared with infor-
useful combination of the two approaches can mation from the real system and any differ-
sometimes be achieved by splitting a complex ences assessed in the context of the intended
model into submodels which can be studied, application of the model. Limitations of the
developed and validated independently model can thus be established in terms of that
before being combined. application and, once a model has been 'cali-
Careful consideration of the purpose of a brated' in this sense, it can be used within its
model together with a full and detailed known limitations until new evidence from the
understanding of the real system are essential real system shows that the model is in some
in defining model boundaries. For example, way inadequate for the intended use and
in modeling an electrical generator considera- should be modified. Figure 1.1 is a flow chart
tion has to be given to the representation of which illustrates some aspects of the iterative
the network to which the generator is con- process of model development, testing and
nected. If the details of the current distri- modification.
bution in different parts of the network are of In considering the use of measured system
interest, or if the network incorporates other response data for the testing and calibration of
generators, it may be necessary to include a simulation model it should be noted that the
much more within the boundaries of the model development process may sometimes
model than would be the case if the network involve inverse modeling techniques in which
EXPERIMENTATION MODELING AND SIMULATION
t
Define the measurements
necessary to characterize
the system behavior and
Define the model
boundaries r-
~
define the available inputs Define the assumptions and
approximations to be used
rl Design appropriate
experiments
I
Apply physical laws and
principles to formulate
the model equations
Carry out experiments and Estimate values of any poorly-
r+ make measurements defined parameters for the l+-
I chosen model structure
Create a database of Formulate the simulation model
measured responses
I
Implement the simulation
model on a computer
~
Carry out simulation experiments
equivalent to the experiments
on the real system
J
results inconsistent l r~lts incQllsistenl
Compare model and real system
for chosen data sets
resylts satisfactoOl
results inconsistent r~lts mk!lWlistent
Test predictive capabilities
of the model using additional
data sets
results satisfadQot
Apply the model in this form
until new evidence suggests
need for changes
Fig. 1.1 Block diagram of the modeling process showing the links between experimentation and the
stages of model development and checking. In situations in which active experimentation is impossible,
current observations and historical data from the system under investigation may be useful for model-
testing purposes.
8 The principles of modeling
estimates of the structure and parameters of a dinated use of experimental techniques and
model are obtained from measured data. This computer simulation can often provide
involves the use of system identification tech- insight which would be difficult, or impos-
niques, as outlined in Chapter 9 in the special sible, to obtain by conventional experiments
context of model validation. Although details alone. Examples of the use of simulation in
of system identification methods fall beyond this way for the testing of scientific hypo-
the scope of this text, it is appropriate that theses can be found in many different fields.
those using simulation techniques should be Problems in biomedical engineering, espe-
aware that such methods exist and that they cially the application of control engineering
should have some understanding of their concepts to the investigation of physiological
potential. In addition to providing a basis for control systems, provide interesting examples
one approach to model testing, they provide of this. The respiratory control system, which
an approach to the problem of establishing involves the pulmonary gas exchange pro-
values for model parameters which are not cesses already mentioned, is one complex
known from theory. Examples of fields in physiological system which has been studied
which system identification techniques are extensively using simulation techniques to
recognized as established tools for model test hypotheses. Respiratory control was in
development and calibration include the fact one of the first biomedical engineering
development of ground-based aircraft flight problems to which continuous system simu-
simulators and the development of simula- lation methods were applied. Although this is
tions for biomedical applications such as a biological problem, the main features of this
anesthesia. application are typical of those encountered
In the modeling of engineering systems it is in the modeling of physical and chemical
often possible, and usually desirable, to estab- processes in other application areas.
lish a close link between elements of a model The respiratory control system can be
and the corresponding components of the real divided conveniently into a 'plant', which
system under investigation. However, even in comprises the gas transport and gas exchang-
some engineering applications, this cannot be ing processes of the lungs, circulation and
done initially because of uncertainties in the tissues, and a separate 'controller', which
quantitative information available about the involves various chemoreceptors and other
system. One objective in applying system sense organs which provide the neural 'drive'
identification techniques in such cases may be to the muscles responsible for breathing.
to help select one representation from a Figure 1.2 shows a much simplified block
number of candidate models. This is closely diagram of the system as it is understood at
linked to the idea of model calibration and present.
emphasizes yet again the fact that the devel- Many early simulations of the respiratory
opment process involves a sequence of steps control processes could provide reasonably
and may require several iterations as the accurate predictions of some aspects of the
model is tested and refined. system's performance but some difficulties
The type of application for which a model remained. One area which presents particular
is being developed always has a considerable problems is the respiratory center, which
influence on the form of the model itself. For combines and processes signals from the
example, in a research environment, mathe- various sensory receptors within the system.
matical models may be studied using simu- Indeed, most studies have avoided many of
lation, together with experimental results, as the complications of the respiratory center
part of a process of testing hypotheses about because detailed quantitative data about
a real system. Results obtained from the coor- some aspects of this part of the system are
Procedures for continuous-variable models 9
BRAIN
COMPARTMENT
I_ feedback pathways
I from
I chemoreceptors
~-PT:V --: 1
CONTROLLER r4- - -
I
I
-f - - - - - - - 1-I
I I
feedback pathways
I I drive signals to
muscles of chest,
from I diaphragm etc.
mechanoreceptors
in chest etc. I
I
I
LUNGS
circulation
---
COMPARTMENT
REPRESENTING OTHER
BODY TISSUES
Fig. 1.2 Highly simplified block diagram of the respiratory control system, which is responsible for the
regulation of breathing. Two types of feedback pathway are shown in this diagram: (1) feedback to the
respiratory centre from chemoreceptors which sense concentration levels of oxygen and carbon dioxide
in the blood and tissues and (2) feedback from other receptors associated with the lungs and chest, such
as the mechanoreceptors responsible for sensing the extension and velocity of the muscles which are
active during breathing.
10 The principles of modeling
still lacking. However, a number of models trast to a more detailed representation in
have been developed which have provided a which the neuromuscular elements are
means of assessing competing hypotheses included. In making decisions about model
concerning interactions between the neural boundaries such as these it is essential to keep
and chemical components of respiratory in mind the intended application of the model
control. Many problems are still unresolved, and the measurements likely to be available
especially in terms of the response of the for model testing and validation. One impor-
system to exercise, and simulation techniques tant factor is that some system variables, such
continue to have an important role in this as gas concentrations at the mouth and gas
field of bioengineering research. flow rates, are accessible for direct measure-
The respiratory system provides other ment on a continuous basis. Other variables,
useful illustrations of some of the modeling such as blood gas concentrations, are less
and simulation problems which have already readily measured.
been touched on in this introductory chapter. While some fundamental physiological
The first of these concerns model boundaries. problems concerning respiratory control
Physiological systems are highly interactive, remain unsolved, respiratory system models
and it is seldom possible to consider one are now beginning to receive attention in the
element in isolation. The study of the respira- context of clinical research and in the design
tory system requires understanding of the gas of medical hardware to provide breathing
transfer and exchanging processes in the functions artificially. Even quite simple simu-
airways, lungs, blood and tissues. The ques- lation models of the respiratory control
tion of boundaries is therefore a complex one system can exhibit maintained oscillations
which must be reviewed constantly through- which appear as periodic patterns of breath-
out the model development process. What ing similar to well-known clinical conditions.
level of complexity is appropriate for repre- Although the above example, which has
sentation of the airways and gas exchanging provided a convenient preliminary illustra-
regions of the lungs? How many lung com- tion of some of the processes of development
partments should there be? Are these lung of a simulation model, relates to a biological
compartments in series or in parallel? What system, the principles are quite general.
level of complexity is needed to represent the Decisions also have to be made about bound-
circulation and body tissues? How many aries in approaching the modeling of physical
lumped tissue compartments should be systems. The model structure is often no more
included in a model intended for respiratory obvious in complex engineering systems than
control system studies? How should breath- it is in this biomedical example. Similarly, the
ing be represented? Do we include the action availability of measurable output quantities
of the respiratory muscles or do we draw the for the purposes of model testing and valida-
model boundaries so that ventilation (the tion may often present problems in engineer-
flow of gas into and out of the upper airways) ing applications which are just as difficult to
is a model input? In the last case the neuro- solve as in the biological case. One vitally
muscular drive may be represented in an important point is that, whatever the field of
implicit fashion by a quantity which alters the interest, the development of a simulation
amplitude and frequency of a periodic wave- model must be undertaken in the context of a
form representing ventilation, and a consider- particular application and must be an itera-
able simplification results in terms of the tive process involving a rigorous testing and
structure and complexity of the model in con- evaluation stage.
References 11
REFERENCES Mathematical Modelling with Computers,
Prentice- Hall, Englewood Cliffs, NJ.
1. Neelamkavil, F. (1987) Computer Simulation 4. Jeffers, J.N.R. (1982) Modelling, Chapman &
and Modelling, John Wiley, Chichester. Hall, London.
2. Pritsker, A.A.B. (1984) Introduction to Simulation 5. Spriet, J.A. and Vansteenkiste, C.c. (1982)
and SLAM II, 2nd edn, John Wiley, New York. Computer-aided Modelling and Simulation,
3. Jacoby, S.L.S. and Kowalik, J.S. (1982) Academic Press, London.
AN INTRODUCTION TO SIMULATION
METHODS 2
2.1 THE NEED FOR SIMULATION Analog principles have been very widely
used for a very long time for the investigation
Simple mathematical solutions can be found of complex nonlinear phenomena. For
for linear ordinary differential equations example, active electrical networks which are
having constant coefficients. Such equations analogs, or models, of some other equivalent
form the basis of many well-known and very physical system can be established readily.
useful lumped parameter models. Intro- Conclusions reached on the basis of simula-
ductory courses in dynamics and in electrical tion experiments on the electrical network
circuit theory provide many illustrations of should then apply to the original system of
models of this type. For example, transients in interest. The use of an analog of this kind
simple electrical circuits involving linear ele- makes it relatively easy to set up the modet
ments such as ideal resistors, capacitors and and to carry out experiments which might be
inductors can be studied very easily using impossible, or at least difficult and time-con-
standard analytical tools for the solution of suming, on the real system. Wind tunnel
ordinary differential equations. Concepts such models and models used in ship testing tanks
as the complementary function, particular provide two other examples of this kind of
integrat D-operator and Laplace transform approach based on the study of analogous
are of considerable importance for problems systems.
of this type, but such analytical tools are of no With conventional digital computers,
use in the case of many practical problems advantage may be taken of the many
involving nonlinear models. For example, general-purpose computer-based numerical
analytical methods are of little assistance if an tools which are now available for continuous
inductive element in an electrical circuit has a system simulation. Although, in principle, any
core of magnetic material and displays general scientific programming language may
significant hysteresis. Similarly, in dynamics, be used for the numerical solution of differen-
linear mechanical systems involving elements tial equations, the use of special-purpose
with mass, spring stiffness and viscous resist- simulation software tools can be very helpful
ance provide a basis for much elegant mathe- in many practical situations and can eliminate
matical analysis which can provide valuable many of the computational difficulties associ-
understanding. Replacement of a viscous ated with writing a simulation program.
damping element by an element with static These special-purpose tools are particularly
friction immediately makes the problem non- appropriate when used in an interactive mode
linear and eliminates the possibility of apply- and allow the user to concentrate attention on
ing standard mathematical techniques. In the model instead of having to struggle with
cases where analytical methods are impracti- computational methods and the technical
cal the only possible approach is through problems of getting meaningful numerical
numerical techniques of some kind or by the results. Simulation is essentially a process of
use of an analog method. experimentation with mathematical models
14 An introduction to simulation methods
and in any experimental situation good inter- train performance with different loads and
action between the experimenter and the gradient profiles.
experiment itself is of vital importance. Good The development of the electronic analog
simulation tools provide good facilities for computer during the period from about 1945
user interaction. to 1955 provided a vitally important stimulus
Simulation software tools for distributed to the field of simulation and introduced tools
parameter models, involving partial differen- that could be translated into hardware which
tial equations, are less well developed than was sufficient reliable and accurate for use in
those for lumped parameter descriptions, large-scale engineering studies. Examples of
involving ordinary differential equations. the use of electronic analog simulation
Although some ,simulation tools do exist for methods can be found in many accounts of
problems involving some special classes of aerospace projects in the 1960s, in the nuclear
partial differential equations, this book deals industry at that time, in the field of chemical
mainly with simulation applications involv- process control and in many other areas of
ing lumped parameter models. engineering. The essential idea in the elec-
tronic analog computer is that electrical volt-
ages or currents can be used to represent
2.2 METHODS OF SIMULAnON
other physical variables in the system under
Simulation may be performed using special- study. The electronic operational amplifier
purpose hardware, or using a general-purpose can form a computational block, which can
computer. Both of these approaches are impor- either take the algebraic sum of these electri-
tant, and the simulation methodology adopted cal quantities or act as an integrator, again
depends very much upon the application area. operating on the electrical analogs of the real
physical variables of the original problem.
Other electronic elements can act as units to
2.2.1 TECHNIQUES BASED ON
multiply variables by constants or to form the
SPECIAL-PURPOSE SIMULA nON
product of two variables. Similarly, there are
HARDWARE
active electronic analog elements which can
The earliest developments in continuous- carry out tasks such as division of one vari-
system simulation were based mainly upon able by another, the determination of a square
special-purpose simulation equipment such root, the representation of some arbitrary
as a mechanical differential analyzer. This function of a given variable or simple logical
type of mechanical computer system involved operations such as comparing the magnitude
operations which were based entirely on of one quantity with another. The intercon-
mechanical principles. For example, the oper- nections between the operational elements are
ation of integration could be based on the fact established by the user, and it is these inter-
that angular position is the integral of angular connections that determine the set of equa-
velocity. Although the mechanical type of dif- tions represented on the analog computer.
ferential analyzer was cumbersome and Considerable ingenuity was shown in the
difficult to use, it was applied to the success- design of commercially available analog com-
ful solution of some very important problems puters to ensure that the user interface was
both prior to and during the Second World suitable and convenient. Much emphasis was
War. Early published accounts of practical also placed on the presentation of results in
uses of this type of simulation computer graphical form. One very important point
include descriptions of applications involving about the analog computer is that such a
the prediction of target position in gun- system is inherently parallel in terms of its
aiming mechanisms and the calculation of internal organization. Since every operation,
Methods of simulation 15
such as integration, multiplication or sum- pure delays and other operations which were
mation, is carried out using independent difficult to implement using analog hardware
hardware elements, the individual processes alone. The principles of hybrid computation
required to form a solution are performed led to the development of some very power-
simultaneously. This gives the electronic ful systems, many of which were-used in
analog computer important benefits in terms dedicated applications within the defense
of speed of operation since the overall speed industries in the 1960s and 1970s.
of the simulation is independent of the size of Although analog and hybrid computers
the problem. The main practical limitation in have ceased to be important as general-
terms of speed of operation of an analog com- purpose simulation tools, small analog
puter is imposed by the finite bandwidth of systems can still be very useful in the solution
the active electronic components, such as of real-time simulation problems. This is
operational amplifiers. In principle, when especially true of simulation studies which
designed with high-bandwidth components, involve real hardware and a simulation
an analog computer can provide performance, model working together in an integrated
in terms of speed, which is still hard to match fashion. Such situations arise frequently in
in simulation applications by even the most the testing of engineering systems. One very
powerful digital computers available today. good example is the evaluation of a control
Problems with analog simulation tech- system which incorporates digital processors
niques, as implemented using 1960s tech- within a feedback structure. The control hard-
nology, lay in the high cost of the hardware ware and software might be tested initially in
and the fact that such computers had no the laboratory using a simulation of the plant
applications other than simulation. Other prior to being installed on site. In this type of
difficulties, which were significant in some situation, interfacing the control system to an
applications areas, included the fact that the analog computer simulation can allow checks
programming of analog simulation problems to be made in a realistic way and may help to
was time-consuming, relatively complicated [Link] errors in the system prior to com-
and required specialist skills. In addition, the missioning on site. Similarly, in laboratory
relatively low accuracy of analog computers work of various kinds, small analog comput-
was an important factor which counted ers (or equivalent special-purpose electronic
against this technology for some types of hardware based on operational amplifiers)
application. The low reliability of early can provide a flexible and low-cost way of
analog computer hardware also raised impor- applying control or carrying out signal pro-
tant issues in terms of maintenance and oper- cessing operations such as multiplication of
ating costs. signals or filtering. Analog methods are
In order to try to overcome the acknow- discussed in more detail in Chapter 14.
ledged limitations of analog computers, the Some developments in the 1970s, such as
manufacturers of simulation computers those at the Technical University of Vienna,
decided to try to combine the best features of where a very powerful multiuser hybrid sim-
conventional serial digital computers with the ulation facility was developed [2], demon-
best features of analog systems [1]. This led to strated the potential of hybrid principles.
what became known as the hybrid computer. However, by that time, major developments
The digital part of a hybrid system was gener- were taking place in the fields of digital hard-
ally used to control the operation of the simu- ware and software. Conventional digital
lation on the analog part of the system, to set computers were becoming ever cheaper and
up and check the analog elements and to more powerful and new concepts were
provide capabilities for function generation, leading to specialized digital systems which
16 An introduction to simulation methods
had architectures designed more specifically Project Whirlwind computer at MIT was
for simulation. This allowed digital machines developed for real-time simulation of flight
to be used in real-time simulation problems, vehicles. Calculations in these machines were
such as flight simulation of aircraft or missile all based on fixed-point arithmetic but, by
systems, in place of analog or hybrid comput- modern standards, were extremely slow. It
ers. Such specialized digital systems could be is worth noting that, although the first gener-
significantly less costly than a large conven- ation of electronic analog computers could
tional digital computer of power sufficient to simulate aircraft and missiles in real time,
deal with the specialized demands of this digital computers could not approach this
type of application and also easier to use and level of performance until about 10 years
less expensive to maintain than most hybrid later. Subsequent developments in digital
systems. computer technology led to very significant
Dedicated special-purpose digital systems increases in speed and the introduction of
continue to be used for some simulation high-level programming languages, such
applications. Considerable efforts have been as Fortran, made digital computers more
made in recent years to produce software accessible to many potential simulation users.
which is appropriate and allows users to For lumped parameter models, traditional
exploit fully the potential of such specialized continuous system simulations are essentially
systems. For example, Applied Dynamics numerical solutions of sets of ordinary differ-
International have produced a series of ential equations. The programming task can
special-purpose digital systems for real-time be approached in a number of ways, includ-
applications and other numerically intens- ing coding from first principles in a general-
ive simulation tasks. Such special-purpose purpose high-level language such as Fortran,
hardware systems require special-purpose Pascal or C. An alternative is to write the
high-level software to allow full advantage program in a general-purpose language but to
to be taken of the speed capabilities which make as much use as possible of an appropri-
are available [3]. Many highly specialized ate library of standard numerical methods
and fully dedicated real-time simulation subprograms such as the NAG [4], IMSL [5] or
applications have been based on computer 'Numerical Recipes' [6]. One important
systems such as these. Other real-time tasks, benefit from using a well-established sub-
such as those involved in piloted flight simu- routine library is that the numerical methods
lation, have also been approached using routines are generally well documented and
multiprocessor architectures based upon much experience has been accumulated in
more conventional computer hardware. their use for many types of application.
These and other aspects of this specialized Provided that the available routines are used
area of simulation are discussed in more with care, and within their published limita-
detail in Chapters 14 and 15. tions, they offer a reasonably efficient and
effective means of generating simulation pro-
grams using a general-purpose high-level pro-
2.2.2 TECHNIQUES BASED ON
gramming language. This approach has been
GENERAL-PURPOSE COMPUTER HARDWARE
used in some special-purpose simulation soft-
Some of the earliest developments in digital ware packages. These consist of a library of
computer technology were associated with routines commonly required for simulations
problems which would now be regarded as applications, and the user must supply a main
simulation applications. The US Army program which defines the model and experi-
EN lAC computer was used in the late 1940s ment, with appropriate calls to the numerical
to solve ballistic trajectory problems, and the subprograms as and when required.
Methods of simulation 17
Simulation languages differ from the simu- provided for interaction with the simulation
lation packages described above in that they problem while the computer was running.
involve coding at a higher level. A simulation The presentation of results was, at best, in the
language should allow the user to write a form of simple graphs, the particulars of
simulation program using a natural form of which had usually to be defined prior to the
application-oriented language and provide a start of a simulation run.
more user-friendly form of environment than Developments in the power of personal
is possible for the nonspecialist programmer computers and workstations has changed the
using general-purpose programming lan- situation considerably in terms of the user-
guages and associated numerical libraries. computer interface and has made truly inter-
The development of specialized simulation active digital simulation possible for many
software for general-purpose digital com- application areas. The availability of the
puters began as a conscious attempt to windows type of environment and of mouse-
represent, on a digital machine, the basic driven software has led to a revitalization of
operations which were possible on an analog block-oriented methods as an alternative to
computer. Examples of early attempts to the equation-oriented approach. The use of
produce software capable of carrying out modern post-processing techniques also
such digital simulations include the DAS offers far more flexibility for the presentation
program, which was followed by MIDAS. of graphical results on a personal computer
Brennan and Linebarger [7] have provided than was available in any of the early simula-
an interesting review of these very early tion languages running on large mainframe
developments in simulation software. These computers. The move from mainframe
tools were generally block oriented and the systems and multiuser minicomputers to
user built up a program to represent, on single-user personal systems has thus
the digital machine, the function of the equiv- brought major benefits in terms of the possi-
alent analog computer diagram showing the bilities for efficient interaction between the
interconnections between the operational user and the simulation problem.
units. It is now generally accepted that a good
MIMIC was an early language which simulation environment should allow models
departed from the block-oriented tradition to be set up and modified quickly and easily
and adopted an equation-oriented format and should permit efficient interactive experi-
which was followed, with variations, by mentation. It should be capable of releasing
many other languages. This approach simulation model developers and users from
involved defining the simulation model much of the detailed coding, and thus allow
directly in terms of the differential equations them to use their time effectively by concen-
of the underlying mathematical model. In trating attention on the model and its behav-
1967 the Society for Computer Simulation, ior. While incorporating many powerful
which at that time was known as Simulation facilities, user-friendly software should gen-
Councils Inc., published a proposed standard erally keep advanced features hidden from
for continuous simulation languages [8], the user until they are actually needed. A
which has formed a basis for many subse- good software system should thus be seen as
quent developments. This standard, which is an appropriate tool by users at many differ-
based upon equation-oriented principles, is ent levels of experience and expertise.
discussed in more detail in Chapter 6. The effectiveness of an interactive simula-
Much of the early software was batch ori- tion depends critically upon whether or not it
ented and relatively inconvenient to use in a is compatible with the typical timescale of
simulation environment. Few facilities were human activities. User-friendly simulation
18 An introduction to simulation methods
software should provide fast keyboard in testing hypotheses or as a means of
responses, fast program translation and fast expressing all of the available knowledge
solution of the model equations. Any exces- about the properties of a complex dynamic
sive delays at the translation stage or at system in a concise and consistent fashion.
runtime can greatly reduce the productivity of The process of building up a simulation
the computer system and of its user and have model can also highlight gaps in existing
an adverse effect on his or her creativity. For knowledge and lead to a search for additional
this reason modern stand-alone workstations information.
and personal computers (especially when It is useful to distinguish between a simula-
equipped with a numeric co-processor) often tion run, which is a term used to describe a
provide a much more effective environment single solution of the model equations, and a
for interactive simulation than multiuser time- simulation study, in which several simula-
sharing computer systems. tion runs are performed and further analysis
In addition to general-purpose simulation is carried out on the basis of the multiple run
languages and packages, a significant number results. Modern commercial simulation soft-
of special-purpose software tools have been ware systems are increasingly being designed
developed for particular applications areas. to provide facilities for simulation studies and
These include software for simulation of elec- have powerful post-processing and graphics
trical and electronic circuits, software for presentation options available or can be
chemical process simulation, software for linked to appropriate computer-aided analy-
simulation of electrical power systems and sis and design packages. Simulation studies
software for simulation of hydraulic and may, for example, be concerned with model
fluid-power control systems. None of these optimization, with sensitivity investigations
specialized aspects of simulation software is to establish the effects of a change in a para-
considered in detail in this book. Further meter, or perhaps with an investigation of
information about these applications can be conditions under which the model becomes
found in more specialized texts (e.g. refs 9 unstable. Although such studies are often
and 10). numerically intensive in nature, and require
high speeds of solution for any interactive
investigation, results obtained by manipulat-
2.3 A REVIEW OF SIMULA nON
ing a simulation model can often provide evi-
APPLICA nONS
dence which would be difficult to obtain by
It has already been pointed out that simula- conventional experimentation alone.
tion may be viewed as a process of experi- Figures 2.1 and 2.2 show flow charts illus-
mentation with mathematical models. This is trating two ways in which models can be
true in the context of both scientific and en- used in the testing of hypotheses [11]. In
gineering applications. In both these fields, Fig. 2.1 a new finding from an experiment on
simulation techniques are especially useful the real system leads to the formulation of a
when the models are nonlinear in form or are hypothesis. Investigations based on a simula-
to be used in real time. In such cases there is tion model and the experimental results,
seldom any alternative to the use of a numeri- taken together, may provide support for the
calor analog type of solution. hypothesis or lead to its modification or rejec-
Simulation applications are very varied, tion. Rejection of the hypothesis may lead to
and the precise approach adopted depends modification of the model or to a new
on the objectives and purpose of the investi- hypothesis which is, in turn, investigated in
gation. For example, models used in scientific the same way. In Fig. 2.2 the hypothesis is
investigations are often developed as an aid based on findings from experiments on the
A review of simulation applications 19
Experiment
t
Examination of experimental
results in terms of
qualitative model
~
Formulation of hypothesis
~
Development of quantitative
mathematical model
+
Correlation of experimental findings
and hypothesis
using mathematical model
I I
+ t
~ Hypothesis supported I Hypothesis rejected I
Fig. 2.1 Flow chart illustrating processes involved in the formulation of a hypothesis from new experi-
mental or observational findings.
simulation model. The model is, in this case, even using a desk-top personal computer of
being used predictively, and new experi- relatively modest power. Further enhance-
ments are needed to provide supporting evi- ment of the interactive capabilities of simula-
dence from the real system. Again, the tion systems is still desirable in order to
process may lead to modification of the create an environment in which complex
model or to a new hypothesis. Practical inves- models can be developed and optimized in a
'tigations can often involve a combination of more natural and user-friendly way. This also
the approaches outlined in these two flow implies careful integration of a number of dif-
charts. In an iterative process of this kind ferent systems, such as databases for experi-
there is a need for a man-machine interface mental records and simulation results,
which allows a close and natural form of graphics packages for the presentation of
interaction to be established between the sim- results, analysis software for any associated
ulation model and the user. Modern simula- numerical operations, as well as the simula-
tion software and hardware can provide an tion software itself. Developments in expert
environment in which many problems can be systems and other areas of artificial intelli-
handled efficiently in an interactive fashion, gence may also be of importance in providing
20 An introduction to simulation methods
Simulation experiment
on mathematical model
Assessment of model
behavior
Formulation of hypothesis
New experiment
I I
t t
I Hypothesis supported ~ Hypothesis rejected I
Fig. 2.2 Flow chart illustrating processes involved in the formulation of a hypothesis from examination
of model behavior.
a basis for intelligent help systems and deci- extract new information which may be incor-
sion aids to guide the less experienced user. porated into an improved model and thus
Simulation may also be required to provide simulation results which are in closer
compare possible experimental strategies or accord with reality. In this type of application
to help with the design of experiments on the simulation again forms an element which is
real system to maximize the information central to the overall scientific method.
content of extracted data. From optimized In engineering design applications simula-
experiments of this kind it may be possible to tion has a role which is often predictive.
References 21
Simulation may be used because no possibili- assessed in this way is by means of a nonlin-
ties exist for experimentation on the real ear simulation model. In the case of the heli-
system. This may be because the system is still copter application this nonlinear simulation
at the design stage or because experiments study would normally include assessments
may not be permissible because of operating carried out using real-time methods in a
requirements or safety constraints. ground-based piloted flight simulator facility.
The field of automatic control systems pro- Simulation models have another type of
vides many good examples of engineering application that is important in engineering
applications in which simulation has a central and in many other fields. This concerns edu-
role within the design process. Control cation and training. Many systems are both
systems are normally designed on the basis of dynamic and complex, and it is often difficult
a linearized description which is valid only to gain an understanding of them from a
at a single given operating condition. written description or from an experimental
Assessment of performance, on the other investigation. A well-planned simulation
hand, frequently requires investigation of the exercise carried out in parallel with experi-
performance of the complete control system mental studies can provide detailed under-
over a range of operating conditions. This is standing where the other methods, used
necessary in order to assess the overall alone, may be more limited in potential. For
robustness of the control system design and example, physiological systems are both
to be able to make predictions regarding per- complex and dynamic. A simulation model of
formance limits. a physiological system, which has been tested
The design of flight control systems for a adequately and is accepted as a credible rep-
helicopter is an interesting example in which resentation, can therefore be used as an edu-
the dynamics of the vehicle changes consider- cational tool. Students can alter the model,
ably with the flight condition. A control often in ways which would be impossible
system designed on the basis of a linearized using conventional experimental methods,
vehicle model to meet particular performance and immediately observe what happens.
requirements in terms of lateral and longitu- Many examples of simulation are now to be
dinal dynamic response at hover could have a found in medical education. The MACPUF
very different performance for fast forward simulation is one of the best known [121.
flight. To be able to undertake an investigation
of system performance in such cases it is not
usually enough to test the control system with REFERENCES
a series of linearized models of the system
1. Bekey, G.A. and Karplus, W.J. (1968) Hybrid
being controlled. Clearly, a new linearized Computation, John Wiley, New York.
description could be used for each flight con- 2. Kleinert, W., Solar, D. and Berger, F. (1983)
dition and the designer could obtain some Status report on TV Vienna's hybrid time
information about the robustness of the com- sharing system, in Proceedings of the 1st European
plete system from linear analysis applied at Simulation Congress, September 1983, Aachen
each point. However, the results would Onformatik-Fachberichte 71), (ed. W. Ameling),
provide information only about the dynamics Springer, Berlin, pp. 193-200.
3. Van Schieveen, H.M. (1986) The System 100-
of the complete system for small perturbations an efficient hardware/software architecture
about each nominal flight condition. Often for real-time and time-critical simulation, in
one is also interested in the response for large Proceedings of the 2nd European Simulation
disturbances or large command signals, and Congress, September 1986, Antwerp, Simula-
the only way in which a design can be tion Councils, Ghent, pp. 413-19.
22 An introduction to simulation methods
4. NAG (1987) Fortran Library Manual, Mark 12, 9. Bowers, J.e. and Sedore, S.R. (1971)
NAG, Oxford. SCEPTRE. A Computer Program for Circuit and
5. IMSL (1987) Library Reference Manual, Version Systems Analysis, Prentice-Hall, Englewood
1.0, IMSL, Sugar Land, TX. Cliffs, NJ.
6. Press, W.H., Flannery, B.P., Teukolsky, S.A. 10. Spriet, JA and Vansteenkiste, G.e. (1982)
and Vetterling, W.T. (1986) Numerical Recipes, Computer-aided Modelling and Simulation,
Cambridge University Press, Cambridge. Academic Press, London, Chapter 6.
7. Brennan, R.D. and Linebarger, R.N. (1964) A 11. Pack, A.I. and Murray-Smith, D.J. (1972)
survey of digital simulation: digital analogue Mathematical models and their applications in
simulator programs. Simulation, 3 (6),22-36. medicine, Scottish Medical Journal, 17, 401-409.
8. The Simulation Software Committee of 12. Dickinson, e.J. (1977) A Computer Model of
Simulation Councils Inc. (1967) The SCi con- Human Respiration, MTP Press, Lancaster.
tinuous system simulation language (CSSL).
Simulation,9 (6), 182-203.
PROBLEM ORGANIZATION FOR
CONTINUOUS SYSTEM SIMULATION 3
It has already been emphasized that the where yes) is the Laplace transform of yet),
choice of state variables is never unique. F(s) is the Laplace transform of f(t) and all
Hence the state variables in this approach initial conditions on yet) and its derivatives
may be numbered in any order. are zero. The quantity yes) can then be manip-
ulated algebraically to give
3.4 TRANSFER FUNCTION DESCRIPTIONS (Ms 2 + Rs + K)Y(s) = F(s) (3.11)
Linear models, expressed in terms of ordinary and thus
differential equations, can also be written
very conveniently as transfer functions, yes) 1
= (3.12)
defined as the ratio of the Laplace transform F(s) MS2 +Rs+K
of the model output variable to the Laplace
transform of the input when all initial condi- In general a transfer function G(s) may be
tions are zero. This type of description is par- written in the form
ticularly widely used in the analysis and
design of engineering control systems. An
yes) = G(s) = A(s) (3.13)
U(s) B(s)
understanding of Laplace transform methods
is an essential prerequisite to the use of trans- where yes) is the Laplace transform of the
fer functions, and a summary of some of the output variable, U(s) is the Laplace transform
relevant aspects of Laplace transforms may be of the input and where A(s) and B(s) are poly-
found in Appendix A. Most introductory nomials in s. In all cases the output transform
texts dealing with control systems engineer- is the transform of the input times the transfer
ing (e.g. refs 1 and 2) provide a detailed function.
account of transfer function models and their Many important properties of a transfer
derivation. An outline of the relevant theory function, and of the system which it repre-
is provided in Appendix B. It is assumed, in sents, depend upon the denominator B(s). The
the remainder of this chapter and in other sec- roots of the characteristic equation B(s) = 0
tions of this book, that the reader is familiar largely determine the form of the output
with Laplace transforms and with the funda- when the transfer function is subjected to a
mentals of transfer function descriptions. given input. These roots are the poles of the
The assumption that all initial conditions transfer function. The number of poles is, of
are zero using a transfer function type of course, equal to the 'order of the system as
Bond graph representations 27
discussed above in the context of state- the loss within the resistor. Equivalent com-
variable descriptions. ponents in mechanical and hydraulic systems
are translational or rotational dampers and
flow resistance respectively.
3.5 BOND GRAPH REPRESENTATIONS
'R
Fig. 3.2 Four generic one-port compne~s for a but it is not always possible to assign this pre-
bond graph: resistor (R), inertia (I), capacItor (C)
and source (S) components. ferred causality in practical applications.
--
R
assigned causality in order to indicate that the
variable which they store is taken from the
junction to which they are bonded. Causality at
Vr
~i C s~
i
·l
1
also shows the order in which equations are
to be solved. In the case of an effort or flow vel
source, causality is clearly implicit. Energy C
storage elements (C and I) have a preferred
causality involving effort output for C com- Fig. 3.4 An example of a common flow junction in
ponents and flow output for I components, a bond graph.
Bond graph representations 29
R
ir ic
e
-}
el R C S \. 0
el' C
Clearly a zero junction imposes the same A gyrator, on the other hand, couples domains
effort on each connected component, while a a and b according to an equation of the form
one junction imposes the same flow on each
(3.15)
component. Both types of junction conserve
power, and in a zero junction the flow vari- This use of coupling components may
ables must sum to zero algebraically. Con- be clarified by considering two specific exam-
versely, at a one junction, the effort variables ples. Take first the hydraulic cylinder and
must sum to zero algebraically. This means piston shown in Fig. 3.6. Here the piston force
that in electrical circuit applications, for F and velocity V can be related to the pres-
example, the one junction is a statement of sure P and volumetric flow rate f through the
Kirchhoff's voltage law, while a zero junc- equations
tion is a statement of Kirchhoff's current
law. F=PA (3.16)
f=VA (3.17)
F,v
a
~I F
V
' TF
p
~
: I
a
V )T,.o V
'GY
T ~
I .n
Fig. 3.6 Examples of bond graph coupling components: (a) a hydraulic cylinder and piston (a
transformer) and (b) a DC electric motor (a gyrator).
3.5.5 SIGNALS IN BOND GRAPHS in Fig. 3.8, the pressure in the second tank
would not influence the flow from the first
A signal is somewhat similar to a bond but
carries a variable, in one direction, from one tank to the second. Since the flow from the
first tank depends only on the pressure at the
point to another and can be amplified. The
output of the first tank a bond is not used and
difference between a bond and a signal may
a signal is employed instead.
be understood best by means of examples.
Consider a system involving two coupled
tanks, as shown in Fig. 3.7. In this system the
3.5.6 SIMULA nON PROGRAM
flow out of the left-hand tank depends on the
DEVELOPMENT FROM BOND GRAPH
pressure difference between the two tanks.
MODELS
The bond graph in this case involves bonds
but no signals. However, if the two tanks In principle it is possible to translate a bond
were connected in a different way, as shown graph model into a simulation program either
0
P1
l
\. 1
P2
f
~o
P1 ~-+L
f
P1 0 ~
l
" 1 ) 0
P
2
1f
manually or using software [7]. There are In the block diagram approach to the
continuous system simulation programs for description of a system each element is
which the dynamic model may be described described by a single block. The arrow enter-
in bond graph form. Examples include ing the block represents the input variable
CAMP-G [8] and ENPORT [9]. For example, and the arrow leaving the block represents
in applying the CAMP-G program the user the output variable. The block contains a
creates a bond graph on the screen and the mathematical expression, usually in transfer
software allows a complete causal analysis to function form, which relates the output to the
be carried out. A CAMP-G output file may input. Signal flow diagrams incorporate
then be created which can be used as an input exactly the same information as the block
file for ACSL, which is a widely used simula- diagram. In this case input and output vari-
tion tool discussed in later chapters. ables are represented by nodes and the line
connecting two nodes is the equivalent of the
block in the block diagram approach. A given
3.6 BLOCK DIAGRAM AND SIGNAL system can be represented at many different
FLOW GRAPH REPRESENT AnONS levels using a block diagram or signal flow
graph. For example Fig. 3.9 shows valid block
Complex systems can be described very conve- diagram and signal flow graph representa-
niently by means of block diagrams or signal tions for the mass-spring-damper system of
flow graphs. Both these forms of graphical rep- equation (3.1). In these diagrams the input
resentation are described in Appendix C. The variable is F(s), the Laplace transform of the
importance of block diagrams and signal flow force f(t), and the output is Y(s), the Laplace
graphs is very considerable for simulation transform of the displacement y(t). No other
since these provide a simple means of express- variables appear in the diagram, which
ing the structure of a complex model. Many describes only the relationship between the
modern simulation packages are block chosen input variable and the chosen output
diagram oriented and a thorough understand- variable. Figure 3.10 shows another possible
ing of such diagrams is therefore essential. level of representation for the same model. In
32 Problem organization for continuous system simulation
course, be noted that the diagrams of Figs 3.9
and 3.10 are exactly equivalent and Fig. 3.9
Frs) 1 Y!s)
1----+ may be derived from Fig. 3.10 using standard
rules for block diagram or signal flow graph
manipulation and reduction, as outlined in
Appendix C.
Frsl + 1 1 VIs)
M S
Fig. 3.10 Detailed block diagram and signal flow graph representations of the mass, spring and damper
system of Fig. 3.1.
Methods for transfer function simulation 33
3.7.1 THE DIRECT CONSTRUCTION the system which is being modeled and is
APPROACH defined from the following equation
Most transfer functions of practical im-
yes) = yes) . E(s) (3.22)
portance can be manipulated into a general U(s) E(s) U(s)
form involving a ratio of two polynomials in
s and a gain factor. Consider the transfer The transfer function of equation (3.21) may
function now be split into two parts as follows
Fig. 3.11 A direct construction form of signal flow graph for denominator terms in a linear nth-order
single-input single-output transfer function description.
U(s) 1 E(s) l 1
S K Yes)
,,
Fig. 3.12 The complete signal flow graph, obtained by the direct construction method, for a linear nth-
order single-input single-output transfer function model.
elements, in order from the right, the com- added to make this set of equations into a
plete set of state equations is as follows: complete state variable representation of the
given transfer function is to add a single alge-
braic equation relating the output y to the
state variables. From the block diagram it is
clear that this equation is
y = K(aox 1+ a1x2 + ... + am-1X n- m + x n- m+l ) (3.27)
Equations (3.26) and (3.27) together form a
(3.26) complete state-space description for the given
transfer function and therefore provide a very
This set of equations is similar in structure to simple basis for implementation of the given
equations (3.9) above. All that has to be transfer function within a simulation. The
Methods for transfer function simulation 35
only disadvantage of this direct construction Using partial fractions it is then possible to
approach is that the properties of the result- express the right-hand side of this equation as
ing state-space model can be highly sensitive a sum of terms each of which involves a
to small numerical inaccuracies in coefficient single pole. The resulting equation is
values, especially in the case of high-order
models. yes) al az a"
- = - - + - - + ... + - - (3.29)
U(s) S+/3l S+/32 s+/3"
Each of the terms on the right-hand side of the
3.7.2 THE PARALLEL CONSTRUCTION equation has the same form, and the corre-
APPROACH sponding signal flow graph involves a parallel
structure as shown in Fig. 3.13. Once again the
The parallel programming approach for the state variables are taken to be the outputs of
derivation of a set of state equations from a each integrator. The state equations describing
given transfer function, and thus for simula- the signal flow graph are then of the form
tion of that transfer function, is appropriate
when the given transfer function has a Xl = -/31X l + u(t)
denominator in factored form as shown
below
(3.28) (3.30)
U(s) 1 Y(s)
Fig. 3.13 A signal flow graph, obtained by the parallel construction method, for an nth-order single-
input single-output transfer function model.
36 Problem organization for continuous system simulation
An additional algebraic equation is then 1
needed to relate the output y to the n state
variables. This has the form
1
(3.33)
X n_, = -13 n_jX,,_j + (a n-j+1 - 13 n-j+1 )x"_j+1
it is possible to construct a submodel block
diagram or signal flow graph which repre-
sents a single factor of the complete transfer
function. This submodel diagram has the
form shown in Fig. 3.14. Since the order of the
denominator is always equal to or greater
X, =-/3,x, +(a 2 -/32)X 2 + ... +(a" -/3 ")x" +u (3.37)
than the order of the numerator of the trans-
fer function for physically realizable systems, As with the parallel construction approach
the number of denominator factors is always the iterative method has properties which are
equal to, or greater than, the number of superior to the direct construction method in
numerator factors. Therefore, at the end of the terms of coefficient sensitivity. A cascaded
process of grouping these factors in pairs, arrangement of first- or second-order transfer
there may be some additional terms involving functions tends to be more robust to coeffi-
Methods for transfer function simulation 37
and
E(s) 1
= (3.41)
That is
Fig. 3.15 A signal flow graph for a single pole for E(s) =U(s) - [6s- 1E(s) + l1s- 2E(s) + 6s-3E(s)]
use in the iterative construction method.
(3.42)
cient errors than the equivalent structure and
involving multiple feedback loops.
Y(s) = S-1 E(s) + 9s-2 E(s) + 20s-3 E(s) (3.43)
3.7.4 AN EXAMPLE OF BLOCK DIAGRAM Equations (3.42) and (3.43) may be expressed
CONSTRUCTION FROM A TRANSFER in block diagram form by a model involv-
FUNCTION ing three integrators connected in cascade
Consider the following transfer function: as shown in Fig. 3.17. Taking the outputs
of integrator blocks as state variables
G(s) = Y(s) = S2 + 9s + 20 (3.38) allows the following set of simultaneous
u(s) S3 +6s 2 + l1s+6 first-order differential equations to be
established:
G(s) = y(s) =
U(s) 1 yes)
Fig. 3.16 A signal flow graph, obtained by the iterative construction method, for the special case of a
single-input single-output transfer function model with an nth-order numerator and nth-order
denominator.
38 Problem organization for continuous system simulation
Fig. 3.17 A block diagram for the example of section 3.7.4 by the direct construction method.
+
, x2
-6
Y(s)
s
,
S
x3
Fig. 3.18 A block diagram for the example of section 3.7.4 by the parallel construction method.
U(s) + 1 1 1
S S 4 S
2 3
Fig. 3.19 A block diagram for the example of section 3.7.4 by the iterative construction method.
systems are distributed in space as well as where z(x, t) is the state of a point at position
time and the mathematical models are de- x and time t. The parameter v is the velocity
scribed by partial differential equations. of transmission.
In situations in which the system reduces Applying the unilateral Laplace transfor-
essentially to a single input variable and a mation to this equation gives, assuming zero
single output variable of interest at one point initial conditions, a new equation of the fol-
in space, as often applies in control systems lowing form
applications, it is possible to derive a transfer
function type of description which involves aZ(x, s)
sZ( x,s ) = v-----'-- (3.56)
a distributed parameter. Such distributed para- ax
meter transfer function models may then be
where
used in the same way as conventional transfer
functions for lumped parameter models. The
two types of distributed parameter element Z(x, s) = f o
exp(-st)z(x, t)dt (3.57)
considered in this section are the pure time
delay and the distributed time delay. The solution of this equation is
az
-=v-
az (3.55)
Thus the transfer function relating the output
at ax transform yes) to the input transform U(s) is
Modeling of distributed parameter elements 41
This allows one of the terms in the above
Y(s)
--=exp (X2
-s-- XI-) (3.61) equation to be eliminated, giving
u(s) v
That is
-
Y(s)
= exp( -sT) (3.62)
Z(X, s) = A exp[ - [xl (3.67)
-~X2l
with an attenuation which varies with fre-
quency. Examples include the conduction of Y(s) = Z(x 2I s) = Aex p[ (3.69)
heat and the transmission of electrical signals
in a medium which is not loss free. The Hence, the transfer function relating Y(s) to
partial differential equation describing this U(s) has the form
situation has the form
(3.70)
(3.63)
----- SYSTEM - - --
INPUTS - -- -- MODEl - - --
OUTPUTS
I I
I I
I
I I
I I
L..,..
-+
------- -----
- - - - - -- CO-SYSTEM
- - --
SENSITIVI TY
FUNCTIONS
Fig. 3.20 A schematic diagram illustrating the generation of sensitivity coefficients using the co-system
approach.
equations with constant coefficients is facili- proved by Kokotovic [16] for multiloop struc-
tated by the application of a technique, devel- tures, such as those arising in the direct con-
oped first by Kokotovic [16], which is known struction method for the simulation of transfer
as the 'sensitivity points' method. functions (section 3.7.1), that the co-system not
A single-input single-output linear system only has the same structural form as the
may be described by the equation model itself but the sensitivity coefficients are
signals appearing at certain predetermined
yes) = G(s) (3.78) points in the co-system block diagram.
U(s) A general direct construction method block
where yes) is the Laplace transform of the diagram is shown in Fig. 3.21 for an nth-order
system output variable y(t) and U(s) is the system. Let this system be described by a
transform of the input u(t). The sensitivity of transfer function
the output to changes of any parameter q of
the model is then given by
1 1
U(s) K 5 S VIs)
,
I
Fig. 3.21 A signal flow graph for an nth-order transfer function (direct construction).
U(s) K
1 1
5 5 Yes)
Fig. 3.22 Sensitivity co-system structures for denominator and numerator coefficients of the nth-order
linear system of Fig. 3.21. Note that the sensitivity functions are the variables at the nodes of the
co-system.
",---ti
"2---tI
1---- V
V1
Vii! Of
dV1
l~1
aq
Y1
Vii! Clf
bVe
Vr h
iii
~ve •
~
~Vr _ _ _ _ _ _--.l
aq
Fig. 3.23 A block diagram illustrating the generation of parameter sensitivity functions for a general non-
linear element.
Fig. 3.24 A schematic diagram of a closed-loop control system involving a digital processor together with
analog to digital and digital to analog converters (ADC and DAC respectively) within the feedback loop.
units be replaced by appropriate idealized digital control systems can be found in any
elements. The process of analog to digital con- text on sampled-data theory and digital
version is a process which involves sampling control (e.g. refs 19 and 20). It should be noted
of the continuous variable which is input to that the analog to digital and digital to analog
the converter. This sampling process is con- converters are normally assumed to operate
ventionally represented by an ideal switch in a synchronous fashion with input and
which closes for a short period of time at output operations at the control computer
regular intervals. The time between samples taking place at the same time. This, and other,
is the sampling period of the analog to digital simplifying assumptions, which are widely
converter and the sampled value is assumed used in analytical treatments of digital control
to provide the discrete input to the control system modelling, may be avoided if thought
algorithm, which is implemented on the necessary in a simulation-based approach.
processor in terms of a software program. The The control algorithm within the digital
digital to analog converter converts the dis- processor is normally represented by a differ-
crete numerical result obtained from the ence equation. This provides a means of
control algorithm into an equivalent analog expressing the output of the digital controller
variable. It is conventionally represented in an at some time t = kT in terms of the input as
idealized model by a zero-order hold, which provided by the analog to digital converter at
gives an output which changes periodically that time, and by the output and input at pre-
but maintains a constant output level between vious time instants t = (k - l)T, t =(k - 2)T, ... ,
each output event. Figure 3.25 shows the etc. Here T is the sampling period and k is an
block diagram of an idealized model of this index defining the sample number being con-
kind. Further details of the modeling of sidered. For example, a very simple algorithm
- IDEAL
SAMPLER
Fig. 3.25 An idealized block diagram model of the digital control system of Fig. 3.24 using standard
sampled-data conventions and symbols.
References 49
which provides proportional control in which 8. Granda, J. (1990) CAMP-G Users Manual,
the error signal l(s) sampled by the analog to Mitchell & Gauthier Associates, Concord, MA.
digital converter and simply multiplied by a 9. Rosenberg, RC. (1990) ENPORT Reference
Manual, Rosencode Associates, Lancing, MI.
gain factor G to produce the output that pro- 10. Doebelin, E.O. (1980) System Modelling and
vides the continuous signal O(s) from the Response, John Wiley, New York, pp. 193-201.
digital to analog converter would involve a 11. Spriet, J.A. and Vansteenkiste, G.c. (1982)
difference equation of the form Computer-aided Modelling and Simulation,
Academic Press, London.
O(kD =G1(kD 12. Bryce, G.w., Foord, T.R, Murray-Smith, D.J.
If the control algorithm was now modified to and Agnew, P.w. (1976) Hybrid simulation of
water-turbine governors. Simulation Councils
involve integral contrbl, in which the error Proceedings, 6 (Part 1), 35-44.
signal is integrated numerically before being 13. Nelms, RM., Sheble, G.B., Newton, S.R and
returned in analog form through the digital to Grigsby, L.L. (1990) Simulation of trans-
analog converter, the algorithm would have a mission line transients with a distributed-
more complex form and would involve previ- parameter model using a personal computer.
ous values of the discrete variables 0 and I as Simulation, 55 (2), 103-107.
well as current values. Further consideration 14. Karnavas, W.J., Sanchez, P.J. and Bahill, A.T.
(1993) Sensitivity analyses of continuous and
is given to sampled-data system models in
discrete systems in the time and frequency
Chapter 5. domains. IEEE Transactions Systems, Man and
Cybernetics, 23 (2),488-501.
REFERENCES 15. Tomovic, R (1963) Sensitivity Analysis of
Dynamic Systems, McGraw-Hill, New York.
1. Palm, W.J. (1986) Control Systems Engineering, 16. Kokotovic, P.V. (1964) Method of sensitivity
John Wiley, New York. points in the investigation and optimisation of
2. Golten, J. and Verwer, A. (1991) Control System linear control systems. Automation and Remote
Design and Simulation, McGraw-Hill, London. Control, 25 (12), 1512-18.
3. Paynter, H.M. (1961) Analysis and Design of 17. Wilkie, D.F. and Perkins, W.R (1969)
Engineering Systems, MIT Press, Cambridge, Generation of sensitivity functions for linear
MA. systems using low order models. Transactions
4. Karnopp, D.C., Margolis, D.L. and Rosenberg, oftheIEEE, AC-14 (2), 123-30.
RC. (1990) System Dynamics: A Unified 18. Vuskovic, M. and c':iric, V. (1966) Structural
Approach, 2nd edn, John Wiley, New York. rules for the determination of sensitivity
5. Rosenberg, RC. and Karnopp, D.C. (1983) functions on nonlinear non stationary
Introduction to Physical System Dynamics, systems, Sensitivity Methods in Control Theory
McGraw-Hill, New York. (ed. L. Radanovic), Pergamon, Oxford,
6. Thoma, J.U. (1990) Simulation by Bond Graphs, pp.154-65.
Springer, Berlin. 19. Isermann, R. (1989) Digital Control Systems,
7. Karnopp, D. (1984) Direct programming of 2nd edn, Springer, Berlin.
continuous system simulation languages 20. Franklin, G. and Powell, J.D. (1980) Digital
using bond graph causality, Transactions of the Control of Dynamic Systems, Addison-Wesley,
Society for Computer Simulation, 1 (1),49--60. Reading, MA.
THE PRINCIPLES OF NUMERICAL
MODELING 4
The solution of a differential equation can be This equation suggests that evaluation of the
obtained by mathematical operations which
term involving the integral is only possible if
involve integration. In order to obtain a
x(to) and x(to + h) are both known. Numerical
solution using a digital computer the analytic integration techniques provide approx-
processes of integration must be replaced by a imations which allow this difficulty to be
numerical method which can yield an overcome and permit the next value of the
approximation to the true solution. A con- variable x[x(to + h)] to be determined using
tinuous variable x which is a function of the
only the current value of that variable [x(to)]:
independent variable t (i.e. time) is rep~ Many different integration methods eXIst
sen ted within a digital simulation by a senes which are suitable for use in the solution of
of numbers, or samples, Xo' Xli XU ... , X". These ordinary differential equations and thus for
samples define the variable X in terms both of
continuous system simulation. Some of these
magnitude and sign at the times to, tl , t2 ,···, tn· are very simple in terms of the numerical
For many purposes it may be assumed that methods used and involve relatively crude
these samples are equally spaced in time and
and simple approximations. Other methods
it is clear that if the sampling rate is high the
are potentially much more accurate in ter~s
information loss owing to sampling can be of the numerical methods used but thIS
small. However, the step size, h, does affect
additional numerical accuracy is achieved at
the size of the error in the numerical some additional cost.
approximation and the choice of h must
One of the penalties associated with some
depend upon the dynamic characteristics of of the more complex integration methods is
the system under consideration. an increased tendency to display numerical
The solution of a differential equation instability. This means that the numerical
involves integration of a variable which can algorithm itself tends to become unstable for
involve a function of that same variable.
a chosen integration step size although the
Consider the first-order differential equation
underlying mathematical model is completely
stable. Even if a technique is stable for a
dx (t) = ax(t) + bu(t) (4.1) chosen integration step size there may well be
dt
problems in terms of the transient behavior of
Here the quantity on the left-hand side of the the numerical algorithm which outweigh the
equation, which on integration gives x(t), ~s advantages of high steady-state accuracy. A
itself a function of x(t). If the value of X IS second, and perhaps more obvious, penalty
known at time to then the value at time to + h is incurred in using complicated integration
given by algorithms is concerned with computation
52 The principles of numerical modeling
time. Complicated algorithms, on a conven- this approach it is assumed that the quantity
tional serial computer architecture, require to be integrated, the derivative function,
more processor time for each integration step. remains constant over the integration inter-
A simple algorithm, although less accurate, val. In other words, the value of the
may well be necessary in applications which derivative throughout the interval h in
are time critical and involve relatively fast equation (4.1) is the value of the derivative at
dynamics. time to. In terms of equation (4.2) this gives
A number of different factors must be taken
into account in the selection of an integration x(t + h) = x(t) + h( ~) (4.3)
method for the solution of a given problem.
One obvious question to be considered is and this shows clearly that the new value of x
whether the integration steplength should at time (t + h) is calculated from the value of x
be constant or variable. Some integration at time t and the derivative at time t. For the
algorithms select their own integration step case of equation (4.1) this means that for a
size since a value of h which is satisfactory in given input variable u(t) the solution of the
one part of the time interval may not be good differential equation starting from an initial
elsewhere. For other types of application, value x(to) follows the sequence shown below:
such as real-time simulation, there may be
good reasons for avoiding such methods.
One broad classification for integration
techniques used in simulation involves the x(to + 2h) = x(to + h)
timing of the sampling operation which is +h[ax(to + h) + bu(to + h)]
inherent in the integration process. On the one
hand, an algorithm which requires only x(to) x(to + 3h) = x(to + 2h)
and u(to) to calculate the right-hand side of +h[ax(to + 2h) + bu(to + 2h)]
equation (4.2) and thus provide an estimate of
x(to + h) is known as an 'explicit' method. On etc.
the other hand, a technique in which x(to + h) The processes involved can be illustrated by
itself appears on the right-hand side of the taking a particular example and following
expression to be evaluated inevitably presents through the calculations. Consider the case
more difficulties and is termed an 'implicit' where a =-2 and b =2 and the input is a unit
approach. Implicit algorithms are not nor- step function, u(t) = 1, with an initial value
mally used in continuous system simulation. x(to) = o. The corresponding transfer function is
One example of an implicit method is the
backward Euler integration algorithm in X(s) 2 1
--=--=--- (4.4)
which the quantity x(t + h) is approximated U(s) s+2 1+0.5s
using values of x(t) and dx(t + h) / dt, where h
For the unit step function input it may be
is the integration interval.
shown, by means of the final value theorem,
that the theoretical steady-state response of
this system is 1.0. It is also evident from the
4.2 FIXED-STEP INTEGRA nON
form of the transfer function, as shown in
ALGORITHMS: ONE-STEP METHODS
equation (4.4), that the system has a time
constant of 0.5 s and represents a system with
4.2.1 INTRODUCTION
simple first-order lag characteristics. An
One simple integration algorithm which example of such a system is the simple
provides a useful introduction to fixed-step resistor-capacitor network of Fig. 4.1 with the
integration methods is the Euler method. In state variable x(t) chosen to represent the
Fixed-step integration algorithms: one-step methods 53
0.9
,...
~o.s
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0
Fig. 4.2 Results from a simulation run involving a unit step input voltage change for the model defined
by equation (4.4) using Euler integration method. (a) Results obtained using an integration steplength of
0.1 s over a period of 5 s. (b) Results for an integration steplength of 0.05 s over a period of 2.5 s. The exact
responses, obtained analytically using Laplace transform methods, are shown by continuous curves in
both cases, whereas the simulated results are shown by means of dashed lines.
1.6 a
1.4
-.e
)(
,
1.2
I
:, .•
, ,"
I
,,"
,, \ ,:'
1.0
, ,, ",
",
0.8
0.6
0.4
0.2
0.0
0 4 6 10 12 14 16
time t
2.6
b
2.4 ,
i •
) ••
~
2.2
fi " I:"
ii
~ ~\ "
I, "
"""
j ,I
.e
)(
2.0
I
:"~
"
"
",I""'I iI,"' "
I:
:\" i"\ :1I1\: ,n, !!
1.8
:\ ,I
I,
"
',
' I
:\ , : 'I
:, I
I\,
I,
,I ', I I,
'' !\ ,, ,, ,," ,I II , ,I
I,
1.6 ,' II
,, ' I
, \ ,I I
:, 1 ,
I ' I
I,
,, , ,i , ,,, ,,, , :,,
' I
; , I
I
I I
,, ,
I
1.4 , I I !, ,I
:, , :
I
, I
,
,i i ,,,
I
: \ ,i I
,,:
I
I
I I
, I
, ! ,,, ,, :, ,:
I I
1.2 , \ I I
,1 :, I :,,
I , :
I I
I :I
1.0 ,I , I I
, I
, I ,,
I I I I
,I : ,, I ,
, :
1,, :, i:
I
0.8
,,i ,, I I
, 1 ,, I
,, I
,
,,
0.6
,I '' ,, :I ,, ,,
! I I
,'
I'
.::, .., . ,,, .''
, I I
I
,, ,, I
\ ,, I , ,
I
I
i
0.4 I
I'
""
,,
, , ,, ,,' , ' ,, I I, I'
I' ,,
" "
I'
" " ,, ',
, I
I,
I'
,: I
H
,
0.2 I'
I' ' I
" . "'
I
""•
" " ,I
'().2
: I
"I'
~
I
~,
I
".,
"
.", "
"" ":i
""
"
"
I,
"I,
I,
I,
,
I
,
I
~
I
I "
~
I
I
,,
'().4
I
,
'().6
0 2 4 6 8 10 12 14 16 18 20
time t
Fig. 4.3 Results similar to those of Fig. 4.2 obtained using the Euler integration method with integration
steplengths of 0.8 s (a) and 1.01 s (b). Dashed lines show the simulated responses in both cases, while
exact solutions (calculated using Laplace transform methods) at the time instants corresponding to the
start of each integration interval are shown by the continuous lines.
56 The principles of numerical modeling
If h is small and if x(t) and its derivatives are value is not known. The only known quan-
all known, the truncated series involving tities are the approximate value and an
(n + 1) terms can be used to calculate x(t + h). estimate of the error or a bound on the
When only terms up to and including maximum size of the error. In the case of
h"xi"l(t)/n! are present in the series the numbers dose to unity the absolute and
representation is said to be a Taylor series of relative error values are very similar. For
order n. numbers which are not close to I, the absolute
Examination of the Taylor series in and relative errors can be very different. For
equation (4.8) shows that Euler's method example, for a true value of 345 678 and an
involves a Taylor series of order 1. Clearly the approximate value of 345 123 the absolute
inclusion of additional terms in the Taylor error is 555 but the relative error is only
series is one way of improving the accuracy 0.0016 (0.16%). Similarly, if the true value is
of the numerical integration process. This is, 0.00003 and the approximate value is 0.00002
however, of limited practical value for the absolute error is 0.00001. However, the
general simulation use since we need to relative error in this case is 0.5 (50%).
obtain the derivatives associated with the Errors in numerical calculations are of three
terms in h, h2, h3 , ••• , etc. It is clearly not types, which are known as inherent errors,
appropriate to carry out this differentiation truncation errors and round-off errors.
analytically in the general case, and esti- Inherent errors are errors in the data upon
mation of the higher derivatives can present which the program depends. In the context of
major difficulties. simulation these might be errors in the values
Although Taylor series solutions are not of constants describing elements within
themselves of great practical value in the mathematical model upon which the
providing algorithms which are of use for simulation program is based. Errors of this
simulation, they do provide a basis for the kind may arise because of simple mistakes in
evaluation and comparison of other methods recording values to be used, because of errors
that are of practical importance. In this in the measurement of some physical quan-
respect the Taylor series provides a tity within the system represented by the
convenient yardstick: some integration model or because a number, such as TT, has
methods will agree with the Taylor series not been represented with sufficient accuracy.
solution only as far as terms in h, others Errors also arise in part because of the fact
through terms involving h2 or beyond. that the integration process involves an
approximation which can be regarded as
equivalent to truncation of a Taylor series
4.2.3 ERRORS IN FIXED-STEP INTEGRAnON
after some finite number of terms. In general,
METHODS
the truncation error at each step is of the
An error is defined here to be the true value order of the first term dropped. For example,
minus the approximate value and may be in the Euler method the truncation error can
expressed either as an absolute error or as a be shown to be of the order of h2 • In more
relative error. The absolute error is defined to complex methods, such as the fourth-order
be the difference between the true value and Runge-Kutta approach discussed in the next
the approximation. The relative error is the section, the truncation error is of the order of
absolute error divided by the approximate h5 or better. For small values of the integration
value. It might appear at first to be better to steplength, there are clearly very large
define the relative error as the absolute error benefits in terms of accuracy in using these
divided by the true value, but it must be more complex techniques, in comparison
remembered that in most situations the true with the Euler method.
Fixed-step integration algorithms: one-step methods 57
The third type of error, round-off error, Euler's method can be extended in a
arises because quantities are represented number of ways. One approach of this kind is
within the computer using a finite number of known as the 'improved Euler method'. This
binary digits. Different compilers handle the is a two-stage method in which an estimate of
rounding of numbers in different ways, and the next point is first made using the Euler
detailed discussion of this topic is beyond the method. The derivative is calculated for
scope of this text. Many books on numerical this estimated point. The average of this
methods provide further information on derivative value and the derivative at the
round-off problems and truncation errors beginning of the step is then used to
(e.g. ref. 1). recalculate the next point. The process can be
illustrated using the example previously
discussed in the context of Euler's method
4.2.4 EXPLICIT RUNGE-KUTTA METHODS
involving the differential equation
Explicit Runge-Kutta methods are closely
dx
related to the Taylor series approach and - = -2x(t) + 2u(t) (4.10)
provide a convenient solution to the problem dt
of determining the necessary derivatives For a unit step input [u(t) = 1 for t > 0)] and an
without using any form of analytic integration interval of 0.1 s the first estimate
differentiation. These methods have three of x(t + 0.1) is given by
main properties which may be summarized
xe(t + 0.1) = x(t) + 0.1[2 - 2x(t)] (4.11)
as follows:
that is
1. Explicit Runge-Kutta methods are one-
step methods and require information to xe (t + 0.1) = 0.8x(t) + 0.2 (4.12)
be available from one preceding point
Using this value for the estimate of x at time
only.
t + 0.1 and substituting in equation (4.10) gives
2. Runge-Kutta methods do not require the
an estimate of the derivative value at this time
evaluation of any derivatives of the
function I in the differential equation
( dX)
dt
= 2 - 2[0.8x(t) + 0.2]
dx e
- = I(x,u,t) (4.9) =1.6-1.6x(t) (4.13)
dt
All that is needed is the value of the The average derivative value is then given by
function I(x, u, t) itself.
3. Runge-Kutta methods include a number
of different algorithms which are dis-
o.s[(~) + (~)J = 1.8 -1.8x(t) (4.14)
1.0
0.9
0.8
-
~
)( 0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0
M ~ M 0.6 M W U U U U w
time t
Fig. 4.4 Results similar to those of Figs 4.2 and 4.3 using the 'improved Euler' integration method with an
integration step of 0.1 s. The corresponding analytic solution is shown for comparison. The dashed line
shows the simulated response while the continuous line is the analytical result.
It should be noted that the calculations second-order case. Let the subscript k denote
shown here for the Euler method and the evaluation at time t k• Then tk+! denotes time
improved Euler method are particularly (t k + h) where h is the integration st~plengh.
straightforward because the input u(t) has Truncating the Taylor series after the second
been chosen to be a constant. In the more term gives
general case where the input is time varying,
the equations to be solved repeatedly at each Xk+l = Xk +hXk +0.5h2xk (4.16)
integration step are slightly more complicated. Thus, in terms of the notation of equation
However, the example chosen does illustrate (4.9), that is
the two methods very effectively.
Xk+l = Xk + h/k + 0.5h 2 (df )
The improved Euler method can be shown (4.17)
to involve an algorithm which agrees with the dt k
/32 = 0.15875964
/33 =-3.05096516
Y3 = 3.83286476
dx b
-=U (4.25)
dt Fig. 4.5 The difference between the Euler (a) and
trapezoidal (b) integration algorithms.
the trapezoidal integration method gives, for
the single step from tk to tk+!
4.3 FIXED-STEP INTEGRAnON
ALGORITHMS: MULnSTEP METHODS
f
tk+l
I
\ I
\ I
\ I TRUNCATION
4.4 PROBLEMS OF INTEGRATION \ I
" TOTAL / /
STEP-SIZE SELECTION AND NUMERICAL
" //
INSTABILITY
' ........... _----""'//
The process of choosing an appropriate
integration step size should begin by careful
consideration of the dynamics of the
mathematical model being simulated. A trial
ROUND-OFF
value, which is small in comparison with time
constants in the model or the period of any STEP SIZE
oscillatory components expected in the model
response, must be selected. A test solution can
then be generated from the simulation Fig. 4.6 An illustration of the dependence of round-
program using this initial choice of integration off and truncation errors on integration step size.
62 The principles of numerical modeling
step size usually involves an initial choice 4.5 VARIABLE-STEP INTEGRATION
which is larger than ideal. Successive ALGORITHMS
iterations involving smaller integration steps
In conventional fixed-step methods there is
generally move one towards the optimum.
always a need to estimate the precision of the
The use of double-precision quantities
integration calculations. A tolerance is
can often be helpful in reducing the effects
usually defined and the integration step size
of round-off errors within simulation
must be chosen to ensure that the errors do
programs.
not exceed the level specified by that figure.
In addition to the problems associated with
The determination of a suitable integration
numerical errors, which have been outlined
step size can present considerable difficulties
above, it is found that limits exist in terms of
in some cases. In particular, in nonlinear
the integration intervals for which a given
problems it may be appropriate to use a small
integration algorithm will be stable. Numer-
step size for one part of the operating range
ical instability can occur quite independently
while a much larger step may be sufficient
from questions of inaccuracy because of
elsewhere. Methods have therefore been
truncation or round-off and depends on the
developed for the automatic adjustment of
form of algorithm chosen and the step size
integration step size.
adopted. In simple fixed-step methods there
One widely used algorithm for variable
are a number of simple 'rules of thumb'
step integration is the Fehlberg algorithm.
which are quite widely used to assess the
The Fehlberg method of order 4 involves an
upper limit of the integration interval. For
algorithm of the Runge-Kutta type. The
example, if instability occurs at a particular
variable-step feature is incorporated by
value of integration step size, he' the critical
determining the difference between the
sampling frequency for the simulation, fe' is values of x(t k + h) calculated by fourth-order
given by the relationship
and fifth-order procedures. This difference is
taken as an estimate of the truncation error,
(4.27) and if this exceeds a given tolerance the step
size is halved. If the difference is much less
than the tolerance the step size is doubled.
The ratio of this critical sampling frequency
Many variable-steplength algorithms are to
to the highest frequency of interest in the
be found in the literature, and a number have
model, fmax, is of special interest. The ratio
been incorporated within widely used simu-
must always be greater than 2 and would
lation packages. A useful review of variable-
generally be chosen to be less than 10. In the
step integration techniques is provided in a
case of complex integration algorithms the
paper by Crosbie and Hay [5].
value chosen is likely to be in the upper half
of the range. One obvious question is how
the value fmax is found and this can clearly
4.6 PROBLEMS OF 'STIFF' MODELS
be difficult in practical problems. One
approach, which is useful in dealing with Systems in which responses contain some
models which are linear or can be linearized components which are very fast compared
for particular operating conditions, is to use with others are called 'stiff' systems. In the
transfer function methods and assess the case of linear models this corresponds to a
natural frequencies of the model and the situation where the characteristic equation
frequencies associated with each simple has roots which cover a very wide range of
pole. The highest of these frequencies is then values. A range of 1000 to 1 for these roots
taken as fmax' certainly constitutes a stiff system.
Discontinuities and their effects 63
The difficulties encountered in the both the very fast dynamics and the very
simulation of stiff systems arise because of the slow dynamics in a single simulation
nature of the transient response of a system of experiment is essential. One can usually make
this kind. Since the response contains some appropriate simplifying assumptions, which
components which are very much faster than can allow the stiffness of the model to be
others, the integration step size must be small reduced significantly, without adversely
enough to suit the most rapidly changing affecting the accuracy of the results. Bond
component. However, this small step size may graphs have been found to provide a basis on
be unacceptable in terms of the computing which to consider alternative models, and
time needed to obtain the complete solution these techniques can therefore sometimes be
which shows the full-time history of the slow of value in the elimination of irrelevant parts
components of the response. Even some of a model [6].
variable-step algorithms may not be capable Situations in which the complications of
of reducing the computing time to an stiffness cannot be avoided include certain
acceptable value in such cases. types of nonlinear problem and models
One obvious solution to the problem of involving time-varying coefficients. In such
stiffness is to redefine the model since this cases use can be made of special integration
type of numerical problem is often avoidable algorithms which have been developed for
by the use of appropriate approximations. stiff systems. A number of such algorithms
Physical insight and practical engineering are available, including one by Gear which is
knowledge can often be used at the initial widely used and is included as an option
model development stage to eliminate in many continuous system simulation
unnecessary features which can lead to languages [7].
problems of stiffness. For example, in the
simulation of a chemical process the
dominant time constants may well be of the 4.7 DISCONTINUITIES AND THEIR
order of many minutes. If a particular EFFECTS
instrument acting as a sensor for some
As has been seen in earlier sections of this
variable within this system has dynamic
book it is conventional to describe a con-
characteristics which can be described by a
tinuous system model by an equation of the
transfer function having a time constant of
form
0.1 s, it may be appropriate to consider
modeling the sensor simply as a gain constant dx(t)
and eliminating the sensor time constant, cit = f[t,x(t),u(t)] (4.28)
which is negligible in comparison with the
other process time constants. Such a course of where x(t) denotes the system state vector (of
action would be particularly appropriate if dimension n), f is the vector of derivative
interest is focused on the long-term behavior functions and u(t) is the input. Discontinuities
of the system. On the other hand, if one was can arise in two ways. The first of these
interested in examining the short-term behav- involves situations in which there is more
ior of the sensor following a sudden change than one set of derivative functions for the
in the variable being sensed, it would be more system model. Switching can occur from one
appropriate to regard the time constants of derivative function to another when a
the process as infinite and consider only the particular state variable in the model reaches
instrument dynamics in the simulation study. a predefined threshold value. Such transitions
There are relatively few situations with a may cause considerable problems for numer-
model of this kind where the inclusion of ical integration routines. The second way in
64 The principles of numerical modeling
which a discontinuity can occur is when there The zeros of the switch functions cPJt, x(t)]
is an instantaneous change in one or more of define events which involve discontinuous
the problem state variables. behavior.
The distinction between these two types of The same type of approach, based on the
situation may perhaps be understood more location of the zeros of a switch function, may
completely by considering some examples of be used also for the case of an instantaneous
simple cases in which discontinuities can change in a state variable. The first problem
occur. An example of the first type, in which is essentially one of establishing whether or
the derivative suddenly changes value, can not an event has occurred within the most
arise when simulating a system in which the recent integration step. If an event has
input u(t) is discontinuous, such as a square occurred the second stage is to determine the
waveform. A small change in the independent exact time of that event. This two-stage
variable (which we are assuming to be time) process can be described as a 'detection
can cause a very large change in the value of problem' and a 'location problem'.
the derivative owing to the sudden change in The detection problem can be approached
u(t) in going from one level to another. The by making use of information about the switch
resulting response in terms of time history of function cP; and its derivatives evaluated at the
x(t) can be highly sensitive to small errors in start and end of the integration step. If the start
timing of such sudden transitions in the of the current integration step is at time t1 and
derivatives. The second type of discontinuity the end of the same step is at time t21 it is
is typified by the problem of simulating an possible to define two auxiliary quantities 51
electrical circuit that involves diodes. These and 521 the signs of which can be used to
devices can conduct when operated with a establish the presence or absence of an event
positive voltage bias, but become non- in most common situations. The technique,
conducting when the voltage across the which has been used by Carver [8], and others,
device is negative. Since currents and voltages involves the following two equations:
are inevitably state variables in the
51 = cPi (t 1 )cPi (t 2 ) (4.30)
mathematical description for an electrical
system such as this, the discontinuities
5 = dcPi (t ) dcPi (t ) (4.31)
involve sudden cl}anges in quantities which 2 dt I dt 2
are state variables of the model. Similar
situations arise in mechanical systems which It can be shown that three main conditions
involve hydraulic elements such as valves. arise. These are as follows:
The simulation of the impact of snooker balls
1. If 51 > 0 and 52> 0 it may be assumed that
provides a further good example of a system
no critical points lie in the integration step
involving discontinuous behavior in terms of
under consideration and the integration
the state variables.
routine may be allowed to continue.
2. If 51 < 0 and S2 > 0 it should be assumed
4.7.1 NUMERICAL METHODS FOR that a single event lies within the interval.
DISCONTINUOUS PROBLEMS 3. If S1 > 0 and S2 < 0 it should be assumed
that there are two events within the inte-
The condition governing the switch from one
gration step in question and the step size
set of derivative functions to another may
should be reduced in order to modify the
often be described in terms of the solution of
situation and get a single event within the
a set of algebraic equations
interval. These three cases are repre-
eMt, x(t)] = a i = l,2,,,.,m (4.29) sented graphically in Fig. 4.7.
Discontinuities and their effects 65
One situation which is not adequately ~i
covered by the conditions given above is
shown in Fig. 4.8. Two zeros occur in the
integration interval, but these are not
detected because the values of <Pi and its first
derivatives at ends of the interval satisfy
condition 1. It should also be noted that in
some cases, such as that shown in Fig. 4.9, a
condition 51 < a and 52 < a could also define an
event within the interval.
til
tz
·1
G(s) PS)
a
-/.-j G,,(s)
H G(s)
~-
b
Fig. 5.1 A linear continuous transfer function model and the equivalent sampled-data approximation.
(a) The block diagram representation of a linear single-input single-output system having transfer function
G(s). (b) The conventional block diagram representation of the corresponding sampled-data approximation
involving a sampler and hold block Gj,(s). This provides the basis for a digital simulation of G(s).
(5.12)
G(S)=~ (5.7)
s+ 10 The response for a unit step function input
applied at time zero can be calculated very
For the case in which a zero-order hold is easily by hand, assuming zero initial condi-
used, the overall transfer function is tions. At the first evaluation (t = 0.0) the
quantities y[(k-1)T] and u[(k-l)T] are zero
Y(z) _ z -1
u(z) z
z[ 10]
s(s + 10)
(5.8) so that y(kD found from the difference equa-
tion is also zero. At the second evaluation
Linear systems modeling using sampled data 69
u[(k-l)T] is unity while y[(k-l)T] is the calculated response from the discrete approx-
value of y(kT) found at the first evaluation imation would not be exact and the error
(i.e. zero). This gives a new value of y(kT) of would be dependent upon the sampling
0.1813, and this is therefore the response at period.
t = 0.02 s. This value is then assigned to the The programming of a difference equation
quantity y[(k -l)T] for evaluation of the model using a general-purpose high-level
response at 0.04 s, and so on. Table 5.1 shows programming language, such as Fortran,
the step response of the approximate model Pascal or C is very straightforward. The solu-
calculated in this way, together with the tion time can be made very short and this
exact response of the continuous model type of approach has proved valuable for
found by means of Laplace transforms. It can real-time applications.
be seen that in this case the agreement is
almost exact, the only differences being asso-
ciated with round-off effects in the fourth
decimal place. For other forms of input the 5.3 MODELING OF LINEAR SYSTEMS
USING SAMPLED-DATA APPROXIMATIONS
Table 5.1 Results for the simulation of the first- 5.3.1 CASCADED ELEMENTS
order transfer function of equation (5.7) for a unit
step function input. The table shows the numerical The modeling of cascaded linear systems by
results y(kT), for T = 0.02 s, together with the means of sampled-data approximations
response of the continuous model, y(t), found involves a straightforward extension of the
using Laplace transform methods approach described in section 5.2 for a single
transfer function. For example, the discrete
k kT y(kT) y(t)
model of a system involving several blocks in
1 0.02 0.1813 0.1813 cascade (as shown in Fig. 5.2a) can be
2 0.04 0.3297 0.3297 obtained by attaching a sampler and data
3 0.06 0.4513 0.4512 reconstructor to the input and a sampler at
4 0.08 0.5507 0.5507 the output of each element. The equation
5 0.10 0.6322 0.6321
defining the output of the first element Gj(s)
6 0.12 0.6989 0.6988
7 0.14 0.7535 0.7534
then becomes, in the discrete model
8 0.16 0.7982 0.7981
9 0.18 0.8348 0.8347
10 0.20 0.8647 0.8647
11 0.22 0.8892 0.8892
12 0.24 0.9093 0.9093
Similarly, for the outputs of the other ele-
13 0.26 0.9258 0.9257
14 0.28 0.9392 0.9392 ments we have
15 0.30 0.9502 0.9502
16 0.32 0.9593 0.9592
17 0.34 0.9666 0.9666
18 0.36 0.9727 0.9727
19 0.38 0.9776 0.9776
20 0.40 0.9817 0.9817
21 0.42 0.9850 0.9850
22 0.44 0.9877 0.9877
23 0.46 0.9900 0.9899
24 0.48 0.9918 0.9918
25 0.50 0.9933 0.9933
(5.14)
70 Sampled-data models and operator methods
-1 GII(S)
Uls)
T
~
Fig. 5.2 Block diagrams of a linear system involving n single-input single-output cascaded blocks. (a) The
continuous model. (b) The sampled-data equivalent.
R(s) Yes)
G(s)
R(s) R(z)
T
- U(z)
-
Y(z)
T
V(z)
b
Fig. 5.3 Block diagrams of a linear system involving feedback. (a) The continuous model. (b) The corre-
sponding sampled-data representation.
The discrete transfer function for the linear y(kT) = f[Yl (kT)] (5.19)
portion of the system is in this case given by
the equation If, on the other hand, the nonlinearity pre-
cedes the linear element, as shown in Fig. 5.5a
Y1(Z) = Z[Gh(s)G(s)] (5.17) and b, the discrete output of the nonlinear
U(z) element is given by
___
U__ ~.I_ __S_)__~I " ·I~_:nei __~-Y
a
Yl(11)
u(l1)
....-. y(J:1)
u
~.- GII(s) G(s)
nonlinear
...-.
element
Fig. 5.4 Block diagrams of a single-input single-output system with a linear element followed by a static
nonlinearity. (a) The continuous model. (b) The equivalent sampled-data representation.
'1
~
y
nonlinear
element -I G(s) ~
a
,)(k1)
.--.-
u(l1) y(l1)
U T
--/ GII(s) nonlinear
element
G(s) • •
Fig. 5.5 Block diagrams of a single-input single-output system with a static nonlinearity followed by a
linear element. (a) The continuous model. (b) The equivalent sampled-data representation.
Hence, the overall transfer function may be and thus the difference equation takes the form
written as
y(kT) = a1y[ (k -l)T] + a2Y[ (k - 2)T]
Y(Z) = Z[Gh(s)n(kT)G(s)] (5.22) + ... + bou(kT) + b1u[(k -l)T] + ...
U(z)
That is (5.24)
Y(Z) = n(kT)Z[Gh(s)G(s)]
In this case the a and b coefficients depend
(5.23)
U(z) both upon the linear elements and on n(kD.
Nonlinear systems modeling using sampled data 73
Thus, unlike the previous cases where The complete closed-loop transfer function is
the coefficients are constant, these quan-
tities must be evaluated at each sampling Y(z) _
R(z)
z[ + n(kT)G1 (s)
1 n(kT)G1(s)G 2 (s)
] (5.26)
instant.
and the transfer function relating the feed-
5.4.2 CLOSED-LOOP MODELS back variable V(z) to the input R(z) is given
by
The treatment of a closed-loop model which
incorporates a single nonlinear element is
essentially the same as that adopted for non-
V(z) _
R(z)
z[ +n(kT)G1(s)G 2 (S)]
1 n(kT)G1 (s)G 2 (s)
(5.27)
-
y(kTJ
~z) -- T
b
Fig. 5.6 Block diagrams of a feedback system with a static nonlinearity in the forward path. (a) The con-
tinuous model. (b) The equivalent block diagram for the sampled-data representation.
74 Sampled-data models and operator methods
avoid this problem it is conventional to insert x(to) the exact solution of the state equation
a delay of one sample period in the feedback can be shown to be
loop. This does introduce an error, but this I
x=Ax+Bu (5.34)
x[kT] = eAT x[ (k -1)T] + u(kT) f eA1kT -7 Id T
tk-1)T
y=Cx (5.35)
(5.40)
where x(t) is a vector of n state variables, y(t)
However, it is clear that
is a vector of p output variables and u(t) is a
vector of m input variables. A, Band Care
constant matrices of dimension n x n, n x m f
kT
eAlkT-7ldr = f
T
eAIT-TldT (5.41)
and p x n respectively. For an initial condition (k-l)T 0
Operational methods for difference equation models 75
Gis) H G(sJ
y
Fig. 5.7 Block diagram, equivalent to that of Fig. [Link], for a sampled-data representation of a linear multi-
input multi-output system. Double lines have been used to indicate that inputs and outputs to the ele-
ments of the model involve a number of variables. The block G,,(s) represents an array of zero-order hold
elements. The block G(s) is a matrix of transfer functions. It should be noted that the matrix G in the
state-space representation of equations (5.42) and (5.45) is not the same as G(s) in this diagram.
and the differeHce equations representing the equation would have advantages. At each
given continuous model are as follows: time step the responses of all the cascaded
blocks could then be generated sequentially
x(kT) = Fx[ (k -l)T] + Gu(kT) (5.42)
since the output of one element is the input
y(kT) = Cx[(k -l)T] (5.43) to the element which follows it. An
important class of methods has been devel-
where the matrices F and G are given by oped which have origins in operational
calculus and which allow such sequential
(5.44) calculations [1,3]. These operational
and methods are quite widely used for real-time
applications.
G= f
T
eA(T-r)dT (5.45)
o 5.6.1 INTEGRATION OPERATORS
respectively. These two difference equations The process of integration may be repre-
are easily programmed to provide digitally sented by a transfer function
simulated responses for the continuous
model. The matrices need be calculated only (5.46)
once at tpe start of the simulation.
This state-space approach may be extended For the case of the trapezoidal integration
without difficulty to the case of nonlinear algorithm the difference equation represent-
equations. ing the integrator transfer function may be
written as
5.6 OPERATIONAL METHODS FOR
T
DIFFERENCE EQUATION MODELS y(kT) = y[ (k - l)T] + - { u(kT) + u[ (k - l)Tl)
2
It should be noted that in the methods so far
(5.47)
discussed in this section, calculation of the
outputs of a number of different block Taking z-transforms gives
diagram elements involves a separate differ-
T
ence equation for every block considered. Y(z) = z- l y(z) + 2{ U(z) + z- l U(z)} (5.48)
Each of these difference equations relates the
chosen output variable to the excitation and rearranging this gives the transfer
applied to the system input. function
This is not very efficient in terms of com-
puter time, and a method in which each GA (z) = y(z) = T z + 1 (5.49)
element is described by a separate difference U(z) 2 z-l
76 Sampled-data models and operator methods
The same approach may be used for many 5.7 SIMULATION USING OPERATIONAL
other integration methods. However, tech- METHODS
niques such as Runge-Kutta methods or pre-
dictor-corrector methods do not yield a 5.7.1 SIMULATION OF A SINGLE
discrete transfer function so readily since TRANSFER FUNCTION
more than one evaluation is needed of the
A digital simulation model of a given transfer
input and output variables at each sampling function element may be obtained from inte-
interval. Discrete transfer function models
gration operators by carrying out a well-
can also be found for continuous integrators
defined sequence of simple steps. Consider
of order 2 or higher. Rosko [1] provides
the case of a simple second-order linear
derivations for integrators of various orders
system shown in block diagram form in Fig.
for a number of the most widely used opera-
5.S. This has a transfer function
tor methods. Table 5.2 shows operators for a
number of cases for integrators of order 2
G(s) = a2 s + a]s + ao
1-4 inclusive. In some of these it may be (5.50)
b2s2 + b]s + bo
seen that the discrete transfer functions
for the integrals of order 2 and higher The steps involved are as follows:
follow directly from the first order case,
but in others the relationship is not so 1. Divide the numerator and denominator
straightforward. by the highest power of s in the denom-
Table 5.2 Commonly used integration operators for integrators of order 1-4 [1]
Operator
U(s) 2 Y(s)
a 2 s +ats+ao
b2 s 2 +b 1s+bo
U(s) U(z) , . - - - - - - - - - - - - - ,
Y(z)
T
---!/'.- A2 + Alz -1 + AoZ -2
B2 + B1Z -1 + B_7- 2
IT""
Fig. 5.9 Block diagram of the digital simulation model equivalent to that of Fig. 5.8.
78 Sampled-data models and operator methods
U(s) Y(s)
G1(s) Gis) G3(s)
·1 ·1 ·1
a
p.-
Y(z)
-/-1 GlA(z)
P-1 G2/a(z)
P-1 G3,.{z)
Fig. 5.10 Block diagrams for a linear system involving three cascaded transfer functions. (a) The continu-
ous case. (b) The diagram for the equivalent digital simulation model.
IOnJ (5.62)
so that, by the first difference substitution
GI(S)=-(lJ
1+ - G2A () z -_~ -
Y(z)
(5.67)
s z -1 1'; (z)
GAl(Z) =
1+(~
1O(~
Z- 1) =
1
10Tz
z(1+T)-l
U(S)·I S 1~ 1 HL..-_;-J~(S)
z-l)
Fig. 5.11 Block diagram of the continuous linear
(5.63) cascaded system used as an example.
Simulation using operational methods 79
That is can be seen that there is reasonably close
agreement between the simulation results
Y(z) T and the corresponding analytical values. For
--=-- (5.68)
~(z) 1-z-1 a smaller value of T, such as 0.001 s, the
agreement would be very close indeed.
and this gives a difference equation It should be noted that the procedure used
here, involving a direct substitution for the
y(kT) = y[ (k -1)T] + TYI (kT) (5.69) integration operator, differs significantly from
the direct z-transformation method described
Results of a simulation based on the differ- in section 5.3 since insertion of a data-
ence equations (5.65) and (5.69) are shown in reconstruction element is no longer needed.
Table 5.3 for a unit step input with a value of One advantage of the operator method is
T of 0.2 s, together with the response calcu- therefore the fact that it avoids any need to
lated analytically using Laplace transforms. It find the z-transform of the combined data-hold
Table 5.3 Results for the simulation of the cascaded system of Fig. 5.11
for a unit step function input. The table shows the numerical results for
the variables y(kT) and y,(kT), for T = 0.2 s, together with the response
of the continuous model, y(t), found by means of Laplace transform
methods
R(s) yes)
+ 10 1
• s+1 s
1.6
~
~
..... 1.4
1.2
1.0
0.8 -
0.6
0.4
0.2
0.0
0 2 4 6 8 10 12 14 16 18 20
time
Fig. 5.13 Simulated response y(kT) in response to unit step change of reference for the feedback system
of Fig. 5.12. A time increment, T, of 0.01 s was used in this case.
values determined using Laplace transforms, Figure 5.15 shows the error in the calculated
is shown in Fig. 5.14. response variable for the same case as that in
A modification suggested by Tagawa [6] Fig. 5.14. It can be seen from Figs 5.14 and
attempts to overcome the effect of the delay 5.15 that Tagawa's modification leads to a
in feedback by introducing a predictive term simulation result which is closer to the theor-
in the difference equation representing the etical response for the initial part of the time
summing element in the block diagram. For history but that errors are larger elsewhere.
the case of an input signal u(t) and a feedback The overall benefits could possibly be greater
signal y(t), with negative feedback, this gives for other values of T or in other types of
a discrete model for the summer which has application.
the form
5.7.4 SIMULATION OF TIME-VARYING AND
e(kT) = u(kT) - {2y[ (k - I)T] - y[ (k - 2)T])
NONLINEAR MODELS BY OPERATOR
(5.73) METHODS
Although the concept of a transfer function,
As an example of this procedure, consider upon which the operator approach depends,
the closed-loop system of Fig. 5.12. Applying cannot be applied strictly to time-varying and
the procedure outlined above leads to the fol- nonlinear system descriptions, it is possible to
lowing three difference equations: approximate the system dynamics in these
cases by what is known as a 'frozen-system'
e(kT) = r(kT) - {2y[(k -1)T] - y[(k - 2)T]} transfer function. This approach is described
in detail by Rosko [1].
(5.74) An alternative approach to dealing with
time-varying and nonlinear systems by
YI (k T ) _ 10Te(kT)+ YI[(k -1)T]
_____ ;::........0_ _---"-
(5.75) means of operator substitutions is to derive a
l+T detailed block diagram of the given model in
which integrators appear individually. Each
y(kT) = TYI (kT) + y[ (k - 1)T] (5.76) integrator may then be replaced by the
0.010
...0 0.009
t:
V 0.008 .
0.007
0.006
0.005
0.004
0.003
0.002
0.001
0.000
·0.001
{).002
0 4 10 12 14 16 18 20
time
Fig. 5.14 Plot of the error in the simulated response shown in Fig. 5.13. This was calculated from the dif-
ference between y(kT) and the exact response y(t), determined analytically using Laplace transforms, at
each time instant t =0, kT, 2kT, ... , etc.
0.035
g oms
V
0.030
0.020
oms
0.010
0.005
0.000
-().005
-().01O
·M15
-().020
-().025
-(),O3O
0 4 6 10 12 14 16 18 20
time
Fig. 5.15 Plot of error in the simulated response for the feedback system of Fig. 5.12 determined using
Tagawa's modification for a value of T of 0.01 s.
Assessment of performance by frequency response methods 83
chosen operator and the difference equations integration, has an amplitude (1/ wYand
derived for all the linear and nonlinear ele- phase (-7T/2)p.
ments in combination. One advantage of this The accuracy of an integration method may
approach is that Tagawa's modification can be assessed by investigating the deviation of
be incorporated wherever there is a feedback the frequency response of a single integrator,
pathway in the model and this may enhance or a multiple integrator, from the theoretical
the overall accuracy [6]. The speed advan- response. This comparison is carried out most
tages of operator methods, in comparison conveniently by determining the ratio of the
with more complex integration methods, are transfer function I'(jw) of the numerical inte-
retained. Care must be taken, however, in grator to the corresponding transfer function
using an approximate approach of this kind, I(jw) of the ideal integrator. That is, by con-
and verification of the simulation is essential. sidering the quantity
This can be carried out by comparing results
obtained using this approach with selected
A (J'w ) = 1'(jw ) (5.77)
results from use of an established, but slower, I(iw)
simulation package.
The larger the range of frequencies over
which A(jw) has a magnitude close to unity
5.8 ASSESSMENT OF SIMULAnON and a phase of zero the more accurate is the
PERFORMANCE BY FREQUENCY integrator.
RESPONSE METHODS Consider the case of the Tustin integration
operator given by the transfer function
The performance of linear continuous system
simulations can be assessed very conve-
1'(z) = T z + 1 (5.78)
niently using frequency response methods.
2 z-l
The accuracy of the simulation may be deter-
mined by investigating the extent to which
Since z = exp(sT) and s = jw for frequency
the frequency response of the simulation
response evaluation this gives
model deviates from the frequency response
of the linear model from which it was T jwT 1
derived. This approach is most widely used I'('w)= e + (5.79)
J 2 e}wT -1
with transfer function descriptions and is par-
ticularly popular in control system applica- = T cos wT + j sin wT + 1
tions, reflecting the widespread use of (5.80)
2 cos wT + j sin wT - 1
frequency-domain methods in control engi-
neering [7-9].
After simplification the resulting expressions
for magnitude and phase are
5.8.1 FREQUENCY RESPONSE OF DISCRETE
n
INTEGRATORS l1'(jw )1 = T cot wT Phase=-- (5.81)
2 2 2
The transfer function of an integrator (1/ s)
has a frequency response l/jw. The ideal inte- Similar processes can be carried out for
gration operator (1/ s) thus has a frequency other types of integration operator. Figure
response characteristic that has an amplitude 5.16 shows relationships between IA(jw) I
which is inversely proportional to frequency and wT for a number of different integration
and a phase lag of n/2. In general, an integra- operators and Fig. 5.17 shows the correspond-
tion operator (1/ sY, where p is the order of ing phase plots.
84 Sampled-data models and operator methods
1.4
/
/
/ 600
I/A
II -.l. /
/ Ve /
~
VB
3' // L 40
0
;::>
< ...--V >-- 1-1- /
- 1.0
-- r-
....... 1' C V
/
V
""
0
20 v
\ /
V
1\ L
~ l.--"
A,C
00
0.2 0.5 1.0 3.0
0.6 ColT
0.2 0.5 1.0 3·0
Fig. 6.2 Listing of the SLIM program VDPOLS. SLI for simulation of a Van der Pol oscillator.
shows the information which appears on the the simulation run illustrated in Fig. 6.3.
computer screen during the running of this Figures 6.5 and 6.6 show typical graphical
simulation. No input file is needed so no file results produced using the post-processing
name is provided. This table also includes the program SLIMPLOT. In Fig. 6.5 the results
first few values of the output as it appears on correspond to selection of channel 3 of the
the screen. results file as the y variable and channel 2 as
Figure 6.4 shows a subset of the results the x variable in an x-y plot. In Fig. 6.6 the
as they appear within the output file same results are shown as two separate time
VDPOLS . OUT which is created as a result of histories.
(a) SLIM
ENTER NAME OF SIMULATION FILE
[Link]
FILENAME FOR INPUT DATA? (HIT RETURN IF NONE)
OUTPUT FILENAME?
[Link]
RESULTS TO BE OUTPUT TO TERMINAL? (YIN)
:t..
3
400
.400000E-OI
.OOOOOOE+OO .500000E+00 -.500000E+00
.400000E-Ol • 479296E+00 -.535330E+00
.800000E-OI • 457164E+00 -.571343E+00
.120000E+00 .433579E+00 -.608069E+00
.160000E+00 .40B509E+00 -.645539E+00
• 200000E+00 etc. etc .
etc.
Fig. 6.3b The screen display at the post-processing stage using the SLIMPLOT program.
Fig. 6.3 Copy of screen display for the Van der Pol oscillator simulation run. Prompts from the computer
appear, together with responses from the user (which are shown here underlined).
100 The principles of equation-oriented simulation software
3 0 -1
.400000E-01 .100000E-05 .100000E-03 .lOOOOOE-03
3
400
.400000E-Ol
.OOOOOOE+OO .500000E+OO -.500000E+00
.400000E-Ol .479296E+00 -.535330E+00
.800000E-Ol .457164E+OO -.571343E+OO
.120000E+OO .433579E+OO -.608069E+00
. 160000E+OO .408509E+OO -.645539E+OO
.200000E+00 .381925E+00 -.683783E+00
.240000E+00 .353796E+00 -.722824E+00
• 280000E+00 • 324088E+00 -.762679E+00
• 320000E+00 .292771E+00 -.803358E+00
. 360000E+OO . 259809E+OO
.400000E+OO • 225171E+OO
.440000E+OO • 188826E+00
.480000E+OO .150743E+00
.520000E+00 • 110897E+00
.560000E+00
.600000E+OO
.640000E+00 .598834E+00
.680000E+00 .618052E+00
.720000E+00 .636853E+00
.655397E+00
-.167311E+Ol .673833E+00
-.164579E+Ol .692295E+00
-.161773E+Ol • 710907E+00
• -.158891E+Ol . 729785E+00
• 158800E+02 -.155934E+Ol .749043E+00
.159200E+02 -.152898E+Ol .768785E+00
.159600E+02 -.149783E+Ol • 789118E+00
.160000E+02 -.146585E+Ol .810145E+OO
.160400E+02 -.143301E+01 .831971E+00
Fig. 6.4 Part of file VDPOLS. OUT created during the Van der Pol oscillator simulation run of Fig. 6.3. The
first three values in the file are all integers and convey diagnostic information of little interest to the user.
The next four items provide the communication interval, the absolute and relative error values specified
in the simulation program and the minimum interval value. All this information is provided automati-
cally at the start of each print file but is not used by the SLIMPLOT program to produce graphical output.
The subsequent items in the file all correspond to TYPE statements in the simulation program and all are
used by SLIMPLOT. The first three of these represent the number of channels being output (three in this
case, corresponding to the variables T, Xl and X2 in the program), the number of samples in the simula-
tion run (400 in this case) and the communication interval (0.04 in this case). The values which follow all
represent the response information for the simulation run, with the three columns corresponding to
values of T, Xl and X2 at each communication interval. Note that the value of the communication interval
appears twice in the header data.
Figure 6.7 shows a slightly modified version Although Van der Pol's equation is rela-
of part of the initial segment of the simulation tively simple in form, it shows some very
program in which the parameter value for the interesting features in terms of its properties.
coefficient, IL, and initial conditions in the It can be seen from the results of Figs 6.5 and
system equations are read from a data file. 6.6 that, for the parameter values and initial
Figure 6.8 represents this input data file, and it conditions chosen, the system response is
may be seen from this that the data are identi- oscillatory and shows a maintained oscillation
cal to that incorporated within the program of of constant amplitude, after an initial tran-
Fig. 6.2. In running the modified program the sient. This is because for small values of y (y
procedure is exactly as described above except less than unity) the damping in this second-
that the input file name, VDPOLS . IN, is given order system is negative. As y increases, the
in response to the prompt from the SLIM value of the factor (if - 1) eventually becomes
program. Results in this case are identical to zero and then positive. One would therefore
those found previously. expect intuitively that the oscillation should
Examples of other simulation languages 101
VRN DER POL OSCILLRTOR SIMULATION
r-~l
--1
P~ogam; [Link]
oot.a f -c L!:? : non!:?
Fig.6.S Simulation results for the Van der Pol oscillator simulation in the form of an x-y plot. Here X, is
the x-axis variable and x, is the y-axis variable.
build up in amplitude until a point is reached values of the parameter IL, including nega-
at which the nonlinearity starts to limit the tive values. A further useful exercise is to
amplitude because the damping has changed modify the simulation program slightly to
in sign. The response thus shows typical 'limit allow study of a system described by the
cycle' behavior. equation
The limit cycle feature of the response
shows up very clearly in the phase plane ----t
d
dt
2
+ IL (y2 - 0.1 y4 -1) + y = 0 (6.4)
type of plot involving X2 plotted against Xl'
On the one hand, if separate simulation runs It is left to the reader to modify the program
are performed for a number of different sets VDPOLS . SLI to investigate the nature of the
of initial conditions it is found that a closed model behavior in this second case.
path is defined in the phase plane to which
all the trajectories tend regardless of the
6.7 EXAMPLES OF OTHER SIMULATION
starting point. This is typical of systems
LANGUAGES
which have stable limit cycle behavior.
Unstable limit cycle behavior, on the other
6.7.1 CSSLIV
hand, would give rise to a locus in the phase
plane from which all the trajectories would CSSL IV [9] is a simulation language which
tend to diverge. Therefore, it is interesting has a form which is very strongly linked to
to investigate the way in which the response the 1967 recommendations for continuous
of the system varies with different sets system simulation languages. CSSL IV is
of initial conditions and with different translator based and converts the source
102 The principles of equation-oriented simulation software
VAN DER POL OSCILLATOR SIMULATION
l
I
l "
o ."
-1
-2
\ ..
'" I
P~ogam: [Link]
Data f{.Lel none
KEY
-- COLUMN 2
-- -' COLUMN 9
Fig. 6.6 The results of Fig. 6.5 displayed as separate plots of Xl and X, versus time. Here the continuous
line represents Xl (column 2 of the output file and the second variable in the TYPE statement in the
DYNAMIC segment of the program) while the dashed line represents X, (column 3 of the output file and the
third variable of the TYPE statement).
Fig. 6.7 Listing showing a modification to part of the initial segment of the SLIM program VDPOLS. SLI
to allow data for the initial conditions (XIO and X20) and the program variable AMU (the parameter J.L of
the mathematical model) to be input from a data file.
program written in the CSSL IV language into generates the equivalent of about 50-250 lines
a set of Fortran 77 programs which are subse- of Fortran. CSSL IV has powerful macro
quently compiled. Each line of CSSL IV code facilities.
Examples of other simulation languages 103
1.0 tinuous system simulation. This is achieved
0.5 by the provision of a DISCRETE section in
-0.5
which difference equations can be coded to
Fig. 6.8 Listing of the data file to be used with the represent discrete elements within a mixed
modified program of Fig. 6.7 for the Van der Pol continuous and discrete model. The DIS-
oscillator simulation. CRETE section is equivalent to a DERIVATIVE
section in a purely continuous simulation
CSSL IV provides for a clear separation of program.
the model definition from the experiments to Simulation programs written using ACSL
be carried out using that model. A powerful have two parts which are quite separate. The
runtime monitor is provided which allows first of these is the model definition, while
the user to control the simulation and carry the second is the set of commands required
out all the necessary interactive processes to control the simulation at runtime. The
such as changing parameters, collecting data benefit which this type of approach gives is
and presenting results. Experiments may be that the model may be defined once
interrupted using a 'break key', which allows and analyzed repeatedly using runtime com-
a simulation to be stopped while remaining mands which are specified either interac-
within the CSSL IV environment. tively or in a batch mode of operation.
In addition to the standard facilities Figure 6.9 is a listing of an ACSL program
expected in simulation languages which for the simulation of the system described
adhere to the 1967 recommendations, CSSL by Van der Pol's equation. Many similar-
IV incorporates a number of useful addi- ities may be found when this program
tional features. As well as special-purpose is compared with the SLIM version in
facilities associated with the simulation of Fig. 6.2.
continuous and mixed continuous-discrete As with CSSL IV, the model definition in an
systems, CSSL IV incorporates analysis tools ACSL program is translated into Fortran. The
such as a steady-state finder, FFT (fast translated model is then loaded with the
Fourier transform) processing and lineariza- runtime library which reads and interprets
tion routines. the set of runtime commands. A variety of
CSSL IV is available on all major computer forms of output may be produced, in both
systems. CSSL IV simulation models are tabular and graphical form.
highly portable and the accuracy is claimed to ACSL allows the user considerable scope in
be machine independent. terms of interactive usage. Simulation com-
mands may be issued individually and data
may be modified at the time of the simulation
6.7.2 ACSL
experiment. Unlike many other simulation
ACSL (Advanced Continuous Simulation languages, ACSL imposes no limits for inter-
Language) [10] is an equation-oriented simu- nal tables used by the translator and thus has
lation language that has a form which corre- no arbitrary limits on the number of state
sponds closely to the CSSL specification [6] in variables, the number of symbols or the
many respects. It has a number of features number of labels in the simulation program.
which exceed the requirements of the 1967 As with most other modern simulation lan-
specification very considerably and it is guages, ACSL has a macro capability which
undergoing a continuous process of upgrad- allows the set of special ACSL statements to be
ing and improvement. ACSL, as with most extended. Macros are similar to a subroutine
other commercial packages currently avail- or function in conventional high-level pro-
able, has facilities for mixed discrete and con- gramming languages. Once defined the macro
104 The principles of equation-oriented simulation software
PROGRAM Van der Pol Oscillator simulation
INITIAL
CONSTANT mU=1.O,yic=O.5,ydic=O.O,tend=15.0
CINTERVAL CINT=O.04
END $ "of INITIAL"
DYNAMIC
DERIVATIVE
ydd=-mu*(y*y-l)*yd-y
yd=INTEG(ydd,ydic)
y=INTEG(yd,yic)
END $ "of DERIVATIVE"
TERMT ([Link])
END $ "of DYNAMIC"
TERMINAL
END $ "of TERMINAL"
END $ "of PROGRAM"
Fig. 6.9 Listing of an ACSL simulation program for the Van der Pol oscillator. Note the explicit definition
of the INITIAL, DYNAMIC and TERMINAL segments and form of syntax used. Comments are denoted in
this listing by means of a $ symbol.
Fig. 6.10 DESIRE program listing for the Van der Pol oscillator. In this case the program has been written
in such a way that multiple runs are being performed for a range of values of the parameter 'mu'.
DESIRE has a command-mode interpreter statement 'drun' which calls for a simulation run. In this
example, use has been made of the similar statement 'drunr', which is an extension of the 'drun' type of
statement for use in repetitive simulation experiments in which initial conditions must be reset before the
start of each separate run.
SLIM program. Once again, typical applica- The absence of facilities to model a pure
tions involve much smaller programs. As has time delay (transport delay) is one important
been explained in section 6.7, some simula- example of a feature available in many simu-
tion languages, such as ACSL, do not impose lation packages which has not been included
limitations of this kind. within the version of SLIM provided with this
While the restrictions in terms of program book. The lack of any facility to represent a
size are not thought to represent limits which function in tabular form is another example
are likely to be approached by the beginner, it of a potentially useful feature which has not
is important to understand that languages been included. Problems which require the
such as CSSL IV, ACSL and DESIRE provide use of such facilities cannot be tackled using
many features which are not available within SLIM unless the problem can be reformulated
SLIM. Some of these restrict the type of in some way to allow it to be solved within
system which can be modeled using SLIM the limitations of the language. Other lan-
and others limit the options available to the guages provide a much more extensive set of
user at runtime. The subset of facilities incor- facilities and therefore provide much more
porated in SLIM has been chosen to allow the flexibility at the model definition stage.
beginner to gain experience rapidly, without In comparison with ACSL and other similar
having to deal with the complexities associ- languages SLIM has limited capabilities in
ated with certain types of feature which are terms of runtime options. There is no real sep-
included in more comprehensive languages. aration of the model definition and the exper-
A careful examination of the user manual for iment in SLIM. Changes of parameters
each specific language allows the potential require editing of the simulation program
user to establish the precise differences in itself or changes of an associated data file.
every case. However, flexibility is available, in terms of
References 107
output display facilities, in SLIM because of 7. Control Data (1970) MIMIC - A Digital
the fact that the graphical output options Simulation Language Reference Manual, revision
involve a separate post-processing program. D, Publication No. 44610400, Control Data,
Minneapolis, Ml.
The reader can gain some appreciation of 8. Control Data (1971) Continuous System Simu-
the facilities of languages other than SLIM in lation Language, Version 3, User's Guide, Control
the applications chapter which include some Data, Sunnyvale, CA.
comparisons of SLIM with other implemen- 9. Simulation Services (1984) CSSL-IVReference
tations of the same simulation problems. Manual, Simulation Services, Chatsworth,
SLIM can be extremely effective when learn- CA.
ing about the principles of simulation and 10. Mitchell, E.E.L. and Gauthier, J.S. (1976)
Advanced continuous simulation language
when used for the solution of problems for
(ACSL). Simulation, 25, 72-78.
which its facilities are appropriate. It is 11. IBM (1968) System/360 Continuous System
important, however, to understand its limita- Modeling Program (360A-Cx-16X) Application
tions and to have an appreciation of the addi- Description, IBM Publication GH20-0240-2.
tional facilities of the software products 12. Speckhart, F.H. and Green, W.L. (1976) A
which are available commercially. Guide to Using CSMP, Prentice-Hall, Engle-
wood Cliffs, NJ.
13. Ricketts, l.W. (1977) MIMESIS - a continuous
REFERENCES system simulation language. University of
Dundee. PhD thesis.
1. MathWorks (1993) MATLAB User's Guide, 14. Ricketts, l.W. and Dickie, A.A. (1978)
MathWorks, Natick, MA. MIMESIS - a machine independent CSSL for
2. Billmann, L., Mirab, H. and Winkler, U. (1992) minicomputer systems, in Proceedings of
CACSD - CASE tools. Measurement and the UKSC Conference on Computer
Control, 25, 137-43. Simulation, 1978, Chester. IPC, Guildford,
3. Cellier, F.E. and Rimvall, CM. (1989) Matrix pp. 193-203.
environments for continuous system model- 15. Kom, G.A. (1989) Interactive Dynamic System
ling and simulation. Simulation, 52 (4), 141-49. Simulation, McGraw-Hill, New York.
4. Hindmarsh, A.C (1983) ODEPACK, a system- 16. Kom, G.A. and Kom, T.M. (1993) Simulating
atized collection of ODE solvers. IMACS dynamics and neural networks under DOS,
Transactions on Scientific Computing, I, 55-64. OS/2 and UNIX, in Proceedings of the 1992
5. Petzold, L.R. (1982) DASSL: A Differential/ EUROSIM Conference, EUROSIM'92, 28
Algebraic System Solver, Report No. SAND82- September-4 October, 1992, Capri (eds
8637, Applied Maths Division, Sandia F. Maceri and G. Iazolla), North-Holland,
National Laboratories, Livermore, USA. Amsterdam, pp. 103-15.
6. SCi Software Committee (1967) The SCi 17. Kom, G.A. (1991) Neural Network Experiments
Continuous-System Simulation Language. on Personal Computers and Workstations. MIT
Simulation, 9 (6), 281-303. Press, Cambridge, MA.
THE PRINCIPLES OF BLOCK DIAGRAM-
ORIENTED SIMULATION TOOLS
7
f J x
Fig. 7.1 Block diagram of the Van der Pol oscillator system.
114 Principles of block diagram-oriented simulation tools
iC2 iC1
COli
unity 1 - - - - '
b
Fig. 7.2 Description of the Van der Pol oscillator using the BIOPSl interactive simulation program.
(a) The BlOPSl structure diagram for the simulation. (b) The corresponding table defining the program. In
the structure diagram, block types are indicated in the upper right corner of each block, and in
this example involve INT (integrator), MUL (multiplier), GAl (gain factor) and CON (constant). The block
name is shown in the middle of each block. Circles indicate a SUM (summer) or SUB (subtractor) block, and
in these cases the names are shown next to the symbol. Note that the numerical values (0.5 and -0.5) for
the initial conditions of the two integrator blocks appear only in the Parl column of the program listing.
x1
Fig. 7.3 A SIMULINK diagram for solution of the Van der Pol oscillator problem. The meaning of each of
the SIMULINK icons becomes obvious when this diagram is compared with Fig. 7.1. Note the use of the
'scope' icons to indicate that the variables xl and x2 are to be displayed.
r ... . ,
0 INT2
1 0 INTl
L4~l . ..
MU
.[D ...
1 ... 00e--+-oo
·,
··
1__ ... __ ... ____ ................. _ ....... _____ ...... __ ... __ ~
1_00e--+-Oo
Fig. 7.4 An XANALOG block diagram for the Van der Pol oscillator problem. Note that in this diagram
the dotted lines indicate quantities which have been inverted (multiplied by -1).
dt dt C j
f
- + Mdi2- +1- ij dt (8.1)
described by a linear fourth-order ordinary
Similarly, for the closed path involving the
R 11 current i2(t), application of Kirchhoff's voltage
M
law gives an equation of the following form:
O=l2. R2 + L2-+M-+-
di2 di 1 l2 dt
j J. (8.2)
-1O>-~ 1 - - -.....
dt dt C2
di2
dt =-
1.
T2 12 -
2
W 2
f'
12
d di1
t - K2 dt (8.4)
where
(8.6)
(8.7)
w~ 1
- X3 - X4
1- K j K 2 T2 (1- K j K 2 )
K 2e(t)
+-----==---.:....:...- (8.8)
Lj (1-K j K 2 )
gives the four state equations The four state equations above may be
incorporated into any equation-oriented sim-
(8.5) ulation language program or may be simu-
Fig. 8.3 Listing of part of a SLIM program (COUPLE. SLI) for simulation of the coupled-tuned electrical
circuit.
lated easily by writing a program using a three examples of appropriate parameter sets
general-purpose high-level programming lan- with the corresponding data file names. Any
guage. Figure 8.2 is a block diagram which other situation may be investigated easily by
corresponds to this set of state equations. creating a similar data file containing the nec-
Figure 8.3 shows the listing of the DYNAMIC essary parameter values. Results from this
segment of a SLIM program (COUPLE. SLI) SLIM program are shown in Fig. 8.4 as plots
for this system model for a case where the of the currents i,(t) and i2 (t) versus time for
input, e(t), is a step change of voltage, and the two different sets of conditions. These results
response of interest involves the primary and show very clearly the coupled nature of the
secondary currents. The complete program is dynamics and the periodic exchange of
included on diskette. The initial segment of energy between the primary and secondary
this program involves a set of parameters circuits and the dependence of the behavior
which are read in from a file. Table 8.1 shows of the system on the time constants T, and T2 •
Table 8.1 Parameter values and SLIM data file names for coupled-tuned circuit
example
6 ~
7\\,
'
o '-7~
o 20 40 60 80 t
a (s)
10
5
.-
\
0 )
-5
-10
0 20 40 60 80 t
b (s)
Fig. 8.4 Simulated responses of the coupled-tuned electrical circuit for a step input of applied voltage.
The responses shown are the currents (measured in amperes) in the primary circuit (continuous line)
and secondary circuit (dashed line) for parameter values corresponding to the data file COUP1. IN (a) and
COUP3. IN (b). Note the much lower level of damping in the responses of (b) compared with those
shown in (a). This difference is because of the choice of parameters T, and T2 in the two cases considered,
as indicated by the figures in Table 8.1.
Block diagram-oriented simulation pro- Ideas for further things to tryout on your
grams based on XANALOG and SIMULINK own using this model include investigations
may be derived from Fig. 8.2 and can have an of other parameter combinations involving
identical structure. Results obtained using variation of Wj and W 2 and K and K2 • Those
j
these tools are exactly the same as those with an appropriate understanding of electri-
shown above. The powerful runtime cal circuits may also want to use the simula-
command facilities of XANALOG and tion model to look at the effect of changing
SIMULINK provide a more user-friendly the sign of the mutual inductance coupling
interface for interactive experimentation than the left- and right-hand circuits. This would
is possible with SLIM. Similar runtime involve changing the structure of the equa-
command facilities and high levels of interac- tions which form the simulation model and
tion could, of course, be achieved with other would thus require modification of the source
equation-based languages such as ACSL and program. Since this is a linear model, compar-
with other block diagram-based tools. isons may be made of the simulation methods
A control system simulation 121
with analytic approaches based, for example, useful illustration of continuous simulation
on Laplace transforms. techniques being used to approach an import-
ant class of practical engineering problem.
A number of controller transfer functions
8.3 A CONTROL SYSTEM SIMULATION:
can be considered for this type of application.
SPEED CONTROL OF A WATER TURBINE
Detailed consideration of possible forms of
Control systems investigations form one of controller are not of primary interest in the
the most common application areas for con- present context, but any readers interested in
tinuous system simulation tools. This this aspect, or in the underlying model, may
example is concerned with a simulation of a find further details elsewhere (e.g. refs 1 and
speed-control system for a water turbine used 2). For the present purposes, the controller
for the generation of electricity. Figure 8.5 is a considered is of the 'temporary droop' type
block diagram of a highly simplified descrip- implemented in continuous form. This type
tion of a system of this kind. The turbine, of controller has a transfer function of the
which is of the impulse type, is represented form:
by a linear transfer function, which also incor-
porates dynamics of the pipeline system.
Control of the turbine is accomplished
through feedback of a signal proportional to
turbine shaft speed and comparison with
the desired (reference) speed signal. Any dif- In this controller transfer function only the
ference between the desired and actual speed parameters (J" and JL may be regarded as
provides an error signal, which is processed adjustable quantities. The other parameters Tx
by the controller block to provide the signal and Ty) are fixed quantities and are not avail-
applied to the turbine actuator input, to able for controller tuning.
change the water flow in such a way that the Appropriate parameter values for the plant
error is continuously minimized. The con- and controller are as follows:
troller may be implemented either in analog
(continuous) or digital form. The model is Ta = 7.0 s (inertial time constant)
linear throughout, apart from backlash in the Tw = 1.1 s (water time constant)
mechanical linkages associated with the Ts = 0.2 s (actuator servo time constant)
turbine inlet actuator. The purpose is to allow
the performance of the model of the closed- All model variables, such as turbine speed,
loop system to be investigated. This is a are expressed as normalized (unit) quantities.
TL
+
N ref + V 1 v 1 N.
CIs) N.L.E.
1 + sT. sT.
Fig. 8.5 Block diagram of closed-loop system for automatic control of the speed of a water turbine.
122 Simple examples using common simulation tools
Figure 8.6 shows the DYNAMIC segment of a program listing. The simulation experiment,
SLIM program (TURB. SLI) for the simulation in this case, involves investigation of the
model and the full source code may be found response of the control system to a step change
on the diskette. Parameter values for the of the reference speed. It can be seen from Fig.
nominal set of conditions are provided in the 8.7 that the response is stable but oscillatory in
C
C ******************Start of Dynamic Segment*************
DYNAMIC
C
C ******************Start of Derivative Segment************
DERIVATIVE
DT1= (V-T1) *( 2.0 ITW)
T1=INTEG(DT1,T10)
TO=T1-TW*DT1
TA=TO-TL
DNS=TA/TIA
ANS=INTEG(DNS,ANSO)
E1=REF-ANS
DDY1=(E1-SIG*Y1-( (SIG+AMU)*TX+TY) *DY1) I (TX*TY)
DY1 =INTEG (DDY1,DDY10)
Y1=INTEG(DY1,Y10)
Y=Y1+TX*DY1
DV= (Y-V) ITS
V=INTEG(DV,VO)
DERIVATIVE END
C
C Values of t, ref, tl and ans for current communication
C interval output to results fil
TYPE T,REF,TL,ANS
C
C Test for end of simulation run
IF (T-TMAX) 10,10,12
10 DYNAMIC END
C
C *************************Terminal Segment*******************
12 STOP
END
Fig. 8.6 Part of a SLIM program (TURB. SLI) for the simulation of the water-turbine speed-control system
with a 'temporary droop' type of controller transfer function. No backlash is included in the simulation
model in this case.
<II
Z 0.2
0.1
o 20 40 60 t
(5)
Fig. 8.7 A typical simulated response of the water-turbine speed-control system to a step change in the ref-
erence speed. The choice of controller parameters used in this case gives a transient response which is oscil-
latory in nature. The adjustable parameters of the controller (a and fL) could be altered, through further
simulation experiments, to try to establish conditions giving a more satisfactory transient performance.
A control system simulation 123
.... One interesting extension to this simulation
=
IL.
.... program involves the inclusion of mechanical
=
= backlash between the servomotor and the
turbine inlet. This is an important feature of
real mechanical systems of this kind and is
known to have a destabilizing effect on the
INPUT overall control system. Backlash is a double-
valued form of nonlinearity which has a
steady-state input-output characteristic of the
form shown in Fig. 8.8. It generally arises
because of 'slack' in mechanical linkages and
gears. With backlash present, movement of the
input in one direction produces a proportional
movement of the output, but any reversal of
Fig. 8.8 Typical input-output characteristics of a the direction of the input will cause the output
backlash element. to stop before following the input again.
The representation of backlash and other
nature for the controller parameter values similar double-valued nonlinearities which
used. A program TURBl . SLI, which allows involve a form of 'memory', such as hyster-
input from a data file, is also included on the esis, can be difficult. One approach is to use
diskette supplied with this book. This version principles first established for the modeling of
of the simulation program, together with the backlash elements using analog computers
appropriate data files in the same format as the [3]. This involves the use of an integrator with
test file TURB1. IN, provides a convenient feedback to provide the memory element.
means of experimenting with controller para- Figure 8.9 shows a block diagram for a repre-
meter values. sentation of backlash by this type of method.
--,
: OR
'--D
.-..
,
I
: ,
,
_J
,,
,
:..... -_ .. -
,
I
Fig. 8.9 Block diagram illustrating one possible technique for representation of backlash within a simula-
tion model. Variables shown on this diagram are consistent with the variables of Fig. 8.5 and with vari-
able names used in the SLIM program listing of Fig. 8.10.
124 Simple examples using common simulation tools
Essentially, the system works by comparing the integrator output to remain fixed in value
the integrator output, which is the output for a period of time, and this corresponds to
variable for the backlash element, with the the flat top and bottom sections on the
input. The integrator input is controlled diagram showing the input-output character-
through comparator elements and logic istics for the backlash. Although presented
blocks, and when the input reverses direction here in block diagram form, the implementa-
(at the extremes of travel) the integrator input tion of a backlash model of this kind is quite
becomes zero. The integrator input remains simple using the equation-oriented methods.
zero until the input variable has changed by Figure 8.10 shows the DYNAMIC segment of a
an amount equal to the width of the backlash SLIM program TURB2. SLI, which is a
element. It is then switched back to the version of the turbine speed-control system
output of the summing element. This causes model with backlash incorporated. This
C
C ******************start of Dynamic Segment**************
DYNAMIC
C
C ******************Start of Derivative Section***********
DERIVATIVE
C
C Section for representation of backlash element
IF(AB-0.00l)40,45,45
45 DB=V-YB
IF(DB+AB)60,50,50
60 FB=DB
GOTO 100
50 IF(DB-AB)80,70,70
70 FB=DB
GOTO 100
80 FB=O.O
100 YB=INTEG(10.0*FB,YBO)
GOTO 110
40 YB=V
C
C End of section for backlash
110 DT1=(YB-T1)*(2.0/TW)
Tl=INTEG(DTl,TlO)
TO=TI-TW*DTI
TA=TO-TL
DNS=TA/TIA
ANS=INTEG(DNS,ANSO)
El=REF-ANS
DDYl=(El-SIG*Yl-«SIG+AMU)*TX+TY)*DYl)/(TX*TY)
DY1=INTEG(DDYl,DDYI0)
Y1=INTEG(DYl,YlO)
Y=Yl+TX*DYl
DV=(Y-V)/TS
V=INTEG(DV,VO)
DERIVATIVE END
C
C ********************End of Derivative section**************
C
C Values of t, ref, tl, ans and yb for current communication
C interval output to results file
TYPE T,REF,TL,ANS,YB
C
C Test for end of simulation run
IF(T-TMAX)10,10,12
10 DYNAMIC END
C
C *************************Terminal Segment******************
12 STOP
END
Fig. 8.10 Part of a SLIM program (TURB2 . SLI) for the water-turbine speed-control system with mechani-
cal backlash included between the servomotor and the turbine inlet valve.
Simulation of a neural encoder model 125
program and the corresponding data file tice. For those with an appropriate level of
TURB2 . IN are available on the diskette. understanding of control systems analysis
The oscillation observed on the response techniques, the simulation result for the criti-
shown in Fig. 8.7 is influenced by the backlash cal value of the backlash parameter may be
parameter ab as well as by the parameters of compared with predictions from theory,
the controller. A well-damped step response based upon describing function analysis. It is
in the absence of backlash (a b = 0.0) can interesting to consider the reasons for any dif-
become a maintained oscillation if the back- ferences between the simulation result and
lash is increased sufficiently (Fig. 8.11). the value predicted by theory. It should be
Figure 8.12 is an ACSL program listing for remembered that the describing function
this example. Note that a special statement approach involves some important simplify-
BCKLSH is available in ACSL to represent ing assumptions and approximations. The
backlash effects and use has been made of significance of some of these can be investi-
this facility. Figure 8.13 shows results gated easily using the simulation.
obtained from this ACSL program. In modern
block-oriented tools such as XANALOG and
8.4 SIMULATION OF A NEURAL
SIMULINK special facilities are also available
ENCODER MODEL
for the simulation of nonlinear elements such
as backlash. Figure 8.14 shows an XANALOG Sensory receptors of the peripheral nervous
block diagram for this problem. system have an important role in the control
Further investigations which could be systems which are used in the maintenance of
carried out using this model and the SLIM posture and in the control of movement. One
simulation programs include determination very important receptor in the neuromuscu-
of the minimum value of the backlash para- lar control system is the muscle spindle.
meter ab which gives rise to a maintained Muscle spindles are embedded in skeletal
oscillation in terms of turbine speed. This muscles and provide neural signals that are
'limit cycling' type of behavior is clearly an sensitive to muscle stretch and to the rate of
undesirable situation in a speed-control sys- change of stretch. The precise role of the
tem and one that should be avoided in prac- muscle spindle as an active element in the
0.2 ",
t \
\
..
.1 \
\ ,
.<"
"
, I
l \
0 .. I
-'
-0.2
o 20 40 60 80 t
(s)
Fig. 8.11 Simulated step response of the water-turbine speed-control system with backlash which is
almost sufficient to produce a continuous oscillation. The dashed line shows the response in terms of the
output speed (N,), while the continuous line is the output Yb from the backlash element. This is typical of
the response of a nonlinear closed-loop system when operating in the vicinity of a stable limit cycle.
126 Simple examples using common simulation tools
PROGRAM [Link]
"Simulation of governed hydro-turbine system with"
"a temporary-droop governor and with backlash between"
"servo-motor and turbine inlet actuator"
INITIAL
"Data for plant and governor transfer functions"
CONSTANT tw=l.l, ts=O.2, tia=7.0
CONSTANT ab=O.06 $ "backlash parameter"
CONSTANT tx=16.0, ty=O.3, sig=O.03, mu=O.25
"Data for reference input and disturbance input"
CONSTANT ref=O.l, tl=O.O
"Data for experiment duration"
CONSTANT tend=99.9
ALGORITHM IALG=4 $ "Runge-Kutta second order"
CINTERVAL CINT=l.O
"Initial conditions"
tl=O.O
v~o.
yl=O.O
dyl=O.O
ybO=O.O
END $ "of INITIAL"
DYNAMIC
DERIVATIVE
"Turbine"
dtl={yb-tl)*{2.0/tw)
tl=INTEG{dtl,O.O)
to=tl-tw*dtl
"Load"
ta=tO+tl
dns=ta/tia
ns=INTEG{dns,O.O)
"Governor"
el=ref-ns
ddyl={el-sig*yl-{{sig+mu)*tx+ty)*dyl)/{tx*ty)
dyl=INTEG{ddyl,O.O)
yl=INTEG{dyl,O.O)
y=yl+tx*dyl
"Servo motor ll
dv={y-v)/ts
v=INTEG{dv,O.O)
yb=BCKLSH{ybO,ab,v)
END $ "of DERIVATIVE"
TERMT ([Link])
END $ "of DYNAMIC"
TERMINAL
END $ "of TERMINAL"
END $ "of PROGRAM"
Fig. 8.12 ACSL program listing for the water-turbine speed-control example.
overall control system is not yet fully under- modeled in an approximate way by the
stood but much experimental work has been charging and periodic resetting of a leaky
carried out on the properties of the muscle integrator when its output reaches a thresh-
spindle and many modeling investigations old level. This is an interesting example of a
have been undertaken (e.g. refs 4 and 5). The discontinuous problem and provides a simu-
'encoder', which is the element responsible lation model which has some very interesting
for the conversion of stretch within the nonlinear properties when subjected to peri-
central section of the muscle spindle to a odic inputs representing mechanical stimuli
neural signal in pulse frequency modulated applied to the muscle. Further information
form, provides an interesting simulation about the use of leaky integrator and thresh-
example involving a highly nonlinear old models for neural encoders may be found
problem. The generation of action potentials in the biomedical engineering literature (e.g.
in the nerve leading from the receptor arises refs 6-9). It is also of interest to note that these
through membrane properties which can be 'integrate and fire' models are now creating
Simulation of a neural encoder model 127
i
c
TURBTD2 AB=O . 0 ~ TURBTD2 AB=O.02
[;;
· d;~-+r
(\
[I~
· 1\
~d:+-r4
J' VV'
(f)
~ ~
Z
(f)
~
· d:~-+4
ru
-
c
Fig. 8.13 Responses obtained using the ACSL simulation (a) with backlash (a b =0.02) and (b) without
(a b = 0.0). Note, once again, the destabilizing effect of the backlash.
LTD
K -------""'
[Link]+oo
Nr b
[Link]. SERVO
1
MOTOR TURBINE I-~a.
C:(S) :s'I'a.
GOV_TEMP.MDL
~- SERVO. MOL
5.00e-02
~-
[Link] TNER~rA _ MDLI
.
Fig. 8.14 An XANALOG block diagram for the water-turbine speed-control system simulation. In this
case, the diagram consists entirely of blocks representing submodels. Each of these submodels involves a
further more detailed block diagram based upon combinations of the standard XANALOG block
diagram elements.
some interest in the context of artificial neural sive description of the muscle spindle recep-
networks (e.g. ref. 9). tor or of any other physiological sense organ.
The model which forms the basis of this However, it does embody some elements of
example is extremely simple in form and is widely used physiological models of recep-
not intended to be in any way a comprehen- tors and has properties which are similar in
128 Simple examples using common simulation tools
many respects to those of more detailed phys- of an encoder of this form, including non-
iologically based descriptions. The basic linear effects associated with the response to
equation forming the model is the first-order periodic inputs.
equation Figure 8.15 shows a listing of the DYNAMIC
segment of a SLIM program (ENCODE. SLI)
. vp K (8.10) for this example. This program is included
v =--+-e
P T T on the diskette and provides a further illus-
tration of ways in which the facilities of the
Here the variable v/t) represents the mem- DERIVATIVE section and the DYNAMIC
brane potential, e(t) is the input representing segment can be used to provide quite
mechanical extension of the receptor region complex simulation functions. Figure 8.16
and K and T are constants. When the variable shows some typical results for a constant
vp(t) reaches a threshold level VT the integra- input e(t) and illustrates the 'integrate and
tor is reset to some initial state. When this fire' action of the model very clearly. The
occurs an output pulse is generated, repre- detection of the threshold level and the reset-
senting an action potential in the nerve ting of the integrator is achieved by means of
leading from the receptor to the spinal cord. an IF statement within the DERIVATIVE
The purpose of the model is to allow investi- section and a further IF statement between
gation of some of the fundamental properties the end of the DERIVATIVE section and the
VG=AR*E
C
C ***************Start of Dynamic Segment**************
DYNAMIC
C
C *************Start of Derivative Section*************
DERIVATIVE
VP=INTEG((VG-VP)/TA,VR)
C
C Test whether VP has reached threshold level VT
C
IF(VP-VT)10,20,20
C
C Generate output spike if threshold has been reached.
C If not continue integration process.
C a
20 SP=1.0
10 DERIVATIVE END
C
C *************End of Derivative Section*****************
C
C Values of t, vg, vp and sp for current communication interval
C output to results file
C
TYPE T,VG,VP,SP
C
C Test for end of simulation run
C
IF(T-TMAX)25,50,50
C
C Test for spike. If a spike has been generated the integrator is
C effectively reset by by-passing the DYNAMIC END statement.
C In this case control is passed to line having label 5,
C if end-time not yet reached.
C
25 IF(SP-0.5)30,5,5
30 DYNAMIC END
C
C **************Terminal segment**************************
50 STOP
END
Fig. 8.15 A listing of part of the SLIM program (ENCODE. SLI) for the neural encoder model.
A constrained pendulum simulation 129
......
)0. )
-Ll
/
B t
I
!
~ i!
4~ i
i
~
0'
0 2 4 6 B x1cF t:
a (a)
)a.{ .::,/
b
Fig. 8.16 Membrane potential time histories for two levels of input in the neural encoder simulation
model. These show the 'integrate and fire' characteristics of the system very clearly. Each time the mem-
brane potential (v p) reaches the threshold (10 mV) the integrator is reset (to zero in this case) and an
action potential 'spike' is generated. For the results shown in (b) the input applied was four times the
level used in case of (a). The frequency of the resulting spike train thus depends directly upon the level of
the input applied.
end of the DYNAMIC section. If the integrator facilities for the handling of discontinuities,
output has reached the threshold level an and for the control of integrator operation,
output spike is generated and the integrator which are not available in SLIM. Use of the
is reset by a path which involves bypassing facilities of the DYNAMI C segment to repro-
the DYNAMIC END statement and re-entering duce the operation of the ENCODE. SLI
the DYNAMIC segment. The model can program in other equation-oriented languages
produce interesting and rather complex in exactly the same way may not always be
phase-locking effects when subjected to peri- appropriate. Care should be taken to ensure
odic inputs. Both these situations may be that all of the available features of the lan-
accommodated in the ENCODE. SLI program guage are being used to the best advantage.
with minor editing to allow the sinusoidal
excitation to be used instead of the constant
8.S A CONSTRAINED PENDULUM
input. Figure 8.17 shows an XANALOG block
SIMULATION
diagram for the encoder simulation. Results
in this case are essentially identical to the This example has been proposed by
SLIM results. Breitenecker [11] as an interesting simu-
It should be noted that many equation- lation which can be used for the evaluation
oriented simulation languages have special and comparison of simulation tools. It is
130 Simple examples using common simulation tools
SPIKE
0.0019-+-00
VG
COMPARE
VT
1.00e+0 1
Fig. 8.17 An XANALOG block diagram for the neural encoder model. Note that this diagram includes a
'compare' block which has been implemented as a special-purpose submodel. This XANALOG diagram
is itself a submodel having input variable VG and output variable SPIKE. The block elements associated
with VG and SPIKE are special 'import' and 'export' blocks for use with submodels. Note also the use of
the 'peep' icon, which resembles an eye looking through a circular hole: this provides a means of display-
ing variables within an XANALOG submodel at runtime.
d 2 ({1 g. d d({l
- - = - - SIn ({I - - - (8.11) Fig. 8.18 Diagrams illustrating the constrained
de I m dt pendulum' simulation example.
Simulation of a simple digital control system 131
where d is a damping factor. When the pen- Figure 8.20 shows results from the simula-
dulum reaches an angle 'P = 'Pp its point of tion for both sets of parameters and initial
rotation is at the pin. The length of the pen- conditions given above.
dulum is changed to Is and the angular veloc-
ity d'Pl dt is changed at 'P = 'Pp by a factor Ills'
The form of the equation of motion remains
8.6 SIMULATION OF A SIMPLE DIGITAL
as shown in equation (8.11). When the pen-
CONTROL SYSTEM
dulum swings back so that it again has
length 1the angular velocity is again changed Figure 8.21 is a block diagram of a simple
at 'P = 'Pp' this time by a factor ljI. closed-loop digital control system of the type
Parameters suggested by Breitenecker [11] introduced in section 3.10. The difference
to provide a common basis for software com- between the reference signal and the plant
parisons are as follows: output forms an error, which is sampled by
the analog to digital converter. The digital
m = 1.02 kg 1= 1 m lp= 0.7m processor carries out some form of numerical
g = 9.81 m S-2 operation on the error samples and provides
an actuating signal to the plant input through
Simulation of the motion of the pendulum the digital to analog converter. Figure 8.22
may be considered for a variety of different shows part of the corresponding SLIM
initial conditions, different damping factors program listing for the simplest possible
and for different pin angles 'Pp' Assuming situation in which the control computer
zero initial angular velocity, two sets of con- samples the error signal and outputs the
ditions suggested by Breitenecker are as sampled error values periodically as input to
follows: the plant. The complete program file, named
DIGCON. SLI, and the necessary input data
1. 'Po = 1T16, d = 0.2, 'Pp = -1T112
file DIGCON. IN, can be found on the diskette.
2. 'Po = -1T16, d =0.1, 'Pp = -1T112
Examination of the DIGCON . SLI shows that,
In both cases the response may be considered because SLIM does not have special facilities
over a time period of 10 s. for mixed continuous and discrete system
Figure 8.19 shows part of a program listing simulation, the sampling period is defined as
of a SLIM program (CONS PEND . SLI) on the a multiple of the communication interval
diskette for simulation of this problem. Since parameter, and the facilities of the DYNAMIC
this example involves a discontinuity the segment and DERIVATIVE section are used to
SLIM program makes use of the facilities emulate the discrete action of the controller.
of the DYNAMIC segment to reset the integra- The sampled variable is held constant in a
tion operations at appropriate points. zero-order hold type of action by a loop
The program is therefore similar to the within the DYNAMIC segment. Any discrete
ENCODE. SLI program in some respects. calculations representing the action of a
When the pendulum has reached the critical control algorithm within the control com-
angle and is about to lengthen, or shorten, a puter must also be performed within the
jump occurs out of the DYNAMIC segment. If DYNAMIC segment. Values of variables of
the critical angle has not been reached, opera- interest in the simulation may be written to
tion continues in the normal way within the the output file in the usual way, at times set
DYNAMIC segment and the DERIVATIVE by the communication interval.
section. Comments included within the Figure 8.23 shows results obtained for three
program listing should help to clarify the different sampling periods. It can be seen that
details of this operation. the system shows an increasing tendency to
132 Simple examples using common simulation tools
C Setting flags to indicate phase of motion
IF(XIO-PHIP)4,4,2
2 MARK=-l
MARKl=-l
AL=ALI
GOTO 5
4 MARK=2
MARKl=2
AL=ALI
7 T=Tl
DYNAMIC
DERIVATIVE
DERIVl=X2
DERIV2=-(G/AL)*SIN(Xl)-(D/AM)*X2
Xl=INTEG(DERIV1,XIO)
X2=INTEG(DERIV2,X20)
C Check whether angle of string has reached critical angle PHIP
C If it has not simulation continues without any change of
C parameter AL and without any change of X2 (angular velocity).
C If critical angle has been reached the parameter AL is changed
C (pendulum length) and also variable X2 (angular velocity). The
C changes to be made depend upon the values of PHI immediately
C before the critical angle is reached.
IF(XI-PHIP)20,20,30
C This case applies where the length is ALS (short).
20 MARK=2
GOTO 32
C This case applies where the length is ALI (long).
30 MARK=-l
32 DERIVATIVE END
C output values of variables Xl and X2 to file (and screen)
TYPE T,Xl,X2
IF(T-TMAX)33,33,l2
C Test whether MARK and MARKl flags show that the pendulum
C has reached the critical angle PHIP from either direction
33 IF(MARK-MARKl) 40, 37, 35
C This case applies if pendulum shortens at critical angle.
35 Tl=T
X20=X2*ALI/ALS
MARKl=2
MARK=2
X10=PHIP
AL=ALS
GOTO 45
C This case applies if pendulum lengthens at critical angle.
40 T1=T
X20=X2*ALS/ALI
MARK1=-1
MARK=-1
X10=PHIP
AL=ALI
GO TO 45
37 DYNAMIC END
45 GOTO 7
12 STOP
END
Fig. 8.19 Listing of part of a SLIM program (CONSPEND. SLI) for simulation of the constrained pendulum.
oscillate as the sampling period is increased, The DIGCON. SLI program can be modi-
and eventually it becomes unstable. Sampled fied easily to allow for more complex con-
data theory (e.g. ref. 14) predicts that for this troller action. For example, if one wanted to
system instability occurs when the sampling simulate the system with a controller which
period is greater than 0.549 s. This is consis- implemented a difference equation of the
tent with the results presented in Fig. 8.23, form
and detailed investigations using the simula-
O(kT) = O[(k - l)T] + I(kT) - O.5I[(k - l)T]
tion program can confirm the theoretical
result more precisely. (8.12)
o 2 4 8 t
a (5)
2
1
o
-1
-2
-3
o 2 4 6 8 t
b (5)
o 2 4 6 8 t
c (5)
Fig. 8.20 Responses obtained from CONS PEND. SLI for two sets of initial conditions and different
damping factors. In (a) and (b) the responses shown are the pendulum angular displacement and angular
velocity respectively, for the given set of conditions involving a damping factor of 0.2. Note the discon-
tinuity associated with the string reaching the critical angle and the fact that this is eliminated when the
amplitude of the oscillation becomes sufficiently small. Part (c) is for the second set of conditions, and in
this case the continuous line shows angular displacement of the pendulum while the dashed line is the
angular velocity.
134 Simple examples using common simulation tools
REF.
INPUT
R(s) + I(s) • DIGITAL O(si 4 - OUTP UT
ADC ~ f-+ OAC --
-
PROCESSOR S + 2 YIs)
c
c *******************Start of Dynamic Segment***************
DYNAMIC
1=1+1
C
C *****Test time count for sample instant. Sampling instants
C *****occur at times given by t=MULT*SAM.
C
IF(I)15,25,15
15 IF(I-MULT)20,25,25
C
C *****Take sample of error variable, generate DAC output
C *****and reset sampling timer count to zero.
C
25 AIN=R-Y
OUT=AIN
1=0
C
C *****Write results for current communication interval to
C *****output file (and to screen if required).
C
20 TYPE T,R,AIN,Y
C
C *****************Start of Derivative section**************
DERIVATIVE
DERIV=-(Y-AK*OUT)/TAW
Y=INTEG(DERIV,O.O)
DERIVATIVE END
C
C*******************End of Derivative section*****************
C
C Test for end of simUlation run
C
IF(T-TMAX)10,10,12
10 DYNAMIC END
C
C *********************Terminal Segment**********************
C
12 STOP
END
Fig. 8.22 Listing of part of a SLIM program (DIGCON. SLI) for simulation of the digital control system of
Fig. 8.21. In this case the control computer samples the error signal and outputs the sampled error as the
control input to the plant.
the changes to the simulation program would included on the diskette as the file
all be made in the initial part of the DYNAMIC DIGCONl . [Link] how the discrete input
segment. Figure 8.24 shows the relevant part and output variables are stored for one
of the listing and the complete program is sample period and updated. The transfer
Simulation of a simple digital control system 135
I
I
I
0.8 ~
I
I
L
I
I
I
0.4 , I
\.
----------------------------------------
OL-~
o 2 4 6
______ ~
t
a (s)
~
0.8
0.4
!
-.,I
I -~
2 4 6 t
b (s)
1•
.5
o
I
5 t.. _ _ ...t I
-. t -____________~-L
o 2 4 6 t
c (s)
Fig. 8.23 Responses of the digital control system simulation for three different values of the controller
sampling period. (a) The system output (continuous line) and the error value after sampling by the
analog to digital converter (dashed line) for a sampling period of 0.1 s. (b and c) Equivalent results for
sampling periods of 0.28 sand 0.56 s. Note the instability displayed in the case of (c).
136 Simple examples using common simulation tools
C
C *******************Start of Dynamic Segment***************
DYNAMIC
I=I+I
C
C *****Test time count for sample instant. Sampling instants
C *****occur at times given by t=MULT*SAM.
C
IF(I)15,25,15
15 IF(I-MULT)20,25,25
C
C *****Take sample of error variable and generate DAC output
C *****from present value of controller input AIN and the
C *****values of controller input and output at sampling instant
C *****immediately before the current one (AINI and OUTl).
C
25 AIN=R-Y
OUT=OUTl+AIN-0.5*AINI
c
C *****Re-define previous values AINI and OUTI for use at the
C *****next sampling instant (i.e. the next time around the loop)
C
OUTI=OUT
AINl=AIN
C
C *****Reset the sampling time count variable
C
I=O
C
C *****Write results for current communication interval to
C *****output file (and to screen if required).
C
20 TYPE T,R,AIN,Y
C
C *****************Start of Derivative Section**************
DERIVATIVE
DERIV=-(Y-AK*OUT)/TAW
Y=INTEG(DERIV,O.O)
DERIVATIVE END
C
C*******************End of Derivative Section*****************
C
C Test for end of simulation run
C
IF(T-TMAX)10,10,12
10 DYNAMIC END
C
C *********************Terminal Segment**********************
C
12 STOP
END
Fig. 8.24 Part of the listing of the SLIM program DIGCON1. SLI for simulation of the digital control
system involving the controller transfer function of equation (8.13). With a sample period of 0.347 s this
gives dead beat control characteristics for a step change of reference.
function of the controller of equation (8.12), [14]. In such a situation, the plant output
expressed in terms of z-transforms, is as should exhibit zero steady-state error and
follows: should rise to its final value, in response to a
step change of system reference input, in one
z - 0.5
O(z)
= 1-1-0.5z- sampling period. Figure 8.25 shows results
1
=--- (8.13)
1(z) Z-1 z -1 from the simulation program which are con-
sistent with theory for this special case.
In the special case when the sample period is Some other equation-oriented simulation
0.347 s, the controller should, from sampled languages, such as ACSL, include special facil-
data theory, act as a 'dead beat' compensator ities which can be very useful for the simula-
Simulation of a simple digital control system 137
>-
12
1.
.8
.6
.4
.2
0
0 2 4 6 t
a (5)
...0
......
CD
1. -
.8
.6
.4
.2
°0 2 4 6
t
(5)
b
Fig. 8.25 Responses obtained from the dead beat control system simulation. (a) The output, y, in response
to a unit step change of reference while (b) shows the corresponding sampled error. Note the fact that
this form of controller gives zero steady-state error and that the final value of output is reached after one
sample period.
tion of digital control systems. In ACSL, for This example offers any reader interested
example, DISCRETE sections representing the in automatic control systems many opportu-
difference equations or z-transfer function of a nities for experimentation. It is clear from the
digital controller may be inserted within the listing of the SLIM program DIGCON1. SLI
DYNAMIC segment. Such DISCRETE sections that, with some minor changes to the
are thus similar to DERIVATIVE sections, but DYNAMIC segment, it would be very easy to
communicate with the continuous parts of the replace the dead beat compensator by some
simulation at regular predetermined times. other form of controller. Similarly, any other
Figures 8.26 and 8.27 show XANALOG and form of plant transfer function could be
SIMULINK block diagrams for this digital used in place of the one given in Fig. 8.21,
control problem and illustrate some more of with only some simple changes to the
the specialized blocks and icons available with DERIVATIVE section of the program being
these simulation tools. necessary.
138 Simple examples using common simulation tools
1.00e+OO
1.000+OO -2.00e+OO
Fig. 8.26 XANALOG diagram for simulation of the digital control system.
+ Z-O.5 4 0
Z- 1 S + 2
unit step digital con. plant y
1__-~
Fig. 8.27 SIMULINK diagram for simulation of the digital control system. Note that a discrete transfer
function block in SIMULINK incorporates a built-in sampler at its input and a zero-order hold at its
output. Hence, when discrete blocks are mixed with continuous blocks, as in this example, the outputs of
the discrete blocks are held constant between sampling times.
9.1 THE NEED FOR TESTING The processes of external validation present
many more problems than those of internal
The capabilities and user-friendliness of verification. Although the importance of both
simulation packages have clearly contributed these aspects of testing has been recognized
greatly to the growth of computer-based in some previous texts on modeling (e.g. refs
simulation methods in a very wide range of 1 and 2), it is disturbing to find that most
activities and application areas. This growth applications papers in journals and confer-
means, of course, that simulation techniques ence proceedings pass over questions of
are being used increasingly by nonspecialists external validation in a superficial fashion or
and the risk of misuse is therefore becoming make no mention of it at all. Those appli-
more and more significant. cations studies in which some mention of
Many simulation models are developed to external validation is made seldom provide
provide insight about the dynamic behavior enough detail of the methods used. Ideally
of a complex system, to allow design the processes of internal verification and
decisions to be made, or to provide a basis for external validation should be integral parts of
predictions of future behavior. Simulation the iterative process of model development
results are being used to take important discussed in Chapter 1 and shown in
decisions which may have a bearing on safety diagrammatic form in Fig. 1.1. Therefore, tests
or investment in terms of human and other of validity should be applied many times
resources. It is thus vitally important to be during the development of a model. The
able to test a simulation and to have some documentation associated with the model
objective measures of performance, which can finally adopted for application to the problem
then be incorporated within the documen- in hand need not include every test carried
tation describing the development and out at every stage of development but must
providing information for the user. This include enough detail to allow tests of
involves testing the computer program to external validity to be repeated in full.
establish that it is an accurate implementation Although validation appears to be a
which incorporates the equations of the somewhat neglected topic in most appli-
chosen mathematical model. It also involves cation areas, there are a few safety-critical
testing the mathematical model to ensure applications in which the importance of
that, in the context of the intended model validation is fully accepted. The
application, the mathematical description is aerospace industry provides many good
itself appropriate. These two aspects of illustrations of this since accurate externally
testing involve separate and quite different validated models are needed for the design
activities: the first may be referred to as of aircraft flight-control systems, and for the
'internal verification' and the second as design of simulators for pilot training and
'external validation'. other activities such as investigations of
142 Internal verification and external validation
aircraft handling qualities. The nuclear 'verification' and 'external' with 'validation'
power industry provides other examples in [4]. Internal verification is thus defined as the
which safety requirements make it essential process of proving that a computer simula-
to have a high level of confidence in tion is consistent with the underlying model
simulations. to a specified degree of accuracy, while
The acceptance of model validation external validation involves demonstrating
methods in safety-critical fields means that that the mathematical or conceptual model
there are techniques currently available has an acceptable accuracy over the range of
which can be adopted with benefit in other conditions relevant for the application. These
application areas. Properly validated simu- definitions emphasize clearly that the
lation models, used in an appropriate way, processes which lead to internal verification
can lead to major cost benefits in the design of should provide proof of the internal con-
engineering systems and in other fields of sistency and accuracy of the simulation
application. Inaccurate models, or models program, whereas the assessment of external
applied beyond their range of validity, are of validity of the model involves questions of
little value and can lead to major problems of judgment to a greater extent.
interpretation or understanding and to
engineering design errors, which may prove
9.3 INTERNAL VERIFICATION
difficult and expensive to correct.
Criteria for internal verification involve the
following two aspects:
9.2 TERMINOLOGY AND DEFINITIONS
1. Internal consistency of the simulation
Many different words are commonly used to
program with the mathematical model on
describe aspects of the complete validation
which it is based. The program and the
process. Terms such as assessment, calibra-
model must be shown to involve no
tion, certification, evaluation, qualification,
contradictions in terms of mathematics,
testing, validation and verification are all
logic or concepts.
found. Similarly, words such as accuracy,
2. Algorithmic validity of the simulation
confidence, credibility, fidelity and perform-
program so that all the numerical algo-
ance are used widely in describing the
rithms and associated software routines
qualities of a model. In certain cases these
are shown to be appropriate and provide
words have been given specialized meanings
solutions having a specified numerical
in specific fields of application.
accuracy.
Guidelines on terminology were published
in 1979 by the Technical Committee on Model Internal verification is important at every
Credibility of the Society for Computer stage of model development. All changes
Simulation [3]. Although these have not been within a model must be considered carefully
adopted as a standard by everyone using in terms of internal verification, however
simulation and modeling techniques, the minor they may be. This must involve very
committee has provided very useful careful line-by-line checks of code which
definitions which allow possible ambiguities form the simulation program. In addition,
to be eliminated. One very important there are a number of additional checks and
recommendation was that there should be a comparisons that can be helpful. Checks of
strong distinction made between the words well-understood special cases, which can
'verification' and 'validation'. As indicated in sometimes even be derived using pencil-and-
section 9.1 this has since been extended paper calculations, are particularly appro-
slightly to associate the word 'internal' with priate and these can be divided conveniently
Internal verification 143
into checks of static, or equilibrium, states obtained from the equations suggest that a
and dynamic checks. mistake may have been made in coding the
simulation program. Static checks do not
provide any basis for algorithmic verification
9.3.1 STATIC CHECKS
and are concerned only with questions of
One important example of a check carried out internal consistency.
under static conditions involves a comparison
of equilibrium states in the mathematical
9.3.2 DYNAMIC CHECKS
model and in the simulation. Equilibrium
states are steady-state solutions which can be Dynamic checking, as the term suggests, is a
found by setting all the derivatives on the procedure that goes beyond testing only
left-hand side of the state equations to zero or under static conditions. In general, it involves
by using one of the 'steady state' finders comparing results from the simulation with
which are available in many modern known analytical solutions of the same
simulation software products. Although equations for special cases. Since analytical
algebraic problems may arise in some cases dynamic solutions for nonlinear problems are
with calculations of this kind, the approach is usually not available, this may sometimes
a useful one which can highlight problems require consideration of a slightly modified
that could be much more difficult to identify problem in which nonlinearities have been
under non-steady-state conditions. An temporarily removed for the purposes of
alternative approach to finding steady-state testing. Comparisons of results from a
conditions is to allow the simulation to run simulation program modified in this way
until a steady state is achieved. Whichever with solutions obtained by standard
method of approach is used, comparisons of techniques of linear mathematics, such as the
equilibrium states in the simulation and in Laplace transform, can certainly be of
the real system can also provide a useful basis assistance in locating possible problems both
for external validation. in the structure of the program and in the
It is possible to carry out a partial numerical techniques being used. Dynamic
verification of consistency between the checks therefore address questions of
program code and the mathematical algorithmic accuracy as well as consistency of
formulation by carrying out a static check the simulation code with the underlying
under conditions which do not correspond to mathematical model.
an equilibrium state. This technique was Algorithmic accuracy can be assessed most
introduced first for verification of analog easily in relatively small and simple models.
simulations, but it is equally valid for use Therefore, it is very beneficial to apply
with modern digital simulation tools and is verification of this kind at the submodellevel
especially useful in carrying out checks on and to build up a library of well-tested and
block diagram-oriented simulations. In a documented submodels. These submodels
static check calculation appropriate nonzero can then be incorporated with some
initial conditions are imposed on the confidence in other larger models, provided
integrators, and values of all the other they are appropriate for the intended
variables in the simulation are found. These application.
values are then compared with values In using a well-established simulation
calculated directly from the differential package most of the emphasis in terms of
equations for the same imposed initial state. internal verification is placed on questions of
Any differences in values found from the internal consistency and on the suitability
simulation and the corresponding quantities of options chosen by the user of the package.
144 Internal verification and external validation
The latter may include, for example, 1. Theoretical validity in the sense that the
the integration method, the integration model shows overall consistency with
step length and the communication interval accepted theories or is based upon a
for plotting of output variables. An satisfactory theoretical foundation.
inappropriate choice of anyone of these 2. Empirical validity, with adequate
quantities may cause major problems in agreement shown between the behavior
terms of the eventual interpretation of the of the model and that of the real system
behavior of the simulation. For example, the represented by the model.
choice of an inappropriate integration 3. Pragmatic validity in the sense that the
steplength in a fixed-step integration algo- model satisfies the requirements of the
rithm could cause numerical instability in application.
the simulation program, which could be 4. Heuristic validity in terms of the potential
misinterpreted as a system instability by the of the model for hypothesis testing and
user. In this particular situation, which is a explanation.
common one, simple tests based on changing
the integration steplength, or integration The most important of these criteria, in the
method, can help to establish the true nature context of most practical modeling and simu-
of the instability. If small changes of lation applications, are the ones concerned
integration steplength cause significant with empirical and theoretical validity. These
changes in the model response, the insta- need to be considered at every stage of
bility may well be numerical in nature and development of a model.
be unrelated to any properties of the real
system. Observations of the frequency of the
9.4.1 THEORETICAL VALIDITY
oscillations can also be revealing and a
frequency which is higher than the expected Assessment of theoretical validity involves
bandwidth of the system being modeled checking that the chosen mathematical
may also be indicative of a numerical description does not in any way violate
problem. However, more fundamental important physical laws or principles. It
checks of algorithmic accuracy should not is, for example, very easy to define a
be forgotten, even when using a well- mathematical model in transfer function form
established commercially supported pack- which would incorporate a pure time delay.
age. A responsible user should always give Such a description would be very appropriate
careful consideration to the need for some as part of a model of an engineering process
checks of accuracy, especially when a involving loss-less transport of material in a
particular feature of the package is being pipeline. What comes out of the far end of the
used for the first time by that user or is being pipe is the same as the input but delayed in
employed in an unconventional way. Such time by an amount which depends on the
algorithmic checks normally are of the transport velocity. It would be equally easy to
dynamiC kind. define a mathematical model in which the
delay parameter had a negative value but this
would, of course, be impossible in terms of
9.4 EXTERNAL VALIDATION
the physical system since it would imply that
Criteria for external validation involve the output was available before the input was
assessment of the accuracy and suitability of a applied.
model with respect to the intended Many other situations can arise in which
application. They may include the following mathematical statements have no physical
aspects [5]: relevance, or have restricted validity in terms,
External validation 145
say, of the frequency range in which they are in the real system because of safety
applicable. Situations of this kind must be constraints, but any information available
recognized at the model development stage. about such limits may provide useful insight.
If they involve modeling errors appropriate Comparisons could also be made of the
changes must be made to the structure or frequencies of oscillation of unstable modes
parameters of the model. If they only impose in the real system response and in the model
restrictions in terms of the range of conditions response. Criteria for model acceptability
for which the model will be valid, it may not might then involve agreement to within some
be necessary to make a change in the model given percentage of the relevant system
but decisions must be made about the eigenvalue.
conditions under which the approximation is The third approach involves consideration
acceptable and this must influence the of the response of the system to an
eventual use of the model. Documentation for appropriate input and can be treated in at
the model must then include some statement least two ways:
about the theoretical validity of the model
1. analysis of small perturbation responses
and the range of conditions over which its use
about a selected operating point or about
is permissible.
a series of operating points; and
2. analysis of large-amplitude responses.
9.4.2 EMPIRICAL VALIDITY
Small-perturbation response data may be
Empirical validity is concerned, in the expressed either in the time domain, as a
broadest sense, with comparisons between simple time history, or in the frequency
the behavior of the model and the behavior of domain in the form of separate plots of
the real system. This can include comparison magnitude versus frequency and phase
of chosen system and model variables under versus frequency. Frequency-domain results
equilibrium conditions, comparison of may be obtained from sinusoidal testing or
stability limits and comparisons in terms of from response data generated using other
the dynamic response of selected variables to forms of broadband test input appropriately
chosen input perturbations. transformed using digital signal processing
Evaluation of system and model behavior software. Large-amplitude responses are
under equilibrium conditions can provide normally expressed as time histories as they
much useful insight in the assessment of often involve nonlinear responses for which
external validity, as in the case of internal the frequency domain has no advantage.
verification. Differences in the equilibrium or Responses predicted for a given test input
'trimmed' states of the system and model should, ideally, match the corresponding
may well give indications of problems which system responses to within a given per-
would be much more difficult to identify centage of the full response range in the real
from dynamic responses. Criteria for system. Difficulties can arise with this type of
acceptability of the model could be expressed approach in some cases if the system and
in terms of the percentage difference between model responses are very sensitive to small
appropriate variables of the model and differences of initial condition or of input
system at each equilibrium point considered. amplitude. Such small initial differences may
Stability boundaries for the simulation become exaggerated as time increases because
model should correspond to the stability of integration effects within the system, and
limits for the real physical system. Evaluation direct comparisons of time histories over an
of these limits by actual testing and extended period following the application of
experimentation may not always be possible a test disturbance can therefore sometimes be
146 Internal verification and external validation
of limited value. In the case of a single output can be very helpful as a first step towards
variable, criteria based upon the integral of external validation. For example, the absence
the error squared or the integral of the of an oscillatory response in a model for a
absolute error between the system response system which exhibits strongly oscillatory
and the model response can be of value. A behavior provides an immediate clue which
more general criterion function, which can be can be useful in reconsidering the structure
useful for the assessment of models using and parameter values previously chosen.
direct time history comparisons in the case of Models are not unique, and in most
a model with a number of output variables, situations of practical importance there will
has the form be a number of candidate models which give
an adequate match for a given set of
f
J = eT (t)W(t)e(t)dt (9.1) experimental responses. Models should
therefore be assessed for a variety of
experimental data sets and, before being
where e(t) is the vector representing the time accepted for the intended application, should
response error function and W(t) is a matrix be shown to be capable of matching the
which provides appropriate weighting. experimental results to an acceptable level of
Optimization of criterion functions such as accuracy for all of the test conditions.
those discussed above provides one approach Methods of parameter sensitivity analysis
to the refinement of models in which can help the modeler considerably in gaining
discrepancies have been found. However, an understanding of the potential of a chosen
care must be taken to consider the possibility mathematical description to describe fully the
of errors in model structure as well as errors behavior of a system under both steady-state
in model parameters when carrying out this and dynamic conditions. However, detailed
type of optimization. Unrealistic parameter sensitivity analysis for a complex model is a
values in a model are often indicative of a computationally intensive task and can also
deficiency in the model structure, but it must provide problems in terms of presentation
be recognized that such structural deficiencies and efficient interpretation of results. The
may be difficult to correct. Once corrected, use of parallel-processing techniques and
however, the model structure should be improvements to the user interface tradi-
capable of providing a good fit to system tionally used in simulation work could well
response data over an expanded range of produce benefits in this respect in the future.
conditions. Sometimes, lack of appropriate experi-
Models can also be assessed using particu- mental response data from the real system
lar features of a time history or a small part of may make it impossible to compare model
the total frequency range rather than the responses with measured data. This can arise,
complete record. For example, quantitative for example, in engineering design appli-
comparisons may be made of particular cations when a simulation model may be
features of a time history, such as the shape developed at an early stage in a new project
of a response overshoot for a step input test, before any real hardware or prototype system
or of a frequency response in the vicinity of a exists. Very often, in such circumstances, it
resonant peak. Comparisons may also be may be possible to use previous modeling
useful in terms of individual features of a experience from the appropriate application
response, such as the actual frequency of an area and validated submodels to minimize
oscillatory response or the magnitude or time risks and uncertainties inherent in the new
of occurrence of an overshoot. Even quali- model. It may also be possible to assess the
tative comparisons involving such features overall model results by making comparisons
External validation 147
with data from other validated simulations, validation. Such submodels may be available
even if these other representations are not from previous design studies involving
appropriate for the intended application. similar systems. Complete external validation
Good examples of this kind are often to be of a model, initially validated only at
found in the development of real-time subsystem level, is clearly desirable as soon
simulation models where a simplified as data becomes available from the full
description, capable of running in a fast system.
timescale on a given computer, may be One important, but often neglected, aspect
assessed by making detailed comparisons of the processes of external validation is.
with a much slower but more detailed model concerned with the specification of tests for
which has previously been subjected to a full the generation of time history data for
validation process. This is, of course, a validation purposes. Often experiments have
weaker form of external model validation to be designed specially for model testing and
than is possible by the direct approach, and it may involve monitOring of quantities which
must be used with care. Careful study must are not normally accessible for measurement.
be made of the model which is being used as In engineering applications this may involve
reference, especially in terms of the intended significant augmentation of the standard
area of application, the inherent assumptions instrumentation for the system. There is, in
and the recorded limitations. most cases, an inverse relationship between
the initial level of understanding of a given
system and the complexity of the instrumen-
9.4.3 MEASURED RESPONSE INFORMATION
tation needed for a comprehensive empirical
FOR EMPIRICAL VALIDATION
validation of a model for that system.
Empirical validation requires the availability The information content of the experi-
of data describing the behavior of the real mental response data used for model
system which is being modeled. Such data validation must be maximized, subject to
may involve one or more of the following: constraints imposed by safety. Careful con-
sideration must be given to system band-
1. measured time histories from tests
width, the frequency content of test signals,
carried out on isolated subsystems;
sampling rates, signal amplitudes and signal
2. measured time histories from tests
to noise ratios. Factors such as these are also
carried out on the complete system; and
important in system identification methods as
3. the opinion of an independent person
discussed in section 9.4.4.
who has expert knowledge and practical
One important source of validation
experience of the system in question.
information which should never be neglected
The use of measured response data from is the independent expert. Much useful
tests carried out on isolated subsystems is insight can be provided by someone who is
especially valuable for engineering appli- familiar with the real system, and simple tests
cations involving the evaluation of competing which can be very revealing are often
design options through the use of simulation suggested by such experts. Those responsible
models. As pointed out in the previous for the development of simulator systems for
section, measured response information for training purposes have found that valuable
the complete system under investigation information can be gained very quickly and
clearly cannot exist at the design stage, but effectively in this way, although there may
use may be made of subsystem models which well be differences in the way such infor-
have already been subjected to rigorous mation is used in correcting deficiencies in a
processes of internal verification and external training simulator and in a mathematical
148 Internal verification and external validation
model being validated for other purposes. In are mostly concerned with the estimation of
the case of a training simulator it may well be parameters for small perturbation models for
acceptable to make empirical adjustments particular operating conditions. From meas-
until the experienced operator is satisfied that ured system response data, identification and
the simulation is realistic. In a model being parameter estimation tools can provide para-
developed for other purposes a more thor- meter estimates, together with associated
ough understanding may be needed of all the estimates of variance for the parameter values.
factors underlying the alleged deficiency. Potential problems in the use of system
identification methods for model validation
can be assessed, for a given model structure
9.4.4 SYSTEM IDENTIFICATION AND
and available set of measured variables, by
PARAMETER ESTIMATION METHODS IN
carrying out a preliminary analysis of
THE EMPIRICAL VALIDATION PROCESS
identifiability. In models which are theor-
The factors which are important in the design etically identifiable all of the parameters may
of experiments for empirical validation are be estimated uniquely given ideal test data. In
essentially the same as the factors which models which are theoretically unidentifiable
influence experiment design for system some parameters may be estimated uniquely,
identification. System identification may be whereas others could have a range of possible
defined as a process of determining, on the values consistent with the test data. Many
basis of measured input and response data, a different techniques are available for the
model to which the system under test is assessment of identifiability [7], and some of
equivalent. these are applied in the case studies in later
Some system identification techniques, chapters.
which have been developed largely in the System identification and parameter
context of engineering applications involving estimation methods provide a particularly
automatic control, are often closely linked to important form of external validation in fields
the development of models of data and may of application in which linearized dynamic
not involve models which have any models provide a natural form of system
underlying theoretical structure. In this respect description. One important example comes
system identification methods are very similar from the aerospace industry, in which
to the time-series modeling techniques used by linearized state-space models of aircraft
statisticians and operations researchers, and provide the mathematical descriptions
indeed the underlying theoretical foundations needed for some aspects of flight control
are very closely related. These techniques are system design and handling qualities studies.
now being considered increasingly important In this particular type of application, system
for the development and testing of dynamic identification techniques have also been
models used in continuous system simulation found to be of value in the assessment of
[6]. Frequently, insufficient information is nonlinear models. At each chosen operating
available to allow a mathematical model to be condition the parameters of a linearized
developed entirely from theory, and infor- model obtained by system identification can
mation must be extracted about parameter be compared with the corresponding
values or possible model structures from theoretical parameters found by linearizing
observations on the real system. This process the underlying nonlinear model derived
of inverse modeling from measurements is the using physical principles. Comparison of the
essence of system identification. values of the theoretical and estimated
Aspects of system identification which are parameters for a number of different oper-
of particular relevance for model validation ating points can be revealing. The tendency of
Robustness issues in external validation 149
an estimated parameter to increase or further consideration of it goes beyond the
decrease in value over a range of conditions scope of this book.
involving a series of equilibrium points can
provide considerable insight since a similar
9.5ROBUSTNESS ISSUES IN EXTERNAL
trend in the estimated and theoretical
VALIDATION
parameters tends to support the validity of
the theoretical nonlinear model over the Empirical validity is, of course, one of the
operating range involved [8]. In the case of most important aspects of the external
the aircraft model this might involve validation process, and the robustness of a
comparing theoretical and estimated para- model which has been subjected to the
meters for a [Link] of small perturbation tests processes of empirical validation must be
for a series of trimmed flight conditions at demonstrated if the model is to be credible. In
different forward velocities. general, robustness of a model which has
The importance of system identification been subjected to external validation may be
methods for model validation in aircraft checked further by carrying out tests using
applications has been given further emphasis data sets not employed at any stage in the
recently in a report concerned with the development and refinement of the model. If
application of identification techniques to the differences between the model predictions
helicopters and other forms of rotorcraft [9]. and the corresponding measured variables
This report stresses the need for accurate are all sufficiently small the model may be
mathematical models for the development of accepted as a candidate for the intended
high-performance flight-control systems and application. If the model does not satisfy this
suggests that system identification techniques requirement it must be reconsidered and all
may become mandatory in the future for theory and data used at any stage of model
model validation in ground-based flight development must be reassessed. However, it
simulators. should be noted that even if the model gives a
Although system identification is not good prediction for all the additional tests
carried out exclusively for model validation considered it will not be a unique repre-
purposes, the factors which are important for sentation and other models could give similar
obtaining reliable and robust estimates of results.
parameters in an identified model have The data sets to be used for assessment of
an important bearing on the successful model robustness should be similar in
validation of a given model. As has already character to data sets used at the model
been pointed out, the success of the empirical development or refinement stages in terms of
validation process is closely linked to the the amplitude and frequency distribution of
design of experiments used to test the model. inputs and responses. However, results
Similarly, in system identification and obtained using a number of forms of test input
parameter estimation, the reliability of may give differences which are under-
estimates is critically dependent upon standable in terms of physical reasoning and
experimental design. In both cases the may provide useful insight, which can lead to
information content of the experimental further model improvements. Additional tests
response data must be maximized subject to or new experiments or observations may also
practical constraints. Much has been written be needed in order to resolve difficulties
about system identification techniques and identified at this stage.
about methods of experiment design for Robustness issues are also of considerable
system identification (e.g. refs 6, 7 and 10). importance for the parameter estimates which
This is a specialist topic in its own right and are incorporated into a model at the initial
150 Internal verification and external validation
development stage. This is really inseparable model leads to three possible outcomes.
from the question of robustness, in terms of These are as follows:
model validation, since they both form part of
1. The data sets used for model testing
a larger iterative process. However, it should
cannot be explained by any combination
be noted that the variance of each parameter
of model parameter values and model
estimate provides a useful measure of
structures considered. The candidate
the relative reliability. Whenever system
models have therefore all been falsified. It
identification methods are used as part of the
is necessary to reconsider the model
initial development stage, careful checks
structure and assumptions once again
should be made of the variance values
from first prinCiples.
provided by the chosen parameter estimation
2. One or more models produce a satis-
method, and residuals (the differences at each
factory match to the available system data
point in time between [Link] system
but the uncertainty level is excessive for
quantities and the correspondmg model
certain parameters. Such models may not
variables) should be examined to assess the
be of much predictive value and re-
possible presence of systematic erro:s caused
formulation of the model may again be
by an incorrect model structure .whIch co~ld
necessary.
show up as parametric dIscrepanCIes.
3. A model gives satisfactory agreement
Repetitions of the estimation process. for a
with experimental test results from the
variety of different record lengths or, m the
real system and involves parameter values
case of frequency-domain identification
which are considered plaUSible. In this
methods, for a number of different frequency
case the model is of acceptable validity in
ranges can be useful. Estimated parameters
terms of the particular tests performed.
should show low sensitivity to record length
The model may be used for the intended
and frequency range used.
application until such time as evidence is
In general it is important that identification
found which falsifies the model.
tools which are being used for the
development and validation of con~ius
dynamic models for simulation app!Icattons 9.7 DOCUMENTAnON OF THE
should incorporate a comprehenSIve and VALIDA nON PROCESS
well-engineered user interface which .exploits
When carrying out the complex and time-
the maximum use of simple graphICS. The
consuming tasks of internal verification and
tools available should provide information on
external validation it is essential that the
goodness of fit, confidence intervals f?r
complete process is documented in detail and
estimated quantities, sensitivity to changes ~n
is thus repeatable. The SCS Technical
model structure and sensitivity to changes m
Committee on Model Credibility recom-
test condition [11]. Good facilities are also
mends [3] that the documentation associated
needed for management of the very large
with the assessment of credibility and
amounts of data required for identification,
applicability of a model should include at
and for comprehensive documentation of the
least four basic elements as follows:
results.
1. a clear statement of the purpose of the
model;
9.6 POSSIBLE OUTCOMES OF THE
2. descriptions of the model in both con-
EXTERNAL VALIDA nON PROCESS
ceptual and mathematical terms;
Application of the methods outlined above 3. a statement concerning the range of
for assessing external validity for a given conditions for which the model has been
References 151
tested and the range for which there is 4. Murray-Smith, D.J. (1990) A review of
overall agreement with the real system methods for the validation of continuous
behavior; and system simulation models, in Proceedings of
the 1990 UKSC Conference on Computer
4. a description of all tests used for internal
Simulation (ed. K.G. Nock), United Kingdom
verification and external validation, Simulation Council, Burgess Hill, pp. 108-11.
together with relevant results and a 5. Murray-Smith, D.J. and Carson, E.R (1988)
justification of the tests chosen. The modelling process in respiratory
medicine, in The Respiratory System (eds D.G.
Model documentation of this kind provides Cramp and E.R. Carson), Croom Helm,
a basis for formal certification of a model in London, pp. 296-333.
application areas where this is required. A 6. Unbehauen, H. and Rao, G.P. (1987)
potential user of a model documented in this Identification of Continuous Systems, North-
way can review the evidence provided, Holland, Amsterdam.
consider the context in which it was 7. Beck, J.V. and Arnold, K.J. (1977) Parameter
Estimation in Engineering and Science, John
developed, assess its suitability for the
Wiley, New York.
intended application and take note of stated 8. Bradley, R, Padfield, G.D., Murray-Smith, D.J.
limitations. and Thomson, D.G. (1990) Validation of
helicopter mathematical models. Transactions
of the Institute of Measurement and Control, 12,
186-96.
9. AGARD (1991) Rotorcraft System Identification,
REFERENCES
Advisory Report No. 280, AGARD, Neuilly
1. Shannon, RE. (1975) System Simulation. The sur Seine.
Art and the Science, Prentice-Hall, Englewood 10. Sinha, N.K. and Kuszta, B. (1983) Modeling
Cliffs, NJ. and Identification of Dynamic Systems, Van
2. Spriet, J.A. and Vansteenkiste, G.c. (1982) Nostrand Reinhold, New York.
Computer-aided Modelling and Simulation, 11. Murray-Smith, D.J. (1991) Modelling and
Academic Press, London. Robustness Issues in Rotorcraft System
3. SCS Technical Committee on Model Identification, AGARD Lecture Series No. 178,
Credibility (1979) Terminology for model Rotorcraft System Identification, AGARD,
credibility. Simulation, 32, 103-104. Neuilly sur Seine.
CASE STUDY I - A TWO-TANK LIQUID
LEVEL CONTROL SYSTEM
10
(10.1)
H
M=pV (10.2)
pump
H2
1 av1 avo
H -!.~ _._ .... -. '_. __ .-... -.-.-.-.-.-.- ._- ~
through the insertion of plugs into one or is the cross-sectional area. Similarly for tank 2
more of these holes. The system is equipped we can write
with a drain tap, under manual control, and
the flow rate from one of the tanks can be A 2 dH 2 -Q -Q
dt - v1 vo
00.13)
adjusted through this. The other tank has an
inflow provided by a variable-speed pump,
which is electrically driven. Both tanks are where H2 is the height of liquid in tank 2 and
equipped with sensors which measure the Qvo is the flow rate of liquid out of tank 2.
pressure at the base of each tank and thus Considering the holes connecting the two
provide an electrical output voltage propor- tanks and the drain tap all as simple orifices
tional to the liquid level. allows the flow rates to be related to the
liquid heights by the following two equations
10.2.1 A NONLINEAR MATHEMATICAL
00.14)
MODEL
Following the approach used in section 10.1 and
the equation describing tank 1 in Fig. 10.2 has
the form Qvo = Cd ,a2 (2g(H 2 - H 3 )t 2
00.15)
(10.16)
Reorganization of equations (10.22) and
(10.23) gives a second-order state-space
(10.17) model as follows
(10.18)
(10.27)
(10.20)
and
Similarly
·
W ... -.
o
TANKS _ ALG
'.'
.'
"'iI]Ga~1.
H3
3_0009-02
H2 L.......... Gao
~ .....-_H3~ GaD
1_0309-+-02 INT2
~
W
.' :"
TANKS_ALG
.. •
I
.......•••••..•.••.....•••••••.••••...••..•..................••.......•.......... I
Fig. 10.4 XANALOG block diagram for simulation of the two tank system. Note the use of submodels
and their associated icons.
parameter set are shown in Fig. 10.5. These internal verification is concerned with check-
show the changes of depth h1 and h2 versus ing that the structure of the simulation
time for a number of different initial condi- program is consistent with the mathematical
tions. Are these results meaningful? Is the model. This involves working backward from
mathematical model adequate and does the the statements in the program, especially
simulation program represent the model to a those within the derivative section, to ensure
sufficient degree of accuracy? In order to that when translated back to the form of dif-
answer these questions in a satisfactory way ferential equations they are the same as those
one must first consider carefully how this of the original model. Checks should also be
simulation can be verified and how the model made of the parameter values used in the
can be validated. program or in the parameter input file to
ensure that they correspond exactly to the
parameter set of the model itself.
10.3.1 INTERNAL VERIFICATION OF THE
The second stage of internal verification is
SIMULATION PROGRAM
concerned with numerical accuracy. In the
The simulation program for the coupled-tank case of fixed-step integration methods, com-
system is a very simple one. The first stage of parisons can be made of results obtained with
Simulation of the nonlinear coupled-tank system 159
0.06
~I r-~_,
0.18 ~
0.14 f.
~
0.1 0 ~ ,- - , - - "__
L/ . . -- . . - . . - . . __
-------------..,
-~ -____________________________________
0.06 L
------- ~
~
0.18
Fig. 10.5 Simulation results (from SLIM) for the two-tank system model for three sets of initial condi-
tions. The input flow rate was the same in all the cases, with a value of 16.67E-06 m 3 S-I (1000 cm3 min-I).
The continuous line represents the depth HI and the dashed line denotes H 2 • Note particularly the case
where the two tanks are both initially empty. Can you explain the behavior of the system, on physical
grounds, during the first part of the response?
a number of different sizes of integration ent settings of the relative and absolute error
step length and with different integration limits and with different values of the
techniques. This provides the user with some minimum integration step to be allowed. In
understanding of the sensitivity of results to both cases, some comparisons can be made
the step length and of the overall suitability of using a number of different values of the
the numerical methods chosen. In the case of communication interval to ensure that inter-
variable-step integration algorithms, tests can esting events in the simulation model are not
be carried out to compare results with differ- being hidden from the user simply because of
160 Case study I
an inappropriate choice of this parameter to checking a mathematical model which can
which determines the interval between provide a basis for any definitive statement
output samples. about the overall validity of that model.
Statements about model validity must be
made in the context of an intended applica-
10.3.2 EXTERNAL VALIDATION OF THE
tion. In the case of the coupled-tank system,
SIMULATION MODEL
the computer simulation is to be used as a
The discussion of validation in Chapter 9 basis for the design of an automatic control
shows clearly that there is no single approach system which will ensure that a given level is
Q, H, H, H, H,
measured model measured model
(em 3min-') (em) (em) (em) (em)
C\j-
IE
£.
21
18
15
12
o~-+
o 100 200 300 400 500
__
600
t
(s)
Fig. 10.6 Simulation results for H, (dashed line) and the corresponding measured response of the real
system for an experiment involving doubling of the input flow rate from Qvi =1000 cm3 min-' to QVi =2000
cm3 min-' at time t = 50 s. Note the initial steady-state difference between the measured and simulated
responses.
Simulation of the nonlinear coupled-tank system 161
maintained in one of the tanks. There is par- change in operating conditions. Figure 10.7
ticular interest therefore in the accuracy of the shows measured response data for small-
model in predicting steady-state conditions perturbation step tests carried out about one
and in predicting the form of small transients chosen operating point. Time constants esti-
about any given steady operating point. Such mated from measured step responses such as
comparisons are very easily made in the case these, from which the error plots of Fig. 10.7
of small-scale laboratory equipment of this were derived, can also be compared with
kind, and agreement between steady-state values determined from the linearized model
measurements and steady-state model predic- in the form of equations (10.27) and (10.28).
tions is generally quite good for upper parts The discrepancies in the model exposed by
of the operating range. Table 10.2 shows some the steady-state tests, and the large perturba-
typical results obtained from measurements tion responses, are believed to arise mainly
on the real system and corresponding tests on because of the limitations of equations (10.14)
the simulation model. Differences between and (10.15) in describing the relationships
the steady-state liquid levels in the simulation between output flow and the liquid level in
model and in the real system, for smaller each tank. These equations apply to an ideal
values of input flow rate, are significant and simple orifice and the actual physical effects
vary slightly with operating point. Figures at the tank outlets do not agree exactly with
10.6 and 10.7 show some comparisons in this simplified model.
terms of dynamic tests. In Fig. 10.6 discrepan- With closed-loop control added to the real
cies between the model and system results system, and to the simulation model, the
are shown for a test involving quite a large agreement can be shown to be significantly
~NE
<J I,)
4
3
0 500 1000
sec
~ E
I,)
3
2
o 500 1000 t
sec
Fig. 10.7 Plots showing differences between model predictions and measured values for HI and H, for
small perturbation step tests carried out at one chosen operating point. The step input was applied at time
t =250 s. Note the initial steady-state errors between the simulation model variables and the corresponding
measurements and the increased steady-state error value for H, after the transient has died away.
162 Case study I
closer. This is important since the equipment tively simple nature of the system, and the
is intended to be used for investigations of variables which are accessible for measure-
closed-loop control. In simulation studies ment in the real hardware, make this
involving control system design applications an interesting but straightforward system
there is always particular interest in the for the application of external validation
overall robustness of the control system and methods.
the effect which modeling errors and uncer- Possibilities for using the simulation
tainties may have on the performance of the program TANKS. SLI as a basis for further
closed-loop system. Although control systems investigations on your own fall into two main
are normally designed using linearized areas. One obvious topic to consider would
models, simula,tion studies carried out on a involve investigation of the effect of changing
proposed closed-loop system using a non- the form of relationship used to describe the
linear model of the plant can often be highly discharge nonlinearities. A second area for
illuminating. Such an investigation may independent study would involve adding
reveal problems with the proposed design feedback control. Control system design
which would otherwise only come to light studies could be carried out using the linear-
during the commissioning and testing of the ized form of the mathematical model and dif-
real system. ferent controllers could then be compared in
terms of their sensitivity to changes in operat-
ing points or to changes in plant parameters,
10.4 DISCUSSION
such as the coefficients of discharge.
This case study provides an illustration of a
relatively simple nonlinear system which can
be modeled in a classical way using physical REFERENCE
laws and principles. The simulation model is 1. Wellstead, P.E. (1981) Coupled Tanks Apparatus:
easily implemented using either equation- Manual, TecQuipment International, Long
oriented or block-oriented tools. The rela- Eaton, Nottingham.
CASE STUDY II - AN AIRCRAFT
AUTOMATIC LANDING SYSTEM 11
longitudinal
lateral
axis
rudder
/
______-"7 elevators
Fig. 11.1 The longitudinal (roll), lateral (pitch) and normal (yaw) axes of a conventional fixed-wing air-
craft, together with the control surfaces which are used to apply control forces and moments to the
vehicle.
I
,~ .... ----- -
l' ------ . runway
centre
line __ ,7-~90_HZ=:JIiT<::::':' ______ ~ ••.transmitter
( 150 Hz ___
" ------
Fig. 11.2 Radio transmission patterns associated with the localizer transmitter within an instrument
landing system (ILS).
Simulation of an aircraft directional-control system 16S
ISO-Hz signal will be greater. The received tion to investigate a nonlinear closed-loop
signals provide the input for the indicator system model which could not be studied
display device and can also provide a refer- using analytical methods alone. No consider-
ence for the directional control system. The ation has been given to the question of model
linking of the ILS receiver output and the validation in this particular study. In an
autopilot reference input is carried out using application such as this, in which a simplified
a 'coupler' unit. model is used for the analysis and design of a
The glide-path transmitter is located near flight-control system, validation might well
the point of touchdown on the runway and be carried out by comparing the performance
transmits on a given frequency in the band of the reduced model with a more versatile
329.3 MHz and 33S.0 MHz. The radiated and complex model, which had been the
signal pattern is similar to that of the localizer subject of separate investigations in terms of
transmitter. When on the correct descent path model validity. An incremental approach has
the aircraft receives the 90-Hz and ISO-Hz been adopted in which the simulation model
signals with equal strengths and the hori- of the lateral beam guidance system is built
zontal pointer takes up a central position on up from a simpler linear model of a direc-
the ILS display device. In this case the ILS tional-control system. This allows verification
receiver can be linked to the longitudinal atti- to be carried out through comparisons with
tude-control system of the aircraft and pro- results from linear systems theory.
vides information which ensures that the
aircraft follows the glide path.
11.2 MODELING AND SIMULATION OF
This chapter is concerned with a case study
AN AIRCRAFT DIRECTIONAL-CONTROL
which involves the lateral beam guidance
SYSTEM
system, which uses the localizer signal to
align the aircraft with the runway. The In general, an uncontrolled aircraft has no
second major component of the landing inherent tendency to return to its initial
system, which ensures that the aircraft heading and roll attitude after some distur-
descends on the correct glide slope, is not bance from equilibrium. Thus, in the absence
considered in this simulation study. The two of any kind of flight-control system, the pilot
systems are very similar in terms of their must make corrections continually in order to
block diagram structures and, since they are maintain a required heading and attitude.
essentially independent in their operation, Directional-control systems can provide a
they can be investigated separately. In devel- means of ensuring that a reference heading is
oping a mathematical model for the lateral maintained without pilot intervention.
guidance system it is necessary first to One particular form of directional-control
describe the aircraft's directional control system will be examined in detail. In this
system. Preliminary investigation of this sub- case the system being controlled is the coor-
system, which is of central importance for dinated aircraft, which incorporates coup-
lateral guidance, provides an interesting ilJus- ling between the aileron and rudder. Such
tration of simulation methods for control coupling allows coordinated maneuvers to
system studies. be made, and only the overall transfer func-
Aircraft models are well understood and tion between the aileron input and the roll
well documented. Therefore, in this particu- angle need be considered. A block diagram
lar study, emphasis has been placed on the for the control of heading is shown in Fig.
process of building up a simulation program 11.3. It may be seen from this that the system
to represent a system provided in block involves three feedback loops. The outer-
diagram form and on the use of the simula- most loop involves the feedback of the
166 Case study II
variable of primary concern, which is the where KA is a gain factor and TA is a time con-
heading angle. This feedback is provided stant associated with this highly simplified
through the action of a directional gyro- description of aircraft lateral dynamics.
scope. The other two loops provide sec- Similarly, in terms of the roll angle, ~, for
ondary feedback signals which are needed to an aileron deflection angle, oa' we have a rela-
ensure that the overall system has an accept- tionship of the form
able dynamic response. The block diagram
shows that the second pathway involves roll KA
(11.2)
angle feedback, via the vertical gyroscope, s(1+sTA )
while the innermost loop has roll rate feed-
The control-surface actuation system can be
back. It should be noted that the aircraft
described approximately by a transfer function
description adopted here is a highly simpli-
fied one, and that a detailed description of a 1
(11.3)
practical flight-control system would involve 1+57
a significantly more complex block diagram.
Further details of aircraft flight mechanics where 7 is a time constant.
modeling and of flight-control systems may The three transfer functions given above,
be found, for example, in the texts by and which appear in the block diagram, may
Blakelock [1] and McLean [2]. The reference be converted readily to a state-space form of
quantity is the heading angle commanded model using the methods discussed in Chapter
by the pilot or provided as an electronic 3. This state-space description has the form
signal from some other system within the shown in equations (11.4), (11.5) and (11.6).
aircraft. The various gain constants associ-
dp 1 KA
ated with the feedback pathways in this -=--p+-o (11.4)
system (K d , Kv and K) are chosen at the
dt TA TA a
-------------------------
, /
Fig. 11.5 Results from a simulation experiment on the aircraft directional-control system showing time
histories of the heading angle (dashed line) and the roll angle (continuous line) in response to a
demanded step change of heading of 0.15 rad. The value of the directional gyroscope gain factor (Kd ) in
this case was 2.0, while the parameters K/ and Kv were both 2.5.
~
-s: x10~r-
6
o~-
o 2 4 6 8 10 12 14
t
Is)
Fig. 11.6 Results obtained from simulation of step test on directional-control system for a directional
gyro gain factor of 0.5. Values of the other parameters read in from the input data file were the same as
for the case given in Fig. 11.5.
gyroscope gain factor (0.5) produces a very in Fig. 11.7, where a value of Kd of 8.0 gives a
sluggish response (Fig. 11.6) in which the highly oscillatory response, which is again
required heading is not attained within the unsatisfactory for a system of this kind.
15-s duration of the simulation experiment. Although the mathematical tools of linear
The performance, in this case, is clearly in- systems analysis could be applied very
ferior to that of the system with the higher readily to a linear system such as this, the use
value of gain. However, increasing the gain of a simulation can provide opportunities for
factor can produce other problems, as shown additional understanding, particularly in
Simulation of a lateral beam guidance system 169
~
-s..
X101
6
4
2
0
-2
-4
-6
0 2 4 6 8 10 12 14
t
(8)
Fig. 11.7 Results similar to those of Figs 11.5 and 11.6 for a directional gyro gain factor of 8.0.
terms of parameter sensitivity analysis of the Figure 11.8 is a block diagram of the lateral
kind illustrated by these results. Linear beam guidance system which must be inter-
theory also provides opportunities to verify preted in conjunction with Fig. 11.9. This
the simulation program. In this case, for shows the geometry of the situation and
example, the simulation could be used very defines some of the variables used in the
readily to establish stability limits in the block diagram.
directional-control system and selected sets The variable «/Ifefl which represents the
of results from the simulation could be com- heading to be followed, would be set by the
pared with results obtained by applying pilot and, for simplicity, this reference value
standard stability tests from control theory, is taken to be zero in the equations which
such as Nyquist's criterion. follow. An error angle, A, is related to the
lateral displacement, y, and the range, R, by
11.3 MODELING AND SIMULA nON OF A the equation
LATERAL BEAM GUIDANCE SYSTEM
(11.11)
The lateral beam guidance system is based
upon the directional-control system described Hence, for small angles
in section 11.2. In the lateral guidance
application the directional reference input, (11.12)
«/Ie' is derived from the ILS receiver output.
aircraft
&
lateral
auto ilot
,,
,"
R -~ ",,"
L.. "'aircraft
r'-"=_:~7.;?tC1 ___-,.~ ___ centre
line
Also, it follows from the geometry of Fig. 11.9 where Gc is the gain of the coupler.
that the rate of change of displacement, Combining the block diagram of Fig. 11.8
dy / dt, is related to the aircraft velocity, Uo' by with the block diagram of the directional
the equation flight-control system gives the complete
structure of the lateral beam guidance system,
dy =Uo sin(1/I - I/Iref ) (11.13) as shown in Fig. 11.10.
dt The parameter set for the lateral beam gUid-
Thus, for a small angle 1/1 and for a zero value ance system includes all of the parameters of
of reference heading (I/Ir)' it follows that the directional flight-control system of section
11.2 and, in addition, the coupler gain para-
dy 0:: U ./, (11.14) meter values for the outer loop of Gc = 8.0. In
dt 0."
some cases the coupler involves a transfer
The signal from the localizer course correc- function of a more complex form, incorporat-
tion receiver is approximately proportional in ing proportional plus integral control in a
magnitude to the angle A and, for the case transfer function of the form
where Are! is zero (see Fig. 11.8), the output of
the coupler is given by Gc(s) = G{1 + ~I ) (11.16)
COUPLER
+
Fig. 11.10 Block diagram showing the detailed structure of the lateral beam guidance system.
Simulation of a lateral beam guidance system 171
in many other control system applications. 1t' DYNAMIC
DERIVATIVE
is known as a proportional plus integral con- YD=UO*PSI
troller and can ensure zero steady-state errors Y=INTEG(-YD,YO)
R=INTEG(-UO,RO)
in the presence of a steady external distur- ALAMB=Y/R
PSIC=GC*ALAMB
bance, such as a cross-wind in this instance. PSID=PHI*G/UO
The parameter K] is taken as zero in the exam- PSI=INTEG(PSID,PSIO)
PHIC=(PSIC-PSI)*AKD
ples considered here. EV=PHIC-PHI
Figure 11.11 shows the DYNAMI C segment ES=EV*AKV-PHID*AKR
DAD=(ES-DA)/TAW
of a SLIM program for simulation of the DA=INTEG(DAD,O.O)
Q=(AKA*DA-PHID)/TA
lateral beam guidance system. The complete PHID=INTEG(Q,PHIDO)
program (AILS. SLI) is provided on the PHI=INTEG(PHID,PHIO)
DERIVATIVE END
diskette. The parameter values for the aircraft RP=R/I000.0
TYPE T,PHI,PSI,Y,RP
and directional-control system are assigned in IF(T-90.0)10,10,12
the initial region together with the coupler 10 DYNAMIC END
12 STOP
gain parameter which is input through an END
ACCEPT statement along with the aircraft
Fig. 11.11 Listing of the DYNAMIC and TERMINAL
velocity. Once again, initial conditions, in segments of a SLIM program (AILS. SLI) for sim-
terms of the initial range, initial deviation ulation of the lateral beam guidance system.
from the runway center-line, initial heading Variable names in the program are, in most cases,
and aircraft velocity, are also input through similar to the corresponding variables of the math-
ACCEPT statements in the initial region. Initial ematical model. The complete program may be
conditions for other variables of the direc- found on the diskette.
tional flight-control system are all set to zero.
The range is scaled to provide values in kilo- program with additional statements to allow
meters instead of meters. This is done mainly for the presence of the outer feedback loop
to make the labeling of the range axis more involving the coupler.
convenient. The structure of the derivative Figures 11.12, 11.13 and 11.14 show the
region is essentially the same as that for the response of the system in terms of plots of the
directional flight-control system simulation lateral deviation, y, versus time for three
10
-5 \~
o 20 40 60 80
t
Is.)
Fig. 11.12 Simulation results from AILS. SLI showing response of the lateral beam guidance system in
terms of the deviation, y, versus time for a coupler gain factor G, of 8.
172 Case study II
different values of the coupler gain (Gc = 8.0, characteristic of a truly unstable system. This
Gc = 16.0 and Gc = 32.0). In each case the behavior is a result of the nonlinear nature of
initial position of the aircraft is 6000 m from the lateral beam guidance system. Examin-
the point of touchdown on the runway and ation of Fig. 11.9 suggests that the system
15 m from the center-line in terms of lateral becomes steadily more sensitive as the air-
error. For the two lower values of coupler craft approaches the transmitter because the
gain the trajectory is clearly of a satisfactory angle A increases for a given lateral error, y,
form and the final lateral error is zero or has as the range, R, becomes smaller. This is
a very small value. With the largest of the reflected in Fig. 11.8 by the presence of a
three gain factors the results show a trajec- block involving a gain factor inversely pro-
tory which is clearly unstable. More detailed portional to R. The loop gain in the outer
examination of this case (Fig. 11.14) shows loop of the feedback system therefore in-
that the initial part of the trajectory involves a creases as the range becomes smaller, and
damped oscillation and only the later part of this explains the type of behavior shown in
the response shows a diverging oscillations Fig. 11.14. Although stable for large values
10
o
-5
o 20 40 60 80
t
Is)
Fig. 11.13 Simulation results for the lateral beam guidance system model for G, = 16.
>-S
30
20
10
-10
0 20 40 60 80
t
Is)
Fig. 11.14 Simulation results for the lateral beam guidance system model for G, = 32.
References 173
of R, the system becomes unstable as R the performance of the model when subjected
decreases. In this case, linear theory would to external heading disturbances. The modifi-
provide relatively little information about the cations to the SLIM program provided are
complex nature of the system and could not relatively minor for both these cases [3].
be used to predict the form of the response
for the situation considered in this investiga-
11.4 DISCUSSION
tion, which involves a large change in R.
Another fact of considerable importance is This case study involves the modeling and
that the performance of both the directional simulation of a relatively complicated non-
flight-control system and the lateral beam linear multiloop feedback control system.
guidance system depend upon the approach Mathematical models of fixed-wing aircraft
velocity Uo• The performance is influenced are well established and their limitations are
considerably by this quantity and also by the well understood and documented. Little con-
initial conditions. Unlike a linear system in sideration is therefore given in this example
which the nature of the system response is the to questions of model validity. Emphasis is
same for different initial states, the lateral given instead to the control systems aspects
beam guidance system model has a behavior of the simulation model and to the use of sim-
which depends to a considerable extent on ulation experiments to reveal the complex
the initial lateral deviation, aircraft velocity and highly nonlinear nature of the system
and initial range. Investigation of the model behavior. Similar insight could not be
behavior for different sets of initial conditions obtained from analytical studies based on
is left to the reader as an exercise using the linearized models. As with the system in
lateral beam guidance model file provided on Chapter 10, a nonlinear computer simulation
the accompanying diskette. model provides a basis for very valuable
In a practical system of this kind extensive testing of the control system prior to any
use is made of gain-scheduling techniques, practical implementation in terms of real
with the coupler gain factor being adjusted hardware. This is particularly important in a
according to range information. This gain- safety-critical application such as aircraft
scheduling feature of the real system over- flight control.
comes the major variations of performance
with range. This has not, of course, been
included in the simple model being consid- REFERENCES
ered here. However, the simulation does 1. Blakelock, J. (1991) Automatic Control of Aircraft
provide opportunities to investigate the very and Missiles, 2nd edn, John Wiley, New York.
interesting nonlinear dynamics of the system 2. McLean, D. (1990) Automatic Flight Control
with constant coupler gain and could easily be Systems, Prentice-Hall, Hemel Hempstead.
3. Murray-Smith, D.J. (1983) Use of an aircraft
extended by the reader to include some form lateral beam guidance system simulation in
of scheduled coupler gain. Similarly, the intro- the teaching of control engineering, in
duction of integral control in the coupler is an Proceedings of the 1st European Simulation
exercise which could be carried out without Congress: ESC83, 1983, Aachen (ed.
difficulty by any reader interested in assessing W. Ameling), Springer, Berlin, pp. 337-42.
CASE STUDY III - RESPIRATORY
GAS-EXCHANGE PROCESS SIMULATION 12
COMMON
DEAD SPACE
ALVEOLAR
-
BY - PASS
TISSUE COMPARTMENT
-- --
Fig. 12.1 A simplified model of the pulmonary gas-exchange processes for carbon dioxide.
V(t)~
(12.9)
0.4 Sin(; )
Provision can also be made within the model
equations for a bypass pathway through and that the time required to fill the dead
which a proportion of the blood leaving the space (7") is 0.6 s. The three stages of the
tissue compartment avoids both of the alveo- breath cycle are therefore as shown in
lar compartments. Such a bypass is mainly of Fig. 12.2 for this particular case.
interest for the simulation of abnormalities
and has not been included in the set of equa-
tions presented here.
12.3 THE COMPUTER SIMULAnON
The ventilation, which appears in the equa-
tions describing the alveolar compartments Let us now start considering the task of
and in the relationships defining the three implementing this simulation model using an
stages of the breath cycle, is assumed to be equation-oriented continuous system simu-
sinusoidal in form and is given by the lation package such as SLIM or ACSL. The
equation model equations are already in a convenient
form for simulation and involve a set of three
V(t) =Vrnax sin(2rrft) first-order ordinary differential equations and
additional algebraic relationships. However,
where V max is the amplitude of the sinusoid there is a complication in that the equations
and f is the frequency of breathing. The ven- involve products of variables (e.g. the
tilation level may be defined in terms product of the ventilation and [PD(t) - PA1 (t}]}.
of the tidal volume, Vp This is the volume In addition, there are some rapidly varying
of gas inspired and expired during each terms in the equations owing to the presence
breath cycle, and it is a simple process to of the quantities 51 and 52' which switch
relate this to the amplitude of the sinusoid between values of unity and zero according
above, since to the stages of the breathing cycle. The
f
2 IU
~
VT =tidal volume = V(t)dt (12.10) .... -
::;)-:=
o
o > 0.4 --
f
T
Fig. 12.2 The three stages of the breath cycle for
V(t)dt =VD (12.11) the case of a 4-s breath cycle period, tidal volume
of 0.5 I and a dead space volume of 0.1 1.
180 Case study III
generation of the switching quantities ward, apart from the additional lines of code
requires an extra integration function in addi- needed to generate the switching quantities.
tion to those needed to represent the actual One unusual feature, which could be of value
compartments. Stiffness is not a problem in in demonstrating the effect of changes in
this model, and no difficulties arise over the quantities such as ventilation or blood flow
choice of an integration method. A number rate, is the inclusion of facilities for applying
off standard integration techniques have been step changes of parameters during a simula-
found to give satisfactory results. tion run, and this requires additional state-
The INITIAL section of the simulation ments within the DERIVATIVE segment.
program can be used to establish parameter In simulation languages such as ACSL a
values and other quantities such as the initial simulation user can readily alter constants
conditions for the partial pressures in the which appear in the program using runtime
lung and tissue compartments. The commu- commands without having to edit the simu-
nication interval must be chosen to provide lation program or change data files. As
an adequate number of samples within the explained in Chapter 6, this is a very useful
minimum expected breath cycle period. In facility and is ideal for a teaching environ-
this connection it must be remembered that ment since it provides a means of separating
the objective of this simulation model is to the model itself from the experiments to be
demonstrate gas-exchange phenomena for a carried out on it. Students can use the simu-
wide range of normal and abnormal condi- lation very effectively by employing such
tions. Parameters of the model will be set to a commands to carry out a variety of tests on
variety of different values, and Table 12.1 the model for as many parameter sets and
gives appropriate ranges for the parameters initial conditions as is necessary.
most likely to be adjusted by the user. Figure 12.3 shows relevant sections of a
The structure of the DYNAMIC and DERIV- SLIM listing of a program implementing the
ATIVE sections of the program is straightfor- gas-exchange model described above. The
Fig. 12.3 DYNAMIC segment of a SLIM program for simulation of the gas-exchange model. Variable
names in the program have been chosen to match the corresponding variables of the mathematical model
as closely as possible.
source of the SLIM program may be found on overall results credible? What is the range of
the diskette (LUNGS. SLI). The set of parame- conditions, if any, for which the simulation
ters used (LUNGS. IN) may be regarded as can be used to make useful predictions?
typical of a 'normal' adult human subject. External validation studies have been
Some selected results obtained using this carried out involving comparisons of simula-
SLIM program are shown in Figs 12.4, 12.5, tion results for a particular case of the model
12.6 and 12.7. Similar results may, of course, described above, with experimental records
be obtained from XANALOG and SIMULINK from a human subject. The essential differ-
simulation models, and Fig. 12.8 shows the ence between the simulation used in this case
XANALOG block diagram which is equiva- and the model equations given above is that,
lent to the equation-oriented program listings. in order to allow direct comparisons to be
Questions should always be asked about made with the experimentally measured
the validity or otherwise of any simulation record, the simulation cannot be driven using
model. Clearly a model structure of the kind a sinusoidal ventilation pattern. Instead, ven-
shown in Fig. 12.1 is a highly simplified rep- tilation data logged by computer during an
resentation of a very complicated biological experiment on a human subject must be used
system. Are the simplifying assumptions as an input to the simulation. The model has
used in deriving the model equations justi- been successfully validated in this way for
fied? Does the simulation program accurately normal human subjects and has provided a
represent the model equations? Are the basis for the successful development of a
182 Case study III
35
~
o 50 100 150 200
t
(5)
Fig. 12.4 Response of the simulation model to a step change of carbon dioxide level in the inspired gas
mixture from 0 to 20 mmHg. The continuous curve shows the partial pressure of carbon dioxide in each
of the alveolar compartments, while the dashed line shows the corresponding partial pressure within the
tissue compartment. The carbon dioxide input is applied after 30 s. Note the cyclic nature of the alveolar
partial pressure record, which is caused by the periodic breathing pattern.
-l
a,.. E 38
- E
34
o 50 100 150
t
(5)
Fig. 12.5 The alveolar compartment partial pressure in response to simulated breath holding for a period
of 45 s (from t = 30 s to t = 75 s).
noninvasive method for estimating cardio- piratory 'plant' in simulation studies of the
pulmonary quantities such as the effective control systems responsible for the regulation
lung volume and the total blood flow and control of breathing [11].
through the lungs [8, 9]. The model has been
found to be of more limited value for abnor-
12.4 DISCUSSION
mal subjects, but can nevertheless provide
valuable qualitative insight [10] and has pro- This case study is of particular interest as an
vided a basis for the representation of the res- example of the use of simulation techniques in
Discussion 183
QDOT1
PA2 PT
QDOT2
LUNGPA2_MDL
Fig. 12.8 XANALOG simulation model diagram for the gas-exchange model. Note the use of submodels
for the three compartments, for the ventilation (DV) and for the calculations associated with the switch-
ing quantities (51 and 52).
may reduce significantly the number of experi- in terms of amplitude and frequency, on the
ments required to obtain meaningful results. carbon dioxide concentration within the
Application areas in which simulation-based various compartments. It can also be used to
investigations involving lung gas-exchange investigate the sensitivity of the model vari-
models are of potential value include the ables to the parameter k which determines
development of automatic control systems for how the ventilation is distributed between the
the administration of anesthetic gases, the two alveolar compartments (e.g. Fig. 12.6).
development of improved tests of lung func- Similarly, the blood flow rate parameters Ql
tion and the noninvasive monitoring of car- and Q2 can be varied to simulate abnormalities
diopulmonary quantities. The model described in the pattern of lung perfusion (e.g. Fig. 12.7).
in this study is compartmental in nature and The effects of other parameters such as the
only some of the compartments involve quan- tissue compartment volume V TC and the
tities which can be measured directly in the volumes of the alveolar compartments, VAl
human subject. Issues of external validation and V AU are also of interest.
are therefore of special interest.
The SLIM program LUNGS. SLI can be used
to investigate the behavior of this gas- REFERENCES
exchange model for a wide range of different 1. Vander, P.J., Sherman, J.H. and Luciano,
conditions. It can be used, for example, to D. (1970) Human Physiology, McGraw-Hill,
establish the effect of the ventilation pattern, New York.
References 185
2. West, J.B. (1965) VentilationfBloodFlow and 7. Murphy, T.W. (1969) Modelling of lung gas
Gas Exchange, Oxford University Press, exchange - mathematical models of the lung:
Oxford. the Bohr model, static and dynamic ap-
3. Emery, B. and Pack, A.1. (1971) An experi- proaches. Mathematical Biosciences, 5,427.
mentally verified model of the gas exchange 8. Ferguson, D.R, Mills, RJ., Moran, F., et al.
properties of the lung, in Computers for (1973) Estimation of the parameters of a
Analysis and Control in Medical and Biological lung model with clinical applications, in
Research, IEE Conference Publication No. 79, Identification and System Parameter Estimation
Institution of Electrical Engineers, London, (ed. P. Eykhoff), North-Holland, Amsterdam,
pp.130-34. pp.213-19.
4. Pack, A.I., Emery, B., Moran, F. and Murray- 9. Bache, RA., Gray, W.M. and Murray-Smith,
Smith, D.J. (1974) Computer models of gas OJ (1981) Time-domain system identification
exchange processes in ventilation, in applied to non-invasive estimation of cardio-
Ventilatory and Phonatory Control Systems (ed. pulmonary quantities. lEE Proceedings, Part 0,
B. Wyke), Oxford University Press, Oxford, 128,56-64.
pp.227-47. 10. Bache, RA. and Murray-Smith, D.J. (1983)
5. Murray-Smith, D.J. and Pack, A.1. (1977) Structural and parameter identification of two
Techniques of computer simulation applied to lung gas-exchange models, in Modelling and
respiratory gas exchange, in Non-invasive Data Analysis in Biotechnology and Medical
Clinical Measurement (eds D.E.M. Taylor and Engineering (eds G.c. Vansteenkiste and P.c.
J.5. Whamond), Pitman Medical, Tunbridge Young), North-Holland, Amsterdam, pp.
Wells, pp. 186-202. 175-83.
6. Mills, RJ., Middleton, S., Moran, F., et al. 11. Greer, W., Jordan, M.M. and Murray-Smith,
(1974) Simulation in the teaching of concepts OJ (1982) A structural approach to the respir-
of respiratory gas exchange. International atory control simulation, in Computers in
Journal of Mathematics in Education, Science and Medicine (ed J.P. Paul), Macmillan, London,
Technology, 5, 381-87. pp.339-46.
CASE STUDY IV - A SIMULATION
MODEL OF ACTIVE SKELETAL MUSCLE
13
r I
I
J
Surface
of cell Transmitter I I
I I
I
I
Surface
of cell
Muscle I I
I Mechanical I Mechanical
cell properties properties
-
-
-
-
- -
- - - - -
-1- I
- membrane
I
I
of muscle
I
of load
- Surface
of cell
I
I
I
I
I I
2 3 4
Fig. 13.1 A schematic diagram illustrating the control sequence which links electrical events in the
motoneuron to the mechanical response observed in skeletal muscle.
188 Case study IV
mechanical responses will occur. If the stimu- directly to length. Viscous elements are simi-
lus is repeated in a periodic fashion, at a suffi- larly defined as those in which the tension is
ciently high frequency, the twitches will related directly to velocity. Hill has summar-
merge to produce a contraction which fluctu- ized [2] the results of experiments made to
ates about a mean level which depends upon determine the relationships between length
the applied frequency. The condition in and tension changes. The experiments show
which the decay between individual twitch the following:
responses disappears and a steady contrac-
1. active muscle contains an undamped
tion is maintained is known as tetanus. When
elastic element;
tetanic conditions have been achieved any
2. active muscle contains an apparently
further increase in the frequency of the stimu-
damped element in series with the
lus can produce no further change in output.
undamped elastic one; and
It is possible to investigate the sequence of
3. resting muscle contains the elastic
events which follows the arrival of a pulse
element, but only to a minor degree the
and leads to mechanical activity in the
apparently damped element.
muscle. The complicated system must be
broken down into smaller units which can be The undamped element is termed the series
identified structurally and then isolated for elastic element and the damped element is
independent study. Figure 13.1 shows the the contractile component. Investigations of
control sequence which is commonly sup- the properties of muscle over a wide range of
posed to exist between the input motor nerve lengths indicate that models of resting or
impulses and the mechanical responses of the active muscle must contain a second elastic
muscle. The arrival of a nerve impulse causes element in parallel with the first two. This
the release of a chemical transmitter. The third element is the parallel elastic element.
muscle cell membrane is acted on by the These three elements have been incorporated
chemical transmitter to produce changes in
the mechanical properties of the muscle. The
conceptual boundary between the block rep-
resenting the muscle cell membrane and the
block representing the mechanical properties
of muscle remains somewhat indistinct.
Many models describing the detailed mech-
anisms of muscle contraction have been pro-
posed and, although of great value to Parallel Series elastic
elastic element
physiologists, they do not describe the inter- element
action of a given muscle within a multicom- Contractile
element
ponent system such as the neuromuscular
control system. In this case, the model may be
more empirical in form and can be based on
data obtained from simple measurements.
The concept of muscle as a system of ele-
ments whose mechanical properties could be
described in terms of elasticity and viscosity
was first given clear formulation by Hill in
1922 [1]. In the lumped parameter analogies Fig. 13.2 A diagram showing the three elements of
developed by Hill elastic elements are a simple mechanical analog of active skeletal
defined as those in which tension is related muscle.
Models of active muscle 189
in a conceptual model widely used for many element, and the parallel elastic element can
years to describe the properties of skeletal therefore usually be neglected under these
muscle. This model is shown diagram- conditions. The resulting two-element models
matically in Fig. 13.2. and the associated characteristic curves have
In order to characterize a muscle in terms been the subject of many investigations.
of the three-element model structure dis-
cussed above, at least five curves must be
13.2 MODELS OF ACTIVE MUSCLE
obtained experimentally [3]. These comprise
the force-velocity and force-length curve of As outlined above, the overall mechanical
the complete active muscle; the force- properties of active skeletal muscle can be
extension curve of the series elastic element; described by a set of five characteristic curves.
the force-length curve of the parallel elastic Of these the force-velocity curve is generally
element; and a curve describing the temporal believed to have special significance and has
properties of the muscle activity following the been studied in considerable detail by a
arrival of a motor nerve impulse (the 'active number of different investigators. Each point
state' curve). For values of muscle length on the curve is normally obtained from a
much less than the normal rest length, the measurement of the initial gradient of the
stiffness of the parallel elastic element is very length-time curve of the muscle when con-
much less than the stiffness of the contractile tracting actively in response to a maximal
0_06 1\
,..
\
"I
\J
m
a 0.04 \ ~
~
>t 1\
.ro!
U
0
\
roI 0 . 0 2
>
~ "'" ~
~ I--
I----
0.00
0.0
Force eN)
Fig. 13.3 The force-velocity characteristics of an active frog muscle (based on data of Ritchie and Wilkie
[3]). The circles show experimental points and the continuous curve shows the fit obtained using a
rectangular hyperbola of the form of equation (13.1) with parameters a = 0.143 and b = 0.0175.
190 Case study IV
nervous stimulus applied under conditions of and measuring the tension developed for a
constant tension. Hill's equation [2] give's a number of different values of the overall
good fit to the force-velocity curve for most length. The curve can be approximated quite
skeletal muscles during active contraction. closely by a parabola having the general form
This equation, which has the form of a rectan-
gular hyperbola, may be written in the (13.2)
following way:
In this equation p is the muscle tension, 1 is
(p + a)(v + b) =(po + a)b (13.1) the overall length of the muscle, 10 is the
normal body length of the muscle and Pc is
In this equation P is the force developed by the tension at body length 10 • Figure 13.4
the muscle, v is the speed of shortening, Po is shows the tension-length equation fitted to
the muscle tension at zero velocity for each experimental points found by Ritchie and
value of length and time and a and b are con- Wilkie [3]. The values of k, 10 and Pc are 2060,
stants. Figure 13.3 shows a set of experimen- 0.032 and 0.481 respectively. It should be
tal points for frog muscle obtained by Ritchie noted that, although the theoretical curve
and Wilkie [3] and the corresponding force- defined by equation (13.2) is symmetrical
velocity curve with a = 0.143 and b = 0.0175. about length 10, only the part for 1 less than 10
The tension-length curve of the complete is shown owing to the absence of published
active skeletal muscle can be determined experimental data for 1greater than 10-
experimentally by imposing tetanus on the The length of the series elastic element
whole muscle under isometric conditions cannot be measured directly, but by means of
,...
""
z
0.5
0.4
V
/ --
~
V
0
~ 0.3 I
V
0
~
~ 0.2
/
/
~
0
~
~ 0.1
/
/
~
0.0
J
0.016 0.024 0.032
Fig. 13.4 The tension-length characteristics of an active frog muscle (based on data of Ritchie and Wilkie
[3]). The circles show experimental points and the continuous curve is the parabolic relationship of equa-
tion (13.2) with parameters k = 2060, 10 = 0.032 and Pc = 0.481.
Models of active muscle 191
a series of indirect measurements a curve can the rising phase of the active state curve is
be obtained to relate the extension of this not readily determined with precision, but it
element to the applied tension. This tension- is clear that the contractile element becomes
extension curve gives enough information to an active source of tension a very short time
permit the calculation of corrections to the after the arrival of a nerve impulse at the
tension-length curve of the contractile ele- end-plate. The tension developed in the con-
ment itself. The tension-extension curve of tractile element remains constant for a
the series elastic element (Fig. 13.5) may not period after the application of the stimulus
be expressed in a simple mathematical form and then falls slowly towards zero. Figure
as readily as the other curves, and a table 13.6 shows the active state curve and the cor-
look-up form of function generation is re- responding experimental points found by
quired for the representation of this element Ritchie and Wilkie [3]. It has been assumed
of the model. that the int~al rise in the active state is
The active state curve shows how the instantaneous.
degree of activity of the contractile element The experimental data of Ritchie and Wilkie
changes with time following the application [3] includes results of a series of isoton\: (con-
of a nervous stimulus. The precise form of stant tension) after-loaded twitch experilnents.
0.30
,.. 1
0.20
..,~
D
0
J.f
I
IJ
0
fit
0.10
0.0
0.000 0.001
-- V--
Extensic:>n
V
0.002
<:Ill>
0.003
Fig. 13.5 The force-extension curve describing the series elastic element for the frog muscle used in the
experiments by Ritchie and Wilkie [3]. The circles show experimental points, while the continuous curve
is the corresponding function to be fitted for the simulation model.
192 Case study IV
1.0
~
W
+l
!Ij
+l
m 0.8 ~
W
>
.,.j
0.6
\ \
+l
\
0
<
0.4
0.2
\ 1\
0.0
~ "-...
0.0 0.2 0.4 0.6
'I'i:me (8)
Fig. 13.6 The active-state curve for the frog muscle, with experimental points found by Ritchie and
Wilkie [3] indicated by circles.
To produce an isotonic after-loaded twitch, effect of the parallel elastic element may be
the muscle is held initially under isometric neglected for values of I less than 10, the
(constant length) conditions by a mechanical tension in the contractile element must be
stop which prevents excessive extension of the equal to the tension in the series elastic
inactive muscle arising from the applied force. element. Equation (13.2) may be used to
When the muscle is subjected to a nervous determine the tension at zero velocity for the
stimulus tension is developed, and when the chosen length, assuming the muscle to be in a
active tension developed by the muscle is state of maximum activity. If, on the active
equal to the applied tension contraction starts. state curve, the intensity of the activity at any
The after-loaded twitch therefore involves instant of time is expressed as a fraction of the
both isotonic and isometric conditions at dif- maximum activity, the tension Po at zero
ferent stages in the response. velocity for a given length I must be described
by the equation
13.2.1 SIMULATION OF ACTIVE MUSCLE
(13.3)
FOR CONSTANT LENGTH CONDITIONS
Under isometric conditions the overall length From equation (13.1) the velocity of contrac-
of the muscle is rigidly constrained. Since the tion of the active muscle is given by
Models of active muscle 193
b(po - p) the contractile and series elastic elements
v = ----'----=--- (13.4) must at all times be equal to the applied
p+a
tension. An increase of activity causes the
Under isometric conditions, any change of contractile element to shorten and, since the
length of the contractile element 111c must tension in the series elastic element is con-
cause an equal and opposite change of length stant, the change of length of the complete
111, in the series elastic element. The nonlinear muscle must be equal to the change of length
force-extension curve of the series elastic of the contractile element. Figure 13.10 is a
element gives the tension change appropriate block diagram of the model for isotonic
to the length change 111, and this must be the conditions.
overall tension p developed by the muscle. Figure 13.11 shows the block diagram for
Figure 13.7 is a block diagram of this simple the simulation of an after-loaded isotonic
model. Figure 13.8 shows an XANALOG sim- twitch. Here the muscle is allowed to develop
ulation diagram for this model, and the tension isometrically until the tension reaches
results of a simulation of the isometric twitch the required isotonic value, after which time
response obtained using the XANALOG soft- isotonic conditions are applied. It should be
ware are shown in Fig. 13.9. noted that some intermediate variables have
been introduced on this diagram. These vari-
ables are used in subsequent sensitivity
13.2.2 MODELS OF ACTIVE MUSCLE FOR
analysis studies. Therefore, the system being
CONSTANT TENSION CONDITIONS AND
simulated is initially described by the block
SIMULATION OF THE AFfER-LOADED
diagram of Fig. 13.7 with a change to the
ISOTONIC TWITCH
structure of Fig. 13.10 when the tension
Under isotonic conditions the overall tension reaches the required level.
of the muscle is controlled by external Figure 13.12 shows the XANALOG simula-
means, and the tension developed in both tion diagram for the after-loaded isotonic case,
ActiV"e
[Link]
h
Series a
Musc1e [Link]
e1entent
1ength
1
Fig. 13.7 A block diagram for the nonlinear model for active skeletal muscle for conditions under which
the muscle length is held constant (isometric conditions). Note the use of multiplier, divider and non-
linear function blocks in addition to the standard elements representing integration, summation and
multiplication by a constant. In this diagram the active state (a function of time) is represented by an
external function generator.
b
ACTSTATE_T
TABLE
'----I K
.1..43e-01
4.8119-01 P_DELTAL
TABLE
SERIES
~-ELASTIC"rJ
ELEMENT
[D...... [Link]
3.20e-0 2
Fig. 13.8 XANALOG simulation diagram for the muscle model under isometric conditions. This is equiv-
alent to the block diagram of Fig. 13.7. Note the presence of the INVERSE block which has been created by
the user to provide a form of divide operation when used in conjunction with a multiplier block. Note
also the use of TABLE look-up functions for the active-state function and series elastic element. The
active-state block has time as its input, as shown by the 'clock' icon. The principal output variable is the
muscle tension, p, as shown by the 'eye' icon.
·... ·_· ..._·_·l-·... ·... ·...·_-.-_·__·__·j· ...·...·...·... ·__·j·... ·...·... ·...·-· ...·i· ... ·...·... ·... ·...·...·_,_· ... ·..·... _· ... ·_j_·...· ... ·... ·... ·... ·_j
i i i i i . i .
i i i i i i i I
0.5 .-!i~j·
I , I I I I I I
! ! -! ~ ~ i ! !
1 ! ! I
.-+~
!
I
....-.-.-.-+.-.-.---.-+--.-.-.-.-.+.-.--....
I f I
-.~
!
I
..... !I .. _......... _!.I~
~-
i i i i . i i
. i i i I ; i
0.4 ....-~. .....-.+.-.--.--.-...;.-.-.-.....-.-.•-............._.....--
i i i i i
. i j i
I i i I i i
_. __ ._-+-_............._. .-..... -+.;~ ..;.-.-.-.-........ _j.
i i ji i i j
i i j. i i j
j i i i i ;
i i i i i i I I
-_.- .-.j--.-.-.--;--.-.-.-.--r-. ·-T-·-·-·-·-·-·r-·-·-·-·--j-·-·-·--·-·j--·-·-·-·-·-i
i i i i i i j i
! ~
._..J...•.__.__.-L_. __.__.L.. __.__
! !I ._.•.
!
_._._.!.•.___•.
! _~.J
. ..._._._._._!!
0.2 I I I I I I I
i i j i i j i
i i j I i i i j
.-i~ _.~-i
I I I I I I I I
i
I
!
I
!
I
!I !
I
!
I
!
I
'
I
0.1 ---t--·-----!----·-+-·--···-4-
I I I I
.-~+_!
I I I I
! ! ! ! ! ! i !
--.L.-----i---l---._-!..._..___... -___.L.__.__.--!______.J
!
I
!
I
!I iI iI . i !
I
! ! ! i! !
0.0
0.0 0.2 0.4 0.6 0.8
[Link] (8)
Fig. 13.9 Response of the muscle model to a single nerve impulse under isometric conditions, obtained
using the XANALOG simulation model of Fig. 13.8.
Models of active muscle 195
[Link]
state
h
Applied.
tension p
-~+{
[Link].l
length Ii
Fig. 13.10 A block diagram for the muscle model under conditions of constant tension (isotonic conditions).
Active
[Link]
h
Applied. 0
tension P PA
P. .-J "'-1r----....:.P- - - - + -.. .
A,-__....L.._ _- - ,
Series
e l ...... t i c S e r i e ...
ele_ent [Link]
e.1..[Link]
Initia.l
length Ii
Fig. 13.11 A block diagram for the muscle model under 'after-loaded' isotonic conditions.
and the results of XANALOG simulation conditions also increases. Experimental results
experiments based on this are shown in by Ritchie and Wilkie [3] are also shown, and
Fig. 13.13. These show length changes versus it may be seen that the agreement between the
time for various different values of applied experimental responses and the simulation
tension. It can be seen that as the tension is results is least satisfactory near the peak of the
increased the time taken to establish isotonic twitch for the smaller values of load.
196 Case study W
b
ACTSTATE_T
TABLE
.1. .. 4 : 3 8 - 0 1
[Link]
TABLE
SERIES SERIES
,....--_......:ELASTIC ..._ _ _ _~-=[Link] . .---------1
ELEMEH LEII'JENT
[Link] 6L-t ~p
po
.
.:
: 2 .. 06e+03 4 .. 818-01.
t..................•
LO L-LI
[!]. [Link]
3 .. 2 0 e - 0 2
Fig. 13.12 XANALOG simulation diagram for the after·loaded isotonic case of Fig. 13.11. Note the use of
the switching blocks to provide the change of structure in the model at the transition from isometric to
isotonic conditions.
014~3 IN
L
0.008
'(
0 V
V
..
0/ ~+0436 N
,...
S
\01
~
0.006
VI V
~
C
m
j +/
1+ V
.c:
t)
.c: 0.004
:x.-
v
I
~ x
~
C
m .I x k1515
x/,
N
l fl xL V
~
0.002
I~ ;/v • •
.....-••
j., I~ ~ N
V
0.2j80
0.000 J
0.00 0.10 0.20 0.30 0.40
Ti:rne (8)
Fig. 13.13 Simulated responses (continuous lines) and experimental results (circles and crosses) for the
first phase of the after-loaded isotonic twitch response, during which the muscle is contracting. Four sets
of results are given for different values of applied tension, ranging from 0.0143 N to 0.298 N. The experi-
mental results are by Ritchie and Wilkie [3].
(13.12) ar ap (13.15)
---+YI =0
aq aq
(13.13)
(13.14)
az ad
-+-+Yz =0 (13.17)
aq aq
198 Case study IV
Table 13.1 Functions appearing in the sensitivity
Jpo _ h Jz = 0 (13.18) co-system equations
Jq Jq
Parameter 1'1 1'2 1'3 1'4
(13.19) a -1 0 0 0
b 0 0 -m/r 0
k 0 0 0 -['
(13.20) p, 0 -1 0 0
10 0 0 0 2kf
(13.21)
tension-length diagram of the series elastic
element [4].
Figure 13.14 shows the block diagram of
(13.22)
the general co-system. Figure 13.15 shows the
relative sensitivity coefficients for a number
The functions 1'1' 1'2' 1'3 and 1'4 which appear in of different parameters for a single value of
the co-system equations above are given for applied tension. Figures 13.16 to 13.20 show
each model parameter in Table 13.1. The the sensitivity coefficients for each of the
functions 1'5 and 1'6 are gradient functions parameters of Table 13.1 for a series of values
of a piecewise-linear approximation to the of applied tension.
Ac:ti"e [Link]
h
o
r Ie---Ys
Fig. 13.14 Block diagram of the sensitivity co-system for the muscle model under after-loaded isotonic
conditions.
Parameter sensitivity analysis 199
0.010
P 01
cOp? ~ ---
v::: V
~ - b01
-
Ob
l
/
V
~
~
~
~
'\
'\
- I---
-- t---
ko1-
ok
~
...........
\:--
\ a01
-
oa
\ al0
10 01
-0.010
0.0 0.1 0.2 0.3
T.i:me (s)
Fig. 13.15 Relative sensitivity functions for four parameters of the muscle model for an after-loaded iso-
tonic simulation with an applied tension of 0.2980 N.
The sensitivity coefficients provide an accu- give a better approximation to the experimen-
rate means of estimating the effect of small tal response curves in the vicinity of the peak
changes of parameters on any of the variables of the twitch for the two smaller values of
of the model. The sensitivity functions of Fig. applied tension. However, the same change
13.15 show very clearly the relative impor- does not improve the model response for the
tance of the parameters la, Pc, a and b in the two larger values of tension. Therefore, it is
simulation, and show that near the peak of likely that the errors in the model response
the twitch, where the error in the model are because of a number of factors which
response is greatest, the response is par- include structural features as well as para-
ticularly influenced by the parameter 10 , metric uncertainties. It is of interest to note
Examination of the response curves of Fig. that, in almost all cases, the parametric sensi-
13.13 and their sensitivity functions shows tivity falls as the applied isotonic tension is
that an increase of 10 of the order of 0.3% can increased.
0.06
~I
+l
-.l
),
0.04
v -I--- ~Ot43
I
IN
---;
V-
-.l
+l
-.l
>
V '/ 0.0436 N
~
~
rn
/ V
--
0.02
Ul [Link]
/V
N
-- -
~
l.---
l~ 0.2980 N
0.00
- -1
0.0 0.1. 0.2 0.3
Ti:rne (a)
Fig. 13.16 Parameter sensitivity functions for parameter a in the muscle model under after·loaded iso-
tonic conditions for a number of levels of applied tension.
~Ii@ 0.5
),
+l
.,.j
0.4
>
.,.j
+l N
.,.j
rn 0.3
~ [Link] N
~
Ul
0.2
Time (a)
Fig. 13.17 Parameter sensitivity functions for parameter b in the muscle model for the after-loaded iso-
tonic case.
~I
.1.2
/
V
./
--
0.0.143 N
>t
+l
/
.~
>
V
.~
+l O.B /
Y
0.0436 N
oM
00 V---
C ~
m
Ul
0.4 / /
/
V
/ ,/
0.0
---V- V ./
- ---
0 • .15.15
0~29BO
N
Ti:me (5)
Fig. 13.18 Parameter sensitivity functions for parameter k in the muscle model for the after-loaded iso-
tonic case.
~I
>t 0.0436 N
~ O . 02 r---+---j---+-::=--===+--=;;;;;j
0.0.143
> N
.~
0
I I
• .15.15 N
+l
.~
00
C
m
Ul 0_29BO N
o . 01 t----flf+----iC-.----+--::",........::::::=f-----/
0_00
0.0 0.1 0.2 0.3 0.4 0.5
Ti:me (5)
Fig. 13.19 Parameter sensitivity functions for parameter p, in the muscle model for the after-loaded iso-
tonic case.
t)~( ~\ 0
0.6
I \
0.0.143 N
r-
), ~
.,.,+J
I"
/'
.>
."+J
."In
0.4 / V
0.0436
----
~
N
~
rn
C
/ /
/
V
0.2
V
/ V
/
---
0 . .15.15 N
/ / ~
--- ---
~
0.0 ~ / ~ 0.2980
1
N
Time (a)
Fig. 13.20 Parameter sensitivity functions for parameter 10 in the muscle model for the after-loaded iso-
tonic case.
0.0.10
,...
S 0.008
rI ~
"'"
~
~ P .. = 0.0436 N
'"
tl'
C 0.006
~
~
.d
t)
.d
+J
~
tl'
C
~
0.004
0.002
I( \
I{
......
'\
= 0 • .15.15
'"~
~
0.000
0.0
\ .1.0 2.0
~
" 3.0 4.0
Time (a)
Fig. 13.21 Simulation results for the muscle model: twitch responses under after-loaded isotonic condi-
tions for two levels of applied tension.
Simulation of the complete twitch response and repetitive stimulation 203
13.4 SIMULATION OF THE COMPLETE of the simulation model may be understood
TWITCH RESPONSE AND REPETITIVE from an examination of the block diagram of
STIMULATION Fig. 13.11, which shows that when the active
state is zero the variable Po must be zero.
Although results in Ritchie and Wilkie's Hence, the velocity must take on a constant
experiments [3] were available only for the value which depends upon the applied
rising phase of the isotonic twitch, it is of tension. These results show general qualita-
some interest to examine the response of the tive agreement with experimental records of
model over the complete period following Jewell and Wilkie [5] which show an approx-
excitation, including both contraction and imately linear phase in the relaxing response.
relaxation. Figure 13.21 shows simulation It should be noted that there is evidence
results for this case and indicates that there is which suggests that Hill's equation for the
a linear relationship between length change force-velocity characteristics of the contractile
and time during the relaxing phase after the element is not strictly valid during the relax-
active state has decayed to zero. This behavior ing phase. Hill's equation has been used in
0.016
I I
f > 11 p/s
.-
0.012
V /"
~
r/
'" r ~
,..
~
I I /f = 1 p/s
'oJ
Gl
~ 0.008
V0..... I
~
~
! "'"
""""""
.c
u
.c
~
01
0.004
I
C ~
'"
~
1-1
~
0.000
0.0 1.0
Ti:me
2.0
(s)
"'" 3.0
Fig. 13.22 Simulation results for the muscle model: the response to repetitive neural stimulation under
isotonic conditions. The response to a single nerve impulse is shown, together with responses to stimula-
tion at 1 pulse per second and at over 11 pulses per second. The applied tension in this simulation experi-
ment was 0.0436 N.
204 Case study IV
I
f :> 11 p /5
0.48
/
~
(
.........
,...
Z t~ i ~ ~
1\ \
\
'oJ
k
..,0
\
~ f = :3
\
00
k
W 0.24
8
r-l
W
1\
~
U
00
~ \
0.00
) 0--
0.00 0.25 0.50 1.00
Ti:rne (5)
Fig. 13.23 Simulation results for the muscle model: the response to repetitive neural stimulation under
isometric conditions. The response to a single nerve impulse is shown, together with responses to stimu-
lation at 3 pulses per second and at over 11 pulses per second. The value of muscle length used for this
simulation experiment was 0.032 m.
this simulation model because of lack of must be regarded as acting as a single motor
experimental data for negative values of unit. A typical muscle has hundreds of motor
velocity. The simulation program can be units which, under normal conditions, are
modified easily to accommodate other stimulated in an asynchronous fashion. A
force-velocity curves. closer approximation to a representation of
Figure 13.22 shows the results of simulat- muscle operating under normal conditions
ing the application of repetitive neural stimu- could perhaps be obtained by a parallel
lation at various frequencies under isotonic arrangement of a number of models, but with
conditions. Figure 13.23 shows similar a random or prescribed statistical pattern of
results for the isometric case. The results stimulation.
show very clearly the summation of
twitches, which is a well-known feature of
13.5 DISCUSSION
muscle responses and the limiting condition,
known as tetanus, in which the individual This case study is of particular interest in
twitch responses fuse together to produce a terms of modeling and simulation methods
smooth contraction. for a number of reasons. Firstly, it provides
It should be noted that because of the an example of a nonlinear model which has
method of stimulation used in the experi- been developed on an empirical basis from
ments by Ritchie and Wilkie [3], the muscle measured response data. Secondly, it pro-
References 205
vides an illustration of function generation their most economical speed. Journal of
through curve fitting (the overall tension- Physiology, 56, 19-41.
length curve and the force-velocity curve) 2. Hill, A.V. (1938) The heat of shortening and
dynamic constants of muscle. Proceedings of
and through use of a table look-up approach the Royal Society, B126, 136-95.
(the series elastic element). It also provides a 3. Ritchie, J.M. and Wilkie, D.R. (1958) Dynamics
useful illustration of the generation and use of muscular contraction. Journal of Physiology,
of parameter sensitivity functions in model 143,104-13.
development. 4. Murray-Smith, D.J. (1970) A Model of Active
Muscle, Control Systems Group Report No.
CSG/70/002, Department of Electronics and
Electrical Engineering, University of Glasgow.
REFERENCES
5. Jewell, B.R. and Wilkie, D.R. (1960) The
1. Hill, A.V. (1922) The maximum work and mechanical properties of relaxing muscle.
mechanical efficiency of human muscles and Journal of Physiology, 152, 30-47.
REAL-TIME SIMULATION 14
il
,.\ ('1"
o~-
r~
Fig. 14.2 Circuit diagram of an inverting opera-
Hence tional amplifier acting as a summing unit.
Thus
and thus
and so
if
A coefficient unit provides a means of multi-
i, plying an analog variable by a constant value.
Such units provide a means of setting coeffi-
v, cients which arise in the model equations. In
vo terms of hardware, the simplest form of
coefficient unit is a simple electrical poten-
0
1-1 0
tiometer, as shown diagrammatically and
symbolically in Fig. 14.6. In general, the
Fig. 14.4 Circuit diagram of an inverting opera- coefficient setting, k, is restricted to a range of
tional amplifier acting as an integrator unit. zero to one.
Fig. 14.6 Circuit diagram and equivalent symbol for a coefficient unit.
Analog simulation techniques 213
14.3.5 MULTIPLIERS AND OTHER The interpretation of results from a scaled
NONLINEAR UNITS simulation is very straightforward. For
example, for a problem variable, y, represent-
Although summing units, integrators and
ing a displacement with a range of -50 m to
coefficient units meet all the needs for simula-
+50 m the scaled value of y would be
tions involving linear models, there is often a
need for other operations in the case of non-
ysc = Y/ 50. If at some instant during a simula-
tion run the value of the variable Ysc is 0.2 MU,
linear models. Electronic analog multipliers
the corresponding value of the physical vari-
provide a means of multiplying two variables
able y must be
together, and there are ways of using multi-
pliers in conjunction with high-gain ampli-
y=50ysc =0.2x50=10m
fiers to produce 'units which can divide one
variable by another or carry out square-root
The processes involved in amplitude
operations. Details of such elements, as well
scaling are best illustrated by an example.
as others which are available for more
Consider the unforced mass-spring-damper
specialized tasks, may be found in texts
system shown in Fig. 14.7, where a typical set
dealing specifically with analog simulation
of parameters are M = 20 kg, R = 2000 N m- I s
methods (e.g. refs 3 and 4).
and K = 12000 N m-I . Appropriate maxima for
the displacement y and velocity dy / dt are
14.3.6 SCALING OF VARIABLES 0.05 m and 0.5 m S-I respectively. The second-
order ordinary differential equation defining
Variables within an analog simulation are
a linearized mathematical model for this
restricted in their range by the physicallimi-
system is
tations of the active elements forming the
operational units. Summers, integrators and
d2 d
multipliers, for example, all have a maximum M --.JL + R J.- + Ky = 0 (14.5)
output level at which saturation occurs and at df dt
which they cease to behave as ideal elements
with well-defined characteristics. This second-order equation may be rewritten
Amplitude scaling of variables is a normal- as two first-order equations
ization process which is essential in analog
simulation. In an amplitude-scaled simulation (14.6)
all variables are converted to scaled quantities
which lie between +1 and -1. These scaled, or (14.7)
'per-unit' quantities, are normally expressed in
machine units (MU). In an analog computer,
or in a special-purpose electronic analog simu- where XI = Y and X 2 = dy / dt.
lation developed from electronic hardware In order to scale these two equations in
elements, a scaled value of +1 MU or -1 MU terms of amplitude we must have available
corresponds to the maximum permissible pos- the maximum value of the rate of change of X2
itive or negative output level of the electronic as well as the information given above.
units. Using modern integrated circuit hard- Estimation of this maximum acceleration is
ware this normally corresponds to a range of quite simple in this case since, from equation
+10 V to -10 V. For highest accuracy, the com- (3.5), it is clear that
puter variables should always have values
which are as close as possible to the permis-
sible maximum (1 MU positive or negative).
214 Real-time simulation
K
L-.;M~1
The condition for the maximum acceleration for the full range of conditions over which the
could logically be taken as the situation corre- simulation is to be used. To some extent
sponding to y and dy / dt both taking their scaling always involves some guesswork, but
maximum values. This would give the process may be made much more reliable
if estimates are based upon some preliminary
results obtained from a non-real-time digital
simulation.
For the maxima discussed above we now
This value of 80 m S-2 is, in fact, likely to be a have scaled variables which are as follows
generous estimate of the maximum accelera-
tion since, for physical reasons, it is unlikely X
XI
---
that the displacement and velocity will have Ise - 0.05
their maxima in this system at the same time
instant. It might well be more appropriate to The procedure normally adopted in scaling
examine the two contributions from the terms calculations involves first of all setting up an
involving y and dy / dt separately to give a unsealed analog simulation diagram. This is
maximum acceleration value of 50 m S-2. based upon the state equations, and is very
Possibly, a value between 50 and 80 m S-2 easily constructed in this case, as shown in
might be the most appropriate choice, say, Fig. 14.7. To incorporate scaling the state
60 m S-2. The precise values of maxima used equations must be rewritten in normalized
are not critical but should be chosen to allow form without altering the balance of the terms
Analog simulation techniques 215
on the left- and right-hand sides. This gives and
the following equations:
That is
and
. X d
-.-lrnax
xlsc - - -dtx lsc (14.10)
X 1max
That is . X d
-.-2max
x2sc - - -dtx2sc (14.11)
X 2max
. x2max
X lsc = X2sc -.- (14.8)
x1max Since the input to each of the integrators must
be the time derivative of the output of that
and integrator, it is necessary to introduce addi-
tional coefficients immediately before each
. R X K Xl max
integrator unit in the simulation diagram.
X 2sc = - M X 2sc -.X 2max
- - M X lsc -.- (14.9)
2max X 2max Indeed, one of the important rules of thumb
of analog computer programming is that a
In addition, it is necessary to relate the two coefficient unit must appear between all inte-
derivatives through appropriately scaled grators in a scaled patch diagram. The part of
equations. It is known that the scaled diagram corresponding to equa-
tions (14.10) and (14.11) is shown in Fig. 14.8.
Since coefficients can only be set to values
smaller than 1, it may be necessary to use a
and combination of a coefficient setting with a
summer or integrator gain factor if the
. dX 2 required values are to be obtained. For
X2 = -
dt example, a coefficient unit cannot be set to a
value of 2.5, but such a coefficient could be
and hence it follows that established in a simulation by using a coeffi-
cient unit set to a value of 0.25, along with a
summer or integrator input with a xl0 gain
factor. In the mass-spring-damper problem
Fig. 14.8 Scaling of variables at the inputs and outputs of the integrator units in Fig. 14.7.
In this case the coefficient units C1 and C, must be set to values equal to x1m,j X1m" and x2m ,j X'm"
respectively, in accordance with equations (14.10) and (14.11).
216 Real-time simulation
the substitution of values gives the patch One further important aspect of analog
diagram shown in Fig. 14.9, and it may be computer scaling concerns timescaling. Analog
seen that coefficient units and gain factors simulations may be made to run in a variety
have been used together in this way at the of different timescales simply by altering the
two integrator inputs. It should also be noted gains associated with all of the integrator
that there are two additional coefficient units units by the same factor. For example, increas-
in this diagram which did not appear in ing all of the integrator gains by 10 in a given
Fig. 14.7. These provide a means of setting simulation problem speeds up the simulation
initial conditions on the two integrators. by 10 times. The technique usually employed
Note that, if integrator units can accommo- for timescaling is to define a constant {3 such
date more than one input, the scaling calcula- that
tion may be considerably simplified, unless
dx 2 / dt is needed as an output variable. This
{3 = computer time = ~
simplification arises from the fact that there is
problem time t
no need to use the maximum value of the
acceleration in the scaling calculations if the
integrator involving X2 as output can carry out If this timescale factor is greater than I, the
summation as well as integration. computer is running slower than the real
Fig. 14.9 Complete scaled analog simulation diagram for the mass-spring-damper problem of Fig. 14.7.
It follows, from equations (14.8) to (14.11) with the parameter values given and for the estimated maxima,
that the coefficient unit settings are C] = 1.0000, C2 = 0.1200, C3 = 0.S333 and C4 = 0.5000. Note the use of
gain factors of 10 and 1000 at the integrator inputs to provide coefficient setting which are not greater
than 1.0000. Settings for the initial condition coefficients (Cs and Co) are the initial values of the scaled
state variables and will depend upon the experimental condition to be simulated.
Digital techniques 217
process, while if f3 is smaller than 1 the simu- mode. Modes are normally controlled manu-
lation is faster than real time. From the equa- ally by the user.
tion above it follows that
1
Sf' [f(t)] = F(s) = J
f(t)e- st dt t >0 (A1)
Unit step
s
1
~ (0)]
Y [ :;t] = s" F(s) - s"-I/(O)
M[S2Y (S)-SY(0)-
The transfer function of a linear system is One very important point about a transfer
defined as the ratio of the Laplace transform function model is that, when it is used as an
of the system output variable to the Laplace element within a description of a larger
transform of the system input variable, with system, there are constraints on the type of
all initial values taken to be zero. It is there- interactions which are allowed with other ele-
fore inherently an input-output form of ments of the system. More specifically, it is
system description and does not provide any essential to ensure that the device which the
information about the internal system struc- transfer function represents does not load the
ture or about any intermediate variables. source that provides the input signal.
For the mass-spring-damper system of Similarly, the output of the device must not
equation (A.9) in Appendix A the transfer be loaded by the element that follows it and
function is obtained from the transformed receives its output signal. Thus, in electrical
equation with zero initial conditions. This terminology, it may be said that an element
gives the following equation described by a transfer function ideally has
input from a source with zero (or low) inter-
Ms2y(S) + RsY(s) + KY(s) = F(s) (B.1) nal impedance and supplies its output to an
element which has infinite (or high) input
The transfer function is therefore given by impedance. If these conditions are not satis-
fied the transfer function model will be incor-
G(s) = yes) = 1 (B.2) rect. The process of ensuring that a transfer
F(s) MS2 +Rs+K function description satisfies appropriate
interaction constraints is essentially part of
The denominator of the transfer function, the process of establishing suitable bound-
when set equal to zero, is called the charac- aries for the sub model which the transfer
teristic equation, since the roots of this function represents. Further details of trans-
equation characterize the time response fer function modeling techniques may be
whatever the input forcing function. The found in many textbooks on dynamics and
roots of the characteristic equation are called control (e.g. ref. 1).
the poles of the system transfer function.
Roots of the equation formed by equating
the numerator polynomial to zero in the REFERENCE
same way are known as zeros of the system 1. Thaler, G.J. (1989) Automatic Control Systems,
transfer function. West Publishing, St Paul, MN, Chapter 2.
APPENDIXC
BLOCK DIAGRAMS AND SIGNAL FLOW
GRAPHS
yet) = Iru(t)
yes) = kU(s)
yet)
yet) = UI(t) + u,(t) + u,(t)
Yes)
yes) = UI(s) + U,(s) + U,(s)
c) Integrator
~J +1 t
yet) = yeO) u(r)dr
o
d) Traosfer function
Fig. C.l Examples of linear block diagram elements. Where appropriate both s-domain and time-domain
versions are given. The time-domain versions are also applicable in block diagrams involving nonlinear
functions and other nonlinear elements, such as multipliers and dividers.
and
(e.S)
(e.6)
The complete feedback loop may therefore be
the equations (e.4), (C.S) and (C.6) may be replaced by a single block having a transfer
combined to give function
Signal flow graphs 241
a) Combination of cascaded blocks
----8-
c) Elimination of a feedback loop
--1
Fig. C.2 Examples of simple block diagram operations.
U 6 Y
a) 0 ) 0 Y(s) = G(s)U(s)
Y1
Y,(s) G,(s)U(s)
c) Y2 Y,(s) = G,(s)U(s)
Y,(s) = G,(s)U(s)
Y3
(0.13)
f
(k+l)T
~
during each sample interval by the action of a
x(k + 1) e" x(k) + BU(k{J eMda) (0.10) zero-order hold device. This is essentially the
same as the hold function assumed in the
analysis of section 0.3. A zero-order hold is
where the simplest type of hold device and simply
maintains its output at the level of the preced-
a=(k+1)T-T (0.11) ing sample until the next sampling instant
T
uft) / ulll .. jlERO-ORDER yltl
-----., ee-----_. HOLD
Absolute error 56,96 Communication interval 96,131, Direct construction method 32-5,
ACSL xii, 31. 88, 92, 103-4, 106, llO, 219 38-9
111,125,126,127,179 Comparator 124,210 Discontinuity
Action potential 126,128 Compartmental analysis 5,176--9 detection problem 64-6
Active state (of muscle) 187 Complementary function 13 location problem 64--{i
Adams-Bashforth integration Computer state 63-4,131-8
method 61, 217-18 analog 13-17,123,209-17,225 time 64
Adams-Moulton integration method digital 2,13,15-18,208-9,221-2 Discrete even simulation xi, 2,
61 hybrid 15-16,223 48-50,131-8
AD-10 209 Contractile element (of muscle) 188, DISCRETE section 103,137
AD-100 209, 221-2 189-205 Discretization of continuous model
ADSIM 222 Control systems engineering 3, 13, 243-4
Advisory systems 231 21,48-50,121-5,131-8, Distributed time delay 33,41-3
Aircraft 21. 148-9, 163-73, 207, 224 153--{j2,223 Documentation 139
Algebraic loop 66, 88 Co-processor, numeric 18 DOS xii, 92, 247-9
Analog computer 14-16,123, COSMOS 227 Double precision 62
209-17,225 CRAY-l 208 Dynamic check 143-4
Analog simulation 13-17,209-17, Criterion function 144 DYNAMIC segment 88-91, 92, 95,
225 CSMP /358 85,88 104, 128-9, 131
Analog to digital converter 2,49-50, CSMP-1130 85,110 DYNAMO 233
131 CSSL-III 88
Array processor 208,221-2 CSSL-IV 88,101-3 EASY5 112
Artificial intelligence 19,230-3 CSSL'67 recommendations 87-8, Education 21,223, 225--{j
Automatic landing system, aircraft 91 Effort store 27
163-5, 169-73 Effort variable 27
Auxiliary inputs 6,147-9 D-operator 13 Electrical circuit simulation 5, 18,
Damping, viscous 13 117-20
Backlash 123-5 DAS 17 Empirical validation 145-50,160-2,
BCKLSH function 125 DASSL 87 181, 223
Biological systems xii, 8-10, 21, Delay Emulation 207
125-9,175-85,187-205,225 distributed 33,41-3 Engineering design 20-21
Biomedical engineering xii, 8-10, 21, pure 33,41,88,110 ENIAC 16
125-9,175-85,187-205,225 Delta operator 80 ENPORT 31
BIOPSI 110, 113 DERIVATIVE section 88,92,95, Equation
Block diagram xii, 23, 32-40, 110, 128-9, 131 difference 67-84,132
239-41 DESCTOP 105 ordinary differential xii, 5, 13, 14,
Bond graph 23,27-31,66,233 Design of experiments 20,147,149 23
Boxer-Thaler operator 76 DESIRE 104-5 partial differential 14,41,42
Differential analyzer 14 Equation-oriented methods xii,
CAMP-G 31 Differential equations 85-106
Causality 28 ordinary 4,5,13,14,23 Errormeasures 56--7,61-2,96,220
Chemical process simulation 18,63, partial 14, 41-3 ESL 227
225 Digital control system 2, 48-50, Euler's integration method 52--{j,
CINTERVAL 96,167 131-8 218
Coefficient unit 212 Digital to analog converter 49-50, 13 Experimental design 20,147,149
252 Index
Expert knowledge 147,223 Integration methods continuous-variable 2
Expert system 19,231 Adams-Bashforth method 61, discrete-event 2
Explicit integration methods 52-63 217-18 discrete-time 243-6
External validation 144-50, 160-2, Adams-Moulton method 61 distributed parameter 14,40-43
181, 223 Euler's method 52-6 linear 3
explicit methods 52-63 linearized 21,148,156-7
Final value theorem, of Laplace Fehlberg method 62 lumped parameter 13,14,16,
transforms 236 first difference method 76-83 23
First-difference operator 76-83 fixed step methods 52-62 nonlinear 3,13, 18,63-6,70-74,
Flight control system, aircraft 21, for stiff systems 62-3 121-33, 153-62, 169-74,
148-9, 163-73 Gear's method 63 175-85, 187-205
Flight simulator 8, 16,207,224 implicit methods 52 of data 3,148
Flow store 27 improved Euler method 57-8 of a digital computer 2,48-50,
Flow variable 27 multistep methods 60-61 131-8
Fluid-power system simulation 18 one-step methods 52-60 of systems 3
Fortran xiii, 16, 69, 85, 86, 92, 111 operational methods 75-83 reduced form 24
Frame based systems 231-2 predictor-corrector method 60-61, state space 26,74-5
Frequency response of integrator 83 76 statistical 3
Friction Ralston's method 59 stiff 62-3
static 13 Runge-Kutta methods 56-60,61, time invariant 3
viscous 13 62,76 time series 148
Fuzzy knowledge 231 trapezoidal method 60 time varying 3,63-6,81-3,176-9
variable step methods 62-3, 218 MATLAB 86,111-12,114
Gamma operator 80 Integration operators 75-83 MATRIX 86,111
Gas-exchange, pulmonary 5,8-10, Integration step size selection 51, Matrix packages 86-7
175-85 144,219 Medical applications xii, 8-10, 21
Graphical interface 109,229 Integrator 110,211-12 Menu 109
Gyrator 29-31 Integrator frequency response 83 MERROR 96
INTERACT 222 MICRODARE 105
Halijak operator 76 Interactive simulation 17,19, MIDAS 17
Hardware-in-the-Ioop testing 207, 207-26,228-9 MIMESIS 91
222-3 Internal verification 142-4 MIMIC 17
Helicopter 21, 149 International Association for MINTERV AL 96
Hold Mathematics and Computers Mode 217
triangular 49, 67, 68-9 in Simulation (IMACS) 227 Model
zero-order 67,68 Iterative construction method 36-7, applicability 5
Hold mode 217 40 assumptions 5
Hybrid computer 15-16,223 boundaries 4, 6, 10
Hybrid simulation 15-16,223 Junction, bond graph 28-31 calibration 6, 8
Hydraulic system simulation 18, discrete-time form 246-8
121-5, 153-62 Knowledge representation 231-2 evaluation 10
Hypothesis testing 2, 8, 18-20 linearization 21, 156-7
Hysteresis 13,110 Laplace transform 13,26,32,156, of mechanical system 24, 26, 32,
235-6 121-5,213-14,223
Icon 109 Lateral beam guidance system, of systems 3
Identifiability 148 aircraft 165, 169--73 optimization 18
Identification 8, 148-50 Limit cycle 125 qualitative 232-3
ILLIAC-IV 208 Linearization 156-7 sampled-data form 67-84,243-6
Implicit integration 52 Lung 5,8-10,175-85 sensitivity 18
Implicit loop 66, 88 simplification 6
Improved Euler integration 57-8 MACPUF 21, 225 state space form 26,74-5
IMSL subprogram library 16 Macro facility 103 state-variable form 24-6, 32-40,
Inertia 27 Madwed operator 76 74-5
Inheritance 231-2 Man-machine interface 109,228-9 statistical 3
Initial condition mode 217 Mass-spring-damper system 24-6, testing xii, 6-10, 141-2
INITIAL segment 88-91, 92 213-14 transfer function form 3, 23, 26,
Instability, numerical 51,61-2, Mathematical models 32-40,76-81,165-73,237
219 'black box' 3 validation xii, 6-10, 43,141-51,
INTEG function 95 compartmental 5 160-2, 181
Index 253
Model-method-experiment concept Potentiometer 212 Simulation of transfer function
228 Predictor-corrector integration by direct construction 32-5,38-9
Mouse 17,109 method 76 by iterative construction 36-7, 40
Multiple instruction-multiple data Project Whirlwind 16 by parallel construction 35-36,
(MIMO) 221,229-30 PSI 110 39-40
Multiplier 110, 213 Pure time delay 33,41-3 Simulator 8, 16,21,223-6
Multiprocessor systems 208,221-2 SIMULINK 86, 111-12, 114, 119-20,
Multirate methods 220-1 Qualitative simulation 232-3 125,137-8,181,210
Multistep integration methods Single instruction-multiple data
60-61 Real-time simulation xii, 15, 16, 18, (SIMO) 221
Muscle spindle receptor 125 207-26 Single instruction-single data (SISO)
Reduced form 24 221
NAG subprogram library xiii,16 Relative error 56,96 Single step integration methods
Neural encoder 125-9 Respiratory control system xii, 8-10, 52-60
Neuromuscular control system 125 175-6 Skeletal muscle 187-205
Normalized variables 213-16, Robustness issue, in model SLIM language
219-20 validation 149-50 examples 97-101,118-20,125-9,
NOSORT 90 Round-off error 61-2 129-31. 131-7, 157-9, 167-9,
Numerical instability 51,61-2 Rule-based systems 19,231 171-3, 179-85
Numerical methods Runge-Kutta integration methods hardware requirements 247
discontinuities 63-6 76,87,218 language features 91-2,106,228
integration 51-63 Runtime commands 103 output facilities and graphics
round-off errors 61-2 96-7,251
stability 61-2 S component 28 software installation 247-8
subprograms 16 Sampled-data models 67-84,243-6 syntax 92-6
truncation errors 61-2 Sampling process 244-6 SLIMPLOT program 96-7,100,
'Numerical Recipes' subprograms Scaling calculations 247-9
16 amplitude 213-16,219-20 SmallTalk 232
time 216-17, 225 Society for Computer Simulation 17,
Object-oriented systems 231-2 Sensitivity analysis 18,43-8,98-202 87, 142, 150
ODEPACK 87 Sensitivity coefficient 43, 44-8, Software engineering 233
Operational amplifier 14,209-12 196-202 SORT 90
Operational methods 76-83,218 Sensitivity co-system 44-8, 196-8 Sorting 89-90,103
Operators 76-83,218 Sensitivity function 43, 44-8, Source 27
Optimization 146 196-202 Stability of integration methods
Series elastic element (of muscle) 61-2
Package, simulation 16,85,87 188, 189-205 Standards 17,87-91,227
Pade approximation 42-3 S-function 112 State space 26,74-5
Parallel construction method 35-6, Signal (bond graph) 30 State variable 24-6,32-40,74-5
39-40 Signal flow diagram 23, 32-40, Static check 143
Parallel elastic element (of muscle) 241-2 Static friction 13
188, 189-205 SIMNON 110,111 Steady state finder 143
Parallelism 208,210,221-2,229-30 SIMsystem 222 Stiff models 62-3
Parallel processing 146,208,221, Simulation Submodel 6,28,146,227
229-30 analog 13-17,217-19,225 Summing unit 110,211
Parameter sensitivity analysis 18, block diagram-oriented xii, 17, Superposition principle 3
43-8,146 110-16 SYSMOO 227, 228
PARMACS 230 discrete-event xi, 2, 131-8 SYSTEM 100, 222
Parsimony principle 6 equation-oriented xii, 85-106 SYSTEM BUILD 111-12
Partial system testing 207, 222-3 hybrid 15-16,223 System dynamics methods 233
Particular integral 13 interactive 17,19,207-26,228-9 System identification 8,148-50
Pascal 16,69,86 languages 17,87-106,227-8
Pendulum 129-31 packages 16,85,87 Tanks, coupled 153-62
Peripheral array processor 208, qualitative 232-3 Taylor series 43,53,56-60
221-2 real time xii, 15, 16, 18, 207-26 TERMINAL segment 88-91,92
Physiological systems xii, 8-10, 21, run 18 Testing xii, 6-10, 141-2
175-85, 187-205 standards 17,87-91,227 Tetanus (of muscle) 188
Pipelining 208, 221 study 18 Time delay, distributed 33,41-3
Post-processing 229 Simulation Councils Inc. 17,87 Time delay, pure 33,41-3
254 Index
Training simulators 8, 16, 21, 223-5 empirical 145-9 Water turbine speed control 121-5,
Transfer function 3,23,26,32-40, external 144-50, 160-2, 181 223
76-81,218,237 robustness 149-50 Whirlwind Project 16
Transformer 29-31 terminology 142 WIMP interface 109
Transputer 222 Validity Windows environment 17,229
Trapezoidal integration 60 empirical 145-9,223
Triangular hold 67, 68 heuristic 144 XANALOG 111,112-13,114-16,
Truncation errors 61-2 pragmatic 144 119-20, 125, 127, 129, 137-8,
Tustin operator 76, 80, 83 theoretical 144-5 157-8, 181, 193-5, 210
Twitch (of muscle) 187 Van der Pol's equation 97-101,113
Variable-step integration 62-3, XERROR 96
User interface 109,228-29 218
VAX/VMS 222 Z-transform analysis 67-84,136-7,
Validation Verification 142-4 245
documentation 150-1 Viscous damping 13 Zero-order hold 49,244-6
Chapman & Hall specifically disclaim all warranties, express or
implied including, but not limited to, implied warranties of
merchantability and fitness for a particular purpose. With respect to
the disc and the programs reproduced on the disc, Chapman & Hall
shall have no liability with respect to any loss or damage directly or
indirectly arising from the use of the disc or any program reproduced
on the disc. Without limiting the foregoing, Chapman & Hall snall not
be liable for any loss of profit, interruption of business, damage of
equipment or data, interruption of operations or any other commercial
damage, including but not limited to, direct, indirect, special,
incidental, consequential or other damages.