£.
a. 2.
.9L.5~.:...§v . .p.
23.!fifl has
.n a“. 1 Jun.
faunas.

‘ .. «can.
. J.

r.
Rafi?» ‘
has": ,

. . it}.

’,‘.o
'1‘? 3’
F. . n a
. 4
, huh
{‘I I
Iii-1i
sum;

3&qu

Z 2 a.

:5! 37.? I!

it. A. 99.!!9
x in

inst: .. am?
19:3. "1.1%“ r ‘
in...

010.“!!! 1.3..
z]...
3%
:I.I..4..!l'
. a I?” u
59:03.. J
4‘9):

5.33. c ,
v 5‘“: L
abilitilfsi
2.9.1:}: .. .

V .1525:
29” 103%.:
2!. .1
I.
ablaztil....:.
I ....\..!.).l::tia
.:.mq.,..i....Mii..§..2...i
lilac-‘3‘ if}?.!..b.aa1
..:|‘.e’..¢t.io.9£¢.onitflls.
all ahflilta).»v.
.3 f.
tufts...

5. 3!!
.. .a
I"?! ”out.

:LHL.
3155

95 . £352}!
t...N¢l§?,€I!..
.. ‘

a I! .nflsaliii.lu::v1§! Ir

 

3.§i

 

 

LIBRARY

Michigan State
Unlversity

 

 

 

This is to certify that the

dissertation entitled

ANALYSIS AND DESIGN OF CLASSES 0F TWO-TIME-SCALE
MULTIPORT MODEIS WITH DERIVATIVE CAUSALITY

presented by

A hmed A bdelaziz Omara

has been accepted towards fulfillment
of the requirements for

. Mechanical Engineering
degree in

flC/rwm

 

 

Major prof or

(JP/2% a
/ T

MSU is an Afflrman’w Action/Equal Opportunity Institution 0-12771

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

11m macs-p.14

ANALYSIS AND DESIGN OF CLASSES OF TWO-TIME-SCALE
MULTIPORT MODELS WITH DERIVATIVE CAUSALITY

By

Ahmed Abdelaziz Omara

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSPHY
Department of Mechanical Engineering

2000

ABSTRACT

ANALYSIS AND DESIGN OF CLASSES OF TWO-TlME-SCALE
MULTIPORT MODELS WITH DERIVATIVE CAUSALITY

By

Ahmed Abdelaziz Omara

Eliminating derivative causality by inserting parasitic elements produces modified
systems that are two-time-scale. The computational accuracy and simulation
efficiency of the solution of the modified systems are strongly influenced by the
choice of the inserted parasitic element parameters. Computational accuracy
can be traded off with simulation efficiency. Some important aspects of the
nature and behavior of the modified two-time-scale systems have not been
investigated. In addition, there is no known method that shows how to choose
the values of the parasitic parameters to achieve certain desired accuracy at a
highly efficient computational price. Tools based on the singular perturbation
method, on physical systems modeling, and on numerical simulation are used in
this research to analyze and design some classes of these two-time-scale

modified systems.

Analysis revealed some important aspects of the nature and behavior of
the systems under investigation. Understanding those aspects is extremely
important for the general understanding of these systems. A choice for the
perturbation parameter a that helps to reveal the time-scale properties of the

modified system is given; a suggested preferred standard singular perturbation

form (PSSPF) is shown; a systematic method to transform the modified systems
to the PSSPF is derived; models of the fast sub-systems are obtained and

validated; and proof of stability of such a modified system is shown.

The results of the analysis efforts are used to set procedures to design the
fast sub-systems. The design objective is to achieve certain desired accuracy at
highly efficient computational efforts. The validated models of the fast sub-
systems are used in the design procedures. Examples that demonstrate the
usefulness of these procedures are presented. Analysis and design work shown

here include classes from both linear and nonlinear systems.

To my parents

ACKNOWLEDGEMENTS

For my parents, completing this degree comes in partial return on the
many investments they have made down a long road; and in partial fulfillment of
some of their dreams. Their relentless support and dreams contributed a major

thrust to my drive throughout my life.

I would like to acknowledge the continuous support I received from Dr.
Ron Rosenberg, the chairman of my Ph. D. committee, since I joined MSU. In
particular, I indeed appreciate the quality time and thoughts that he regularly
shared with me. It is a fact that my main motive to join MSU was joining Dr.
Rosenberg’s program. The other important and necessary role for shaping and
guiding these research efforts was Dr. Hassan Khalil’s. I am indebted to him for
the help he has given to these research efforts. He has never hesitated to assign
a slot from his time to strengthen these efforts. If there were a co-chairman post
for my Ph. D. committee he would have earned it without any doubt. Between
him and Dr. Rosenberg, I enjoyed a very precious access to some unique blend
of expertise. I would also like to thank my other committee members, Dr.
Mukherjee, Dr. Newhouse, Dr. Rohde, and Dr. Radcliffe for their important role in

guiding and sharpening these research efforts.

Last but not least, i would like to thank my wife, Manal, for her patience

and understanding during some extremely busy times.

TABLE OF CONTENTS

CHAPTER 1
INTRODUCTION .................................................................................................. 1
1.1 Mathematical Modeling of Mechatronic Systems ........................................ 1
1.2 Multiport Modeling of Mechatronic Systems ............................................... 2
1.3 Some Equation Formulation Issues in Mechatronic Modeling .................... 7
1.4 Problem Statement and Research Objectives ............................................ 9
1.5 Dissertation Outline ................................................................................... 11
CHAPTER 2
ANALYSIS OF A CLASS OF TWO-TlME-SCALE MULTIPORT MODELS ........ 13
2.1 Singular Perturbation Method and Two-Time-Scale Systems ................... 13
2.1.1 Some Important Issues Related to Using the SP Method .................. 13
2.1.2 Additional Remarks on the Standard SP Form ................................... 16
2.2 Classes of Bond Graph Models Under Investigation ................................ 17
2.3 Related Work ............................................................................................ 20
2.3.1 The Algorithm Under Investigation ..................................................... 20
2.3.2 Other Related Work ............................................................................ 23
2.4 Determining the Appropriate Perturbation Parameter ............................... 25
2.4.1 Physical-System-Modeling-Based Analysis ....................................... 26
2.4.2 Simulation-Based Analysis ................................................................. 28
2.5 Finding the Standard Singular Perturbation Form ..................................... 37
2.5.1 Using Transformation Procedures to Obtain PSSPF .......................... 38
2.5.2 A Transformation to Obtain the PSSPF .............................................. 50
2.5.2.1 Derivation and Proof of Existence ................................................ 53
2.5.2.2 Example of Using the Transformation Matrix ............................... 55
CHAPTER 3
DESIGN OF A CLASS OF TWO-TIME-SCALE MULTIPORT MODELS ............ 61
3.1 Introduction ............................................................................................... 61
3.2 A Model for the Fast Sub-System and Proof of Stability ........................... 61
3.3 Validating the Fast Sub-System Model ..................................................... 66
3.3.1 Accuracy of Slow and Fast Poles Locations .................................. I ..... 67
3.3.2 Conditions of Two Critically Damped Fast Poles ................................ 75
3.3.2.1 Advantages of Critically Damped Fast Poles ............................... 79
3.3.3 Effect of Selected Dependent Element Value on Model Performance 81
3.4 An Algorithm to Improve the Model Performance ..................................... 84
3.5 Design Procedures for the Fast Sub-System ............................................ 88
CHAPTER 4
ANALYSIS AND DESIGN OF A CLASS OF TWO-TIME-SCALE NONLINER
MULTIPORT MODELS ....................................................................................... 96
4.1 Introduction ............................................................................................... 96

vi

4.2 Class of Nonlinear Systems Under Investigation ...................................... 96

4.3 Obtaining the PSSPF ................................................................................ 98
4.4. A Model for the Fast Sub-System and Proof of Stability ........................ 102
4.5 Effect of Inaccurate Effective Field Value on Model Performance .......... 106
4.6 Design Procedures for the Fast Sub-System .......................................... 111
4.7 An Example ............................................................................................ 114
CHAPTER 5
CONCLUSIONS ............................................................................................... 123
5.1 Summary of Main Contributions .............................................................. 123
5.2 Suggested Directions for Future Research ............................................. 126
5.2.1 More than One Dependent Element ................................................. 126
5.2.2 Analysis and Design for Eliminating Algebraic Loops ....................... 127
REFERENCES ................................................................................................. 128

vii

CHAPTER 1
INTRODUCTION

1.1 Mathematical Modeling of Mechatronic Systems

Mathematical modeling is becoming more and more an important step in
design, analysis and development of engineering systems. See, for example,
Chang and Rohde, 1996. There are two main reasons for this phenomenon. (1)
Global competition produced a need to reduce time-cycle and related cost of
product design and development. (2) The continuous improvement of computer
support, in terms of both the cost and the capability, makes CAE more effective.
Therefore, mathematical modeling is increasingly perceived in the industry as an
investment that is necessary to survive in a competitive environment. There are
many research and development efforts focused on developing new tools to aid
engineers in the abstraction, mathematical formulation, and analysis processes

of product design.

Engineering systems have become increasingly complex and inter-
disciplinary. Typically, a variety of components are linked together, including
control sub-systems and imbedded "intelligence" to form a system. In 1967
Yasakawa Electric America (a Japanese company) coined the term
“mechatronics” to refer to systems that combine mechanisms, electronics, and
control (T omkinson et al. 1996). The term was coined by combining the words

”mechanics" and "electronics" (Auslander et al. 1996). Currently there are

several common definitions for the term, with perhaps more underway. However,
there is a common understanding that continues to be associated the term. It is
the interaction among different physical sub-systems such as mechanical,
electrical, magnetic, hydraulic, pneumatic, and thermal, with electronics, control,
and imbedded intelligence. For such systems, the multiport approach is seen as
a suitable modeling environment (Cellier, 1992; Cellier 1995; Kamopp et al.,

2000).

1.2 Multiport Modeling of Mechatronic Systems

Physical sub-systems with one or more ports are called multiports
(Karnopp et al. 2000). Ports are places at which power can flow among the sub-
systems. As indicated above, mechatronic systems include signal coupling in
addition to the power coupling. Bond graph is a modeling environment that
consists of elements/sub-systems linked together by power bonds. A bond graph
shows power coupling and describes the exchange of energy in a system. Each
bond in a bond graph, in general, implies the existence of a pair of signals
flowing in opposite directions. The two signals are the flow signal and the effort
signal. As multiports essentially transmit finite power when interconnected,
which requires the existence of a flow variable and an effort variable, bond graph
is a very efficient way to model multiports. Figure 1.1 shows a simple hydraulic

system and its bond graph model (from Rosenberg et al. 1983). On the other

hand, block diagrams are efficient in modeling signal coupling. Therefore, block

diagram can be combined with bond graph to model mechatronic systems.

 

 

 

 

 

Q.

P2 I .__I

 

 

 

 

 

 

 

Figure 1.1. A hydraulic cylinder and its bond graph model

The main reasons why the bond graph is valuable as a modeling environment for

mechatronic systems are

1) It provides a unified graphical representation for multiport modeling
regardless of the physical domain type. This makes the bond graph
especially valuable for systems containing transducers.

2) It provides tools to perform graphical analysis that provide insight to
the user efficiently.

3) It enables algorithmic equation formulation that is suitable for computer
implementation. This leads to an efficient and error-free formulation of
system equations for large and complex models (Rosenberg et al.,

1983; Cellier, 1992; Cellier 1995; Kamopp et al., 2000).

The literature includes many examples of using bond graphs to gain
insight into the dynamics of systems and to reveal some of their structural
properties. To mention a few examples, Sueur and Dauphin-Tanguy (1989,
1991) present necessary and sufficient conditions for structural
controllability/obsewability of linear systems represented by bond graphs.
Rosenberg and Zhou (1988) proposed a power-based tool for visualizing the
dynamic response of engineering systems in bond graph form. They showed
that the power response of dynamic systems could be used as a basis for model
order reduction under some conditions. In particular, they used the RMS power
levels in a bond graph model to determine the weakest power connections in the
model. These bonds became candidates for elimination to reduce the model
order. The tool is suitable for handling both linear and nonlinear systems.

Rosenberg (1987) discusses several applications of causality in understanding

the dynamic structure of models, such as revealing the presence of algebraic
couplings. Zeid and Rosenberg (1986) show that, by transforming a bond graph
model of a linear time invariant system to a canonical form, the gyrobondgraph,
inspection of the resulting graph structure reveals various properties of the
system eigenvalue spectrum. Dauphin-Tanguy and Borne (1985) show that the
fast and slow dynamics of bond graph models can be estimated by determination
of causal loop-gains. They used this result in model order reduction of multi-
time—scale systems. The powerful analysis tools, great insight, and the
computer-aided enabling environment offered by the bond graph give it a
decisive advantage as an engineering modeling environment. Consequently,
numerous research efforts are devoted to addressing issues that arise in the
bond graph environment. This effort is among them; it is dedicated to addressing

the derivative causality issue.

The literature also provides many examples of algorithmic formulation of
system equations from bond graph models. A variety of algorithms to obtain
formulations for different physical domains have been proposed (Karnopp et al.,
1990). That includes obtaining the preferred formulation for certain classes of
problems such as the Lagrange Formulation for non-linear mechanics (Karnopp
1977, 1978). One can proceed algorithmically from a bond graph model to a set
of fundamental equations or to different sets of secondary equations. A set of
equations representing a system is called a fundamental set if it preserves all of

the structural information about the system. Otherwise, the set of equations is

called a secondary set. Structural information depends only on the system’s
topology and the nature of its components. Structural information does not
depend on the numerical values of the parameters of the system elements. As
an example, consider an electrical circuit. A set of fundamental equations is
obtained by applying the volt-current relation for every element in the circuit,
Kirchhoff’s current law, and Kirchhoff’s voltage law. The state equations, the
nodal equations, and the mesh equations are three different sets of secondary
equations. While it is usually straightforward to proceed from a set of
fundamental equations to different sets of secondary equations, the case is
different if you try to proceed from a set of secondary equations to another set of
secondary equations. Karnopp (1983) points out that the bond graph is unique in
providing a mechanism for demonstrating causal interactions among different
components and sub-models; therefore, studying alternative equation
formulations becomes easier. He describes how different causal patterns can
lead to different equation formulations for the same dynamic system. Because
the procedures offered in the bond graph environment are algorithmic, they can
be computerized. ENPORT is one example of a software package that
implements bond graph procedures (Rosenberg and Karnopp 1983). The author
views the bond graph as a basic physical model representation that can lead to a

variety of equation formulations and analysis methods.

1.3 Some Equation Formulation Issues in Mechatronic Modeling

A typical equation formulation goal is to obtain the standard explicit state
space form

%’f—=f(x,u,o (1.1)

Where, X is the vector of independent state variables, I is time, and U is the
input vector. However, systems could have two different types of couplings
effects that may prevent this form from being attainable. The first is the so-called

“algebraic coupling.” The resulting equations have the general (functional) form

£=f(X,H,U,t) (1.2)
dt
0=g(X,H,U,t) (1.3)

where H is the coupling vector. The second type of coupling that may prevent
the realization of the explicit state space form is “energy-variable coupling.” The

resulting equations have the general (functional) form

dX. dX
—_L : X'a d aUat 1.4
dr f( "7 > ( )
X4 =(—0-(Xi,U,t) (1'5)
Where, X, is the vector of independent state variables, X d is the vector of

dependent state variables. Of course it is possible for both types of couplings to
exist in the same model. The energy-variable coupling is known as ”derivative

causality“ in bond graph terms (Rosenberg et al. 1983). This research effort is

related to alleviating the effect of the derivative causality for some classes of
dynamic systems. A definition for these classes will be provided later; equations

(1.4) and (1.5) define only a functional form.

Equations (1.2) and (1.3) and equations (1.4) and (1.5) are two implicit
sets of differential equations. They can be converted to explicit sets by algebraic
manipulation if the system is linear. If the number of coupled elements is greater
than two or three, the algebra may become cumbersome. If the system is
nonlinear the algebraic manipulation might be very difficult or even impossible
(Rosenberg et al., 1983; Karnopp et al., 2000). When one applies an integration
formula to a set of differential equations, the result will be an explicit set of
difference equations only if both the integration formula and the set of differential
equations are explicit. Othenlvise, the result will be an implicit set of difference
equations. Therefore, integrating an implicit set of equations such as the set
given by equations (1.4) and (1.5), is an implicit-integration problem, regardless
of the type of the integration formula used. Implicit integrators are generally
computationally more involved than explicit integrators (Bennett, 1995).
However, the computational cost is not the major reason why we usually do not
prefer the presence of derivative causality. The following section includes a

summary on why we do not prefer the presence of derivative causality.

1.4 Problem Statement and Research Objectives

As stated above, bond graph models with derivative causality produce
mathematical models with algebraically coupled derivative terms. Some times it
is preferable to eliminate the derivative causality, and consequently to enable an
algorithmic derivation of an explicit ordinary differential equations (ODEs) model.
The main reasons for preferring explicit ODEs are:

1) to take advantage of many control algorithms that are available for
systems expressed in the explicit ODEs form (Turner & Zeid, 93;
Krishnaswami, 93);

2) to take advantage of the standard CAD control packages that can
algorithmically linearize, analyze, and synthesize controllers for dynamic
systems (Turner & Zeid, 93; Krishnaswami, 93);

3) to take advantage of easy and error-free automated formulation of
dynamic multiport equations (Rosenberg, 71; Karnopp, 79; Zeid, 88);

4) to obtain a more parallelizable set of equations, for computation on
parallel processors (Turner & Zeid, 93; Krishnaswami, 93).

In addition, this work shows that the derived set of explicit ODEs might require
less computational effort to integrate, compared with the original implicit set,
even if the explicit state-space order is increased as compared to the original

model.

There are two common ways to eliminate the effect of the derivative
causality. The first is to obtain an equivalent expression for the behavior of the
variables. This approach guarantees that all responses’ behavior match exactly
those of the original model; except for the internal details of the coupled sub-
system. This approach needs some algebraic manipulation that might be difficult
or impossible if the system is nonlinear. For these cases where the algebraic
manipulation is not preferable/possible the second way becomes an attractive
option. In the second way, some parasitic elements are inserted to break the

coupling, and consequently eliminate the derivative causality.

