(Asif Mahmood Mughal) Real Time Modeling, Simulati (B-Ok - Xyz)
(Asif Mahmood Mughal) Real Time Modeling, Simulati (B-Ok - Xyz)
(Asif Mahmood Mughal) Real Time Modeling, Simulati (B-Ok - Xyz)
Modeling and simulation of dynamical systems has been taught all around the globe
with many methods, including standard modeling techniques of differential equa-
tions. These differential equation models are later transformed into frequency
domain for analysis. In the latter part of the last century, time domain analysis
began to be taught at undergraduate and graduate levels. Many books are currently
serving the purpose of introducing modeling and simulation of linear time invariant
systems. In the current scenario, new degrees are evolving such as systems engi-
neering, mechatronics engineering, and biomedical engineering, where concepts of
different conventional degrees (electrical, mechanical, etc.) merge together. Stu-
dents are always curious in taking elective courses which widen their depth of
knowledge. Accordingly, the courses and relevant text are being introduced con-
tinuously for meeting this demand. A modeling technique that covers the interdis-
ciplinary modeling problems in a much better way than any other technique is the
“bond graph method.” This scheme has been taught since the mid-1970s, and the
best available text is System Dynamics: Modeling, Simulation, and Control of
Mechatronic Systems by Dean C. Karnopp, Donald L. Margolis, and Ronald
C. Rosenberg, now in its fifth edition. Most of the literature or lectures are drawn
from this book, as all respectable authors themselves were students of Prof. H. M.
Paynter, who introduced bond graph techniques in the late 1960s. Recent develop-
ments, especially software in modeling and toolboxes in MATLAB as well as
Maple, allow students to use different modeling paradigms interchangeably. In
addition, a modeling software 20-Sim originally used for bond graph now also
provides other modeling techniques as well. A simple course only on bond graph
modeling may interest students from different backgrounds. The purpose of this
book is to introduce modeling, simulation, and control perspectives to mixed-
background students of different undergraduate degrees at final (senior) year or as
a first/elementary course at graduate level. This book introduces state space model-
ing derived from two modeling techniques, namely, Lagrangian formulation and
bond graph method, followed by analysis of system and usage of model in control
system design. This is a one-semester text derived into eight chapters for lectures
vii
viii Preface
and other texts as case studies in different applications. I myself taught this course
five times at the graduate level and improved the text accordingly in order to meet
the requirements of students from various engineering disciplines and concentration
areas. This book in no capacity is a replacement of System Dynamics, which is
always the unique reference book for students and teachers in terms of examples,
exercises, and descriptions of problems. I myself prefer that students should attempt
the examples and problems of this book for better understanding and practice to get
a good grip on the subject.
The book presented here, Real Time Modeling, Simulation and Control of
Dynamical Systems, provides a better and comprehensive study text to cover over
one semester. The focus is to introduce state space modeling and Lagrangian and
bond graph techniques such that students can work on advanced problems easily
after this course. The book also introduces MATLAB commands and most of the
examples and problems are developed in 20-Sim software. Analysis and simulation
chapters use MATLAB and Simulink in order to broaden students’ knowledge of
various software utilities. Real-time modeling means covering the system in dif-
ferential equations or equivalently state space formulation. The analysis of the state
space model is presented in time domain with concepts of eigenvalues and gains
mainly. Frequency-based analysis is briefly discussed. The input-based responses
are also simulated in time domain with state space and not with frequency domain.
The present control techniques are well understood with pole placement rather than
frequency domain-based analysis; hence the focus of time-based responses are
given as compared to frequency-based responses, which are mostly popular in
current textbooks. The chapters are organized systematically in the natural way
so students keep on building from previous knowledge. There are no recommended
prerequisites for engineering students as they already cover engineering mathemat-
ics in detail before taking these electives. However, universities that offer graduate-
level engineering programs to non-engineering graduates should check for engi-
neering mathematics knowledge before recommending this course or make an
“advanced engineering mathematics” course as a prerequisite to it. The chapter
descriptions are given as follows.
Chapter 1: The first chapter is very important in this text in order to explain
details of any system, its dynamics, what a model is, and why to model. This
chapter discusses in detail the model types and classifications of models in engi-
neering practices. Furthermore, state space formulation is introduced from ordinary
differential equations for nonlinear as well as linear time varying and invariant
systems. Concepts of linearization and transfer function realizations are also
discussed. The state space realization developed in this chapter will further be
used in Lagrangian and bond graph techniques in the following chapters.
Chapter 2: Lagrangian modeling is discussed in this chapter from basic to
advance level. A single-point rigid body model in Cartesian coordinates is initially
discussed and converted into polar coordinates. Later the same model is expanded
in vector notation for any generalized coordinates (Cartesian, spherical, cylindrical,
etc.). The basic physics behind work-energy relationship is discussed with
Preface ix
conservative and nonconservative forces. At the end of the chapter, state space
formulation is obtained and elaborated with a few examples.
Chapter 3: This chapter introduces preliminary concepts of bond graph model-
ing techniques by drawing analogies of different energy systems. Concepts of
power and energy variables, word bond graphs, causality, and flow of power are
introduced with arguments building toward elements of bond graphs.
Chapter 4: Elements of bond graphs, which are 1-port, 2-port, and junctions, are
discussed with their causality assignments and the governing or constitutive laws
and examples of different energy systems. Examples within different energy sys-
tems with corresponding analogy are described. Junctions are discussed with signs
of power flow as well as effort/flow equations. Modulated transformer and gyrator
(2-port elements) are discussed in text, and modulated sources (1-port) are intro-
duced in the problem section. It is assumed that students are well versed in basic
knowledge of physics for electrical, mechanical, translational, and rotational
systems.
Chapter 5: Writing analytical equations in state space formulation using con-
stitutive laws is explained with different examples of electrical, mechanical, and
mechatronics systems. As per existing literature, a generalized rule for all kinds of
systems is not available. A simplified procedure for all kinds of systems is to be
evolved into a basic methodology of bond graph as independent energy systems
(electrical, mechanical, hydraulics, etc.) must not be violated. This procedure is
emphasized and discussed to develop bond graph of multidisciplinary energy
systems. This chapter develops capability in students to fully model different
systems in bond graph and write its equations.
Chapter 6: Different concepts of advance bond graph modeling are summarized
and combined in this chapter in a logical sequence so the reader can understand why
specific advancement or extra elements are needed. It covers algebraic loops,
derivative causality, and fields, as well as some advanced mixed system with partial
bond graph portions. A further combination of elements in different requirements is
discussed, and, at the end, the topic of vector bond graph is touched. This chapter
helps to develop a significant grasp on the bond graphs method by exploring
different key concepts so students can utilize them in their projects. This is the
last of four chapters covering specific bond graph techniques.
Chapter 7: Analysis and simulation of state space systems are discussed with
details in time domain and some analysis in frequency domain. Most engineering
students are familiar with responses in transfer function approach, and in this
chapter, these responses are discussed with state space method. Generalized solu-
tions of state space equations followed by free and forced responses are discussed as
homogenous and non-homogenous solutions. Concepts of internal and external
stabilities are explained with free and forced responses, thus discussing damped
and amplified systems. Forced responses are studied with impulse, step, sinusoidal,
and decaying sinusoidal signals in order to elaborate the bounded input and
bounded output of control system designs with stability planes and resonance.
Chapter 8: A complete paradigm of control system design based upon state
space model is discussed in this chapter. It introduces a detailed simulation diagram
x Preface
for control system design followed by its simplification based upon design require-
ments. Open- and closed-loop observers are also discussed with complete compen-
sator design. Performance parameters and tuning of PID and control gains are also
discussed with the help of MATLAB tuning and constraint optimization blocks,
which help to understand the requirement of control design. The basic idea is to
familiarize the readers with different types of control systems, keeping the design
scope for advance courses.
Chapter 9: This chapter consists of different projects and research work for
multidisciplinary energy systems with special focus on bioengineering systems.
Recent advances in bond graph allows the research in the area of biomedical and
bio-mechatronics systems. Bond graph models of biomechanics systems are
discussed briefly in order to introduce their applications in this field. The selected
projects of this part can also be taught in the classroom for description, analysis, or
assessment purposes depending upon the time of the semester.
The book is written for a semester workload in real-time modeling, simulation,
and control system design. This can serve as a beginning course for control system
specialization and may be taken with linear system theory during the same semes-
ter. I prefer and teach this course with a semester project requiring students to
develop and simulate a model of any electromechanical or mechatronics system by
utilizing bond graph technique. It is highly recommended that the book System
Dynamics should be employed as a reference book and guideline for the semester
project from students. The lecture-wise breakdown is given in the table below for a
16-week semester for a 3-credit-hour course.
I am grateful to my Allah (( )ﺃﷲthe Lord of the heavens and earth) for His countless
blessings in my life. I begin with Salat-o-Salam (peace and blessing) upon His final
messenger and prophet Muhammad (ﺣﻀﺮﺕ ﻣﺤﻣﺩ ﺻﻠﻰ ﺍﷲ ﻋﻠﻴﻪ ﻮ ﻋﻠﻰ ﺍﻠﻪ ﻮ ﺻﺣﺑﻪ ﻭ ﺍﻣﺘﻪ
)ﻮﺒﺎﺮﻚ ﻮﺳﻠﻢ ﺗﺴﻟﻳﻤﺎﱟ ﻜﺜﻴﺮﺍaccepting that I am not capable of achieving any objective but
with the blessings of Allah ()ﺃﷲ, the Great and the Majestic.
I dedicated this book to Mr. Qaiser Yasin Khan, who was my supervisor at work
in my early days of professional life. He became a true source of inspiration for me
during these years due to his excellent mental acumen, professional integrity, and
visionary guidance. This dedication is a small appreciation for a very few qualities
I tried to imitate from his legendary professionalism.
I am very thankful to Prof. Dr. Afzaal Malik, who introduced me to bond graph
modeling during a graduate course which developed my interest in this area and
lead to several publications. I am also thankful to D. C. Karnopp, D. L. Margolis,
and R. C. Rosenberg, who responded to my queries and even sent me scanned
explanations of handwritten and drawn solutions. I am highly indebted to my Ph.D.
supervisor Dr. Kamran Iqbal from the Dept. of Systems Engineering, University of
Arkansas at Little Rock (UALR). He taught me several subjects of controls and let
me introduce bond graph modeling to his undergraduate course during my Ph.D.
work. Any words of appreciation will be less for his contributions to enhance my
skills, and he is always motivating me for new endeavors. I also like to especially
acknowledge Tallal Saeed, Madiha Zoheb, Tayyaba Qaiser, Sana Javed, Maryam
Iqbal, Hassan Javed, Muhammad Furqan, Mona Jaffar, and Mariam Javed, whose
work became a significant contribution in this text. I am greatly thankful to Marta
Moldvai for her precious suggestions and improvement of this manuscript and for
providing me great encouragement and support during the review and publication
process. I am beholden to my family that they all contributed in their capacity
xi
xii Acknowledgments
xiii
xiv Contents
Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
Chapter 1
Dynamical Systems and Modeling
A system is an isolated part of the Universe that is of interest due to a specific reason
for an appropriate time. This isolated reality can be used for discussion, study,
analysis, improvement, betterment, protection, or any other objective. The isolated
part in which we take interest actually becomes the system and everything else is
surroundings, environment, or the rest of the Universe. For example, an astrophys-
icist wants to study the gravitational forces of the sun and its planets, so this isolated
interest becomes a solar system for him. The rest of the galaxies, and even relatively
small events on the earth and planets, become surroundings for him. An oceanog-
rapher analyzing the waves and currents in the Pacific is only interested in his
isolated reality and the rest of the happenings in the world are collectively an
environment. Similarly, an engineer working on a comparatively small problem
whether designing, studying, or discussing is a system for him and everything else
in the world becomes an environment or surroundings to him till he stops working.
A car is a system, an isolated part of reality for a buyer. But this car itself is a
combination of many other systems such as a combustion system, braking system,
fuel injection system, electrical system, steering systems, etc. A buyer at a time is
interested in all including the color of car, financial value in the market, fuel
consumption, etc. He is overall looking for an appropriate driving system that
suits his needs. But several different types of engineers are interested in several
subsystems of the same car. This car in turn is a collection of subsystems. A
subsystem is a complete system when one develops an interest in it and other
subsystems of the same car become an environment or surroundings for it. This
example is similar to the human body system with several subsystems like the
respiratory system, digestive system, reproductive system, central nervous system,
etc. A body needs all subsystems at a time but an orthopedic surgeon usually only
focuses on the skeletal system, and a gastroenterologist is only interested in the
digestive system.
Systems engineering approaches clearly distinguish between what the matter of
interest is and what is not important in the current time. This leads to focus on
achieving required objectives clearly and thoroughly. When boundaries are drawn,
then improvement of things inside boundaries can easily be targeted. The Systems
Engineering approach leads to quantify the objectives and their solutions with
well-established analytical methodologies. Modern day bioengineering is an exam-
ple as it is solving a gap between myriads of experimental data and well-developed
analytical theories. Generally the systems engineering approach quantifies the
system’s dynamics through mathematical representation in an appropriate time
frame. A dynamical system is one that changes its state with time. In reality there
is not a single system that exists in the entire universe that does not change its state.
The system that appears static is either not considered within an appropriate time
frame, or the change is out of the scope of interest. The North Star has been there for
centuries and is responsible for helping with navigation. But the rotational axis of
the earth changes during its movement in an elliptical path and about every 5000
years, the North Star passes his duties to another star for the same purpose. The
appropriate time frame has to be defined properly when system dynamics are
studied. Now let’s consider an example of a book in a local library, which has
never been issued and never been read. So in a library system the status of this book
may be constant as it does not come into a scope of interest. But changes are there as
pages are turning pale and binding is getting loose. For defining any dynamical
system, two things are very critical: one is an appropriate time frame to notice its
change, and the second is the state that is being changed. As the Universe moves
and changes its state in space in every moment, so there is only one thing that is
persistent in that Universe, and that is the occurrence of change. For an engineer, a
system that changes its state (points of interests) with a significant rate of change is
a dynamical system. This window of time can be a few thousand years for an
astrophysicist and a few hundred years for the geologist who is considering the
formations of canyons and gorges. Similarly, a financial engineer plans the econ-
omy for few years to monthly basis. A mechanical engineer is normally interested
in minute scale and a communications engineer worries about milliseconds of
delay. This discussion and examples of dynamical systems can go longer and longer
as this approach is applicable to each and every subsystem of the Universe at all
levels of time scales. But to define an isolated part of reality, it is important to
choose the right size time scale and the right set of state variables.
State variables are a set of defined variables that are required to completely describe
a system for its purpose. These are variables of interests and the minimum set of
these variables that completely define the system are assembled together for the
description of the system. The combined set of these variables, whether minimum
or redundant, is called state vector. For a control systems engineer, the time-
*
dependent state vector is represented as x ðtÞ and each independent variable is
defined as x followed by a numeral x1(t), x2(t), etc.
1.3 Modeling 3
Each system interacts with its surroundings, universe, or its environment through
the input or the output. Whatever is getting in by a system from the environment is
the input and whatever is passing out to the environment is the output. Input
describes the changes of a system posed from its surroundings whether these are
big, small, desirable, or inappropriate. Whatever a system is contributing to its
surroundings, whether asked for or not, is the output of this system and usually an
input to some other system. A car takes air from its environment for combustion and
this is an input to the car and it runs on the road at a desirable speed. The desirable
speed is an output of the car, along with the gases passed from its exhaust. Those
gases are also the input to the ecological system of the city providing air to the car.
The desired input that is also a function of time is usually denoted by u(t) and a
*
vector u ðtÞ for multiple inputs. Similarly required outputs or objectives are denoted
*
by a vector y ðtÞ.
1.3 Modeling
gravitational pull from the stars to the car, he thus makes a very complex and
unnecessary model.
Causality: The system may be dependent upon different inputs or requirements that
must be considered appropriately. A shuttle in space can never rely upon air supply
for combustion. Similarly, an engineer cannot provide a water-cooled model for
cooling of machinery in a desert. All dependencies must be considered properly,
including time-dependent scenarios such as linking of events, as one cannot occur
until the other starts.
Redundancy: This refers to variables or states that can be avoided because some
other variable is doing the similar and required work. But reliable systems have
some redundancy because if one fails to deliver, the other may help to achieve the
objective. It is the same as driving a car in a city where fuel is readily available or
driving on a long route highway where fuel may or may not be available easily. A
cautious approach is to fill the tank whenever you get a chance for uncertain travel
conditions. Redundant systems are generally more reliable but they require more
cost to operate or consume more energy. An engineer should not always avoid
redundancies but tradeoff between cost and reliability for the better solution.
Redundancy also exists in natural systems where required, such as two eyes, two
ears, two nostril paths, two lungs, but one tongue.
Experimental Validation: If a model behaves like a real system as conceived,
designed, and analyzed through experimental trails then this model is more suitable
than any other model that is untested. But usually experimental validation is not
possible for all systems and is often a costly requirement for developing prototypes
followed by the actual product delivery. A car model can be tested on roads
rigorously before delivery to customers, but in the fashion design industry a
model may wear new trendy clothes and catwalk.
• From the intrinsic evidence of this creation, the Great Architect of the Universe
now begins to appear as a pure mathematician (Sir James Jeans)
• This universe is created by an Entity who knows all the Mathematics behind
it. Who knows how to solve, when, and where to solve the right variable at right
time by the right method (Author of this book—Asif Mahmood Mughal)
Mathematical modeling plays a thorough and pivotal rule in understanding and
exploring natural phenomenon. A Navier-Stroke’s equation can be applied to
determine where the wind will take a falling leaf or understand the flow of blood
in the human arteries. Tidal waves in oceans, light rays from sun, northern lights
(aurora) or power transmission through different medium, all can be explained with
electromagnetic and wave equations. A golf ball shot straight to the hole, an off
swing delivery to clean bold, a hawk’s dive to catch a fish, an airplane flying to its
destination are a few examples of balancing aerodynamics forces. Pade’ approxi-
mation is used to estimate delays in communication channels as well as delays in a
feedback loop in the central nervous system of the human body. Mathematical
models are applied as financial engineering tools to predict slumps in markets
which act as forcing functions for small companies to either quit or merge into
large business. The supply and demand of daily needs, upcoming trends in people
purchases, and data mining to boost sales are examples of economical modeling of
globalization. When a system is modeled in the form of mathematical equations, it
will be easy to understand, interpret, or change accordingly for a variety of
perspectives. There are different kinds or classifications of mathematical models.
Lumped System Model: An ordinary differential equation with time as the depen-
dent variable is often known as the lumped system model; it changes its state with
time or it may change itself with time. A particular definition of lumped system is a
system with only one variable of change (dependant variable) and a finite number of
states (number of elements in state vector).
Distributed Model: A model with two or more independent variables of change is a
distributed model. It is a system that may change with time as well as with spatial
coordinates. A good example of this system is weather, which is always changing
with time and location. A distributed model occurs with more than one independent
variable and the partial differential equation is required to represent this system. A
system with infinite number of states is also known as distributed system, e.g., a
sample time delay in a system which may have infinite real numbers to fill between
two consecutive samples.
Autonomous System Model: A dependent variable does not appear explicitly in the
modeling equation and only appears in solutions. So variables can be defined as
xðtÞ ¼ x and in general it is time invariant system. Most time invariant models are
autonomous such as the application of Newton’s second law of motion on ordinary
objects, F ¼ m €xðtÞ ¼ m a. Where acceleration and applied force are functions of
time, but variable of time t does not appear when values or equations for force or
acceleration are quantified.
6 1 Dynamical Systems and Modeling
Time Varying or Time Invariant System: A model that changes with time is a time
varying model; it means the mathematical equation that is representing the system
also changes with time. So for a one specific instant a modeling equation is
applicable but after some time another modeling equation is governing the same
system. Large air carriers, space shuttles, and other high fuel consuming objects are
considered time varying because of their rapid change of mass during their journey.
Newton’s second law of motion is represented as FðtÞ ¼ mðtÞ €xðtÞ for these time
varying systems. But if we treat mðtÞ ¼ m, constant mass during its operation then
this is a time invariant system, though force and acceleration are changing with
time. A system which is considered unchanged in physical properties for some
specific time is known as the time invariant system.
Linear and Nonlinear Systems: If an input excites a system to perform some actions
under specific conditions and this can produce the same results for the same input
and conditions in every other time, then it is a linear system. In simple words, the
principle of superposition should apply to the system. Let us consider a system y
ðtÞ ¼ xðtÞ þ uðtÞ where two different results from two different values and inputs
are given as y1 ðtÞ ¼ x1 ðtÞ þ u1 ðtÞ and y2 ðtÞ ¼ x2 ðtÞ þ u2 ðtÞ and if their homogenous
and additive combination can be obtained for x3 ðtÞ ¼ α x1 ðtÞ þ β x2 ðtÞ and u3 ðtÞ
¼ α u1 ðtÞ þ β u2 ðtÞ as y3 ðtÞ ¼ α y1 ðtÞ þ β y2 ðtÞ, then the system is linear. If this
principle cannot be applied to a system then that system is nonlinear. An example is
if either of x(t) or u(t) are trigonometric functions or other algebraic functions where
superposition is not applicable.
Continuous or Discrete Systems: If a system is defined at each and every instant of
time, so that its solution can be obtained for each instant of time from the smallest
fraction of a microsecond to several thousands of years, then it is a continuous
system. If in order to define a system or to obtain its solution, we have to sample the
time instants according to our convenience, then this is a discrete system. A discrete
system is silent between two intervals and choosing the sampling interval as
smallest as possible makes it closest to a continuous system. A car runs on road
continuously but a microcontroller-based electronic fuel injection system of a car
processes the combustion in discrete samples.
Causal and Non-causal System: A system that only depends upon current and
previous values to give output is a causal system. A causal system holds a memory
for some past events in accumulation with current scenario to produce results.
Causal systems do not depend upon the future or predict or anticipate anything
for output, unlike a non-causal system that does depend upon future inputs to have
current output results. All physical systems are causal in nature such as mass
springs, dampers, capacitors, and inductors which hold some previous stretches or
values in addition to applied forces to give present output. Non-causal systems can
only exist in virtual worlds and are used to solve computer algorithms or estimation
problems.
1.6 Mathematical Model 7
All physical systems in the Universe are causal, continuous, nonlinear, and time
varying (and distributed) in reality. So why study Acausal (non-causal), discrete,
linear, and time-invariant lumped systems? The only answer is ease to handle the
problem. The laws of nature do not require future inputs to solve for current
problems, as causal systems use only previous history and current inputs to solve
for present output and future responses with some valid assumptions. But we
require developing an algorithm which predicts future states and satisfy the current
requirements, so a computer program may be non-causal. If we like to solve an
exact solution of a small problem, it is possible to do so but for large problems
continuous systems require memory storage and computation speed that may not be
available to us. Thus it requires to discretize the system with such close sampling
interval that it seems to behave like a continuous system. Any continuous system
problem solver in computer is actually a discrete problem solver with very small
sampling time. The principal of superposition does not hold in nature in an exact
sense; whatever we get today may not be equal to what we get tomorrow in
actuality. Somewhere some depreciation or addition has been done. But we study
linear systems because superposition principal makes us to solve analytical models
with great ease. We can easily assume that under certain conditions the solution of
linear and nonlinear will be so close that their relative error is meaningless to
us. Similarly a dynamical system is the one that changes with time but to monitor
change in each and every parameter, make the model with almost infinite states
very difficult to solve analytically. Anything that comes on the axis of time is bound
to change but monitoring change of a state is relatively very easy as to monitor
change in physical properties. A book may be losing its shine or binding but if never
read or issued is in a constant condition for a librarian. Changing values of resistors
on a minute scale due to aging and environment may be useless when designing
a circuit. A car’s motion on the road is represented by the time invariant system
F ¼ m a where m is constant for the car, though there is continuous consumption
of fuel and neglecting variation of total mass by assuming mðtÞ m for maximum
fuel and loading capacity. This assumption is not valid for a rocket taking a shuttle
to the space.
:::
z ðtÞ þ aðtÞ €zðtÞ þ bðtÞ z_ ðtÞ þ cðtÞ zðtÞ þ d ðtÞ δðtÞ ¼ 0 ð1:1Þ
In this equation a(t), b(t), c(t) and d(t) are time varying physical parameters, and it is
a linear time varying model of a system. If we consider these parameters as constant
a, b, c and d then this becomes a linear time invariant (LTI) system. This is a causal
system and requires initial conditions to solve for its solution. It is difficult to solve
higher-order different ODEs and solution of first-order ODE is simple. It is better
and convenient to solve three ODEs than to solve a single third-order ODE. The
Eq. (1.1) can be represented as three first-order ODEs. Let x1 ðtÞ ¼ zðtÞ,
x2 ðtÞ ¼ z_ ðtÞ ¼ x_ 1 ðtÞ, and x3 ðtÞ ¼ €zðtÞ ¼ x_ 2 ðtÞ and Eq. (1.1) is now represented in
matrix form for the three first-order ODEs
2 3 2 3 2 3 2 3
x_ 1 ðtÞ 0 1 0 x1 ðtÞ 0
4 x_ 2 ðtÞ 5 ¼ 4 0 0 1 5 4 x2 ðtÞ 5 þ 4 0 5 δðtÞ ð1:2Þ
x_ 3 ðtÞ c b a x3 ðtÞ d
The third equation in Eq. (1.2) is same as Eq. (1.1) with a change of variable from
z(t) to x3(t), but in order to solve for first order x3(t), it is required to solve x1(t) and
x2(t). The solution of a single first-order differential equation and three first-order
equations simultaneously is the same algebraically. Methods of linear algebra
provide solution to a matrix differential equation (1.2).
Let us define two outputs as z_ ðtÞ þ k δðtÞ and eðtÞ zðtÞ þ hðtÞ €zðtÞ. So two outputs
are in a vector form and with definition of state variables are given as
* y 1 ðt Þ x2 ðtÞ þ k uðtÞ
y ðt Þ ¼ ¼
y 2 ðt Þ eðtÞ 2
x1 ðtÞ þ3hðtÞ x3 ðtÞ
x1 ðt Þ ð1:4Þ
0 1 0 k
¼ 4 x 2 ðt Þ 5 þ uðtÞ
eðtÞ 0 hðtÞ 0
x 3 ðt Þ
1.7 State Space Method 9
The Eqs. (1.2) and (1.4) represent the model from input to output as state space
representation using state, input, and output vectors in Eqs. (1.3) and (1.4). These
equations can be represented generally for any linear system as
_
* * *
x ðtÞ ¼ AðtÞ x ðtÞ þ BðtÞ u ðtÞ ð1:5Þ
* * *
y ðtÞ ¼ CðtÞ x ðtÞ þ DðtÞ u ðtÞ
_
* * *
x ðtÞ ¼ f x ðtÞ, u ðtÞ, t
ð1:6Þ
* * *
y ðtÞ ¼ g x ðtÞ, u ðtÞ, t
Both linear and nonlinear systems have first-order derivatives for all state variables.
If state and input vectors can be separated through superposition then the system is
linear as in Eq. (1.5), otherwise they are nonlinear. The time variance depends upon
the existence of variable t explicitly and state equation and output equation in
Eqs. (1.5) and (1.6) are functions of time. The system is LTI if matrices in Eq. (1.5)
are constant.
_
* * * *
x ðtÞ ¼ f x ðtÞ, u ðtÞ
ð1:8Þ
* * * *
y ðtÞ ¼ g x ðtÞ, u ðtÞ
For LTI system the order of matrices A, B, C, D depends upon the number of states,
* *
inputs, and outputs. The state vector x ðtÞ is n 1 for n state variables, u ðtÞ is a
*
p 1 vector for p inputs, and y ðtÞ is q 1 vector for q outputs of the system. This
definition of these vectors determines the order of matrices as:
The matrix A is always square and the order of matrix D represents the direct input–
*
output relation. In Eq. (1.5) A is 3 3 matrix, and only one input so u ðtÞ is 1 1 and
*
so B is 3 1. For two outputs, y ðtÞ is 2 1 vector and C is 2 3 matrix and finally
for two outputs and single input D is 2 1. For nonlinear systems there are only
n 1 state equations. The vector field F is either real space ℝ or complex space ℂ
depending upon the system; and for all physical systems, these matrices belong to
10 1 Dynamical Systems and Modeling
real space ℝ. These matrices may belong to complex space ℂ for a virtual system or
computational purposes, but cannot be physically implemented through real
devices.
1.8 Linearization
The nonlinear state space system is linearized to obtain linear state space equations.
The mathematical principle behind the concept is Taylor’s series approximation for
* *
defining any function within a specified range. Let x e and u e be equilibrium state
and input at equilibrium point respectively and let Δx, Δu be the small perturbations
from the equilibrium where the system is linear. The nonlinear functions in n 1
vector for state space with n state variables and p input variables is given as
_
* * * *
x ðtÞ ¼ f x ðtÞ, u ðtÞ ð1:9Þ
* * *
As we perturb the state x within small region Δx from the equilibrium x e at given
*
inputs u e such as
* * *
x ðtÞ ¼ x e þ Δx ðtÞ ð1:12Þ
* * *
u ðtÞ ¼ u e þ Δu ðtÞ ð1:13Þ
_
* * * * * *
x ðtÞ ¼ f x e þ Δx ðtÞ, u e þ Δu ðtÞ ð1:14Þ
_
* * * * * *
x ðtÞ ¼ f x e þ Δx ðtÞ, u e þ Δu ðtÞ ð1:15Þ
1.8 Linearization 11
* *
dx ðtÞ * * * ∂f * * * ∂f * * *
*
ffi f xe ; u e þ * xe ; u e Δx þ * xe ; u e Δu
dt ∂x ∂u
þ higher order terms ð1:16Þ
* * *
We know that f xe ; u e ¼ 0 and we can ignore second-order and higher-order
*
terms by assuming negligible effect within the perturbation region Δx ðtÞ and
*
Δu ðtÞ. We also observe that the equilibrium point is a constant so
* * *
*_ d
*
x ð tÞ d x e þ Δx ðt Þ d Δx ð t Þ *
dΔx *_
x ðtÞ ¼ ¼ ¼ ¼ ¼ Δx : ð1:17Þ
dt dt dt dt
Where state matrix A and input gain matrix B are computed Jacobian matrices at
equilibrium points
2 3
∂f 1 ∂f 1
∂f 1
* ∂x1 ∂x2 ∂xn
∂f * * 6 ∂f 2 ∂f 2 ∂f 2 7
6 7
A ¼ * xe ; u e ¼ 6 ∂x1 ∂x2 ∂xn 7 ð1:19Þ
∂x 4⋮ ⋮ ⋱ ⋮ 5
∂f n ∂f n ∂f n * *
∂x1 ∂x2
∂xn
x ¼ xe
* *
u ¼ ue
2 3
∂f 1 ∂f 1
∂f 1
∂up
* 6 ∂u ∂u2 7
∂f * * 6
1
∂f 2 ∂f 2 ∂f 2 7
B ¼ * xe ; u e ¼ 6 ∂up 7 ð1:20Þ
6 ∂u1 ∂u2
7
∂u 4⋮ ⋮ ⋱ ⋮ 5
∂f n ∂f n ∂f n * *
∂u1 ∂u2
∂up
x ¼ xe
* *
u ¼ ue
It can also be linearized at the same equilibrium point using Taylor series expansion
in a similar way. The Output-state matrix C and output–input gain matrix D are also
computed as Jacobian matrices as follows
12 1 Dynamical Systems and Modeling
2 3
∂y1 ∂y1
∂y1
∂x1 ∂x2 ∂xn
∂y * * 6
* ∂y2 ∂y2 ∂y2 7
6 7
C¼ xe ; u e ¼ 6 ∂x1 ∂x2 ∂xn 7 ð1:22Þ
4⋮ ⋱ ⋮ 5
*
∂x ⋮
∂yq ∂yq ∂y *
∂xqn x ¼ xe
*
∂x1 ∂x2 *
*
u ¼ ue
2 3
∂y1 ∂y1 ∂y1
∂u
6 ∂u ∂u2 7
∂y * * 6
1 p
* ∂y2 ∂y2 ∂y2 7
D ¼ * xe ; u e ¼ 6 ∂up 7 ð1:23Þ
6 ∂u1 ∂u2
7
∂x 4⋮ ⋮ ⋱ ⋮ 5
∂yq ∂yq ∂yq *
∂u1 ∂u2
∂up * x ¼ xe
* *
u ¼ ue
The relation between outputs to input is called transfer function and mathematically
it is a rational function G(s) of a single variable. This single variable expresses the
mutual relationship of different states and also the relationship from input to states
and then to states to output. Transfer functions are representatives of a system in
frequency domain in the form of algebraic equations. The mathematical relation-
ship in frequency domain between states X(s), input U(s), and output Y(s) is given as
There is a single transfer function for single input to single output (SISO) system;
but for multiple inputs multiple outputs (MIMO) system, transfer function G(s) is a
matrix of q p order. Each entry in the matrix transfer function represents the
relationship of a specific output to a specific input. There is a relationship between
transfer function and state space matrices, which is known as realization of a system
in frequency or time domain. State space is realization in time domain and taking
the Laplace transform of Eq. (1.5) for zero initial conditions provides realization in
frequency domain as G(s). A system with zero initial condition is known as a
relaxed system. Taking Laplace transform of Eq. (1.5) as LTI yields
1.9 Transfer Function 13
For relaxed system xðt0 Þ ¼ xð0Þ ¼ 0 and so transfer function (TF) G(s) is
*
G ðsÞ ¼ C ½s I A1 B þ D ð1:26Þ
For the system in Eqs. (1.2) and (1.4), there are two outputs and an input, so there
are two transfer functions in a matrix. The first and second rows of C matrix
correspond to Y1(s) and Y2(s) respectively.
2 3 2 31 2 3
Y 1 ðsÞ
s 1 0 0
* 6 U ðsÞ 7 0 1 0 6 7 6 7 k
G ðsÞ ¼ 6 7¼
4 Y 2 ðsÞ 5 40 s 1 5 4 0 5 þ
e ðt Þ 0 hðtÞ 0
c b a d
U ðsÞ
ð1:27Þ
The order of numerator m is always less or equal to the order of denominator n for a
proper rational function. If m ¼ n then the corresponding entry in D is non zero and
if m < n then corresponding entry in D matrix is 0. In Eq. (1.27), the numerator of
Y 1 ðsÞ Y 2 ðsÞ
U ðsÞ has order equal to denominator which is 3, but for U ðsÞ the order of numerator
should be less than 3. The order of denominator is always equal to the order of
square matrix A and determinant of ½s I A is denominator for all transfer
*
functions in G ðsÞ.
Example 1.1: State Space Equations Formulate a state space representation for the
following set of ODE and obtain transfer function matrix.
€θ þ 2θ_ þ 3θ 2x ¼ τ þ F
τ
2 ð1:29Þ
x_ 5θ ¼ F
5
14 1 Dynamical Systems and Modeling
Solution The set of state variables consists of θ, θ_: , and x. The total order of the
system is three by summing highest derivative terms in all variables. Let us define a
* *
state and input vectors x and u as
x1 ¼ θ u1 ¼ τ
x2 ¼ θ_ ¼ x_ 1 u2 ¼ F
ð1:30Þ
x3 ¼ x
This is state space form given as in Eq. (1.7) as LTI system with
2 3
2 3 0 0
0 1 0 6 1 1 7
A ¼ 4 3 2 2 5 and B ¼ 6 4 275. The output in this equation is not
5 0 0 1
1
5
explicitly mentioned and so in this case, we can measure all the state variables at
output. The output equation will be
* *
y ðtÞ ¼ I 33 x ðtÞ ð1:33Þ
2 3 2 3
1 0 0 0 0
In this case, C ¼ 4 0 1 0 5 and D ¼ 4 0 0 5 in output equation. Matrix D is
0 0 1 0 0
often a null matrix and null matrices are not written in state space representations.
MATLAB command ss2tf gets a TF from state space representation and vice
versa for tf2ss command. Using symbolic math, the transfer function can directly
be computed as
syms s
A¼[0 1 0;-2 -2 2;5 0 0]
B¼[0 0;1 1/2;-1/5 1]
C¼eye(3,3)
G¼C*inv(s*eye(3,3)-A)*B
simplify(G)
1.9 Transfer Function 15
*
This generates 2 3 G ðsÞ matrix of transfer function for three outputs and two
inputs. All denominators in TF are same, whereas there are six different numerators
2 3
*
s 0:4 0:5s þ 2
1
G ðsÞ ¼ 3 4 s2 0:4 0:5s2 þ 2 5 ð1:34Þ
s þ 2s2 þ 2s 10
0:2s2 0:4s þ 4:6 s þ 2s þ 4:5
2
Example 1.2: State Space and Transfer Function Formulate a state space repre-
sentation and transfer function matrix of
:::
z ðtÞ 2€zðtÞ 4z_ ðtÞ þ 8zðtÞ ¼ f 1 ðtÞ þ 8f 2 ðtÞ; ð1:35Þ
x1 ðtÞ ¼ zðtÞ
u1 ðtÞ ¼ f 1 ðtÞ
x2 ðtÞ ¼ x_ 1 ðtÞ ¼ z_ ðtÞ yðtÞ ¼ mðtÞ ð1:36Þ
u2 ðtÞ ¼ f 2 ðtÞ
x3 ðtÞ ¼ x_ 2 ðtÞ ¼ €zðtÞ
Note: If there is nonzero term in D then the corresponding numerator has the same
order as denominator
Example 1.3: Nonlinear State Space Equations Formulate a state space represen-
tation for the following nonlinear system
where output is ηðtÞ ¼ ðy_ ðtÞÞ2 þ sin ðpðtÞÞ cos ðzðtÞÞ. Later find the transfer
*
function of the system by linearizing at xe ¼ ½ π=4 1 4 0 T and
*
ue ¼ ½ 0 0 T .
Solution The state and input vector is
x1 ðtÞ ¼ yðtÞ
x2 ðtÞ ¼ x_ 1 ðtÞ ¼ y_ ðtÞ u1 ðtÞ ¼ uðtÞ
ð1:40Þ
x3 ðtÞ ¼ zðtÞ u2 ðtÞ ¼ uðtÞ
x4 ðtÞ ¼ x_ 3 ðtÞ ¼ z_ ðtÞ
The state space representation for four state variables and two inputs is given as
* *
x_ 1 ðtÞ ¼ f 1 x ðtÞ, u ðtÞ, t ¼ 2 ∙ cos ðx1 ðtÞÞ 2x3 ðtÞ
* *
x_ 2 ðtÞ ¼ f 2 x ðtÞ, u ðtÞ, t ¼ 3 cos ðx1 ðtÞÞ 3ðx2 ðtÞÞ2 þ sin ðu1 ðtÞÞ þ 1
pffiffiffiffiffiffiffiffiffiffiffiffiffi ð1:41Þ
* *
x_ 3 ðtÞ ¼ f 3 x ðtÞ, u ðtÞ, t ¼ 2 sin x1 ðtÞ ðx3 ðtÞÞ
* *
x_ 4 ðtÞ ¼ f 4 x ðtÞ, u ðtÞ, t ¼ u1 ðtÞ u2 ðtÞ
syms x1 x2 x3 x4 u1 u2
f1¼2*cos(x1)-2*x3
f2¼3*cos(x1)-3*x2*x2+sin(u1)+1
f3¼2*sin(x1)-sqrt(x3)
f4¼u1-u2
a11¼diff(f1,x1)
b12¼diff(f1,u2)
c13¼diff(g1,x3)
C¼[c11 c12 c13 c14]
D¼[d11 d12]
x1¼pi/4; x2¼1; x3¼4; x4¼0; u1¼0; u2¼0;
eval(C)
Problems
P1.1 Input–output
Explain the input–output of different subsystems of a car, how is the input of
one system derived from the output of the other system?
P1.2 Redundant Systems
Redundancies increase reliability of the system but also increases costs. Are
expensive systems always reliable? Explain your answer with examples,
models, and system.
P1.3 State Space of Linear Systems
Solve to find state space representation and TF of the following systems
2 €yðtÞ 4 y_ ðtÞ þ 6 yðtÞ þ 2 φ_ ðtÞ ¼ 8 δðtÞ ρðtÞ
(a) € ðtÞ 3 yðtÞ ¼ δðtÞ
φ
hðtÞ ¼ 5 φðtÞ þ yðtÞ þ 2 δðtÞ
(b) €I ðtÞ þ RL I_ ðtÞ þ LC
1
I ðtÞ ¼ EðtÞ
€yðtÞ þ ω yðtÞ ¼ FðtÞ
2
(c)
lðtÞ ¼ yðtÞ
P1.4 Linearization for State Space Representation
Find state space representation of nonlinear system and linearize at relaxed
conditions
(a) l θ€ðtÞ þ g sin θ ¼ f ðtÞ
(b) m v_ ðtÞ ¼ Fg b v2 and output is h(t) where vðtÞ ¼ h_ ðtÞ
y_ 1 ðtÞ ¼ a y1 ðtÞ b y1 ðtÞ y2 ðtÞ
(c)
y_ 2 ðtÞ ¼ k y1 ðtÞ y2 ðtÞ l y2 ðtÞ
References
1. Friedland, Bernard. 2005. Control System Design—An Introduction to State Space Methods.
Mineola: Dover Publications.
2. Chen, Chi-Tsonga. 1999. Linear System Theory and Design, 3rd ed. Oxford: Oxford University
Press.
3. Heidi, Christian, André Ran, and Freak van Schrage. 2007. Introduction to Mathematical
Systems Theory—Linear Systems, Identification and Control. Basel: Birkhäuser Verlag.
4. Rugh, Wilson J. 1996. Linear System Theory. Upper Saddle River: Prentice Hall.
5. Franklin, Gene F., J. David Powell, and Michael L. Workman. 1997. Digital Control of
Dynamic Systems, 3rd ed. Richmond: Addison-Wesley-Longman.
6. Gershenfeld, Neil. 2011. The Nature of Mathematical Modeling. New York: Cambridge Uni-
versity Press.
7. Ogata, Katisuhiko. 1995. Discrete—Time Control Systems, 2nd ed. Upper Saddle River: Pren-
tice Hall.
8. Hespanha, Joao. 2009. Linear Systems Theory. Princeton University Press: Princeton.
Chapter 2
Lagrangian Modeling
The basic laws of physics are used to model every system whether it is electrical,
mechanical, hydraulic, or any other energy domain. In mechanics, Newton’s laws
of motion provide key concepts to model-related physical phenomenon. The
Lagrangian formulation of modeling derives from the basic work–energy principle
and Newton’s laws of motion. The basic law states that the force acting on a body is
directly proportional to its acceleration equated with constant mass.
dv d2 x
F ¼ ma ¼ m ¼m 2 ð2:1Þ
dt dt
dx_ dx
€xdx ¼ dx ¼ d x_ ¼ x_ dx_ ð2:4Þ
dt dt
ΔK þ ΔU ¼ 0 ð2:7Þ
So
ðB
ΔU ¼ W ¼ F dx ð2:8Þ
A
dU
F¼ ð2:9Þ
dx
This shows forces required to change the potential energy, whereas the force
required to change the kinetic energy is determined from the definition of K.
1 1
K ¼ mv2 ¼ mx_ 2 ð2:10Þ
2 2
So
∂K
¼ mx_ ð2:11Þ
∂x_
L¼KU ð2:14Þ
This is Lagrangian formulation for a system with motion only in single axis.
Now consider a motion in three axes and every equation now to be represented in
three dimensional vector notations.
*
* d2 x
F ¼m ð2:17Þ
dt2
*
The force F has three components, Fx, Fy, and Fz in x, y and z axes respectively. The
potential energy U is a scalar field in three dimensions and the negative gradient of
this field provides a force vector in three axes. Now Eq. (2.9) represents as
T
* ∂U ∂U ∂U
F ¼ ∇ U ðx; y; zÞ ¼ ð2:18Þ
∂x ∂y ∂z
22 2 Lagrangian Modeling
This provides the three components of a force through kinetic energy Eq. (2.19) as
T !
* d ∂K ∂K ∂K
F ¼ ð2:20Þ
dt ∂x_ ∂y_ ∂z_
This gives us Lagrangian formulation of rigid body (particle) for each axis in
Cartesian coordinates as
d ∂L ∂L
¼0
dt ∂x_ ∂x
d ∂L ∂L
¼0 ð2:22Þ
dt ∂y_ ∂y
d ∂L ∂L
¼0
dt ∂z_ ∂z
x ¼ f ðr; θ; φ; tÞ
y ¼ f ðr; θ; φ; tÞ ð2:23Þ
z ¼ f ðr; θ; φ; tÞ
The functions in Eq. (2.24) follow the chain rule for respective derivatives because
x, y, z are differentiable with respect to time as well as with respect to each
generalized coordinates. Now we follow the convention that x_ ¼ dx dt ; a dot on top
of the variables represents its differentiation with respect to time for both Cartesian
and generalized coordinates. The change with respect to time in Cartesian coordi-
nates is defined as
∂p1 ∂p ∂p
p_ 1 ¼ q_ þ 1 q_ þ 1 q_
∂q1 1 ∂q2 2 ∂q3 3
∂p2 ∂p ∂p
p_ 2 ¼ q_ þ 2 q_ þ 2 q_ ð2:25Þ
∂q1 1 ∂q2 2 ∂q3 3
∂p3 ∂p ∂p
p_ 3 ¼ q_ 1 þ 3 q_ 2 þ 3 q_ 3
∂q1 ∂q2 ∂q3
We observe that
! !
d ∂p ∂p d ∂pi
p_ i i ¼ €pi i þ p_i ð2:28Þ
dt ∂qj ∂qj dt ∂qj
24 2 Lagrangian Modeling
∂p_ i ∂pi
¼ ð2:30Þ
∂q_ j ∂qj
As pi is a function of qj then the partial derivatives with respect to time are given as
! !
d ∂pi X3
∂ ∂pi
¼ q_ ð2:31Þ
dt ∂qj k¼1
∂q k ∂qj k
We know that in complex or real fields, second-order cross partial derivatives are
equal for any position variable with respect to other generalized coordinates
!
∂ ∂pi ∂ ∂pi
¼ ð2:32Þ
∂qj ∂qk ∂qk ∂qj
We know that change in work done is the dot product of force applied and
corresponding change in position given as
* *
δW ¼ F δp ð2:37Þ
*
Equivalently, force vector F has three components FP1 , FP2 ; and FP3 , such as
X
3
δW ¼ FPi δPi ð2:38Þ
i¼1
X
3 X
3
δW ¼ m€pi δPi ¼ FP i δ P i ð2:40Þ
i¼1 i¼1
Let
X
3
∂ Pi
Fq j ¼ FP i ð2:42Þ
i¼1
∂ qj
mX 3
K¼ p_ 2 ð2:45Þ
2 i¼1 i
There are three Lagrange equations which generate from Eq. (2.44) for each axis
ðj ¼ 1, 2, 3Þ as
X ! 2 !
3
d ∂ p_ i 2 ∂ p_ i
m ¼ Fq j ð2:46Þ
i¼1
dt ∂ _
q j 2 ∂q j 2
Each value of j represents a Lagrange equation in its axis. The purpose of keeping
δW in Eq. (2.44) is that we need to equate change in energy through work done by
other physical interpretations. The Lagrange equation deals with the kinetic energy
of an object in motion in Eq. (2.44). Work done is also represented by change in
potential energy as given in Eq. (2.8). There may be non-conservative forces that
are acting on the body or internal energy of the system, but due to the law of
conservation of energy these forces must balance each other. The negative gradient
of potential energy provides forces acting upon the body which may be causing it to
move or stop. So if we define potential energy U in Cartesian coordinates, forces in
Cartesian coordinates are represented as
∂U
FP i ¼ ð2:47Þ
∂pi
X 3
∂U ∂Pi ∂U
Fq j ¼ ¼ ð2:48Þ
i¼1
∂pi ∂qj ∂qj
X ! 2 !
3
d ∂ p_ i 2 ∂ p_ i ∂U
m ¼ ð2:49Þ
i¼1
dt ∂ _
q j 2 ∂q j 2 ∂qj
Again using the definition of kinetic energy from Eq. (2.45) to represent Eq. (2.49)
in energy variables
2.3 Work and Energy Formulation 27
!
d ∂K ∂K ∂U
¼ ð2:50Þ
dt ∂q_ j ∂qj ∂qj
Now using a Lagrange variable L from Eq. (2.14) in three axes system and using
Eq. (2.51) in each coordinates we get the following equation for a conservative
system in each qj coordinate
!
d ∂L ∂L
¼0 ð2:52Þ
dt ∂q_ j ∂qj
In a conservative field Eq. (2.7) holds for change in both potential and kinetic
energy, but in a non-conservative field the change is also due to non-conservative,
external, or internal forces in each coordinate system. Enc is a non-conservative
energy due to all of these forces acting upon the body, so the energy equation
becomes
ΔK þ ΔU ¼ Enc ð2:53Þ
In this case the Lagrange equation (2.52) does not satisfy the due non-conservative
force acting upon a body. Let Fncqj is the sum of all non-conservative forces acting
in qj direction. The generalized Lagrange equation in generalized coordinates is
!
d ∂L ∂L
¼ Fncqj ð2:54Þ
dt ∂q_ j ∂qj
Depending upon the space in which coordinates are defined, the value of j can be
changed. For a motion with two degrees of freedom there will be only two Lagrange
equations but in a motion with three degrees of freedom there will be three
Lagrange equations to describe the system. In certain cases, there is motion in
both Cartesian coordinates and generalized coordinates; then accordingly a direc-
tion of motion in Cartesian is also treated as a variable qj.
28 2 Lagrangian Modeling
The Lagrange equation itself is a nonlinear or linear state space representation. All
non-conservative forces are treated as input to the system whether these are
controllable or exogenous. The generalized coordinates and their first-order deriv-
atives constitute the state vector.
Example 2.1: A Simple Pendulum The simple pendulum is a body with mass
m attached at length l of massless string from the origin (or ground) to move
between two end points as shown in Fig. 2.1. If no external force is applied then
the movement will depend upon the starting point and eventually comes to end at
mass exactly below the origin. In Cartesian coordinates x and y mass move in a
semicircle, which can also be related through polar coordinates l and θ, which are
fixed length of a string and angle from the origin respectively. There is only one
degree of freedom, so there is only one generalized coordinate q1 ¼ θ. Conven-
tionally these simple examples are represented with θ instead of q1. Now we define
the position and velocity of Cartesian coordinate as a function of a generalized
coordinate as
xðθÞ ¼ l sin θ
ð2:55Þ
yðθÞ ¼ l cos θ
x_ θ; θ_ ¼ lθ_ cos θ
ð2:56Þ
y_ θ; θ_ ¼ lθ_ sin θ
1 1 1
K ¼ mv2 ¼ m x_ 2 þ y_ 2 ¼ ml2 θ_ 2 ð2:57Þ
2 2 2
−
2.4 State Space Representation 29
1
L ¼ ml2 θ_ 2 mglð1 cos θÞ ð2:59Þ
2
€θ þ g sin θ ¼ 0 ð2:62Þ
l
State space formulation of the system represent Eq. (2.62) with definition
of state
variables given as x1 ¼ θ and x2 ¼ θ. _ We have state vector x ¼ θ and the
*
θ_
nonlinear state space representation as
x_ 1 ¼ x2
g
x_ 2 ¼ sin x1 ð2:63Þ
l
T
The system can be linearized at relaxed equilibrium point
" # θe θ_ e ¼ ½0 0
0 1
with only one matrix A ¼ g 0
l
Example 2.2: Inverted Pendulum on a Moving Cart A cart of mass M is moving by
an applied force F acting along x-axis. An inverted pendulum of mass m is attached
on the center of cart with a rod of length l (Fig. 2.2).
The motion of the cart is defined from origin and the angle is measured from
vertical y-axis. The system can either be represented with horizontal and vertical
variables of cart and bob or by using horizontal displacement of cart and angular
position of bob. There are only two degrees of freedom; an angle can be measured
from the horizontal but in order to be conversant with remaining literature we chose
the angle θ measured from y-axis. So for 2-DoF system, there are two generalized
coordinates x1 ¼ x, the cart’s displacement and θ the angle of a bob from vertical.
The position of bob (x2, y2) can be given in terms of x1 and θ as
30 2 Lagrangian Modeling
( ) M
x2 ¼ x þ l sin θ
ð2:64Þ
y2 ¼ l cos θ
x_ 2 ¼ x_ þ l θ_ cos θ
ð2:65Þ
y_ 2 ¼ l θ_ sin θ
1
K 1 ¼ Mx_ 2 ð2:66Þ
2
1 2
K 2 ¼ m x_ 2 þ y_ 2 2 ð2:67aÞ
2
1 h 2 2 i
K 2 ¼ m x_ þ l θ_ cos θ þ l θ_ sin θ ð2:67bÞ
2
1
K 2 ¼ m x_ 2 þ 2lx_ θ_ cos θ þ l2 θ_ 2 ð2:67cÞ
2
So the total kinetic energy of the system from Eqs. (2.66) and (2.67c) is given as
1 1
K ¼ K 1 þ K 2 ¼ ðM þ mÞx_ 2 þ mlx_ θ_ cos θ þ ml2 θ_ 2 ð2:68Þ
2 2
1 1
L ¼ K U ¼ ðM þ mÞx_ 2 þ mlx_ θ_ cos θ þ ml2 θ_ 2 mgl cos θ ð2:70Þ
2 2
we represent two equations of motion sing the variable L and two generalized
coordinates (x, θ) as
d ∂L ∂L
¼F ð2:71Þ
dt ∂x_ ∂x
and
d ∂L ∂L
¼0 ð2:72Þ
dt ∂θ_ ∂θ
∂L
¼ ðM þ mÞx_ þ mlθ_ cos θ ð2:73aÞ
∂x_
∂L
¼ mlx_ cos θ þ ml2 θ_ ð2:73bÞ
∂θ_
∂L
¼0 ð2:73cÞ
∂x
∂L
¼ mgl sin θ ð2:74Þ
∂θ
These are nonlinear equations for the system of a moving cart with an inverted
pendulum. In order to make a state space system, first we need to solve variables
with double derivatives simultaneously. Each generalized coordinate of the
system is second-order so a total order of the system is 4. We have both €x and €θ,
so we need to solve these equations simultaneously. Eqs. (2.75) and (2.76) are
represented as
l
ml2 €θ ¼ F ðM þ mÞ€x þ mlθ_ 2 sin θ ð2:77Þ
cos θ
l€θ
€x ¼ x_ θ_ tan θ þ g tan θ ð2:78Þ
cos θ
32 2 Lagrangian Modeling
Inserting value of €x from Eq. (2.78) in Eq. (2.75) and then inserting the value of
ml2 €
θ from Eq. (2.77) in Eq. (2.76) simultaneously solve the equations for state space
T
representation of a state vector x θ x_ θ_ as
M þ m 1 _ _ F
€x ¼ m cos θ mx_ θ sin θ þ mg sin θ þ mlθ tan θ þ
2
cos θ cos θ
ð2:79Þ
€ 1
θ¼
l
1
M m
þ m cos θ ðM þ mÞ x_ θ_ tan θ þ g tan θ þ mlθ_ 2 sin θ þ F
cos θ cos θ
ð2:80Þ
Problems
x ¼ r cos θ sin φ
y ¼ r sin θ sin φ
z ¼ r cos φ
P2.5 The position of bob of the inverted pendulum on a moving cart can be
controlled by applying external torque τ at the hinge changing the
Eq. (2.78) as
Linearize the complete system by the assumption given in P2.4 and obtain
state space representation to monitor all state variables at output
independently.
References
1. Friedland, Bernard. 2005. Control System Design—An Introduction to State Space Methods.
Mineola: Dover Publications.
2. Zak, Stanislaw H. 2002. Systems and Control. Oxford: Oxford University Press.
3. Brizard, Alain J. 2008. An Introduction to Lagrangian Mechanics. Singapore: World Scientific
Publishing.
4. Wellstead, P.E. 1979. Introduction to Physical System Modelling. London: Academic Press.
5. Amirouche, Farid. 2006. Fundementals of Multibody Dynamics—Theory and Applications.
Basel: Birkhäuser.
Chapter 3
Introduction to Bond Graph Modeling
Fig. 3.1 A half arrow representing a bond between two energy ports
variables in each energy domain with different names. The product of effort and
flow is power P(t) in every energy domain.
In bond graph representation shown in Fig. 3.2, the effort variable is written above
or to the left and the flow variable is written in under or to the right side of the half
arrow.
In electrical systems, effort is voltage and flow is current and their product is
power either in watts, VA or VAR, each of which is a representation of power in
Watts or Joules per second. In mechanical energy for translational motion, effort is
force and flow is velocity and their product is power either in Watts or horsepower
(hp). Conceptually voltage is the driving force in electrical systems and force is the
excitation in mechanical systems. Current flows in the electrical circuits and a
mechanical body moves with a velocity. Similarly, we understand efforts and flow
variables for other types of energy systems. Table 3.1 shows the effort and flow
variables for different energy systems
3.3 Energy Variables 37
In some cases it is not possible to generate a set of power variables whose product is
power. In Table 3.1, one set of power variables for a thermal system shows
temperature and heat change rate as effort and flow variables respectively. The
product is not power rather power times temperature, which conveys that a power
transfer is taking place through change in temperature. This type of bond where
product of effort and flow is not purely power is called pseudo bond graph.
Different effort and flow variables can be associated with the same energy
systems depending upon their application and integration with other systems.
Each type of energy systems have already established modeling techniques through
their respective laws of physics, but a bond graph modeling scheme is vital when
more than one energy systems are modeled together.
Power is the rate of change of energy in the system, so the change must be occurring
from either effort or flow and in certain cases from both variables. There are two
variables that associate effort and flow variables with change in the energy of the
system. These variables are called energy variables. One variable is called gener-
alized momentum p(t) and change in momentum is an effort at the port. Another
variable is generalized displacement q(t) and the change in generalized displace-
ment is the flow variable. In electrical systems, change in flux linkage is voltage/
effort so flux linkage λ(t) is the generalized momentum in electrical systems. The
rate of change of electrical charge is the current flowing in the conductor, which
provides us generalized displacement as electrical charge q(t). In a mechanical
system, momentum of motion (translational or rotational) is the generalized
momentum. As we know, mechanical momentum for translational motion is
the product of mass and velocity given as pðtÞ ¼ m vðtÞ and change in momentum
p_ ðtÞ ¼ m v_ ðtÞ ¼ m aðtÞ ¼ FðtÞ ¼ eðtÞ is a respective force. Similarly, change in
angular momentum is the torque, which is an effort variable in rigid body rotational
motion. The generalized effort and generalized momentum relationship is
defined as
ðt
qð t Þ ¼ f ðtÞdt ð3:3bÞ
0
dEðtÞ
Pð t Þ ¼ ð3:4aÞ
dt
ðt
EðtÞ ¼ PðtÞdt ð3:4bÞ
0
Now we can define energy function with either of the energy variables by substitut-
ing the relationship of effort Eqs. (3.2a and 3.2b) or flow Eqs. (3.3a and 3.3b) in
Eq. (3.5). By substituting Eq. (3.2a) in Eq. (3.5), we get energy function as a flow
variable integrated with momentum.
ðt
Eð t Þ ¼ f ðtÞ dpðtÞ ð3:6aÞ
0
ðp
EðtÞ ¼ f ðpÞ dp ð3:6bÞ
0
Using the relationship of generalized flow and displacement energy from Eqs. (3.3a
and 3.3b) for Eq. (3.5) give us energy function as an effort integrated over
displacement.
ðt
EðtÞ ¼ eðtÞ dqðtÞ ð3:7aÞ
0
ðq
EðtÞ ¼ eðqÞ dq ð3:7bÞ
0
Both Eqs. (3.6a and 3.6b) and (3.7a and 3.7b) represent energy E(t) as a function of
time, whereas momentum and displacement are autonomous functions of time.
Power is a product of flow and effort variables whereas energy is an integral of the
power. The energy is obtained by integrating either effort variable with respect to
displacement or flow variable with respect to momentum. So in a bond there are
actually four variables related together; two are power variables, which are inte-
grals of their respective energy variables. In many cases both effort and flow are
often not measureable, and power is determined through integration of momentum
or displacement variables. This relationship was explained as a tetrahedron of states
3.5 Word Bond Graph 39
between power and energy variables in [1]. Energy function in Eqs. (3.6a and 3.6b)
and Eqs. (3.7a and 3.7b) is an accumulation of integration at any instantaneous
value, so both energy variables are associated with energy-storing elements of a
physical system. Energy variable is an integral of a power variable and stores
energy with other power variables for different systems shown in Table 3.2.
The input to the system is excitation or cause to excite the model which produces
effect or response as output. The half arrow representation of a bond determines a
power input, which is flowing through one end to half arrow end. The source of
power attached with one end provides power to the other end of system (half arrow
side) conventionally given in Table 3.3. The relationship from output (effect,
response) to input (cause, excitation) is called transfer function. In real-time
modeling, input is excited with impulse, step or at different frequencies to check
the dynamic response of the output.
Word bond graph is a first step to represent the block diagrams of a model towards
complete bond graph modeling. These bond graphs explain the structure of the
system, and the flow of power and interconnectivity of the components. These are
like signal flow diagrams of the electrical circuits or free body diagrams. We first
draw an electromechanical system as word bond graph and later expand it into a full
structured bond graph model. The word bond graphs of various energy ports are
given in Table 3.4
40 3 Introduction to Bond Graph Modeling
Force
Torque
It is important to understand that the power variable is the product of effort and
flow variables, which means either effort or flow is going into the bond (port) or
coming out of the port. When two components are connected through a bond, if
one component is taking effort as input then it must be giving flow as output or
vice versa. Causal stroke determines the output of bond as either effort or flow
variable. Conventionally we place causal stroke at the side of the component,
which is providing flow to the other component. So causal stroke end specifies
where flow is being determined. Assigning causality through a stroke helps to
write algebraic relationships between ports by determining which port is provid-
ing effort or flow. We know that the law of conservation of energy forms the
basic principle of bond graph modeling by equating powers at both sides.
Accordingly, if one variable is known at one end, then the other end measures
other power variables.
Two systems, A and B, interconnect with equal power at both ends as shown in
Fig. 3.3, so there are a variety of ways in determining effort, flow, and power
transfer between two systems. There are two choices for power transfer and two
choices for assigning causality through a stroke at each end. Overall there are four
possible representations with respective understanding as given in Table 3.5.
If we connect two elements and power transfer is not taking place, but either effort
or flow is providing a transfer, then this action is represented by active bond. Active
bond shown in Fig. 3.4 is denoted by a full arrow with either effort or flow marked
as a main variable and the other variable assumed as zero. In reality, there is a very
low power being transmitted but cannot be equated with the other end and preced-
ing systems. This situation arises when some other power system is driving some
3.7 Active Bond 41
Belt pulley ( ) ( )
Belt Pulley
( ) ( )
Cam rod ( ) ( )
Cam Rod
( ) ( )
Crank rod ( ) ( )
Crank Rod
( ) ( )
Electrical transformer
1( ) 2( )
Transformer
1( ) 2( )
Universal coupling
1( ) 2( )
Cam Rod
1( ) 2( )
components and thus represented by active bond. A voltage amplifier takes voltage
input and then converts it by multiplying with a fixed gain value. The power at the
input of the amplifier is not equal to the power at the output because the
42 3 Introduction to Bond Graph Modeling
( )
( )
Fig. 3.4 Active Bond Graphs, (a) Active Effort and flow is zero, (b) Active flow and effort is zero
( )
Sensor (any power or energy variable) ( )
Sensor Req.
( ) Variable
amplification is also receiving power from another source as shown in Table 3.6 for
different sources. Graphically, we represent active bond without any causality
assignment; either effort or flow is the key input through this bond and the other
variable is assumed at very low value and approximated as zero.
3.8 Connection of Bonds 43
The bonds or ports are connected with each other and thus form a complete bond
graph. A block diagram representation followed by a word bond graph is usually the
first step. At each bond, power is equal to the proceeding or succeeding bond, and
thus complete variables can be solved.
Example 3.1: Bond Graph of Multienergy System Construct a word bond graph of
the block diagram shown in the system of an overhead storage tank filled through
electric power by a driving shaft coupled with AC motor and hydraulics pump.
Calculate the following
(a) What is the driving current if the system is consuming 1 kW of power?
(b) What is the constant pressure at the pump when a 400 L tank fills in an hour?
(c) Drive shaft is reducing the speed of motor by one-tenth, and torque coupling of
the shaft with AC motor is 5 Nm/A. What is torque and the speed of drive shaft
at the pump?
(d) Drive a numeric relation between pressure of the pump and the torque of the
drive shaft (Fig. 3.5).
Solution A simply connected word bond graph of a system is given in Fig. 3.6, it
shows the coupling of the interconnected systems:
(a) In a AC 220 V system, driving current of the motor is
PðtÞ
iðtÞ ¼ ð3:8aÞ
V ðtÞ
1000
iðtÞ ¼ ¼ 4:55 A ð3:8bÞ
200
phi
T
Fig. 3.5 Block diagram of a water storage system. AC motor driving a shaft, a shaft is rotating a
pump which is storing water in over head tank
Fig. 3.6 Word bond graph of a water storage system. AC motor driving a shaft, a shaft is rotating
a pump which is storing water in over head tank
44 3 Introduction to Bond Graph Modeling
ΔVol 0:4 m3
Q_ ¼ ¼ ¼ 11:1 103 m sec ð3:9Þ
3
P 1000
Phyd ¼ ¼ ¼ 9 M Pa ð3:10Þ
Q_ 11:1 103
(c) The torque of drive shaft is linked with current of the motor as
P 1000
ω1 ð t Þ ¼ ¼ ¼ 44 rad= sec ð3:12Þ
τ1 22:73
ω1 ðtÞ
ω2 ð t Þ ¼ ¼ 4:4rad= sec ð3:13Þ
10
P 1000
τ 2 ðt Þ ¼ ¼ ¼ 227:27 Nm ð3:14Þ
ω2 4:4
(d) The pressure of the pump and the torque of the motor is related as power
variables
ω1
Phyd ¼ τ1 ð3:15Þ
Q_
References 45
Problems
Name Figure
Gear
Timing Belt
Thermistor
Tachometer/encoder
Fork (a combination of
sources)
References
1. Karnopp, Dean C., Donald L. Margolis, and Ronald C. Rosenberg. 2012. System Dynamics—
Modeling and Simulation of Mechatronics Systems, 5th ed. Hoboken: Wiley.
2. Borutzky, Wolfgang. 2010. Bond Graph Methodology—Development and Analysis of Multidis-
ciplinary Dynamic System Modeling. London: Springer.
Chapter 4
Elements of Bond Graph
Components and elements bond the system, which define the relationship between
two consecutive ports. These bonding elements follow constitutive laws for equat-
ing power between bonds. These elements are analogues in different energy
systems as we discussed effort and flow variables in Chap. 3. If a bond or port
has only one element at either end and does not require any other element to define
it, then it is a 1-port element. If an element requires two bonds for defining
relationship, then it is a 2-port element. If three or more bonds connect for a
constitutive law, then these elements are called junctions. It is important to note
that these 1-port and 2-port elements and junctions are different than vector bond
graphs, which will be discussed later.
There are five 1-port elements, two of which are sources, one of which is dissipa-
tive, and two of which are energy storing elements.
4.1.1 Sources
Source of generalized effort and source of generalized flow are both 1-port ele-
ments. Voltage sources in electrical systems and applied force or applied torque in
mechanical systems are sources of effort. Graphically, source of effort connects at
the starting end of the bond with a causal stroke on the other side (Fig. 4.1).
Generally, sources provide power to the systems, so the half arrow side is away
from the elements where source is placed. Similarly, source of flow represents the
current source in an electrical system or the source of fixed velocity in a mechanical
system. Figure 4.2 represents the source of generalized flow with causality assign-
ment near the source.
In source of effort, effort is fixed and flow is arbitrary in the bond; in source of
flow, flow is fixed and effort is arbitrary in the bond. Causality assignments for both
sources are fixed and cannot be changed as flow is uniquely determined away from
Se and at the source in Sf. In different literature, it is often confusing when sources
are mentioned particularly for energy systems, for example SF is also mentioned as
a source of force, which is an effort variable but also as a source of flow. SV is often
mentioned for the source of velocity, which is flow variable as well as for source of
voltage, which is an effort in electrical systems. In this text, Se is a generalized
source of effort that represents source of voltage, force, torque, etc in all energy
systems. Similarly, source flow Sf represents current, velocity, angular velocity or
volume flow rate, etc. A source in bond graph represents the ideal source with a
given value whether constant or alternating, sinusoidal, decaying or time-based
function. If nonlinearities, dissipation, or dependencies upon other power variables
is to be considered then these sources can be modeled as a complete system. For
example, a 12-V voltage source of a car battery with power rating of 100 Ah is
dependent upon the current supplied. In a bond graph model, a simple SE with
assigned value of 12 V will model the system but it will not take the power
consideration of 100 Ah in the simulation. In order to represent a non-ideal source,
a detail bond graph model for source is required.
Resistor is an element that defines the relationship between generalized effort and
generalized flow. Nonlinear constitutive law of a resistor is given as
or
f ðtÞ ¼ φ1
R feðtÞg ð4:2Þ
element which is defined by a direct proportion between effort and flow. Linearly
these relations are given as
or
eðtÞ
f ðt Þ ¼ ð4:4Þ
R
The remaining two 1-port elements are energy storing elements, which are impor-
tant in order to obtain differential equations for the system. These two elements
define a relationship between a power variable to an energy variable. The capacitor
is a first 1-port energy storing element which directly holds a relationship between
generalized effort and generalized displacement. As generalized displacement is
integral of generalized flow then the relationship between effort and flow of the
bond with capacitor has an integral or derivative given as
ð t
eðtÞ ¼ φ1
C f ðtÞ dt ð4:5Þ
0
d
f ðt Þ ¼ φ fe ð t Þ g ð4:6Þ
dt C
50 4 Elements of Bond Graph
f ðtÞ ¼ C dedtðtÞ
or
It is clear from expression that if effort is given then flow is arbitrarily determined
or vice versa and thus affects the causality assignment accordingly as given in
Table 4.2. Generalized capacitor is preferred with integral causality capacitor over
differential causality.
In electrical systems, the capacitor is capacitance {in Farads} defined through
the relationship of applied voltage and the charge developed by relation
qðtÞ ¼ C V ðtÞ. In mechanical translational systems, the capacitor is spring con-
stant k N m1 by applied force and displacement by which spring compresses or
stretches by defining relation FðtÞ ¼ k xðtÞ. By the concept of generalized
capacitance, spring constant k is known as a stiffness parameter, which is the
inverse relation as compared to electrical systems. However, in rotational mechan-
ics, both rotational spring constant k N m rad1 and torsional constant
c rad N1 m1 are used.
d
eðtÞ ¼ φ ff ðtÞg ð4:11Þ
dt I
or
Equivalently
In case of energy storing ports, the concept of derivative (differential) and integral
causalities prompts an important step in modeling. Algebraically, effort or flow can
be determined at either side of bond but in modeling it requires a preference to
integral causality over derivative causality. In the case of sources, causality assign-
ment is fixed as either effort or flow determines at respective source end and hence
its causality cannot change. In the case of resistor, which is an energy dissipating
element, both causality assignments make an algebraic relationship between both
power variables. Both choices are fine in modeling because it does not induce or
52 4 Elements of Bond Graph
deduce existing information of the model. In the case of capacitor and inertia,
defining relationship between effort and flow is either an integral or derivative of
one variable to determine other variable. In the case of integral causality, it means
existing/past information is being added to formulate a value. Whereas, in case of
derivative causality, existing function is differentiated to deduce the other variable;
and in this case no new information is added and it makes the system dependent
upon already available values. Algebraically, integrals add on and produce new
information whereas derivatives use the existence of the functions or variables.
When an integral causality is used, it means a new state variable can be determined
and when derivative causality is used, it means a state variable is dependent upon
some other state variable and it is a redundant system. Physical systems can be
redundant but it increases the order of the system to solve without adding more
information. When assigning causality, it is preferred to avoid derivative causality
in order to avoid having a redundant equation to solve; dependency on other
variables as well as it may produce constraints in the model.
These elements connect two bonds together, ideally by using the law of conversa-
tion of energy for power transfer. In real scenarios, systems are not efficient to
transfer complete energy into other systems without loss, but even these losses can
be modeled with bond graphs by appropriately addressing the dissipation or loss.
There are two basic 2-port elements, namely transformer and gyrator, ideal for
power transfer linearly from one system to other.
4.3.1 Transformer
Ideal power transformer TF represented with two bonds in Fig. 4.3 means that
power at one bond completely transfers to other without any loss by implying
Transformer is a device that relates effort at bond one to the effort of bond two and
flow of bond two with flow of bond one by a transformer ratio or modulus m.
e1 ðtÞ ¼ m e2 ðtÞ
ð4:16Þ
m f 1 ðt Þ ¼ f 2 ðt Þ
n e1 ðtÞ ¼ e2 ðtÞ 1
n¼ ð4:17Þ
f 1 ðt Þ ¼ n f 2 ðt Þ m
In both equations, we can see effort at one end relates with effort at the other end of
the transformer. Accordingly, flow relates with flow in both ends of the transformer,
and variables m and n are transformer ratios. In order to avoid confusion and
repetition, we will use only m as a transformer ratio in the text. In electrical systems,
electrical transformer relates voltages and current between primary and secondary
coils by turns ratio as modulus or transformer ratio. There are numerous examples
of the transformer in mechanical systems, such as a lever mechanism which relates
forces F1 & F2 at both ends and velocities V1 & V2 at both ends, and the transformer
modulus is the ratio of lengths m ¼ L1 =L2 from the center point as shown in Fig. 4.4.
The direction motion of lever will depend upon the magnitude of moment arm and
accordingly the lever will move up and down from either end.
Similarly, in gears, power transmission is equal in rotation and torques τ1 & τ2
and angular velocities ω1 & ω2 are related with ratio of teeth between first gear and
second gear (Fig. 4.5).
There are two choices of causality assignment for a transformer; simply if effort
is known at one end of the first bond then automatically effort is automatically
determined at the same end of second bond. Similarly, flows determine accordingly,
one at transformer end and one away from transformer end. Causality assignments
by determining flow has two choices in this case as given in Table 4.4
4.3.2 Gyrator
The second 2-port element is the gyrator represented in Fig. 4.6, which evolved
from the concept of the mechanical gyroscope in which applied force affects the
velocity of precession, and force in the rotational axes provides a translational
motion towards the direction of force. By the law of conversation of momentum and
energy, powers are equal at both sides but relating forces with velocities.
It is difficult to quote a physical device which directly acts like a gyrator in the
same energy system. But there are numerous devices, sensors, and actuators which
behave like gyrators, such as the dc motor which performs a gyrator action between
torque at the shaft and the current of the armature and angular rate with the voltage.
Another example is a speaker with a magnetic coil that has relation of sound
pressure to the current in the coil and velocity with the voltage. These relationships
demonstrate the use of a gyrator in the physical sense. In most of the cases, gyrator
action performs the power transformation between two forms of energy. A gyrator
relates the effort at first bond with the flow of second bond through gyrator ratio or
modulus of gyrator r. So, it also inversely relates the flow of the first bond with the
effort of the second bond.
e1 ðtÞ ¼ r f 2 ðtÞ
ð4:18Þ
r f 1 ðtÞ ¼ e2 ðtÞ
The causality assignments of the gyrator also have two choices, i.e., either effort or
flow determines at the gyrator end as given in Table 4.5.
4.4 Junctions
These are multiport elements that require connections of three or more bonds or
subsystems governing by total power conservation in all bonds represented in
Fig. 4.7. The net power at a junction is conserved, which means that total power
in is equal to total power out at a junction. There are two types of junctions known
as 0-junction and 1-junction. These junctions represent the series or parallel con-
nections between the subsystems. The idea of series or parallel connections is not
universal when it comes to electrical or mechanical systems and also varies with
generalized efforts as voltage or force, generalized flow as current or velocity. Even
adding basic 1-port elements (resistor, capacitor, and inductor) in series and parallel
yields different configurations; because resistor in series is the sum of all resistors,
but capacitors in series are the sum of their reciprocals just like resistors in parallel.
4.4 Junctions 55
5 4
Same is the case for series connection in electrical systems equivalent to the parallel
connection of mechanical systems analytically. Though 0-junction and 1-junction
connect subsystems in either series or parallel configurations, we cannot classify
either junction as a series or parallel junction. Series and parallel connections
conserve power by either fixing the effort or flow at the junction.
The power at this junction J, which is either 0-junction or 1-junction, can be
written as
It is important to understand that this equation from bond graph depends upon
direction of bond as well as causality assignment. In this bond all arrows are coming
towards junction so the sum of all incoming power is equal to zero. In most cases, if
some bonds are supplying power, there are other bonds that are consuming power as
well so that the directions of arrows are reverse accordingly for these bonds.
Figure 4.8 shows a junction with three ports inward to the junction and two ports
outward with power relation given in Eq. (4.20).
4.4.1 0-Junction
0-Junction is a flow junction in which flow adds at the junction while effort remains
constant (Fig. 4.9).
4.4.2 1-Junction
1-Junction is an effort junction or common flow junction, where effort adds up and
flow remains constant at the junction shown in Fig. 4.10. This bond represents
a series combination in electrical systems using Kirchhoff’s voltage law.
4.4 Junctions 57
It is the common flow junction, so flow remains constant in all three bonds and the
sum of efforts are equal to zero. This is same as sum of voltages equal to zero in a
loop, or sum of forces at a point is equal to zero in translational mechanics and sum
of torques equal to zero at a point in rotational mechanics. There is a single choice
of causality assignment because one flow is determined at the bond and rests are
equal to this flow. Effort at this bond is equal to the sum of the remaining efforts
algebraically as shown in Table 4.7. In this junction, we can also observe that f1(t) is
an input to the 1-junction and flows at bonds 2, 3, and 4 are equal to the f1(t).
The effort of bond 1, e1(t) is determined by equating to the sum of remaining three
efforts. This defining bond like bond 1 in this case is often referred as a strong bond
of 1-junction and others as weak bond of this particular 1-junction. This is also clear
with the half arrow direction in this specific case where power is only input from
bond 1 and the distributed power is output from rest of the bonds.
Example 4.1: Equations of Junctions Write the equations of the following junc-
tions in Fig. 4.11
58 4 Elements of Bond Graph
Solution
(a) Bond 1 and Bond 3 are inputs to the junction and bond 2 and 4 are outputs of
the junction so the equation will be
This is a common flow junction and by causality assignment we know that f1(t)
is determined before the bond so remaining flows are
f 2 ðtÞ ¼ f 1 ðtÞ
f 3 ðtÞ ¼ f 1 ðtÞ ð4:24Þ
f 4 ðtÞ ¼ f 1 ðtÞ
Equations (4.24) and (4.26) are solution equations for this 1-junction.
(b) The 0-junction has two inputs and two outputs and this is a common effort
junction. Bond 2 is a strong bond of this 0-junction and so e2(t) is determined
at the junction and remaining efforts are equal to it as given in Eq. (4.27)
e1 ðtÞ ¼ e2 ðtÞ
e3 ðtÞ ¼ e2 ðtÞ ð4:27Þ
e4 ðtÞ ¼ e2 ðtÞ
The modulus of a transformer and gyrator may vary depending upon other variables
in the system or due to a combination of few variables. This modulus now connects
through an active bond, which means it does not represent power transfer, rather it
only supplies a variable criterion for transformer and gyrator ratio. The modulus is
usually represented as a function of state variables in the system or may be
dependent upon input or output. In electrical transformers, output voltage is usually
4.6 Modulated Sources 59
Fig. 4.12 Bond graph of modulated 2-port elements (a) A modulated transformer, (b) a modu-
lated gyrator
fixed at a desired level and the number of turn ratio is changed depending upon load
requirement and input current. A door with a hinge joint is a modulated transformer
because it depends upon the angle at which force applied, and this angle and the
moment arm make a torque which revolves the door. Figure 4.12 shows the
graphical representation of the modulated transformer and gyrator.
Remember that there is no causality assignment for active bond graph and the
modulus of transformer or gyrator represents without any causality assignment. In
these cases, a modulus is a function of variable which is dependent upon for
defining relationship. For example in order to relate the torque/angular rate of
revolving door with applied force and translational velocity, we relate the modulus
as m(θ) in a door with hinge joint. In this case θ is generalized displacement on the
side of the transformer. Problem 4.5 elaborates more on this concept.
If a source represents a time function or is dependent upon some other variable, then
these sources are represented as modulated sources. These sources take an input
through an active bond and then supply effort or flow by a normal bond to the
system. Dependent sources, like voltage-dependent voltage source, current-
dependent voltage source, etc take an input from their dependencies through active
bond and supplies output in the form of effort or flow to the system.
Problems
P4.1 Find the second gyrator equation and modulus r as per the definition given in
Eq. (4.18) if
f 1 ðtÞ ¼ s e2 ðtÞ
P4.2 Find the constitutive equations at all bonds for the following two systems and
explain why there is a differential causality with transformer but not with
gyrator?
60 4 Elements of Bond Graph
P4.3 A piston and cylinder component is shown in the figure below, which is
converting mechanical translation into hydraulics system. Define a 2-port
element for the system and explain its modulus with constitutive law between
both systems
P4.4 Define constitutive laws for each element given in Tables 3.3, 3.6 and
Problem 3.1 and specify the type and causality assignments for each element.
P4.5 A hinge movement of a revolute joint is shown in figure below which trans-
forms translational energy into rotational energy. The force applied resolves
into two components, and torque at hinge makes a moment arm with a
component. Modulated transformer represents this scheme in bond graph
methodology, and modulus is dependent upon the length of the moment arm
and angle at which force is being applied. Draw a word bond graph of this
scheme and formulate a constitutive with definition of modulus as the func-
tion of variables associated with the model.
.
4.6 Modulated Sources 61
P4.6 There are modulated sources of efforts and flow which relate time as a
mathematical function to effort or flow. An example is AC electrical power
source with a fixed frequency component. These sources are represented as
active bond graph at input with either effort or flow as a main output power
source. Write the constitutive equations of following models in each
subsystem with modulated sources and 2-port elements.
(a)
(b)
(c) Note that bond 7 is the effort of bond 6 only through active bond,
represents the effort sensor in the respective bond. The bond 6 remains
same before and after sensor in all respects of causality and power
transfer.
P4.7 Find a constitutive law between bond 1 and bond 3 with all possible choices of
causality and simplify the each model
References
1. Karnopp, Dean C., Donald L. Margolis, and Ronald C. Rosenberg. 2012. System Dynamics—
Modeling and Simulation of Mechatronics Systems, 5th ed. Hoboken: Wiley.
2. 20-Sim Software. www.20sim.com
3. François E. Cellier, and N. Ángela. 2005. The Modelica Bond Graph Library. In Proceedings of
the 4th International Modelica Conference, 57–65. Hamburg, Germany, 7–8 Mar 2005.
Chapter 5
Analytical Formulation by Bond Graph
Modeling
Bond graph modeling possesses more applicability towards systems with different
energy domains. It is easier to model electrical systems with electrical theories
using KCL or KVL. It is easier to model mechanical systems using Newtonian
mechanics, Lagrangian mechanics, etc. Similarly, for electromagnetic systems
there are Maxwell laws to deal with both electrical and magnetic forces; in
hydraulics systems, there are Bernoulli’s /Stroke’s law to deal with pressure and
forces. But now, with the advent of modern technologies, there are systems which
are now amalgams of electrical, mechanical, hydraulics, electromagnetic systems,
etc. Bond graph modeling provides help to these modeling schemes, and dealing
with systems in their energy domains provides an easier and technically mature
approach. Bond graph provides a unified approach to model a system comprised of
two or more energy systems. On the other hand it is difficult to provide a universal
methodology that is applicable to all systems. In electrical and hydraulic systems,
nodes appear where generalized flow distributes between them whereas, in mechan-
ical systems, nodes appear where forces act together both in series or parallel
configurations. As we discussed in previous chapters, that series and parallel
combination is not a universal phenomenon for modeling. So there arises a problem
in formulating a unified procedure for modeling at system level. A generalized set
of procedures are written for modeling the systems that is applicable for modeling
systems within same energy domain as well as modeling in multi-energy domain. It
needs careful attention when a series or parallel combinations of different energy
systems are modeled as a single system.
1 0 1
1 R
Se 1 0 1 I
Step 3: There are nodes where flow is different but a zero is already placed, so skip
to next step.
Step 4: Inserting 0 and 1 junction between 0 and 1 such that 0 and 1 alternate
(Fig. 5.3).
Step 5: Insert elements (Fig. 5.4).
Step 6: Insert bonds (Fig. 5.5).
Note: Elements are joined with inserted 0/1 junctions between the bonds.
Step 7: Delete reference node (Fig. 5.6).
Step 8: Deleting extra bonds and simplifying: There are two 0-junctions with less
than three ports and a 1-junction with less than three ports. 1-0-1 can be
simplified as 1 because both Se and C have the same flow, similarly R2 and
66 5 Analytical Formulation by Bond Graph Modeling
0 1 0 1 0
1 R
Se 1 0 1 I
0 1 0 1 0
1 R
Se 1 1 I
Se 1 0 1 I
I has the same flow. R1 is in parallel with two 1-junctions, so this 1-junction can
be removed (Fig. 5.7).
Step 9: We know that Se ¼ V DC and C have the same flow but different effort
(current in series components) and R2 and I have the same flow and different
effort. Flow is distributing between three links in series represented by a
0-junction, where effort is the same in these links. This combination logically
conforms to the series and parallel connections in the circuit and verifies the
bond graph representation.
Step 10: Causality assignment is easy in this case, as we assign causality to source
of effort Se and integral causality to C and I, the rest will follow in line.
Step 11: Finally numbers are assigned to each bond in the model (Fig. 5.8).
Se 1 0 1 I
Step 10/11: There are two choices of causality assignments depending upon which
capacitor is chosen for differential causality (Figs. 5.11 and 5.12).
Example 5.3 Dependent Sources
Draw a bond graph of the following system with a voltage-controlled current
source for measuring current passing through resistor R6 and voltage drop across R2
(Fig. 5.13).
68 5 Analytical Formulation by Bond Graph Modeling
Solution This model seems complicated but in reality it is a very simple bond
graph model to construct. The complete system is an electrical network with an AC
voltage source and a dependent voltage control current source. We leave measure-
ments for a while for constructing a bond graph representation. We can directly
reach till Step 8 as we know the points where effort is distributing and where flow is
distributing.
1. An AC voltage source is in series with resistor R1 and in series with parallel
combination of capacitor C1 and inductor L1. A 1-junction connects a voltage
source MSe, R1 and two ports to 0-junctions in this part. A 0-junction connecting
this 1-junction, a capacitor C1 and Inductor L1 represent the parallel combination
in series.
2. The 1-junction in step 1 is in parallel with resistor R2 and a series combination of
R3 and TF. A 0-junction has two 1-junctions on both sides; one connecting with
1-junction of Step 1 and a second connecting with a series combination of R3
and TF.
3. The other side of the transformer is also connected with a 1-junction, which
connects R4 and a 0-junction.
4. The 1-junction of Step 3 is followed by a 0-junction, which connects a capacitor
C2 and a voltage-controlled current source MSf.
5. The output of modulated source of flow MSf is connected with a series combi-
nation of R5, R7, and a parallel combination of L2 and R6 represented by
1-junction connecting a 0-junction like Step 1.
5.1 Modeling Procedure 69
R1 R2 R3 R4 C2 R5
R R R R C R
MSe 1 0 1 TF 1 0 MSf 1 0 R
R6
Sine
0 C
R I
C1
R7 L2
I
L1
R1 R2 R3 R4 C2 R5
R R R R C R
MSe 1 0 1 TF 1 0 MSf 1 0 R
R6
Sine
0 C
R I
C1
R7 L2
I
L1
R1 R2 R3 R4 C2 R5
R R R R C R
m
2 7 9 12 14 17
1 6 8 10 11 13 15 16 19 20
MSe 1 0 1 TF 1 0 MSf 1 0 R
R6
Sine 3 18 21
4
0 C
R I
5 C1
R7 L2
I
L1
Fig. 5.16 A complete bond graph of an electrical network given in Fig. 5.13
Fig. 5.17 A complete bond graph of an electrical network given in Fig. 5.13 with two sensors
5.2 Sensors
The current passing through R6 and voltage drop in R2 can be measured by placing
sensors in 20-Sim software and these constitute the required measured output
vector. A sensor for each requirement can be placed using 20-Sim software, and
there is no need to assign a new bond number because it does not change causality
or mathematical formulation. These sensors are used for simulation purposes for
monitoring desired response or feedback to the control mechanism (Fig. 5.17).
Solution The difference between mechanical and electrical system is the logic of
series and parallel combination. In mechanical systems, velocities are measured
with respect to ground or reference. 1-junctions are placed for different velocity
nodes in the system.
Step 1: There are three nodes in this system with different velocities, zero velocity
at ground, velocity between spring, mass, & damper and at the other end of mass
where force is pulling the mass.
Step 2: In this system there is one uniquely determined effort so we skip to Step
3 and flow distribution varies with different nodes.
Step 3: There are two places where effort is distributing and the flow is same, so we
place 1-junctions at two points.
Step 4: Insert zeros between 1-junctions
Step 5: Add 1-port elements to the zeros (Fig. 5.19)
Step 6: Add half-arrow bonds between the junctions (Fig. 5.20)
72 5 Analytical Formulation by Bond Graph Modeling
Step 7: We now delete the ground represented as the source of flow which is zero or
reference velocity in the system (Fig. 5.21).
Step 8: Now there are two zero junctions and a 1-junction which have only two
ports, we squeeze these together.
Step 9: Power flow can be verified that a single source of effort is pulling with three
components with a same velocity and different distribution of applied force.
Step 10: A source of effort, integral casualties for capacitor (spring), and inertia
(mass) leave a single choice for resistor (damper).
Step 11: Numbers are assigned for four bonds in the system (Fig. 5.22).
The other end is also a 1-junction, which is connected with gyrator, inertia, and a
transformer.
3. The crank rod mechanism is represented as a transformer as it relates the torque
to translational force as well as angular velocity of wheel to translational
velocity of rod.
4. The other end of this transformer is connected with a spring, this is represented
by a 0-junction where force remains the same for a rod pulling the spring
followed by a spring pulling a belt.
5. Belt pulley is directly connected to the mass so it is not transforming variable so
there is no need to connect any element to a fixed belt pulley.
6. The transformer of Rack-Pinion system is now connected to a 1-junction where
speed or mass remains the same under two different forces. One force is applied
through DC motor via mechanical assembly and the other is gravitational force.
7. A position sensor is a generalized displacement sensor connecting to PID control
system which provides feedback to modulated source.
Step 4: There is no need to place 0 or 1 junction in between existing junctions.
Step 5: There are five elements that are directly connected to the respective 0 and
1-junctions. One modulated source of effort and force under gravity connects
with the first and last 1-junctions respectively.
Step 6: Half-arrow direction is from modulated source of effort to the rest of
components except only for a second source F ¼ mg which also directs power
from source to mass m.
Step 7: There is no ground or reference point in electromechanical loop.
Step 8: Already simplified
Step 9: Flow of power is from modulated source of effort to rest of the system.
Step 10: First we assign causality to both source of efforts, integral causality to
three generalized inertia element and a generalized capacitive element. Starting
for 1-junction from either end, the causality of a gyrator and both transformers
74 5 Analytical Formulation by Bond Graph Modeling
State space equations from bond graph models can be directly obtained in matrix
form by writing algebraic equations at each 0-junction and 1-junction. In order to
start writing state space equation, we must define state and input vectors. All
sources constitute the input vector in any order of choice. State vector is defined
by energy storing elements i.e., generalized capacitance and generalized inertia.
Generalized displacement of each capacitor and generalized momentum of each
inertia element define the state variables. The differential of generalized displace-
ment is flow in the capacitor and differential of generalized momentum across the
inertia is effort. Algebraic solution of efforts at 1-junction and flow at 0-junctions
leads to state space formulation of the model.
Example 5.6 Equation Derivations-1
Write state space equations of the following bond graph (Fig. 5.25)
Solution In this bond graph we know that there is only one source and two energy
storing elements. So we define the state vector and input vector as
* q2 ð t Þ *
x ðt Þ ¼ , u ðtÞ ¼ ½Se ð5:1Þ
p7 ð t Þ
associates with bond number 2. Similarly p7(t) represents the generalized momen-
tum related with bond 7 of 1-port inertia. Source of effort Se constitutes an input
vector of the system and in this case there is no need to represent numbers with
it. Now we can write algebraic equations for each 1-junction and 0-junction of the
model. 1-junction is a common flow and 0-junction is a common effort, so the three
set of equations are
1-junction
0-junction
1-junction
e 1 ð t Þ ¼ Se ð5:5Þ
q2 ðtÞ
e2 ðtÞ ¼ ð5:6Þ
C
f 2 ðtÞ ¼ q_ 2 ðtÞ ð5:7Þ
e7 ðtÞ ¼ p_ 7 ðtÞ ð5:8Þ
p7 ðtÞ
f 7 ðtÞ ¼ ð5:9Þ
I
76 5 Analytical Formulation by Bond Graph Modeling
There are two state variables so we should get two equations involving differential
of state variables as a function of state variables and input. Equation (5.3b) now
forms as
Whereas
e4 ðtÞ
f 4 ðt Þ ¼ ð5:11Þ
R
q2 ð t Þ
e4 ðtÞ ¼ e3 ðtÞ ¼ Se ð5:12Þ
C
So
1 q ðtÞ p ðtÞ
q_ 2 ðtÞ ¼ Se 2 þ 7 ð5:13Þ
R C I
q2 ðtÞ p ðt Þ
p_ 7 ðtÞ ¼ Se R 7 ð5:16Þ
C I
Where
p2 ðtÞ
e3 ðtÞ ¼ f 3 ðtÞ B ¼ f 2 ðtÞ B ¼ B ð5:22Þ
m
p2 ð t Þ
p_ 2 ðtÞ ¼ FðtÞ k q4 ðtÞ B ð5:23Þ
m
p2 ðtÞ
q_ 4 ðtÞ ¼ ð5:24Þ
m
Fig. 5.27 Bond graph of an electromechanical network of Example 5.5 without controller
p12 ðtÞ
There are two transformers and a gyrator with the following governing equation as
e4 ðtÞ ¼ μ f 5 ðtÞ
ð5:27Þ
e5 ðtÞ ¼ μ f 4 ðtÞ
f 8 ðtÞ ¼ α f 7 ðtÞ
ð5:28Þ
e7 ðtÞ ¼ α e8 ðtÞ
e10 ðtÞ
e11 ðtÞ ¼
r
ð5:29Þ
f 11 ðtÞ
f 10 ðtÞ ¼
r
Fig. 5.30 A hybrid electromechanical network with control system in the feedback loop
Now replacing with state and input variables yields these equations as
Where
p3 ðtÞ
f 2 ðtÞ ¼ f 3 ðtÞ ¼ ð5:32Þ
L
and
p6 ðtÞ
f 5 ðtÞ ¼ f 6 ðtÞ ¼ ð5:33Þ
J
80 5 Analytical Formulation by Bond Graph Modeling
p3 ðtÞ p ðtÞ
p_ 3 ðtÞ ¼ FðtÞ R μ 6 ð5:34Þ
L J
This implies
Where
p3 ðtÞ
f 4 ðtÞ ¼ f 3 ðtÞ ¼ ð5:38Þ
L
So
p3 ðtÞ
p_ 6 ðtÞ ¼ μ α k q9 ðtÞ ð5:39Þ
L
f 11 ðtÞ
q_ 9 ðtÞ ¼ α f 7 ðtÞ ð5:41Þ
r
Where
p6 ðtÞ
f 7 ðtÞ ¼ f 6 ðtÞ ¼ ð5:42Þ
J
p12 ðtÞ
f 11 ðtÞ ¼ f 12 ðtÞ ¼ ð5:43Þ
m
So
With state variables and a constant acting force the Eq. (5.45) becomes
e10 ðtÞ
p_ 12 ðtÞ ¼ þW ð5:46Þ
r
k
p_ 12 ðtÞ ¼ q9 ðtÞ þ W ð5:47Þ
r
The output at bond 12 can be measured easily if it is flow in this bond or effort at
this junction. The flow in bond 12 relates to velocity of a mass. But in order to
measure position, further integration of flow is required. In this bond p12 is
generalized momentum and dividing by mass m will yield a velocity and integrating
this velocity will result in displacement of mass. This integration increases the order
of the system as one more equations need to be solved. Now let us define a new state
variable q12 i.e., generalized displacement of bond as
p12 ðtÞ
q_ 12 ðtÞ ¼ ð5:49Þ
m
So state space formulation given in Eq. (5.48) is now represented with new state
variable as
2 3 2 3 2 3 2 3
p_ 3 ðtÞ R=L μ=J 0 0 0 p3 ðtÞ 1 0
6 p_ ðtÞ 7 6 μ=J 07 6 7 6 7
6 6 7 6 0 α k 0 7 6 p6 ðtÞ 7 6 0 0 7 " #
6 7 6 7 6 7 6 7 F ðtÞ
6 q_ 9 ðtÞ 7 ¼ 6 0 0 α=J 1=r m 7 6 7 6 7
0 7 6 q9 ðtÞ 7 þ 6 0 0 7
6 7 6
6 7 6 7 6 7 6 7 W
4 p_ 12 ðtÞ 5 4 0 0 k=r 0 0 5 4 p12 ðtÞ 5 4 0 1 5
q_ 12 ðtÞ 0 0 0 1=m 0 q12 ðtÞ 0 0
ð5:50Þ
The output matrix for measuring both position and velocity is given as
2 3
p3 ðtÞ
" # 6
6 p6 ðtÞ 7
7
* 0 0 0 1 0 6 7
y ðtÞ ¼ 6
6 q9 ðtÞ 77 ð5:51Þ
0 0 0 0 1 6 7
4 p12 ðtÞ 5
q12 ðtÞ
82 5 Analytical Formulation by Bond Graph Modeling
If more variables need integration then order may be increase relatively. It depends
upon sensors or required output to obtain a state space formulation. The state space
in Eq. (5.48) fulfills the modeling requirement but for a specific output, the system
needs more integration and hence the increases in system order. This can be
achieved separately in simulations by obtaining velocity from this state space and
further integrating it for position. A zero column in Eq. (5.50) also indicates the
redundancy in the system caused by the sensor, which in this case is not adding new
dynamics to the system but causing increase of order for specified variable output.
The major step in analytical formulation is to draw a bond graph representation
appropriately with logical flow of power and causality assignments. When the bond
graph has been drawn through steps given in the procedure then it will be easy to
write equations. Equations at junctions are algebraic manipulations which at the
end yield a state space representation of a model. This state space representation
can be subjected to further analysis and design. In Example 5.7, we obtain two state
space representations with 4th- and 5th-order model. The advantage of the
5th-order model is that we can obtain required output directly from model equa-
tions, but when linear analysis or controller design are performed, 4th-order model
will prove itself a better and simple representation of the system. It is also an
advantage that higher order differential equations are automatically represented in
the system of equations through efforts and flows at respective junctions. However,
these equations can be combined to obtain a single higher order system if needed.
20-Sim software can directly generate system equations by drawing bond graph
properly. The software then checks the necessary conditions of the model. If there
are conflicts of causality, loops, or improper combination of elements, the software
generates errors accordingly. The highlighted errors must be rectified and when
there are no errors left, the modeling equations can be generated. The 20-Sim
software can generate modeling equations in the presence of some warnings, but
it is suggested that warnings must be taken care of for equations and properly
understood in order to obtain a proper model of the system. By utilizing the
frequency domain analysis tools, 20-Sim can also generate symbolic representation
in state space formulation as well. Linearization in 20-Sim is only possible for
single input and single output system.
Example 5.8 Equations through 20-Sim
Verify equations for a system given in Example 5.6 using 20-Sim Software
Solutions The state space formulation of the system is given in Eq. (5.25) as
5.5 20-Sim Software Tips 83
2 3
" # B " # " #
p_ 2 ðtÞ 6 m k 7 p ðtÞ 1
¼6 7 2
4 1 5 þ ½FðtÞ
q_ 4 ðtÞ q4 ð t Þ 0
0 ð5:52Þ
m
" #
* p2 ðtÞ
y ðtÞ ¼ e4 ðtÞ ¼ ½ 0 k
q4 ðtÞ
F\p.e ¼ F\effort;
Dynamic equations:
C\p.e ¼ C\state/C\c;
C\p.f ¼ m\state/m\i;
B\p.e ¼ B\r * C\p.f;
m\p.e ¼ F\p.e(C\p.e + B\p.e);
System equations:
Static and state equations represent input and state vectors respectively.
Dynamic equations represent the state equations but their number may be equal
or more than the order. By representing the system equations in terms of derivatives
on the left hand side rather than integrals on the right hand side, we get the state
space representation. Backslash with variables show parametric assignment and
forward slash is used for division. C\state and m\state represent the q4(t) and p2(t)
respectively. Looking at the first system equation and substituting C\p.f accord-
ingly, we obtain the 2nd equation in state space formulation as in Eq. (5.24). Now
substituting m\p.e in the 2nd system equation we get the 1st equation in state space
representation i.e., Eq. (5.23). These equations are the same as the derived equation
in Example 5.6. Using “Model Linearization” from frequency domain tools, we can
get state space representation directly in 20-Sim environment as well as transfer
function of the system (with default numerical values of parameters). In addition
there are 14 removed equations after algebraic manipulation to obtain dynamic
equations, which are not needed in the process.
84 5 Analytical Formulation by Bond Graph Modeling
Problems
P5.1 Draw a bond graph of each system and represent its equation in state space
formulation and verify it using 20-Sim software.
P5.2 Use state space formulation separately for both cases of Example 5.2 using
20-Sim and verify equations of 20-Sim with equations manually calculated of
both representations.
P5.3 A good practice for drawing bond graph and generating an equation is to study
the related examples and material in the literature. Ref. [1], “System Dynam-
ics—Modeling, Simulation, and Control of Mechatronic Systems” by D C
Karnopp, D L Margolis & R C Rosenberg is a wonderful resource for bond
graph methodology in detail. Draw bond graph for each example and exercise
of Chap. 4, of this text and obtain equations manually and verify it using
20-Sim software.
P5.4 All examples in this chapter and P5.1 have been implemented in 20-Sim
software with both bond graph methodology and block diagram. It is possible
for many examples to directly obtain equations from block diagrams rather
than converting these into bond graphs. Implement these examples/P5.1 in
20-Sim using block diagrams and obtain equations and compare with state
space representation in this form.
P5.5 Obtain transfer functions of the systems in the examples in Chap. 5 and P5.1
and compare with the transfer function obtained from 20-Sim linearization
tool. Note that the numerator of the transfer function depends upon the output
chosen for the purpose whereas the denominator remains the same regardless
of the output variable.
References
1. Karnopp, Dean C., Donald L. Margolis, and Ronald C. Rosenberg. 2012. System Dynamics—
Modeling and Simulation of Mechatronics Systems, 5th ed. New York: Wiley.
2. 20-Sim Software, www.20sim.com
3. Asif M. Mughal, and Kamran Iqbal. 2013. Bond Graph Modeling and Optimal Control
Design for Physiological Motor Control System. International Journal of Modeling and Sim-
ulation, 33(2):93–101.
Chapter 6
Advance Bond Graph Modeling
q_ 3 ¼ f 4 ¼ f 5 þ f 7 ð6:1Þ
The flow (current) in the first 1-junction divides into limiting resistor and load
resistor where more simplification leads to
e5 1
q_ 3 ¼ þ p ð6:2Þ
R2 L 7
e5 ¼ e4 ¼ e1 e2 e3 ð6:3Þ
1 1 1
_q 3 ¼ V R f 2 q3 þ p7 ð6:4Þ
R2 C L
Now flow in bond 2 through resistor R is the same as the flow in bond 3, which is
being solved. This creates an algebraic loop.
1 1 1
q_ 3 ¼ V R q_ 3 q3 þ p7 ð6:5Þ
R2 C L
Now moving the term RR2 q3 on the left hand side solves this loop as
R 1 1 1
q_ 3 1 þ ¼ V q3 þ p7 ð6:6Þ
R2 R2 C L
p_ 7 ¼ e6 e8 ¼ e4 f 8 RL ð6:7Þ
6.1 Algebraic Loops 87
We get
1 RL
p_ 7 ¼ V R q_ 3 q3 p ð6:8Þ
C L 7
The cumbersome state equation now requires careful algebra in order to get the
state space representation as
1 R2 1
q_ 3 ¼ q þ p þ V ð6:10Þ
ðR þ R2 ÞC 3 ðR þ R2 ÞL 7 ðR þ R2 Þ
1 R þ 1 þ RR2 RL 1
p_ 7 ¼ q3 p7 þ V ð6:11Þ
1 þ RR2 C 1 þ RR2 L 1 þ RR2
More than one algebraic loop may result in more complex algebraic operations in
order to obtain state space equations. In this model if current limiting from the
source is taken care at just after source and before capacitor, the model will be
simplified with a zero junction and a one junction without algebraic loop. This can
further be simplified as a source of flow with a one junction as shown in Fig. 6.3 but
it will introduce a differential causality in the system. However, it must be under-
stood that this arrangement may not be possible with all the modeling assumptions
e.g., if there is a transformer in the system or the voltage drop across capacitors
needs attention, etc.
Fig. 6.3 Bond graph simplification of system in Fig. 6.1 but causing derivative causality
88 6 Advance Bond Graph Modeling
1
Se ¼ f 2 R þ q ð6:13Þ
C1 4
The effort at 0-junction can be expressed as either q4 not at q5 due to differential
causality at port 5. We now express the state equations by equating flow in
1-junction i.e., f 2 ¼ f 3 and
f 3 ¼ f 4 þ f 5 ¼ q_ 4 þ q_ 5 ð6:14Þ
1 1
q_ 4 ¼ Se q q_ 5 ð6:15Þ
R R C1 4
q4 ðtÞ q5 ðtÞ
¼ ð6:17Þ
C1 C2
Now differentiating Eq. (6.17) with time implies to
C2
q_ 5 ¼ q_ ð6:18Þ
C1 4
The second state equation is linearly dependent upon first state equation as given in
Eq. (6.18)
=C 1
C2
1 1
q_ 5 ¼ Se q ð6:20Þ
1 þ C2 =C1 R R C1 4
The A matrix has one column zero and that is due to the linear dependence, which is
the result of derivative causality. This simply means that one state equation is
enough to solve the system and the other equation is simply an algebraic multiple of
another equation. The derivative causality introduces an algebraic loop in the
systems which are required to be solved by manipulating variable in the formula-
tion. Let us consider a circuit in Fig. 6.6, where current flow in C2 is dependent upon
C1 and both are related to different state variables.
The bond graph of the circuit given in Fig. 6.6 is given in Fig. 6.7 with causality
assignment at C1, however, causality assignment in this case may be changed either
to C1 or C2. In this case the independent state variable is q2 and the dependent state
variable is q4.
In this algebraic loop, we can only get one state equation for q2 and the other
state q4 can be obtained at output. Problem P6.5 elaborates more for the complete
state space representation which requires differentiation of source.
90 6 Advance Bond Graph Modeling
1 q 1 C2=C1
q_ 2 ¼ V 2 ð6:22Þ
R C1 1 þ C2 =C1
C2
q4 ¼ C2 V þ q2 ð6:23Þ
C1
6.3 Fields
Basic 1-port elements may need higher representation in modeling, e.g., a spring is
attached with both ends exerting force or capacitors banks in electrical circuits,
different configurations of inductor as well as modulated transformers and gyrators
with own feedbacks. These basic elements can be represented as Fields, in which a
1-port element is connected with two or more ports with the same physical law. The
constitutive algebraic law will be represented in matrix form for generalized vari-
ables. Let us consider two similar examples in mechanical and electrical systems
with two generalized capacitances with two sources of efforts as shown in Fig. 6.8
and Fig. 6.9:
The respective bond graph for both systems is same and shown in Fig. 6.10
Both generalized capacitors have differential causality and it makes it more
difficult to obtain the equations. The effort, flow, or displacement functions in C1
are dependent upon C2 and vice versa. This situation can easily be represented with
Field capacitors with two ports as shown in Fig. 6.11.
6.3 Fields 91
2-I element, etc. In 20-Sim, it is necessary to change in code to represent the II-field
or resistive fields given as
The main advantage of this port-scheme is significant when there are different
energy systems with transformer or gyrator actions through fields, e.g., an electro-
magnetic relay which consists of three primary systems: electrical, magnetic, and
mechanical (Fig. 6.13).
A normally open electromagnetic relay converts electrical energy into magnetic
energy, which closes the circuits till the voltage is being supplied. As soon as
voltage across the windings is released, a mechanical spring pushes the mass and
thus opens the relay. The magnetic force equates the rotation, damping action, and
exerting force in the mechanical spring. A gyrator is needed to represent conversion
of electrical energy into magnetic energy, and a gyrator modulus is the number of
turns in the coil. In order to avoid detailed magnetic circuits, we can represent the
magnetic resistance on a zero junction followed by a C-field capacitor connecting
the spring and mechanical assembly together. A simplified bond graph of the
system is represented as (Fig. 6.14):
Remember this is a simplified version of bond graph for an electromagnetic
relay. A detailed discussion and better model for this electromagnetic relay is given
in ref [1]. This system has three state variables, two displacement variables q6, and
q7 and a momentum p8.
6.4 Series Motor 93
2 3
1 R 1 1 R 1
2
3 6 C6 N 2 Rmag þ þ 0 7
q_ 6 ðtÞ 6 C67 N 2 Rmag 7
6 7
6 7 6 1 7
4 q_ 7 ðtÞ 5 ¼ 6 0 0 7
6 I 7
p_ 8 ðtÞ 6 7
4 1 1 B5
C76 C7 I
2 3 213
q6 ðtÞ
6 7 6N7
4 q7 ðtÞ 5 þ 6 7
4 0 5 ½V ðtÞ ð6:26Þ
p8 ðtÞ 0
A block diagram of a series motor with armature and field voltages and current is
shown in Fig. 6.15 with ideal bond graph.
The bond graph of a network with a series motor in connection with loads
through a transformer is shown in Fig. 6.16. The input to the transformer is the
rotational effort (torque) and angular rate of motor shaft which drives the load of
mass m with damping in rope Bout. There are field Rf, armature Ra, and output
resistances R of series motor respectively, inductance at output, spring and damper
at shaft end are given by I, k, and B respectively. In this circuit an IC field represents
a capacitor and inductor at respective bonds as well as parasitic capacitance and
inductance in the windings. In IC field, effective parasitic capacitance and induc-
tance along with individual elements are represented in a vector notation given as
94 6 Advance Bond Graph Modeling
Fig. 6.15 A series motor with armature and field sources and bond graph
Fig. 6.16 A network with a series motor driving shaft and pulling a rope loaded with mass
2 3
1
a62
f6 6 I6 7 p6
¼6
4
7 ð6:27Þ
q2 1 5 q2
a26
C2
If a62 ¼ a26 ¼ 0 then the effect of parasites, windage, or air gap are ignored, and it
represents a simple inductor element at port 6 and a simple capacitive element at
port 2. Otherwise, a62 ¼ C162 is a parasitic capacitance faced at field windings from
armature and a26 ¼ L126 is a mutual inductance faced at armature due to field
windings. Both energy storing variables p6 and q2 in IC field will also be state
variables in the system equations.
6.5 Resistive Fields 95
In some literature, other single port elements with gyrators or transformers are
replaced for simplicity, e.g., a gyrator followed by an inductor can replace both
gyrator and inductor as equivalent capacitor. Also, a capacitor followed by a gyrator
is simplified as an inductor. A transformer followed by a capacitor or inductor will
respectively replace with capacitor or inductor. But these simplifications often lead
to confusion relating to second ports of the two-port elements and needs more
attention while formulating constitutive laws as well as state space equation. It is
therefore not necessary to follow these simplifications. A resistor followed by a
transformer or gyrator remains a resistor with appropriate constitutive law (P6.10).
Electrical networks or mechanical networks for fluid flow and other applications
have many similar ports representing the same constitutive laws between elements as
shown in Fig. 6.19. There may be interactions between elements that can also be
represented such as mutual inductance, two-way spring action just as C-fields and I-
fields. The constitutive laws remain the same but appear in matrix notations.
The figure shows a vector bond graph. Let us consider it is a 3 3 system for
each single port element. This implies, C, I, and R are 3 3 matrices with diagonal
elements representing values of each individual elements and non-diagonal values
represent interaction elements given as
2 3 2 3 2 3
k11 k12 k13 m11 m12 m13 R11 R12 R13
6 7 6 7 6 7
C ¼ K ¼ 4 k21 k22 k23 5 I ¼ M ¼ 4 m21 m22 m23 5R ¼ 4 R21 R22 R23 5
k31 k32 k33 m31 m32 m33 R31 R32 R33
If non-diagonal elements are zero, then diagonal elements only represent the self
values. In this system, it should be noted that the values of any element [aij]
represent the interaction of ith element with jth element in a matrix. For example,
k12 is the combined compliance of two capacitances C1 and C2 just like field
98 6 Advance Bond Graph Modeling
capacitance. The equations of the system have the same methods except vectors,
and matrices will be used. In this example, the state vector consists of two vector
* *
elements, q2 and p5 , where each element in state vector can be represented as
2 3 2 3
q21 p51
* 6 7 * 6 7
q2 ¼ 4 q22 5 p5 ¼ 4 p52 5
q23 p53
ð6:28Þ
" * #
* q2 ðtÞ
x ðtÞ ¼ *
p5 ðtÞ
The causality assignment of vector bond graph is the same; the methodology of
obtaining equations at each junction also remains the same. The first set of equa-
tions is obtained as
* * *
e3 ðtÞ ¼ e1 ðtÞ e2 ðtÞ ð6:29Þ
* *
e3 ðtÞ ¼ T e4 ðtÞ ð6:30Þ
* * *
f 4 ðtÞ ¼ f 5 ðtÞ þ f 6 ðtÞ ð6:31Þ
" * #
* * *
q_ 2
The solution to be found is x_ ¼ * , where q_ 2 and p_ 5 are two independent
p_ 5
vectors. From the figure and set of equations, we can get following equations as
* * *
q_ 2 ¼ f 2 ðtÞ ¼ f 3 ðtÞ ð6:32Þ
* * *
q_ 2 ¼ f 3 ðtÞ ¼ T 1 f 4 ðtÞ ð6:33Þ
6.8 Vector Bond Graphs 99
* * *
1
q_ 2 ¼ T f 5 ðtÞ þ f 6 ðtÞ ð6:34Þ
Where
* *
f 5 ðtÞ ¼ M1 p5 ð6:35Þ
and
* *
f 6 ðtÞ ¼ R1 e6 ðtÞ ð6:36Þ
* * *
e6 ðtÞ ¼ e4 ðtÞ ¼ T 1 e3 ðtÞ ð6:37Þ
* * *
e3 ðtÞ ¼ F ðtÞ K 1 q2 ð6:38Þ
* * *
f 6 ðtÞ ¼ R1 T 1 F ðtÞ K 1 q2 ð6:39Þ
*
Similarly for second variable vector p_ 5 we have
* * *
p_ 5 ¼ e5 ðtÞ ¼ e6 ðtÞ ð6:41Þ
*
The state vector and input F ðtÞ are 6 1 and 3 1 vectors respectively for the
overall sixth order system given in Eq. (6.43).
100 6 Advance Bond Graph Modeling
Problems
P6.1 Find the state space representation of loop circuits given in Fig. 6.20. Is
it possible to obtain different state space by reversing the loop (i.e.,
counterclockwise)?
P6.2 Draw a block diagram of the bond graph of a system given in 6.16 with given
details. Find state space representation to measure the position of mass
attached to the rope and maximum force acting on the rope.
P6.3 Formulate a state space relation of the bridge circuit given in Fig. 6.17. Is
there any difficulty in obtaining the output equation for measuring voltage
across load resistor? Is it possible to obtain the transfer function of the
system?
P6.4 Find the bond graph diagram of the following circuits in Fig. 6.21 and
explain how many choices of causalities are possible. Write state space
equations of the system of all possible choices. Are the state space repre-
sentations different or equivalent? Explain
P6.5 In order to represent state space with the two variables of systems given in
Fig. 6.6, an effort source can be represented as V ðtÞ ¼ V m sin ðωtÞ. Now
differentiation of this is possible, which may lead to state space representa-
tion with two variables. Formulate a state space representation and give a
transfer function from each state as output to the input V(t).
P6.6 Capacitor bank in parallel are used for amplification in voltage in certain
applications. A simple circuit is shown in Fig. 6.22 ith n number of capac-
itors. Derive a nth order state space representation and discuss the structure
of A and B matrices.
Problems 101
P6.7 Formulate complete state space representation of series motor from Fig. 6.16
by measuring voltage and current of the load at the output.
P6.8 Formulate the state space representation from the circuits with fields
in Fig. 6.23.
P6.9 There are three resistive elements in the Bridge circuit given in Fig. 6.17.
Find out efforts or flow measurements in each; give a matrix of mixed field
form accordingly.
102 6 Advance Bond Graph Modeling
P6.10 If capacitors in problem P6.6 are replaced with resistors and voltage across
each resistor, and flow in each branch needs to be measured. Provide a
generalized matrix to give the measurement and specify how many elements
in vector to measure effort and how many elements for measuring flow in a
mixed form causality assignment.
P6.11 Simplify the Bond Graph given in P6.1 by eliminating the bond 5 (from both
parts).
References
1. Karnopp, Dean C., Donald L. Margolis, and Ronald C. Rosenberg. 2012. System Dynamics—
Modeling and Simulation of Mechatronics Systems, 5th ed. Hoboken: Wiley.
2. 20-Sim Software. www.20sim.com.
3. Zimmer, Dirk., François E. Cellier. 2006. The Modelica multi-bond graph library. In Pro-
ceedings of the 5th International Modelica Conference, Vienna, Austria, 4–6 September 2006,
559–568.
4. Craig, K., P. Voglewede. 2010. Multidisciplinary Engineering Systems Graduate Education:
Master of Engineering in Mechatronics. In Transforming Engineering Education: Creating
Interdisciplinary Skills for Complex Global Environments, IEEE, Dublin, 6–9 April 2010, 1–14.
5. Cetinkunt, Sabri. 2006. Mechatronics, 1st ed. Somerset: Wiley.
Chapter 7
Simulation and Analysis of State Space
Systems
This equation is non-homogenous due to the term bu(t), so the solution requires
computation of the homogenous and particular solution given as
The homogenous solution xh(t) obtains by initial conditions and system is solved if
uðtÞ ¼ 0 and general homogenous solution of Eq. (7.1) is given as
The particular solution xp(t) is solved by calculus due to function u(t) in the
equation and total solution of Eq. (7.1) is given as
ðt
xðtÞ ¼ e at
xð0Þ þ e eaτ b uðτÞdτ
at
ð7:4Þ
0
We can see from Eq. (7.2) that at t ¼ 0 system starts responding from homogenous
solution and if uðtÞ ¼ 0 then only the homogenous solution contributes to the final
shape of the response. If the system is relaxed i.e., xð0Þ ¼ 0 then the system is only
excited by the given input u(t). A system excites to either non-zero initial conditions
or given input, otherwise the system remains at rest. These two excitation sources
are known as free response which are due to initial conditions only and forced
response due to given input to the system. The output state space response of the
first-order system is
8 9
< ðt =
yðtÞ ¼ c eat xð0Þ þ eaðtτÞ b uðτÞdτ þ d uðtÞ ð7:5Þ
: ;
0
_
* * *
x ðtÞ ¼ A x ðtÞ þ B u ðtÞ ð7:6Þ
* * *
y ðtÞ ¼ C x ðtÞ þ D u ðtÞ
This equation is solved by linear algebra and the general form of the solution in
vector field is
ðt
* * *
x ðt Þ ¼ e At
x ð 0Þ þ e At
eAτ B u ðτÞdτ ð7:7Þ
0
*
In this case A, B are matrices of n n and n p orders respectively, where x ðtÞ,
*
u ðτÞ are vectors of n 1 and p 1 orders respectively. The output response
solution is
8 9
< ðt =
* * * *
y ðtÞ ¼ C eAt x ð0Þ þ eAðtτÞ B u ðτÞdτ þ D u ðtÞ ð7:8Þ
: ;
0
*
In the output response, y ðtÞ is a vector of q 1 order and C, D are matrices of q n
and q p orders respectively. The general formulation of state space equation of
any order solves by integrating the parameters with given input over the time. The
simulation diagram of the system in Fig. 7.1, illustrates the output Eq. (7.8) solving
state space formulation of Eq. (7.6).
The figure shows the real-time simulation diagram in MATLAB/Simulink of
*_
state space equations; an integrator integrates the vector x ðtÞ during the simulation.
*_
This vector x ðtÞ is composed of two terms, first is the multiplication of state vector
*
A with output of integrator which is the state vector x ðtÞ and the term is input vector
*
u ðtÞ multiplied with input gain matrix B. Similarly, the output vector is composed
*
of two terms, first is the multiplication of output gain matrix C to state vector x ðtÞ
*
and a feedforward component from input u ðtÞ by multiplying feedforward (input–
output) matrix D. If the system is of first order then all matrices and vectors are
7.1 Free Response of First-Order System 105
treated as 1 1 elements to represent Eq. (7.1). This system can also be simulated
as a transfer function approach, which is only for a relaxed system and for which
initial conditions cannot be specified in simulation. In this block diagram, the
integrator takes initial conditions as a vector and starts processing from these values
to determine the total response of the system.
The free response of a system is obtained with zero input and exciting through
initial conditions only. The homogenous solution of an nth order differential
equation either in algebraic formulation or in state space formulation provides the
free response of the system. The free response of the first-order system is given as
xðtÞ ¼ eat xð0Þ and with output it is given as
If the initial condition is zero then the output remains silent for all the time but if
there is non-zero initial condition then the output will start from that value and
either increase or decrease exponentially depending upon the sign of parameter a.
In a first-order system, this parameter a or matrix [a] constitutes the only eigenvalue
106 7 Simulation and Analysis of State Space Systems
1.2
a<0
Normalized y(t)
0.8
a>0
0.6
a=0
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
time (normalized)
Fig. 7.2 Normalized responses of a first-order system with change in parameter {eigenvalue} a
0.8 2
1.5
0.6 1.5
y(t)
1
0.4 1
0.5
0.2 0.5
0 0 0
0 5 10 0 5 10 0 5 10
time (sec)
of the system. If a > 0 then the response will increase exponentially and if a < 0
then the response will decrease exponentially. The response will remain constant if
a ¼ 0 at initial condition. In a first-order case we can simply observe that output
gain matrix C just amplifies or attenuates the response. Free responses in Fig. 7.2
show the dependency on a, which is the only eigenvalue of the first-order system.
The output is normalized to maximum value in order to show the comparison of
responses for exponential growth and decay for the same eigenvalue with sign
inversion. The responses for three cases are also shown in Fig. 7.3 to show the level
of amplitude obtained in 10 s time.
7.2 Eigenvalues of Higher Order System 107
An nth order system will have n eigenvalues; the free response of such a system
depends upon all the eigenvalues. A system matrix can be constituted such that all
eigenvalues appear explicitly in state matrix, but generally this is not the case. The
eigenvalues are also called poles of the system, which are determined from the
denominator of transfer function. A transfer function is a realization of the relaxed
state space system in frequency domain as we already discussed in Chap. 1.
A system with p inputs and q outputs has p q transfer functions from each output
to each input.
*
G ðsÞ ¼ C ½s I A1 B þ D ð7:12Þ
Each transfer function shows the relationship between a single input to a single
output represented as
DðsÞ ¼ js I Aj ¼ 0 ð7:14Þ
The response of each eigenvalue is similar to the response of the first-order system.
If an eigenvalue is less than zero then it contributes a decaying response and if an
eigenvalue is greater than zero then it grows exponentially. In higher order real
systems, complex eigenvalues may appear in a conjugate pair. If DðsÞ ¼ 0 has any
complex root then its conjugate pair is also a root of the polynomial D(s). This is
due to the quadratic solution of the polynomial; a third-order polynomial may have
three real roots or a real root and two complex roots (in which both must be
conjugate to each other). In a matrix notation, eigenvalues may appear in diagonal
of [A] matrix if the matrix is diagonal or upper triangular but for real systems,
complex eigenvalues do not appear in complex form in real matrices. Finding roots
of the polynomial D(s) or obtaining eigenvalues through matrix decomposition
presents the complex eigenvalues. The response of nth order system realizes as the
combination of n first-order systems if all poles are real, whereas we can also
observe the response of nth order system as a combination of second-order system if
roots are complex conjugate pair with the remaining first-order system. A transfer
function given in Eq. (7.13) can also be written as
108 7 Simulation and Analysis of State Space Systems
Qm
s sj
G ðsÞ ¼ K Qj¼1
p : ð7:15Þ
i¼1 ðs si Þ
Where there are m zeros that are roots of N(s) and n poles that are roots of D(s).
Each real root has the simple form s ¼ si where si is a real number. A complex root
has a form
si ¼ σ þ j ω ð7:16aÞ
si * ¼ σ j ω ð7:16bÞ
Which yield as
D 2 ðsÞ ¼ s2 þ a s þ b ð7:18Þ
It is not possible to have a real system with complex coefficients in either denom-
inator or numerator as it cannot be implemented with real devices. So poles or zeros
always exist in conjugate pairs which lead to polynomials with real coefficients.
A second-order system consists of two poles and the location of these poles in
complex plane “s” describes the behavior of response. A solution of higher order
system in Eq. (7.8) may appear as a combination of two linear solutions of the first-
order system
Note here y1(t) and y2(t) are two first-order solutions due to two eigenvalues (poles)
of the second-order system λ1 and λ2. The solution will be
The values or a and b forms two types of responses, undamped and damped
responses. These are also characterized as undamped oscillation or damped oscil-
lation of the second-order system.
If a ¼ 0 then the system is undamped and two eigenvalues are conjugate pairs
given as
pffiffiffiffiffiffiffi
λ1, 2 ¼ b ¼ jω ð7:23Þ
pffiffiffi
In this case we refer ω ¼ b as undamped natural frequency of the system. The
solution of the system is given as
Solution The transfer function is obtained from MATLAB by using the following
command
A¼[0 -1;1 0]
B¼[1;1]
C¼[1 1]
110 7 Simulation and Analysis of State Space Systems
D¼0
[Num,Den]¼ss2tf(A,B,C,D)
2s
GðsÞ ¼ ð7:26Þ
s2 þ1
The system has two eigenvalues which can be found by either command roots
(Den) or eig(A) to obtain roots of characteristic equation or eigenvalues of [A]
computed as
0 þ 1i
0 - 1i
We can use initial conditions to find the values of constants in either of two
solutions. Algebraically in this problem we plug initial conditions to the equation
* *
y(t) and y_ ðtÞ or we can use free response simulation as y ðtÞ ¼ eAt x ð0Þ. We can
also directly simulate the system in SIMULINK given in Fig. 7.1 with zero input
and obtain solutions for different initial conditions given in Fig. 7.4
* *
We can see that responses two x 01 and x 02 are 90 out of phase because either
*
cos(t) or sin(t) will remain in a solution. Using initial conditions x 03 both terms
share their parts in the solution of free response. All three initial conditions yield a
sinusoidal signal in response without any amplification or attenuation over time,
this is known as undamped oscillation of the system. The undamped oscillatory
response only appears when real parts of second-order poles are zero or we can say
eigenvalues (poles) lie on the imaginary axis of s-plane.
2
x01
x02
1 x03
y(t)
-1
-2
0 1 2 3 4 5 6 7 8 9 10
time (sec)
Fig. 7.4 Free response of undamped system with different initial conditions
7.3 Free Response of a Second-Order System 111
and if parameter a 6¼ 0 then the system has damped oscillations based upon the
value of a. Let us now consider
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
a 4b a2
λ1, 2 ¼ σ jω ¼ j ð7:28Þ
2 2
The term σ ¼ a=2 is known as damping frequency and also referred as attenuation
ratio. The term ω is the frequency of the system which is measured in rad=s and related
with time period T as ω ¼ 2π=T and frequency ω ¼ 2πf . The characteristic
equations of the second-order system often realize with two parameters of damping
ratio ζ (zeta) and natural frequency ωn.
A system is underdamped with two distinct real roots because there are no imag-
inary parts in the complex roots. This case arises when
a2 4b > 0 or a2 > 4b and ζ 2 > 1 or ðζ > 1Þ.
The system will have two distinct real roots given in Eq. (7.22)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
a a2 4b
λ1, 2 ¼
2
112 7 Simulation and Analysis of State Space Systems
In this case we may get both positive, both negative, or a positive and negative root
but for overdamped systems both roots are negative. With negative roots as
t ! 1, yðtÞ ! 0 so quickly that it does not allow the system to oscillate and
there is no damped natural frequency in the system. If roots are very small as
compared to sluggish model parameters (a, b), there may be overshoot or under-
shoot depending upon the initial conditions of the system.
Example 7.2 Overdamped Oscillations
Find the free response of the following system with given initial conditions and
explain why this is an overdamped system.
*_ 1 0 * 1
x ðtÞ ¼ x ðt Þ þ uð t Þ
0 2 1
*
yðtÞ ¼ ½ 1 1 x ðtÞ
* 1 * 1 * 1 * 1
x 01 ð0Þ ¼ x 02 ð0Þ ¼ x 03 ð0Þ ¼ x 04 ð0Þ ¼
1 1 1 1
2s þ 3
G ðsÞ ¼ ð7:32Þ
s2 þ 3s þ 2
A second-order system is said to be critically damped if it has two real double roots
and parametrically a2 4b ¼ 0 or a2 ¼ 4b. This case deals with a class of second-
order system with unity damping ratio ζ 2 ¼ 1 or ðζ ¼ 1Þ. The general solution to
this type of system with both same roots λ ¼ a=2 is given as
This system will lead to a very small constant value because of term t e =2 t and
a
2
x01
1.5 x02
1 x03
Response y(t)
0.5 x04
-0.5
-1
-1.5
-2
0 0.5 1 1.5 2 2.5 3 3.5 4
time (sec)
The zero-poles-gain form of transfer function also shows that there is one zero at
s ¼ 32 and two poles at s ¼ 1, 1. This form can be obtained in MATLAB by
the following command:
sys
sys¼ss(A,B,C,D)
sys1¼zpk(sys)
A second-order system may have both complex roots, as we consider in Eq. (7.26).
The system will oscillate like undamped systems but the difference is the magni-
tude of damping will decrease with time due to a damping ratio ζ. The case arises
when ða2 4bÞ < 0 or ζ 2 < 1 or ð0 < ζ < 1Þ. The eigenvalues are represented as
qffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
a 4b a2
λ1, 2 ¼ ζωn jωn 1 ζ ¼ 2
j ¼ σ jω ð7:35Þ
2 2
7.3 Free Response of a Second-Order System 115
In each of the solutions there are two constants {(c1, c2), (k1, k2), (K, δ),} which can
be found from the given initial conditions. (c1, c2) and (k1, k2) are simple amplitude
gain constants whereas K and δ are amplitude and phase difference (leading or lagging
depending upon the sign of phase difference) respectively. For underdamped systems
as t ! 1, the response is yðtÞ ! 0 after oscillating in sinusoidal fashion.
Example 7.4 Underdamped Oscillations
Find the free response of the following system with given initial conditions and
explain why this is an underdamped system.
*_ 1 0 * 1
x ðt Þ ¼ x ðtÞ þ uðtÞ
0 2 1 :
*
yðtÞ ¼ ½ 1 1 x ðtÞ
* 1 * 1 * 1 * 1
x 01 ð0Þ ¼ x 02 ð0Þ ¼ x 03 ð0Þ ¼ x 04 ð0Þ ¼
1 1 1 1
Solution Figure 7.7 shows that initial conditions vary in amplitude of response as
well as phase with respect to zero. These responses die down to zero as t ! 1 while
2 x01
x02
x03
1
Response y(t)
x04
-1
-2
0 0.5 1 1.5 2 2.5 3 3.5 4
time (sec)
oscillating at frequency ω. The first and second initial conditions are both positive
or negative so these constitute higher initial magnitude of response while 180 out
of phase from each other, similarly third and fourth initial conditions have mixed
signs and so these have lesser amplitude but also 180 out of phase with each other.
1 1
τ¼ or ð7:39Þ
σ λ
Each eigenvalue of a damped system will lead the system to zero; the higher the
magnitude the faster the system reaches zero or simply higher value of negative real
part, faster the system reaches to steady-state or zero.
The combination of damping ratio and natural frequency determines the oscillatory
behavior of the system. The roots (eigenvalues) of a second-order system given in
Eq. (7.27) are
7.5 Damped Natural Frequency 117
N ðsÞ
G ðsÞ ¼ ð7:42Þ
s2 þ 2ζωn þ ωn 2
2
x01
1.5
x02
1 x03
x04
0.5
y(t)
-0.5
-1
-1.5
-2
0 0.5 1 1.5 2 2.5 3 3.5 4
time (sec)
2
Overdamped
1.5 Critically Damped
Underdamped
+ Envelope of underdamped
1
- Envelope of underdamped
0.5
y(t)
-0.5
-1
-1.5
-2
0 0.5 1 1.5 2 2.5 3 3.5 4
time (sec)
Fig. 7.9 Free response of overdamped, critically damped and under damped {with envelopes on
starting from positive or negative initial conditions} systems
observe that at ω ¼ ωn , there is a peak and the amplitude of this peak depends upon
the damping ratio ζ; the higher the ζ, the lower the amplitude of |G( jω)|. This
corresponds to Fig. 7.9 that higher ζ will underdamp the system and lower ζ will
tend towards overdamping. In other part of the Bode diagram which shows the
relationship of ∠ðGðjωÞÞ with frequency, a steep shift from 0 to 180 occurs at
ωn for ζ ¼ 0 and as ζ ! 1, the curve of ∠ðGðjωÞÞ becomes more smooth and the
phase of 180 occur at higher ω values. Bode plots of closed loop feedback
system dynamics helps to design PID (or lead-lag) controllers using the phase
crossover and gain crossover margins (amplitudes).
Closed loop systems have to exhibit as a damped system, but open loop or not
controlled systems may have damping behavior which is characterized as stable or
unstable behavior due to some sort of internal amplifications. The poles of transfer
function or eigenvalues of state matrix may be such that it produces negative
damping or amplification. A second-order section of a transfer function can have
four choices of roots like damping system; a set of two distinctive roots, a real
double root or a conjugate pair with or without real part. In case of amplification we
have ζ < 0 but this parameter ζ is known as a damping ratio and so it cannot be
attributed to amplification. A damped system is a real system with stable behavior
and an amplified system is a real system with unstable behavior so these terms are
7.6 Amplification Not Damping 119
not used for describing amplification systems. The behavior of these types of
systems remains similar but instead of damping down, the amplitude of response
will be amplified. Now let us consider a second-order system with different
combination of poles.
A system will be undamped exactly like given in Eq. (7.25) and Fig. 7.10 if there is
no damping terms ζ or a ¼ 0 or σ ¼ 0. The response in Fig. 7.10 is classified as
undamped and unamplified response. Now there are three other choices of roots
with σ 6¼ 0. A system has two distinct roots of a2 > 4b, a system may have a real
Bode Diagram
10 ζ=0.0
Magnitude (dB)
ζ=0.1
0
ζ=0.2
-10 ζ=0.3
ζ=0.4
-20 ζ=0.5
ζ=0.6
0
ζ=0.7
-45
Phase (deg)
ζ=0.8
ζ=0.9
-90
ζ=1.0
-135
-180
10-3 10-2 10-1 100 101
Frequency (Hz)
double root if a2 ¼ 4b i.e., no imaginary part and a system may have a conjugate
pair if a2 < 4b.
The response shown in Fig. 7.11 of Table 7.2 is similar to responses of Table 7.1
but the difference is amplitude of response increases with time whereas, in a
damped system it decreases with time.
1
z =0
0.8
0.6
0.4
normalized y(t)
z =1
0.2
-0.2
-0.4
-0.6
-0.8
-1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
time (normalized)
0.8
0.6
Normalized y(t)
0.4
0.2
-0.2
-0.4
-0.6
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
time (normalized)
Poles in s-plane
Imaginary Axis
Real Axis
Fig. 7.13 Location of poles in s-plane; poles exist in conjugate pairs except on real axis.
response depicting the stable and unstable nature of the system. The position of
imaginary axis shows the oscillating frequency of the system and its amplitude
either grows out if it is on the right side or damped down if it is on the left side of the
s-plane. As we know, complex roots exist in pairs and so both roots have the same
response starting from the same initial conditions on both sides of the real axis.
Same responses can be seen at same distant locations of imaginary axis in upper and
lower sides of real axis. The upper and lower sides of real axis have nothing to do
with the stability of the system, they just describe the amount of oscillatory nature
present in the system. If a complex root z ¼ a þ jb has jaj
jbj then there is quicker
7.8 Forced Response 123
damping (or amplification) and lesser oscillations observed over a period of time.
On the other hand if b
jaj then oscillations are more visible in the response and
envelope of response on both sides of s-plane show the speed of convergence or
divergence of the system. A conjugate pair of poles on imaginary axis shows an
undamped behavior with a constant amplitude sinusoidal behavior at different
values of ω. In simple words, internal stability of the system depends upon the
location of the pole in the right or left half of the plane. A mirror image of the pole
in the left or right half plane stabilizes or destabilizes the system as given in
response in Fig. 7.14. A pole in the left half of plane makes the system oscillate
with a frequency ω and dampened down towards zero; on the other hand a pole of
same magnitude with an opposite sign of real part makes the system oscillate at
same frequency and grows the response out of bound.
In Eq. (7.7) we represented the response of the system in homogenous terms and
non-homogenous terms
124 7 Simulation and Analysis of State Space Systems
ðt
* * *
x ðt Þ ¼ e At
x ð 0Þ þ e At
eAτ B u ðτÞdτ ð7:45Þ
0
*
Where the term eAt x ð0Þ determines the free response and the internal stability of
the system. The same equation is used in output equation given as
8 9
< ðt =
* * * *
y ðtÞ ¼ C eAt x ð0Þ þ eAðtτÞ B u ðτÞdτ þ D u ðtÞ: ð7:46Þ
: ;
ðt 0
*
The term eAðtτÞ B u ðτÞdτ is a non-homogenous term and it determines the
0
*
forced response of the system due to input vector u ðtÞ. Considering a relaxed
*
system i.e., x ð0Þ ¼ 0, the only excitation is due to the input given to the system
which causes its effect in output response of the system given by
ðt
* * *
y ðtÞ ¼ C eAðtτÞ B u ðτÞdτ þ D u ðtÞ ð7:47Þ
0
The relation between each element of input vector to each element of output vector
is also given by a transfer function. Each element of input contributes to the system
and the response then varies accordingly. In order to understand the type of
response excited by different inputs, we now consider a stable single input–single
output system. If we compare free response term with forced response term, the
parameter eAt as solution of differential equation is common in both parts. We have
understood that the stability of the systems depend upon the eigenvalues of the
A matrix and so if the system is unstable internally then the force response will also
be unstable (if it is not being controlled specifically). There are three responses
which are often studied:
• Impulse response
• Step response
• Sinusoidal response or decaying sinusoidal response
*
We now simulate Fig. 7.4 but instead of taking uðtÞ ¼ 0, we now take x ð0Þ ¼ 0
in integrator.
where
(1 )
=k ataþk
f k ð t aÞ ¼ ð7:49Þ
0 otherwise
ð
1
δ ð t aÞ ¼ 1 ð7:50Þ
0
We can also note that Laplace transform of impulse function or Dirac’s Delta
Function is given as
So the impulse response brings the response of internal dynamics of the system at
the output. The Dirac’s Delta function has digital equivalent known as Kronecker
delta in digital control systems.
Example 7.5 Simulating Impulse Response:
Consider the Examples 7.1–7.4 and plot impulse responses for relaxed systems.
Solution A matrix is different in all examples for undamped (UND), overdamped
(OD), critically damped (CD), and underdamped (UD) systems. We simulate the
system using the Simulink diagram on Fig. 7.15 for relaxed systems.
" # " #
0 1 1 0
AUND ¼ AOD ¼
1 0 0 2
" # " #
1 1 1 1
ACD ¼ AUD ¼
0 1 1 1
126 7 Simulation and Analysis of State Space Systems
Impulse Response
1
Undamped
0.8 Over Damped
Critically Damped
0.6 Under Damped
0.4
Response y(t)
0.2
-0.2
-0.4
-0.6
-0.8
-1
0 1 2 3 4 5 6 7 8 9 10
time (sec)
Fig. 7.16 Impulse response of Undamped, Overdamped, Critically Damped and Under Damped
System
Impulse response is created by using “Signal Builder,” block for a Dirac’s delta
function (Fig. 7.16).
We note that impulse response shows the internal dynamics of the system as
shown in Fig. 7.17. The impulse response of unstable systems will also show the
amplification in the magnitudes growing out of bound. A system hit by a huge
impulse force excites the internal dynamics of systems which later settle down just
like excitation with initial conditions.
7.8 Forced Response 127
Step Response
3
2.5
1.5
1
y(t)
0.5
-0.5
-1 Undamped
over damped
-1.5 Crtically damped
Under damped
-2
0 1 2 3 4 5 6 7 8 9 10
time (sec)
Fig. 7.17 Step response of an undamped, overdamped, critically damped, and underdamped
system
A step response of the system is evaluated after giving a constant step input usually
starting at t ¼ 0. A step function is given as
( )
1 t ti
uð t t i Þ ¼ ð7:54Þ
0 otherwise
1
Lfuðt aÞg ¼ eti s ð7:55Þ
s
The transfer function of the system with a step input ðt aÞ can be represented as
GðsÞ ti s
Y ðsÞ ¼ G ðsÞ U ðsÞ ¼ e ð7:56Þ
s
The output y(t) is like a free response of a system with an additional pole at s ¼ 0
starting with a delay of ti seconds. The parameter eti s in the transfer function
represents the delay in the system. The application of a step response analysis arises
128 7 Simulation and Analysis of State Space Systems
0.5
y(t)
-0.5
ωn =ωE /2
ωn =ωE
-1
ωn 0= 2ωE 2 4 6 8 10 12 14 16 18 20
time (sec)
when a system is hit by a constant force over a period of time of operation like
combustion, constant power supply, constant force exerting on a spring, etc.
Example 7.6 Simulating Step Response
Find a step response of the all four systems in previous Example 7.5 with step
input uðt 1Þ and discuss responses.
Solution A step function of MATLAB Simulink can be used in Fig. 7.15 to
determine the step responses of the undamped, overdamped, critically damped,
and underdamped systems. Figure 7.18 shows the responses of the systems.
An undamped system oscillates with twice the magnitude as three poles now lie
on y-axis starting at ti ¼ 1 seconds. The remaining three systems settle at a constant
value but not at zero steady-state. The constant input makes the system constant at
non-zero position depending upon the location of poles in s ‐ plane. An
underdamped system may get lower magnitude while settling close to zero when
the real part is smaller, on the other hand if real part is larger than it may also settle
at a higher constant value.
The dynamic response of the system with a forcing sinusoidal function as an input has
a variety of parameters to be discussed. An input with a frequency component excites
7.9 Resonance 129
the system and the system has a combination of frequencies in its output response. A
forced input with magnitude F and excitation frequency ωE can be given as
We can observe that the system will have two additional poles on imaginary axis at
jωE and the combination of these poles with poles of G(s) affects the response of
the overall system. A transfer function G(s) has natural frequency ωn and the
location of poles due to ωn and jωE defines the response of the system. A
second-order system may be easily studied in order to check the different responses.
7.9 Resonance
When excitation frequency approaches the natural frequency i.e., ωE ! ωn then the
magnitude of the system approaches infinity. An undamped system operating at ωE
¼ ωn becomes an unstable system because of double pole at the same location on
the imaginary axis. Remember these poles are in conjugate pair, a double pole on
imaginary axis in upper half plane {above x ‐axis} and a double pole on imaginary
axis in lower half plane {below x-axis} of s-plane. A double pole on the imaginary
axis is an unstable system, so the undamped system becomes an unstable system
when excited at the natural frequency. But when an undamped system is excited by
sinusoidal input with a different frequency response y(t) contains both frequencies:
natural frequency ωn as well as excitation frequency ωE. The magnitude of the
response also varies due to variation in ωE (Fig. 7.19).
The responses of overdamped and critically damped systems are similar under
excitation from sinusoidal frequency because the poles of both systems do not have
imaginary parts. The system oscillates at the excitation frequency ωE after an initial
overshoot in magnitude. The normalized response plot of overdamped and critical
damped system shown in Fig. 7.20 depicts the initially amplified cycle and then
almost constant amplitude for both systems in both excitation frequencies ωE and
2 ωE. There is no concept of resonance in overdamped and critically damped
systems because of damping and no dominating natural frequency component in
the poles. The response of an underdamped system creates a scenario of practical
resonance because the real systems have some damping in them. If a real part of
130 7 Simulation and Analysis of State Space Systems
0.8
0.6
0.4
Normalized y(t)
0.2
-0.2
-0.4
-0.6
Over damped
-0.8
Crtically damped
Over damped 2 ωE
-1
0
Crtically 1
damped 2 ωE2 3 4 5 6 7 8 9 10
time (sec)
Fig. 7.19 Excitation response of overdamped and critically damped systems with frequency ωE
ωE = ω ωE < ω ωE > ω
1 1 1
σ=ω
0 0 0
-1 -1 -1
0 0.5 1 0 0.5 1 0 0.5 1
1 1 1
σ>ω
0 0 0
-1 -1 -1
0 0.5 1 0 0.5 1 0 0.5 1
1 1 1
σ<ω
0 0 0
-1 -1 -1
0 0.5 1 0 0.5 1 0 0.5 1
time (Normalized)
Fig. 7.20 Excitation response of an underdamped system with frequency ωE
7.10 Decaying Sinusoidal Response 131
ωE = ω
1
Response of Undamped System y(t)
-1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ωE < ω m=-2
1 m=-1
m=-0.5
0
-1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ωE > ω
1
-1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
time (Normalized)
Fig. 7.21 Excitation response of an undamped system with decaying sinusoidal input
eigenvalues of A matrix is very small than imaginary part, it means there is less
damping in the system. This system may also resonate or oscillate with high
amplitude when excited by frequency equal to the imaginary part {natural fre-
quency} of the system. Due to damping effect, the system settles at oscillating
frequency after following a transient response in the first few cycles (Fig. 7.21).
Where m < 0 is the decaying ratio, and if m ¼ 0, the input signal becomes pure
sinusoidal and for m > 0 the input signal will be amplified with time; ωE is the
excitation frequency. Figure 7.22 shows that the undamped system settles at
oscillating frequency after transient response and the higher the magnitude of
m, the quicker the settling time (Fig. 7.23).
132 7 Simulation and Analysis of State Space Systems
small ωE Large ωE
1 1
Over Damped
0.5 0.5
0 0
-0.5 -0.5
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
m=-2
Critically Damped
1 1 m=-1
m=-0.5
0.5 0.5
0 0
-0.5 -0.5
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
time (normalized) time (normalized)
Fig. 7.22 Excitation of critically damped and overdamped system with decaying sinusoidal input
ωE = ω ωE < ω ωE > ω
1 1 1
σ=ω
0 0 0
-1 -1 -1
0 0.5 1 0 0.5 1 0 0.5 1
1 1 1
σ>ω
0 0 0
-1 -1 -1
0 0.5 1 0 0.5 1 0 0.5 1
1 1 1
σ<ω
0 0 0
-1 -1 -1
0 0.5 1 0 0.5 1 0 0.5 1
time (normalized)
Step Response
Impulse Response
1
3
0.8
Amplitude
0.6
2.5
0.4
0.2
0 2
0 2 4 6 0 2 4 6 8 10
Time (sec) Time (sec)
Sinusoidal Response
Sinusoidal Unbounded response
4 1
Decaying sinusoidal
2 0
0 -1
-2 -2
-4 -3
0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1
Time (sec) Normalized Time
Fig. 7.24 Excitation responses of impulse, step, sinusoidal, and decaying sinusoidal functions as
bounded input and a sinusoidal function as unbounded input for system in Eq. (7.71)
134 7 Simulation and Analysis of State Space Systems
There are different variants of responses due to the combination of slower and
faster eigenvalues of the system. Eigenvalues farther from origin on the x-axis are
faster {large value of damping or amplification} and those closer to origin are
slower to settle. It should also be noted that the third- and higher order systems
usually have a combination of different poles and thus the response of linear
systems also adds up accordingly.
A system is internally stable if the eigenvalues are in the left half plane, but a
stability related with input of the system is known as external stability. This
stability is related with forced response and also known as Bound Input Bounded
Output {BIBO} stability of the system. This stability determines the response of the
system with bounded input. The forced response of the system is given as
ðt
* * *
y ðtÞ ¼ C eAðtτÞ B u ðτÞdτ þ D u ðtÞ ð7:61Þ
0
ðt
* * *
y ðtÞ ¼ gðt τÞ u ðτÞdτ þ D u ðtÞ ð7:64Þ
0
A bounded input will yield response as
ðt
*
y ðtÞ ¼ gðt τÞ Udτ þ D U ð7:65Þ
0
7.11 External Stability 135
*
As D < 1 is a real bounded matrix so the response y ðtÞ will be bounded if
ðt
* * *
y ðtÞ ¼ gðt τÞ u ðτÞdτ þ D u ðtÞ < 1: ð7:66Þ
0
*
As we know that both u ðτÞ and D are bounded so
ðt
* * *
y ðtÞ ¼ gðt τÞ u ðτÞdτ þ D u ðtÞ ð7:67Þ
0
ðt
*
y ðtÞ ¼ gðt τÞ Udτ þ D U ð7:68Þ
0
t
ð
*
y ðtÞ ¼ jU j gðt τÞ dτ þ jD U j ð7:69Þ
0
t
ð
*
So y ðtÞ will be bounded if and only if gðt τÞ dτ is bounded
0
t
ð
*
y ðtÞ < 1 iff gðt τÞ dτ < 1:
ð7:70Þ
0
_
* 1 0 * 1
x ðtÞ ¼ x ðtÞ þ uð t Þ
2 1 1 ð7:71Þ
*
yðtÞ ¼ ½ 1 0 x ðtÞ þ ½2 uðtÞ
We can see that unstable pole has been canceled with a zero and the remaining
system is actually a stable first-order system. If we observe C matrix which shows
that only first state of the system is contributing in the output and first state of the
system is an independent stable eigenvalue of the system {decomposition of states}.
This system is internally unstable but externally stable when a bounded input is
given to the system. We now plot the bounded responses of the system with respect
to four bounded inputs with an unbounded output. Figure 7.25 shows the
corresponding responses computed from the following MATLAB commands
Total Response
2
Free Response
Forced Response
1.5 Total Response
1
y(t)
0.5
-0.5
0 2 4 6 8 10
time (sec)
Fig. 7.25 Total response of a system as combination of free response and forced response
7.12 Total Response of the System 137
*
If a system is not relaxed then x ð0Þ 6¼ 0 and given initial condition excites the free
*
response of the system, a non-zero input vector u ðtÞ excites the forced response of
the system. Total response is the linear combination of both responses, a superpo-
sition of homogenous and particular solution of differential equation. The final
stability is the combination of internal and external stability of the system and
profile of response also depends upon combination of both. We studied four types
of free responses of second-order system, and four types of bounded input forced
responses, which show that different combinations of responses are possible
depending upon the location of eigenvalues as well as the type of input given.
The same eigenvalues of the system or poles of characteristic contribute to free and
forced response so usually there are similar types of responses that combine to
produce an overall effect. Some examples are the combination of a resonant forced
response and undamped system, the combination of an overdamped system with
decaying sinusoidal input, etc.
Example 7.8 Simulating Total Response
Simulate and plot the free response, forced response, and total response of the
system
1 2 1
A¼ , B¼ , C ¼ ½ 1 1 , D ¼ ½1 ð7:74Þ
2 1 1
* 1
x ð 0Þ ¼ ð7:75Þ
1
uðtÞ ¼ et sin ð2tÞ ð7:76Þ
Solution The total response of the system is a combination of free response and
forced response. Simulating Fig. 7.1 with plant dynamics given in Eq. (7.73), initial
conditions given in Eq. (7.74) and decaying sinusoidal input given in Eq. (7.75),
free response and forced response given in Fig. 7.4 and Fig. 7.15 respectively yields
the total response by giving non-zero initial conditions in Fig. 7.15 (modified
diagram of Fig. 7.1). We can observe that free response from magnitude of 2 for
y(t) and zero magnitude in forced response. The total response also starts with
magnitude 2 and settles earlier due to out of phase free and forced responses
(Fig. 7.26).
Problems 139
0.5
Response y(t)
-0.5
-1
0 1 2 3 4 5 6 7 8 9 10
time (sec)
Fig. 7.26 Comparison of responses for the eigenvalue of same magnitude and frequency in right
and left half of the s-planes
Problems
P7.1 Complex eigenvalues or complex poles always appear in pairs but do not
appear in matrix A or coefficients of characteristic equation DðsÞ ¼ 0. Ana-
lyze third- and fourth-order polynomial equations and show that possibilities
of different complex roots {poles} in both equations. How you relate this to
state matrix A?
P7.2 A harmonic oscillator, Equation of a mass spring system without a damper are
given as m€yðtÞ þ kyðtÞ ¼ 0. Find a solution of this oscillator using Laplace
transform and find its natural frequency. Represent this system in state space
formulation and simulate it using the following parameters: 100 N force pulls
the system for 0.995 m in relaxed position. What will be the response of the
system for position and velocity if we pull the system to following situations
and release? {Use g ¼ 9:81 m=s2 }
(a) 0.25 m only with negligible speed
(b) 0.25 m with 1 m/s
P7.3 A harmonic oscillator will keep on moving forever if there is no damping in it
but all real systems have damping and in the presence of damping the equation
of motion is given as
P7.5 Find the transfer function of the system given in P2, P3, and P4 using
MATLAB. Find the natural frequencies, damping ratios of all systems and
use MATLAB for Bode plot of all the systems in one figure. Compare
damping ratios of the system in bode plot and relate with responses given in
Fig. 7.8. What happens if the numerator of the system changes but the
denominator remains the same in Bode plot?
P7.6 Find transfer functions and plot the simulated free responses of the following
system with different initial conditions. Also plot the envelopes of oscillatory
systems and find out how the initial or maximum point of the envelope is
calculated.
(a)
_
* 0 2 * 1
x ðt Þ ¼ x ðt Þ þ uð t Þ
2 0 2
*
yðtÞ ¼ ½ 2 1 x ðtÞ
* 1 * 0 * 1
x 01 ð0Þ ¼ x 02 ð0Þ ¼ x 03 ð0Þ ¼
0 1 1
(b)
_
* 2 0 * 2
x ðt Þ ¼ x ðt Þ þ uð t Þ
0 3 1
*
yðtÞ ¼ ½ 1 2 x ðtÞ
References 141
* 10 * 1 * 2 * 2
x 01 ð0Þ ¼ x 02 ð0Þ ¼ x 03 ð0Þ ¼ x 04 ð0Þ ¼
1 10 1 1
(c)
_
* 2 1 * 1
x ðt Þ ¼ x ðt Þ þ uð t Þ
0 2 1
*
yðtÞ ¼ ½ 1 1 x ðtÞ
* 1 * 1 * 1 * 1
x 01 ð0Þ ¼ x 02 ð0Þ ¼ x 03 ð0Þ ¼ x 04 ð0Þ ¼
1 1 1 1
(d)
_
* 2 3 * 1
x ðt Þ ¼ x ðt Þ þ uð t Þ
3 2 1
*
yðtÞ ¼ ½ 2 2 x ðtÞ
* 2 * 1 * 1 * 2
x 01 ð0Þ ¼ x 02 ð0Þ ¼ x 03 ð0Þ ¼ x 04 ð0Þ ¼
1 2 2 1
Normalize the time and response amplitude to the maximum values and compare
your plots without normalization; why does normalization sometimes produce
better understanding? Is there any other parameter which may be used as a
normalization factor for response? Find the eigenvalues of all systems and find
out whether the systems are stable, unstable, or marginally stable based upon
eigenvalues and free responses.
References
1. Kreyszig, Erwin. 2006. Advanced Engineering Mathematics, 9th ed. Hoboken: Wiley.
2. Karnopp, Dean C., Donald L. Margolis, and Ronald C. Rosenberg. 2012. System Dynamics—
Modeling and Simulation of Mechatronics Systems, 5th ed. Hoboken: Wiley.
3. Friedland, Bernard. 2005. Control System Design—An Introduction to State Space Methods.
New York: Dover.
4. Chen, Chi-Tsong. 1999. Linear System Theory and Design, 3rd ed. New York: Oxford Univer-
sity Press.
Chapter 8
Introduction to Control Systems
Resource Objectives
Actuation System
Controller
Steering
Measurement
Sensors
Computation Guidance
from the behavior of plant m(t) may be different than what we actually want to
achieve as objectives y(t). Reference output is the desired output d(t) at sensors with
the error calculated between actual measurements and desired (reference) outputs
as error signals e(t), which become input to control systems. We can summarize the
input and output of the closed loop system as
Inputs: There are two types of inputs to the closed loop systems.
Controllable inputs u(t): An output of a control system that drives the plant.
Exogenous inputs: These are either reference inputs or disturbances in the system
that are due to external sources or internal dynamics, e.g., delays or approximation
8.1 Representation of Controller 145
errors of modeling, etc. These inputs are not controllable by the control system and
add the burden to achieve desired objectives in required time.
Reference inputs: These inputs are the desires for the system; a control system tries
to achieve these but cannot control these inputs, which is why these are also termed
as exogenous inputs. Reference input r(t) and reference output d(t) are both
reference inputs to the closed loop feedback system.
Outputs: There are two types of outputs of the plant, one that can be measured
through sensors known as measurements m(t), and the other, which are desired
objectives of the plant y(t). In a simplified system, measurements are usually taken
as desired objectives which yield yðtÞ ¼ mðtÞ.
U ðsÞ
K ðsÞ ¼ ð8:1Þ
Eð s Þ
But we also know that each transfer function can also be represented as a state space
formulation, and so K(s) can also be represented as:
_
* * *
x c ðtÞ ¼ Ac x c ðtÞ þ Bc e ðtÞ
K ðsÞ ¼ * ð8:2Þ
* *
u ðtÞ ¼ Cc x c ðtÞ þ Dc e ðtÞ
Here, subscript “c” denotes the controller or compensator in the state vector and
state space matrices. If dðtÞ ¼ 0 then measurements are directly fed into the control
system to get the compensation commands i.e., u(t). The simulation block diagram
in Fig. 8.3 shows that a state space plant can be represented as an integrator with
matrices as gains. The objectives are required outputs y(t) with C as output gain
matrix and D as input–output gain matrix or feedforward matrix. Measurements
may be different than the actual required outputs; combinations of state vector and
input vector to represent measurement are given by Cm and Dm matrices. A
controller may be represented as a transfer function as well, such as a PID block
in Fig. 8.3. PID is one of the most commonly used and technologically accepted
controller design methods.
146 8 Introduction to Control Systems
Fig. 8.3 Simulation diagram of control system as PID controller with different measurements and
objectives. Actuators and sensor measurement transfer functions are taken as unity with no
reference input
In linear systems it is also very common to write model equations in terms of the
error dynamics by defining state vector appropriately. A stable linear system will
settle at zero steady-state if initial conditions are zero and there is no constant term
being added. The error model shows that error goes to zero during the process. In
this case zero is the desired reference value which can eliminate the signal d(t) from
the simulation. The error model of the plant dynamics are given with the same state
space matrices while state variables are defined as the error in the states.
Example 8.1 Linear Model of an Inverted Pendulum
State space model of a simple pendulum (Example 2.1) is given as
θ_ ¼ θ
ð8:3Þ
€θ ¼ g sin θ
l
If we assume θref ¼ π2, write a state space formulation in the error dynamics
Solution Let
x1 ¼ θ θref ð8:4Þ
8.3 Estimator or Observer 147
So
x_ 1 ¼ θ_ ¼ x2 ð8:5Þ
x_ 1 ¼ x2
g ð8:6Þ
x_ 2 ¼ cos x1
l
A set of equations that estimates the state vector during process is called an
estimator or observer. In digital control theory the word “estimator” is more popular
or appropriate and in continuous time system “observer” is commonly used. In
reality all states of the system are not explicitly available because taking measure-
ments of each and every state through sensor is almost impossible. So in order to
effectively take control of all states it is required to observe states through some
mechanism. The first question arises that if we know model, initial conditions, and
starting time as well as controlled input to plant we can simultaneously run this
model to get the values of these states. True, these are called open loop observers as
given in Fig. 8.4.
State space formulation for open loop observer is given as
* *
x_ e ðtÞ ¼ A xe ðtÞ þ B uðtÞ
ð8:7Þ
* *
ye ðtÞ ¼C xe ðtÞ þ D uð t Þ
* *
Observer output ye ðtÞ should reach the actual output of plant y ðtÞ simultaneously as
both systems are starting at the same initial conditions as well as excited with the
same input vector. In reality it is not possible due to the fact that actual plant, which
also have uncertainties, facing environmental disturbances, approximation errors
and transmission delays, etc. does not behave exactly as same as the analytical open
148 8 Introduction to Control Systems
Fig. 8.4 Open loop observer to estimate the states of the plant simultaneously by taking same
inputs and same state matrices
loop observer. A closed loop system is required to get the actual states of plant and
passes to the observer, which estimates the correct values of states instantaneously
for regulation or control commands. A closed loop observer needs measurements or
error with respect to measurements and input calculated by controller to give the
required estimates of the states.
A closed loop observer shown in Fig. 8.5 is called closed loop due to the
feedback loop from the error between outputs of the plant and the observer. The
task of the closed loop observer is to minimize the error between actual outputs of
* *
the plant and observed outputs and as eðtÞ ! 0 then y o ðtÞ ! y ðtÞ, thus matching
the observer with plant. The equation of the closed loop observer is given as
* *
n *
o
* *
x_ o ðtÞ ¼ A xo ðtÞ þ B u ðtÞ þ L y ðtÞ yo ðtÞ
ð8:8Þ
* * *
yo ðt Þ ¼ C xo ðtÞ þ D u ðtÞ
*
An input vector u ðtÞ may be composed of feedback controller output and reference
input. The output of closed loop observer regulates the plant, which is actually a
controller gain matrix multiplying either states or estimates of states. A control
system is actually a combination of both closed loop estimator and regulator to
provide controlled input. This is also referred to as compensator—a combination of
regulator and estimator as given in Eq. (8.2) and below:
8.3 Estimator or Observer 149
Fig. 8.5 A closed loop observer to estimate the real states by minimizing error between actual and
estimated output with same input source u(t)
_
* * *
x c ðtÞ ¼ A c x c ð t Þ þ Bc m ð t Þ
ð8:9Þ
* * *
u ðtÞ ¼ Cc x c ðtÞ þ Dc m ðtÞ
M ðsÞ
G ðsÞ ¼ ð8:10Þ
U ðsÞ
150 8 Introduction to Control Systems
The controllers are given as transfer function output of the control system to the
input of the controller which is also the output of plant
U ðsÞ
K ðsÞ ¼ ð8:11Þ
M ðsÞ
It must be clear that both K(s) and G(s) are different transfer functions where the
order of denominator is always greater and equal to the order of numerator and
these are not reciprocal to each other. Furthermore, observer can also be classified
as full state observer and reduced order observer. A full order observer estimates all
states of the plant while the reduced order observer only estimates those states that
are not being measured through sensors.
Linear control theory describes the details of controller design, compensator design,
and full/reduced order observers design for robust and optimized performance of
the plant. A simple controller may be a PID controller which is designed to stabilize
the plant with specific parameters. A PID controller is equivalent to lead-lag
compensator design mathematically and both are basically second-order systems.
The performance specifications of second-order systems are applicable to the
designs in frequency domain for the control systems.
In this chapter, we will explore the simple PID control system design to understand
the basic performance specifications of the closed loop system. Figure 8.7 shows
the typical second-order transient response for a variable y(t) to achieve the final
desired value yf(t).
This is the time it takes for the response to decay and remain within the specified
tolerance band of the steady-state value. It is denoted by ts and the value is given
when the response enters the tolerance band and remains in the band afterwards. An
usual approximation of settling time for a second-order damped system is given in
terms of damping frequency σ or damping ratio ζ and natural frequency ωn.
4:6 4:6
ts ¼ ð8:12Þ
σ ζωn
This is the time required to raise the response from 10 % of final value to 90 % of the
final value. Some text give from 0 to 90 % of time but this can only be a true
measure if there is no undershoot in the response. The rise time is approximated by
using natural frequency of second-order system given as:
1:8
tr ð8:13Þ
ωn
152 8 Introduction to Control Systems
8.4.1.4 Overshoot
This is the relative measure that system raised above the final value before decaying
and settling down. It is given as the percentage of peak relative error (maximum)
and steady-state value
ðymax yss Þ
Mp ¼ 100 % ð8:14Þ
yss
8.4.1.5 Undershoot
This is the value of a response when a system initially goes in the opposite direction
of the required value. It is the opposite of overshoot and given as:
ðyss ymin Þ
Undershoot ¼ 100 % ð8:16Þ
yss
In many cases undershoot does not occur, especially if all initial conditions are
either positive or negative. The definition of responses can be a little different if a
response starts from a higher value and settles at lower value. In MATLAB, a
command S¼stepinfo(y,t,yfinal) calculates the performance parameters of the
required functions.
Designing a control system from a model is often a challenging task due to the
variety of parameters to be satisfied. These parameters are related with physical
abilities of the plant or the system under control. It may not be possible to achieve a
very small rise time if the natural frequency of the system is high, and higher
damping in the system increases settling time and overshoot, etc. It is difficult to
fully control and change the natural characteristics of the plant through transfer
functions, pole–zero cancellations, or lead-lag compensation. However, control
objectives can be achieved significantly with appropriate control designs. There
are different design methods to appropriately shape the dynamic response of the
plant by feedback loop. A feedback function is designed to operate in a loop as
given in different configurations in Figs. 8.1, 8.2, 8.3, 8.4, 8.5, and 8.6. Simplest of
all approaches is a unity feedback with PID controller given in Fig. 8.3. The
8.5 Controlling the Response 153
simplest and most commonly used controller design is known as the PID control
design or the Proportional–Integral–Derivative control design. There can be other
methodologies as well to design a controller in the feedback loop.
The parameters in the three terms are Kp, KI and KD as the design values in the PID
controller. Kp is a simple proportionality between the output of PID and input to the
PID e(t) simply representing increasing or decreasing gain based upon magnitude of
error. When error is integrated in order to observe settling of the response, then
integral parameter KI is tuned for appropriate gain. The sudden changes and rate of
change of error is observed and tuned through derivative KD control based upon
error of the plant output. The PID controller gain is also given as:
Fig. 8.8 Feedback control loop with unity feedback, output reference, and plant configuration in
state space as well as equivalent transfer function
154 8 Introduction to Control Systems
U ðsÞ KI
K ðsÞ ¼ ¼ Kp þ þ KD s ð8:20Þ
Eð s Þ s
This is the second-order transfer function required to solve for the numerator or
denominator from this equation. This will be achieved by taking inverse Laplace
transform term-by-term and then finding an actual transfer function. There is often a
requirement to use the controller as either PI or PD due to simplification and lesser
design specifications. The integral parameter KI controls the steady-state error
which may occur due to the proportional term Kp. Damping overcomes the over-
shoot, and the speed of the response regulates by the derivative parameter KD.
8.5.2 Compensation
s þ zc
K ðsÞ ¼ K c ð8:21Þ
s þ pc
The subscript “c” is for compensator and if zc < pc then this is known as lead
compensator. Lead compensator is equivalent to the derivative controller term KD,
which is used to decrease the overshoot of the response and increase bandwidth or
damping in the system. If the pole of the compensator is lagging the zero of
compensator zc > pc then the system is known as lag compensator, which is
equivalent to the integral control term KI whose primary objective is to decrease
the steady-state error. A second-order compensator is actually equivalent to the PID
control design, which may have both lead and lag properties for full design
specifications. A compensator is the combination of both regulator and observer
design of the system in state space representation given as
_
* * *
x c ðtÞ ¼ Ac x c ðtÞ þ Bc m ðtÞ
K ðsÞ ¼ ð8:22Þ
* * *
u ðtÞ ¼ Cc x c ðtÞ þ Dc m ðtÞ
A state space compensator can directly take measurement output and generate
regulating commands for the actuators. PID and lead-lag compensator adopts
frequency-based techniques for selection of the PID gain matrices or zero/pole/
gain of lead-lag compensator such Bode Plot, Root locus method, and Nicholas
chart.
8.6 Shaping the Dynamic response 155
Regulator and estimator are designed in a state space formulation in order to change
the closed loop poles of the system. A regulator is achieved by designing a matrix
K taking state input to regulate the feedback command
* *
u ðtÞ ¼ K x ðtÞ ð8:23Þ
The matrix K is of order p n for p inputs and n states and it is selected in such a
way that closed loop poles of the system are placed at the desired locations. The closed
loop poles are active when the feedback command is regulating the system and the
dynamic profile of the system is determined by the location of the poles. In Eq. (8.22),
Cc ¼ K and Dc ¼ 0 for simple regulator design of the closed loop system, where
Bc ¼ L is an observer gain matrix as shown in Fig. 8.5. The observer gain matrix is
chosen such that closed loop poles for observer are stable and faster than regulator
poles. The matrix Ac is a feedback matrix from both regulator and observer poles.
The dynamics response can be shaped either for selection of PID gains, or lead-lag
compensator using MATLAB optimization tools. The performance specification
met using different methodologies of designs. Problem 8.6 describe the simulation
technique in detail for pole placement and shaping the dynamic response
accordingly.
Example 8.2 PID Tuning
Simulate the open loop system and tune the PID gains in order to obtain the desired
characteristics as given in Fig. 8.8 using step input.
* 1 4 * 1 *
x_ ðtÞ ¼ x ðtÞ þ u ðtÞ
4 1 3 ð8:24Þ
*
yðtÞ ¼ ½ 2 2 x ðtÞ
tr < 800 ms ts < 4 s for 10 % of yss Over Shoot < 9 %
Get the information by using stepinfo command and plot the open loop and closed
loop response. Is it first peak in the profile which is used to measure overshoot?
Solution The open loop response of this system given in Fig. 8.7 with following
parameters
S¼stepinfo(y.signals.values,y.
time,1.65,‘SettlingTimeThreshold’,0.1)
RiseTime: 0.2535
156 8 Introduction to Control Systems
SettlingTime: 2.9753
SettlingMin: 1.2576
SettlingMax: 2.5016
Overshoot: 51.611
Undershoot: 13.766
Peak: 2.5016
PeakTime: 1.9015
PID block tunes itself for the different values and options and after tuning to step
input response in order to meet the desired objective, the following parameters are
achieved for yss ¼ 1 as shown in Fig. 8.9.
RiseTime: 0.71328
SettlingTime: 3.5662
SettlingMin: 0.7687
SettlingMax: 1.063
Overshoot: 6.3006
Undershoot: 1.8029
Peak: 1.063
PeakTime: 4.0859
0.8
0.6
y(t)
0.4
0.2
-0.2
0 2 4 6 8 10
time (sec)
Fig. 8.9 Step response of a system with tuned PID controller starting at t¼1 s
8.7 Gain Scheduling in Multivariable Control 157
Controlling more than one output using one or more inputs makes a multivariable
control problem. Each output–input loop is treated differently and some loops may
be ignored due to insignificant contributions in the loop. The input to the plant is
dependent upon the combination of different loops. The gains in each loop need to
be adjusted so that the overall system remains bounded and stable.
Example 8.3 Tune a PID controller using signal constraints block of the following
system:
" # " #
* 2 3 *
0 *
x_ ðtÞ ¼ x ðtÞ þ u ðt Þ
1 3 2
" #
* 1 0 *
y ðt Þ ¼ x ðt Þ
1 2
Solution This is a multivariable system due to two outputs, and PID is a SISO
system so both outputs tune separately. The gain multiplies in both loops in order to
effectively control the plant. The signal constraint can also optimize the both gains
g1 & g2 for PID controllers K1(s) & K2(s) respectively for components of same
input. It is important to note that PID blocks K1(s) & K2(s) are tuned independently
to meet performance specifications before tuning of gains g1 & g2 for input. The
inputs u(t) are given as:
The signal constraint block optimizes for performance specification of net input to
the plant. This signal constraint block optimizes outputs as well, and in many cases
it can also optimize more than one variable at a time (Fig. 8.10).
Fig. 8.10 Simulation block with two PID loops and gain scheduling for input
158 8 Introduction to Control Systems
Problems
P8.1 Draw the block diagram of the systems given in Examples 8.1 and 8.2 by
treating each element differently. There will be n-integrators for n states and
each element of matrices will be a gain in the loop.
P8.2 In Fig. 8.2, if we augment the actuators and measurement into the plant state
space, what will happen to the order of plant, which matrices will be changed
and how this will affect control loop? Take the following transfer functions
and augment with the plant of Example 8.1.
1 2
ActðsÞ ¼ Sensor ¼
s þ 1000 s þ 2000
P8.3 Find a PID gains for the system in Problem 8.2 with the same performance
specifications given in Example 8.2.
P8.4 Simulate the open loop observer of the problem in Example 8.1 and note the
error. Introduce the small random noise in feedback path with noise power
0.0001 W and simulate again to observe changes in error.
P8.5 Tune the block of Problem 8.4 for different rise time and overshoot of the
input. What difficulties are encountered in order to tune PID blocks and input
optimization? What are the minimum rise time, peak amplitude, and steady-
state error achieved in your solution? Is the solution unique?
P8.6 Pole Placement using Transfer function: Simulate Fig. 8.6 for the following
systems:
1. Example 5.5
2. System in Fig. 5.28
3. Problem 7.6
4. Example 8.1
5. Example 8.2
Choose K(s) appropriately such that the zeros of K(s) cancel all the unstable
poles and the poles of K(s) are in the left half plane with minimum order of
transfer function.
(a) Simulate with K(s) as a transfer function and state space and observe if
there is any difference.
(b) Identify the maximum bound of poles of K(s) and explain what happens
if bound is violated.
(c) Each pole placement design will give different performance specifica-
tions; re-iterate the design until better results are obtained.
(d) Tune the performance by constraint optimization as well as PID tuning
and compare your results.
References 159
References
1. Friedland, Bernard. 2005. Control System Design—An Introduction to State Space Methods.
Mineola, NY: Dover Publications.
2. Chen, Chi-Tsong. 1999. Linear System Theory and Design, 3rd ed. New York: Oxford Univer-
sity Press.
3. Heij, Christiaan, André Ran, and Freek van Schagen. 2007. Introduction to Mathematical
Systems Theory—Linear Systems, Identification and Control. Basel: Birkhäuser Verlag.
4. Hespanha, Joao. 2009. Linear Systems Theory. Princeton, NJ: Princeton University Press.
Chapter 9
Recent Applications of Bond Graph
Modeling
Bond graph modeling is used in diverse engineering applications due to its universality
between energy systems. Every year the number of research papers and publi-
cations in the domain of modeling and engineering using this technique is on
rise. Lagrangian modeling is a popular method, especially in robotic manipula-
tors, but with the help of bond graphs these models can be extended with
electromechanical systems of actuators and sensors. In the last decade, numerous
papers in conferences and journals have appeared for bond graph modeling and
simulation. It is being taught at undergraduate and graduate level courses in
different universities showing the developing interest in this technique. Bioen-
gineering or biomedical engineering is a key area of interest for research and
development for its applications. Bond graph modeling made its way to this area
by providing an alternative paradigm for modeling and simulation. In this
chapter, we will focus on the recent applications of bond graph in biomedical
engineering followed by some applications in a few other fields.
Hill-type muscle models are used by physiologists and researchers for the biome-
chanical modeling and development of the musculoskeletal control strategies.
Muscle structures exhibit mechanical, electrical, and chemical properties simulta-
neously. Figure 9.1 represents the block diagram of two Hill-type muscle models;
their respective bond graphs are given in Figs. 9.2 and 9.3, generated with 20-Sim
software.
In these models, we represent general compliance with C, which is reciprocal of
stiffness k, in series and parallel elements component Cs and Cp respectively.
The mass of muscle m is represented by I (inertial component), and viscous
damping of contractile element is represented by B. The source of effort of force
acting on the muscle is Se, which is externally applied to the muscle. The state
variables in Hill’s 1st muscle models are generalized displacement (q) at port 3 and
port 5, and momentum ( p) at port 2. Bond graph can be used to write linear as well
as nonlinear state equations. Equation (9.1) represents the state equation of the 1st
Hill model.
a b
CE SE CE
SE
Fout Fout
PE PE
q3 q5
p_ 2 ¼ SE
Cp Cs
p2
q_ 3 ¼ ð9:1Þ
m
p2 1 q5
q_ 5 ¼
m B Cs
Similarly, we can write state equations for alternate Hill muscle model. In this
model, there are two inertial components m1 and m2 at both ends of muscles with a
total of four state variables. These variables are generalized momentums at port
2 and 8 and generalized displacement at port 5 and 7; output is the force at port
8, which is the rate of change of momentum at port 8.
164 9 Recent Applications of Bond Graph Modeling
q5
p_ 2 ¼ SE
Cs
P 2 P2 1 q5
q_ 5 ¼
m1 m1 Bs Cs
ð9:3Þ
P8
q_ 7 ¼
m2
q q
p_ 8 ¼ 5 7
Cs Cp
γd II
Chain SR
9.1 Bond Graph Models of the Physiological Elements 165
Nout
m
ksp
an intrafusal bag model with sensory region bag are given in Eq. (9.4) with three
state variables, generalized displacement at port 3 and 6, and generalized momen-
tum at port 7 (Fig. 9.6).
1 q q
q_ 3 ¼ SE þ Γ 6 3
B Cs Cp
1 q q p ð9:4Þ
q_ 6 ¼ SE þ Γ 6 3 7
B Cs Cp m
q6
p_ 7 ¼
Cs
The combined effect of the bags and chain are shown in Fig. 9.7 for the model
represented in Fig. 9.4. In this model, mechanical input Se is from muscle force for
fascicle length and velocity; the neural or excitation inputs, i.e., dynamic (γ d) and
static (γ s) fusiomotor inputs are represented by the active bonds in this model via
modulated source of efforts.
1 q11 q17
q_ 11 ¼ Se þ γ s ð9:5Þ
B1 Cp Cs
1 q12 q18
q_ 12 ¼ Se þ γ d ð9:6Þ
B2 Cp Cs
166 9 Recent Applications of Bond Graph Modeling
1 q13 q19
q_ 13 ¼ Se þ γ d ð9:7Þ
B3 Cp Cs
1 q q p
q_ 17 ¼ Se þ γ s 11 17 þ 25 ð9:8Þ
B1 Cp Cs m1
1 q q p p
q_ 18 ¼ Se þ γ d 12 18 25 26 ð9:9Þ
B2 Cp Cs m1 m2
1 q q p p
q_ 19 ¼ Se þ γ d 13 19 25 26 ð9:10Þ
B3 Cp Cs m1 m2
Two neural outputs, i.e., primary and secondary afferent, are given as:
9.1 Bond Graph Models of the Physiological Elements 167
1
Fout 1 ðIaÞ ¼ ðq þ q18 þ q19 Þ ð9:13Þ
Cs 17
1
Fout 2 ðII Þ ¼ ðq þ q19 Þ ð9:14Þ
Cs 18
Golgi Tendon Organ (GTO) responds to tensions in the muscle fibers and sends
proprioceptive feedback to the muscle for activation. GTO is usually represented by
static gains, but in bond graph, detailed models of GTO can be implemented. In a
simple approach, GTO can be modeled as only an elastic spring element or
combination of elastic element with viscous damping of negligible mass. In
Fig. 9.8 Se is the source of force from muscle with viscous and damping elements.
This is the simplest one-degree model of GTO. The output equation for this
model is given as
1
x_ gto ¼ Se kgto xgto ð9:15Þ
Bgto
I
l
θ k
x
168 9 Recent Applications of Bond Graph Modeling
of support, with mass m, inertia I, and a distance k from the center of mass to a joint
was discussed in [1]. The inertial model has joint angle θ and angular velocity θ_ as
state variables for a 2nd order nonlinear system which was linearized at stable
equilibrium (upright) position, and the dynamic state space representation of the
inertial subsystem is given as
€θ ¼ ðτ þ mgk sin
θÞ
ð9:16Þ
I þ mk2
The input torque to the inertial model comprises muscle force transformed into joint
torque. The outputs of this model are angle and velocity that are transformed into
fascicle length and fascicle velocity as inputs to the muscle spindle structure. The
bond graph of an inverted-pendulum model has been discussed in detail in several
textbooks.
The complete musculoskeletal systems with inertia, muscle, muscle spindle and
GTO can be modeled using bond graphs. In this model we also need two-port
elements called transformers for conversion from translational to rotational motion
or vice versa. Inertial equations of joint mechanics are usually represented in
angles, i.e., flexion of limbs (segments) at joints. Muscle equations are mechanical
translational equations, so a transformer is needed at every joint to convert muscle
forces into joints torques. Figure 9.10 shows the block diagram representation of a
system with comprising components.
In Fig. 9.10, muscle output force is transformed into torque for inertial block and
added up to external torque for limb movement. Transformer ratio (r) is the length
of a limb. Similarly angle and angular velocity outputs from inertial block are sent
to muscle spindle by transforming these to fascicle length and fascicle velocity. The
output angle and velocity from inertial block are represented by the vector bond
graphs. The vector bond graphs are used when multiple variables are connecting
one port to another identically. These are also used for representing assembly of
subsystems. The negative ratio (r ) is used here because if force to moment is
positive then the resulting angular motion results in muscle shortening. The muscle
block takes proprioceptive feedback and produces appropriate forces for stable
Ia & I I
CNS
9.2 The Complete Musculoskeletal System Model 169
movement. The muscle block has two components, one connected to central
nervous systems (CNS) and the other producing controlled forces to inertial
block. GTO also responds to the muscle forces and their feedback changes the
neural excitation. The muscle output forces are transformed into joint torques
adding up with external or disturbance torques and then cause the movement in
inertial block. The musculoskeletal feedback system can be combined into a single
state space representation for analysis and robust design of a controller structure.
This dynamic system takes effect of all physiological variables interacting with
each other in the closed loop response. The state variables of closed loop system are
xM, xS, xG, and xI which represent the dynamics of muscle, muscle spindle, GTO
and inertia respectively. The combined musculoskeletal model for the closed loop
feedback system is given in Eq. (9.17).
2 3 2 3 2 3 2 3
x_ M AM AMS AMG 0 xM BM
6 x_ S 7 6 0 AS 0 ASI 7 6 xS 7 6 0 7
4 x_ 5 ¼ 4 AGM 0
AG AGI 5 4
þ
xG 5 4 0 5
U ð9:17Þ
G
x_ I AIM 0 0 AI xI 0
The order of the system is dependent upon the order of each subsystem. The muscle
structure is 4th order, inertia is 2nd order systems, GTO is 1st order system, and
muscle spindle is 6th order when using two intrafusal bags and a chain system
represented in Fig. 9.4. The musculoskeletal system shown in Fig. 9.10 is 13th order
system. However, if a simplified 3rd order muscle spindle structure is used then the
overall system is reduced to 10th order. Only the muscle structure is directly
affected by the input through CNS, and so in B matrix, there are zeros for other
subsystems. The interaction matrices define the relationship between subsystems:
AMS, AMG, represent interaction between muscle structure with muscle spindle and
GTO respectively. Similarly AIM, ASI, AGI represent the interaction of inertial
subsystem with muscle, muscle spindle, and GTO. Transformer actions from
muscle to inertial block and from inertial block to muscle spindle block are
represented by T. There are three exogenous or disturbance inputs in the system,
external (disturbances) torque added to the control torque, neural static (γ s), and
dynamic (γ d) fusiomotor inputs to the muscle spindle. Figure 9.11 presents the
complete musculoskeletal system with Hill’s 1st and alternate muscle models
respectively. In this case there is a differential causality at series compliance Cs
at a zero junction, which adds one redundancy in the system, i.e., one dependent
equation.
Figure 9.12 shows the complete musculoskeletal bond graph with Hill’s alter-
nate muscle model and this doesn’t have any differential causality. The nonlinear
model is given in Eqs. (9.18), (9.19), (9.20), (9.21), (9.22), (9.23), (9.24), (9.25),
(9.26), and (9.27) for the 10th order system. In this system there is only one
fusiomotor input with simple 3rd order muscle spindle model of, 2nd order
nonlinear inertial model, 1st order GTO model and 4th order muscle model. This
is an open loop model, which requires a feedback controller to regulate the sensory
commands.
170 9 Recent Applications of Bond Graph Modeling
Fig. 9.11 Bond graph model of musculoskeletal system with 1st Hill’s muscle model
Fig. 9.12 Bond graph model of musculoskeletal system with alternate Hill’s muscle model
q5
p_ 2 ¼ Se ð9:18Þ
Cs
p2 1 1 q 1 q10
q_ 5 ¼ 5 ð9:19Þ
m B Bgto Cs Bgto Cgto
9.2 The Complete Musculoskeletal System Model 171
p8
q_ 7 ¼ ð9:20Þ
m1
q5 q7 1
p_ 8 ¼ þ τext ð9:21Þ
Cs Cp l
1 q5 q10
q_ 10 ¼ ð9:22Þ
Bgto Cs Cgto
1 p8
x_ 17 ¼ m g k sin ðx18 Þ þ l ð9:23Þ
I þ m k2 m1
x_ 18 ¼ x17 ð9:24Þ
1 1 1 q25 q29
q_ 26 ¼ Γ þ x17 þ x18 ð9:25Þ
Bsp r1 r2 Csp Csr
1 1 1 q25 q29 p
q_ 29 ¼ Γ þ x17 þ x18 30 ð9:26Þ
Bsp r1 r2 Csp Csr m2
q29
p_ 30 ¼ ð9:27Þ
Csr
Fig. 9.13 Bond graph model of extended musculoskeletal structure with muscle spindle
actuator for the prosthesis. The effect of field and armature winding are modeled to
represent the energy dissipation as the winding losses and inductive behavior of
winding coil. The electrical energy is converted into mechanical energy as a
function of field current. Planetary gear assembly is used to convert the rotational
mechanical energy into translational energy, which is the function of 21:13 gear
ratio. As the prosthetic mechanism is actuated through DC motor, a string, attached
to the shaft of DC motor, is used to carry energy up to each phalange of the
prosthetic hand. This string works similarly to the tendons in human hands. Because
each finger has different parameters (length, weight etc), the length of each string is
different with respect to the finger to which it is attached. These strings are fastened
at the shaft just like the bracelet-like tissue that holds the tendons at the carpals in a
real hand. The block diagram represents only one phalange, for compactness of
block diagram, by three cascaded circles representing proximal, middle and distal
extremities. There are five phalanges in the actual model of a hand. The bond graph
model is a precise mathematical description of the physical system in the sense that
it leads to the state space description of the system.
The bond graph model of an anthropomorphic hand with five phalanges is
presented in Fig. 9.15. The bond graph of the anthropomorphic hand has thirteen
energy storage elements, and all have integral causalities.
The anatomy of the hand is efficiently organized to carry out a variety of complex
tasks. These tasks require a combination of intricate movements and finely con-
trolled force production [8–9]. The human finger consists of three rigid bodies or
bones called phalanges. The phalange closest to the palm is called the proximal
interphalange (PIP) followed by the Middle intermediate phalange (MCP), and the
174 9 Recent Applications of Bond Graph Modeling
Fig. 9.15 Bond graph model of an anthropomorphic hand of five phalanges implemented with DC
motor and pulleys arrangement. The notations used in the bond graph are given as: Ia Inductance of
armature winding, Ra armature winding losses, r Gyrator modulus representing the field current,
Imm mass of DC motor, Rf frictional losses of motor bearings, k stiffness of the rotating shaft,
m gear ratio of planetary gear box, N {1, 4, 7, 10, 13} radius of the pulleys between metacarpals
and phalanges of five fingers, N {2, 5, 8, 11, 14} radius of the pulleys between proximal and
phalanges of five fingers, N {3, 6, 9, 12, 15} radius of the pulleys between middle and distal of five
fingers, M f1 5g mass of five phalanges, k f1 5g stiffness of string attached to five phalanges,
R, R f1 4g cumulative damping of pulleys in each phalange
last one, the distal interphalange (DIP). Joints have unique characteristics to
provide compliant motion in different directions. The function of tendons is to
transmit forces from muscles to bones in order to actuate joints. The tendon
lubrication system reduces the friction during tendon excursion. In many cases
tendons cross more than one joint and impart rotation to these joints. The extensor
mechanism of the fingers contains complex tendinous networks. The nonlinear-
elastic approach considers the ligament to respond like a nonlinear spring. By
adding a dashpot in series or in parallel to the spring, time-dependent behavior is
achieved. The finger kinematic chain consists of three rigid bodies, connected to
each other by joints. First, a single rigid body with two connection points for the
joints is modeled. Bond graph notation is used to create the model of the rigid body.
The actual finger is made such that the phalanges can only rotate around one axis,
9.4 Movement Coordination Problem for Human Fingers 175
and the center of gravity can only move in a plane perpendicular to this rotation
axis. Muscles provide rotational and translational motion. Joints between the bones
are modeled as revolute joints with viscoelastic behavior due to the presence of soft
cartilage. Viscoelastic behavior is modeled by providing C and R elements between
the 0 junctions of all the bones. The nature of elements C and R depend on the
properties of the biological material. Similarly, rotational motion is allowed about
one axis only and limited motion is permitted about the other two axes. Consider the
distal interphalangeal joint of the little finger: when a force is applied at distal
phalanx, a rotational motion is experienced by the phalanx.
The joint tries to oppose rotational motion (bending) beyond a certain limit in
order to avoid any injury. That opposition is caused by joint stiffness, similar to
stretched spring and friction. Accordingly, a joint can be modeled via a capac-
itance (C) and resistance (R). Although force has been applied at the distal
phalanx of the little finger, it is being transferred to adjacent phalanx through
DIP, so a joint performing force transformation also acts as a transformer. The
translational momentum of the phalanx can be considered to be concentrated at
the center of its mass. Since the motion of the center of gravity of a rigid body, or
the rotation with respect to a fixed axis, is defined with reference to an inertial
frame, kinetic energy stores directly to the 1-junction representing their velocity
or angular velocity. The body’s kinetic energy is a function of its velocity. In
six-links Biomechanical model, three joints DIP, PIP, & MCP of both little and
ring finger will be considered. The complete system with bones’ inertia, joints’
stiffness, and friction is modeled with the help of bond graph. Five 2- port
elements have been used for conversion from translational to rotational motion.
Torque is distributed across joints through 2-port elements called transformers
and 3-port elements called 0-junction and 1-junction as shown in Fig. 9.16. There
are twelve energy storage elements, so the order of the system will be 12; leading
to 12 state space equations. Resistances R, are the viscous-elastic damping,
Compliances C are stiffness and inertia I are the masses of the corresponding
segments.
Fig. 9.16 Bond graph model of a movement coordination model of little and ring finger
176 9 Recent Applications of Bond Graph Modeling
Heat transfer mechanisms play an important role in physical systems. There are lots
of physical systems in which their study has been conducted in different ways. A
baby incubator is a special type of biomedical device that provides an ideal
environment for the baby and operates on the mechanism of thermodynamics
because premature babies are unable to keep themselves sufficiently warm due to
their immature regulation systems. They are also very weak and get infected very
quickly. There are four types of thermodynamic processes in baby incubators:
conduction, convection, evaporation and radiation. Each plays the important role
of keeping the infant warm. There are two types of evaporation that occur in infant
incubators. The first evaporation occurs at the surface of the infant due to thermo-
regulation. To overcome this problem a second type of evaporation is produced in
an incubator by introducing air from the outer environment through a mechanical
process. This air is warmed with the aid of a heater; water in the tank which is
placed in the incubator is heated with the help of this warm air to produce vapors
and to warm the inner side of incubator so that the required temperature is provided
through latent heat for the infant. Moreover, there is no method for controlling and
providing a constant evaporation, which is produced by air coming from the outer
environment. We use the pseudo bond graph approach to model the second type of
evaporation process in a baby incubator. We show that an incubator is an open loop
system and equations for an incubator with the help of the first law of thermody-
namics can be obtained to further develop a controller design as shown in [11]. The
major cause of hypothermia is the heat loss in small infants that happens due to the
evaporation from the skin of infants. In order to overcome this problem outside air
is entered into the incubator chamber which is warmed through a heater. With the
help of a water tank in the incubator, vapors produce the latent heat which warms
the infant body and saves it from this heat loss. These vapors then go to the chamber
of the incubator and produce the latent heat which hence provides the required
temperature and humidity to the infant. To maintain a constant temperature and
humidity in the chamber of the incubator, a microcontroller is placed that takes the
information from a temperature and humidity sensor, therefore monitoring as well
as controlling the whole mechanism of air and heat flow in the incubator according
to the set point. According to Pseudo Bond graph methodology shown in Fig. 9.18
we modeled the air as a source of flow, heating unit (heater) as a heating capacity
and processes of coveting water into vapors as a gyrator. Additionally we introduce
a capacitance as an addition of water vapors into warmed air, and at the output the
heat flow in the chamber is modeled as a source of flow.
This output of this pseudo bond graph model is a value from heat flow sensor
which is an input to the controller in order to regulate the air flow (Sf).
178 9 Recent Applications of Bond Graph Modeling
Due to the recent advancements in the “Robotics,” new and customized robotic
parts imitating human body movements are developed everyday. These kinds of
developments also help where the presence of human beings is either injurious or
deadly. The field of physiotherapy has reached new milestones owing to the
introduction of robotic arms and limbs, which replace these amputated parts of a
human body [6]. Their need is not only widespread but also desired to make life
easier and more comfortable. Their use can be varied to fit our own purposes and
requirements. This example discusses modeling of a miniaturized robotic gripper,
which moves in a vertical and then horizontal direction to grab any object in
consideration [12]. The arm has two aluminum bars at its extreme ends which
perform the griping action when instigated. The complete structure on which this
arm is placed along with the arm is shown in Fig. 9.19.
The assembly is such that the gripper arm is initially at rest at an initial vertical
position, with the two gripping rods of the arm stretched apart. The first motion is
the movement of the whole gripper bar from vertical initial point to a vertical final
point. After that has been achieved, the bars that perform the gripping action move
towards each other. In a modeling approach the whole system in considered as a
single plant. A Block diagram of the system is shown in Fig. 9.20. The purpose of
9.8 Bilateral Master–Slave Telemanipulation 179
the first motor is to convert the electrical energy through voltage source into
rotational mechanical energy. The requirement is to move the bar vertically or
translational upwards, so a transformer represents this motion by converting rota-
tional mechanical energy into translational. After the transformation of the energy,
this bar reaches the vertical final point and closes a limit switch. This limit switch
performs the action of a gyrator because the mechanical energy is converted into
electrical energy by relating the velocity (flow) to voltage (effort). After reaching
final height, the two gripping rods receive input energy to start moving towards
each other, a motion which is again represented by a transformer followed by a
shock absorber to dampen any jitters. Motors are represented as gyrators and output
is monitored before shock absorbers. The bond graph of the system is given in
Fig. 9.21.
Human Hybrid
Master Slave Environment
Operator Controller
Teleoperation
There is hybrid feedback controller between master and remote slave, which
interacts with the environments. The operator commands are usually in the form of
force or position input. Similarly, the environment can also exert the force on the
slave or change its position. Either combination can also be applied in the entire
model and also in the master–slave dynamics. The operator force command is
transferred to the slave force command through master and the position commands
work similarly. However, force and rate of change of position (velocity) are the
variables of interest in our systems and these are described as power variables in the
bond graph methodologies. The hybrid model between master and slave is
presented in Fig. 9.23. The two-port model presented in Fig. 9.23, can be modeled
in various types of hybrid systems depending upon the input and output consider-
ations. From master force in to slave force out system impendence matrix is
calculated; similarly velocity in from master to velocity out from slave admittance
matrix is to be calculated. Hybrid matrix gives the relationship from master force to
9.8 Bilateral Master–Slave Telemanipulation 181
x´1 x´2
Zh Ze
Fh + fm Master fs + Fe
Slave
x´1 x´2
Bh Be
Mh Mm Bm x´m Ms Bs Me
βx´m
+
Fs Bl
Fh Kh αFs + Ke Fe
Ks
-
Master - Slave
the slave velocity and vice versa is the chain matrix. In this paper teleoperation
model for complete network presented in [13] is modified by incorporating the
human and environmental impedance characteristics as shown in Fig. 9.23. In that
figure master–slave configuration is replaced by the dependent voltage and current
sources. The voltage-dependent voltage source at master is proportional to the slave
force (voltage) and current-dependent current source at slave end is proportional by
the velocity (current) at master [14]. Fh is the force delivered by the human operator
through the human operator impedance to the master.
The human operator impedance consists of mass Mh of the human operator arm
acting on master; Bh and Kh are the viscous and stiffness coefficients of human
muscles exerting force. The Mm and Bm are the mass and viscous coefficients of the
master arm. Similarly for the slave, there is slave mass Ms, viscous coefficient of
slave Bs and viscous coefficient of loss Bl at slave end. Ks represents the stiffness of
the slave and Bl is the due to the losses at slave end and these are sometimes
negligible. Me, Be and Ke are the mass, viscous coefficients and stiffness of the
environment or the object slave is interacting with the system. This object at
environment also exerts some contact or reaction forces Fe to the slave arm as
shown in block diagram in Fig. 9.24. The objects at environment variables have a
large range to choose depending upon the specific applications of the robot. If a
182 9 Recent Applications of Bond Graph Modeling
dynamic robot is designed to pluck a flower or to move high metal, these variables
are different keeping the human operators, master arm, and slave arm variables in
the same range. The velocity of the master arm is shown as x 0 m and this is treated
as the hybrid variable at slave end as x0 s ¼ βx0 m. The two variables of velocity
must be theoretically the same but there is some difference due to the delays. β ¼ β
ðtÞ is a function that will correspond to the delays between master and slave
velocities as x0 s ¼ x0 mðt T Þ, where T is the time delay. When T ¼ 0 or the
delay time is short enough as compared to responses or is not of precise importance,
β ¼ 1 can be assumed. The force at slave end is controlled and acted through
voltage (force) control voltage (force) source. The Fm ¼ Fs ðt T Þ ¼ αFs where
α ¼ αðtÞ. The slave force is available to the master for transmission but not the
contact or reaction force of environmental object, Fs ð1 þ αf ÞFe . Therefore, in
this paper we have used the total impedance of the environment as
Zte ð1 þ αf ÞZ e , where αf is force gain. Figure 9.25 presents a bond graph
model from human operators to the environmental object with impedance charac-
teristics. Due to analogies between mechanical and electrical systems, an electrical
network similar to the networks is adopted for the bond graph simulation. As the
system consists of both electrical and mechanical components, generalized vari-
ables are used in the bond graph model. The power variables are “e” and “f” for
effort and flow. Momentum and Displacement are represented as “p” and “q.” The
bonds at 10 and 11 are active bonds for effort (voltage or force) and flow (velocity
or current). These active bond relations are described in Eq. (9.28).
e9 ¼ α e11
f 9 ¼ f 11 ¼ 0
ð9:28Þ
e12 ¼ e10 ¼ 0
f 12 ¼ β f 10
Problems 183
This is so because current flow in the 9 and 11 bond remains the same and
modulated by an external power source, but effort is changed which is represented
by the voltage controlled voltage source, and hence in active bond with source of
effort between 9 and 11. Similarly for 10 and 12, the effort remains the same and
modulated by an external power source, but flow changes, which is represented by
the source of flow between 10 and 12. Fh represents the force exerted by the human
operator and Fe represents the contact or reaction force by the object at environ-
ment. If the contact force is negligible then SEe ¼ 0 can be assumed, especially
when a relative big slave arm is pulling the small object. This assumption will not
make any changes in the bond graph. The losses at slave arm are represented by Bl.
In case no loss condition is to be required for simulation, bond 14 will be replaced
as bond 13 and bond 15 should be deleted.
Problems
both position and velocity sensors at bond number 22 and obtain state space
equations. Convert this problem into a feedback control system design with
two outputs (i.e., position and velocity at bond number 22).
References
A D
Active bond, 40–42, 58, 59, 61, 165, 182, 183 Damped, 109–117, 119, 121–123, 126,
Actuators, 40, 42, 54, 143, 146, 154, 158, 129–133, 140, 151, 152
161, 173 Decaying, 48, 107, 116, 124, 131–134,
Algebraic loops, 61, 85–87, 89, 100 137–139, 152
Amplification, 42, 100, 110, 118–119, 123, Derivative, 7, 9, 19, 23, 24, 31, 49, 51–52, 83,
126, 134 87–90, 153, 154
Analytical modeling., 4–6 Discrete systems, 6
Anthropomorphic hand, 172–174 Distributed model, 5
Autonomous system model, 5 Dynamical system, 1–17, 35
B E
Bilateral, 179–183 Eigenvalues, 105–108, 110, 112–118,
Bridge circuit, 96, 100, 101 121, 124, 131, 134–136,
138–141
Electromechanical network, 72–74, 78, 79
C Energy variables, 26, 33–37, 42,
Capacitor bank, 100, 101 45, 49, 162
Cartesian, 22–24, 26–28, 32 Error model, 146–147
Causality, 4, 40–42, 48–62, 64, 66, 67, 69, 70, Estimator, 147–150, 155, 183
73–74, 82, 85, 87–91, 95, 96, 98, Evaporation, 177–178, 183
100–102, 169, 173 Excitation, 36, 39, 40, 104, 121, 124, 128–133,
Combination of elements, 82, 96–97 138, 150, 164, 165, 169
Compensation, 145, 152, 154 Exogenous inputs, 144–145
Competence, 3 Experimental validation, 4
Complex poles, 121, 139 External stability, 134–137
Computational systems, 6–7
Constitutive equations, 59–61, 91
Control design, 150–154 F
Controller, 78, 82, 118, 143–146, 148, Fields, 9, 20, 21, 24, 27, 45, 85, 90–97, 101,
152–157, 169, 177, 180 103, 161, 173, 174, 178
Control system, 2, 79, 125, 143–158, 183 Forced, 104, 123–129, 134–139
Critically, 113–114, 116, 118, 120,
125–133, 140
G N
Gain scheduling, 157 Natural frequency, 109, 111, 112, 117, 129,
Golgi tendon organ (GTO), 83, 167–169 131, 139, 140, 151, 152
Gyrator, 52, 54, 55, 58–59, 64, 72, 73, 78, Non-causal system, 6
90, 92, 96, 97, 100, 174, Nonlinear systems, 6, 9, 15, 17, 32, 147, 168
176, 177, 179
O
H Observer, 147–150, 154, 155, 158
Harmonic oscillator, 139, 140 Oscillations, 109–113, 115, 116, 123, 133
Hinge movement, 60 Outputs, 3, 6–9, 11–17, 39–42, 56–59, 61, 68,
Human fingers, 173–175 70, 81, 82, 84, 89, 93, 95, 100,
104–107, 121, 124, 125, 127, 129,
134–136, 143–145, 147–150, 153, 154,
I 157, 163, 164, 166–169, 176,
Impulse, 39, 124–127, 137 177, 179, 180, 183
Inertial, 162, 163, 167–169, 175 Overdamped, 111–114, 116–118, 120,
Infant incubator, 177–178 125–133, 138
Integral causality, 50–52, 64, 66, 69, Over Shoot, 155
72, 73, 173
Internal stability, 121–123, 134–136
Inverted pendulum, 29, 30, 32, 146, 167, 168 P
Performance specifications, 150–152, 155,
157, 158
J Physiological elements, 161–168
0-Junction, 54–58, 63–66, 68, 69, 72–75, 80, PID tuning, 155, 158, 171
88, 92, 96, 169, 175 Piston and cylinder, 60
1-Junction, 54–58, 63–66, 68, 69, Plethysmography, 176, 183
71–73, 75, 77, 78, 80, 85, Pole placement, 155, 158
87, 88, 96, 175 1-Port capacitor, 49–50, 74
Junctions, 47, 54–58, 63–66, 68, 69, 1-Port elements, 47–51, 54, 71, 90, 97
71–75, 77, 78, 80–82, 85, 87, 88, 2-Port elements, 52–54, 59, 61, 63, 64, 175
92, 96, 98, 162, 169, 175 1-Port inertia, 50–51, 75
1-Port resistor, 48–49
Power variables, 35–40, 48–51, 180, 182
L Proportional–integral and derivative (PID)
Lagrangian, 19–33, 63, 161 control, 73, 146, 150, 152–154, 156, 157
Law of conversation of energy, 20, 54 Pseudo bond graph, 37, 177, 178, 183
Linearization, 10–12, 16, 17, 82–84, 147
Lumped system model, 5
R
Redundancy, 4, 82, 169
M Redundant systems, 4, 17, 52, 88
Marginally stable, 121, 141 Reference inputs, 143–146, 149
Master–Slave, 179–183 Resistive, 92, 95–96, 101
Movement coordination, 173–175 Resonance, 129–131
Multi-energy system, 63 Response, 7, 35, 39, 70, 103–141, 143,
Multivariable control, 157, 183 150–156, 169, 182, 183
Muscle spindle, 164–169, 171, 172 Rise time, 151, 152, 158
Muscular structures, 162–164 Robotic gripper, 172, 178–180, 183
Musculoskeletal, 162, 164, 168–172, 183 Rotating pendulum, 32
Index 187
S T
Second order system, 31, 85, 103, Telemanipulation, 179–183
107–115, 117–121, 129, 138, Thermodynamic, 177
150–152, 154 Time varying/time invariant system, 5–7
Sensors, 42, 43, 45, 54, 61, 70, 72, 73, 82, Total, 7, 14, 20, 30, 31, 54, 103, 105, 137–139,
143–147, 150, 158, 161, 164, 163, 182
165, 177, 183 Transfer function, 12–16, 19, 39, 84, 100, 105,
Series motor, 93–94, 100 107, 109, 112–114, 117, 124, 127, 129,
Settling time, 131, 151, 152 135, 136, 140, 143, 145, 146, 149, 150,
20-Sim, 70, 82–84, 91, 92, 162 153, 154, 158
Simple pendulum, 28, 146 Transformer, 41, 52–54, 58–60, 64, 68, 69, 73,
Simplicity, 3, 97 78, 87, 90–93, 96, 97, 100, 102, 168,
Sinusoidal, 48, 110, 115, 116, 123, 124, 128, 169, 175, 179
129, 131–138
Sources, 35, 39–43, 45, 47–48, 51, 59–62, 64,
66–69, 71–75, 78, 84, 87, 89–91, 94, 96, U
97, 100, 104, 144, 149, 162, 164, 165, Undamped systems, 109–110, 114, 116, 117,
167, 177, 179, 181–183 120, 128–129, 131, 138, 139
Spherical, 22–24, 32 Under damped, 111, 114–116, 118, 120,
Spring, 50, 70–73, 76, 77, 79, 88, 90, 92, 93, 125–133, 140
97, 128, 139, 161, 164, 167, 174, 175 Undershoot, 112, 151, 152, 156
Stable systems, 119, 121, 156 Unstable systems, 119, 121, 126, 129, 156
State space method, 8–9
State variables, 2, 8–10, 14, 16, 33, 52, 58,
74–76, 81, 88, 89, 92, 146, 163, 165, V
168, 169 Vector bond graphs, 47, 97–102, 168
State vector, 2, 5, 8, 9, 28, 29, 32, 74, 78, 83,
98, 99, 104, 145–147
Steady state value, 116, 128, 150–152 W
Step, 39, 43, 51, 64, 82, 124, 127–128, 155, Word bond graph, 39–43, 45, 60, 64
156, 179, 183 Work energy principle, 19
Systems engineering, 1, 2