"‘V-‘\v.vvv

 

’4‘-.'('~L”".j.' 2" $2.. 'J. "1: «‘3': .‘ m . A A f .. T-;’§\.u: \ .. ‘ A. _ ’1“: .
‘ I..." 3“...“ . 2': ’I ’ " - '- ’7 A? %‘ ;\'¢\ . "(1 -. ‘ A _ 39‘ $ ~ $X¢ zJ§ o
"I. ’ . . .

. t.
, jig? :{h “I...“

I;

55“..”
‘ 7.

'

q '- 1;.
Wm

v

§
I

.- ? W

V. akin/3y ‘-
. '3“ "0-7 ‘9.
Q‘110 L._ft_v '1?

i
A‘;:': ‘ \‘§\:‘;3
r \ e V. “ i mgfrnif
0.5:. . ”M75“

‘.‘ '. '. ‘ ‘~‘V nk‘3ffik "L J’V’Q‘M
2 5;? 5335. -,,,.L :4 '
n ‘J. L 5:? E ‘ A
. 2 In} ' p ’ 9‘ -~ ‘1:ch . $4.3. .Sé. ”3:1; ‘.'.‘
' V :4: “6.319.”. 29"?“ fit - {V g?” \‘.::“r "9.x *{g , kiln-3 “M: iv“) 5,0.-
‘M I :O'J? fit. 5: g V. i:\‘-‘.‘ p‘ v ‘ Q.
..'-U’ '\ A’VJJHA _ {1" I V?“ ‘ A :‘ ‘v‘.
‘1 I I.

.7: J
“~21 y

‘- ""
a. u w
‘F-‘AN' a?“ n
'3" 16 at. "
.23 .3, .-. .
i ,tfqé ’2‘”
2:» ~ > r

'3»: k321i
(1%" 7%.
ifékh' " «'13.,
«m
’6)“- 5‘455} V

”11%;:1‘ . gm
.1" I 1')”. 3.1!}
3253}: ”fig?

-"" I

,3
..

V
‘lvlef'.
5‘ I
N ..
g3 35. .
I

.
."

kI‘m":
I‘l

l' :‘élé':