Karnopp and Margolis (1979) showed a method to eliminate derivative
causality in bond graph models. It breaks the coupling between the coupled
energy storing elements by inserting parasitic elements into the original model.
Consequently, a modified model that approximates the original model is
produced. This modified model is supposed to be a two-time-scale (stiff) model
so that the dynamic behavior of the original state variables would not be greatly
affected (Karnopp et al., 1979; Zeid, 1988; Zeid, 1989). This approach is the
approach being investigated in this research work. The choices of the inserted
parasitic elements’ values directly affect both computational efficiency and
approximation accuracy; approximation accuracy can be traded off with the
system stiffness. In general, there are some important aspects of the nature and
the behavior the modified two-time-scale system that have not been investigated

in detail. In addition, there is no known method or algorithm to show the user

10

how to choose the values of the inserted parasitic parameters to achieve high

computational efficiency while satisfying certain desired accuracy.

This research work has two main objectives:

1) Analyze the class of two-time-scale modified systems that results from
using the above-mentioned approach to eliminate the derivative causality.
The purpose of this analysis is to better understand the important two-
time-scale properties of that class of two-time-scale systems.

2) Forrnalize the above-mentioned approach of using parasitic elements to
eliminate the derivative causality in a design context. The purpose of this
formalization is to enable the user to achieve the greatest possible
simulation efficiency for specified approximation accuracy, by making the
appropriate choices for the values of the parasitic elements parameters.

A description of the classes of dynamic systems for which this research effort will

be conducted will be defined later.

1.5 Dissertation Outline

CHAPTER 2 starts by briefly discussing the usefulness of the singular
perturbation method for analysis and design of two-time-scale systems. Then,
the classes of bond graph models for which these research efforts will be
conducted are defined. This is followed by a brief review of related work with

more emphasis on the bond-graph-based methods for eliminating derivative

11

causality. The last two sections of CHAPTER 2 show how to determine the
appropriate perturbation parameter a for the classes of modified systems under
investigation; and how to use this perturbation parameter to transform them to
the standard singular perturbation form. A preferred standard singular
perturbation form (PSSPF) is shown and the transformation needed to obtain it
directly from the modified systems equations is derived. CHAPTER 3 carries on
the results of CHAPTER 2, and shows some procedures to design the fast sub-
system. First a model for the fast sub-system is derived and validated; and the
stability of the modified system is demonstrated. An algorithm to improve the
model performance is shown, before the design procedures are given. An
example demonstrating the usefulness of the design procedures is included.
Nonlinear systems are studied in CHAPTER 4. Therein, a class of nonlinear
systems is defined. Then, obtaining the PSSPF is shown and a model of the fast
sub-system is derived. Stability of the modified system is also proved. The
chapter concludes by giving some procedures to design the fast sub-system of
the class of nonlinear systems under investigation, with an example
demonstrating the usefulness of the design procedures. CHAPTER 5 gives a
summary of the main contributions and suggests some future research

directions.

12

CHAPTER 2
ANALYSIS OF A CLASS OF TWO-TIME-SCALE MULTIPORT MODELS

2.1 Singular Perturbation Method and Two-Time-Scale Systems

It was stated in CHAPTER 1 that the algorithm, for eliminating the
derivative causality, under investigation breaks the coupling between the coupled
energy storing elements by inserting parasitic elements into the original model. It
was also stated that a modified model that approximates the original model is
produced; and that this modified model is supposed to be a two-time-scale (stiff)
model so that the dynamic behavior of the original state variables would not be
greatly affected (Karnopp et al., 1979; Zeid, 1988; Zeid, 1989). The singular
perturbation method has proved to be very useful in analysis and design of multi-
time-scale dynamic systems (Kokotovic et al., 1976; Khalil, 1986; Kokotovic et
al., 1986). It decomposes a two-time-scale system into a slow and a fast sub-
systems by assuming a small singular perturbation parameter approaching zero.
Tools to analyze the slow and the fast sub-models and the interrelation between
them are available. For these reasons the singular perturbation (SP) method will

be used (among other tools) in this research work.

13

2.1.1 Some Important Issues Related to Using the SP Method

There are two important issues related to the use of SP method. These
two issues can be defined as (1) the determination of the “appropriate”
perturbation parameter (s). (2) Transforming a singularly perturbed model to
the standard singular perturbation form. Solving the first issue is a pre-request
for solving the second one. This sub-section will briefly discuss these issues;
later in this chapter, they will be dealt with so that the singular perturbation

method can be of useful use to the problem under investigation.

Although the singular perturbation theory has seen significant
developments, and has many applications; theory users have a major concern.
The time-scale separation is not always, easily detectable. In many situations, a
good understanding of the physical process Is needed to detect the separation,
and to choose the perturbation parameter(s) (Kokotovic et al., 1986) and (Sueur,
1991). Choosing the “appropriate” perturbation parameter is usually crucial to
the usefulness of singular perturbation analysis. This task becomes more difficult
when the system involved is complex. In addition, the perturbation parameter of
a two-time-scale system does not have to be unique; there might be several

possible perturbation parameters.

The general form of the standard singular perturbation model is

X:f(X,Z,£,t), X(t X0, XER", (2.1)

0)=

l4

£2 = g(X,Z,€,t), 200) = 2°, 26 R’" (2.2)

where; e is a small positive scalar; a superscript dot denotes a time derivative;

and f and g are functions that are assumed to be sufficiently several times

continuously differentiable. This model is the most studied SP model in the
literature and the first one to be used in control and systems theory (Kokotovic et
al., 1986). Consequently, numerous analysis and design tools are available for
two-time-scale systems that are in the standard form (2.1) and (2.2). Of course,
expressing the model in the standard form is not going to change any thing in the
system’s time-scale properties; but it is going to reveal them. In contrast, a
singularly perturbed model that is not in the standard form might hide its time-
scale properties. Even though, there are procedures to help transform a two-
time-scale system from a non-standard form to the standard form; these

procedures are general and their outcome is not unique.

Zeid (1988 and 1989) suggested that the modified systems that result
from removing derivative causality by inserting parasitic elements could be
transformed to a standard singular perturbation form. He also noticed that once
a standard form is reached the singular perturbation theory tools could be
conveniently used. However, there was no method defined to show how to
transform the system and consequently, there was no attempt made to use the

singular perturbation analysis tools.

15

The objectives of the work reported in this chapter are

(1) To introduce a choice for the singular perturbation parameter (a) that
helps reveal the two-time-scale properties of the class of modified
systems under investigation.

(2) To show how the general procedure outlined in Kokotovic et al. (1986) for
transforming singularly perturbed models to the standard form can be
used to obtain a preferred standard singular perturbation form (PSSPF)
and to gain insight into the problem under investigation.

(3) To define the transformation needed to obtain that PSSPF.

2.1.2 Additional Remarks on the Standard SP Form

When a is set to zero, equation (2.2) reduces/degenerate to the algebraic
equafion
o = g(X‘,Z,o,r) (2.3)
The bar in )7 and 7 indicates that the variables belong to a system with a = 0.
Model (2.1) and (2.2) will be said to be in the standard form iff (2.3) has k 21
distinct real roots

2 = (Trix), i: 1,2,...,k. (2.4)

This assumption ensures that each root of (2.4) will correspond to a well-defined
n-dimensional reduced model that can be obtained by substituting (2.4) into (2.1).

The ith reduced model would be

16

)7 = f(X‘,¢Ti(Y,z),o,r), in = 0) = X(t = 0) (2.5)

Notice that the initial conditions of X— are the same as those of X . The model
(2.5) is called a quasi steady state model because when a is small Z will
converge rapidly to the solution of (2.3), (Kokotovic et al., 1986). It has been
shown in the literature (Kokotovic et al., 1986; Khalil, 1986; Kokotovic et al.,

1976) that, for a small 5 and a stable system
X(t) — in) = 0(5), r 2 0 (2.6)
and, in addition,

2(r) - in) = 0(5), t > 0 (2.7)
2.2 Classes of Bond Graph Models Under Investigation

The classes of bond graph models under investigation in this research

effort are the classes that can be represented using the basic bond graph set
S1={g,g,§,S,,Sf,TF,GY,O,1}. The basic set S1 includes multiport capacitances,

multiport inertances, multiport dissipations, sources of effort, sources of flow,
displacement-modulated transformers, displacement-modulated gyrators, 0-
junctions, and 1-junctions. The general homogeneous mathematical form for the
derivative causality that occurs in these classes of bond-graph models can be
expressed in equations (2.8) to (2.11). This form assumes that the model does

not include algebraic loops that are not separable from the derivative causality.

l7

Xi:Jrr(Xr)*Zr+er(Xr)*Xd (2-8)

2, = J,,(X,)*Z, (2.9)
2, = (p,(X,) (2.10)
X, =(p,(z,) (2.11)

Where, the J matrices are the constraints imposed by the junction structure

between sets of ports; they are in general nonlinear. X , is the independent
energy variables vector, X, is the dependent energy variables vector, Z, and Z,
are the coenergy vectors for X, and X, respectively, U is the sources input
vector; and (a, and (p2 are generally nonlinear functions, Rosenberg (1971), van
Dijk (1994). Throughout this document the dimension of X, will be one. The

above general formulation allows non-linearity in both the fields and junction
structure. If the fields and junction structure are both linear, the fields’ relations
and the junction structure matrices will become linear. Consequently, the

junction matrices will no longer be function of the independent state vector X ,.

The general mathematical form of these classes of linear systems can be written

as follows.
X, =1, *2, +1, *X, (2.12)
z, = 1,, *z, (2.13)
z, = s, * X, (2.14)
X, = s,“ *z, (2.15)

18

Rosenberg (1971) has also shown that 1,, and 1,, are one another’s negative

transpose. If 2, and Z, are eliminated, the linear general mathematical form will

simplify to
X,=J,,*S,*X,+J,,*X, (2.16)
S,*X,=J,,*S,*X, (2.17)

Zeid’s proposed algorithm is supposed to be implemented in the bond graph
environment. Actually, facilitating the automated derivation of an explicit set of
equations is one of the main reasons for using his algorithm. However, since this
work includes analysis and design of the modified systems that are produced
from this algorithm, it is necessary to show the general mathematical form for
these modified systems. If the original system is in the form (2.16) and (2.17) the
homogeneous modified system will be in the general form (2.18) to (2.20).

X, =J,, *S, *X, +J,, *(K, *q+R*(J,, *S, *X, —S, *X,)) (2.18)

X,=K,*q+R*(J,,*S,*X,—S,*X,) (2.19)

q = Jdi *5) :1: X. -54 :1: X. (2.20)
Where, K p and R are the inserted parasitic stiffness and resistance

respectively; q is the flow variable in the inserted parasitic elements. In these

research efforts we will study systems with one derivative causality per storage
field. Of course using the algorithm is not limited to linear systems only. First
linear systems that can be expressed in the general form (2.16) and (2.17) and
whose modified systems are in the general explicit form (2.18) to (2.20) will be

studied. Then, in CHAPTER 4, a class of nonlinear systems will be studied.

19

2.3 Related Work

Because the method of eliminating the derivative causality under
investigation here is bond-graph-based, this section will discuss only the bond-
graph-based methods to eliminate the derivative causality. Examples of the
other methods include the kinematics-based work of Allen (1981) to obtain
explicit ODEs; the work of Gordon et al. (1998) to convert systems of ADEs to
ODEs; and the work of Borutzky and Cellier (1996a and 1996b) which partition
(tear) the resultant ADEs from the bond graph using symbolic manipulation for a
more efficient solution. Presentations of the bond-graph-based methods will be

very brief except for the method that is under investigation in this research work.

2.3.1 The Algorithm Under Investigation

This research work is a continuation of the method for eliminating
derivative causality demonstrated by Karnopp-Margolis (1979) and generalized
and automated by Zeid (1988). In 1979 Karnopp and Margolis demonstrated a
method to eliminate the derivative causality in bond graph models of planar
mechanisms. The method eliminates the derivative causality by inserting
parasitic damping and parasitic stiffness. Zeid (1988) showed an algorithm to
generalize and automate the method demonstrated by Karnopp and Margolis.

Zeid summarized his proposed algorithm as shown in Figure 2.1 to Figure 2.4.

20

I

I . .le
/\ / I

Figure 2.1. Handling an inertance that is attached to a 1-junction
I
0

i
\/°\.:i/ \T’ I

Figure 2.2. Handling an inertance that is attached to a O-junction

I

21

C R
\/0/\ l:\>/ 0 ——’I 1
Figure 2.3. Handling a capacitance that is attached to a O-junction
C
_L

/‘ KT /‘1/\_"

Figure 2.4. Handling a capacitance that is attached to a 1-junction

_.I
I

_ I‘E— Q___’l 3

As the figures above indicate, a dependent inertia is treated by adding a parasitic
damping and a parasitic stiffness; the parasitic elements are in series with
respect to one another, and they are both in parallel with respect to the

dependent inertia. On the other hand a dependent stiffness is treated by adding

22

a parasitic damping and a parasitic inertia; the parasitic elements are in parallel
with respect to one another, and they are both in series with respect to the
dependent inertia. The dependent stiffness problem and the dependent inertia
problem are the dual of one another. Therefore, it is sufficient to analyze only

one of them. It has been decided to do the dependent inertia problem.

2.3.2 Other Related Work

Joseph et al. (1974) presents a method for causality assignment (relaxed
causality) for nonlinear systems with algebraic loops and/or derivative causality.
The method will still produce a set of ADEs for systems with algebraic loops
and/or derivative causality. However, he claims that it is more suitable for
computerized formulation of nonlinear systems than the widely used Sequential

Causality Assignment Procedures (due to Karnopp and Rosenberg).

Granda (1984) shows different alternatives to deal with derivative causality
and algebraic loops. The proposed alternatives can be summarized as being
deleting elements from the original model, adding elements to the original model,
or algebraic manipulation. Each of these alternatives could be the right choice,
depending on the problem under investigation and depending on the objective of
the process. He emphasized the importance of alternatives that can be

implemented on a computer.

23

Bos et al. (1985) describes a method that uses formula manipulation to
remove differential causality from mechanical systems. It states that the

computational demand of the formula manipulation is “excessively large.”

Karnopp (1992) proposed another approach to eliminate derivative
causality in bond graph models of mechanical systems. The approach defines
an l-field or an IC-field using generalized momenta or generalized coordinates
similar to those used in Lagrange’s equations. Karnopp (1992) showed that this
method would require a matrix inversion, at each time step if nonlinear
constraints such as modulated transformers were involved. In addition, this
method would be more difficult to automate compared with the method under

investigation in this research.

Margolis (1993) acknowledges that when the goal of analysis is to
simulate the system's transient response, the derivative causality issue has to be
dealt with; and shows a way to derive input/output transfer functions for bond

graph models with derivative causality.

Edstrom (1999) shows that by inserting a parasitic resistance into a class
of bond graph models with derivative causality; then deriving the explicit set of
equations for the modified system; and setting the parasitic resistance to zero (or
infinite); a set of explicit ODEs can be obtained for the original system that

contained a derivative causality. Results were proved for a very limited class of

24

linear bond graph models that did not include transformers, gyrators, or multi-port
elements. The main importance of this paper is highlighting the relationship
between mathematical models derived from bond graphs with derivative causality
and the standard singular perturbation model; although it was addressing a less
general (compared to the modifications under investigation in this work) type of

modifications and no details have been investigated.

2.4 Determining the Appropriate Perturbation Parameter

As stated earlier, there might be more than one possible perturbation
parameter for a two-time-scale system. In addition, there is no systematic way to
identify the most "appropriate" one for a problem under investigation. Some
insight into the problem under investigation is needed to enable the user to
acquire this information. The objective of this section is to gain that insight and to
determine the appropriate perturbation parameter s for the problem under
investigation. To the author’s knowledge, there are no published efforts that
discuss this issue in detail. Zeid (1988) suggested (without analysis) using two
different perturbation parameters, C and R with C being "very small“ and R

being "very big”. This suggests that RC is 0(1). On the other hand, some

simulation results, which were obtained within this research work and will be
shown later, show that the modified system can include two fast poles. These
results showed also that the separation between the slow and the fast sub-

systems increases as RC decreases. This suggests using RC as the

25

perturbation parameter. This section will show another possible choice for the
perturbation parameter; later it will be demonstrated that this choice helps reveal

the time-scale properties of the systems under investigation.

2.4.1 Analysis Based on Physical System Modeling

Figure 2.5 shows the simplest version of the problem under investigation.

It shows a case in which the original system consists of only two “coupled”

masses, [1 and I 2.

 

 

 

Figure 2.5 A local problem

Let us assume that these two coupled masses are decoupled by inserting a
parasitic spring, whose compliance is C and a parasitic damper whose
resistance is R. This is exactly equivalent to using the method under

investigation to eliminate the derivative causality. Assume that the translation of

26

the two masses are x1 and x2 respectively. The equations of the modified

system can be written as
1 R
5612—— x2—x1+— XZ—Xl
IlC( ) Il( )

1 R
"2=—— 2— 1—— xz—xi
x 12C” x) 12( )

5 =(21—5cz)

This system can be written in the state space form with the state vector

j 5 ;1_
x1 11 11 C11
. . R —R l
X = x2 and system matrix A: — — —
12 12 CI2
6 1 —1 0
_ I

 

 

It can be shown that the characteristic equation of this system is

222+ £+£ +i —l—+—L =0
11 12 C II 12'

where, 11 is a root of the characteristic equation; and/l = 0 root represents the
slow “original system” eigenvalue. To find the eigenvector that corresponds
to/i = 0, set ,1 = 0 in the equation

AV = 2V
Where, V is the eigenvector that corresponds tol. It is easy to show that for
,1 = O ,

V1

27

This mode represents a solid body mode in which the two masses have the

same velocities; this is an accurate description of the original system.

The other two eigenvalues of the modified system are the eigenvalues of

the fast sub-system and they are the solution of

12+2[%+I£2)+%(71I+715]=0 (221)
As the best separation will occur when the two fast poles are equally fast which
means that the sluggishness of neither of them would dominate. We examine
here the condition at which the two fast poles of (2.21) are critically damped. It is
easy to show here that for roots of (2.21) to be critically damped, the following
condition (2.22) must be satisfied

4
1 1

II 12

RRC=

 

(2.22)

Equation (2.22) will be examined subsequently, together with the results from the

next sub-section.

2.4.2 Simulation-Based Analysis

In another effort to discover the appropriate perturbation parameter for the
problem under investigation, simulations of some simplified versions of the
problem were conducted. In one attempt, Example 1 was used; Figure 2.6
shows a schematic diagram and a bond graph of the original system of Example

1. Where, the original system refers to the system with the derivative causality.

28

As shown, the original system includes two inertances, L1 and L2, that are
coupled through a transformer, TF, whose modulus is m. Figure 2.7 shows a
schematic and bond graph model for modified system of Example 1. Where, the
modified system refers to the original system after inserting parasitic elements
into it to eliminate the derivative causality. The terms original system and
modified system will be used throughout this thesis to convey the above-

mentioned meanings.

TF RI

 

 

Figure 2.6 (a). Schematic diagram for the original system of Example 1

—L|1}-—, TF}-——,1|—712

I I

11 R1

Figure 2.6 (b). Bond graph model for the original system of Example 1

29

L1 L2

lllli

TF F“

 

Figure 2.7 (a). Schematic diagram for the modified system of Example 1

11 :l—I—r 1&0

Figure 2.7 (b). Bond graph model of the modified system of Example 1

The following set of parameters was used for Example 1 in this attempt:
I 1 = I 2 = 2000
R1 = 2000
m =1
Eigenvalues of the modified system were calculated at different values of the

product RC that were chosen. For every chosen value of RC, the values of R

30

and C were changed till a pair of two critically damped poles for the fast sub-

system are obtained. For these values of R and C for which two critically

damped fast poles occur, the product RZC was recorded vs. its corresponding
value of RC. Again, the critically damped condition was chosen because it was
thought to provide the best separation between the fast sub-system and the slow
sub-system. To obtain suitable values for the fast poles some trial and error
alongside with knowing, that generally reducing RC will produce faster poles in
the fast sub-system were used. Figure 2.8 summarizes the results of this

investigation for Example 1. The horizontal axis in Figure 2.8 is the product RC
and the vertical axis is RZC normalized by dividing it by the right hand side of

equation (2.22).

 

 

100.0%
98.0%
96.0%
94.0%
92.0%
90.0%

RRC Normalized

88.0%
86.0%

 

 

 

 

 

 

 

Figure 2.8 Relation between normalized RRC and RC, Example 1

31

To further the understanding of the results shown in Figure 2.8 the
procedure described above in this section were implemented again for another
example (Example 2). Figure 2.9 shows a schematic diagram and a bond graph
model of the original system of Example 2; and Figure 2.10 shows a schematic
diagram and a bond graph model of the modified system of Example 2. The
difference between Example 1 and Example 2 is that Example 2 has an
additional energy-storing element (capacitor, C1). Consequently, Example 2 is
more general than Example 1 and its order is higher than that of Example 1 by

one.

L1 Cl L2

TF F“

 

Figure 2.9 (a). Schematic diagram for the original system of Example 2

32

C1

1

—l>|1|.__> TFH1Ih—AR1

I i

l1 I2
Figure 2.9 (b). Bond graph model for the original system of Example 2

L1 C1
1 I I I I
TF RI

Figure 2.10 (a). Schematic diagram for the modified system of Example 2

IIIll

33

Figure 2.10 (b). Bond graph model for the modified system of Example 2

The same procedure described above for Example 1 was applied for Example 2
for three different cases; every case had a different set of original system
parameter values. The three different sets of parameters were chosen such that:

1) One set produces an over damped original system, one set produces a
critically damped original system, and one set produces an under damped
original system.

2) The real part of the fastest pole in all cases is equal to —0.25. This will
make comparing the three different cases with one another and with
Example 1 easier and more meaningful.

The used three different sets of original system parameters along side with their

corresponding eigenvalues are as follows:

34

Case one:

C1 = 0.0107 11 = 12 = 1250 R1 = 1000
ll=-0.25 x12=-0.15 m=1
Case two:

Cl=0.008 11:12:1000 R1=1000
21: —0.25 22 = —0.25 m = 1
Case three:

Cl=0.004 11:12:1000 Rl=1000
21: —O.25 + 10.2501 112 = —O.25 — i0.2501 m = l

The results of these three cases are summarized in Figure 2.11

 

 

[—o—c1=o.oro7 +cr=oooe +C1=0004 ]

1 00%

 

90%

l

80% —

Normalized RRC

 

70%

 

 

 

 

 

 

Figure 2.11 Relation between normalized RRC and RC, Example 2

35

Figure 2.8 and Figure 2.11 indicate that the product RZC converges to a

constant as RC gets smaller. In other words, RZC converges to a constant as
the time-scale separation between the slow sub-system (original system) and the

fast sub-system increases. Moreover, Figure 2.8 and Figure 2.11 indicate that

the constant (to which RZC converges) could be related to the original system
parameters only. These conclusions agree with and perhaps explain equation
(2.22) that was derived for the two coupled masses example. After examining
Figure 2.8 and Figure 2.11 it is easier to understand that equation (2.22) was
obtained because the original system had a zero eigenvalue; therefore any value

for the two critically damped fast eigenvalues will provide enough time-scale

separation between the original sub-system and the fast sub-system forRZC to

converge to . The observation that RZC converges to a constant as e

 

l 1

ll 12

tends to zero means that there is another option for choosing a in addition to
£=RC; it is s=%. These results will aid the efforts to find the standard

singularly perturbed model for the class of systems under investigation. These

efforts are reported in the next section.
Although examples used in this section are somewhat restricted, they

were useful in demonstrating some important aspects of the general problem for

these reasons:

36

1)

2)

3)

They contain the important aspects of the problem under investigation;
they contain a coupling between two energy-storing elements in the
original system; and they contain the same structural modifications under
investigation in this work. The modifications do not depend on the details
of the original system.

As the time-scale separation between the slow sub-system and the fast
sub-system in modified system increases; most of the details of the slow
sub-system (the original system) are not going to affect the fast-sub-
system. Later, it will be shown that, only the inertia field will affect the fast
sub-system as the time-scale separation increases.

As the time-scale separation between the slow sub-system and the fast
sub-system in the modified system increases; the effect of the fast sub-

system on the slow sub-system tends to vanish.

2.5 Finding the Standard Singular Perturbation Form

The results obtained from section 2.4 are more of clues that provide

insight, which is typically needed before being able to use the singular

perturbation method effectively. In addition, until this point of this presentation, it

has not been confirmed yet that choosing 8:115 and assuming that RZC

converges to a constant as 8 tends to zero, will be a better option for the
problem under investigation. A subsequent analysis based on these results will

prove their usefulness. An important part of this analysis is formulating the

37

modified system in the standard singular perturbation form. This issue is

investigated in this section.

2.5.1 Using Transformation Procedures to Obtain PSSPF

Even though there are some general procedures to help transform a two-
time— scale singularly perturbed model from a non-standard form to the standard
form (Kokotovic et al. 1986), the result of these procedures are not unique; and
you have to decide what a is and put the model in a required form before you
can start using those procedures. This section shows how the transformation
procedure outlined in Kokotovic el al. (1986) can be best used for the classes of
models under investigation. It demonstrates that by showing how to use the
procedure in transforming two simple examples; then it characterizes the

preferred standard singular perturbation form (PSSPF).

This procedure assumes that the systems are linear; they can be

summarized in three steps as follows:

$21.3

Put the system in the form
81? = F (£)v, v e Rn+m (2.23)
Where, n is the order of the slow sub-system and m is the order of the fast sub-

system.

Let F(£) = FO + €F1(8)

38

For two or more time-scales to appear we should have detF(0)=0.
lfdetF(0)¢0 then all the components of v would be fast. Assume that the

dimension of the null space (N) of F0 is equal to n. To be able to find the
transformation, rank of F0 must equal to the number of the fast variables. The

null space (N) of F0 is defined as
_ n+m _
N(F0)_{yeR ’FOy—O}

Step 2:

To select a fast variable vector 77 choose m linearly independent vectors in R’H'm

orthogonal to the null-space N(F0) and arrange them as the rows of an

mx(n + m) matrix Q.

Step 3:

To select a slow variables vector x, find P (the left null space of F0) such that

PF

0:0

The change of coordinates

x = PV, 7) = Qv; T 2 IS] (2.24)

will transform the system (2.23) into the standard form

if = A1 1(6)): + A12 (5)77 (2.25)

err=eA21<e>x+An<e>n (2.26)

39

as long as the transformation matrix T is nonsingular. In addition, as 13
becomes small, A11(0) will contain the it slow eigenvalues and £A22(0) will
contain the m fast eigenvalues in the original time-scale I (vs. the fast time-scale
r=—;-), Kokotovic et al. (1986). The next paragraphs will demonstrate how to

use this general procedure by explaining how to apply them to two simple

examples.

Following the Sequential Causality Assignment Procedures (Rosenberg et
al. 1983) it can be shown that the homogeneous original system of Example 1

can be written as

1131 = _ leL1 43112. (2.27)
dr Il dt

102. 21

12 m“ (2.28)

and that its homogeneous modified system can be written as

2
dpl m (R+Rl) mR m
——=———— 1+— 2—— 2.2
d: 11 p 12" Cq ( 9)

dp2 mR R l
——=— 1—— 2+— 2.
dt 11 p 12” Cq ( 30)

dq m 1
—=— l—— 2 2. 1
dt 11p 12" ( 3)

Where i1 and i2 are the currents in the inductors L1 and L2 respectively; q is
the charge in the parasitic capacitor; pl = I 111 ; p2 = I 2i2. Next we will apply the

transformation procedure to convert the modified system to the standard form.

40

First, the modified system equations (2.29) - (2.31) can be expressed in the

 

matrixform
F—m2(R+Rl) mR —m
pl 11 12 C
1’. p2 = _R_':'. :5 -1. (2.32)
dt 11 [2 C
q 1': :_l 0
II 12

 

 

Inspecting equation (2.32) shows if we proceed from here to calculate F0, the
rank of F0 would equal to one. Since it was stated above that the number of fast
variables should be equal to the rank ofFO, and we already decided that there

should be two fast variables; this indicates a potential inconsistency, and means

that some scaling is needed to avoid this potential inconsistency. Knowing that

s = %, and that RZC becomes constant as 5 becomes small; an experienced

observer can notice that he/she needs to scale the variable q to Rq to be able to

proceed correctly. As R is equal to i, this scaling is an e-based scaling, which
a

is exercised only to avoid potential inconsistency such as the one detected

above. Implementing this scaling and substituting a for RZC will produce

1' 1

 

—m2(R+Rl) mR —mR
pl 11 I2 a
R —R R
1 p2 = ._ _ _ (2.33)
dt 11 12 a
Rq m_R j 0
11 12

 

 

Consequently, we can obtain

41

 

—m2(l+£R1) _111_ -mT
pl 11 12 a
d m -1 1
g__ 2 : — — — 2.34
cit p 11 12 a ( )
Rq m —l
— — 0
ll 12

 

 

Comparing (2.23) and (2.34), we can see that

-

:12. 2 :2
11 12 a
m —-1 1
F = — — —— 2.35
0 11 12 a ( )
fl :_1. 0
11 12

 

 

and that the rank of F0 is two, which is the correct rank.

In the second step, to select the fast variables, find the right null space (y) of

F0. Calculating y from F0 y =0 leads to

 

 

where yl is an arbitrary constant. Notice that the dimension of this null space is

one, which is equal to the number of the slow variables. The next step is to

choose two linearly independent vectors that are orthogonal to y. Obviously, the
choices of these vectors are not unique. This research work recommends

choosing the two orthogonal vectors to be

42

_ fl _
I l 0
:—l and O (2.36)
I 2
l
0

 

 

This choice is recommended because the two fast variables this choice is going
to produce are physically meaningful and moreover, understanding and
characterizing their behaviors are useful to this research work. With this

recommended choice the two fast variables are

m l

— l—— 2 2.37
11 p 12 p ( )
Rq (2.38)

The first fast variable (% p1———112— p2) can be seen as the error introduced into the

second equation (2.28) of the original system model due to inserting the parasitic
elements. The second variable (Rq) is a scaling of the charge in the parasitic
capacitance; it is basically, the integral of the first fast variable. From design
perspective, it is important to design the modified system such that these two
variables are fast; so that the effect of the parasitic elements on the original
system behavior will vanish quickly. Therefore, these two variables are of
particular importance to this investigation.

In the third step, to select the slow variables, find P (the left null space of F0)

such that PFO = 0. This will produce

P=[Pl mPl 0] (2.39)

where, P1 is an arbitrary constant, it will be chosen to equal to one; so that

43

P = [1 m 0] (2.40)
With this choice, the slow variable in this example is going to be
p1 + mp2 (2.41 )

This slow variable has a physical meaning. If we were to algebraically
manipulate the original system equation (2.27) to sum the coupled derivative
terms together, this algebraic manipulation would have produced the derivative of

this variable (2.41) as the sum. With P and Q available, the transformation

matrix T can be written as

P 1 m 0
T: =_"l;lo (2.42)
Q 1112
0 01

Consequently, the system matrix of the transformed system is

 

 

 

 

 

 

-m2(R+Rl) mR -mR
11 12 a
3:7" £5 j .1: T_1,or
ll 12 a
m_R LI: 0
11 12 j
- _ 2R1 _
m2 0
Il+m12
~ —m3R1 1 R m2 1
A: 2 2 — —-+— —— —+—-
11(11+m 12) 1112(11+m13 11 12 a 11 12

0

 

L

R

0

 

.1

Therefore, this system can be written in the standard singular perturbation form

5: = A11(8)x + A12(£)n

577 = 9‘21“.” + 1422(8)”

(2.25)

(2.26)

with

3.". 1___1_ 2
x=pl+mp2. 77= 11” 12p 1
Rq

A _ —m2Rl A1 _ —m312R1 0
11 (11+m212)’ 3 IIl+m212I ’

 

 

l—m3Rl 2 —m14122R1.e,_ m2 +_1_ _1 23+:
A21: 11(11+m212) rand A22: 1112(11+m213 11 12 a II 12
0 1 0

Examining the eigenvalues of the original system shows that it is really equal to

the one indicated by Al 1 as it was stated. It was also stated above that

l1122(0) contains the eigenvalues of the fast sub-system as 8 gets small.
1:

43.; _1 12,;
1122(0): 11 12 a 11 12. (2-43)
1 0

Using equation (2.43), the eigenvalues (11) of the fast sub-system, as 6 gets

small, can be calculated from
2 2
2+ Ln—+—l— R +-1— -m—+—l— R2 =0 (2.44)
11 12 a ll 12

Next we will demonstrate how to apply the procedure to Example 2.
Example 2 is more complex than Example 1; its original and modified systems

were given above in Figure 2.4.5 and Figure 2.4.6 respectively. Following the

45

SCAP (Rosenberg et al. 1983), the homogeneous original system of Example 2

can be formulated as

@2121
d1 11
dr C1 11 d1
12sz
12 II

(2.45)

(2.46)

(2.47)

and its homogeneous modified system can be written as

5111.12.21
d1 11

m
2.__
Cq

(2.48)

(2.49)

(2.50)

(2.51)

1) After scaling q to Rq, this system can be written in the matrix form as

 

 

 

 

0 i 0 0
f 1‘ 11
q] :_1 —m2(R+Rl) mR -—mR
i p = Cl 11 12 a
dt p2 0 m_R :5 5
_Rq_ 11 12 a
0 1R: 1 0
_ 11 12
Therefore,

46

—

(2.52)

 

 

10 0 0 0
0 _m2 1: I
ll 12 a
F0 = 0 1'1 _-_1_ _1_ (2.53)
11 12 0!
0 fl 3—1 0
_ II 12 -

 

 

2) To select the fast variables, find the right null space (y) of F0. Calculating y

from F y =0 leads to

 

 

 

 

0
_ yl . '1' _ O _
12 0 1
y= m12 =11 +y2 m12 (2,54)
——y2 0 —
11 11
0 - _0_ _ 0 _

 

 

Where yl and y2 are arbitrary constants. Notice that the dimension of this null
space is two, which is equal to the number of the slow variables. To select the
same meaningful fast variables we selected for Example 1 the two vectors that

are orthogonal to the right null space should be selected as

— -

 

 

 

 

0 _
L". i0
11 0
and 2.55
__1 O < )
12
l
L 0 _ _ J

Therefore, the two fast variables remains to be

m l

— 1-— 2 2.56
11p 12” ( )
Rq (2.57)

47

3) To select the slow variables, find P (the left null space of F0) such that

PFO = O. This will produce

P11 P12 P12 0
P: m (2.58)
P21 P22 mP22 0

Where, P11, P12, P21, and P22 are arbitrary constants to be chosen so that
the transformation matrix is not singular. We still desire to get the meaningful
slow variable that was obtained in Example 1. To do that select

P21: 0, P22 =1 (2.59)
This will select one of the two slow variables as

pl + mp2 (2.60)

The other meaningful slow variable would be the state (ql) in the original

system. To select it select
P11=1, P12=0 (2.61)

With P and Q available, the transformation matrix T can be written as

q

)—

 

l 0 0 O
P 0 l m 0
T = = m -1 (2.62)
Q 0 — — 0
ll 12
O 0 0 1

 

One can easily notice that the lower right 3x3 sub-matrix in the transformation
matrix of Example 2 (2.62) is the same as the transformation matrix of Example 1

(2.42). The system matrix of the transformed system is

48

 

 

 

 

 

 

 

 

 

 

0 i 0 0
11
—_1_ ~m2(R+R1) mR -mR
3:1 Cl 11 12 a T—l,or
o m_R j 5
11 12 a
0 3'5 j. 0
I. II 12 .1
’ 1 m12
2 2 O
11+m 12 11+m 12
—_1 -m2Rl -m312Rl 0
A: C1 11+m’12 Il+m212
—m —m3R1 —m‘12R1 _ i _1_ i£+i
IlCl 11(11+m212) 11(11+m212) 11 12 a 11 12
_ 0 0 R 0
‘11 fl 1—-l—p2
x: , 17: up 12 ,and
p1+mp2
Rq
—m412R1£ m2 1 —l m2 l 1
mummy 77+}? 7; THE
A..(e)= (2.63)
1 0 j
Consequently,
I m2 1 —1 m2 1 _
‘ 71+}? "07 7W5
A..(°)= (2-64)
1 0

 

 

It is remarkable that A22(0) for Example 1 (2.43) is the same as that for Example

2 (2.64). In other words, the fast sub-systems in both examples will become

same as 8 decreases. More details will be given on that in CHAPTER 3.

49

As explained while demonstrating how to use the procedure on the two
examples, there are some variables that are preferable, to be states in the
transformed standard form, to others. The standard form that includes these
preferred variables will be called the preferable standard singular perturbation
form (PSSPF). This PSSPF has the following characteristics:

1) The slow variables consist of the states that are not part of the coupled
field in the original system, and the variables whose sum of derivatives are
obtained if we were to sum the coupled derivatives (in an equation-by-
equation basis) in the original system as shown in the examples.

2) One of the fast variables is the error introduced into the original system’s
algebraic equation due to inserting the parasitic elements.

3) The other fast variable is a scaling of the charge (Rq) in the inserted

parasitic capacitance.

2.5.2 A Transformation to Obtain the PSSPF