(521';
?.h:-¢KA
mt. Mia

. . l:'.'( i:{ I
”‘1‘ 1:00,}; " r

- A( ‘er ‘ I 3' ...
”*2 Egg-+23% w:
’7‘ f’ f. *5 ' .
\ fi‘ig‘ fl$?: g
} J #1. “I;

x»
a~'f: ' . h.
;' h’h'fi:

 

. - '. A “O

”(:53 31.54? 97; .3
figs". :I 3' 4.0 ,'

‘.'.‘_

lidOBqu

MICHIGAN STATEU

HWMMWMW

300563 9558

 

 

 

 

 

 

 

 

 

 

 

 

 

mum [II

This is to certify that the

dissertation entitled

Improving Design Assessment and Simuiation
of Large-Scale Dynamic Systems

presented by

Tong Zhou

has been accepted towards fulfillment
of the requirements for

Ph. D degree in Engineering

Major professor

{éfiwm/
gym/Jo" /

MSUi: an Affirmatiw Action/Equal Opportunity Institution 0- 12771

 

LIBRARY
Michigan State

 

 

 

 

 

 

University i
MSU RETURNING MATERIALS:
Piace in book drop to
“saunas remove this checkout from
“ your record. FINES will

 

 

be charged if book is
returned after the date
stamped beiow.

 

 

 

IMPROVING DESIGN ASSESSMENT AND SIMULATION

OF LARGE-SCALE DYNAMIC SYSTEMS

By

Tong Zhou

A DISSERTATION

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

DOCTOR OF PHILOSOPHY

Department of Mechanical Engineering

1988

1’. L: (0

(0'7 (4)

ABSTRACT

IMPROVING DESIGN ASSESSMENT AND SIMULATION

OF LARGE-SCALE DYNAMIC SYSTEMS

By

Tong Zhou

When accurate models of complex large-scale engineering systems are
made, they may involve large numbers of coupled nonlineam'ciifferential
and algebraic equations. It can be both time-consuming and difficult for
the design engineers to get a good assessment of the performance of the
design, especially in a way leading to design improvements. The
objective of this research is to find methods that improve the design
assessment and simulation of complex. large-scale dynamic models,

thereby reducing time spent to gain increasing insight by the engineers.

With respect to simulation efficiency.tfim:problem of implicit R
fields is very important to the bond graph modeling of engineering
systems. Simulation efficiency in such problems may be humeasedinz
reducing the number of iteration variablefix :% new algorithm for
determining the minimum number of iteration variables required in a
model with implicit R fields is presented. The algorithm generates exact
results for implicit R fields containing l-port R nodes and weighted
junction stnunnues (l, O,’TF). An extension is made for causal IRFs

containing multiport R nodes to get the minimum number<fiiiteration

variables. As a further extension the basis order properties of general
junction structures (1, O, TF, GY) were derived by using gyrographs and

a maximum matching algorithm.

To increase insight a new approach to model order reduction and
simplification has been developed in the bond graph framework. A pilot
version has been implemented in software. Power responses measured in
various ways are displayed directly on the bond graph by color coding.
Such a display gives engineers an easy way to see the power
distributions in a large-scale dynamic system under various sets of'
operating conditions. By interpreting this power-based data suitably,
possibilities for model order reduction and simplification can be
identified without detailed analysis or equation manipulation. The power
response data is available for nonlinear models as well as for linear
ones, so the power-based reduction approach may be used for nonlinear

time-varying systems.

The results obtained contribute to the design assessment and
simulation problem by increasing simulation efficiency for certain
classes of models and by increasing the insight about system behavior
available to the design engineers. Further directions of work are

discussed based on the results presented.

Dedicated to

My parents, Wenjin Zhou, MD., Hangxiang Lee, MD..
My wife Shuling Jiang. and

My son Fan

iv

ACKNOWLEDGMENTS

The author is very grateful to Dr. Ronald C. Rosenberg for his
guidance,enmeuragement, persistence, and friendship in his role as
major advisor, critic, taskmaster, and friend. The author also would
like to thank Dr. Hasan Khalil, Dr. Clark Radcliffe, and Dr. Steve
Dragosh for their'contribution throughout his graduate study and
research in both class teaching and serving on his doctoral cxnmnittee.
Many thanks go to Mr. Dave Withers of the IBM Information Systems for
his interest and support in the research and for serving on the

author's guidance committee.

The support of the IBM Thomas J. Watson Research Center is greatly
appreciated. The research was also supported by the use of the excellent
computing facilities and the aid of the cooperative personnel at the A.
H. Case Center for Computer-aided Engineering in the College of

Engineering at Michigan State University.

TABLE OF CONTENTS

LIST OF TABLES ....................................................

LIST OF FIGURES ...................................................

Chapter 1 INTRODUCTION ...........................................

1.1 Problem Statement ........................................

1.2 Research Objectives ......................................

1.3 Dissertation Organization ................................
Chapter 2 EFFICIENT COMPUTATION OF IMPLICIT DISSIPATION FIELDS

2.1 Problem Definition .......................................

2.2 General Solution by Iterative Methods ....................

2.3 IRFs with One-Port R Elements .............................

2.4 IRFs with Multiport R Elements ............................

2.4.1 Multiport resistances in physical systems ............

2.4.2 Extension of the Algorithm to IRFs
with multiport R nodes ...............................

2.4.3 Assigning causality to an IRF Containing
One Multiport R Node .................................

2.5 IRFs with General Junction Structure .....................
2.5.1 General Junction Structure ...........................

2.5.2 Gyrographs and Maximum Matching ......................

2.5.3 Basis Order Algorithm for General Junction Structures

2.5.4 Examples .............................................

2.5.5 Remarks on the Algorithm .............................
Chapter 3 MODEL ORDER REDUCTION METHODS ..........................
3.1 Aggregation Methods ......................................

3.2 Modal Analysis ...........................................

3.3 Lyapunov Function Method .................................

vi

12
19
25

42

59
66

77
79

8O
82
85

3.4 Perturbation Methods .....................................
3.4.1 Regular perturbation method ..........................
3.4.2 Singular perturbation method .........................

3.5 Component Connection Methods .............................
3.5.1 Component connection method ..........................
3.5.2 Diakoptics ...........................................

3.6 Component Cost Analysis of Large Scale Systems ...........

3.7 Model Order Reduction in Bond Graph Models ...............
3.7.1 Sequential causality assignment method ...............
3.7.2 Reciprocal system method .............................

3.7.3 Finite mode distributed system models ................

Chapter 4 POWER PORTRAITS OF DYNAMIC SYSTEM MODELS ..............

4.1 Graphical Representations of Dynamic System Models ........
4.2 Uniform Performance Measure --- Power ....................

4.3 Power/energy Visualization ...............................

b

.3.1 Computation of power variables .......................
4.3.2 Power measures .......................................
4.3.3 Scaling and color-coding .............................

4.4 Power Method for Model Order Reduction and Simplification
4.4.1 Simplification of physical effects ...................
4.4.2 Decomposition and simplification of subsystems .......
4.4.3 Procedure for model simplification ...................

4.5 Examples .................................................
4.5.1 A radar pedestal unit ................................

4.5.2 An Euler-Bernoulli beam with a vibratory load ........

Chapter 5 SUMMARY AND CONCLUSIONS ................................

5.1 Implicit R-field Simulation ..............................

5.2 Power Distribution and Display ...........................

86
87

89
89
91
92

94
95

96
98
98

105
105

110
113

113

117

137
139

APPENDIX
Appendix Standard Bond Graphs and Gyrobondgraphs .............. 142
1. Standard bond graphs
2. Gyrobondgraphs
3. Gyrographs

BIBLIOGRAPHY 145

viii

 

LIST OF TABLES

Table page
2-1 Possible variable selection ............................. 38
4-1 Powers in various energy domains ....................... 104
4-2 Parameters of pedestal model ........................... 121
4-3 Effect indices ......................................... 122
4-4 Eigenvalues of the original and reduced models ......... 125
4-5 The power listing of five-mode model ................... 130
4-6 The powers on cut bonds ................................ 132
4-7 The power listing of two-mode model .................... 134
4-8 Eigenvalues of the original and reduced models ......... 134
4-9 Computer times ......................................... 134
A-l Standard node sets ..................................... 144

ix

LIST OF FIGURES

Figure page

2-1 A mechanical system with dissipative coupling ............ 7
(a) Physical system
(b) Bond graph model
(c) Causality assignment

2-2 Multiport Field Structure ................................ 13
2-3 Key vectors .............................................. 13
2-4 A system containing two IRFs ............................. 18

(a) The identified IRFs
(b) Causality of the IRFs

2-5 An example of an IRF ..................................... 24
(a) Bond graph of the R-field
(b) Causality assignment

2—6 Four-way-valve amplifier ................................. 26
2-7 Reduced form of hydraulic amplifier ...................... 26
2-8 Power transistor and its bond graph model ................ 27
2-9 Key vectors in an IRF .................................... 28
2-10 Functional diagram of an TRF ............................. 32
2-11 Variable flow diagram .................................... 33
2-12 Loop formed by an improper vortex e change ............... 36
2-13 Inclusion of l-port R ncdts ........ .. .................. 37
2—14 An electronic circuit and its bond graph model ........... 39
2-15 Acausal IRF with multiport R node ........................ 42
2-16 An example of self-loop and R-block ...................... 43
2-17 An R-block of type 1 ..................................... 45
2-18 An R-block of type 2 ..................................... 46
2-19 An R-block of type 3 ..................................... 46
2-20 An R-block of type 4 ..................................... 47
2-21 A multiport R node without structural coupling ........... 51

2-22

2-23

2-24

2-25

2-26

2-27

2-28

2-29

3-1

3—2

Example 1 ................................................
Example 2 ................................................

Graphs of an R-field .....................................
(a) The standard bond graph

(b) The gyrobondgraph

(c) The gyrograph

Causality in different graphs ............................
(a) In a standard bond graph

(b) In a gyrobondgraph

(c) In a gyrograph

A graph and one of its matchings .........................

An alternating path ......................................
(a) Before reversing markings
(b) After reversing markings

An alternating cycle .....................................
(a) Before reversing markings
(b) After reversing markings

Example 1 ................................................
(a) Bond graph and gyrobondgraph
(b) Gyrograph

(c) Maximum matching in GJ

(d) Maximum matching in CC
(e) Causality corresponding to the F-minimum basis
(f) Causality corresponding to the F-maximum basis

Example 2 ................................................

(a) Gyrograph

(b) Maximum matching of OJ

(c) Maximum matching of CG

Reciprocal system method .................................
Finite modal bond graph model ............................
Graphical representations of mathematical modeling .......
Powers in system models ..................................
(a) schematic diagram

(b) block diagram

(c) bond graph

Power responses of the system in Figure 4-2 ..............

Powers on external bonds .................................

Powers on cut bonds ......................................

xi

54
56

62

64

66
69

69

75

76

95
97
99
103

4-6

4-7

4-8

A radar pedestal unit .................................... 120

Bond graph for the radar pedestal unit ................... 120
Power display on the bond graph of the pedestal .......... 124
Reduced-order model of the pedestal ...................... 125
Responses of the original and reduced models ............. 125
A beam-coupled system .................................... 128
Modal bond graph model for the beam coupled system ....... 128
Power distribution in the modal bond graph model ......... 131
Reduced model with two modes ............................. 133
Comparison of input and load velocities .................. 135

Chapter 1
INTRODUCTION
1.1. Problem Statement
Suppose we have an energy-based structural model (i.e., a bond

graph model) of a large-scale, complex, nonlinear electro-mechanical
system, including its controls. In addition we have a set of nominal
constitutive equations for the energy/power elements and a set of input
functions for the source elements. A simulation can be run to estxflolish
the nominal system response under the conditions noted. We assume that
the response so obtained is satisfactory in terms of required system
performance. This situation represents the nominal design (Rosenberg,

1985).

Now we ask the following question: How much change can be tolerated
in a set of parameters associated with the constitutive equations and
still have the system response meet the performance criterion
acceptably? It is necessary to define an acceptable performance
criterion C quantitatively to attack the question meaningfully. We

assume this has been done in terms of a set of system output variables

y(t).

The most obvious method of attack is to choose a set of values for
a parameter vector P, run a simulation, and assess the resulting y(t)
with respect to the performance criterion C. There is a statistical
version of the problem in which statistics are assigned to the P vector,

and the statistics of y(t) with respect to C are generated. The

0

(—
statistical approach is of great value in making cost/benefit analyses
associated with manufacturing and assembly decisions (Prakash, et. a1,

1985).

One deterministic version of the problem is based on worst-case
analysis. We try to find.the lindts on the values for P that push the
response y(t) to the acceptable limits of C. In another version we try
to guarantee that choosing P values with certain boundswdjl meetCZ
acceptably. The first version is the best we can hope to do with a
deterministic approach, given C in an acceptable/unacceptable sense (as
contrasted with a cost function sense). Since setting error tolerance on
performance of many classes of electroqmufimndcal equipment is a common
practice, we will assume that C is stated in an acceptable/unacceptable

sense .

Progress in treating the problem described above will benefit
engineers in discovering and eliminating failure modes characterized by
exceeding dynamic tolerances within the design cycle. Therefore we can
expect to achieve greater reliability and manufacturability of new

products.

1.2 Research Objectives

The particular research objectives are derived from the problem
described above. Given a fairly general class of nonlinear dynamic
systems that are physically (i.e., energetically) based, we expect to
use a numerically-implemented approach, rather than functional analysis,
as the principal tool in improving practical assessmenttxmhniques.

Included in numerical computation are tools for studying thereffects of

3

variations in model topology as well as parameter set values<n1the
response y(t) relative to C. The principal goal is to find‘ways both to
reduce computational effort relative to deterministic modeling and to

guide the search for model alternatives automatically.

The main problem of interest, which is related to the efficiency of
simulation, but more concerned with the solution of the system:
equations, is that of coupled.nonlinear algebraic equations that arise
when coupled bond graph R elements exist. The implicJJ: (coupled)
nonlinear algebraic equations are often difficult to solve, and they
typically must be solved iteratively at 6%Kfll integration step. The

primary objectives in this sub-problem are:

l. to develop an automatable method to identify and solve the
coupled nonlinear algebraic equations correctly; and
2. to study possible improvements in the existing solution process

such that the computation cost can be reduced.

In order to make the design process more efficient and
computationally economical, one avenue of approach is to reduce the
order of the system model, since the computation effort typically
increases much faster than the system dimension does, and high-dimension
systems lead to more complicated control design. Since no general
theoretical method is available for assessing the response of nonlinear,
large-scale complex engineering systems. a need arose to find a new
method to consider the possibility of reducing the order of the dynamic
system model. One of the research objectians is to provide a general
purpose nunwrical approach to model order reduction and simplification.

The application of this approach should not be restrictxui'by the

4

linearity of the dynamic systems. Also, it should be easy to use and
should provide engineers additional information about the system under

study.

1.3 Dissertation Organization

In pursuit of these objectives the following subjects will be
discussed in the thesis. The dissertation consists of two parts:
improved solution methods for implicit dissipation fields and
power potraits of dynamic system models. Although both problems are
related to the improvement of simulation and design efficiency, they are
rather different in the nature. Therefore, we will discuss then:

separately in the following chapters.

Chapter 2 deals with the implicit dissipatitnl field problem. It
presents results concerning the efficient cmmnmetion of duacoupled
algebraic equations by an iteration process. A procedure for determining
the mirdJmun number of (independent) iteration variables for implicit R
fields with simple and weighted junction structures is given. For
implicit R fields with general junction structures an algorithm for

finding the minimum and maximum flow input variables is described.

In Chapter 3 a brief review of some recent work in the area of
model order reduction is given. This chapter provides a background
relative to the general model order reduction problem and analytically-

oriented methods.

Chapter 4 introduces the power concvpt and its application in bond

graph models. A new tool for exploiting tile pntwn' attrilnites of

5

a dynamic system is presented. Some examples are given to illustrate the

potential scope and utility of this method.

In Chapter 5 the results for improving the design assessment tools

are summarized and some future research directions are identified.

The computer simulations used in this research were run using the
ENPORT-7 bond graph/block diagram software. The program was version 7.1
run cum a PRIME 750 under PRIMOS. In several cases noted in the text the
algorithms developed in this research were implemented within the ENPORT

framework on a pilot basis.

6

Chapter 2

EFFICIENT COMPUTATION OF IMPLICIT DISSIPATION FIELDS

2.1 Problem Definition

One of the major objectives of improving design assessment is to
automate the process of formulating and solving the state equations
associated with bond graph models of engineering systems (See Appendix
A). A second objective is to provide timely and insightful feedback to
the designer. There are several sources of difficulty in accomplishing
the major objective. The one we wish to focus on here is that of solving
coupled nonlinear algebraic equations that arise when the R elements are
connected in particular ways in the model.7fluainmlicit equations are
often difficult to solve, and they typically must be solved several
times in each integration step. Consequently, it is helpful to be able
to inform the modeler in detail of the existence of such coupling.
Furthermore, increasing the efficiency with which such solutions awe

obtained can dramatically decrease the overall solution time.

The existence of algebraic loops in the equations of a physical
system may not be detected until the sorting or reducing process starts
in most traditional simulation approaches.lhu:their existence can be
verified even before equation formulation when the bond graph approach
is employed. To illustrate the case let us first consider a physical
device shown in Figure 2-1a. In Figure 2-lb. the corresponding bond
graph model has been built. The 1 element represents the inertial effect
and the compliance effect is indicated by a C element in the mechanical

system. The R elements represent energy dissipative effects. The SF

 

 

 

 

 

 

 

 

 

(a)

1 1
7“ f‘
\ k
C Ra C Rh

( b ) ( c )

Figure 2-1 A mechanical system with dissipative coupling
(a) Physical system
(b) Bond graph model

(c) Causality assignment

8

element indicates an imposed velocity on the left plate (assumed
massless) as an input. It is desired that the state-space equations be

obtained in an explicit form as follows

p1 ‘ 81(P1. 42. V0) (2.1)

QZ a 82(p19 Q2. V0) (2.2)

where p1 is the momentum of mass I and q2 is the deflection of spring C.

Causality can be assigned to the bond graph of Figure 2-lb
according to the general rules (Rosenberg and Karnopp, 1983). After
finishing the first step (assigning required causality to the source SF)
and the second step (assigning the integral causality to the storage
elements C and I), we find that the causality does not fully extend
through the graph. Some acausal bonds (bonds 4, 5, and 6) will be left.
At this stage, we realize that an R-field exists irllflnis system. This

implies that there will be an algebraic loop in the system equations.

Suppose we continue the causality assignment by imposing an

arbitrary causal orientation on one of the two R elements, say, R Then

4'
we extend the causal implication through the graph using the constraint
elements (0 and 1). Now the causality assignment has been completed

(Figure 2-lc). The state vector X and input vector U are identified as

follows

If we define F4 and V6 as auxiliary variables, then the auxiliary

equations are readily obtained as

 

p1 - F4 (2.3a)

a2 - V6 (2.3b)

and the constitutive equations are

F4 ’ 84(V4) B 84(V0' Pl/ml' V6) (2-43)

Assume that both R4 and R6 are linear, that is,

F4 - Rava (2.5a)

v - R' r (2.5b)

After some algebraic manipulations to eliminate the auxiliary
variables F4 and V6, an explicit state-space equation set can be
developed; namely,

p1 - -[R4R6/<R4+R6)mlpl + [Rak/(R4+R6)lq2 + [R4R6/<R4+R6>1v0 (2.6a)

qz - -[R4/<R,+R6>mlp1 - [k/(R,+R6>1q2 + [R,/<R,+R6>lvo (2.6b>

Now suppose that the RA and R6 are nonlinear. Then we may have
difficulty solving the auxiliary equations to get an explicit state
form. In general explicit analytic solutions of nonlinear coupled

equations are difficult, if not impossible, to achieve.

From the development above, we see that the process of causality
assignment is an.aid.in the process of identifying the algebraic loops
in dynamic systems. In Figure 2-1, causality on bonds 4” ES, and 6 can
not be determined after assigning causality to source and energy storage

elements. Furthermore, algebraic loops in the matimnmatical sense are

10

physically related.to the existence of dissipation fields. In Figure 2-
1, this dissipative field is consisting of nodes Ra"Rb’ and associated
junction structure. Such fields are called implicit R-fields (IRFs).
Reading the partial causality-assigned bond graph of Figure 2-1b, we can
easily identify the R fields from other dynamic fields. It is also clear
that nonlinear algebraic loops in system equations may prevent the
subsequent reduction of the equations to an explicit state-space form.

They may make system simulation very difficult to accomplish.

The solution of algebraic loops may not be very important from a
theoretical point of view. However, especially for nonlinear systems
'where an analytical solution is not always possible, solving the loops
is very computer-time consuming. Often when the model involves implicit:
algebraic equations, the approach is to use iterative solution methods.
These methods can be very costly and they nun; introduce difficulties

related to the existence and uniqueness of solutions.

Several researchers have been working in a bond graph environment,
using three different approaches to these problems. Barreto and Lefevre
(1985) try to avoid hmfljcit algebraic loops by modifying the system
model. Their basic idea is to consider the implicit part of the model
before attempting the simulation and to modify the model suitably so it
can be simulated using only explicit methods.7fimfli'proposed ways to
modify the model include: (1) imposing restrictions in the set of values

of admissible solution; (2) changing the model to fit reality.

The second direction is to improve the cxNMMJtational efficiency
within the implicit solution framework. Lorenz and Wolper (1985) made an

observation on the causality assignment in the case of algebraic loops.

11
'Phe algebraic loops are found by a related signal graph and the minimum
number of independent variables of the algebraic loops is the minimum
number of nodes required to break the topological loops. Two general

rules are suggested for assigning causality to the implicit R-fields.

The third direction in which the parasitic physical elements are
added to eliminate the nonlinear coupled R-field has been employed by
Zhou (1984). The method based on bond graph augmentation converts an IRF
into a dynamic subsystem that exhibits the proper static characteristics
at steady-state and employs a two-time-scale integration technique. This
method has a philosophy similar to that of the charge-up method in the
ASTAP program (Zeid, 1985). The charge-up method has proven to be very
reliable, but it is computationally costly. Since all capacitors and
inductors introduced into the original system are assigned a value of 1,
for a typical nonlinear circuit, it would take close to 2OON passes to
arrive at a solution, where N is the number of dynamic elements
introduced. In Zhou's work, an augmentation sequence and a general rule
for parameter selection for arbitrary n-th order subsystems have been
suggested. These have been numerically tested in several cases. The
algorithm appears promising, but it has not been optimized. Granda
(1984) also proposed several approaches to algebraic loops, including

adding a storage element.

In this study we assume that a system bond graph containing one or
more implicit R fields is given. The task is to solve the coupled
algebraic equations efficiently during the system simulation. The
decomposition technique is applied to the IRFs and the relevant
associated vectors. A direct bond-graph-based algorithm for identifying

the minimum number of iteration variables in certain implicit R-field

12
problem has been developed (l-lood, et a1., 1987). An advantage of the
algorithm is that it works directly with the bond graph itself. This
algorithm is also extended to IRFs containing multiport R elements. For
an IRF with a general junction structure, a preliminary result on the
basis order is developed which facilitates the efficient solution for

simulation of system containing such IRFs.

2.2 General Solution by Iterative Methods

A bond graph model may characterized by the diagram in Figure 2-2
(Rosenberg and Karnopp, 1985). The system inputs are defined by the
collection of Se and Sf elements and are referred to as the source
field. Dissipative effects from the R elements of bond graph are
grouped into the dissipation field. Dynamic effects are the result of C
and I elements in the model and are represented in the energy storage
field. Eacil<xf these fields has constitutive equations associated with
it. They are coupled by a power-conserving connective multiport
represented by the junction structure.7flu2Peynter junction structure
consists of O and 1 junction elements, (It is named after H. M. Paynter,
the inventor of the bond graph.) and it is invariant. The transducer
field is the collection of the transformers and gyrators, which may have
varying umnhrli. The key vectors of each field are identified in Figure
2-3, where U is the input vector, a function of time; X is the energy
‘variable vector and Z is the coenergy variable vector, the subscripts i
and (1 denote the independent and dependent parts of the storage field;
the output vector of the dissipation field D0 is the function of the
input vector to the field Di; Ti and T0 are the input and output vectors

of the modulated junction field.

13

Source Field

Se Sf

l I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Independent C fi—l I‘— TF
Storage I l‘—- Paynter Junction —=I Modulated
Dependent c in— k—— Juncnon
R R
Dissipative Field
Figure 2-2 Multiport Field Structure
—L..
__A i, U TV __
cbi 13- I Ii...
2 <1)
(pd ‘Ri‘ 19.. ‘
V Di A Do
(DI

 

 

 

Figure 2-3 Key vectors

14

When causality has been assigned to the bond graph according to the
Sequential Causality Assignment Procedure (SCAP) (Rosenberg and Karnopp,
1983), it is possible to identify each separate implicit R-field within
a graph. As a result of assigning causality, a computing diagram can be
generated, based on a set of key vectors for various fields in the
graph. Ignoring the Zd’ Xd vectors, the general equation structure is

given by

U - ¢u(t) (2.7)
2 - ¢X(X) (2.8 )
Do - ¢L(Di) (2.9)
dX/dt - $112 + 81300 + Sl4U (2.10)
Di - S312 + S33Do + S34U (2.11)

for systems with no 2 X vectors (dependent storages). If implicit R-

d’ d
fields exist in the bond graph, then the matrix S33 is non-zero. If S33
is zero, then no implicit R-fields exist, and a straightforward

procedure for integrating the system equations can be employed. If
implicit R-fields exist in the bond graph, then the matrix S33 would be
non-zero. These fields can be identified and isolated by proper use of
causality data. Referring to the field structuring of system equations,

the particular subset of equations with which we are concerned is

D0 - ¢L(Di) (2.9)

Di - S312 + S33Do + SBAU (2.11)
We seek an efficient solution to these equations at each time step,
given values for the Z and U vectors. A more succinct form for the

equations at a given time is

15

Do - ¢L(Di) (2.9)

Di - 333 D0 + C (2.12)
or

Di - S33¢L(Di) + C (2.13)

which lead to the single implicit vector equation, where C is a constant
vector. Clearly, it would be possible to iterate on D1 to obtain a
solution to the problem. There are two useful ways to reach the

solution.
(1). Single-partitioned iteration vector

If an R element with port j has its causality assigned uniquely,
then the corresponding row in the matrix S33 contains all zeros. The
value of D; is determined by D; — Cj. If an R element belongs to a
impliciizll-field therlits corresponding row in the matrix S33 contains
nonzero terms. The Di and DO vectors are reorgnized by gathering all
implicit input and output variables into the subvectors D: and D: ,
respectively, and all explicit input and output variables into the
subvectors D? and DE , such that D. a [D2 , D? ]t, and D - [DE , DI ]t.
1 o 1 1 1 o o o

The superscripts E and I denote the explicit and implicit field,

respectively. Rearranging the matrix 833, the upper part of the matrix

 

 

S33 15
DE 0 0 DE CE
1 o
- + (2.14)
I I I I
D1 0 833 D0 C

 

 

16

If the dimension of the DE and DE vectors is dE’ the dimension.of'
the D: and D: vectors is d1, then the dimension of the implicit
algebraic equation set to be solved is d1. That is,

I I

Do - Q (Di) (2.15s)
I I I I

Di - 833 D O + C (2.15b)

(3). Multiple-partitioned iteration vector

If the bond graph model contains several implicit R-fields and if
the vectors D: and D: are partitioned suitably, the submatrix S§3 is a

block diagonal matrix. Equation (2.15b) may be rewritten as

DIl ‘ S11 DI1 + C11
1 33 o
(2.16)
DIk = SIk DIk + CIk
1 33 o
where each subvector has dimension equal to mi, and m1 + m2 + ...... +

mk — dI' The Di and Do vectors are sorted into explicit and implicit
subsets; then the D: is partitioned according to its fields. This result
in the nonzero submatrix of S33 having a blcmfl<<diagonal form. Now our

task is to solve the following equation set for each local IRF:

D0 = ¢L(Di) (2.17)

D1 = S33 D0 + C (2.18)

17

Note that all vectors and matrices in Equations (2.17), (2.18) and.the

following are local representations.

The advantage of treating the implicit R-fields separately is; that
the total computation load is smaller. Since the computation cost of an
iterative solution to the IRFs appears to vary approximately as the
square of the size of the iteration vector, the savings by partitioning

into several IRFs can be substantial. Basically, we have a cost of

k k
a( 2 mi) rather than a( 2 mi)2 by partitioning effectively.
i=1 i=1

An example is given in Figure 2-4. Two implicit R-fields each
contains two l-port R.nodes. The causalities associated with the other
four l—port R nodes are determined uniquely by the SE , C, and I
elements. The dimension of the system vectors Di and D0 is 8. If the Di’
Do vectors are not partitioned at all, the computation time is 11.43 CPU
seconds for a simulation. If the explicit part is separated from the
implicit part, then the dimension of iteration vector D: becomes four,
and the solution time is reduced to 6.43 CPU seconds. If the D: vector
is partitioned into two local IRF vectors, then computation time is
reduced to 5.65 CPU seconds. This represents a 50.6 % savings in

simulation CPU time by partitioning.

18

l

.mmmH ocu mo xuwfimmsmo Any
.mmmH cowuwucopm wee Amy

mmmH 03o wewcwmucoo Emumxm < a-~

Apv
_ yn ya “u - ~_
A H H H a H
AlZTchIIZTloiIZTIo
a
_
Auv

ounwwm

19

2.3 IRFs with 1-port R Elements

An additional refinement for increasing solution efficiency has
been developed by Hood, Rosenberg, Withers and Zhou (1987). Tfluainethod
uses earlier work on the basis order for bond graph junction structures
(Rosenberg and Moultrie, 1980). In this section we impose two
restrictions on the R-field problem, which still leaves us with the most
common practical subclass of the general problem. (1) A11 R nodes in
implicit fields are l-port. (2) The bond graph does not contain
transformers (TF) and gyrators (GY) in the implicit fields. The
resulting problem is practically important, since many dynamic models of
engineering systems are contained in this subclass. Extension to more

general problems is discussed subsequently.

Now we ask these questions about the solution of a given IRF.
Subjecn: to the restriction imposed above, each obeys local equations of
the form of Equations (2.17) and (2.18). We ask.‘Wflnat is the smallest
number of iteration variables?" "How can such a setlxafbundT'"How
should Equations (2.17) and (2.18) be used for best computing

efficiency?"

To answer these questions one of the important properties of the
junction structure, basis order, will be used. The basis order is given
by two critical numbers, the number of effort inputs (E) and the number
of flow inputs (F) required at the junction structure (JS) ports in
order to determine all internal effort and flow variables and outputs.
The algorijflnn developed in this section will further reduce the number
of iteration variables of an IRF to the mininumicif its E and F. This

algorithm applies to IRFs that contain only l-port R elements.

20

The bonds of a given IRF can be sorted into one of three mutually
exclusive sets:

1. the boundary bonds, connecting the IRF to the rest of time graph
(these bonds are causal);

2. the R bonds, incident to R nodes, with which are associated the
local Di and Do vectors; and

3. the remaining bonds, which are internal to the local junction

structure (JS).

First we focus on the simple junction structiuna (SJS) which is
composed only of 0- and l-junction nodes. Earlier work by Rosenberg and
Moultrie (1980) has shown that there are two critical numbers associated
'with a JS. These represent the number of effort (E) and flow (F) inputs
required at the JS ports in order to deterndlmaeall internal variables
and the outputs. For completeness we state the rule here for calculating

the numbers:

- N - B (2.19a)

(2.19b)

where NB is the number of bonds of the JS:
NO is the number of O-junctions;
N is the number of l-junctions;

B0 is the number of bonds incident to the O-junctions;

B1 is the number of bonds incident to the l-junctions.

Next,vm observe that a SJS has a pair of separate but related
transformations associated with it. Namely. input effbrts determine

output efforts, and input flows determine output flows. Furthermore, if‘

21

all powers at the SJS ports are oriented outward (or all in), the

associated matrix 833 is skew-symmetric, subject to proper ordering of

the port variables (Ort and Martens, 1973, Perelson, 1975). Now assume

that causality assignment to the IRF has been completed and is

consistent. The proceeding observation allows us to organize Equations

(2.17) and (2.18) in detail as follows.

Sort the D1 and Do'vecbors into a resistance set (r) and a

conductance set (g). The r set has flow inputs to the R ports and effort

outputs whose dimensions are E; the g set has effort inputs to the R

ports and flow outputs whose dimensions are F. Write Equation (2 17) as

er - ¢r(fr) (2.20a)

fg - ¢g(eg) (2.20b)

where er and fr are associated with the rlxnuiset, and f

associated with the g bond set. Write Equation (2.18) as

f - S e + S f + C (2.21a)

e - S e + S f + C (2 21b)

are

Since the SJS transforms efforts to efforts and flows to flows, then Srr

and Sgg must be zero. Consequently, we have

f - S f + C (2.22a)

e - S e + C (2 22b)

Furthermore, we note that the combined set (ff, eg) contains the SJS

outputs,\flfile the combined set (er, fg) contains the SJS inputs. A

22

computational algorithm may be derived from this observation.

Algorithm

1. Assign causality to the source and storage nodes, using the

SCAP.

2. Identify each implicit R-field within the partially causal bond

graph.

3. For

a .

each implicit R-field:

Calculate E and F.

If either E or F is less than one, stop. (There is no
guarantee that there are unique outputs from the inputs for
the SJS )

Obtain a complete, consistent causal orientation for the
IRF. (It will obey the E, F numbers.)

Order the R bonds by resistance (r), then conductance (g)

causality. Define the vectors er. fr’ fg’ and e

. Assume that E is less than or equal to F. Use er as the

iteration vector. Make an initial guess eri for er

Use Equation (2.22b) to find eg

. Use Equation (2.20b) to find f8

Use Equation (2 22a) to find fr'

Use Equation (2.20a) to find er.

Compare eri to er. If the error is within tolerance, stop.
Else return to Equation (2.22b) and repeat the sequence
with an updated guess for er.
Note: If E is greater than F. use fg as the iteration

vector. The equation order is then (2.22a), (2 20a),

(2.22b), (2.20b).

23
It is observed that the minimum iteration set has the size of min
(E, F}. This is always less than or equal to one half the number of R
bonds. Reducing the dimension of the iteration vector has a significant

positive influence on computing efficiency, as noted previously.

The restrictions placed on the problem structure can be relaxed to
a certain extent without changing the algorithm as stated above. A given
IRF can contain R nodes with more than one port, provided each such R
node is a pure r, or a pure g, type. See Equation (2.20a) and (2.20b) in
this regard. In addition, the R-field junction structure can contain
transformers (TFs), since they do not alter the structure of the effort-
to-effort, flow—to-flow transformation properties. Their effects are
combined irnn: 8. See Equations (2.22a) and (2.22b). The formulas for E
and F need to be modified by introducing the term NT’ which is the
number of transformers in the junction structure. The modified formulas

are:

- B - N (2.23a)

F - N + N - N - B - N (2.23b)

An example of an implicit R-field is shown in Figure 2-5a. The
inputs to the IRF are ea, eb, and ec. The goal is to calculate all R-

field variables. We firstind E and F from the data given in Figure:

E = 6 + l - 3 - 3 = 1 (2.24a)

F = 6 + 3 - 1 - 6 = 2 (2.24b)

24
A solution does exist, since both E and F are greater than zero. We
obtain a complete, consistent causal orientation, as shown in Figure 2-
5b. The ordered r bond vectors are defined, based on the causality, as

indicated in the figure. The equations can be written as

e1 - ¢1(f1), (2.25)
f2 - 32(e2), (2.26s)
f3 - ¢3(e3), (2.26b)
f1 - f2 + £3, (2.27)
e2 - -el + (ea + eb), (2.28a)
e3 - -el + (ea + eC). (2.28b)

Since E is less than F, use er - e1 as the iteration variable. The

equation sequence is (2.28), (2.26), (2.27), (2.25).

a 4 6 c

It 51 3‘ External bonds: a, b, C
R-field bonds: 1 - 6
R19//'|\\<\ R3
R2 NB=6.NO-1.N1-3
BO = 3. B1 = 6
(a

.JL.41.1._40_2_.4IF.JL.
I ’1 l
R l .
ij/A‘ 2 R3 er . 1 l t fr l f1 1 t
R, f { fg. f j e = [ e e l

(D)

II
(D

II

Figure 2-5 An example of an IRF
(a) Bond graph of the R-field

(b) Causality assignment

25

2.4 IRFs with Multiport R Elements

2.4.1 Multiport Resistances in Physical systems

Multiport resistances arise in many kinds of physical systems, such
as hydraulic, thermal, and electronic. A well-known example of multiport
resistance can be found in a hydraulic system shown in Figure 2-6a
(Rosenberg and Karnopp, 1983). The four-way control valve is
characterized by four resistances that are formed by the four edges of
the spool lands and corresponding lands in the valve body and modulated
by spool position 2. The bond graph model is depicted in Figure 2-6b.
For a closed-center valve, one can work out the relationships between pm
and Q8 with 2 as a parameter. These relationships are nonlinear since
the pressure drops at individual ports (resistors).are proportional to
flow squared. The algebraic reduction of the R-field can be represented
in.a simplified‘way (Figure 2-7b) and the constitutive laws for the

modulated 2-port R are those shown in Figure 2-7a.

Another example of multiport R element may be found in an
electronic circuit. For a 3-terminal element such as the grounded-
emitter transistor of Figure 2-8a (Rosenberg and Karrmuna, 1983), if we

consider the collector-emitter voltage. e base-emitter voltage,

CE’ eBE’

and the currents i and i the 2-port R representation of Figure

C, 13! E!

2-8b with a mixed causality has the following constitutive laws:

1C = ¢1( eCE’ iB) (2.29)
i (2.30)

eBE ‘ ¢2< ecs, B)

The figure shows a possible connection for a power transistcn: in a

26

 

 

 

 

 

 

 

“XL—Eel

 

 

 

Figure 2-6 Four-way-valve amplifier

.15
+0.“ \

 

 

 

 

 

‘Qmu
-P

 

 

Ia)

Figure 2-7 Reduced form of hydraulic amplifier

27

switching mode. The control current i is switched rapidly between a

B
(low) blocking value and a (high) conducting value. During the blocking,

phase, ic is nearly zero even when e takes on.quite large values.

CE

During the conducting phase e is small even for fairly large currents,

CE
so that the voltage source is effectively applied directly to the load,

the averaged load current can be continuously varied.

 

 

 

 

Figure 2-8 Power transistor and its bond graph model

28

2.4.2 Extension of the Algorithm to an IRF with Multiport R elements

Assume a causal IRF with a number of l-port R nodes and M multiport
R nodes is given as shown in Figure 2-9, where R denotes a set of R
nodes in resistant form (R), C denotes a set of R nodes in conductance
form (C), while RGi denotes a multiport R node with mixed causaliigl. In
addition, U denotes the input variables, while V denotes the output

variables. S characterizes the junction structure. The key vectors are:

to
fr = [ fro’ frl’ ...... , er] .
e - [ e e e ]t'
g go, g1, 000000 ’ gM ’
e - [ e e e ]t'
r ro’ rl’ """" ’ rM ’
t
f = f , f , ...... , f ;
s [ so s1 all]
and in an aggregated form,
to
D1 - [ fr’ eg] ,
p - [e , f it
o r g
Ue vf Ve Uf

 

e e
ro , o
R h—‘f IRF JUNCTION STRUCTURE —E—>lf G

ro go

frllérl fgl-‘leg1 ..... erIFrM fgfilng

RCl RGM

 

 

 

Figure 2-9 Key vectors in an IRF

29
where fr is the flow input vector to the R ports; e is the effort input
‘vector to the G ports; er is the effort output vector from the R ports;
and f8 is the flow output vector from G ports. D1 is the input vector to
R nodes and Do is the output vector from R nodes. The vectors with
subscript o are associated with the l-port R nodes. The vectors with

subscript nlare associated with the m-th multiport R node, where m = 1,

2, ...... , M. The constitutive equations for the l-port R nodes are:
ero - ¢ro( fro) (2.31a)
f - o ( e ) (2.31b)

e = Q ( f , e ) (2.32a)

f - ¢ ( f , e ) m = l, 2, ...... , M (2.32b)

The input Di and the output D0 are linearly relatxxi through the
junction structure of the IRF. The relations between Di and D0 are

represented by the junction structure matrix S for the IRF, namely,

D. = S D + C (2.33)

1 o
where C is a constant vector derived from input vector U. Assuming that
the powers on all bonds of the R-field are directed towards the R nodes,

the local junction structure matrix S has the following form:

30

 

 

Srogo Srogl "' SrogM
S
0 rlgo 1
O Sf SngO OOOOOO rM I1
8 - ‘ Sgoro Sgorl "‘ SgorM
S 0 S
e glro
0
Sngo "' SngM

where S is skew-symmetric.

The connective Equation (2.33) may be written in some detail as:

f = S f + C (2.34a)

e = S e + C (2.34b)

fr0 a Lro(fgo’ fgl’ ...... , ng) + Cr0 (2 35a)

frm - er(fgo’ fgl’ ...... , ng) + Crm (2.35b)
and

ego = Lgo(ero’ erl’ """ ’ erM) + Cgo (4.36a)

egm - Lgm<ero’ erl’ """ ’ erM) + Cgm (2.36b)
where m = l, 2, ...... , M.

We now seek as an iteration set a minimum satisfactory set of
variablesm IIf we construct a functional diagram to represent the
mathematical relations among the variables basemlcna Equations (2.31),

(2.32), (2.35), and (2.36), it appears as in Figure 2-10. The weight of

31

each edge is the dimension of the associated vector. From this diagrann
the algebraic loops can be identified and a minimum iteration set may be
determined by removing a set of edges which break all loops and whose
sum of weights is a minimum. However, for our purpose of identifying
the minimum iteration variables, a simpler variable flow diagram in
Figure 2-11 can be generated such that a vertex represents a vector and
an edge represents a condition of input-output. This diagram enmflnasizes
the computation requirement of which input vectors are needed to
generate a particular output vector without any detail of their
functional relations. For convenience in our initial development let us

ignore the existence of e0, fO temporarily.

The algebraic loops can be identified from this variable flow
diagram. The iteration variable set with minimum dimension can.tu3 found
from its associated loops in the variable flow diagram. The removal of
the verthufi;(iteration variable set) from the variable flow diagram
destroys their associated loops. If the selection can be guided by some
rules, it will facilitate the automatic computation efficiently. For
consistency, we always select.the iteration variables from the input
variables to the junction structure, namely, er and f
Every set consisting of four vertices representing the vectors (er,
ff, eg, and f8)m and the edges among them is defined as a block of the
digraph. Each block connects (M-l) other blocks with the same structure
by 2*(M-l) outgoing edges departing from vertices fg and er, while the

vertices eg and fr connect the (M-1) blocks by 2*(M-l) incoming edges.

All blocks are arranged as shown in Figure 2-11 such tfluat all the

Vertices (\muztors) associated with r-type causality are located at the

 

 

 

 

Figure

 

2-10

32

 

 

 

 

 

W2
1 \ Li"
’ ID
{81 d)
fg2 0

 

 

Functional diagram of an IRF

34

right side, while all the vertices (vectors) associated with g-type
causality are located on the left side of the digraph. The vertex
associated with r-type causality is classified as r-type vertex” vflnile
the vertex associated with g—type causality is classified as g-type
vertex. In order to break all loops in the digraph, it is obvious that
one of the two vertices in each row ( either (e , f ) or (e , f ) )
g g m r r m
must be removed since it is necessary to destroy the loops formed by

these two-vertex sets.

Referring to this digraph with a repeated and symmetrical pattern,
every loop in this digraph contain both g-type and r-type vertices.
Therefore, removal of either all the r-type or all time g-type vertices
breaks all loops. Thus (eg, fg) or (er, fr) can be chosen as a feasible
solution. The size of the selected vertices (vector) set is two times
the dimension of (eg, f8) or (er, fr)' i.e., 2F'or 2E, respectively.

Considering the number of the basis order, E and F. the set of iteration

variables with smaller size can be chosen. Thus the size is given by

N = 2* min { E, F } (2.37)

For convenience in the next development, let us assume the case of
E S F, in which all the vertices on the right side are renmnmxi (If F S
E, then certain arguments are modified suitably ). Although E S F, this
does not guarantee that in every blockLHL lflr s Ng is true. Since the

erms air: the inputs of the junction structure, let us keep these in our

solution temporarily.

Under this condition for a particular block with er removed,

removing fg instead of fr also ensures no loops eumxxriated with this

35

‘particular block, because there are no more outgoing edges from this
block. A size reduction may achieved by removal of one of the two
vertices (fg’ fr)m with smaller size. By applying this step to each

block, N can be reduced from Equation (2.37) to

M

N - min { E, F } +mil min { Nr’ Ng }m (2.38)

Can we reduce N further? Suppose that there exists an i-th block
such thatzli . < N .. If vertex e . is removed instead of e ., at least
g1 r1 g1 r1

one loop will be found in the remaining digraph (Figure 2-12). The
reason is the following. Since N . < N_., the vertex f . has been
g1 11 g1

removed, while fri is kept in the last step. Howevwnr, since E 5 F, we

always can find at least one block, say the j-th block, such that Nrj s

N j. In this block, vertices erj and frj have been removed and vertices
e . and f . remain. It is easy to verify that a loop (e ., e ., f .,
81 8] r1 SJ 83
fri’ eri) exists in the remaining digraph. This shows that sucrl

exchange of vertices e . and e . is not allowed. even though N . S N ..
g1 r1 g1 r1

Thus no further reduction on N can be achieved.

Now consider the existence of l—port R nodes. These are represented

by the block consisting of the vertices e_ . f , e , anui f as shown
10 ro go ro

in Figure 2-133. It is a special case of the general blocks, namely, no

two-vertex loops exist. In the case of E <:l’. removing only er0 can

break all possible loops associated with this block. Therefore,
incorporating this advantage3hnx>the result. the smallest size of the
iteration variables for an IRT‘containing both l-port and multiportll

nodes with assigned causalities is given by

M
N = min 1 E, F l + 2 min i N
m=l

r' kgim (2.39)

36

 

 

 

Figure 2-12 Loop formed by an improper vertex exchange

38

Table 2-1 Possible variable selection

 

 

 

 

F _ E E _ F
M M
( U f ) ( U e
m=0 gm m=0 rm )
N S N + e + f
gm rm g1 g1
Ngm Z Nrm + er1 + fr1
M M
Note here that E - 2 Nrm and F = Z N m' An iteration variable set with
m-O m=O

this property can be identified.

From the analysis above, for each multiport R node four
possibilities exist, which are listed in Table 2-1U 1% simple rule to
determine the minimum iteration variable set for the IRF with multiport

R elements is derived from the variable flow diagram and Table 2—1 as

follows:
M M
1) If F s E, choose ( U f ). else ( U e );
m=0 gm m=0 rm
2) For each multiport R node:
M
a. In addition to ( U f ). if N, s N , choose e ,
m=0 gm 1m gm rm
else e
gm
M
b. In addition to ( U e- ), if N S N , choose f ,
m=0 1m rm gm rm

else f
rm

39

An example

4A simple electronic circuit is shown in Figure 2-143. After
assigning causality to the Sf and Se elements, an IRF is identified
(Figure 2-lhb). By assigning suitable causality to ports 3anu15, a

causally completed bond graph is obtained (Figure 2-14c).

 

 

R1 R2
11
' (b)

R: R2

1‘ 1”
sr l-J—*Rm‘,‘-3— 1t—L‘Rr—g— lF—b— Se
(C)

Figure 2-14 An electronic circuit and its bond graph model

40

The vectors associated with l-port R nodes are

e = [ e
ro

e

1' 2 ]

fro ' [ fl’ f2]

and

e = f - [0]
8° go

The vectors associated with multiport R nodes are

e

r1 3 fr1 = [O]

fgl - [ f3 1

eg1 = [ e3 ]

and

8r2 - [ ea ]
f - [ f

r2
£82 = [ f5 ]
eg2 - [ e5 ]

we may calculate the numbers E and F by

2

E - 2 dim (erm) - 2 + O + l = 3 (2.40a)

m=0

2

F - 2 dim (fgm) = O + l + l = 2 (2.40b)

m=0

The minimum number of sufficient iteration variables is determined by

2
II

M

min { E, F ) + 2 min { Nr . N }

m=l

min { 3, 2 } + min { O. 1 } + min < 1. l } = 3 (2.41)

41

The iteration variables are chosen by the rule as follows: :f3, f5,

and e5. The constitutive equations for the R nodes and the junction

structure are

Do ero: el - ¢l(f1) (2 42a)
e2 - ¢2(f2) (2.42b)
fgo: -----
erlz ....
fgl: f3 = ¢3(fa, e3) (2.43)
er2: e4 = ¢4<f4’ e5) (2.44)
fg2: f5 - ¢5(f4, e5) (2.45)
Di fro: f1 = f3 (2.46a)
f2 - f5 (2 46b)
ego' -----
frlz -----
egl' e3 = - el - eA (2 A7)
fr2 f4 = f3 (2.48)
eg2. e5 = eb - e2 (2.49)

Time iteration variables are fg0 U fgl U ng U eg2' namely, (f3, f5, and

e The iteration process may start with a set of initial guess (£31)

5)'

f5i’ and e51) to (f3, f5, 5).

(2.48), (2.42). (2.44), (2.47). (2.49), (2.45). (2.43). A new set of f

e The iteration sequence order is (2.46),

3 )
f5, and e5 are generated and become updated initial values. Thus the

iteration can be repeatedtnnjl the errors are within specified

tolerance.

42

2.4.3 Assigning causality to an IRF Containing one Multiport R node

From the study of iteration variables for a causal IRF containing a
multiport R node, it is clear that if the multiport R node has its
causality assigned in such a way that it is in either complete R-form or
complete G-form ( i.e., either Nrm - O or Ngm - 0), the number of
iteration variables is minimum and equal to the smaller of E and F. If
such a causal assignment is not possible, we should make the additional
iteration variable set as small as possible. That is, either Nrm or N

gm
is the minimum. How can we determine the minimum number of Nrm or Ngm?
How do we assign causality so that the minimum number of iteration

variables can be used? This subsection addresses these questions.

., R are

Figure 2-15 is a diagram for an acausal IRF where R1. n

l-port R nodesanui the large R is a multiport R node. Ue and Uf are

input vectors to the IRF and Vf and Ve are output vectors.

Ue‘L Vf Ve-[Uf

Junction Structure

R 3 §

 

 

 

 

Figure 2~15 Acausal IRF with multiport R node

43

Suppose a multiport R node has n ports. If we start at one of the
ports to trace out a path consisting of alternating bonds and junction
nodes (0, l, TF), the path may end up at a loport R node or the
multiport R node itself. A path starting at and ending up at the same
multiport R node is called a self-loop. Several paths may share some
common ports and common junction nodes. The collection of the ports of
such coupled self-loops and all the junction nodes in the paths, and all
l—port R nodes adjacent to the junction nodes. is defined as an R-block.
Figure 2-16 shows an implicit R-field containing a multiport R node
which has a causally determined bond 7 and three acausal bonds 4, 5, and
6. A self-loop is identified, that is (Rm - 0a - la - Ob - lb - Rm). The

R-block associated with this self-loop consists of node Rm, 0a, R la,

2 l

0b, 1b, and R3.

3 I7 R
1 I 4 5 y 3
Rf“— O —“Rm"— 1b

Figure 2-16 An example of self-loop and R-block

44

Before we make a further discussion, we need to have the following
assumptions. First, no redundant R ports exist in the IRF, namely, no
two or more l-port R nodes are allowed to be adjacent to a common
junction node, since such l-port nodes may be aggregated ixnua a single
lrport R node. Second, no short loop exists among the multiple R ports,
namely, no two or more ports of a multiport R node are adjacent to a
common junction node. This has the same reason as for the first

assumption.

It is easy to show that a port of the multiport R node, which is
not contained in any R-blocks, can be assigned with either R-form or G-
form causality. For example, bond 4 of the multiport R in Figure 2-16 is
not in any R—block, it can be assigned with either effort in or effort
out. However, the ports in an R—block, for example, bond 5 and 6, may or
may not be assigned freely depending on the local junction structure of

the R-block.

Each R-block has its own local basis order pair. Let Mi denote the
number of ports associated with the multiport R node hithe idfllR-
block. Let E1 and F1 denote the basis order pair for the i-th R-blockg
Cixmni E and F for a SJS or a WJS. not any subset of ports can always be
assigned in, say, R causality. However. if we assume that any subset of
ports in an R-block with given E and F can be assigned desired causality
obeying E and F, we can derive a possible minimum number of iteration
variables and guide the causality assignment from the discussion on four

types of R-blocks.

45

Type 1. M1 is less than or equal to both E1 and Fi'

Since no short loop is allowed in the bond graph, tracing the
causal orientation from a port of the multiport R node will always end
up at a l-port R node. The M1 ports of the multiport R node can be
assigned either all in R-form or G-form. Figure 2-17 illustrates an R-
block of type 1. The local basis order is calculated to be E - 3 and F -
3.‘The number of acausal multiport bonds in the R-block is 2. These two

ports can be assigned either in R-form or G-form.

Type 2. M1 is greater than Ei but less than or equal to Fi'

In this case, Mi ports of the multiport R node can be assigned all
in G-form, but not in R-form. In Figure 2-18. since E1 - l, and F1 - 2,
the twatmflmdports can not be assigned as all R-form. but it may be

assigned as all G-form.

R3 R3
_1/ Rl--1\/
Rm b m b”\
\ ‘l'
/ 4 I \/
0 ""‘ 1 0 la
a a T /\
l \R R.
R2 5 R2
(a) (b)
NB=9. VO=2. 311“”
BO=6, 81:6 211:2
E1=9+2-2-6=3 E1>M1
r1=9+2 2 TH) r1>~1

Figure 2-17 An R-block of type 1

46

 

0,,"""'1a
R
(a)
NB = 6, N0 = 2, N1
BO - 5, B1 - 4 1
E1 - 6 + 2 - 2 - S -
F1-6+2-2-4-
Figure 2-18 An R-block of type 2

l ///,0b
0, — 1,\
R
(a)
NB = 6, N0 = 2 N1
BO = 4, B1 = 3 11
El = 6 + 2 2 - 4 =
Fl = 6 + 2 2 3 =
Figure 2-19 An R-block

of

H99

(b)

H

47

Type 3, M1 is greater than Fi but less than or equal to E1.
The result is opposite to that in case 2. The R-form can be
assigned on all Mi ports, or F1 ports can be in G-form and the remaining

ports must be in R-form. The causality assignment can be seen from the

example in Figure 2—19.

Type 4. M1 is greater than both E and Fi‘

1

Obviously, the M1 cannot be assigned either complete R-form or G-
form. They can be assigned Ei R-form causalities and (Mi-Ei) G-form
causalities, or inversely. Fi G-form and (Mi-F1) R-form causalities. It

is shown in Figure 2-20.

 

(a) (b)

3 = a 9
NB 8. NO 2, N1 _

= = .N = ‘
BO 6. Bl 5. J1 I
‘ = 3 + 3 ._ 4’ ‘ - L: "x :'
L1 D 1 11
F1 = 8 + 2 2 3 = 3 F.1 \ V,

Figure 2-20 An R-block of type 4

48

Let NO denote the number of free ports of the multiport R node,
which are not in any R-blocks. An IRF is decomposed and all R-blocks are
identified and classified according to above types. Let Ei’ F1 and Mi
denote the basis order pair and the number of ports of the multiport R
node in the i-th type-l R-block, respectively, where i - l, ..., I. Let
Ej, Fj and M. denote the basis order pair and the number of ports of the
nmltiport node in the;i¢flxtype-2 R-block. respectively, where j a l,

., J. Let E Fk ang M denote the basis order pair and the number of

k!
ports of the multiport node in the k-th type-3 R-blcmflc, respectively,

where k - 1W ..., K. Let E F1 and M1 denote the basis order pair and

1 v
the number of ports of the multiport node in the Imwfli type-4 R-blockg
respectively, where l - l, ..., L. For a particular multiport R node, we

may try to assign asrmnu'R-form causalities as the structure allows.

This process is called the maximum R-form assignment. The maximum G-form

assignment is just its opposite.

In the maximumiflmform assignment, the maximum possible number of
effort outputs and the minimum possible number of flow outputs from the

multiport R can be determined by

I J K L

Er = 2 N1 + E Ej + Z Mk + E E1 + NO (2.50a)
J L

Fr = 2 (Mj- Ej) + 2 (M1- E1) (2.50b)

where Er is the maximum possible number of ports with R causality and Fr
is the minimum possible number of ports C causality for the multiport R

node.

49

In the maximum G-form assignment. the minimum possible number of
effort outputs and maximum possible flow outputs from the multiport R

node can be determined by

K L

Eg - Z (Mk- Fk) + 2 (MI- F1) (2.51a)
I J K L

Fg = 2 Mi + 2 Mj + Z Fk + 2 F1 + NO (2 51b)

where Eg is the minimum possible number of ports with R causality and F
is the maximum possible number of ports with G causality for the

multiport R node.

The smallest dimension set is found from Fr and Eg' If Fr 5 Eg , then
the maximum R-form assignment would be used and the number of iteration
variables of IRF would be the minimum. Since Eg is the minimum number
of the R-form bonds, it is the smallest possible size of the vector erl’

i.e. DJ in Rule 1; since F} is the minimum possible number of the G-

r1

form bonds, it is the smallest possible size of f i.e. N in Rule 1.

gl’ g1
By comparing E8 and Fr’ we find smallest additional iteration set to the

set determined by min {E, F}. In general. The minimum possible number of

iteration variables is determined by:

N = min (B, F) + min {E , F,} (2.52)
g 1 m
This result may be applied to guide finacmumality assignment for
the multiport R node in computation automation. Let us discuss the

implications of the above result.

50

1. If there are no R-blocks, Fr and E8 are equal to zero. This fact
suggests that the causalities of the multiport R node can freely be

assigned in either complete R-form or G-form.

2. If only type-l R-blocks exist, Fr and E8 are also equal to zero.

(Same as in case 1.)

3. If neither a type-3 nor a type-4 R-block exists, but type-1
and/or type-2 R-blocks exist, Eg = O and the multiport R node can have

complete G-form causalities.

4. If neither a type-2 nor a type-4 R-block exists, but type-l
and/or tqu3-3 R-blocks exist. then Fr = O and the multiport R node can

have complete R-form causalities.

5. If a type 4 R-block exists, then Fr ¢()andiig¢ O.The

multiport R node will not be in complete R-form or G-form for sure.

From the discussion above, if there is no self-loop of the
multiport R node, the multiport R node is not structurally coupled
(Figure 2-21). In this case. from the resulteflxnmn it is easy to see
that thermfltiport R can always be put into either R-form or G-form
completely. The causality can be assigned according to the required
physical laws of the multiport R node. The minimum number of sufficient

iteration variables can be determined by the previous rule, i.e.

N = min i F, F } (2.53)

51

Figure 2-21 A multiport R node without structural coupling

For example, if the multiport R node is in complete R-form, the

constitutive and system equations will be reduced to a simpler form

er0 - ¢ro( fro) (2 54a)

er1 - ¢r1( frl) (2.54b)
fgo = ¢go( ego) (2.55)
f = S f (2.56a)
ro rogo go

fr1 = Srlgo fgo (2.56b)

e (2.57)

= S e, + S e
go goro io rlgo r1

If the multiport node is associated with self-loops. then this
node is structurally coupled. The multiport R node can or cannot be put
in either R-form or G-form completely, depending on the structure of the
IRF. If the multiport R node can not be put into either R-form or G-form
(i.e.,ndxed fann), the coupling among the ports of the multiport R
node is called essential structural coupling. If the causality is
assigned. we will find a causal loop which smmarts at one of the

multiport node and ends up at another port or the same multiport node.

52

A strategy for assigning causality to an IRF containing one

multiport R node is stated as follows:

1. Assign causality to the source and storage nodes, using the

SCAP.

2. Identify each implicit R-field within the partially causal bond

graph.

3. For each implicit R-field:

a .

Calculate E and F.

If either E or F is less than one. stop. (There is no
guarantee that there are unique outputs from the inputs for
the IRF.)

If a multiport R node exists. identify all self-loops and

associated R-blocks. For each R-block:

Calculate the local basis order (E. F)m and Mm'

Classify the type for the R-block.

Calculate the pair of (E . F ) and determine the possible

g r m
minimum number of iteration variables by
M
N = min {E, F) + 2 min {B . F_}
i=m g i m

Assign causality to multiport R node first according to

min (Eg’ Fr)m. Obtain a complete. consistent causal
orientation for the IRF. (It will obev the E, F numbers and
<5. mm.)

Order the R bonds by resistance (r). then conductance (g)

causality. Define the vectors er, fr’ f , and e

. Assume that E is less than or equal to F and Eg is less

than or equal to Fr‘ then R = E + Eg' Use e , e

and f
ro r

rl’ l

53

as the iteration vectors. Make an initial guess e ., e

r01 rli’

, e

and f . for e
r11 ro

rl’ and frl'

. Use Equation (2.36) to find ego and egl'

. Use Equation (2.3lb) to find fgo'

. Use Equation (2.32a and b) to find er1 and fgl'
. Use Equation (2.35) to find f and f .
r0 r1
. Use Equation (2.31a) to find ero'
Compare eroi’ erli’ and frli to ero’ erl’ and frl' If the

error is within tolerance, stop.
Else return to Equation (2 36) and repeat sequence

with the updated guess for e . e and f .
ro rl

rl’

Note:

1. If E is less than F and Fr is less trmn1 Eg' use e and

e
ro’ rl’

fg1 as the iteration vectors. The equation order is then (2.36),

(2 31b), (2.35), (2.31a), (2.32).

2. If F is less tfluni E and E8 is less than Fr‘ use fgo’ f and
er1 as the iteration vectors. The equation order is then (2.35),

(2.313), (2.32), (2 36), (2.3lb).

3. If F is less than E and F is less than E . use i? , f , and
r s so 81

e81 as the iteration vectors. The equation order is then (2.35),

(2.31a), (2.32), (2.36).(2.31 b).

The following two examples show how this strategy works. The first
example in Figure 2-22a depicts an isolated IRF with one multiportll

node. The pair of basis orders are calculated first:

54
E-8+2-3-5-2 (2.58a)

F-8+3-2-6-3 (2.58b)

At this stage we do not know how to assign causality to this IRF
and also we are not sure that the number E - 2 is the number of the
iteration variables. Now, a self-100p (Rm - 0a - la - Ob - lb - Rm) is
identified, and so the associated R-block (Rm, Oa, la, Ob, 1b, R5). The

local basis orders of the R-block is obtained by:

E1-6+2-2°5=l (2.59a)

F1-6+2-2-4-2 (2.59b)

3
RH 1 ‘4le le

c
4_.{ Sf
0,, 6 \Z/
6 o I‘— 1
0a 18 d a a’\
is [\Se 13L? SC
R2 2
2(a) _ :b‘)

Figure 2-22 Example I

55

The number of ports of the multiport R node in this R-block is 2.
Since PH .. 2 > E1 = l, and M1 - F1, this R-block belongs to type 2. We

have F - M - E -=1” and E — M - F = O. It means that this
r l l g l l

multiport R.node may be in complete C-form. but not in complete R-form.

The possible minimum number of the iteration variables is determined to

be:
N - min {2, 3} + min { 0, l} = 2 (2.60)

Based on this result, the multiport R node is assigned to have
complete G-form causality (Figure 2-22b). The associated input and
output vectors are identified and the constitutive and connective

equations are listed as bellow:

er0 e1 = 61(f1) (2.6la)
e5 = ¢5(f5) (2 61b)

erl ---

f ---

go

fgl f2 = ¢2(e2, e3. ea) (2.62a)
f3 = ¢3(e2, e3, ea) (2 62b)
f4 = ¢4(e2, e3, ea) (2.62c)

fr0 f1 = f2 (2 63a)
f5 = fc - f3 - f4 (2.63b)

fr1 --

e ---

go

egl e2 = ea - el (2.648)
e3 = eb + ed - e5 (2.64b)
e4 = e5 (2 64c)

56

. . . . t . .
‘The iteration variable is are. [ e1. e5] , and the iteration

sequence order is (2.64), (2.62), (2.63), (2.61).

The second example shows a case in which a multiport R node is not
able to be put in a complete R-form or complete G-form. T he isolated

IRF is given in Figure 2-23a. The pair of basis orders are:

E - 11 + 4 - 3 - 9 - 3 (2.65a)

F - 11 + 3 - 4 - 7 - 3 (2.65b)

Two self-loops are identified, those are (Rm - (M1 - la - 0b - lb
Rm) and (Rm - lb - Oc - lc - Rm). The associated R-block is (Rm, Oa, la,

Ob, lb, Oc, 1c, R6). The local basis orders are:

E1 = 9 + 3 - 3 - 7 = 2 (2.66a)
F, = 9 + 3 - 3 - 7 - 2 (2 66b)
Sc Se
1
3f 11:1— oc/C‘Sf sr 1,347— 0g}Sf
1. a]. .2 i3 5 is 1 a1[ 2 it; 5 8
lifh"' ""'th‘h" lb-‘<: ‘b Iii‘__' (ll—‘7'th‘7"'lb 9 l3
4 (h;"-‘Sf 4 (Efr'4.Sf
0.41— ,a/m 08..., ,3 .0
6 i A§<ESe: 6.i AQNE Se
R6 R6

(a) (b)

Figure 2-23 Example 2

57

Since M1 -:3 in this case, we have M1 > E1 and M1 > Fl' This R—block is

of type 2 and Er - 3 - 2 — l, Fgl- 3 ~ 2== 1. It tells that the
multiport R node can not be assigned in either a complete R-form or G-
form. The possible minimum number of iteration variables is determined
by:

N - min {3, 3) + min {1, l} = 4 (2.67)
Assigning causality according to a maximum R-form scheme, the causal
graph is shown in Figure 2-23b. Note here that we may also use the

maximum G-form scheme in this problem. The constitutive equations are

given and the connective equations are derived from the causal graph:

Do: er0 ----
er1 e2 - ¢2(f2, f4, f5, e3) (2.68a)
eA - ¢4(f2, f4, f5, e3) (2.68b)
e5 - ¢5(f2, f4, f5, e3) (2.68c)
fgo f1 - ¢(e1) (2.69a)
f6 - ¢(e6) (2.69b)
fg1 f3 = ¢3(f2, f4, f5, e3) (2.70)
Di fro ----
fr1 f2 = fd - f1 (2.7la)
f4 = fb - (fC - f3) - f6 (2.71b)
f5 = 1?C - f3 (2.71C)
ego e1 = e2 (2.72a)
86 = e4 (2.72b)
eg1 e3 - ee + eS - (ea -éa) (2.73)

The iteration variables are f3, f1, f6’ and e3anuithe iteration
sequence order is (2.71), (2.68). (2.72). (2 73), (2.70), (2. 69)..'The
result can be extended to.aIHU?containing several multiport R nodes if

some modifications are made.

58

2.5 IRFs with General Junction Structures

2.5.1 General Junction Structures

The jruuction structure of a bond graph is that portion which
represents the power-conserving features of the modeled system; as such
it can be viewed as a multiport transformation. In an abstract sense,
the junction structure represents the energy topology of a bond graph in
the same way that a generalized linear graph represents the topology of
an electrical network. Many researchers have made investigations into
the properties of junction structures (Karnopp, 1969; Rosenberg, 1971,
1978, 1979; Ort and Martens, 1973; Perelson. 1975). Many system
properties, such as solvability and basis order, are obtained by
studying the junction structures (Rosenberg and Andry, 1979; Roserfiuarg,

1980).

Junction structures are classified as simple (SJS), weighted (WJS),
or general (GJS), according,to whether their junction elements are in
the set (0 and l), (O, l, and‘TF)<nr(O. 1. TF, and CY), respectively.
The properties of SJS are contained in WJS. so we shall not discuss SJS

further.

Thezlnasis order of a WJS is the number of effort variables (E) and
flow variables (F) thatrmun:be specified for a particular WJS to have
all its power variables known. The pair of HZ.EU is a unique property
of a WJS. The principal application of the basis order is to determine,
in advance of assigning causality, the number of independent port effort

and flow variables for a given junction structure. The basis order

59

formula for a given SJS or WJS has been derived by Rosenberg and

Moultrie (1980).

The existence of gyrators in a junction structure generally makes
the basis order pair no longer unique. The transformation between input
and output vectors is not in the form of effort—to-effort and flow-to-
flow, since gyrators transform effort to flow (or flow to effort)
simultaneously on their ports. The results obtained.in.the previous
sections may not be applied to IRFs containing gyrators. Existing
luumvledge about the properties of the general junction structure is not
complete. In order to extend the results for minimum iteration variables
to an IRF with gyrators, we must develop some additional properties of
the general junction structure. Until now there has been no general
theory that treats the CJS systematically. In this section we develop

the basis properties of the CJS.

jFirst we introduce the gyrograph, derived from a GJS in an
abstracted form. Then we focus our attention on the relation between
causality in bond graph and the graph matching concept in the gyrograph.
Following this, we present an algorithm for determining,tflne existence,
the Itqninimum, and the F-maximum of the basis order for a GJS. Finally,

some examples are given.

2.5.2 Gyrographs and Maximum Matching

2.5.2.1 Gyrographs

A gyrobondgraph (GBG) is a bond graph formed from a limited set of

multiport nodes, called the primitive set- One such set i4; {1, CY, Se,

60
I, R}. See Appendix A for details. The junction structure of any bond
graph can be represented in CBC terms by the use of l junctions and CY
nodes. Although the CBC represents a canonical graph form, its essence
is the same as that of the standard (9 node type) bond graph. The
smaller primitive element set may facilitate the study of some
properties of bond graphs. For our purposes, however, it does not seem
to be a suitable graph representation for applying the well-developed
linear graph theory and associated algorithms to explore the properties
of general junction structures. Therefore we consider a more abstract
form of the CBC, called a gyrograph (CC) (Paynter, 1967; Zeid, 1982).
The CC preserves all the information of the CBC and therefore all the
information of the original bond graph. However, it depicts the
topological information and connectedness more clearly for certain

purposes. Now let us briefly define the gyrograph.

Nodes of a CC are represented by squares and circles. Edges of a CC
are lines. The square nodes are used to mark multiport field elements,
such as Se, 1, and R. These nodes are called environment nodes. The
circles represent l-junctions. An edge represents a gyrobond connector
(i.e., 2 bonds and a CY) between a pair of l-junctions if it joins two
circles. An edge represents a bond if it joins a square and a circle.
The former is called a junction edge. The latter is called an
environment edge. The junction structure is represented by the circles
and the edges; its ports are the environment edges. See Appendix A for

details.

In a CC the set of environment nodes (squares) is denoted by VEn'

The set of junction nodes (circles) is denoted by V The set of

J .

environment edges is denoted by E The set of junction edges is

En'

61

denoted by EJ. The nodes of the CC are denoted by V, where V - VEn + VJ.

The edges of the CC are denoted by E, where E = EEn + E]. No edges are
allowed between the environment nodes. since there are no ports are

defined between elements Se, Sf, C, I, R.

The example in Figure 2-24 illustrates the derivation of a
gyrograph from a standard bond graph model. Consider the R-field of
Figure 2-24a, which has been isolated from a system. Parts b arui<: show
the associated gyrobondgraph and gyrograph, respectively. The nodes and
edges are labeled to show the corresponding elements as the

transformations are made from the SBC to the CBC to the CG.

2.5.2.2 Causality and gyrographs

Standard bond graphs and gyrobondgraphs have the same causality
properties and can be assigned by the sequentia1.cnmasality assignment
procedure. For the primitive set ( I. R. Se. 1. and CY ) we have the

following causal properties:

Se has the effort directed outward:

I has two types of causality, integral and derivative;

R has two types of causality, resistance and conductance;

CY has two permissible causal forms. namely, both efforts directed
inward or both efforts directed outward; and the

l-junction has exactly one flow directed inward.

The complete causality assigned to a gyrobondgraph may be
recognized in the associated gyrograph by the marking on the edges. lfiar

edges in E a marked edge indicates a causal orientation such that both

J)

62

 

 

A B c D E F

1’ 01 GY‘—14-—0‘—GY4 0 >1

R1 R2 R3 R4 R5 R6
(60

A 13,0 D E F

lacYg—lthé—lA—GYL 14 CY4——l

J {\ J J J

 

R1 G] R3 G] G] R6
1 C l H l I
R2 R4 R5

(b)

 

 

 

 

(C)

Figure 2-24 Graphs of an R-field.
(a) The standard bond graph.
(b) The gyrobondgraph.

(c) The gyrograph.

63

efforts are directed as inputs to the gyrator (or equivalently, both
flows are directed as inputs to their adjacent l-junctions). An unmarked
EJ edge of a causal CC indicates the opposite causal state. For edges in
E a marked edge indicates that the flow'is input to the circle (1-

E’
junction), and the effort is input to the square (environment node). An.
unmarked EE edge of a causal CC indicates the opposite state.

Figure 2-25 shows causality and its equivalent in three related
graphs: the SBC oi'Figure 2-24, CBC, and its CC. We start with a
complete causality for the SBC in part (a). Part (b) shows the
corresponding CBC causality, and part (c) shows the equivalent causal

information in the causal (i.e., marked) CC. Observe that each circle

node in the CC has exactly one marked incident edge.
2.5.2.3 Causal 665 and the cardinality matching problem

From the pattern of the marked edges in a gyrograph, we shall refer
to an edge-marked CC as a causal CC, since it carllna interpreted as a
causal CBC. A causal CC also leads itself to interpretation in terms of
the cardinality matching problem of standard graph theory. Briefly, the
cardinality matching problem (Syslo, et al. 1983) may be stated as: In a
given graph, find a maximum matching, that is, a matching with as many

edges as possible. The concept of matching is defined next below.

Consider a graph G defined by a set of vertices V and a set of
edges E. By definition, a set of edges M in an undirected graph G = (V,
E) is called a matching if no two edges in M have a node in common. For
example, in Figure 2-26 the edge set {(v1. v6), (v2, v4)L,shown ha

heavy lines, is a matching, because the two edges are not incident to a

64

l—AGYH-IP—OHGYlé—OH—‘l

II II

(a)

yup—HIP
zlAL——4c>m

A B,C D 1»: F
111—— CY‘—-{1|l——CY‘—-{1|‘—GY‘-—il‘—-I can—>1

I 1\ I I I

R1 R GY R6

I. I. I.
1 I

 

Figure 2-25 Causality in different graphs
(a) in a standard bond graph.
(b) in a gyrobondgraph.

(c) in a gyrograph.

65

common node. The set {(v1, v6), (v2, v5), (v3, v4” is another matching.
But {(v2, v3), (v3, v4), is not a matching because the edges share a
node, namely v3. A single edge in a graph without self-loops is

obviously a matching.

Clearly, a maximum matching in a graph with n nodes can not have
more than (n/2) edges. It may have fewer edges. With respect to a given
matching M in a graph C, an edge is said to be matched if it is in M.
Similarly, a node x is said to be matched or saturated if it is an end
node of some matched edge, say, (x, y). A node that is not matched is
called an exposed or free node. In Figure 2-26 nodes v3 and v5 are
exposed and the rest are saturated. A path P - (v1, v2, ...... , vk) is
called an alternating path with respect to a given matching M if the
edges of P are alternately in M and not in M. An alternating path that
begins at an exposed node and ends at another exposed node is called an
augmenting path. The maximum matching algorithm is based on Berge's
Theorem that a matching M in a graph G is maximum if C has no augmenting

path with respect to M (Syslo, et al., 1983).

For our later purpose, we introduce the concept of an alternating
cycle. An alternating cycle is an alternating path that begins at an
exposed node and ends at the same node. For example, in Figure 2-26 the
path (v3, v2, v4, v3) is an alternating cycle with respect to matching M

- {(v1, v6), (v2, v4)). So are the paths (v5, v4, v2, v5) and (v5, v4,

v2, v1, v6, v5).

For a l—junction node, exactly one incident bond must be set with a
flow input to the 1-junction. The implication of this requirement in a

gyrograph is that a circle node must have. exactly one marked incident

66

 

 

Figure 2-26 A graph and one of its matchings

I

edge. In a gyrograph marked with the equivalence of complete causality,
all the circle nodes must have exactly one marked incident edge. Hence
they are all saturated. This property will be utilized to explore the

basis order of a CJS.

2.5.3 Basis Order Algorithm for General Junction Structures

From Equations 2.23a and b we see that the basis variable structure
(E and F pair) for a given WJS is unique. However, for a junction
structure with essential gyrators (Rosenberg, 1979: Breedveld, 1984) the
basis variable structure is not unique. That is. the number of effort
variables (E) and the number of flow variables (F) required as inputs to

determine the general junction structure is not a unique pair.

We here present an algorithm for determining the existence of basis

order properties for an arbitrary CJS. If a basis does exisn:, then we

67

can determine the minimum and the maximum flow basis orders, Fmin and

Fmax’ respectively. We will employ the gyrograph and apply a variation
of the standard cardinality matching algorithm to find the properties of

the basis order.

In the study of the junction structure, we will sometimes focus our
attention to the portions of the gyrograph which contains only jLunction
nodes and edges between them; namely, subgraphs induced by junction node

sets. The following definitions are introduced to assist in development.

Definition 1. External junction node: A junction node which is
adjacent to at least one environment node is called an external junction
node. The set of external junction nodes is denoted as VE'

Deflhfixion 2. Internal junction node: A junction node which is
adjacent to no environment nodes is called an internal junction node.
The set of internal junction nodes is denoted as VI'

Definition 3. Junction gyrograph: A subgraph of a gyrograph in
which alleundionment nodes and their incident edges are removed is

called a.jrnn3tion gyrograph. We use the symbol C to denote a junction

J
gyrograph.

For a given gyrograph to have a basis each junction node must
exactly have one marked incident edge. Such nodes are said to be
saturated. This corresponds to each of the l-junctions in the associated
CBG having an incident bond which gives the junction a strong causal
determination. The basis order properties are derived from the junction

gyrograph CJ and the gyrograph CG.

68

The nodes in a junction gyrograph GJ may be classified as belonging

to either the external set VE’ or the internal set VI. In a marked

gyrograph the nodes may be saturated or free. The collection of free

nodes in the internal set and external set are denoted as Vi and v: ,

respectively. The collection of saturated nodes in the internal set and
the external set are denoted as V: and VE, respectively.

Let us study the junction gyrograph first. Suppose we have a
maximum matching for a given junction GC. There are four cases (tuit all

mutually exclusive) that can exist. We present each case and discuss its

implications.

Case 1. All the nodes of internal set VI are saturated.7flun1a.basis

solution exists.

Case 2. There exists an alternating path P which starts at a free node

in the internal set VI and ends at a saturated node in the external set

VE' Then we reverse the marking of the edges ( i.e., the edges of P

currently in matching M are removed from the matching and those edges of
P not in M are put into the matching). The size of the matching M is
unchanged and so is the number of saturated nodes. However, the starting
node ir1\L[ is now saturated, while the ending node in VB becomes free.
By repeating this procedure we can move all such free nodes fhxnn Vi to
VE. An example of this case is shown in Figure 2-27a. It depicts a
portion of the maximum matching generated by a maximum matching
algorithm. The matched edges are shown in heavy lines. Along the
alternating path (v, x, y, z, u), node v is a free internal node,

whereas u is a saturated external node. This path can be changed to that

of Figure 2-27b by reversing the markings.of the edges. Thus the node v

69

(a) (b)

Figure 2-27 An alternating path.
(a) Before reversing markings.

(b) After reversing markings.

V V

2 y 2 Y

(a) (b\

Figure 2-28 An alternating cycle.
(a) Before reversing markings.

(b) After reversing markings.

70

becomes a saturated node, while the node u becomes a free node.

Case 3. There exists an alternating cycle C, which starts and ends at a

. . f .
free node v in the internal set VI and contains at least one external

node u in VE' Then we can reverse the markings of some of the edges of C
to make v saturated and u free. This procedure does not change the size
of the matching set M or the number of free nodes, but it moves the
position of the free node from the internal set VI to the external set
VE' This is illustrated in an example (Figure 2-28). Suppose tflnat after

a maxinnun matching is generated, an alternating cycle C (x, y, z, u, v,

x) is identified. There is a free node x in V The cycle contains nodes

1'
z and u which belong to the external set VE. If we choose node u to be
the destination for the free node x, we reverse the markings starting at
the marked edge (u, v) along the partial cycle (u, v, x) until edge (v,
x), such that node x becomes saturated and node u becomes free (Figure
2-28b). We may also choose node 2 to be the destination of the free node
x, by reversing another partial cycle (2. y. x) starting at the marked
edge (2, 3y). From this example, we can extract the fact that any

saturated node of V in the alternating cycle can be the destination of

E

an internal free node in the same cycle.

Case 4. There exist one or more free nodes which may not be saturated by
means of alternating paths and/or cycles. This means that the l-junction
associated with the free node can not be assigned by deterministic
causality. Hence no solution exists for this CJS.

The full saturation of the set V implies the completeness of the

I

deterministic causality assignment to the internal junction nodes, since

every node has one flow input. The free nodes in the set VE may get

71

saturated by their adjacent environment nodes. Therefore the existence

of a solution to the CJS can be judged by whether I Vf | - O or not.
(I . | denotes the size of a set.) If an external node is saturated,

then it implies that the associated l-junction has been assigned a flow
input. Thus the adjacent port bond must have effort as input to the l-
junction. If an external junction node is free, then the adjacent port
bond has flow as input to the l-junction. Since a maximum matching
algorithm is applied, the number of free nodes in VE is a minimum. Hence
a maximum matching with all free nodes in set VE gives the minimum
number of flow inputs to the associated CJS, since each external node

has at least one environment edge (i.e., each external l-junction has at

least one port).

Now let us expand the junction gyrograph to include all the
environment nodes. Starting with the maximum matching obtained within
the junction gyrograph, a second maximum matching with respect to the

entire gyrograph may be done such that all the nodes in the set V are

E
saturated. In this process only two cases occur. First, a free node in
the set VE resulting from the first maximum matching can be saturated by
simply marking the incident environment edge. Thus the adjacent
environment node becomes saturated. Since one free external junction

node makes one environment node saturated. Fm. free external junction

in
node will make Fmin saturated environment nodes. In the second case, an
augmenting path starting at an environment node and ending at another
environment node is found. Reversing matching of this path will increase
the size of the matching set M by one and the size of saturated
environment nodes by two. If a maximum matching is found, then the

number of saturated environment nodes will be a maximum. It corresponds

the maximum number of flow input variables; namely, the F-maximum basis.

72

An algorithm for determining the existence, the F-minimum basis,
and the F—maximum basis for general junction structures is stated based!

on the above discussion.

Basis Ogder Algorithm for CJS

Assume an unmarked CC derived from a given CJS.

(1) Identify the sets VE’ VI’ and VEn'
(2) Start with the junction gyrograph CJ induced by the node set VB
and VI'
(3) Apply the Pape-Conradt maximum matching algorithm (Syslo, et.
al., 1983) to the junction gyrograph CJ.
* If ( |V§| = O ) then
a solution exists.
Else
Start at a free node v 5 Vi.
If ( there exists an alternating path which ends at a node

u c VE) then

reverse the edges. Therefore we have
f

f f .
VI = VI - v, | VI | - | VT | - l, and
f f f ,f
VE a VE'flL | VE |-—| \E | + 1

Go to (*) above.
Elseif ( there exists an alternating cycle which contains a
node u 6 VB ) then
reverse the markings of edges starting from the matched
edge incident with u to the edge incident with v. Therefore
we have

= - 7 = V - '
VI VI \, | V | | | l. and

f f f f

VE-VE+u, IVEI-IVEI+1
Go to (*) above.
Else
No solution exists. Stop.
Endif
Endif
(4) Fmin - | v; |. The minimum number of flow inputs to this

general junction structure is equal to the number of free nodes in the

external set.

(5) Apply the Pape-Conradt maximum matching algorithm to the entire
gyrograph CC.

(7) Fm - S The maximum number of flow inputs to this

V .
ax I En |
general junction structure is equal to the number of saturated

environment nodes.

All bases for a given CJS have F ranging between F.. and F .
min max
From the discussion on the maximum matching for the entire gyrograph we

also deduce that
F = F . + 2k (2.74)

where k --(l, 1, 2, ....., K and K is given bv K = (F - F . ) / 2
~ ma min

x
If the F . and the F of a general junction structure obtained by the
min max

above algorithm are the same, it indicates that this CJS has unique

basis order numbers (E, F).

74

2.5.4 Examples
2.5.4.1 Example 1

A bond graph with the environment nodes denoted by EN is given in
Figure 2-29a. Its associated gyrobondgraph is identical, auui the
gyrograph is given in Figure 2-29b. The nodes a, b, c and d, in due
gyrograph correspond to the junctions 1a. 11), Lc and 1d in the bond
graph (gyrobondgraph), respectively. The junction gyrograph CJ consists
of nodes a, b, c, and d. Nodes a, b, and c belong to the external set
and node d belongs to the internal set. Applying the maximum matching
algorithm to this gyrograph yielded no free node in the external set.
Thus Fmin - | v: | - 0. By using the result from the maximum matching in
ij we applied the maximum matching algorithm again to the entire
gyrograph (Figure 2-29d). There was one augmenting path {(A, a), (a, c),
(c, C)). Reversing this path yielded a augmentedlnatching set.P1-= {(A,
a), (b, d), (c, C)) (Figure 2-29e). Two environment nodes are saturated.
We found that Fmax = 2. The corresponding causality assignments for Fmin
and Fmax are shown in part f and g.rknjce that the difference between

the F . and the F in this example is 2.
min max

2.5.4.2 Example 2

An arbitrary gyrograph is given in Figure 2-30a. In this graph, VE

- (e. f. g. h. i, j, k), and V = (A, B, C, D. E. F.

= (a. b, c, d}. V En

I

C}. We first apply the maximum matching algorithm to its junction

gyrograph CJ which is induced by the set VEenKl\7. A maximum matching

I
may be found, namely, M = {(b. f). (C. g). (CL 11), (e, i), (j, h)).
Since Vf = { O I, a solution does exist (Figure 2-30b). The F-

I

75

RA‘—la‘—-CYd—lb —->RB

\ /
\CY\ /CY/

 

 

 

 

 

 

 

CY 1d CY d
\.:/

RC

(a) (b)
RA RA 3
RB d RB - d
RC RC C

(C) (d)

RA‘__41a|-— GYa—qlby—aRB RAF—lalt—- CY4-—1 lb i—bRB

(x. W. ya] (a? 267'
,

CY CY 1d CY

\I

C GY

1

1c 1c

T 1

RC RC

(e) (f)

Figure 2-29 Example 1
(a) Bond graph and gyrobondgraph
(b) Gyrograph
(c) Maximum matching in CJ
(e) Maximum matching in CC
(e) Causality corresponding to the F-minimum basis

(f) Causality corresponding to the F-maximum basis

76

 

 

 

 

 

 

 

En E
A a
B::
(Zn—a b n— i
D/
c
F
C 5 k
(a)
VI
e
f .
i
' , J'
h k
(b) (c)

Figure 2-30 Example 2.
(a) Gyrograph
(b) Maximum Matching of Cl

(C) Maximum Matching of CC

77

minimum basis was found from Fmin - | VE | = 1. Now we apply the maximum

matching algorithm to the entire gyrograph. An augmenting path is found;
namely, P - {(C, b), (b, f), (f, g), (g, c), (c, E)}. Reversing the
markings of the path made the environment nodes C and E saturated. The

number of saturated environment nodes were increased by two. Thus the F-

S

E | = 3. In this case the flow inputs

maximum is found to be Fm - | V

ax
will be applied to the junctions a, b, and c. Note that we may also find
another solution (a, b, d). Also we noticed that Fmax is greater than
F . by 2 in this case.

min

2.5.5 Remarks on the Algorithm

An algorithm for determining the basis order properties of an
arbitrary CJS has been developed through the use of the gyrograph
representation and the graph matching concept. The existence of a basis
order is determined by the absence of free nodes in the internal set of
the jruuztion gyrograph. The number of free nodes in the external set of
the junction gyrograph gives the F-minimum basis order of the CJS. The
resulting marked gyrograph is further used to determine the F-maximum
basis order by increasing the number of saturated nodes in the
environment set. The possible number of flow inputs.to a CJS ranges
between Fmin and Fmax in steps of 2. We conjecture that if Fmin - Fmax

in a CJS, no essential gyrators exist in the CJS. This remains to be

investigated.

The methodology for determining the basis order of a bond graph
with a general junction structure by using a transformed gyrograph may
be applied to the solution of the implicit R-field problem and the

implicit C(I)-fie1d problem. The identification of the basis properties

78
of the IJHMJC variables and determination of feasible inputs can
contribute to improving computational efficiency for large scale

nonlinear dynamic systems.

79

Chapter 3

MODEL ORDER REDUCTION METHODS

A simplified model of engineering system is always highly desired
in dynamic analysis, synthesis, and design. The simplification of a
dynamic system model brings many benefits to engineers. First,
repetitive simulations become easier and cheaper to perform, since the
computational load due to large dimension size and the widely separated
system time constants may be reduced dramatically. It happens in the
investigation of the influence on the system performance as some of
system parameters have been varied. Second, the complexity«of a higher
order model often makes it difficult to obtain a good understanding of
the behavior of the system. Salient features of the system, previously
hidden in a mass of detail, may be revealed. Third. controllers may be
designed for the reducedlmxkflq since most currently available control

design methods only work on small-dimension systems.

The simplification of a system model can be achieved by reducing
the model order, or by eliminating the adverse parameters which cause
the computational difficulties. Because of their importance in system
analysis and the design of controllers, model order reduction methods
have received considerable attention over the past three decades. In the
existing literature there are a number of books dedicated to this topic
(Decarlo and Saeks, 1981; Happ, 1971; Jamishidi, 1983; Michel and
Miller, 1977; Sage, 1978; Siljak, 1978; Kokotovic, Khalil, and O'Reilly,
1986) and numerous research papers. The objective of model order
reduction is to find a lower order model which preserves the dynamics of

more complex high order system in both time and frequency domains. The

80
techniques in the existing literature can be divided into two groups.
The first group of methods attempts to retain the dominant modes of the
original systems. It includes all aggregation methods and perturbation
methods. Another approach is based on applying an identificatdxni
procedure to input-output data obtained by driving the original system
with a specific input, for example, Walsh functions (Kawaji and
Shiotsuki, 1985). Since the latter is not widely applied in practice, we
will restrict ourselves to the first group. All the methods surveyed
below are concerned with time-domain models. Alternatively, the linear
time-invariant systems in state-space form can be represented in
frequency domain. By far the greatest effort in model order reduction
techniques based on frequency domain has been for single-imunit single~
output systems. We will not discuss these techniques here and in this

regard Jamishidi's book (1983) would be an excellent reference.
3.1 Aggregation Methods

The notion of aggregation was introduced in the control literature
by Aoki (1968). The intuitive idea behind the notion of aggregation is
quite simple. Suppose that S1 is a mathematical description of a
physical system using a given set of variables, and S2 is a consistent
description of the same system using smaller set of variables. Then S2
is termed an aggregate model for S1 and the variables of the system S2

are termed aggregate variables.

In the literature there are a number of aggregation methods. The
most basic aggregation method is the exact aggregation which illustrates

the notion of aggregation most clearly.

81

Consider the system

x(t) - Ax(t) + Bu(t) (3.1a)

y(t) - Dx(t) (3.1b)

Let z(t) - Cx(t), where x is a vector of dimension of n and z is a
vector of dimension of r. We want to obtain a new system model with

lower dimension r, i.e.

z(t) - Fz(t) + Gu(t)
= FCx(t) + Gu(t) (3.2a)

y(t) - W2(t) (3.2b)

By the requirement of consistency( or dynamic exactness), for any

u(t) and 2(0) - Cx(O), we need

CA - FC (3.3a)
CB - c (3.3b)
and WC 2 D (3.3c)

If the above equations are satisfied, then the triple (F,G,W) is
said to be a perfect aggregation of duetmiple (A,B,D) relative to the

aggregation matrix C.

Further insight into the nature of the class of matrices for which
dynamic exactness can be achieved is obtairmxilyy realizing that the
aggregation problem as posed for linear system is in fact a problem of

minimum realization. The exactness of the original system and the

82

reduced system only can be achieved when there are pole-zero

cancellation in the transfer function of the original system.

However, the aggregate state variables Zi(t) will not in general
correspond exactly to physical variables (Sandell, et al., 1978).
Therefore, an alternative point of view . which may be more useful in
applications, is to regard z(t) as an approximation to physical

variables. In other words in addition to (3.33) and (3.3b), we desire

z(t) z y(t) = Cx<t> (3.4)

where the matrix C picks out components or linear combinations of
components of x(t) that are to be approximated. Of cause, the choice of
the aggregate matrix C can greatly influence the nature of the

approximation.

There are many other aggregation methods such as controllability
matrix approach (Aoki, 1968), continued fraction method (Chen and Shieh,
1969), chained aggregation (Tee, et al., 1977), aggregation via
covariance equivalent realization (Yousuff, et al., 1985) and model
reduction via balanced state-space representation (Pernebo and
Silverman, 1982). All these methods are based on matrix similarity
transformation, therefore, they are not directly applicable to the

general nonlinear problems.

3.2 Modal Method

83
The modal method (Davison, 1967) is a well-known model order

reduction method. Consider a generalized linear system in the state-

space form

x(t) = Ax(t) + Bu(t) (3.5a)
y(t) - Cx(t) (3.5b)
Assume that A has distinct eigenvalues A“ A2, .H., Anvfith

negative real part and M is the modal matrix with columns consisting of
the corresponding eigenvectors of A. Define a new state vector u(t)
which is transformed via

_1
0(t) ‘ M X(t) (3.6)

then the system equation can be transformed into

n(t) = . n(t) + . u(t) (3.7a)

 

 

 

 

y(t) = [c,,c,, ...... ,cr]Tn(t) (3.7b)

Suppose tflnat the eigenvalues cluster into dominant and undominant
eigenvalues and that the undominant eigenvalues are asymptotically
stable. If there exist some bi's - O, we can say that the corresponding

states are not controllable, and if there exist some cj's = O, we can

84

say that these corresponding states are not observable.lhxm1input-
output consideration, the only important modes are those which are both
controllable and observable. Therefore, we may drop the equations
corresponding to the 1's and j's. Eliminating the uncontrollable and
unobservable modes leads to a minimal realization. By the same
reasoning, we may further approximate the transformed system by dropping
the states which are weakly controllable or weakly observable, that is,

the corresponding bi's and cj's are small.

Since the transformed system matrix has diagonal form, the new
system states are decoupled. We partition n according to the the dynamic
speed, controllability and observability into subblocks as fine

perturbation form:

fi1(t) - J1 01(t) + C1 u(t) (3.8a)
€2fi2(t) - J2 "2(t) + 62 u(t) (3.8b)
fi3(t) - J3 n3(t) + 6363 u(t) (3.8c)
in“) - J. n.(t) + G, u(t) (3.8d)
y(t) = H1n1(t) + H202(t) + H3n3(t) + €4H4fl4(t) (3.8e)

If we set 6i = O, i = 2,3,4, then
Mr) - J. mt) + c. u(t) (3.9a)

,1
y(t) - H1 01(t) - H2J2 C2u(t) (3.9b)

This method is workable theoretically. buttunzadvised in
practice. The computational effort is not trivial for large scale system
since the work involved in eigenanalysis OfiNlTIX.n matrix goes up as
us. The other trouble is from the numerical characteristics of large

scale systems which are ill-conditioned mathematically. Besides, if the

85

state variable of the original model x(t) represents physical variables,
the physical nature of the variables will be lost since n1(t) represents
nmthematical variables which in most case do not correspond to any

physical variables.

3.3 Lyapunov Function Method

Stability is one of the most important properties of’a dynamic
system. There are a number of qualitative analysis methods based on this
property for large-scale systems. Lyapunov's second method is an ideal
mechanism for accomplishing the aggregation plan in the stability
analysis of large-scale dynamic systems (Siljak, 1978). Actually, the
Lyapunov method itself can be viewed as an aggregation process. A
stability prOperty, involving several state variables, is entirely
represented by a single variable --- the Lyapunov function. However,
this approach simplifies the stabiliqyrnxmdem. but sacrifices detailed
information about the size of variatfons of each separate state

variable.

The concepts of vector differential inequalities and vector
Lyapunov function have been developed by Matrosov (1962) and Bellman
(1962) and other researchers. The concept associates with a dynamic
system several functions (say 5) in such a way that each of them
determine the desired stability properties iritflne system space ( of
dimension n > s ) wherever the others fail to do so. These scalar
functions are considered as components of a vector Lyapunov function,
and a differential inequality is formed in terms of this function, using
the orighuflqsystem of equations. As in the case of scalar Lyapunov

function, the stability properties of an n-th order system are

86

determined by considering only the s-vector differential inequality
Lyapunov functions. This can bring about a considerable reduction in the
dimensionality of a stability problem. It should be mentioned
immediately that there is no general systematic procedure for choosing
vector Lyapunov functions and that is the most serious drawback of the

approach.

3.4 Perturbation Methods

‘The other scheme of model order reduction for large-scale system
is perturbation, which is based on ignoring certain interactions of the

dynamic or structural nature in a system.

Perturbation methods are useful for dealing with a system that can
be approximated by a system of simpler structure. Mathematically, the
difference in the response between the actual and approximated systems
is modeled as a perturbation term driving the latter. In principle they
can be applicable to both linear and nonlinear problems. The
perturbations are divided into two classes of regular and singular

perturbations.

3.4.1 Regular Perturbation Method

The regular perturbations (Kokotovic. et al.. 1969) are those that
appear in the right-hand side of a diffinmmmial equathxu the general

formulation is the following form

x1(t) = f11(xl) + ef12(x2) + b1u1(t) (3.103)

x2(t) - ef21(x1) + f22(x2) + b2u2(t). (3.10b)

87
where e is a small positive parameter. The system is connected by the
small (weak) connections ef12(x2) and ef2,(x1) and it can be decomposed
into two completed independent subsystenm by ignoring the weak
connections. The computation of the two independent lower dimensional
problems is fewer than that of the single high dimensional problem. This

effect is enhanced for more than two subsystems.
3.4.2 Singular Perturbation Method

By singular perturbation is meant a perturbation to the left-hand
side of a differential equation (Kokotovic, et al., 1986). Consider a

dynamic system of the form

x(t) - f(x,y,t,e) (slow subsystem) (3.11a)

ey(t) - g(x,y,t,c) (fast subsystem) (3.11b)

where e is a small positive parameter. If we set 6 - 0, then the reduced

slow subsystem becomes

§(c) a f(x,y,t,O) (3.12a)

o - g(x,y,t,O) (3.12b)

If Equation 3.12b has isolated roots y a h(x.t).tfimzlimit model

for the slow system is

dx(t)/dt = f(x,h(x,t),t,0) (3.13)

88

The fast or boundary layer system can be obtained by stretching the

time scale from t to r = (t-to)/c as

A
A

SY - g(x,y,to+er,e) (3.14)
T
Under the condition that Rex(g$) s -p < O, the Tichonov Theorem can

be applied such that

x(t) - {{(c) + 0(a) (3.15a)

y(t) - y(t) + h(x,t) + 0(6) (3.15b)

Since the concept involves essentially an asymptotic approximation,
quantitative design results are difficult to obtain --ii:is hard to
say 'how small is small enough'. Furthermore. singular perturbation is
not a coordinate free concept,which may be the reason faréilack of
modeling procedures and prescriptive (computer-oriented) decomposition
techniques for model reduction via the characterization of the fast and

slow scales (Siljak, 1983).

As Sandell and Athans (1978) state: "From the practical point of
view, the main problem with this method is that a nKMhel of a physical
system is hardly ever given in the standard form with the slow and fast
variables separated and the parameter c conveniently appearing in the
left-hand side of the equations. It is a completely nontrivial exercise,
requiring considerable physical insight. to model a physical system with
slow and fast modes in the framework demanded by the theory!" This
problem lSIMNJHElly more severe for a poorly understood large scale

system.

89

3.5 Component Connection Methods

3.5.1 Component Connection Method

For reasons of efficiency it is often profitable to handle the
component equations of a system as separate entities. This permits one
to store the different component models separately in computer memory
and to analyze them one at a time. By decomposing a large scale system
into a number of smaller subsystems (components) and connecting them
based on their interactions, the component connection model can be used

(Decarlo and Saeks, 1981).

An interconnected dynamical system may be composed of many

components, each of which has a mathematical model of the following form

dxi/dt - fi(xi’ ai) (3.16a)

bi - gi(xi, ai) (3.16b)

where ai is the vector of input signals for duainth component;bi is
the vector of output signals for the i-th component; and x1 is the state

vector of the i-th component.

The interaction between the i-th component and the rest of

components is described by an algebraic equation

where u is the system input vector. From a theoretical point of

90

view, however, it is convenient to lump the n component equations

together, forming a single Composite Component Model. It takes the form

dx/dt - f(x, a) (3.18a)
b - g(x, a) (3.18b)
and
a = L11 b + L12 u (3.19a)
y - L21 b + L22 u (3.19b)

where x is the composite component state vector, a and b are the
composite component input and output vectors. respectively, and 11.and y

are the composite system input and output vectors.

The component connection model divides the system into two sets of"
equations: component equations, characterized by partially decoupled
differential equations, and the connection equation, characterized by

coupled linear algebraic equations.

There are two algorithms based on the component connectiorirmethod.
One is the Sparse Tableau Approach which "stacks" the various component
equations together with the connection equations to form a large, highly
sparse set of simultaneous equations.(fiyen an hunuzvector,1n and a
set of initial conditions, one can solve for a. b, and y by use of
sparse matrix inversion. The other one is the Relaxation Algoridmn,
which builds around a predictor-corrector integration scheme. It can
solve linear and nonlinear systems. The algorithm allows one to use a
different variable order and/or step-size integration routitua:for each

component of the system.

91

3.5.2 Diakoptics

It should be mentioned that research on decomposition-aggregation
methods was conducted by Kron in the 19505. He developed a scheme called
Diakoptics (Kron, 1963) whose main procedure can be summarized as

follows

(a) Tear the system apart into logical groups, each of which can
conveniently be analyzed as one unit;

(b) Set up the equations ofeuufixcomponent unit separately, as if
the other units were non—existent;

(c) Set up a "connection matrix" C showing how the various
components are interconnected;

(d) Using the matrix C with the laws of transformation of tensor
analysis, establish the equations of the interconnected system;

(e) Solve the equations piecewise.

It is rare to find a new research paper on Diakoptics in the recent
literature. Harrison pointed that "Diakoptics which has been successful
in solving large electric networks, turned out to be not as successful

in other type of models" (Harrison, 1972).

Karnopp (1970) also stated "All too often. however, the
mathematical style of Kron's presentation obscured even his basic
philosophy and many workers hntflmzfield of system dynamics elected to
ignore his contribution. If we examine the relationship of the
Diakoptics with bond graphs, a junction structure may be regarded as an)

almost physical representation of the tensors Kron discussed". However,

92
Kron's basic philosophy of decomposition-aggregation is still

significant to the research on the large scale systems.

3.6 Component Cost Analysis of Large Scale Systems

The performance of a dynamic system is often evaluated in terms of
a performance metric V. The performance metric V might represent the
system energy or a norm of the output error over some interval of time.
A question is "what fraction of the overall system performance metric \7
is due to each component of the system?" Based on the notion of
significance of system component in a dynamic system, the component cost
analysis for linear systems has been deveIOped by Skelton, et al.,

(1980, 1983).

Component cost analysis (CCA) consists of the decomposition of‘V
into the sum of contributions Vi associated with each component state

Xi' where the Vi's satisfy the cost—decomposition property,

V = 2 V. (3 20)

It seems equally natural and basic, therefore, to characterize the
system's behavior in terms of contributions from each of the entities in

the system. The CCA algorithm is summarized as:

Step 1. Determine a performance metric

v = lim E II y ll2 II y IIZ = yty (3.21)

t—*00

for the system

93

n n
xi = E Ai.x. +DiW y = E C.x. (3.22)
3'1 J J j=1 J J
Step 2. Compute V. from V.- tr[XCtC].. and
i i 11
0= XAt + AX + DDt ( Ricatti-type equation) (3.23)

Step 3. Rank the component costs in the manner

IV1|Z|V2|Z ..... ZIVI

The 'most critical' component of the system is xl having component
cost V1 and 'the least critical' component of the system isnx having
component cost Vn
Step 4. The least critical components will be deleted from the system.

The accuracy is controlled by the cost perturbatdxni index defined

by

k n

A = 2 V. / 2 V. k s n (3.24)
. i . 1
i=1 i=1

where k is the number of retained components.

3.7. Model order reduction in bond graph models

As has been suggested in the brief discussion about Diakoptics
above, there is in a bond graph model the potential for exploiting
efficient solution techniques. Here we mention two methods, one of which
is based on the Sequential Causality Assignment Procedure (Rosenberg and

KarnOPP, 1983), and is the standard approach for implementing

 

 

94

computational methods based on bond graphs (Van Dixhoorn, 1977; Granda,
1985; Rosenberg, 1984, 1986). The second method is newer and has not

been implemented computationally at this time.

3.7.1. Sequential causality assignment method

There is a close connection between the Sequential Causality
.Assignment Method and the Component Connection Method. When the
sequential causality assignment procedure has been followed for a bond

graph model, the system equations can be expressed in the form

Z1 - fi(X1, Xd) (3.25a)
Zd - fd(xi' Xd) (3.25b)
Do - g(Di) (3.25c)
U a h(t) (3.25d)
dxi/dt = $1121 + $12( Xd/dt ) + 51300 + SlaU (3.25e)
Zd - Slei + + $24U (3.25f)
Di - S3lZi + + S33Do + S34U (3.25g)
V - 84121 + 342( dXd/dt ) + 343Do + $44U (3.25b)

where X1 is the independent energy (state) vector:
Xd is the dependent energy vector:
U is the system input vector;
V is the system output vector;
Z. is the independent co-energy vector;
2 is the dependent co-energy vector;

D is the dissipation field input vector:

D is the dissipation field output vector

95
In fact, this method decomposes a system into the individual
elements, and then connects them by the junction structure matrix. It
solves not only linear problems, but also very general nonlinear

problems.

3.7.2 Reciprocal system method

There has been an attempt to connect the bond graph modeling
approach with singular perturbation theory (Dauphin-Tanguy, et al.,
1985). It defines the notion of a reciprocal system and then applies the
theorjrcaf perturbation. Thereby greater accuracy on the fast time scale
behaviors of the system can be obtained. However, one must construct a
reciprocal bond graph model. The analysis procedure can be illustrated

as shown in Figure 3-1.

81, S2, S3 and S4 are four different system models. The author
claims thattflds method can be numerically implemented without
udifficulty, even if some matrix.elements differ greatly in magnitude,

because no matrix inversion is required.

 

 

 

 

Reciprocal
transformation
Initial $1 > 82 Global system
global A
system
Singular
Reciprocal Perturbation
transformation
Fast S4 < S3 Slow
reduced decoupled
system system

Figure 3-1 Reciprocal system method

 

96

3.7.3. Finite mode distributed system models

Finite mode distributed system models (Margolis, 1984) are
extremely accurate in a chosen frequency range, while requiring only a
fraction of the number of equations required by finite difference
methods. The principal drawback is that the normal modes and frequencies
must be obtained before the modeling process begins. This can be a
tedious, if not impossible, task in itself. However, in many instances
the actual system being modeled is composed of nearly uniform structural
elements(such as beams,plates and membranes) for which the modes are
easily obtained. Another problem usually associated with finite mode
models is in the selection of appropriate boundary conditions for

determining the original system modes.

The decoupled modal equations are obtained as

J
mini + kini = .2 Fjwi(xj) (3.26)
j-l
where kin wizmi , mi is the modal mass, wi is the frequencies for the
unforced system, and m1 = m JD Wi2 dD.

The actual system displacement is computed from

W(x,t) - Z Wi(x)ni(t) (3.27)
i=1

The bond graph representation of Equation 3JNSis shown in Figure
3-2. Each of the retained modes is an I-C pair emanating finnn.a common

l-junction. Each I element is armxkfl.mass mi. while the corresponding

97

compliance, C, is an inverse modal stiffness or l/miwiz. The external

forcing is properly represented by the Se node associated with the
discrete force Fj(t). These forces are "felt" to each mode by
multiplication with the proper mode shape evaluated at the location of
the force. The TF moduli are nothing more than the mode shapes evaluated
at the location of the respective forces, i.e., TFll - W1(x1), TF12 -
W1(x2), ...., TF1k - W1(xk), ...., TFnk - Wn(xk) . Thus the
transformers (TFs) simply convert the actual force into modal forces
while at the same time converting the modal velocities, i71(t), into
actual system velocities at the force location. In this fashion, the
interactions between the continuum system and the lumped system are

established.

Figure 3-2 Finite modal bond graph model

98
Chapter 4
POWER POTRAITS OF DYNAMIC SYSTEM MODELS

All the model order reduction methods surveyed in Chapter 2 have a
common aspect; namely, the starting point is the explicit system
equations or transfer functions. For a physical system of large
dimension, iJzzis not easy work for an engineer to model it by a
mathematical equation representation. Bond graph technology provides a
potentially useful missing link between the physical system and its
mathematical models. It is a computer-oriented modeling language and can
be applied uniformly to many kinds of energy domains. In addition, it
provides a clear picture of system t0pology. A new approach to model
order reduction based on bond graphs is developed such that the model
order reduction can be considered without first having to obtain the
system equations. Let us next introduce the bond graph modeling

technique.
4.1 Graphical Representations of Dynamic System Models

In addition to representing physical system models by explicit
system equations, they may also be represented in graphical forms. Some
examples are schematic diagrams, block diagrams. signal flow graphs, and
bond graphs (Figure 4-1). Each of these representations has its own
strengths and weaknesses. Schematic diagrams depict the configuration of
physical systems in terms natural to particular physical domains, such
as electrical and hydraulic circuits. The mathematical relations of the

components are implied by the associated physical laws. The signal flow

 

 

99

F(t)

l __I

.s {,L...

////// /7

 

 

 

 

 

 

 

 

§(_S.)_..{ 6(8) ! Y(S) : y(s) - C(s) x(s)

 

1 s
x O / :O Y y(s) - (l/s) x(s)

 

 

Figure 4-1 Graphical representations of mathematical modeling

100

graphs are similar in certain respects to block diagrams. Only are the

block diagrams briefly discussed here.

Block Diagrams

Block diagram is a pictorial representation of the local cause and.
effect relationships among components of aimxkfl” Blocks may be
aggregated or made more detailed as need arises. Block diagrams provide
a convenient and useful representation for characterizing the functional
relationships among the various components of control systems. The block
diagram model stresses functional properties of modeled objects and
their signal connections. The block diagram is a signal-based modeling
metiuni. Blocks of the model (i.e., nodes of the graph) are connected by
directed lines which represent the direction of unilateral information
or signal fltnv. Block diagram modeling has been widely used in control
system design and simulation. In general, the block diagram is obtained
by represerming the particular mathematical equations, for example, the
control laws. Block diagrams reveal signal relations clearly, but the

power/energy aspects of a dynamic system are not readily accessible.

Bond Graphs

Power and energy attributes of models can be conveniently accessed
by employing a power-based model representation, namely, bond graph
model (Rosenberg and Karnopp. 1983). The bond graph is a pictorial
modeling representation based on power coupling among components. A
bond, the means of energy transfer between multiports,<xnwmcts two
ports. The bond graph is a power-based modeling method, since each tnnui

connector contains a pair of power variables whose scalar product is the

101
power. A half-arrow on the bond represents the positive reference power
direction. Using a rather small set of ideal elements, one can uniformly
construct models of electrical, magnetic, mechanical” lmydraulic,
pneumatic, thermal, and other systems, or mixed systems. Standard
techniques allow the models to be translated into a set of differential
and algebraic equations by hand or by computer. The bond graph model
depicts the physical effects considered by its modeler and their
topological relations. It is easy in practice to modify the model
structure to include additional effects. The bond graph technique has
found many applications in engineering, biology, and even economics (Bos

and Breedveld, 1985).

While bond graphs are excellent for modeling the 'plant' of a
system, they are not well suited for modeling the controls of the plant.
The signal communication aspect is better modeled by block diagrams. A
need for a system model containing both bond graph and block diagram
elements arises when the system being analyzed consists of subsystems
which are best approached with different formalisms. The mixed bond
graph/block diagram graph, developed by Zalewski and Rosenberg (1986)

will be used in this study.
4.2 Uniform Performance Measure ---_ Power

Power is often a neglected aspect of dynamic system response in
system dynamics. We may use a very simple oscillation system as our
illustration example. In Figure 4-2a a sclmummtic diagram depicts the
configuration and the linear parameters. By using Newton's law we may

derive the equation of motion as follows

102

mx'+bx+kx=F(t) (4.1)

where a superdot denotes a time derivative.

Figure 4-2b is the equivalent block diagram of the system. You may
well ask: "where is the power?" Neither the schematic diagram nor the
‘block diagram can give you the answer directly. However, aspects of the
power response can be accessed and displayed clearly in its bond graph,
Figure 4-2c. Each bond in the graph contains a pair of power variables,
effort e(t) (force) and flow f(t) (velocity), whose product is power. In
this example, the bond graph gives the following constitutive and

connective equations:

p2 - F2 (4 2)
V2 = m-xpg (4 3)
x,| - V2 (4 4)
F, - kx, (4 5)
F3 - bV3 (4 6)
F2-F1-F3-F, (47)

where Fi is force on bond i,

V. is velocity on bond i,

p is momentum on bond i. and

x is displacement on bond i.
Also associated with each bond in this bond graph is the power. Wi,
denoting the power on bond i, is the product of F1 and Vi’ Also the net

energy transfer, E. can be found from W. d1
1 t l

1(13

13(t)

 

l7rr77/777/l

b

 

(a)
-.1. ’5...
's

 

 

 

 

 

 

 

 

 

(C)

Figure 4-2 Powers in system models
(a) schematic diagram
(b) block diagram

(c) bond graph

’m—q

 

 

.Although the power variables (effort e(t) and flow f(t)) in
different energy domains have different physical meanings and units,
their product in every domain is power. As long as the efforts and flows
have appropriate units (e.g.,
products are the same (e g., watt). Table 4-1 summarizes the power
variables of the bond graph in several energy domains. We can see that
power is a uniquely uniform variable in a multi-domain system. Some

properties of physical systems may be investigated efficiently by using

such a uniform basis.

104

in the SI unit system),

 

 

Table 4-1 Powers in various energy domains
Domain Effort. e(t) Flow, f(t) Power, F(t)
Mechanical Translation Force, F(t) Velocity, V(t) F(t)V(t)
Mechanical Rotation Torque, t(t) Angular Vel., (u(t) t(t)m(t)
Electrical Circuit Voltage, V(t) Current, 1(1) V(t)l(t)
Hydraulic Circuit Pressure Volume Flow Rate P(t)Q(t)
1’0) Q0)
Magnetic Circuit Mmf Flux. Rate M(t) on)
W!) 00)
Thermal Temperature Entropy Rate T(t) S(t)
T(t) S(t)
Heat Transfer Temperature Heat Flow Rate T(t)Q(t) *

T0) Q0)

 

* Not a true power.

their power

105

It is natural at this point to raise the question, "How can we
exploit the knowledge about powers and energy transfer in a bond graph
model for better understanding of the system dynamics?" Answers may make
a valuable tool for model order reduction. The first step is to make the

power information readily available and accessible.

4.3. Power/energy Visualization

4.3.1 Computation of power variables

Refering the diagram in Figure 2-2 and 2-3 of Chapter 2, the source
field, energer storage field, dissipatMNiffield, and transducer fieLd
are coupled by a power-conserving connective multiport represented by
the junction structure. The Paynter junction structure consists of O and
l junction elements, (It is named after H. M. Paynter, tine inventor of
the bond graph.) and it is invariant. The modulated junction structure
is the collection of the transformers and gyrators, which may have
varying moduli. The key vectors of each field are identified. The system

equations are defined for each field as follows:

Source field U = ¢S(t) (4.9 )
Storage field Zi = ¢i(xi. Xd) (4.10)

2d = ¢d(Xi, Xd) (4.11)
Dissipation field D0 = ¢L(Di) (4.12)

Transducer coupling structure

 

106

To = ¢t(xi’ t)*Ti (4.13)

The invariant connective structure is represented by the junction

matrix and the associated vectors as follows

v - s * v. (4.14)
O 1
v - [z x D U T 1‘
1 i’ d’ o’ ’ o
1'.
v0 — [x1, zd, D1, v, T1]

where Vi is the input vector to the junction field and V0 is the output

vector of the junction field.

Furthermore, the variables of bonds incident to only 0 and 1
junctions are collected in a single vector, Y, which includes both
effort variables and flow variables. The connective relations in the
Paynter.rnumjon structure are defined by the square matrix P which
consists only of the integers O, l.euKi-l. The P matrix can be

decomposed such that

V a P .* V. + P * Y (4.15a)
o 01 1 0y

Y - P .* V. + P * Y . (4.15b)
Yl 1 YY

Since the coefficients of P are all constants,FHOVided dyarequired

inverse exists, we may find a solution for Y from

Y = I-P ' P . av. a. 6
(( yy) y1) 1 ( 1 )

107

Using (4.16) in (4.15a) we can solve the output vector for the Paynter
junction structure as
-1 ¢
Vo - (Poi + Poy((I - Pyy) Pyi)«~Vi (4.17)
Thus far, we are able to obtain every pair of power variables in a
bond graph model. Therefore we can evaluate the power response on bond i
simply by multiplying the proper pair of the two power variables;

namely,
Wi(t) - ei(t)*fi(t) (4.18)

The powers in a bond graph may be classified into two major types,
namely, external and internal. The external powers are those on the
external bonds which connect physical nodes (Se, Sf, C, I, R), while the
internal powers are those on the internal bonds which are incident only
to 0 or 1 junctions. The power and energy flow on the external bonds are

defined as follows:

Source elements

Q'— Sc (Sf)

W -= e*f (4.198)

"=I:?wm (4.1%)

The net energy transfer on tnelxnmirepresents the energy supplied

to the system during the time period from C1 to t2

108

Reeistagce elements

For passive resistors, the power w is always positive, i.e.
W. - e.*f. 2 0
1 1 1

The energy transfer on the bond represents the energy consumption

during the time period from t1 to t2

t2
Ti Itl ei‘kfi d1 (4 20a)
The total energy absorbed is

E a z T. (4.20b)

Capacitance elements

‘3 f
Cl 'n
__———-> C
f1

For springs and other types of capacitors. the power represents the

time rate of change of energv storage. Ifiit>C>the capacitor receives

 

 

109

energy, while for W < 0 the capacitor releases energy. The net energy

transfer on a bond is computed by

_ t2 _ t2
Ti Itl Wi d1 It: ei*fi d1 (4.21a)

The energy stored is given by

n
E-ETi+E
i-l

0 (4.21b)

Ine a e eme ts

For masses and other types of inertial elements, the power W is the
time rate of change of kinetic energy. If N > O the inertial element
receives energy, while for w < O the inertial element loses kinetic

energy. The net energy transfer on a bond is computed by

-tz =t2~.-
Ti It ”1 d1 It ei fi dr (4.223)

1 1

The energy stored is given by

T. + E (4.22b)

110

Internel_bonds
The power on an internal bond is the product of its power variables
Wi - ei* fi
4.3.2 Power measures

In a dynamic system the power on a bond is a function of time that
is positive when the power is in the half-arrow direction, or negative,
if it opposes the half arrow. The power history can be displayed in a
usual way to assist system analysis. Figure 4-3 is the power responses
versus time of the example system of Figure 4-2 for ixumat F(t) = 1.0,
'where W.Bl is the power on bond Bl, and so on. Note that W.B3 is always

greater than zero since it is dissipated from this system by the damper.

A
The power conservation 2 W.Bi - 0 is true at all time.
i=1

Power responses versus time show how the powers associated with
different types of physical effects in different energy domains vary.
However, for a large scale dynamic system its bond graph may contain
tens or even hundreds of external bonds and their associated physical
effects, plus many internal bonds. While it is possible to display every
power response history for all bonds in a plot, it is very difficult to
abstract useful information from the vast amount of data. An
alternative way is to display aspects of the power response on the bond
graph itself. Since the power can be evaluated over’aa time perirmi, an
averaged power over that period may refhum:the intensity of the local

interacthmnvdthin the global system. The bond provides the perfect

111

am.

 

 

Fri-i.

Ill---.-...sm . m
j.

~-¢ madman 5H EoumAm as» mo noncommeu nosom m-q

a
d

'9‘

'

.zemts

mi:
5.5

i... amé llmmé. -..- Nmé. 1.15.;
ms:

muswuh

Bzmmmm
sm . s

 

I‘ll l!‘

sans

 

mm . a -
.-, in a ms.s-
\~\ a... on .9; \\II//

ss . \~ .o o I
ml \ / \ .t / .. xx / .VN S
T. \. / .. \ / x I /
_letili IE... ..I iv. 313.... I... \ ..iiIIVlilll. \ ll!!! / iii]! “uni r it
all “I 11(1‘ (ll-I‘ll \ [I‘ll/‘\ \ fl!!!pl\\\\
{\sc /II\\ o f/ \s\ f/ \ / .\\
\t .. ...\ u , x .. .
..a .N / lu\\ on /.~ mm B
... N. .o. 4. I\
..x w a ms.s
a i a, mess
xx. mn<um

 

.- mm.“

112
place to display this intensity quantitatively by color, width of the

bond, or some other means.

For the purpose of possible decomposition and simplification of the
system model, a statistical measure is introduced as below. The
magnitude of the power reflects the strength of the interaction on a
bond. For a suitable time period the magnitudes of averaged powers on
all bonds are computed and compared to find the relative interaction

strengths. Some possible power measures are:

arithmetic mean W - (l/T) IT W d1 ; (4 23a)
absolute mean W - (l/T) IT | W ldr ; (4.23b)
and root mean square W = [ (l/T) IT W2 dr ]1(2 (4.23c)

For easmetaf implementation the approximate measures corresponding

to the above definitions are employed:

K
arithmetic mean W = (l/K) 2 Wk ; (4.24a)
k-l
K
absolute mean W = (l/K) Z | Wk| ; (4.24b)
k=l
- K 2 1/9
and root mean square W = [ (l/K) 2 Wk ] “; (4.24c)
k=l

where K is the total number of data stored and Wk is the power at the
step k. 13mg range for k can be subinterval within the computed range.
Each operation may be used to make the best display in different cases.

Powers and energy transfers on all bonds can be computed. The energy of

the external nodes may also be computed.

 

 

113

4.3.3 Scaling and color-coding

Color graphics techniques are being used in the practical worlds of
science and industries as well, to help manageamuiinterpret the vast
amounts of scientific and technical data. Graphics reveal trends and
relationships that would otherwise remain buried underimnnmains of
numerical detail. With the aid of computer, graphic information is
generated, analyzed, and displayed for the power properties in a bond
graph model of a dynamic system. A color graphical display of the power
strength can be shown on a computer screen. which gives user very clear

picture of the power distribution in the bond graph.

In the power display the color of a bond is determined by a
certain scale and classified into a finite groups. There are three
scales available in the computer implementation. They are linear,
logarithmic, and rank-ordered. All bonds are sortxxi into six groups.
each of which is assigned a color from red (highest power) to blue

(lowest power) in the descending order of magnitude.
4.4 Power Method for Model order Reduction and Model Simplification

Since any physical system is energy related, in addition to
the generalized momenta and displacements and the generalized efforts
and flows, the power and energy should also pnxruide insight into the
system dynamics. As discussed above the power on a bond is the product
of its two power variables, effort e and flow f. The ideal elements are
(zonnected by the bonds which are the energy pipes between the connected
elements or subsystems. The time integration of the power is the energy

flow through the bond over the time period. Power can be treated as the

...3

 

114

uniform indication of the interaction between two sides of the bond. It:
can be analogue to arielectrical system in which the current through a
line and the voltage across the line are important attributes, but their

product, power, also provides additional knowledge about the system.

The baskzidea behind the notion of power analysis is that the
importance (weight) of a physical effect or a subsystem to the global
system dynamics may be evaluated by the strength of its interaction with
the rest of the system. The averaged power at the bond could be used as
a measure of the strength, therefore, the significance of a physical

effect or a subsystem can be revealed.

The simplification of a dynamic model can be achieved 111 two
stages. The first is at the subsystem levelanmlthe second is at Hue
component level..At the subsystem level we seek to decompose a bond
graph model into a number of subsystems, based on physical insight,
linear/nonlinear separation, and the topology of the junction structure.
At the component level, the importance of each physical effect is

evaluated based on the power interaction with the whole system.

4.4.1 Simplification of physical effects

A power display on a bond graph as shown in Figure 4-4 gives a
clear picture of the power distribution with respect motflmzphysical
effects. A set of the weakest physical effects can be identified
immediately from the powers on the external bonds. Since the power on an
external bond reflects the strength of the interaction with the restz<3f
system, a measure of the power magnitude may suggest the importance of

the associated effect. Weak physical effects carting removed from the

115

system on a trial basis. To guide the removal of the weakest physical
effects, some criterion must be set according to the nature of the
effects, the power metrics, and accuracy of the approximation. One

potentially useful measure is the effect index, El, defined as follows:

811 - ( w. / a ) s 100 s (4.25)

where W1 is the absolute mean (or root mean square) power of'tfiue i-th
bond, and Whax is the maximum absolute mean (or root mean square) power
among the all bonds. The use of the effect index will be illustrated

through an example in Section 4.5.

Se-l—I-hllr—4-*C :k
2 3

Figure 4-4 Powers on external bonds

 

 

116

4.4.2 Decomposition and simplification of subsystems

Elimination of unimportant physical effects may yield a reduced
order model if the ignored effects include the capacitance or inertance.
A model order reduction in the form of decoupling also can be achieved
by decomposing the system into subsystems by removing the least
significant subsystems whose connecting bonds have the lowest level
powers. In this case, a out set of bonds with the low power levels
indicate the weak connections between the subsystems and the "mother"
system, even though the bonds inside of the subsystems may have high
level powers. For example, assuming that in Figure 4-5, bonds 1, 22, and
up to k are internal and they are cut bonds where the system is
decomposed into k subsystems 81’ $2, and up to Sk' By examining the
powers on these bonds we found that one of them, say bond 1, has very
small power compared with the rest of bonds. It shows that subsystem S1
has weak connection to the other subsystems: therefore we can simplify
the system by removing subsystem S1 with a lxnatrisk of removing
significant effect. The resulting computathmusvdll involve separated

subsystems.

Similar to the effect index, the connection index, Cl, is defined

as follows:
CIi - ( a / z wj ) a 100 a < 1.0 (4.26)
j

where W1 is the root mean square ( or absolute mean) power of the cut

bond of the i-th subsystem, Wj's are the power on all the cut bonds

connected to the junction structure.

—u

 

117

 

 

 

"Mother"
31 System sk

 

 

 

 

 

 

Figure 4-5 Powers on cut bonds

4.4.3 Procedure for model simplification

By using the above suggested effect index and connection index a

general procedure to simplify the system model in a bond graph framework

has been developed. The procedure is described as below:

1. Conduct a simulation under nominal conditions on inputs and

parameters.

2. Identify and remove insignificant subsystems:

a. By applying the connection index criterion, identify the

weakest cut bonds from all cut bonds. thereby the most

insignificant subsystems:

 

 

118

b. Remove the most insignificant subsystems.
3. Eliminate the unimportant physical effects:

a. By applying the effect index criterion, identify the weakest

external bonds, thereby the most unimportant physical effects;
b. Remove the most unimportant physical effects

The developed method would be best illustrated by some examples in

the following section.
4.5. Examples
4.5.1. A Radar pedestal unit

For illustration, let us model a radar pedestal unit and its
control. A sketch of the radar pedestal unit is shown in Figure 4-6. The
purpose of the drive motor and control system is to set the angular
position of the pedestal about its vertical axis as desired. The system
graph for pedestal position control is shown in Figure 4-7. The system
graph includes a bond graph part for the unit and a block diagram

portion for the feedback control.

The bond graph part models a DC drive motor connected through a
shaft and gears to a pedestal unit. The input voltage to the drive motor
is supplied by node SEE. The dominant effects in the electrical part of
the motor are the inductance (IE) and the resistance (RE). The field

current (signal I) is converted to a torque with no back-voltage effect.

 

 

119

IM and RM are the motor inertia and fricthnr The shaft has torsional
stiffness (CS). The gears are modeled as an ideal transformer (TF). The
pedestal load is composed of a rotational inertia (IP) and a friction
effect (RP). Two types of power bonds are used, the electrical set (El,
E2, E3) and the mechanical rotation sets (M1, M2, M3 for the motor),
(81, $2, $3 for the shaft), and (P1, P2, P3 for the pedestal). The bond
graph part, starting with the motor input SEE, can be used to

investigate the open-loop characteristics of the drive system and load.

The position feedback control part of the model includes a transfer
function, l/s, which integrates the velocity of the pedestal, w, to get
its position 9. Node SUM compares the position with the desired value,
REF, the error CV is sent to modulate the input voltage of the drive

motor .

 

 

120

mm m

not
(Inuit-um

 

 

 

ii

Figure 4-6 A radar pedestal unit

 

 

 

 

 

 

 

 

 

 

 

 

REF
FB
SUM TRF N
CV in
V

 

 

SEE AIE ASEM 41M “-5‘0 —->‘TF -—>1P

El [ Ml 31 93 PL F]
El H3 32 P2
E2 M2

[E RE [M RM CS [P RP

Figure 4-7 Bond graph for the radar pedestal unit

 

121

The state equations for the radar pedestal unit system can be

derived as follows:

 

 

 

P2 -(r2/12)p2 ' kq + (Km/11)Pi (4
9 (1/12)P2 ' (m/I3)p3 (4
P3 mkq ' (rs/13)P3 (4
o' <1/13>p3 <4
The system matrix A and the input matrix B are as follows:
-r1/Il O 0 O 'K
A = 0 1/12 0 -m/I3 O
O 0 0 1/13 0
a - [ K6 0 o o 0 1‘
The parameters for this system are given in Table 4-2.
Table 4-2 Parameters of pedestal model
IE = 0.1 henrys rE a 5.0 ohms
2
IM = 0.25 Lg.m rM = 0.3 N.s/rad
kS = 100.0 N/rad IP = 320.0 kg.m2
rP = 50.0 N.s/rad Ke = 1.0
Ml - 20.0 * I TF = 30.0

 

.27a)
.27b)
.27c)
.27d)

.27e)

Therefore the matrix A and B are evaluated as:

 

122

 

-50 0 0 0 -1
200 -1.2 -100 o 0
A - o 4 0 -0.09375 0
0 0 3000 -O.15625 0
0 0 0 0.003125 0
B = [ 1 0 0 0 0 1‘

The simplification by removing the weakest physical effects may be

seen clearly from the power distribution of the radar pedestal unit. An

RMS measure of the power on each bond has been computed and the

resulting list of values has been rank-ordered.
response of the system is listed in Table 4-3. The

lists the RMS power value as a percentage of the

The abstracted power
column at the right

maximum value. Figure

Table 4-3 Effect indices
RMS
ORDER BOND POWER
1 M1 0.1059E+02
2 M3 O.8120E+01
3 S3 0.4030E+Ol
4 Pl 0.4030E+01
5 SI 0.4023E+01
6 P2 0.3739E+Ol
7 M2 O.3239E+Ol
8 P3 0.1506E+01
9 $2 0.2620E+00
10 E3 0.8605E-01
11 El 0.8662E-01
12 E2 0.17llE-Ol

El

100.
76.
38.
38.
37.
35.
30.

.210

000
642
041
041
969
287
575

 

.473
.818
.818
.016

H
OOONL‘

123

Figure 4-8 shows the root mean square powers displayed in six
colors. We found that bonds El, E2, and E3 have the lowest power levels.
However, these are the powers in the controlled field of the motor. The
major power supply of the system is from the armuture. The connector I
separates the two parts. In the main part bond 82 has the smallest
effect index with 2.473%. The physical effect associated with bond 52 is
the shaft compliance. It is expected that the shaft transmits large
torque with small rate of change of the torsional displacement under a
high oscillation frequency. It shows that the compliance of the shaft is
not significant with respect to the general system dynamics and suggests
that the ignoring the element CS would not degrade the model. The
simplified model is shown in Figure 4-9. Since the removal of node CS
causes derivative causality on one of the bonds M2 and P2, the two
inertance IM and IP are combined into a single equivalent inertance IEQ.
The order of the simplified model is reduced by two. The eigenvalues of
these two models are listed in Table 4-4 and their responses are plotted
in Figure 4-10. It is not surprising that the modes with high frequency
are truncated and the responses are almost identical except for the

ripples on the curve of the angular velocity of the original system.

 

124

Hmumovod on» mo :dmuw econ ecu co amaamfim umsom mic ouawwm

mm mm mm Z_.ng : mg n:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

ni. mm. . mm .
is mm m _ a: i . .qm :,
eiluZisii Ali 11:... minimally
>
.i . 3 >9
GEE: Fez _
mm.» mg;
mammmK .
_ zumZ v.2:m
N$uWww.mH w mm“ g
arms; as
gawkd mafi—
smmzwmefl .rI
ems.
.mduweww
maxi
. swag

 

 

125

 

 

 

 

 

 

 

 

 

w SRC
REF
F8 ~ .
SUM TRFN
CV '

 

W
SEE ——‘1-15 -—>SEM —> m —>@ —\n= -—> 1p

5 I N! St 53 P1
53 N3 P2 P3
E2

IE RE RM IEO RP

 

Figure 4-9 Reduced-order model of the pedestal

Table 4-4. Eigenvalues of the original and reduced models

 

 

Original Reduced
-3.862lOE-Ol i j 2.60896E+01
~2.91432£-Ol i j 3.67950E-Ol -2.91340E-Ol : j 3.67800E-Ol

-5.00010E+01 -5.000hSE+Ol

 

 

 

 

sa.m sm.d

 

mHovoa vmosvou can HmcHwauo mcu we mmmcoammm o~-¢

 

mZHF
am.s

a mad*
am.“

musmfim

...,- m: 1.1%.“. a_,.._m3._

s¢.g

 

126

 

 

 

 

 

sm._
m- wad*
mqaum

as.m

 

 

 

 

127

4.5.2 An Euler-Bernoulli beam with a vibratory load

As a second example consider a pinned beam with a vibratory load
attached at position x1 and a driving force F(t) acting at position x2,
as shown in Figure 4-11. The beam acts as a coupling element between a
driving force and the load. Following the modal bond graph approach of
Karnopp and Rosenberg (1968) and Margolis and Tabrizi (1984) the
continuous beam is modeled by a modal bond graph with five modes
retaixuxi (See Figure 4-12). Mode 1 is represented by C1, 11; mode 2 is

represented by CZ, 12; H”. The modal frequencies of the beam are given

by

2 __
= ’ 4
”i (1 " ) J EI/pAL

where “3.15 the i-th modal frequency; E is the Young's modulus; I is

the area moment of inertia; p is the density of the beam material; and L

4
is the length of the beam. In this example, we assume / EI/pAL - 1.0,
2 2 2
so that (.01 - n - 9.87, w2 - (2x) = 39.48, w3 = (3n) - 88.83, w4 -
2 2
(4n) - 157.91, “5 - (Sn) - 2A6.74. The model mass mi is determined as

follows

E
I

L 2 L 2
I pAYi(x) dx = I pA sin (inx/L) dx
0 O

1 2
J pAL sin (inx/L) d(x/L)
O

pAL [ l - cos(inx/L)]/2 d(x/L)

pLA/2 = M/2 (4.29)

 

 

128

33

 

 

-F_l—§>_
L-
$2

 

 

 

Figure 4-11 A beam-coupled system

31
[U PM 39 [Y “CD

127 t“

1X RD SE
l” l“
,aA as
,/’/’ 3 5
,/’ 3\ / n
../ 32127 PR 2 nfllk
rm TFIZ ma m4 ms TF21 TFZZ/TFZB TF24 was

  

M R R L’\

Cl [1 [2 I2 C3 [3 C4 [4

53<1;'fii
%

Figure 4-12 Modal bond graph model for the beam-coupled system

 

 

 

129

where M is the mass of the beam. Let M - 10.0 kg; then 1111 = m2 - m3 - m4

- m5 - 5.0 kg. The modal stiffnesses are computed by

k1 - mlwl 487.08 N/m

k2 - m2W2 = 7793.4 N/m,

k3 - msws 39453.8 N/m,

124677.8 N/m,

w
‘
I
a
a
8
a
I

Rs 3 msws 304403.1 N/m.
The parameter values for the lumped mass and spring are chosen to
‘be m - 2.() 'kg and k.- 200.0 N/m. Thus the load natural frequency is

w = 10 0. For a damping ratio 0.1, we choose b = 4.0 Ns/m .
The mode shapes at position Xj are determined by
Y; = sin (iji/L) i = 1, 2, ..... ; j = 1, 2, .... (4.30)

The mode shapes at the positions of the predescribed force F(t) and

the load are calculated from Equation 4.30 when x = 1405 and x. == 2L/3,

l 2

where L is the beam length. They are
Y1 = [0.5878, 0.9511, 0.9511. 0.5878. 0.0]‘,
Y2 = [0.866, -O.866. 0.0, 0.866. -O.866]t
The RMS power distribution for a unit step force ixunit (E.42) in
the bond graph is given in Table 4-5. The data are coded on the bond
graph in Figure 4-13. Bonds 30 and 20 have zero power because tflua load
is located at a node of mode 5 (x1 = L/S). Bonds 33anu123 have zero
power because the input force is located at a node of mode 3 (i e., x2 =

2L/3). The lower parmzof the bond graph is the modal representation of

 

 

130

Table 4-5 The power listing of five-mode model

RMS

ORDER BOND POWER EI %
1 1 0.6239E-02 100.000
2 42 0.6181E-02 99.079
3 11 0.5808E-02 93.098
4 21 0.5783E-02 92.692
5 31 0.5783E-02 92.692
6 2 0.2155E-02 34.535
7 3 0.2071E-02 33.188
8 l2 0.1931E-02 30.956
9 22 0.1892E-02 30.327
10 32 0.1892E-02 30.327
11 38 0.1587E-02 25.443
12 39 0.1491E-02 23.900
13 41 0.1403E-02 22.487
14 36 0.8567E-03 13.731
15 37 0.8567E-03 13.731
16 4 0.7426E-03 11.903
17 16 0.7200E-03 11.541
18 26 0.7200E-03 11.541
19 7 0.6410E-03 10.275
2O 14 0.5910E-O3 9.474
21 34 0.5850E-03 9.378
22 24 0.5850E-O3 9.378
23 40 0.476OE-03 7.630
24 9 0.4250E-03 6.813
25 17 0.3908E-03 6.264
26 27 0.3908E-03 6.264
27 15 0.3871E-O3 6.205
28 25 0.3871E-03 6.205
29 35 0.3871E-03 6.205
30 8 0.2477E-03 3.970
31 10 0.176OE-O3 2.821
32 29 0.7394E-04 1.185
33 19 0.7394E-04 1.185
34 5 0.7855E-05 0.126
35 13 0.7753E-05 0.124
36 18 0.7753E-OS 0.124
37 28 0 7753E-05 0.124
38 6 0.3685E-06 0.006
39 20 0.0000E+00 0.000
40 23 0.0000E+00 0.000
41 3O 0.0000E+00 0.000
42 33 0.0000E+00 0 000

131

 

Hopes ammuw “Eon .2..on on”. H: Gowusnfluumfiu uwsom Mala «.5me

 

 

 

 

Hump

oummmo Hoz _ mHLH NWLH

as,mm©.m
wa,mam._
wa-mw+.m
ms-mma.4
éa-mma.m .A
ms.mmm.u outr-runwuXSJwvm:
ma-d4m.m
:
wzm
mmxom

 

 

 

.CC _ .. w. 7x33

921
aiwfinlw:

 

Hopes :mmuw vaoo. H3808 ma... ca Coausnfiuumflv uozom 2|.» muswflm

131

adam_o eoz _
as.mm©.m
wa-Mam.u
wa.wc4.m
ms-mmm.+
as-mme.m
SEN; oodlcllélmle
ma.dam.m
.905 “m5<um

 
  

 

 

 

 

 

 

132

the beam and consists of five modal pairs of C and I nodes. Each pair
may be treated as a subsystem. The subsystem 1 composed by C1 and 11,
and so on. These subsystems connect with the mother system by bonds 11,
12, 13, 14, and 15, respectively. These bonds are cut bonds and consist
a cut set. We can apply the connection index CI to judge the
significance of each individual subsystem. First, we compute the

connection indices as follows

CI - 011 / (011 + 012 + 013 + 015 + 015)

- 0.005808 / (0.005808+0.001931+0.000007735+0.000591+0.0003871)
- 0.666

CI - 0.221

CI - 0.000887

CI - 0.0677

CI - 0.0440

By comparison of the obtained connectjxni indices in Table 7, we
found that the third mode represented by subsystem 3 has the weakest
connection. Therefore, removal of mode 3 should not have much effectzcn1
the performance of the whole system. We may try to remove mode 5 and 4,

sequentially.

Table 4-6 The powers on cut bonds

RMS
ORDER BOND POWER CI %
1 11 0.5808E-O2 66.567
2 12 0.193lE-02 22.131
3 14 0.5910E-03 0.887
4 15 0.387lE-03 6.774
5 13 0 7753E-05 4.436

”'3

 

133

If a model simplification that ignores the highest three modes is
used ( See Figure 4-14), then the resulting RMS power distribution in
list format is given in Table 4-7. The eigenvalues of the original and
the reduced models are listed in Table 4-8. To further investigate the
behavior of the simplified model as compared to the original model we
plotted the velocities at the input (F.42) and the load (F.37) locations
versus time in Figure 4-15 ('F' denotes flow, hence velocity in ENPORT).
The curves for reduced model are more smooth. The deviations are quite
small. This example shows that the weak connection of a subsystem to the

mother model is a good indicator to lead to a model simplification.

The computational effort is reduced dramatically by the model
reduction. The data in Table 4-9 were obtained by running the models in
a Prime 750 system using ENPORT7-7.l.5. The computation time is
decreased by reducing the number of state variables and eliminating the
high frequency modes. It is most clear when both the fourth and fifth

modes were eliminated.

I
ID-<——ox-;;1v-;fco
37 Lu

1X RD E
l” 1"
’95 .3.B
;;?/ .-vi)”
A 27 .-z “/33
TFll TFIZ ‘fIF21 TF2?

 

1AA“' 188

t: Lu
1A 18

LC\3 Lf\{\
‘\

Cl [1 C2 12

Figure 4-14 Reduced model with two modes

‘5

 

Table

Tab

1e

134

4-7 The power listing of two-mode model

ORDER BOND

\ONVChUlbU-DNH

1

42
ll
21
31
2

3

12
32
22
38
39
41
36
37
A

16
26
4O
17
27

OOOOOOOOOOOOOOOOOOOOO

4-8 Eigenvalues of the

RMS
POWER

.6248E-02
.6147E-02
.5817E-02
.5792E-02
.5792E-02
.2154E-02
.2051E-02
.1913E-02
.1875E-02
.1875E-02
.1582E-02
.1486E-02
.1398E-02
.8515E-03
.8515E-03
.7340E-03
.7197E-O3
.7197E-03
.4739E-03
.387lE-O3
.3871E-03

%

.000
.379
.101
.698
.698
.469
.822
.626
.010
.010
.316
.782
.374
.628
.628
.748
.518
.518
.586
.196
.196

original and reduced models

 

 

Original Reduced
0.00000E+00 i j 2.46739E+02
-l.40411E-Ol i j 1.57952E+02
-3.75707E-01 i J 8.90238E+01
-4.12466E-Ol i j 3.99433E+01 ~4 17791E-01 i j 3.99438E+Ol
-8.06808E-01 i j 1.17678E+01 -8.l7999E-Ol i j 1.17833E+01
-2.64687E-01 i j 8.24499E+00 -2.64250E-01 i j 8.25372E+00

 

Table 4-9 Computer times

 

 

Modes in CPU time
model (sec)
1,2,3,4,5 27.809
1,2,4,5 20.488
1,2,4 13.591
1,2 2.955

 

 

.moauuooao> pooa new name“ mo comauwano ma-¢ ouswwm

---mu.ll.; “ozmowu
5 mad“ wrap

135

 

\

‘

‘3
‘

 

 

 

ss+wss.s
ssmwm._
mu

 

 

 

 

 

 

 

 

 

 

 

 

ss+wsc.s
us-msm.m
;

 

 

.wzflu<um

136

From the examples above the blue color on the thinnest bonds
indicates the lowest power level and the user can easily spot those
bonds qualitatively, hence the least significant physical effects.
Referring to the power listing, the user obtains the power levels
quantitatively. Such process is very natural and intuitive. The
simplification may be done repeatedly by eliminating the lowest power
bond and associated node, or subsystems. In fact, the beam-load model
experienced three times of reduction, each time a pair of modal C and I

nodes were removed.

Note here that the two examples are linear and time-invariant. It
is for the purpose of making a comparison by using the system
eigenvalues. It is observed that in this power-based method there is no
functional analysis or matrix transformation involved. The method is
IMJlIC on the solution data obtained in simulation, therefore, no
restriction exists with respect to the linearity of time system in the

problem.

It is clear that since this method must use the simulation data to
find the power distribution, we can not escape the initial computation
for solving the original high-order, complex, or ill-modeled system
model. It may prevent the application in the case that the computer does

not have enough storage space for large number of solution data.

 

 

137

Chapter 5
SUMMARY AND CONCLUSIONS

In this dissertation two approaches to the design assessment and
simulation of engineering systems were presented. The principal research
objective of providing tools for improving insight and efficiency in

dynamic system design was achieved.
5.1 Implicit R-field Simulation

‘With respect to the simulation of models containing implicit R-

fields, our research has accomplished the following.

1. For a bond graph containing IRES, the IRFs were identified at
the graph level from causality data. Each IRF was identified separately.
The total computation effort was reduced by this localization of IRFs.

The identification procedure was implemented in software.

2. A new algorithm that finds the minimum iteration set of a given
IRF was presented. This algorithm uses a;finunfion structure property,
the basis order, to determine the minimum number of iteration variables
and applies to IRFs with a simple junction structure and l-port R nodes.
The minimum number is equal to the minimum of the pair of basis order
numbers, E and F. This number is always less than or equal to one half

the number of bonds on the IRF R nodes.

138

3. To extend the above algorithm to an IRF containing multiport R
nodes and a simple junction structure, a rule to determine the minimum
number of iteration variables for a given causal IRF was deduced from an
equivalent digraph. The coupling of port variables by multiport R nodes

increases the size of the minimum iteration variable set from min {E, F}

M
to min {E, F} + 2 min {Nr‘ Ng 1m . The summation term reflects the
m-l

couplings within multiport R-nodes.

4. Next an acausal IRF with a multiport node and a simple junction
structure was considered. The question was "What assignment of causality
will lead to a minhmnnileration variable set?". A strategy for
.assigning causality such that the number of iteration variables for an
IRF containing one multiport R node is a minimum was developed. This

strategy also determines such a set of iteration variables.

5. The algorithm described above utilizes the property of effort-
to-effort and flow-to-flow transformations of simple (and weighted)
junction structures. For general junction structures (containing one or
more essential gyrator nodes) we developed a new result that
gives the basis order. The basis order numbers, E and F, are not a
unique pair. In order to find the range of basis order numbers, we
transformed a general junction structure to an associated gyrogramdm. By
using the gyrograph representation and the graph matching concept, an
algorithm for determining the F-maximum and the F-minimum of basis order
for a general jLuurtion structure was developed. The structure of the
basis set was determined to range from Fmin to Fmax by steps of tn“) in

F.

 

 

139

The author makes the following suggestions for future research on

this topic.

1. The restriction on the number of multiport R node to ones in the
algorilflun for optimum causality assignment in IRFs should be relaxed.

The interactions among multiport R nodes must be considered.

2..App1y the basis order algorithm for general junction structure

to the solution of IRFs with a general junction structure.

3. Implement the algorithms from 1 and 2 above within a bond-graph-
based simulation program, such as ENPORT-7. It will make the software

more powerful and efficient in design assessment and simulation.

5.2 Power Distribution and Display

To increase insight and reduce simulation time an approach to model
order reduction and simplification based on the power distribution
response of dynamic systems was investigated. A new display tool was
developed and implemented in software. A number of examples were
provided to illustrate the use of the tool. To obtain power distribution
attributes the engineering system is modeled by a tnnui graph, a
simulation is made under a set of given test conditions, and time
histories of system variables, including powers. are obtained. Various
measures of the power response can be calculated and displayed directly
on the bond graph by color coding. Such displays give engineers an easy'
‘way to visualize power distribution attributes in large-scale, complex
dynamic systems under a given set of operating conditions. Often insight

may be gained into possibilities for model simplification.

 

140

A major advantage of the approach is that it is associated with the
graph model and does not require the engineer to deal with equations.
The power response approach is not restricted to linear models. In the
nonlinear case there are no convenient data of the class of eigenvalues
to which to appeal for insight. Here the power distribution may offer an
efficient way to make a trial assessment of possible model

simplification.

The following issues remain for future development with respect to

the power distribution method.

1. It will be useful to get the existing tools for power response
assessment, now implemented in pilot version, into wide use so that data
regarding modeling approximations can be accumulated. From such shared
experience it may be possible to find guidelines that are automatable

and will provide the engineer with additional insight.

2. We have observed that different power measures, such as mean or
RMS, can yield different insight into the relative importance of parts
of the model. With respect to possible model simplification we are at
the "cut and try" stage. It would be valuable to have some guidelines as
to how to exploit power level distribution data with respectztua system
response. One such question is "What different measures should be
considered for different physical effects, such as power sources, power

dissipation, and energy storage?"

3. Power distribution properties may open a new research (Lirectitni
in control. For example, are there some ways to control a dynamic system

based.on not only signals but also power properties? To answer this

 

 

141

question we need to study the power properties in a dynamic system
theoretically. The relation of power variables to many important aspects
in control, such as controllability, observability, and stability of a
system, need to be investigated. One may start with a simple linear

time-invariant model.

APPENDIX

142
APPENDIX : Standard Bond Graphs and Gyrobondgraphs

1. Standard bond graphs

A.standard bond graph consists of elements of the standard set
(Table A-l): C, I, R, Se, Sf, 0, 1, TF and CY, which are called
capacitance, inertance, resistance, source of effort, source of flow,
one-junction, zero-junction, transformer and gyrator, respectively. The
elements C, I, R, Se and Sf represent the energy field effects of
lumped-parameter multiport physical components, while the elements 0, 1,
TF, and CY are used to form the junction structure of the bond graph.
The junction structure defines the connection pattern among the field
multiports, and it is power-conserving. The SJS contains only elements 0
and 1, while the WJS contains the elements<3,]” and TF. They share a
common causality structure, in the sense that the effort transforms to
efforfl:‘and flow transforms to flow. The efforts and flows are disjoint.
Their causal properties are clearly established. The CJS contains
gyrator (CY) as well as the elements 0, l, and TF. The properties of CJS

have not been established.
2. Gyrobondgraphs

The gyrator is a fundamental element in constructing power dual
pairs of elements, such as (C, I), (Se, Sf) and (O, 1). Based on the
gyrators, the standard set can be reduced to a smaller working set (I,
1!, Se, 21 and CY). Bond graphs can be transformed into a canonical form,
gyrobondgraph, from which the properties of the original bond graph cant

be deduced. For example, see Rosenberg [9]. hitfds paper, we will use

143
the set (I, R, Se, 1 and CY) as our primitive set for expressing the

gyrobondgraph.

A standard bond graph can be transformed to its associated

gyrobondgraph by the following procedure:

(1) Replace each standard element not in the primitive set by its
equivalent:
(a) an Sf is equivalent to an Se and a CY;
(b) a C is equivalent to an I and a CY;
(c) a TF is equivalent to two CYs in cascade (one of which has
unity modulus), and
(d) a O is equivalent to a 1 and a CY at each port.
(2) Eliminate all pairs of unit gyrators.
(3) Combine all l-junctions that are directly bonded.
(4) Insert 1-junctions to meet the adjacency conditions, namely,
(a) each I, R, and Se is adjacent to a l-junction; and
(b) each CY is adjacent to two distinct l-junctions.
(5) Combine any fragments of the graph that require it into

equivalent 1- or R- fields.

3. Gyrographs

The vmnxi, gyrograph, was first used by Professor H. M. Paynter in
1968. He stated that all reduced gyrobondgraphs contain ideal multiport
junction structures consisting solely of interconnected 1's and GY's,
and all such systems are but specializations of a general ideal
multiport we call gyrostructure. A graph, simplified from its

gyrobondgraph is called the gyrograph.

144

“To simplify'a gyrobondgraph to a gyrograph the following steps are

used

(1). Replace each l-junction by a corresponding circle.

(2). Replace each gyrator together with its bonds by a edge joining

the two adjacent circles.

(3). Replace each external node by a square.

‘Table

SYMBOL

55...£L.,

c¢L_%__.

'r"

R f

A-l Standard node sets

DEFUNITION

e - e(t)

f - f(t)

e - 4K0)
e(t) - u(toh {fodt

f - 4K9)

u(t)-p00» {e-a:

0(30f) a O

8‘ ' M082
m-f' I f2

e, - refz

82 . nfl
f' + f2 - f3 3 0
f1 8 f2 3 f3

NAME

source of effort

source of flow

capacitance

inertance

resistance

transformer

gyrator

common effort

junction

common flow
junction

BIBLIOGRAPHY

145

BIBLIOGRAPHY

Aoki, M., "Control of Large-scale Dynamic Systems by Aggregation," IEEE
Trans. Automatic Control, AC-13, pp. 240-253, 1968.

Bellman, R., "Vector Lyapunov Functions," SIAM Journal of Control, Vol.
1, pp. 32-34, 1962.

B05, A. M. and Breedveld, P. G., "1985 Update of the Bond Graph
Bibliography," J. Franklin Institute, Vol. 319, No. 1/2, pp. 269-286,
1985.

Breedveld, P. C. , "Essential Gyrators and Equivalent Rules for 3-port
Junction Structures," Journal of Franklin Institute, Vol. 318, No.2, pp.
77-89, August 1984.

Burreto,;L” and Lefevre, J., "R-fields in the Solution of Implicit
Equations," J. of the Franklin Inst., vol. 319, No 1/2, pp. 227-237,
1985.

Chen, C. F2, and Shieh, L. 8., "Continued Fraction Inversion by Routh's
Array," IEEE Trans. Circuit Theory, CT-l6, pp. 197-202, 1969.

Davison, E. J., ”A Method for Simplifying Linear Dynamic Systems," IEEE
Trans. Automatic Control, AC-12, pp. 119-121, 1967.

Dauphin-Tanguy, G., Borne, P. and Lebrun, M., "Order Reduction of Multi-
Time Scale Systems Using Bond Graphs, The Reciprocal System and Singular
Perturbation Method," Journal of the Franklin Institute,\kfln 319,
No.1/2, pp. 157, 1985.

Decarlo, R” A., and Saeks, R., Interconnected Dynamical Systems, Marcel
Dekker Inc., New York, 1981.

Van Dixhoorn, J. J. , "Simulation of Bond Graphs on Minicomputers," J.
Dynam. Syst., Meas. & Control, Vol 99, pp. 9-14, 1977.

Granda, JR ;I., "Bond Graph Modeling Solutions of Algebraic Loops and
Differential Causality in Mechanical and Electrical Systems," Proc.
Applied Simulation and Modeling” pp. 188-193, IASTED Conference, San
Farancisco, CA, 1984.

Granda, JR ;l., "Computer Generation of Physical System Differential
Equations Using Bond Graphs," J. Franklin Inst. , Vol 319, No.1/2, pp.
243-255, 1985.

Guardabassi, G., and Locatelli, A., Large Scale Systems: Theory and
Applications, IFAC, Pergamon Press, Oxford, 1976. ‘

Happ, H. H., Diakoptics and Networks, Academic Press, New York, 1971.

Harrison, B. K., "Large-scale Linear Systems, in: System Theory," L. A.
Zadeh and E. Polak (Eds.), McGraw Hill, New York, 1972.

146

Hood, S. , Rosenberg, R., Withers, D., and Zhou, T., "An algorithm for
Automatic Identification of R-fields in Bond Graphs," IBM Journal of
Research and Development, Vol. 31 , No. 3, pp. 382-390, May 1987.

Hrovat, D., Tobler, W., and Tsangarides, M. , "Bond Graph Modeling of
Dominant Dynamics of Automotive Power Trains," Dynamic Systems: Modeling
and Control, ed., M. Donath, ASME DSC-Vol. 1, 1985.

Jamishidi, M., Large-Scale Systems: Modeling and Control, North-
holland, New York, 1983.

Karnoppu I). C. and.Rosenberg, R. G., "Analysis and Simulation of
Multiport Systems," MIT Press, Cambridge, MA, 1968, See Section 4.3.

Karnopp, D. C., ”Power-Conserving Transformations: Physical
Interpretations and Applications Using Bond Graphs,"'.I. Franklin
Institute, Vol. 288, No. 3, pp. 175-201, 1969.

Karnopp, D. C., ”Correspondence on Bond Graphs and 'The Penetrating but
Olympian Tensorial Methods of Kron' ," The Matrix and Tensor Quarterly,
pp. 120-124, June 1970.

Kawaji, S., and Shiotsuki, T., "Model Reduction By Walsh Function
Techniques," Mathematics and Computers in Simulation, Vol. 27, North-
Holland, 1985.

Kokotovic, P. V., Perkins, W. R., Cruz, J. B., and D'Ans, G., "c_
coupling Method for Near-optimum Design of Large-scale Linear System,"
Proc. IEE 116, pp. 887-892, 1969.

Kokotovic, P. V. , Khalil, H., and O'Reilly, J., Singular Perturbation
Methods in Control: Analysis and Design, Academic Press Iku:., Orlando,
Florida, 1986.

Kron, G., Diakoptics, The Piecewise Solution of Large-scale Systems,
Macdonald, London, 1963.

Lorenz, F., and Wolper,;Lq "Assigning Causality in the Case of’
.Algebraic loops," J. of the Franclin Inst., Vol. 319, No. 1/2, pp. 237-
243, 1985.

Margolis, D. L. , and Tabrizi, M., "Reduction of Models of Interacting
Lumped and Distributed System Using Bond Graphs and Repeated Model
Decomposition," J. Franklin.Institute, vol. 317, No. 5, pp. 309-322,
1984.

Matrosov, V. M., "On the Theory of Stability of Motion (in Russian),"
Priklandaia Mathmatica i Mekhanika, 26, pp. 992-1002, 1962.

Meirovitch, L., Analytical Methods in Vibrations, The MacMillan Co., New
Yourk, NY, 1967.

Michel, A. N. , and Miller, R. R., Qualitative Analysis of Large Scale
Systems, Academic Press, New York, 1977.

Ort, J. 11., and Martens, H. R., "The properties of Bond Graph Junction
Structure Matrices," ASME, Journal of Dynamic Systems, Measuremernz, and
Control, Vol. 95, No. 4, pp. 362-367, December 1973.

147

Paynter, H. M., "The Rational Foundations of Multiport Dynamics," July
27, 1967 (unpublished notes)

Perelson, A. S, "Bond Graph Junction Structures," ASME, Journal of
Dynamic Systems, Measurement, and Control, Vol. 97, No. 2, pp. 189-195,
June 1975.

Pernebo, I.., and 1“ Di. Silverman, "Model Reduction via Balanced State
Space Representations," IEEE Trans. Automatic Control, AC-27, No. 2, pp.
382-387, 1982.

Prakash, R., Palmer, E., and Withers, D., "Design assessment: Electro-
Mechanical Systems," Research Report RC-10993, IBM Thomas J. Watson
Research Center, Yorktown Heights, NY, 1985.

Rosenberg, R. C. , "State-Space Formulation for Bond Graph Models of
Multiport Systems," ASME Journal of Dynamic Systems, Measurement, and
Control, Vol. 93 N0. 1, pp. 35-40, March 1971.

Rosenberg, R. C. and Karnopp, D., "A Definition of the Bond Graph
Language," ASME, Journal of Dynamic Systems, Measurement, and Control,
Vol. 94, No. 3, pp. 179-182, Sept. 1972.

Rosenberg, R. C. , "On Gyrobondgraphs and Their Uses," ASME, Journal of
Dynamic Systems, Measurement, and Control, Vol. 100, No. 1, PM). 76-82,
March 1978.

Rosenberg, R. C., "Essential Gyrators and Reciprocity in Junction
Structure," Journal of Franklin Institute, Vol. 308, No. 3, pp. 343-352,
September 1979.

Rosenberg, R. C., and A. N. Andry,Jr.. "Solvability of bond Graph
Junction Structures With Loops," IEEE, Transaction on Circuits and
Systems, Vol. GAS-26, No 2, pp. 130-137, February 1979.

Rosenberg, R. C., and Moultrie, B., "Basis Order for Bond Graph Junction
Structures,’ TIEEE Transaction on Circuits and Systems, GAS-27, No. 10,
pp. 909-920, Oct. 1980.

Rosenberg, R. C., and Karnnop, D. C., Introduction to Physical System
Dynamics, McGraw-Hill, New York, 1983.

Rosenberg, R. C., The ENPORT-5 User's Manual, ROSENCODE Associamai
Lansing, MI, 1984.

Rosenberg, R. C. , "Improving Design Assessment Techniques for
Reliability and Yield by Bond Graph Technology" , Proposed to 18.1
Manufacturing Research, Department oftflufiunncal Engineering, Michigan
State University, 1985.

Rosenberg, R. C. , The ENPORT-7 Reference Manual, ROSENCODE Associates,
Lansing, MI, 1986. -

Sage, A. P., Methodology for Large Scale Systems, McGraw-Hill,Pkmz
York, 1978.

148

Sandell Jr., N. R., Varaiya, P., Athans, M., and Safonov, M. G., "Survey
of Decentralized Control Methods for Large Scale Systems," IEEE Trans.
on Automatic Control, Vol. AC-23, No. 2, April 1978.

Siljak, D. D., Large-scale Dynamic Systems, North-Holland, New York,
1978.

Siljak, D. D., "Complex Dynamic Systems: Dimensionality, Structure, and
Uncertainty," Large Scale Systems 4 (1983), pp. 279-294.

Skelton, R. E. , "Cost Decomposition of Linear Systems with Application
to Model Reduction," Int. J. Control, Vol. 32, No. 6, pp. 1031-1055,
1980. ~

Skelton, R. E., and Yousuff, A., "Component Cost Analysis of Large Scale
Systems," Int. J. Control, Vol. 37, No.2, pp. 285-304, 1983.

Syslo, M. M., Deo, N., and Kowalik, J. 8., Discrete Optimization
Algorithms, Prantice Hall, Englewood Cliffs, NJ, 1983, pp. 320.

Titli, A., and Singh, M. G., Large Scale Systems: Theory and
Applications, IFAC, Pergamon Press, Oxford, 1980.

Thompson, W. T., Theory of Vibrations with Applications, Prentice-ikrll,
Inc., Englewood Cliffs, NJ, 1981.

Tse, E. C. Y., Medanic, J. V., and Perkins, W. R., "Chained Aggregation
of Linear Time Invariant Systems," Proc. JACC. San Fransisco, CA. 1977.

Yousuff, A., and Wagie, D. A., and Skelton, R. E., "Linear System
Approximation via Covariance Equivalent Realization," J. Math. Analysis
and Appl., 106, pp. 91-115, 1985.

Zalewski, Z. , and Rosenberg, R. C., "Simulation of Engineering Models
Containing Bond Graphs and Block Diagrams," ASME Computer in Engineering
Conference, 355-359, Chicago, IL, July 1986.

Zeid, D. A., "Bounds on the Eigenvalues for Certain Classes of Dynamic
Systems," Ph. D Dissertation, Michigan State University, 1982.

Zeid, D. A., "Solution of a set of Nonlinear Algebraic Equations for
General-Purpose CAD Program," IEEE Circuits and Devices Magazine, pp. 7-
20, Sep. 1985.

Zhou, T., "A Parallel Computation Method for Dynamic Systems with
Coupled nonlinear Dissipation," MS Thesis, Michigan State University,
1984.