This section derives the transformation matrix needed to transform a

general modified system from the class of systems under investigation to the

PSSPF and proves its existence.

Assume that a linear system includes a coupled energy field together with

some other energy storing elements that are not coupled with the field. Assume

50

that the vector of the states that are not among the coupled field is X , and its

dimension is r. One way to recognize the elements of the coupled field is by the

appearance of their parameters in J Further, assume that the coupled field

d1”
includes only one dependent energy-storing element. Based on general form for
the coupled field given in equations (2.16) and (2.17), the general form of this

class of original systems can be written as

X, = AWX, + A,X, (2.65)
X, = 1,, *S, *X, +1,, *X, (2.66)
s, * X, = 1,, * S, * X,. (2.67)

Consequently, the general form of the modified systems of this class, and based

on the general form given in equations (2.18) to (2.20), can be written as

 

 

 

A" A, 0 0
FXr- l ”X,.-1
A. J..S.+RJ.J.S. -RJ.S _1-
X, rr 11 1 rd d1 1 rd (1 11! X
i ' = C1 ' (2.68)
‘1‘ X4 0 RJ,,S, -RS, — Xd
_ q 1 C _ ‘1 1
0 1,5, —S, 0
_ J

 

 

Where, —é—= K p. The transformation matrix that transforms this modified system

to the PSSPF is

' I, 0 0 0 '
o 1, -1,, 0
T: (2.69)
0 1,5, —S, 0
_ O O 0 R J

 

 

If the modified system is in the explicit form

51

X=AX+BU
then, the transformed modified system will be

Z=TAT"Z+TBU, Z=TX
Where I i is a square unit matrix whose dimension is the same as that of X i and
I r is a square unit matrix whose dimension is the same as that of X ,. Notice
that because there is only one dependent energy-storing element, the dimension

of S d is one by one and the dimension of X d is one. Further notice that one

does not have to formulate the original system to be able to formulate the
modified system and transform it to the PSSPF. To obtain the transformation

matrix given in (2.69) one needs to know Ji d’ S,, and S d' The matrix

J di = 4;, which relates the flows (velocities) of the dependent and independent

energy storing elements can easily be obtained from the original system bond

graph; and S,- and S d are inverse inertia matrices that are parts of the original

system parameters.
2.5.2.1 Derivation and Proof of Existence
To derive equation (2.69), we will start from the general form for the

modified systems given in equation (2.68). This form is rewritten again here in

equation (2.70).

52

F"

 

 

 

—

A

IT

A.

If

0

0

Ari
11151 + RerJdrSr
111.5.

JdiSi

0

- R1,, 5,

- RS,

_Sd

 

0
, _
1] Xr
_ id X.
C ' (2.70)
i X.
C _ q)
0

 

Changing q to Rq and substituting a for RZC in equation (2.70) will lead to

 

 

 

 

 

A” A, 0 0
IX: R IX,—
d X, _ Air J,,S,+RJ,,,J,,S, TRJI'de 31111 X, (271)
‘1’ X4 0 111,5, —RS, 5 X4 '
Lqu a _Rq_,
0 RJ,,S, —RS, 0
Equation (2.71) means that
0 0 0 0
l
0 J. J .S. J. S —J.
F0: rd d1 1 1d d alzd (2.72)
0 J .S- —S —
d1 1 d a

 

 

 

 

and that the two orthogonal vectors to the null space can be chosen as

53

 

 

 

 

0 '0'
(1,3,1 0
and
- S d 0
0 _1,

O J .S. —S 0
_ d1 1 d
Q ‘I 0 0 0 1 ] (2'73)
To select the slow variables, find P (the left null space of F0) such that PF0 = O.

This will produce

P Pri —PriJid O
P: P" —P J 0 (2.74)
ir ii ii id

Where, P”, Pri’ Pir’ and P11 are arbitrarily chosen so that the transformation

matrix is nonsingular. To get the meaningful slow variables that we want (as

explained in the examples above), select

P. :0, P..=I. (2.75)
1r ll 1

P :1, P.=O (2.76)
17' r 7'1

combining the change of variable q to Rq with P and Q, we can get the

transformation matrix T as in equation (2.69) to be applied directly to the general

form (2.70).

Next we will show that the transformation matrix given in equation (2.69) is

always nonsingular.

S4

det(T)= R det(I,)det(Ii)det (-sd+1d, 5, 11711”)

= i R det(sd+1,T 5,1,d)

d

Since Si and S d are inertia matrices, they must be positive definite. Therefore,

Sd>0,and

1,281.11. d 2 O, with the inequality satisfied when Jid is nonsingular.

This proves that det (T) at O which means that the transformation exists.
2.5.2.2 Example of Using the Transformation Matrix

This sub-section shows an example of using the transformation matrix
defined above to transform a modified system to the PSSPF. Figure 2.12 and
Figure 2.13 shows a schematic diagram of a fluid pipe system and its bond graph
model respectively. Figure 2.14 shows a bond graph model for the modified

system.

55

I3 ,R3
Arrl [1,Rl

Efi / / 7
511—4 D K1/14,R4
s..__. Q‘\ k,

Arr2 [2 ’R2 \

 

 

 

 

 

 

 

Is,R5

Figure 2.12. Schematic diagram of the fluid pipe system

 

 

 

 

1.—
-(P-—

[3

[4

 

 

Carri RI 1' i / 1
S121 I
1 I 0 /1 1
Carl R2 ’2

Figure 2.13. Bond graph for the original fluid pipe system

56

 

 

 

 

 

I I
5"1 / 0 ’I 1 F
Carrl R. [I
512} 7 0 7Jl 1
Car-r2 R2 [2

 

 

 

fl

I

0 7i”
I

11L 7R
I

C

Figure 2.14. Bond graph for the modified fluid pipe system model

For this system

pl
p2
p3
p4

- d

 

 

Xd =[P5I1

57

_ E _1 _ _ _
S,_ 1 , Sd—[—], 1d,_[11 1 1]

 

 

:44; ,_2
dt 11
£13:sz 2
dt 12

 

II 12 13 I4 Is [1 12 I3 I4
5122.. 1 q,_R,_123_iq-12.2_2_2_2)_R{2+2_2_2)
2

dp1_ 1 q,_R,2_%q_R(2,_12_12_2:_2)_R{2+£3_23__1:)
1

 

Ii 12 I3 14 Is It 12 13 I4

dt C It 12 13 I4 Is Ii 12 [3 14 I4
5122-: 2.2_2_2_2
d1 C (11 12 13 I. 15

12-21.222-222
dt 11 [2 13 I4 15

Consequently, the modified homogeneous system can be written in the matrix

form as

58

 

 

 

 

 

 

 

 

 

 

J

 

 

 

 

 

 

 

() () -—— I) () () I) 0
Ir
F-q - 0 0 0 -;;- 0 0 0 0
I 2
m l _R+&+Rs *R+Rs R+Rs R+Rs __ _1
C1 Ii Iz I3 I4 Is C
P‘ 0 L _R+Rs _R+R2+Rs R+Rs R+Rs i -i
ji_ ’72 :z (T2 It 12 [3 I4 Is (3
d; p3 0 0 R+Rs R+Rs _R+R3+R5 _R-l-Rs _fl 1
P4 II 12 I3 I4 Is C
0 0 R+R5 R+R5 _R+Rs _R+R4+Rs _R 1
P5 11 12 13 1. 15 C
-q- 0 0 5 5 -5 -3 -5 _1_
Ii 12 I3 I4 Is C
0 0 _1_ _1_ __1_ _i _i OJ
_ 11 12 I3 I4 Is

Using equation (2.69) the transformation matrix can be written as

r- '1

 

 

10000 0 00
01000 0 00
00100 010
00010 010
T=000010—10
000001—10
00ii__1___1-_i0

II 12 13 14 Is
_00000 0 0R,

With this transformation, the states of the transformed modified system are
I' q, "
qz
p1+ ps
p2 + p5
p3 - ps
[)4 - p5
111,33 P3 P4 P5

’1 I2 ’3 I4 ’5

 

 

— .1

Notice that the original system equations are the following algebraic differential

equations

59

——-=S __
dt fl 1]

143._S,,_L

dt 2

dP1 1 P1 dPs P1 P2 P3 P4
—=—q1-R1—-——-R — —————
dt C1 Ii dt 11: 12 13 I4
sz 1 P2 dPs P1 P2 P3 P4
—=—q2-R2—-—- — —————
dt 2 12 dt [1‘ 12 I3 I4

60

CHAPTER 3
DESIGN OF A CLASS OF TWO-TlME-SCALE MULTIPORT MODELS

3.1 Introduction

Analysis carried out in CHAPTER 2 gave insight into the problem under
investigation; and showed how to transform a modified system from the class of
systems under investigation to the PSSPF. This sets the stage for taking
advantage of several tools that are provided by the singular perturbation method.
CHAPTER 3 will take advantage of some of these tools, together with some
other tools, to design the class of systems under investigation. The design
objective, which is the same as the chapter objective, is to achieve the greatest
possible simulation efficiency for a specified approximation accuracy by making
the appropriate choices for the values of the parasitic element parameters. As
explained earlier, the computational efficiency and the approximation accuracy of

the modified systems can be traded off together.
3.2 A Model for the Fast Sub-System and Proof of Stability

As indicated earlier, if the modified system is correctly designed, it should
be a two-time-scale system, so that the modifications would not greatly affect the
behavior of the original system’s variables. This two-time-scale system should
contain two sub-systems; the slow sub-system and the fast sub-system. The

slow sub-system represents the original system, which is known and there is an

61

obvious interest in not altering its behavior. The fast sub-system results from
inserting the parasitic parameters. It is well understood that the fast sub-system
will affect the behavior of the slow (original) one. Unfortunately, there is not
much known about the fast sub-system. Particularly, there is no published work
that characterizes or models it. Therefore, the first step in this effort of designing
the fast sub-system is to model it. This section will carry on some of the results
obtained in CHAPTER 2 to achieve this task. In addition, the results obtained in
this section directly address the stability of the fast sub-system and

consequently, of the whole modified system.

Work presented in CHAPTER 2 derived a general form for the

transformation matrix T and the matrix FO that is used to select the fast and

slow variables, so that a modified system in the general form (3.1) can be
transformed to the standard form (3.2) and (3.3). These general forms are

rewritten again here, for the reader’s convenience.

 

 

 

 

 

A” A". 0 0
”Xr- 1 -Xrfi
d X, Air J,,S,+RJ,de,S, _RJI'de —Jid X.
_ r = C1 ' (3.1)
dt X" 0 RJm'Sr _RSd — Xd
L q .1 C — q J
0 1,5, —5, 0
X = A11(£)x + A1 2(5)17 (3.2)
877: £712, (8)x+A22(8)17 (3'3)

62

— q

0 0 0 0
l

 

 

 

 

0 J. J .S. —J. S —J.
d d 1 1d d 1d
F0: ’ ‘ a 1 (3.4)
_O Jd1'Si _Sd 0 —
’1, 0 0 0‘
T: F = _0___..I_"____—.{'L_9 (35)
Q 0 1.8.- S. 0
_0 0 0 1,
Kokotovic et al. (1986) showed that
422 (0) = QFOW (3.6)
where,
T“=[V W] (3-7)

As the general form for both F0 and Q have been derived in CHAPTER 2 (they

are rewritten above as (3.4) and (3.5)), the only remaining piece that is needed to

calculate A22(0) from equation (3.6) is W. To find W, we will derive T" using

equation (3.5), and then use equation (3.7) to find W. From equation (3.5)

1;‘ 0 0
_1_ -1 _ Ir —Jr'd
T — 0 G 0 where,G—

O 0 1 111151 _Sd

Using an inverse block matrices formula (see for example Kailath 1980), it can be

shown that

63

 

 

 

 

 

 

 

 

F141 _ 11.11.1151 11.1 -
G‘l — 53-514‘5‘1“ _S" +{J'S‘J‘d and consequently, it can be
(If i
_ Sd -J,,,S,J,d —Sd +JdrSrJ111 _
shown that
r _
0 0
Jid 0
W = "Sq +JdiSiJr'd
l 0
-Sd +1105:er
0 l
_ _I

 

 

Therefore, using equation (3.6) it can be shown that

_(Sd TJLSIJM) 1(84 +J;S,J,,)
422(0) = QFoW = a (3'8)
1 0

and consequently, as 5 becomes small, the fast sub-system dynamics can be

written as

_ 11,, ,0) = (3.9)
8 R 0

Assuming that s is small enough, equation (3.9) models the fast sub-system.

The term

1

(s, + 1,3; 5.1..) E 1" (3'10)

 

is the equivalent of the inertia field evaluated at (reflected on) the dependent

energy-storing element. Therefore, equation (3.9) can be written as

 

A = '1 '1 (3.11)

 

 

where A, refers to the fast sub-system dynamics, and 1,, refers to the

equivalent of the inertia field evaluated at (reflected on) the dependent energy-

storing element. Using equation (3.11), the characteristic equation of A, can be

written as

22 +£2+ 1
I CI

eq (q

 

= o (3.12)

Equation (3.12) is the characteristic equation of a series R-L-C circuit, with

the term 1,, as the inductance of L. This result shows that, as 8 gets small, the

fast sub-system can be approximated by the parasitic capacitance, the parasitic

resistance, and 1,, connected in series. This model should be very useful for

any design work for the two-time-scale modified system; and it will be used in the
design algorithm to be presented later in this chapter. The condition, of having a
small enough a, imposed on the validity of this model is not detrimental at all for
our framework, because it is the same condition that applies for any useful
modified system. Equation (3.12) also proves the stability of the fast sub-system
for a small enough a, which means that the whole two-time-scale modified
system will be stable, if the original system is stable (Kokotovic et al. 1986). The
model of the fast sub-system described above can be shown in a graphical form.

Figure 3.1 shows a bond graph model of the fast sub-system.

65

Figure 3.1. A bond graph model of the fast sub-system

Notice that the parasitic elements R and C are linear by choice. The model
given in Figure 3.1 is derived for linear systems; however, it will be useful for

nonlinear systems too. That will be shown later in CHAPTER 4.

3.3 Validating the Fast Sub-System Model

The model of equation (3.11) for the fast sub-system is very important; but
it should be validated before it is used. As the model is an approximation that is
good when 8 is "small enough," this validation should also help us to get an idea

of what 'small enough 8 " means.

3.3.1 Accuracy of Slow and Fast Poles Locations

In this sub-section we examine the accuracy of the fast and slow poles of
the modified system after using the fast sub-system model to choose the values
of the parasitic elements to place the fast poles at certain specified locations. The
first example used is the electrical circuit that was introduced as Example 2 in
CHAPTER 2. The homogeneous original system of Example 2 that was given in
equations (2.45) to (2.47); and its homogeneous modified system that was given

in equations (2.48) to (2.51) are rewritten here again for the reader's

convenience.

10:1 =5:11 (2.45)
1571 = —éql — n{lef—ll+%) (2-46)
22. = mfl (2-47)
12 11

if = 5.11. (2.48)
%=_éql_flnfl‘fllpr+’y_§p2—gq (2.49)
if=lgpcgpz+gq (2.50)
j—:I=%pl—71§p2 (2-51)

67

Three different cases of Example 2 were examined in this validation work; every
case had a different set of original system parameters values. The three different
sets of parameters were chosen such that:

1) One set produces an over damped original system, one set produces a
critically damped original system, and one set produces an under damped
original system.

2) The real part of the fastest pole in all cases is equal to —0.25, to make
comparing the three different cases with one another easier and more
meaningful.

The used three different sets of original system parameters along side with their

corresponding eigenvalues are as follows.

Case one:

Cl=0.0107 11:12:1250 Rl=1000 m=l
21: —0.25 22 = -—O.15

Case two:

Cl=0.008 11:12:1000 Rl=1000 m=l
21 = —O.25 22 = —0.25

Case three:

Cl=0.004 11:12:1000 Rl=1000 m=1

111 = —0.25 + i0.2501 22 = —0.25 - 10.2501
We assume that we desire the two fast poles to be critically damped and to be at
a certain desired pole location (DPL). Then, equation (3.12) is used to calculate

the values of the parasitic parameters that will place the fast poles at the desired

68

location. These values of the calculated parasitic parameters are used in the
modified model and subsequently, the eigenvalues of the modified model are
calculated to examine how close the fast poles to the specified desired ones and
how close the slow poles to the original system poles. Notice to use equation

(3.12) to calculate the values of the parasitic parameters, we need to know 1,], in
addition to the two desired location of the poles. 1,, Can be calculated using

equation (3.10); it can be shown it is, for this example, equal to

l

l l
_+__.
Ii 12

 

Tables 3.1 to 3.3 show the calculated eigenvalues from the modified system (as
described above) for every desired pole location (DPL) that was examined; the
three different cases of the example are shown. The first column in every table
lists the eigenvalues of the original system. PR is the ratio between the used

desired pole location and the real part of fastest pole in the original system.

 

 

 

 

 

 

m DPL=-1 .25 DPL=-2.5 DPL=~5 DPL=-7.5 DPL=-12.5 DPL=-17.5 DPL=-25
PR=5 PR=10 PR=20 PR=30 PR=50 PR=70 PR=100

-2.2744 -3.7827 -6.6724 -9.4794 -14.9726 -20.3766 -28.3881

-0.5709 -1.6101 -3.7261 -5.9199 -10.4271 -15.0233 -22.01 19

-0.25 -0.3093 -0.2592 -0.2529 -0.251 9 -0.251 4 -0.251 3 -0.251 2
-0.15 -0.1454 -0.148 -0.1486 -0.1488 -0.1488 -0.1488 -0.1488

 

 

 

 

 

 

 

 

 

Table 3.1 Validating the fast sub-system model against Example 2, case 1

69

 

 

 

 

 

 

 

 

 

 

 

 

 

Orig DPL=-1 .25 DPL=-2.5 DPL=-5 DPL=-7.5 DPL=-12.5 DPL=-17.5 DPL=~25
PR=5 PR=10 PR=20 PR=30 PR=50 PR=70 PR=100

-2.4509 -3.9817 -6.91 12 -9.7517 -15.3003 -20.7505 -28.821
-0.4201 -1 .5041 -3.5859 -5.7471 -10.1993 -14.7493 -21 .6789

-0.25 0.1 194i -0.2869 -0.2648 -0.2593 -0.2553 -0.2537 -0.2526
-0.25 -0.2089 -0.2274 -0.2381 -0.241 9 -0.2451 -0.2465 -0.2475

 

 

 

Table 3.2 Validating the fast sub-system model against Example 2, case 2

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Orig DPL=~1 .25 DPL=-2.5 DPL=-5 DPL=-7.5 DPL=-12.5 DPL=-17.5 DPL=-25
PR=5 PR=10 PR=20 PR=30 PR=50 PR=70 PR=100
-2.3995 -3.9485 -6.8893 -9.7343 -1 5.2872 -20.7396 -28.81 19
-0.64 -1 .5547 -3.61 1 1 -5.7658 -10.2128 -14.7604 -21.6881
-0.25 -0.2303 -0.2484 -0.2498 -0.25 -0.25 -0.25 -0.25
0.25i 0.27231 0.2561 i 0.2514i 0.2506i 0.2502i 0.2501 1 0.2501 i

 

Table 3.3 Validating the fast sub-system model against Example 2, case 3

In addition to the results shown in the tables 3.1 to 3.3, Figure 3.1 shows a

summary of the error in the fast poles locations compared with the DPL; and

shows the value of R, both against the normalized fast poles location (PR);

notice that R = -1—. Every error shown in Figure 3.1 is the average of the errors of
8

the two fast poles with respect to the DPL. The three different cases of the

example are shown. The results of this example and of the next example will be

examined and explained together.

70

 

 

 

 

EO—eror,casel +eror.cas82 +eror.case3 +R,caset +R,cas62 —O-—R,case$]

 

 

 

 

 

 

80% I T. I I 35000
in i i I I
5 70%” : ' """" + 30000
g 60% , g g 3 j R 25000
g 50%4»—-\--—' --------- ----- : —————— I -------
-5 I I I - 20000
3 40°/o 7 I I : m
g . . 15000
; 30% + —:~ —~
:3 10% 4 5000
0% 0
0 20 40 60 80 100
Normalized desired fast poles location (PR)
9_99.ctmt4

 

 

 

Figure 3.1. Effect of fast pole locations on the accuracy of the fast

sub-system model for Example 2

The second example that is used to validate the fast sub-system model is
the fluid pipe system (Example 3) that was introduced in CHAPTER 2. The

homogeneous original system equations for this example was given in

CHAPTER 2 as
E _ _fl

dt 11

£4: .. -2

dt 12

71

dt 01' 11 d: 7? 12 13 I4

dP2 I 2 R22_1&_R{2+2_5:_2]

$213+
dt dt

dP4 (”’5 P1 P2 P3 P4 P4
— = — + R +
dt dt

P5_P1 P2 P3 P4

7771' 12 13 14

And the modified system equations were given as

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0 0 -i O O 0 0 0
It
- - 0 0 0 —i O O 0 O
4‘ 12
qz i O _R+R1+Rs _R+Rs R-i-Rs R+Rs _l: _1_
C1 ll 12 13 I4 15 C
P‘ 0 L _R+Rs _R+R2+Rs R+Rs R+Rs 5 __1_
_4_ P2 C2 11 12 13 14 15 C
R+Rs R+Rs R+R3+Rs R+Rs R 1
dt p3 0 0 __ __ __ _
p4 [1 12 [3 I4 15 C
0 0 R+Rs R+R5 _R+Rs _R+R4+Rs _3 i
p‘ 11 12 I; 14 15 C
_q‘ 0 0 5 _ -5. -5 -5 _1_
ll 12 13 14 Is C
0 o _1_ _1_ __1_ _i __L 0
_ In 12 [3 I4 15 _

The following set of parameters values were used:

C1=0.001 C2=0.001
|1=1000 |2=200 |3=1OO |4=100 l5=1OO
R1 =3000 R2=1000 R3=1 000 R4=1 000 R5=1000

With the above given set of parameters, it can be shown that the original system

eigenvalues are approximately the following:

72

~10.0 ~10.0 ~4.9042

~2.6946 -0.96254 ~0.32758
These slow eigenvalues were obtained by calculating the eigenvalues of the
modified system with two very fast poles (~5012.57 & -4962.41; DPL = 5000) with
respect to the original system. The same procedure described above in this sub-
section that was used with Example 2 was used with this example. Again
equation (3.12) can be used to calculate the values of the parasitic parameters

after obtaining I“, that can be calculated using equation (3.10). Using equation

(3.10), I“, for this example is equal to

 

Table 3.4 shows the calculated eigenvalues from the modified system (as

described above) for every desired pole location (DPL) that was used.

 

 

 

 

 

 

 

 

 

 

 

OriL DPL=~50 DPL=~1 00 DPL=~200 DPL=~300 DPL=~500 DPL=~700 DPL=~1000
PR=5 PR=10 PR=20 PR=30 PR=50 PR=70 PR=100
~103.6 ~166.68 ~286.273 ~401 .678 ~626.409 -846.66 ~1 172.3
~18.46 ~58.1 06 -1 38.77 ~223.405 ~398.692 ~578.447 ~852.8
~10 ~13.039 -10.3186 -10.0675 ~10.0287 ~10.01 -10.005 ~10
~10 ~10 ~10 ~10 ~10 ~10 ~10 ~10
~4.9042 ~4.9169 ~4.9072 ~4.9049 ~4.9045 ~4.9043 ~4.9043 ~4.9
~2.6946 ~2.6946 ~2.6946 ~2.6946 ~2.6946 ~2.6946 ~2.6946 ~2.7
~O.96254 -0.9625 -0.9625 ~0.9625 ~0.9625 ~0.9625 -0.9625 -1
~0.32758 -0.3276 -0.3276 -0.3276 ~0.3276 ~0.3276 -0.3276 ~0.3

 

 

 

 

 

 

 

 

Table 3.4 Validating the fast sub-system model against Example 3

73

 

Figure 3.2 (similar to Figure 3.1) shows a summary of the error in the fast poles
locations compared with the DPL; and shows the value of R both against the

normalized fast poles location.

 

 

 
 
 
 
 
 
 
  

 

 

 

90% 60000
in
g 80% ——————————————————————————————————————— -- ------
'5 F 50000
8 700/04” _-_ ———————
a 600/0 _ __________________ . -. ~ .— — - -------------------- " 40000
c
- 50°/ + ----------------------------------------
8 ° 1 -. 30000 a:
§40701"‘*' —————————— ~: ~~~~~~~~~~
"c' 300/04 ——————————— - -, ———————— . ~~~~~ -: —~20000
‘5 20% 4.-. ............... l ————————— ' l-
: I I I 10000
§10%~~—-~ ~~~~~~ J, ~~~~~~~~ ———————— : —————————
0% l i 'r i 0
0 20 40 60 80 100
Normalized desired fast poles location (PR) a m3 m1

 

 

 

Figure 3.2. Effect of fast pole locations on the accuracy of the fast sub-system

model for Example 3

Examining tables 3.1 to 3.4, it is noticeable that across the range of PR
that is shown in the tables, the slow eigenvalues are very close to the original
system eigenvalues except for one value in case two of Example 2 (13t column in
table 3.2), where PR is equal to five and there is no obvious separation between
the fast and the slow sub-systems. This should not be disappointing because,
the rule of thump is that there should be an order of magnitude separation (this

corresponds to PR=10) so that the fast sub-system would not greatly affect the

74

slow one. On the other hand, it is apparent from tables 3.1 to 3.4 that the
locations of the fast poles are somewhat different from the DPL that was used in
equation (3.12) to calculate the parasitic elements parameters. Figure 3.1 and
Figure 3.2 show that as the ratio PR (called normalized desired fast poles
location in the figures) increases over ten the error rapidly decreases. It can also
be noticed that for all cases shown in both figures, for a DPL that is an order of
magnitude, or more, faster than the original system (this corresponds to PR
bigger than or equal to ten) the following holds true:

1) The error in the fast poles locations are at or less than 55%.

2) The locations of the slow poles are very close to those of the original

system (within 5%).
3.3.2 Conditions of Two Critically Damped Fast Poles

In this sub-section we examine the validity of the conditions, provided by
the model under validation, for having two critically damped fast poles. These
conditions are of particular interest; later it will be argued that a pair of two
critically damped poles is the recommended choice for fast poles. The
characteristic equation of the fast sub-system model was given in equation (3.12)
as

Al +£A+ 1
CI

eq eq

 

= o (3.12)

For the roots of this equation to be critically damped, the following relation must

hold

75

R2 _ 1
2 _
41“, CI“,

 

 

This means that

2

R C = 41 (3.13)
eq

Substituting :— for L in equation (3.12) then, the characteristic equation of
R C qu
the fast sub-system can be written as

4+ 4 _
RC RZCZ

 

112+).

which means that

_:_2
RC

P
Where P (identical to xi) is the location of the two fast critically damped poles as
5 becomes small. This result can be written as

P*RC =-2 (3.14)
Equation (3.14) is a condition provided by the model for obtaining a pair of
critically clamped fast poles. To validate equation (3.14), different values of the
product RC were chosen. For every value of RC, the values of R and C were
changed until a pair of critically damped fast poles (P) were obtained from the
modified system. For these values of R and C for which two critically damped
fast poles occur, the location of the two critically damped poles (P) was recorded

versus its corresponding R and RC. Figure 3.3 summarizes the results of this

attempt for Example 2. The‘horizontal axis in Figure 3.3 is the value of R

76

(R :1); the first vertical axis shows the product P*RC. The second vertical
8

axis shows the ratio between the critically damped fast poles location in the
modified system and the real part of the fastest pole in the original system. This
ratio is referred to as MO. The same procedure described in this sub-section
was implemented for Example 3. Figure 3.4 (similar to Figure 3.3) summarizes
the results of Example 3. Figure 3.3 and Figure 3.4 show that as R increase (i.e.
as .9 decreases) the product P*RC converges to minus two; this result confirms
equation (3.14). They also show the product P* RC is close to minus two (within
%15) for an M0 ratio of ten or more (i.e. for fast poles that are an order of
magnitude or more faster than the original system). It is also obvious that the

product converges to minus two rapidly as MO gets bigger than ten.

Validation results shown in this sub-section and the last sub-section show
that the model under validation (equation (3.11)) would be fairly accurate if the
fast poles are an order of magnitude faster than the original system; and that the
accuracy improves rapidly as the fast poles gets faster than an order of

magnitude.

77

 

 

 

 

 

 

 

 

[+P'RC.C1=0.0107 +P'RC. C1=O.m8 +P'RC.Cl=0.004 +MO.C1=0.0107 +140, 0130008 +MO.C1=0.W1
-1.65 fl. 90
l
4.7.-“ --- -- ”LL- _- ' _______________ ".7 80
l I
P'RC ; 1 —— 70
175 L- - \ ————————————— . —————————————— 4 —————————————————
. 1 M0 1 ~— 60
l
1.8» _- . - - -4 - -- -4
g \ I I 50 g
at .
O. 185 — ~ \ ———————— 1 ~~~~~~ -1 —'3_ 40
l
.\ i I —- 30
-1.9 + ~~~~~ :— ~ ~ - : —————————————————
: : ‘— 20
'1.95 ‘1 i. """" I ““““ : """"""""""""" j... 10
I
l I
0.005+oo 1 00 +04 2.00E+04 3.00E+04
R
W_9_99.h.cf‘ml5

 

 

 

 

  

 

     
       
     

 

 

 

-1.6 7 . 100
-155- .............. L _______ :~ 90
-17__ ”_Praowg. ___-Mo--: _____________ :80
° : : -- 70
0'1-75‘ ““““““ j ********** 460
9: ~18 — : —- 50 g
o. ' :
-1.85 ~ ; :— 40
1 —~ 30
‘1-9”“"“ ““““““““ 5 “““““““““ :_20
-1.95 — 3 ««««« T ---------------- ;_ 10
-2 i l 0
0.00E+00 2.00E+04 4.005+04 6.00E+04
R
ex_13.xls.

 

 

Figure 3.4. Conditions for two critically damped fast poles, Example 3

78

3.3.2.1 Advantages of Critically Damped Fast Poles

The second order fast sub-system should be designed to be critically
damped, to avoid unnecessary computational cost and potential excessive
approximation errors. The reason is that if the fast sub-system is over damped,
then it would contain two real poles; one of them is faster than the other. In this
case the faster pole will determine the stiffness of modified system and, the less
fast pole (of the two fast poles) will determine simulation accuracy.
Consequently, there will be an un-necessary additional computational cost. On
the other hand, if the fast subsystem is excessively under damped, there will be
an overshoot and an oscillatory behavior, which is undesirable. To demonstrate
this we simulate Example 1 that was introduced in CHAPTER 2. Figure 3.5
shows simulation results for calculating the state of the independent variable (it)
from Example 1 using three different systems: The original implicit system, the
modified system with two excessively under damped fast pole, and the modified
system with two critically damped poles. The product RC has the same value in
the used two modified systems. It is obvious that the modified system with the
two excessively under damped fast poles has oscillatory behavior that increases
the approximation error introduced by the modified system. Figure 3.6 show the

simulation results for the state (q) of the dynamic parasitic element (the

capacitor) using the above mentioned two modified systems. It is obvious that
the case with two excessively under damped fast poles has an amplitude that is

much bigger than that of the critically damped case.

79

 

 

 

—— i1, implicit ——i1, under damp fast — i1, crit. damp fast

 

 

 

  

 

 

 

2 . . I

S 1 . --------------- I ------------- : ——————————— :

m o 1 _I 1

g 3 0.8 - ---------- : ———————————— : ————————————
'U-g 06" ******** ; — 1* q ~~~~~~~~~~~~~
5“04~~~ ' -

o. > . : .1 .'

3 0.2 4 ......... ; .............. ; ............. ; ______________
E o ' i 1

0 5 10 15 20
Time, sec “mum,

 

 

Figure 3.5. Approximation error for an excessively under damped and a

critically damped fast sub-systems

 

 

 

 

q, crit. damp fast —q, under damp fasfl

 

J

,,,V.._«

9.0.0
NCO-P-
:11

.0
at
l

- -1--+—-—--—-—-—

State of the dynamic
parasitic element

 

 

 

 

O

_-

O
01

1 O 1 5 20
Time, sec RC_P5.)ds, chart2

 

 

 

Figure 3.6. Behavior of the dynamic parasitic element for an excessively

under damped and a critically damped fast sub-systems

80

3.3.3 Effect of Selected Dependent Element Value on Model Performance

If we have a storage field with a single instance of derivative causality, the
choice of the dependent element is arbitrary. The fast sub-system model as
given in equation (3.11) does not suggest which element one should choose as
the dependent element. Therefore, we would like to answer the following
questions. Suppose that we have a storage field whose elements have different
values. Does the value of the selected dependent element (with respect to the
independent ones) have an effect on the model performance? If there is an
effect, what is the choice for which the model would perform best? This section
addresses these questions through simulating Example 1 that has been

introduced in CHAPTER 2.

Example 1 includes only two coupled inertances; its homogeneous original
system was given in equations (2.27) to (2.28) and its homogeneous modified
system was given in equations (2.29) to (2.31). These equations are rewritten

here again for the reader’s convenience.

fl = - mm£1 Ni2 (2.27)
dt 11 d1

£2- : mp—1 (2.28)
12 11

81

_: —————- 1+— 2"— 2'29
dt 11 p 12” Cq ( )

iq___'zp1_ip2 (2.31)

To explore the effect of the value of the selected dependent element we simulate
Example 1 with different combinations of the values of the inertances 11 and I 2
and with both of them alternatively chosen as the dependent energy storing
element. This is done for two different values of the transformer modulus (m).
The parameter R1 of original system was held constant throughout these
simulations at R1=1000. For every combination of II, 12, and m, we chose
certain desired pole location (DPL) for a pair of critically damped fast poles.
Based on this chosen location equation (3.12) was used to calculate the values
of the parasitic elements R and C. These calculated values of R and C were
then used in the modified model to calculate its eigenvalues. The fast
eigenvalues of the modified model were then compared with the desired pole
location (DPL) that was chosen and used in calculation. Figure 3.7 and Figure
3.8 summarize these simulation results for the two simulated cases, m =1 and
m = 2 respectively. The horizontal axis in both figures is the difference between
the values of II and I 2. The two solid lines are the errors of the fast poles
location for the two cases, when 11 is the independent element and when 12 is
the independent element. Every calculated error is the average of the errors of

the two fast poles with respect to the corresponding DPL. The dashed lines in

82

both figures shows the value (#227 It can be noticed from both figures that
m

choosing either 11 or I 2 as the independent element is not going to affect the

. 12
model accuracy only at the pomt ma. Notice that (II/m2) is the
I

l/m2

equivalent of 11 when reflected on I 2. Away from this special point the model
will be more accurate if the dependent element value (or its equivalent) is smaller

than that of the independent element.

It can be concluded from this sub-section that the dependent element
should be chosen as the element of the least value. Notice that all element

values should be evaluated at the same location.

 

 

 

 

 
 

 

 

 

[+m=1, 11 ind +m=1, 12 ind - o - m"2‘L2/L1,m=1]
0.5 . . 1 I 4.5
‘0 ‘ ' I r 4
2 0.4 .-- ‘ +3.5
0 1
D. g | _ 3 1-
§2°$ ‘ ******* “T """ Qsfi
E 0 .2“ --------- .W‘ . ' <
.3 .9. o | : ' F 1.5 E
z 01« ~1 --.1 ”1-2” *1
In I I . "-o-:.._. 4~ 0.5
l l I l . I.
0 1 1 1 1 0
~2500 -1 500 -500 500 1 500 2500
I1 42 ext _mdel_val.xts. chartt

 

 

 

Figure 3.7. Effect of selecting the dependent element, m=1

83

 

 

[+m=2,11 ind —Il-—m=2, 12 ind - at: - m"2'L2/L1, m=2J

 

 

  

 

 

 

50% 5
E 40% _- -------------- E ------------- --+4
3,357,307., ---- ----- - ------ ~3§
'2 § 20%«1 ------ 1x —————————————————————— —2 9'
'2- , E
2 10%--- .............. 1
'“ s a s
0% 1 1 + 1 1 10
0 2000 4000 6000 8000 10000 12000

I1 -I2

0111 _mdol_vol.xls, chat 2

 

 

 

Figure 3.8. Effect of selecting the dependent element, m=2

3.4 An Algorithm to Improve the Model Performance

Last section showed that the fast sub-system modeI (equation (3.11))
could be used to choose the value of the parasitic elements to locate the fast
poles at certain user-defined location. Although, there would be some error in
the poles locations, this error as shown above is negligible for the design
objective of this work. However, this section shows an algorithm to improve the

model results, in case this is of interest for other uses of the model.

The fast sub-system characteristic equation is rewritten here again for the

reader’s convenience.

84

I = 1
“7 (Sd +J1;Si‘]id)

 

(3.10)

12.13.14. 1
CI

‘9 eq

 

= 0 (3.12)

Before the algorithm is given some of the features of the characteristic equation
(3.12) are demonstrated. For the roots of the characteristic equation (3.12) to

have certain damping ratio (d, the following must hold true

2 R2
4; =——CI ,which means that
d 12 661
eq
2 _ 2
R C_4qu§d (3.15)

Equation (3.15) means that the product R2C will always be a constant as 8

become small. The real part of the roots of the characteristic equation (3.12) is

—R _ -R2C _-R2C 1
21 2RC1,q 21 RC

eq “1

 

 

(3.16)

As RZC is constant, equation (3.16) means that the location of the real part of

the fast poles is inversely proportional to the product RC.

The method that is proposed here for improving the model results is based

on successive improvements of the calculated qu using the calculated fast

eigenvalues from the modified system. Assumes that the user has defined some

specific locations ,1, and A, at which he/she would like to place the fast poles to

separate them from the slow poles. Further assume that the two poles A, and A,

85

have a real part that is referred to as Re(1,), a damping ratio that is referred to

as (d, and an undamped natural frequency that is referred to as com .1- The

algorithm can be described as follows:

1)

2)

3)

4)

5)

Use equation (3.10) to calculate I“), for convenience we will refer to it as

I

eq,u'
Calculate the values of the parasitic parameters RC and Cc using the

following equations

 

 

= 3.17
Icq.uCr wn‘d ( )
1R1 = 24,11)“, ' (3.18)

eq,u

Use the calculated parasitic parameters (RC and Cr) from step two in the

modified system model to calculate the eigenvalues of the modified
system, and consequently identify the two fast poles. Let these two fast

poles be A” and 22., , let their corresponding damping ratio be 5.. and

let their natural frequency be a)”. If 4’, is “equal to” {d go to step six,

otherwise proceed to step four.

 

Calculate an updated equivalent inertia (law) as follows

R
1 = r 3.19
242w... ( ’

Go to step two above.

86

6) The final values of the parasitic elements parameters (R f and C I) should

be calculated from the following

 

Rea”)
R ,C , = Rcc,
R6011)

R {C
——’— f = .1 (3.21)
2 I,“

The loop that includes steps two to five is used to improve the accuracy of the

(3.20)

 

damping ratio of the placed poles and then step six is used to improve the
accuracy of their real part. Equations (3.17) to (3.19) and (3.21) come directly

from equation (3.12). Equation (3.20) is suggested by equation (3.16).

To demonstrate the effectiveness of the algorithm, we apply it to Example
3, the fluid pipe system. The desired fast poles were set to have real part of (-
100) and a damping ratio of 0.7. Table 3.5 shows simulation results for
implementing the algorithm with the loop involving steps two to five, repeated
only five times. The first column in Table 3.5 has five labels, each one is used for
its corresponding row; the next six cells are the original system’s eigenvalues
and the last two cells are the real part and the imaginary part of the desired poles
locations respectively. The following six columns in Table 3.5 show the results
for repeating the loop five times followed by the final result obtained from step six
in the algorithm. It can be noticed from the table that the damping ratio of the
placed fast poles converges monotonically to the desired one as we keep
updating the value of the equivalent inertia. It is also remarkable that

implementing step six has moved the fast poles very close to the two desired

87

values. lt is also noticeable that in all the shown results the slow eigenvalues are

very close to original system’s eigenvalues.

One can repeat the loop involving steps two to five as many more times as
needed; and further more can implement another loop involving steps 2 to 6 until
the desired accuracy of pole placement is reached. However, as the purpose of

this work does not need such accuracy; there will be no further demonstration to

 

 

 

 

 

 

 

 

 

 

 

 

 

 

this idea.
Loop, 1 Loop, 2 Loop, 3 Loop, 4 Loop, 5 Step 6
I «1.- 24.5272 22.6945 21 .6016 20.9271 20.501 1 23.4966
9'. 0.79277 0.75653 0.73541 0.72257 0.71455 0.69554
R 5555.276 4905.165 4536.673 4320.147 4165.213 4669.534
0 1 .76E-06 2.005-06 2.165-06 2.27E-06 2.34E-06 1 .645-06
RC 0.009601 0.009601 0.009601 0.009601 0.009601 0.006605
032756 032756 032756 032756 032758 032756 032756
096254 096253 096253 096253 096253 096253 096253
-2.6946 -2.6946 ~2.6946 ~2.6946 ~2.6946 -2.6946 ~2.6946
-4.9042 -4.9056 ~4.9058 49059 -4.906 -4.906 -4.9056
-10 -10 -10 -10 -10 -10 -10
-10 -10.1335 -10.1517 -10.1643 -10.1729 -10.1766 -10.1377
-100 -112.463 400.772 -94.1667 90.2309 -87.7992 ~96.5376
102.02 66.4619 67.1102 66.7695 66.3261 65.9614 99.7226

 

 

 

 

 

 

 

 

 

Table 3.5. Using the algorithm to Improve model results for Example 3

3.5 Design Procedure for the Fast Sub-System

lt was demonstrated in the previous sections that the fast sub-system

model could be used to place the fast poles at certain desired locations. This

88

section shows some design procedure whose objective is to achieve simulation
efficiency while meeting certain accuracy bounds. As discussed earlier, the
computational efforts needed to simulate the modified system increases as its
stiffness increases. We have also discussed that the fast poles should be faster
than the original system to avoid affecting the original system significantly.
Therefore apparently, there is a need to know the location of the fastest pole in
the original system. But the original system as we know is in the algebraic
differential equation form, and it is not convenient to find exactly or directly it

eigenvalues. This issue will be dealt with in the procedure to follow.

Assume that there is an original system whose modified system is
required to be accurate to within certain specified bound, for some specified
inputs. Assume that the error bounds are set for all the states of the original
system and for the state of the dynamic parasitic element. The proposed
algorithm can be summarized as follows:

1) Set i=0, were j is a counter; choose the smallest inertia as the dependent
inertia. Make an initial estimate of the original system eigenvalues by
omitting the dependent energy-storing element from the system.
Consequently, one can determine an initial estimate of the radius of the

original system’s eigenvalues; let us refer to it as p,.

2) Choose the two fast poles of the modified system to be critically damped;
and choose a desired pole location for the two critically damped fast poles

(DPL) such that

89

3)

4)

5)

6)

DPL = 10 * p, .
Use the fast sub-system model (equations (3.10) and (3.12)) to calculate
the values of the parasitic parameters R and C needed to place the two
fast poles at DPL. lfj > 1 go to step five directly
Use the calculated values for the parasitic parameters in step three in the
modified explicit system and calculate the modified system eigenvalues.
Assume that the real part of the two fast eigenvalues (or the slower one, if

they are over damped) is Re f; and the calculated radius of the slow

eigenvalues is pr. If you can not distinguish two fast poles from the rest

Re I
of the poles, assume that = 1. If

 

 

 

pr
Ref Ref
< 5, or > 20 take
pr pr
DPL = DPL*10*-£"—
Re

I

and go back to step three; otherwise proceed to step five. This will ensure
enough separation between the fast and the slow sub-systems and
prevent simulating an excessively stiff system during the design process
unnecessarily.

Set j=j+1, and DPL]. = DPL.

Simulate the modified system and the original system and calculate the
errors of the states for which the user requires bounds on their errors.

Calculate the ratios between the calculated errors and the user-specified

9O

error bounds; assume that these ratios are r,,r2,...,r,,. Find 'max such

that r

max, j =max{rl,r2,...,rn }. Assume that the error whose ratio is

rm,x is "e" and that its corresponding desired error is e .1 .
7) If 0.8 s 'max 51, go to step nine.
Otherwise;

if j>1, calculate

e.-e

)6 = J 1-1
DPL]. — DPLH '

 

DPL 1,, = DPL}. +

 

This assumes that the error is a linear function of the DPL.

lf i=1, take

 

DPL. . .
DPL: f ’. Where fac=2, 1f rmax >1, and fac=0.5 1f rmax <O.8
ac

8) Go to step three.

9) End of design procedure.
Notice that, when the dependent energy-storing element is omitted from the
original system, the system will become an explicit ODE system whose

eigenvalues can easily be determined.
The procedure presented above assumes that the model of the original

system is conveniently available. This assumption is true for several practical

classes of systems. However, the user might prefer to deal only with the

91

modified explicit system. One possible reason for this preference could be
avoiding deriving the implicit original model. In this case, there would be a need
to find a reference other than the original implicit system to calculate the errors of
the modified system with respect to. Since it is known from the singular
perturbation theory that the errors of the modified system should decrease as 8
decreases, we can choose this reference to be the modified system itself but with
an excessively small 6. An excessively small 6 means an excessively large R.
To verify that certain 5 is small enough for the modified system to be an
“accuracy reference”, it is recommended to ensure that a further decrease (say
an order of magnitude) of 8 would not change the simulation results obtained
from the modified system within the number of digits of accuracy required by the
user. Figure 3.9 shows a flow chart for the main steps of the suggested design

procedure.

492

No Yes
Original system

Find an accuracy Simulate original

reference system

 

 

 

 

 

 

 

 

 

 

 

 

 

4 Estimate rho t

1

Take DPL=10'rho

,l

Calculate R 81 C

l

Simulate modified
system 81 calculate errors

Satisfactory
error?

Figure 3.9. Flow chart for the proposed design procedure

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

     

 

 

 

Change DPL End procedure

 

 

 

 

 

 

 

3.5.1 An Example Using the Design Procedure

To demonstrate the usefulness of this procedure we apply it to design the

fast sub-system for Example 3. The following set of parameters was used in this

example:

93

C1 =0.0001; C2=0.0001;

l1=100; |2=100; l3=100; l4=100; l5=100;

R1 =300; R2=300; R3=1 00; R4=1 00; R5=1 00;
Assume that the error bounds are set to be 1% for all the states of the original
system for unit step inputs. Further assume that the error is defined as the ratio
between the maximum absolute error value and the maximum absolute signal

value.

All the inertia elements here have the same value, which is the most
difficult situation to approximate; we will choose l5 to be the dependent. To
estimate the radius of the slow sub-system, omit l5 and then calculate the
eigenvalues of the original system. This will produce the following set of
eigenvalues

-1.0000

~3.1331

~2.4334 ~ 9.4779i

~2.4334 + 9.4779i

-1.5000 - 9.8869i

-1.5000 + 9.8869i
The fastest pole in this set is (31331). Therefore we will choose DPL=~31.
Then we will use the fast sub-system model (equations (3.10) and (3.12)) to
calculate the values of the parasitic parameters that are needed to obtain two

critically damped fast poles at the chosen DPL (~31 ). This will result in the values

94

C = 5.20296-005 R =1240
Using these two values in the modified model, and subsequently calculating its
eigenvalues produce the following set of eigenvalues for the modified system
-1.0000
~1.0037
-1.0920 - 7.5092i
-1.0920 + 7.5092i
-1 .5000 - 9.8869i
-1.5000 + 9.8869i
~22.5286
~44.2837
It is obvious from this set that Re , = —22.5286, and ,or = —1.5. Therefore,

Re 1
P1

 

z 15, and we can proceed. Next we set DPLl = —31 and i=1; and simulate

the implicit original and the modified systems to calculate the errors. The errors
calculated for the seven states of the original system are, in the order of the

states, as follows:

 

 

 

 

 

 

 

 

 

 

 

DPL=~31
61 0.044539
62 0.044539
63 0.029023
64 0.029023
65 0.035824
66 0.035824
67 0.046535

 

95

Inspecting these errors show that some of them are larger than 1%. Therefore,
there is a need to reduce the DPL used in the design. Since j is equal to one,

DPL will be set to 2*DPL,, which means that DPL2 = —62. Using the new value

for DPL in the fast sub-system model to calculate the values of the parasitic
parameters produces
C =1 .30076-005 R =2480

Using these values of R and C in simulating the modified system, the obtained

errors were as shown below

 

 

 

 

 

 

 

 

 

 

 

DPL=~31 DPL=~62
61 0.044539 0.012109
62 0.044539 0.012109
63 0.029023 0.007981
64 0.029023 0.007981
65 0.035824 0.009484
66 0.035824 0.009484
67 0.046535 0.01237

 

Again there are still some errors that are higher the set bounds. Therefore, there

should be additional iterations as described in the design procedure. The results

of all the iterations are summarized below.

96

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

DPL=~31 DPL=~62 DPL=~64.2 DP L=~68.4

R 1 240 2480 2568 2736

61 0.044539 0.012109 0.011309 0.009988
62 0.044539 0.012109 0.011309 0.009988
63 0.029023 0.007981 0.007452 0.006579
64 0.029023 0.007981 0.007452 0.006579
65 0.035824 0.009484 0.008851 0.007818
66 0.035824 0.009484 0.00851 0.007818
67 0.046535 0.01237 0.011551 0.010199

 

MATLAB ode156 stiff integrator was used to simulate both the implicit and the
modified systems.
seconds and for the modified system whose errors satisfied the errors bounds
was 3.63 seconds. One can notice that to double the value of DPL we needed to

double the value of R. This is not a surprise because equation (3.14) shows that

RC should be set equal to ‘1);152}: and it has also been shown that the product

RRC converges to a constant as 8 gets small. It can be noticed also that the

number of iterations that was needed to achieve the design objectives was not

large.

Simulation time for the original implicit system was 2.96

97

 

CHAPTER 4
ANALYSIS AND DESIGN OF A CLASS OF TWO-TIME-SCALE NONLINEAR
MULTIPORT MODELS

4.1 Introduction

CHAPTER 2 and CHAPTER 3 demonstrated analysis and design of a
class of two-time-scale linear systems. This chapter uses the results and insight
obtained in CHAPTER 2 and CHAPTER 3 and further extends them to cover a
class of nonlinear systems. Tools for analysis and design of a class of nonlinear

systems are presented.

4.2 Class of Nonlinear Systems Under Investigation

The class of nonlinearity under investigation here is the displacement-
modulated nonlinear junction structure. An example for this class of nonlinearity
is a planar mechanism that has linear inertia field with displacement— modulated
transformers. The general form of the homogeneous original system of this class

of nonlinear systems is

X, =J,,(X,)*X, +J,,(X,)*x, (4.1)
X, = J,,(X,)*X, +1,,(X,)*Z,+J,,,(X,)*X,, (4.2)
2, =Jd,(X,)*Z, (4.3)
2, = 5,11, (4.4)
X, = 532,, (4.5)

98

Where, the J matrices are the nonlinear constraints imposed by the junction

structure; X ,. is the independent energy variables vector in the coupled field; X .1
is the dependent energy variables vector in the coupled field; Z, and 2,, are the
coenergy vectors for X,. and X d respectively; and X , is the state vector whose

states are not among the coupled field and its dimension is r. As stated earlier,
one way to recognize the elements of the coupled field is by the appearance of

their parameters in J . The state vector of the corresponding modified explicit

di
system is
”X,“
Xi
X = ,
Xd
_ ‘1 1

 

 

where q is the state of the inserted parasitic compliance. The explicit equations
of the modified system will be in the general form
11, = J,,(X,)* X, + J,,.(X,)* X,. (4.6)

X,. =J,,(X,)*X, +J,,(X,)*s, *X, +

1 (4.7)

J..<X.)*<R*<J..-<X,>*S. *X. -5, *X.)+Eq>
X.=R*<J..<X,>*S.*X.-—S.*X.>+-(1.-q (4.8)
qudi(Xr)*Si*X1_Sd*Xd (49)

Where, C and R are the inserted parasitic compliance and resistance

respectively. Throughout this chapter the junction matrices are functions of the

99

vector X ,; however, sometimes this will not be written explicitly for the sake of

clarity.
4.3 Obtaining the PSSPF

CHAPTER 2 showed how to use some general transformation procedure
to obtain the PSSPF for the modified models of the class of linear systems that
was under investigation. Kokotovic et al. (1986) gives some other general
procedure to transform nonlinear systems to the standard singular perturbation
form. This procedure is analogous to the above-mentioned procedure for linear
systems. Therefore, the effort of applying the nonlinear procedure to get the
PSSPF will directly benefit from the results and insight gained while discussing
how to use the linear procedure. This section shows how to use the general
procedure for nonlinear systems to obtain the PSSPF for the general modified
nonlinear systems in the form (4.6) to (4.9). This procedure can be summarized
as follows:

1) Put the system in the form

51? = h(v,8), ve Rn+m (4.10)
2) The fast variables are n=®(v), where C(v) is the m independent

equations that the n +m equations h(v,0)= 0 reduce to, knowing that
rank %(2 = m .

V

3) The slow variables are x = 0(v). Find 0(v) from

100

60'
—h ,0 =0 4.11
a, (v) 1 1

The assumption is made that
_Qg

rank 35 =n+m (4.12)

5?
Next, this general procedure will be applied to the class of nonlinear modified

systems expressed in equations (4.6) to (4.9). While doing so, the insight gained

from CHAPTER 2 will be used.

We start by doing the scaling (change of variables) Rq for q in the

equations (4.6) to (4.9); the usefulness of this scaling was demonstrated in

CHAPTER 2. Doing this scaling produces

X,=J,,*X,+J,,*X, (4.13)

X,. =J,,*X, +1, *5, *X, +J,, *(R*(J,, *S,*X, -S, *X,)+—R:Rq) (4.14)
a

X, =R*(J,,*S,*X,—S,*X,)+5Rq (4.15)
or

R1): R1,, *S, *X, —RS, *X, (4.16)
After this scaling the system is ready to be put in the form

£V=h(v,£)
Based on what was demonstrated in CHAPTER 2, we will take £=% and

RZC = a, where a is a constant, which will produce the following equations:

101

5X, =aI,*X,+aI,*X, (4.17)

8X, =£l,,*X,+£I,, *S, *X, +J,, *((J,,*S, *X, —S, *Xd)+—l—Rq) (4.18)
a

8X. =0...- *S.— *X. ~S. *X.>+;:-Rq (4.19)

gRq' = 1,, * s, * X,. — s, * X, (4.20)

Following step two in the procedure above we obtain

h,(v,0) = 0 (4.21)
h,(v,0)=J,, *((J,, *S, *X, ——S, *X,)+qu)=0 (4.22)
a
h, (17,0) = (1,, *s, *X, -s, *X,)+qu =0 (4.23)
a
h, (v,0) = 1,, * s, * X, - S, * X, = 0 (4.24)
’11,— _Xr-
h, )1.
Where h = and v = ‘
)1, X,
_h‘IJ _Rq_

 

 

 

 

It is noticeable that either equation (4.22) or equation (4.23) is redundant.

Therefore, h(v,0) = 0 reduces to two independent equations and the vector of the

two fast variables 17 is:

1. X *S.*X.-S *X
”=[ d1( r) 1 1 d d] (4.25)

R9

The rank of 317- is the rank of
V

102

_ . .X. _
313:: I 111151 _Sd 0
O O 0 l

 

 

which has the rank of

— S, 0
0 1 '
This rank is two because S, at 0. The slow variables vector is x = 0(v), where

01

0(v) can be found from equation (4.11). Assuming that o=[g-R-]= -0-"-. ,

l r+1

 

 

r+1' _,

equation (4.11) can be written as

 

 

 

 

 

 

 

 

 

 

 

I 801 801 801 801 '
3X,- 3X, aXd aRq Phr (‘40)-
hi(v,0)
=0.
h d (11,0)
. . . . h m0
a‘7r+i a("Hi a0r+i a‘7r+i _ q( )—
_ 6X, 3X, aXd aRq 1
Examining equations (4.21) to (4.24) shows that
h,(v,0) = 0, and
hi (17,0) = Jidhd(v’0)
Consequently, we can make the following choices
80,, = 80,, -O 80,, ___ 30R :0,

ax ' ax.

r I

8X , aRq

103

80‘, _aJ,,X, 80', I 80', J 60',

 

 

 

 

ax ax, ’ ax,‘ "’ ax, ‘1’ 6114

r

Therefore, the slow variables vector can be written as

X7
0 = (4.26)
X1 —J1'd(Xr)Xd
and
aa 1, 0 0 0
__ = 4.27
av —aJidXd I _ J 0 ( )
1 id
_ 8X, _

 

 

has a rank of r+i.
4.4. A Model for the Fast Sub-System and Proof of Stability

This section will define a model for the fast sub-systems for the class of
nonlinear systems under investigation. To do that, the general form of the two
fast variables that was obtained above will be used to derive a general form for

the dynamics of the fast sub-system. The general form of the two fast variables

was given in equation (4.25) as

J.X *S.*X.—S *X
”=[zl]:[ d1( r) qur d d] (4.25)
2

Therefore, we can obtain

 

d"!=dei(Xr)*S,~*X1+Jd1(Xr)*S1*X1_Sd *Xd (4.28)
dt dt

d1;2 .

-—=R 4.29
dt q ( )

104

Using the general form of the modified system as given in equations (4.13) to

(4.14) we can obtain

£71.: f(X,,X,)+J,,S,J,, *(R*(J,,S,X, —S,X,)+£Rq)-
a

 

dt
R
S,[R*(J,,S,X, —S,X,)+ZRq] (4.30)
“:7: = RJ,,S,X,. — RS,X, (4.31)
Where,
did!
f(X,,X,)=—d;—S,X, +J,,s,*(1,,X, +J,S,X,.) (4.32)

is, in general, a nonlinear function that does not include the parasitic parameter

R. Notice that f(X,,X,) includes the term(

%sxg, which is the part of

the dynamics that results from the dependency of the junction structure on the

state vector X ,. As 8 gets small, R gets large and equation (4.30) can be

approximated as

16171—1: JdiSiJui *(R*(J71151X1 —SdXd)+£Rq)_
a

S,l:R*(J,,S,X, — S,X,)+-:-Rq] . (4.33)
The approximation given in equation (4.33) is an 0(5) approximation. It shows

(X,.)

that the additional dynamics term (%Six") that results from the junction

nonlinearity has a secondary effect that can be neglected for a small enough 5.

105

If the junction structure matrix 1,, were a function of X,., then its derivative with

respect to time would have included X,.. But, equation (4.14) shows that X,

includes the parasitic parameter R as a factor for some of its terms. Therefore, if

J ,,. were a function of X, we would not have been able to neglect the term

55-115}, as an 0(8) approximation. Terms in equation (4.33) can be

dt

rearranged to write the equation as

%=_R*(Sd +Jhrisi'lid)*(JmSiX1'ded)-_R'(Sd +Ji¢TiSiJid>Rq
a

(4.34)
Equations (4.34) and (4.31) are the equations of the fast dynamics; they can be
written as follows

d 17, — R R
= _ _ _ 4.35

ct) cq

d772
= R 4.36
dt 11. ( )

 

l

1,, = T (4.37)
(S, +J,,(X,)*S, *JM(X,))

 

Equations (4.35) to (4.37) represent the fast sub-system model for a small

enough 8.

Examining the fast sub-system model shown in equations (4.35) to (4.37)

shows that the model is a function of the slow variable(s) (X ,). Equation (4.26)

shows that the variables of X , are among the slow variables. Because of the

106

time-scale difference in a two-time-scale system, the slow variables will change
only slightly during the dynamic behavior of the fast variables. Therefore, when
solving for the fast variables that are dependent on the slow variables in a two-
time-scale system, the slow variables are taken as constants at their initial
conditions (Kokotovic et al. 1986). In fact, two-time-scale systems are solved by
solving the fast sub-system as just described, and then solving the slow sub-
system. This is true even for frequency excitation as long as the frequency of
excitation is slower than the fast dynamics. Based on this understanding the

equivalent inertia term 1,, in equation (4.37) needs to be evaluated only once at

the initial conditions of the slow variables X ,. As a result, the fast sub-system

model shown here is the same as the fast sub-system model (equations (3.10)
and (3.11)) that was derived in CHAPTER 3 for the linear case, except that the

matrix Ji d here has to be evaluated at the initial conditions of X , instead of

being totally independent of all variables.

1

 

 

5] 3.10
@HLSM.) “’ ( )
__L _R _
1 R201
4,: 4 “1 (3.11)
R

 

 

Consequently, the bond graph model and the physical interpretations of the fast
sub-system for this class of nonlinear systems are the same as those for the

linear systems studied in CHAPTER 2 and CHAPTER 3.

107

The only nonlinearity that appears in the fast sub-system model is in the

term 1,, due the dependency of 1,, on the slow variables X,. Therefore, if the

fast sub-system is stable as a linear system over the entire range of the values of

1“,, this proves that the fast sub-system is stable (Kokotovic et al. 1986).
Consequently, if we prove that 1,, is always positive, this will prove that the

second order fast sub-system, expressed in equations (4.35) to (4.37), is stable.

As S, and S, are inverse inertia matrices, S, must always be positive definite
and J,f,S,J,, must always be positive semi definite. Consequently, 1,, is always

positive and the fast sub-system is always stable for a small 6. If the fast sub-
system is stable the entire two-time-scale system is stable if the slow (here the
original) system is stable (Kokotovic et al. 1986). This proves that the entire

modified system is stable for a small enough 6 if the original system is stable.
4.5 Effect of Inaccurate Effective Field Value on Model Performance

To use the fast sub-system model (4.35) to (4.37) to design the fast sub-

system, I,, has to be evaluated. As explained above, 1,, can be closely

approximated by evaluating it at the initial conditions of X ,. This sub-section will

investigate the effect of using an erroneous value for qu on the fast sub-system

model performance. To investigate that, we will simulate Example 3 (the fluid

pipe system) that was used in section 3.5. It is a linear system whose qu is a

108

known constant. We will use the linear fast sub-system model (3.12) to place the
fast eigenvalues of the modified system of this example at certain locations while

using some erroneous values for leg. The effect of different percentages of
errors in the value used for qu on the model performance will be evaluated.

The same set of parameters that was used in section 3.5 is used here.

Tables 4.1 to 4.4 summarize the results of simulation carried out in this
section. Every table presents simulation results for a different PR ratio. The PR
ratio was defined previously as the ratio between the desired location of the two
critically damped poles (DPL) that is used in calculations, and the fastest real
part in the original system’ s poles. The first column in every table includes two
labels (in the first two cells, they label their corresponding rows); the DPL value
used in the whole table; and the eigenvalues of the original system. The two

labels are % qu, which refers to the percentage of the correct qu that is used;

and Zeta-fast, which refers to the damping ratio of the two fast poles. Every
column from the remaining columns in the tables shows, on its top, the

percentage of the correct qu used in that column, the damping ratio of the two

fast poles that are obtained from the modified model after using that percentage

of qu’ and the eigenvalues calculated from the modified system after using the
percentage qu that is shown on the top of the column. When the two fast poles

are complex conjugate, the real part is shown in one cell followed by the

109

imaginary part in the following cell. Notice that these situations (complex
conjugate fast poles) will correspond to Zeta-fast (damping ratio of the two fast
poles) that is less than one. It can be noticed from tables 4.1 to 4.4, that in all the
presented cases, the original system's eigenvalues were affected only slightly

(less than 5%). That is true over the entire range of the used % qu (150% to

50%) and over the entire range of the used PR ratio (10 to 40). it can also be

noticed that when the correct qu (100% leg) is used, there are still some errors

in the fast poles locations; these errors decrease as the ratio PR increases.

Using an erroneous value for qu could increase or decrease this error;
depending on whether this error in qu is going to counter affect or add to the

approximation error of the model. The result will be obtaining fast poles that are
(or one of them is, if they are over damped) slower than the desired fast poles
location. Figure 4.1 shows the percentage error in the location of the slower of

the two fast poles vs. the error in the value that is used for leg. When the two

fast poles are under damped, their real part is used in the error calculations and
the mark that identifies their location in Figure 4.1 has no background color.
Figure 4.1 presents the four different cases for the PR ratio. It shows that in

general equation, (3.10) slightly overestimates qu, and that its estimate

improves as 6 gets small (means PR gets bigger). It can be noticed also, that in

all cases for errors in la; within 50%, the additional errors in the pole location

due to using the erroneous value for qu is always less than the percentage error

110

in the value that is used for leg. This suggests that if we are designing the fast

sub-system to set the fast poles no slower than certain location and
approximating the calculation of qu , we should increase the used PR ratio, than
what we would use if we had accurate calculation for qu, with a percentage that
is similar to the expected error in qu. Results shown here show that, in the

worst case, when the inaccuracy in qu adds to the model approximation error,

this can be counter affected by increasing the PR ratio; therefore, the price in the

worst case will be slightly additional computational efforts compared to what

would have been required if qu were known exactly.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Met! 150% 140% 120% 100% 60% 60% 50%
Zeta-fast 1 .340 1 .303 1 .225 1 .142 1 .054 0.959 0.910
DPL=~100 -270.626 -250.171 ~208.902 -166.663 -121.432 -72.276 ~62.212
DPL=~100 -54.275 -54.714 -55.944 ~58.107 -63.271 21.3394i 28.3816i

-4.9 ~4.906 -4.906 -4.907 -4.907 -4.906 -4.909 -4.910

-2.7 ~2.695 -2.695 ~2.695 -2.695 -2.695 -2.695 ~2.695

0327 -0.328 -0.328 0326 0326 0328 0328 0326

-1 -0.963 0963 0963 0963 0963 0963 0963

-10 -10.208 -10.224 -10.263 -10.319 -10.405 -10.554 -10.680

-10 -10 -10 -10 -10 -10 -10 -10
Table 4.1. Effect of erroneous leq on model performance, PR=10

111

 

 

%Ie

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

150% 140% 120% 100% 60% 60% 50%
Zeta-fast 1 .279 1 .239 1 .156 1 .066 0.969 0.860 0.801
DPL=~200 -507.341 -465.297 -379.191 -286.273 472.513 -132.496 ~112.467
DPL=~200 -117.725 -1 19.766 -125.864 -138.770 44.329921 76.545821 64.092541
-4.9 ~4.905 -4.905 -4.905 -4.905 ~4.905 -4.905 -4.906
2.7 ~2.695 -2.695 2.695 2.695 ~2.695 2.695 -2.695
0327 0328 0328 0328 0328 0328 0328 0328
-1 0963 0963 0963 0963 0963 0963 0963
-10 -10.045 ~10.048 -10.056 ~10.068 -10.065 -10.113 -10.137
-10 -10 -10 -10 -10 -10 ~10 -10
Table 4.2. Effect of erroneous leq on model performance, PR=20
% 1
84 150% 140% 120% 100% 80% 60% 50%
Zeta-fast 1 .260 1 .220 1 .135 1.043 0.943 0.831 0.768
DPL=~300 -743.992 ~680.252 548.759 -401 .676 252.537 -192.531 -162.526
DPL=~300 -161.100 ~184.839 -196.328 223.405 89.22632i 129.10441 135.3561
-4.9 -4.904 -4.904 -4.905 ~4.905 -4.905 -4.905 -4.905
2.7 2.695 ~2.695 2.695 2.695 2.695 2.695 ~2.695
0327 0328 0328 0328 0328 0328 0328 0328
-1 0963 0963 0963 0963 0963 0963 0963
-10 -10.o19 -10.021 -10.024 -10.029 ~10.036 ~10.048 -10.058
-10 -10 -10 ~10 -10 ~10 -10 -10
Table 4.3. Effect of erroneous leg on model performance, PR=30
% I
64 150% 140% 120% 100% 80% 60% 50%
Zeta-fast 1 .251 1 .210 1 .125 1 .032 0.930 0.816 0.753
DPL=~400 -980.621 -895.1 52 ~71 8.077 ~514.793 ~332.546 -252.542 -212.540
DPL=~400 244.479 249.946 267.021 210.302 130.9763i 178.7838i 185.925i
~4.9 ~4.904 ~4.904 -4.904 -4.904 ~4.904 ~4.905 ~4.905
~2.7 ~2.695 -2.695 -2.695 ~2.695 ~2.695 ~2.695 ~2.695
0327 0328 0328 0328 0328 0328 0328 0328
-1 0963 0963 0963 0963 0963 0963 0963
~10 ~10.011 ~10.011 -10.013 -10.016 -10.020 ~10.026 ~10.032
-10 -10 -10 -10 ~10 -10 -10 -10
Table 4.4. Effect of erroneous l9q on model performance, PR=40

112

 

 

 

 

 

 

+ PR=10 + PR=20 + PR=30 + PR=40
50% . . , ,

 

 

 

 

 

40%

  
 

30% "

~~~~~~~~~~~

I
l
[2,2
I

20% ~

   

“v-__.__.

I
I
I
1
I
I
l

10%-

 

 

0% 1 1 1
-60% -40% ~20% 0% 20% 40% 60%
% Error in l-eq

 

%Error in slow pole location

“13.3

 

 

 

Figure 4.1. Effect of errors in l-eq on fast poles location errors

4.6 Design Procedure for the Fast Sub-System

This section gives a procedure for designing a fast sub-system, from the
class of nonlinear systems under investigation, to achieve simulation efficiency
while meeting certain accuracy bounds. As discussed earlier, the computational
efforts needed to simulate the modified system increases as its stiffness
increases. The design procedure presented here are very similar to those
presented in CHAPTER 3 except that we linearize the original and the modified
systems to estimate the linearized systems eigenvalues. The linearization is
done by evaluating the junction structure matrices at the initial conditions. As

explained earlier, the slow variables change only slightly during the dynamic

113

 

 

 

behavior of the fast variables. Therefore, for fast sub-system design it would be
sufficient to evaluate its slow-variables-dependent parameters at the initial

conditions of these slow variables.

Assume that there is an original system whose modified system is
required to be accurate to within certain specified bound for some specific inputs.
Assume that the accuracy bounds are specified for all the states whose accuracy
are of concern to the user. The design procedure can be stated as follows:

1) Set i=0, where j is a counter. Linearize the original system around its
initial conditions by evaluating the junction matrices at the initial conditions
of X ,. Choose the smallest inertia as the dependent inertia. If you can
calculate the eigenvalues of the linearized system do so. If you cannot
due to some cumbersome algebra, omit the smallest inertia and then
calculate the eigenvalues for the linearized approximated system.
Consequently, one can determine an estimate of the spectral radius of the

linearized original system’s eigenvalues; let us refer to it as p,.

2) Choose the two fast poles of the linearized modified system to be critically
damped; and choose a desired pole location for the two critically damped
fast poles (DPL) such that

DPL = 10 * p, .

3) Use the fast sub-system model (equations (4.35) to (4.37)) to calculate the

values of the parasitic parameters R and C needed to place the two fast

poles at DPL.

114

4) Set j=j+1, and DPL, = DPL.

5) Simulate the modified nonlinear system and the original nonlinear system
and calculate the errors of the states for which the user requires bounds
on their errors. Calculate the ratios between the calculated errors and the

user-specified error bounds; assume that these ratios are r,,r2,...,r . Find

11

rmax such that = max(r1,r2,...,r,, }. Assume that the error whose

r max, j
ratio is rmax is "e" and that its corresponding desired error is e, .
6) If 0.8 s ’max 5 1, go to step eight.

Otherwise;
if j>1, calculate

e.—e

,6 = J 1‘1
DPL, - DPL,_, ’

 

e,—e

fl

1'

 

DPL,,, = DPL, +

This assumes that the error is a linear function of the DPL.

If j=1, take

 

DPLj . o
DPL = f . Where fac = 2, If rm, >1, and fac =0.5 If rm, <0.8
ac

7) Go to step three.
8) End of design procedure.
If the user prefers to deal only with the modified system, the same method

described in CHAPTER 3 can be used to find an “accuracy reference.”

115

4.7 An Example

The slider crank mechanism shown in Figure 4.2 is taken from Rosenberg
et al. (1983). It has a nonlinear junction structure manifested in its nonlinear

modulated transformer (m ). The equations of the original system of the slider

crank mechanism can be written as

p

1‘3 = i (4.38)
d: J

dp p dp

_Q. = _ le—Q-+—x- (4.39)
p p

J- = m—6 (4.40)
M 1

This means that

Jid=_m’ szfl, X,.:po, Xdsz (4.41)

The transformer modulus m is function of 6 and can be calculated as follows

dx dxdd d6 dx

—=——=m— :

d1 d6 dt d1 0'9

 

x(6)= rcosB+\fl2 --r2 sin2 19

2 smflcosfl (442)
\fi2 - r2 sin2 0

Where 6 is the crank angle, M is the cylinder mass, J is the crank inertia, l is

m =-rsin9—r

the length of the massless rod, and r is the radius of the crank. A bond graph

116

model for the modified system is shown in Figure 4.3; the modified model of this

example can be written as

p
d—9=—9 (4.43)
d: J
dp R +R R 1
__6 :—m2——( 1 p +m—p -m—q (4.44)

dt J 0 M x C

dp R R 1

——l= — —— +— 4.45
dt m] 109 M 17,, C4 ( )

l
iq—=mlp6——p (4.46)

   

\\\\\\\\\\

 

 

 

.16) //7////// “I

Figure 4.2 (a) Schematic diagram for the slider crank mechanism

--T-=Ill~--A MTFI—AII—ARI

I ”’8 I

I M

Figure 4.2 (b) Bond graph model for the slider crank mechanism

117

Figure 4.3. Bond graph model of the modified slider crank mechanism

First we put the modified model in the PSSPF. To do that, start by

changing q to Rq and then follow the above-outlined procedure. Changing q to

 

Rq produces
p
d—g- = i (4.47)
dt J
dp R + R R R
6 —m2I l p +m—p —m—Rq (4.48)
dt J 0 M x a
dp R R
—J—‘-=mflp ——p +—Rq (4.49)
dt 1 9 M x a
dRq R R
— = m— — — 4.50
d J pg M Px ( )

1) Put the system in the form EV = h(V.8)

p
d: J l

118

 

 

 

 

 

Edpg -m2 I8R1+l)p +mip —mqu =h
(11 J 6 M x a 2
dpx 1 l +1R h

6—= — -— — =
dt JPB Mpx a q 3

8dRq__ i _ 1 -h
dt JPB Mpx 4

_h1_ _ 19 _
h
h: I12 V: pg
3 Px
_h4_ _Rq_
2) h(v,0)=0zs

and the fast variables are
1 l

7): m7p0 _fipx (4.51)
R4

119

0 3111 m —1
31: [p737 7 H O],hasranktwo.
0 0 O 1
3) F1nd the slow variables x: 0(v ), %—0 h(v ,0)=0

P 80, Ba, 80', 80',
86 8p, 8p, BRq
Bo, Ba, 80', 802
_ 86 3p, 8p, €9qu

 

h(v,0) = 0, leads to

 

 

 

53—0—1 (v,0)+a—0’-‘—h (v,O)—h(vaU,-10) 39—111

0 =0
ape 2 apx 8R Rq 4,(V )

80’ 30’ 30 Ba
392 hi (V10) '1' '51—);2'h2 (V10)+ $12113 (V10) '1' 3R; h4 (V10) = 0

 

 

To obtain the slow variables vector

I 9 I
0': ,
176+me

make the following selections

 

 

 

 

 

Ba_ Ba, _ 0 80', _ fl =
a—Hz ’ 3179 ’ 8p. ’ 3R9
a_0'2 _ am 80'2 :1 802 =m Ba,
36 319 p,, 3p, ’ 3p, ’ aRq
Consequently,
I 0 0 0
3;, = am , has rank two.
V
_ 1 0
39 P; m

120

 

‘l.

(4.52)

Assume that the example has the following set of parameters
71
J=2, M=10, R1 :1 00, 6(1 = 0) = 4

r=0.2, l=3*r

Where, all units are SI units. Further assume that the rest of the initial conditions
are zeros and that, the input torque is

T=2,t$2; T=0 otherwise.

Assume that error bounds are specified for both the speeds that are coupled in
the inertia fields in the original system; and whose coupling is broken in the
modified system by inserting the parasitic elements. Those two speeds are the
angular velocity and the piston velocity. Assume that the errors are specified
such that the ratio between the maximum absolute value of error and the
maximum absolute value of the signal should be less than or equal to 1%. We

will call the error in angular velocity 61, and the error of piston velocity 62.

The value of the transformer modulus at time zero is
m(%) = 0.17577

This is the value that will be used for the fast sub-system design. Substituting
this value in the original model and calculating the eigenvalues of the associated
linearized model shows that the fastest calculated eigenvalue is about -1.34.
Therefore we will start our design iterations by choosing DPL=-13.4. Using this
value to design the fast sub-system in the linearized modified model, gives the

following values for the parasitic parameters

121

 

R = 232.14 C = 642956-004
Using these two values in the linearized modified model gives the following
eigenvalues for the linearized modified model

~15.2594

-11.7448

-1.3406
Proceeding to step five in the algorithm to simulate the modified and the original
systems (both nonlinear) and calculate the errors, we find that

e1 = 0.0090

62 = 0.0457
As the 62 is still larger than the desired bound (0.01) and j=1, we will take DPL
equal to twice its current value and go to step three. The results of this iteration

and the rest of the iterations are summarized in table 4.5.

 

 

 

 

 

 

 

 

 

 

DPL ~13.4 ~26.8 ~36.9 ~48.4 -60.3
R 232.1401 464.2801 639.25 838.47 1044.6
61 0.009 0.007 0.0067 0.0066 .0065
6 2 0.0457 0.0234 0.0164 0.0129 .010

 

 

Table 4.5. Summary of results of design procedure

It is noticeable that the algorithm converges quickly. It can also be realized that
62, which is the error in the speed of the dependent inertia is larger than that of
the independent inertia and also more responsive and consistent to changes in
6. Analysis presented in this work did not investigate this issue. Figure 4.4 and

Figure 4.5 compare simulation results of both the original and the modified

122

systems resulting from the previous design. Comparisons are shown for the two
velocities that were coupled in the original system and whose coupling was
broken in the modified system. The modified system is seen to follow the original
system closely. It is interesting also that simulation time for the original implicit
system was 5.49 seconds while it was 2.86 seconds for the modified stiff system
with the parasitic parameters values corresponding to DPL=-60.3; same
integrator (MATLAB ode159) was used in both cases. This indicates that
following the design procedure prescribed here can lead to computational

savings even as compared with the effort required for the original implicit system.

 

 

 

 

— omega, modified — omega, original

 

 

......................

—————————

22--..

 

 

 

 

0.15 — 1
0.05 —§ —
0 5 10 15 20
Time, sec

setf 1_ex3 attach 1.1116, chant

 

 

 

Figure 4.4. Comparing omega from the original and modified systems

123

 

 

 

 

 

 

 

 

 

 

 

 

 

—velocity, modified —velocity, original
0 . I
§ : :
E ~0.02 '
E 004 ~ : E- E
0 l 1 1
.9 . : :
g 006-6 ~~~~~
5 : : :
E008~~ ~~~~~~~~ ~~~~~ ------------

O. I : :
‘0.1 1 1 I

0 5 10 15 20

Time, sec
ex31_sett 1 .146. 6116112

 

 

 

Figure 4.5. Comparing velocities from the original and modified systems

124

CHAPTER 5
CONCLUSIONS

5.1 Summary of Main Contributions

This research work had two main objectives:

1) Analyze the class of two-time-scale modified systems that results from
using parasitic elements in a specified structure to eliminate derivative
causality. The purpose of this analysis is to better understand the
important two-time-scale properties of that class of two-time-scale
systems.

2) Formalize the above-mentioned approach of using parasitic elements to
eliminate derivative causality in a design context. The purpose of this
formalization was to enable the user to achieve the greatest possible
simulation efficiency for specified approximation accuracy, by making the
appropriate choices for the values of the parasitic element parameters.

These two main objectives have been met. Analysis shown in this document
revealed some important aspects of the time-scale properties of the modified
problem under investigation, which greatly improved our understanding of the
systems under investigation. The design procedure presented in this document

serves to meet the second objective.

Several contributions have been made toward achieving these objectives.

The following gives a brief description and comments on these contributions:

125

 

1)

2)

3)

4)

5)

6)

Determining the choices of 6 and variables scaling (E=%, RZC is

constant, and Rq for q) that help formulate the modified systems in the

standard singular perturbation form and reveal its time-scale properties.
Formulating the modified system in a physically meaningful standard
singular perturbation form (PSSPF). The states of this PSSPF are related
directly to the states of the original system and other physically meaningful
quantities that are of particular interest in the study of these systems.
Using the general procedures outlined in Kokotovic et al. (1986) could
lead to a standard singular perturbation form whose state variables are not
physically meaningful and not of special interest. Methods to
systematically obtain this PSSPF were shown for classes of both linear
and nonlinear systems.
Showing that the modified system can be divided into two distinct sub-
systems and that each sub-system can have its own model. This is an
enabler for any analysis and design work on the modified system, and
offers great insight into it.
Showing that the fast sub-system can be closely approximated by an R-L-
C model.
Showing that as 6 gets small the only parameters from the original
system that continue to affect the fast sub-system are the parameters of
the coupled field; all other effects from the original system diminish.
Showing that the transient of the error introduced into the algebraic

equation can be directly designed, which is a significant improvement on

126

what was known of taking R “very big” and C “very small” to obtain fast
parasitic dynamics in prior published work.

7) Giving procedures to design the fast subsystem to achieve certain
required accuracy at a highly efficient computational price. The design
procedures used the results obtained from the analysis work. The key
part of the procedures is the use of the fast sub-system model. In the
case of linear models it was shown that the user could also place the fast
poles as desired. It was suggested that the designer should target the two
fast poles to be critically damped to provide for efficient computation.
Examples of using the design procedures were given. These examples
show that the proposed design procedures can effectively and efficiently
select the appropriate values for the parasitic parameters, so that the
modified system is just stiff enough to meet the accuracy bounds. They

also showed that the typical number of iterations is small.

The results of these research efforts make the method of derivative
causality elimination that was under investigation more understandable, more
tractable, and much more rigorous as an option for eliminating the derivative
causality. In addition they bring a new set of tools and a very useful mind-set for

attacking certain classes of problems.

127

5.2 Suggested Directions for Future Research

Research results that were presented herein point to some directions of
further research efforts. This section briefly describes two such research

directions.

5.2.1 More than One Dependent Storage Element

Research efforts presented in this document addressed directly the case
in which there is one dependent energy-storing element in the coupled field.
However, there are situations in which there is more than one energy-storing
element in a coupled field. In these situations, it the modified system is designed
to be two-time-scale, then there will be couplings between the different parasitic
fast dynamics. This observation is made based on some preliminary work done
by the author. The nature, the effect, and the chances of manipulating and/or
approximating these couplings are to be evaluated. Another possible solution to
this case would be to design a multi-time-scale system. The number of time-
scales would be equal to the number of the dependent energy storing elements
plus one, for the original system. The process could be a successive design
process in which the designer starts with the fastest time- scale sub-system and
then does the second fastest sub-system, and so on. This solution might raise

computational-efficiency-related concerns, especially as the number of time

128

scales increases. The proof of stability of the fast sub-systems also must be

addressed.

5.2.2 Analysis and Design for Eliminating Algebraic Loops

It was indicated in CHAPTER 1 that the presence of an algebraic coupling
in a resistance field could also lead to an implicit set of differential equations.
Algebraic loops could also occur in the junction structure. As explained earlier it
is often desirable to obtain an explicit set of ODEs. One way to break these
algebraic loops is to insert parasitic dynamic elements. As a result some
parasitic dynamics occurs. It seems that this parasitic dynamics would share
some similar aspects with the parasitic dynamics that was dealt with in the
problem that was under investigation in these research efforts. Consequently, it
is recommended to examine the possibility of treating this problem with the same

techniques that were used in the derivative storage element.

129

REFERENCES

130

REFERENCES

Allen, R. R., "Dynamics of Mechanisms and Machine Systems in Accelerating
Reference Frames," Journal of Dynamic Systems, Measurement, and Control,
Vol. 103, pp. 395-403, Dec. 1981.

Auslander, D., Kempf, C., Mechatronics Mechanical Systems Interface, Prentice
Hall, 1996.

Bennett, B.S., Simulation Fundamentals, Prentice Hall, 1995.

Borutzky, W., Cellier, F. E., “Tearing Algebraic Loops in Bond Graphs,“ Trans. of
SCS, 13(2), 1996a, pp.102-115.

Borutzky, W., Cellier, F. E., "Tearing in Bond Graphs with Dependent Storage
Elements," Proc. CESA 96 IMACS Multiconference, Lille, France, July 9-12,
1996b, pp. 1113-1119.

Bos, A., Tiernego, M., "Formula Manipulation in the Bond Graph Modeling and
Simulation of Large Mechanical Systems,“ Journal of The Franklin Institute, Vol.
319, No. 1/2, pp. 51-65, February 1985.

Cellier, F. E., “Hierarchical Non-Linear Bond Graphs: A Unified Methodology for
Modeling Complex Physical Systems,“ Simulation, Vol. 58, No. 4, pp. 230-248,
April 1992.

Cellier, F. E., "Bond Graphs: The Right Choice for Educating Students in
Modeling Continuous Physical Systems," Simulation, Vol. 64, No. 3, pp. 154-
159, March, 1995.

Cellier, F. E., Continuous System Modeling, Springer-Verlag, 1991.

Chang, D. C., and Rohde, S.M., "The Role of Synthesis, Analysis, and
Simulation in Engineering 6 Complex System: the Automotive Vehicle,"
International Congress of Mechanical Engineers, 1996, Atlanta.

131

Dauphin-Tanguy, G., and Borne, P., “Order Reduction of Multi-Time Scale
Systems Using Bond Graphs, the Reciprocal System and the Singular
Perturbation Method," J. Franklin lnst., Vol. 319, No. 1/2, pp. 157-171, 1985.

Edstrom, K., “Generating State Space Equations from A Bond Graph with
Derivative Causality Using Singular Perturbation Theory," ICBGM, pp. 105-109,
1999.

Gordon, 3., Liu, S.,"A Singular Perturbation Approach for Modeling Differential-
Algebraic Systems," ASME Journal of Dynamic systems, Measurement, and
Control, Vol. 120, 1998, pp. 541-545

Granda, J., "Bond Graph Modeling Solutions of Algebraic Loops and Differential
Causality in Mechanical and Electrical Systems,“ Proc. Applied Simulation and
Modeling, IASTED Conference, SF,CA,1984, pp.183-193.

Joseph, B. J., Martens, H. R., "The Method of Relaxed Causality in The Bond
Graph Analysis of Nonlinear Systems,“ Journal of Dynamic Systems,
Measurements, and Control, March 1974, pp. 95-99.

Karnopp, D and Margolis, 0., “Analysis and Simulation of Planar Mechanism
Systems Using Bond Graphs,” Journal of Mechaical Design, Vol. 101,April 1979,
pp.187-191.

Karnopp, D. C., Margolis, D. L., Rosenberg, R. C., System Dynamics: Modeling
and Simulation of Mechatronic Systems, -John Wiley 8 Sons, 2000.

Karnopp. 0, "Alternative Bond Graph Causal Patterns and Equation
Formulations for Dynamic Systems,“ Transactions of the ASME, Vol. 105, June
1983, pp. 58-63.

Karn0pp, D., “An Approach to Derivative Causality in Bond Graph Models of
Mechanical Systems," J. Franklin Institute 329 (1), 1992, pp. 65-75.

Karnopp, D.C., “Lagrange’s Equations for Complex Systems,“ Trans. ASME, J.
Dynamic Systems, Measurement, and Control, December 1977, pp. 300-306.

132

Karnopp. D.C., “Power-Conserving Transformations: Physical Interpretations and
Applications using Bond Graphs,“ J. Franklin Institute 288 (3), Sep. 1969, pp.
175-201.

KamOpp, DC, ”The Energetic Structure of Multibody Dynamic Systems," J.
Franklin Institute, V. 306, No. 2, 1978, pp. 165-181.

Khalil, Hassan K., Nonlinear Systems, New York MacMillan 1992.

Kokotovic, V., Khalil, Hassan K., O’Malley, R. E., Singular Perturbation Methods
in Control. Analysis and Design, Academic Press, London, 1986.

Kokotovic, V., O’Malley, R. E., Sannuti, P., “Singular Perturbation and Order
Reduction in Control Theory- An Overview,“ Automatica, Vol. 12, pp. 123-132,
1976.

Margolis, 0., "Bond Graph Causality and the Derivation of Analytical Transfer
Functions,“ ICBGM, 1993.

Ogata, K., Modem Control Engineering, Prentice Hall, 1990.

Paynter, H. M., 1961, Analysis and Design of Engineering Systems, MIT Press,
Cambridge, MA.

R. C. Rosenberg, ”Exploiting Bond Graph Causality in Physical System Models,“
ASME Journal of Dynamic Systems, Measurements, and Control, Vol. 109, pp.
378-383, 1987.

Rosenberg, R. 0., "State Space Formulation for Bond Graph Models of Multiport
Systems," Transactions ASME Journal Dynamic Systems, Measurements and
Control, Vol. 93 Series G, March 1971, pp. 35-40.

Rosenberg, R. C., The Enport Reference Manual. 1987, Rosencode Associates
Inc., Lansing, MI.

133

Rosenberg, RC. and D. C. Karn0pp, Introduction to Physical System Dynamics,
McGraw-Hill Book Co., Inc., New York, 1983.

Rosenberg, RC. and J. McCalla, “Power to the User in Formulating Model
Equations,” ASME 1991, DSC-Vol. 34, Automated Modeling.

Rosenberg, R.C., “Exploiting Bond Graph Causality in Physical System Models,“
ASME Journal of Dynamic Systems, Measurements, and Control, Vol. 109, pp.
378-383, 1987.

Rosenberg, RC, and Andry, A.N., Jr., “Solvability of Bond Graph Junction
Structures with Loops,” IEEE Trans. Circuits Syst., Vol. CAS-26, No. 2, pp. 130-
137, 1979.

Rosenberg, RC, and Zhou, T., "Power-Based Model Insight,” Proceedings of
the 1988 ASME Winter Annual Meeting, Symposium on Automated Modeling for
Design, pp. 61-67.

Stein, Jeffrey L., Yih-Tun Tseng, "Strategies for Automating the Modeling
Process," Proceedings of the 1991 ASME Winter Annual Meeting, Symposium
on Automated Modeling, pp. 69-87.

Sueur, 0., "Bond Graph Approach to Multi-time Scale Systems Analysis,“ J.
Franklin lnst., Vol. 329, No. 5/6, pp. 1005-1026, 1991.

Sueur, C., and Dauphin-Tanguy, G., "Bond-graph Approach for Structural
Analysis of MIMO Linear Systems," J. Franklin lnst., Vol. 328, No. 1, pp. 55-70,
1991.

Sueur, C., and Dauphin-Tanguy, 6., "Structural Controllability/Observability of
Linear Systems Represented by Bond Graphs," J. Franklin lnst., Vol. 326, No. 6,
pp. 869-883, 1989.

Tomkinson, D., and Home, J., Mechatronics Engineering, McGraw-Hill Book Co.,
Inc., New York, 1996.

134

 

Van Dijk, Johannes, 0n the Role of Bond Graph Causality in Modeling
Mechatronic Systems, Ph. D. dissertation, University of Twente, Netherlands,
1994.

Zeid, A, "An Algorithm for Eliminating Derivative Causality,“ Proc. ASME Winter
Annual Meeting, Vol. 8, Automated Modeling for Design, 1988, pp. 15-21.

Zeid, A. And Rosenberg, R., “Eigenvalue Spectra and Bounds for Certain
Classes of Dynamic Systems Having Tree Bond Graphs," IEEE Transactions on
Circuits and Systems, Vol. 33, No. 12, Dec. 1986, pp. 1232-1240.

135

 

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

1111111111111111111111111211111