Integrated Simulation
Integrated Simulation
Integrated Simulation
Thèse de doctorat
présentée en vue
de l'obtention du grade de
Docteur en Sciences Appliquées
par
Olivier BRÜLS
Ingénieur Civil Electro-Mécanicien (Mécatronique - Productique)
Aspirant F.N.R.S.
Février 2005
Abstract
A mechatronic system is an assembly of technological components, such as a mech-
anism, sensors, actuators, and a control unit. Recently, a number of researchers and
industrial manufacturers have highlighted the potential advantages of lightweight paral-
lel mechanisms with respect to the accuracy, dynamic performances, construction cost,
and transportability issues. The design of a mechatronic system with such a mechanism
requires a multidisciplinary approach, where the mechanical deformations have to be
considered. This thesis proposes two original contributions in this framework.
(i) First, a modular and systematic method is developed for the integrated sim-
ulation of mechatronic systems, which accounts for the strongly coupled dynamics of
the mechanical and non-mechanical components. The equations of motion are formu-
lated using the nonlinear Finite Element approach for the mechanism, and the block
diagram language for the control system. The time integration algorithm relies on
the generalized-α method, known in structural dynamics. Hence, well-dened concepts
from mechanics and from system dynamics are combined in a unied formulation, with
guaranteed convergence and stability properties. Several applications are treated in the
elds of robotics and vehicle dynamics.
(ii) Usual methods in exible multibody dynamics lead to complex nonlinear mod-
els, not really suitable for control design. Therefore, a systematic nonlinear model reduc-
tion technique is presented, which transforms an initial high-order Finite Element model
into a low-order and explicit model. The order reduction is obtained using the original
concept of Global Modal Parameterization: the motion of the assembled mechanism is
described in terms of rigid and exible modes, which have a global physical interpreta-
tion in the conguration space. The reduction procedure involves the component-mode
technique and an approximation strategy in the conguration space. Two examples are
presented: a four-bar mechanism, and a parallel kinematic machine-tool.
Finally, both simulation and modeling tools are exploited for the dynamic analysis
and the control design of an experimental lightweight manipulator with hydraulic ac-
tuators. A Finite Element model is rst constructed and validated with experimental
data. A reduced model is derived, and an active vibration controller is designed on
this basis. The simulation of the closed-loop mechatronic system predicts remarkable
performances. The model-based controller is also implemented on the test-bed, and the
experimental results agree with the simulation results. The performances and the other
advantages of the control strategy demonstrate the relevance of our developments in
mechatronics.
Acknowledgements
This dissertation is the result of several years of research at the University of
Liège, and I feel deeply indebted to a number of people who directly or indirectly
inspired its realization.
First, I would like to thank my advisor Professor Jean-Claude Golinval for his
precious help and support during those years. My gratitude also goes to Professor
Pierre Duysinx for his technical guidance, and for his enthusiastic encourage-
ments. Every discussion with Professor Michel Géradin had a constructive inu-
ence on the orientation of this research, I am pleased to acknowledge him. Those
three Professors highly contributed to the improvement of the draft manuscript, I
am grateful to them for their advice during its elaboration.
I wish to express my thanks to Professor Wayne J. Book for his cheerful
welcome at the Georgia Institute of Technology in summer 2004, and for making
my stay as protable as it could be from all points of view. During the stay,
I especially appreciated the collaboration with Ryan Krauss, who initiated me
with the Ralf test-bed. Our numerous discussions were highly valuable for the
elaboration of chapter 6. I also acknowledge the technical support of J.D. Huggins
during this experimental work.
The IAP/AMS framework oered great opportunities to cooperate with bel-
gian specialists in mechatronics. The research on the semi-active suspension has
been initiated by the teams of KULeuven-PMA and UCL-PRM, which are grate-
fully acknowledged. In particular, I am indebted to Professor Paul Fisette (PRM)
who took the time to provide the data and the complementary results associated
with this problem. I am grateful to Denis Joassin, whose master thesis contributed
to the modeling of the suspension.
I wish to thank Professor Rodolphe Sepulchre for interesting discussions about
the application of control theories to exible multibody systems.
Professor Philippe Wenger (IRCCyN, France) is also acknowledged for mak-
ing the Orthoglide data available.
I would like to thank the Oofelie community, in particular the team of Open
Engineering, for its invaluable support in the software implementation of the
methods presented in this dissertation. Special thanks go to Professor Alberto
Cardona and Elisabet Lens (CIMEC-INTEC, Argentina), whose active collabo-
vi
1 Introduction 1
1.1 Integrated design in mechatronics . . . . . . . . . . . . . . . . . . 2
1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.2 Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.2.1 The constraint elimination method . . . . . . . . . . . . . 44
3.2.2 The augmented Lagrangian method . . . . . . . . . . . . . 44
3.3 Time-integration . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.3.1 Newmark algorithm . . . . . . . . . . . . . . . . . . . . . . 46
3.3.2 Generalized-α method . . . . . . . . . . . . . . . . . . . . 47
3.3.3 Linear stability analysis . . . . . . . . . . . . . . . . . . . 50
3.3.4 Convergence analysis . . . . . . . . . . . . . . . . . . . . . 54
3.4 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7 Conclusion 229
Bibliography 235
I
Introduction
Reference -
Controller - Ampliers - Actuators - Transmission - Mechanism -
trajectory
6
Sensors
the base, reducing the moving mass and simplifying the transmission mechanisms.
A diculty comes from the complexity in the geometrical workspace and in the
nonlinear kinematic relation between the motion of the actuated joints and the
eector location.
From this early analysis, the advantages of exible and/or parallel designs
appear especially attractive for high-speed applications (e.g. high-speed robots
and machine-tools), and problems where the mass is a critical issue (e.g. large
manipulators, space robots and foldable structures). In order to design a con-
troller for those complex dynamic systems, an appropriate mechanical model is
highly desirable at this level.
For a wide range of practical problems such as car suspensions and tradi-
tional machine-tools, the exibility of the mechanical system is not dominant,
and the complex kinematic behavior can be approximately represented by lin-
ear equations. After a design based on simplied linear models, the engineer
may worry about the validity of the predicted performances. Then a multidisci-
plinary simulation tool may be extremely valuable at a pre-prototyping stage, in
order to achieve inexpensive checks and performance predictions before physical
implementation and testing.
Even though the rigid body and serial topology constraints are discarded,
the sequential design is still a pragmatic approach, which allows a relative inde-
pendence in the design of the mechanism and of the controller. This sequential
approach suers from inherent limitations which can be overcome by the devel-
opment of a concurrent optimization of the whole system. At this level, a mul-
tidisciplinary simulation tool is valuable for the estimation of the optimization
criterion and its sensitivities.
Those modeling, simulation, and optimization concepts are relevant at spe-
cic stages of the design procedure. In this work, we assume that a rst design
of the system is available, which includes the description of the mechanism and
of the actuators. Hence, the objective is to dene the control strategy, and to
optimize the remaining design variables. A procedure to solve this problem may
involve the following steps:
In order to obtain the most appropriate information at each step of the design
procedure, various models should be considered. Two important distinctions
among modeling concepts are introduced in the next two paragraphs.
Mathematically, the dynamic model of a mechanism is naturally expressed
by the equations of motion, which dene an instantaneous relationship between
the external actions (actuator forces, disturbances), the mechanism generalized
coordinates and their time-derivatives (velocities and accelerations). The equa-
tions of motion, supplemented with a solver algorithm, constitute a simulation
model of the system, relating time-domain responses (displacements, strains, etc)
with time-domain actions (applied forces, commands, etc). In system theory, the
same distinction can be done between the state equations, and their formulation
with a simulation algorithm.
For complex multibody systems, a modeling software is helpful to formulate
automatically the equations of motion from a high-level description. Among the
computer modeling methods, symbolic methods allow to build the equations of
motion in symbolic format, whereas numerical methods produce the equations
of motion as complex numerical procedures. The symbolic format has the ad-
vantages of portability and eciency, and it provides interesting insights in the
analytical structure of the equations. However, numerical methods are able to
deal with a more general class of problems, and they are especially suitable to
model the dynamics of a exible mechanism with complex topology in a system-
atic way.
After this clarication, let us further characterize the modeling requirements
in the design procedure, which are directly associated with the objectives of this
research.
In step 2, the control design usually exploits the mathematical structure of
the mechanical equations of motion formulated in step 1. If a linear model is
able to capture the essential dynamics, linear control theory can be eciently ap-
plied. Using linearizations at dierent operating points, it is sometimes possible
1.2. OUTLINE 5
to generalize this paradigm for nonlinear systems. For the critical applications
considered here, a nonlinear control technique is necessary, which advantageously
exploits the structure of a nonlinear model. Since the order of the controller is
often proportional to the order of the model, a low-order model in explicit and
analytical format is desirable at this stage. However, existing modeling meth-
ods for complex exible mechanisms are not able to construct such a nonlinear
model; they lead to high-order and computationally demanding models, which
are denitely not appropriate. An important contribution of this thesis is a sys-
tematic reduced-order modeling method, leading to nonlinear mechanical models,
compatible with the requirements of a control design procedure.
Steps 3 and 4 involve a multidisciplinary simulation model, for which less
stringent simplicity requirements are imposed: the computational load and the
complexity of the model should be balanced with the generality, reliability, modu-
larity and systematic implementation issues. In this thesis, we propose an original
integrated simulation tool for mechatronic systems based on a unifying Finite El-
ement formulation. This reliable method allows a modular model denition for
both the exible multibody dynamics and the dynamics of the control system.
In the presentation, a specic eort is delivered to dene natural connections
between the elds of modeling and control. Standard formulations from both
areas are combined for the developments of our original mechatronic concepts.
Hence, we hope that our point of view can be understood, exploited and developed
by specialists from both elds.
1.2 Outline
For the sake of consistency, the dissertation follows a progression from analy-
sis to design concepts, in contradiction with the sequence of the design procedure.
Thus, the topics of integrated simulation, model reduction and control design are
successively addressed in chapters 4, 5 and 6.
Prior to the presentation of our personal contribution, an extended state of
the art is developed in chapter 2. Preliminary modeling and simulation concepts
in multibody system dynamics are reviewed before specialized discussions on the
simulation of mechatronic systems, the reduced-order modeling and the control
of exible mechanisms. The Finite Element approach proposed by Géradin and
6 CHAPTER 1. INTRODUCTION
Cardona for exible multibody systems [GC01, Car89, CG89, CG88] is at the
core of our developments, and an overview is presented in chapter 3.
The simulation of mechatronic systems is addressed in chapter 4. The block
diagram language is selected for the description of the control system dynamics,
and the generalized-α integration scheme is extended to deal with the dynamics
of the state variables. After a formal presentation of the simulation method, theo-
retical convergence and stability results are established. Three examples illustrate
the generality and the eciency of the formulation: a four-bar mechanism, an
industrial robot, and a car equipped with a semi-active suspension.
Chapter 5 is devoted to the reduced-order modeling of exible mechanisms.
The method is an extension of the linear component-mode synthesis, which ac-
counts for the nonlinear kinematics of the system. Therefore, the original concept
of Global Modal Parameterization is dened, and the reduced equations of mo-
tion are formulated. The variations of the model in the conguration space are
approximated by piecewise polynomial functions. Critical implementation issues
are discussed, in connection with the control design requirements. The method
is illustrated with two examples: a rigid parallel machine-tool, and a exible
four-bar mechanism.
In chapter 6, the theoretical tools presented and validated in chapters 4 and 5
are exploited for the modeling and control design of a large and exible manipu-
lator, which is represented in Figure 1.2. A reduced-order model is validated with
experimental results, and the new insights in the dynamics of the manipulator
lead to the development of an original control law. Numerical simulations of the
closed-loop mechatronic system predict that a compromise can be obtained be-
tween the motion bandwidth and the stabilization of the mechanical vibrations,
which is conrmed by the experimental results. The various advantages of the
control strategy demonstrate the relevance of a global approach in mechatronics.
Our contribution is thus described from the general to the particular: simula-
tion of mechatronic systems, modeling of exible mechanisms, and application to
an experimental manipulator. Finally, chapter 7 draws conclusions and discusses
possible directions for future research.
1.2. OUTLINE 7
Lumped approaches
Wanner [HW91].
14 CHAPTER 2. STATE OF THE ART
mechanisms.
3. how does the simulation algorithm deal with the coupled problem?
not be considered here. At the kernel level, the concepts are built using the syn-
tax of the programming language, and are exploited to establish the element and
analysis libraries. The element library aims at describing every subsystem and its
contribution to the global equations of motion, whereas the analysis module de-
nes the treatments to solve those equations. At the top, the user interface gives
access to the element and analysis libraries in order to construct the simulation
model.
In the following, several modular simulation approaches are successively con-
sidered: the Finite Element method, physical modeling techniques (the bond
graph and the linear graph), variational principles, mathematical approaches
(the block diagram). We shall discuss their ability to represent the dynamics
of a exible mechanism within the mechatronic system. A general comparison is
summarized in Tables 2.1 and 2.2: the answers to the three fundamental questions
(elementary description, coupling strategy and time-integration) are decomposed
according to the software architecture given in Figure 2.1. We encourage the
reader to refer to them during the presentation of the dierent methods.
Physical modeling methods are especially well-suited when the port variables
are scalar quantities. Among the numerous references in the literature, we may
cite [OE97, KMR00] for bond graphs and [KTKH67] for linear graphs. Recently,
a signicant eort has been spent for the standardization of multidisciplinary sim-
ulation techniques, leading to the development of the unifying modeling language
Modelica [EMO98, Til01], which is based on bond graph theory.
Several authors have extended the concept of bond graph and linear graph
for multibody systems, see [Kar97, Fav97] and [McP96, McP98, McP03, SM03a],
respectively. According to the comparison study by Sass [Sas04], the linear graph
is more appropriate in that case, and it is noticeable that the topological analysis
automatically solves the kinematic problem. The extension of the linear graph
method for exible multibody systems was also presented in [SM00, SMH01],
relying on an a priori spatial discretization with assumed-modes. McPhee et
al. [SM03a, SM03b] demonstrated the relevance of the linear graph method in
electromechanics.
engineers prefer the more general block diagram formalism, where the equations
of the subsystems are directly manipulated. The block diagram is a mathemat-
ical language for the modular description of dynamic systems in terms of input,
output, and internal variables; it should not be associated with a particular sim-
ulation method. Depending on the coupling strategy, three dierent approaches
are considered below.
The weakly coupled strategy, available in several commercial software (Simulink,
Vissim, etc.) is often regarded as the standard approach for the simulation of
a block diagram model. Each individual subsystem is usually represented by
ODEs, but algebraic equations come from input/output interconnections. The
time-integration is conducted separately for each block in an asynchronous way,
and the coupling results from the exchange of input/output numerical values. The
sequence followed for the successive treatments of the blocks is obtained from a
causality analysis, itself relying on the input/output causality property of every
block. However, this causality analysis breaks down in case of algebraic loops3 , so
that reliable algebraic loop detection and solver algorithms are essential for con-
sistent simulation results. As represented in Table 2.2, the ODE solver directly
deals with the individual blocks, so that the separation between the element and
analysis library is virtual. For this reason, the weak coupling strategy is not
suitable for globally implicit integration schemes, and it is not recommended for
sti problems.
Partitioned simulation methods can also be formulated mathematically in
terms of block diagrams. Each subsystem can hide a highly complex dynamics
(e.g. the mechanical or the electrical part of a mechatronic system), and a parti-
tioned treatment intends to improve the computational eciency and the mod-
ularity of the simulation. Among those methods, let us mention the multi-rate
integration strategies for ODEs suggested by Gear and Wells [GW84], their gener-
alization by Schiehlen and co-workers to DAEs with algebraic loops [RS98, KS00],
and the methods based on Gauss-Seidel iteration for systems coupled by alge-
braic constraints developed by Arnold et al. [Arn01, HAV03]. Co-simulation
techniques are also partitioned methods, where the subsystems are implemented
in dierent specialized software [VBSV99, VGVDS+ 99]. In this case, a master
3 An algebraic loop is a topological loop where all the blocks are direct-feedthrough, which
means that their outputs are instantaneously aected by their inputs.
2.3. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS 21
linker module manages the task scheduling, and the calling procedures of the
other software.
Finally, the block diagram language leads to a strongly coupled simulation if
the coupling is imposed prior to the integration scheme: an assembly procedure
is necessary to build the coupled equations, which are treated by a monolithic
solver. Sass [Sas04] recently adopted this point of view to establish symbolically
the coupled dynamic equations of an electromechanical system composed of two
subsystems: a rigid mechanism and an electrical network. The strongly coupled
approach is also selected in this thesis, but the assembly is performed numerically,
as in the Finite Element method.
A few methodologies escape to the strict classication between weak and
strong coupling, such as the fastest-rst approach evoked in [GW84], or the sub-
cycling technique proposed by Cardona [Car89] for the dynamic analysis of mech-
anisms with hydraulic actuators. In both cases, smaller time-steps are applied for
the fast subsystem according to a weakly coupled approach, but a strong coupling
is considered at every global time-step.
Element
Coupling strategy Time-integration
description
Element type,
User
Finite Element method
Ports, causality
User Element type, physical
assignment/tree Solver selection
interface properties
Physical modeling
selection
Element Constitutive equations Port variables
library
Analysis ODE/DAE solver
library
Port concept,
Kernel Element concept topological equations,
graph analysis
Solver selection
interface properties coordinates denition
Element Energetic Generalized/redundant
library contributions coordinates
Analysis ODE/DAE solver
library
Relation generalized/
Kernel Element concept
redundant coordinates
User Solvers/scheduler
Partitioned methods
Dof/node
Block diagram concepts,
Kernel Element concept
concept numerical
assembly
- the reliability of the method, since the strongly coupled approach prevents
any diculty due to algebraic loops.
Two important contributions of this thesis are the formulation of the block di-
agram concepts in a Finite Element framework, and the extension of the
generalized-α method for the time-integration of the coupled equations.
This work addresses the control design of exible mechanisms. At this design
stage, it is commonly accepted that an ideal mechanical model should combine
the following properties:
2. low-order, since the order of the controller is often related with the order
of the model,
in both elds of linear system theory and linear structural dynamics. We shall
review them before considering the reduction of a nonlinear exible multibody
system.
Several state-space reduction methods have been proposed for linear time-
invariant systems. The reduction relies on two low-dimensional subspaces: the
former denes a coordinate transformation, and the second, a projection operator
for the state equations. The objective of a reduction technique is to select those
subspaces in order to minimize the accuracy loss. Some methods, well-suited
for large-scale problems, are based on Krylov subspaces [GVVD04], and others
rely on a truncated balanced realization, as presented by Gawronski [Gaw98]. In
order to preserve physical properties of the system, such as dissipativity, congru-
ence transformations4 may be advantageous, especially for passive RLC electrical
circuits and passive mechanical systems [PLS03].
In linear structural dynamics, dedicated reduction techniques exploit the
second-order and lightly damped nature of the equations of motion. Initially,
those methods were not developed for control applications, but for substructured
analysis of complex mechanical systems: the reduced contributions of each sub-
structure are established independently in a rst step, in order to simplify the
analysis of the assembled structure in a second step. A congruence transforma-
tion is dened in the space of generalized coordinates as a mode shape matrix.
Hurty [Hur65], Craig and Bampton [CB68], Herting [Her79] and Craig [Cra87]
proposed various methods to select the truncated modal basis. They are usu-
ally denoted component-mode techniques, and a more detailed presentation will
be given in chapter 5. Interesting combinations of reduction methods in system
theory and in structural dynamics have been investigated by Gawronski [Gaw98]
and De Fonseca [DF00].
The component-mode synthesis leads to an optimal modal representation of
a exible body, which can be exploited in a oating frame formulation, leading
to the concept of superelement [SW83, CG91, GC01]. But again, for a complex
exible multibody system, the resulting model does not meet the requirements
of a control design procedure.
4 In a congruence transformation, the subspaces for the coordinate transformation and for
the projection are identical.
28 CHAPTER 2. STATE OF THE ART
design.
This systematic and general method leads to a highly compact dynamic rep-
resentation of complex exible mechanisms. The resulting model is appropriate
for control design, but it can also be used to build a more structured model, e.g. a
polytopic linear model or a linear fractional model. As a special case, the method
is applicable to parallel rigid mechanisms [BDG06].
of large space structures. Another solution is to avoid extra actuators, and use
the existing actuators for the simultaneous control of the large amplitude motion
and of the vibrations. This is the indirect vibration control strategy, which is
especially suitable for the control of robots and machine-tools. This approach is
cheaper and simpler from a mechanical point of view. However, if one is interested
in controlling the tip-position, the non-collocated conguration of the sensors and
the actuators complicates the design of the control law. This research focuses on
the indirect vibration control problem, which has received a major attention in
the literature.
The control methods reported here rely on a model resulting from spatial
discretization with a nite number of coordinates. Since a exible system has, in
principle, an innite number of dofs, the higher order dynamics is neglected, but
its interaction with the control system may cause spillover instability. Special
care is necessary to prevent the system from such undesirable phenomena.
Most control techniques rely on the structure of a mathematical model, hence,
we successively consider methods based on linear approaches, on quasi-linear ap-
proaches, and nally on nonlinear concepts; of course, each method is conditioned
by the availability of an appropriate model.
There exists a special class of systems, called dierentially at systems, for
which there is a one-to-one correspondence between trajectories of a set of "at
outputs" and full state and input trajectories [FLOR95, VNM98, FLMR99]. This
special invertibility property leads to an ecient trajectory planning tool: a tra-
jectory specied in output space can be lifted to the state and input space through
an algebraic mapping. Ollivier and Sedoglavic [OS01] exploited this tool for the
trajectory planning of a very exible rod, and Thummel et al. [TOB01] for
exible joint robots. Since this formalism requires the analytical structure of the
dynamic equations, its systematic implementation for complex exible mecha-
nisms seems to be a dicult problem.
Feedback linearization is a geometric method to transform a nonlinear system
into a virtual linear system by a state transformation and a feedback transforma-
tion [SL91, Isi95]. Then, linear control theory can be applied for the virtual linear
system. For a rigid robot, this method is equivalent to the classical computed
torque technique. Sometimes, feedback linearization is classied among inverse
dynamics techniques, since it realizes a cancellation of the nonlinear part of the
dynamics. Feedback linearization typically requires a full state estimation, which
is not possible if part of the dynamics is not observable at the outputs. This unob-
servable dynamics is called the zero dynamics and it is interpreted as the system
dynamics when the outputs are forced to zero. By denition, a nonlinear system
is minimum-phase if its zero dynamics is asymptotically stable. Several authors,
such as Wang and Vidyasagar [WV91a], applied those advanced concepts to the
indirect vibration control of mechanisms with exible links, and they concluded
that the full state dynamics cannot be linearized. The choice of the outputs has a
decisive inuence on the stability of the zero dynamics: the collocated joint angles
lead to an oscillatory zero dynamics, whereas the noncollocated tip position leads
to an unstable behavior [WV91a, DLS93]. The output redenition concept can
also be generalized from the linear case [MPK00]. Feedback linearization relies
on an analytical expression of the dynamic equations, so that this tool may be
hardly applied for complex parallel mechanisms. Moreover, some authors [SJK97]
described the pitfalls of the feedback linearization technique, in particular its lack
of robustness, since it may involve useless and expensive cancellation of stabilizing
nonlinearities, using their dangerous destabilizing counterpart.
Singular perturbation theory has also been investigated for the control of ex-
36 CHAPTER 2. STATE OF THE ART
ible multilink manipulators [SB88, MC95, GLL98, MPK00, GS00, SV01, CYC02].
This theory assumes the partitioning of the state variables into slow and fast vari-
ables, so that the dynamics can be represented by two reduced subsystems: the
slow subsystem obtained by neglecting the fast oscillatory motion, and the fast
subsystem obtained by freezing the slow variables. A composite control can then
be applied, combining control laws established separately for each subsytems. In
the context of exible mechanisms, the elastic coordinates are naturally selected
as fast variables. The slow subsystem corresponds to a corrected rigid system
(e.g. including the static contribution of the exible variables), whereas the fast
subsystem has a linear dynamics, parametrically depending on the conguration.
The slow controller can be designed on the basis of well-established joint tracking
schemes for rigid manipulators, and the fast subsystems can be stabilized using
standard methods for linear systems [SB88, GLL98, GS00, SV01] or for linear
parameter-varying systems [MC95, MPK00, CYC02]. This theory is extremely
appealing, but, as pointed out by Book [Boo93], the performances are limited by
the necessary time-scale separation between the rigid and the exible controllers.
The Lyapunov method is a powerful tool for the stability analysis of non-
linear systems. However, for closed-loop systems, the construction of the con-
trol algorithm does not completely follow from this theory and involves engineer
judgement and intuition.
If the transfer function between the input and a chosen output is passive 6 ,
Lyapunov-type control schemes can be more systematically developed with guar-
anteed stability. Usually, this method leads to simple and robust feedback strate-
gies. For a exible mechanism, the passivity property of collocated transfer func-
tions can be easily assessed, leading to a classical joint controller. Dening a
passive output that includes deformation eects is a more dicult problem.
Considering plane mechanisms that are not inuenced by gravity, Ge et al.
[GLZ96] proposed an energy-based feedback which includes deformation eects
and leads to a stable closed-loop system in the sense of Lyapunov. An attractive
property of this strategy comes from the independence of the control design with
respect to the system dynamics, so that no model is necessary. However, a draw-
back comes from its conservatism (i.e. its low performances), and no guideline is
available to optimize the parameters of the controller.
6 An exact denition of passivity in system theory is given by Sepulchre et al. [SJK97].
2.5. CONTROL OF FLEXIBLE MECHANISMS 37
Specic nonlinear control design methods have been developed for systems
with a cascade structure, e.g. the states of one subsystem are the control variables
for the following. A step-by-step recursive design procedure can then be exploited,
such as backstepping or forwarding [SJK97]. Backstepping has been explicitly
applied for exible joint robots [OL99]. Recursivity is also at the basis of the
inertial damping concept developed for the control of macro/micromanipulators.
In the macro/micro conguration, a small rigid robot is mounted at the tip of a
large and exible manipulator. Several approaches have been proposed for the
control of such system, and an extended review is presented by George [Geo02].
Among them, the inertial damping method consists in controlling the motion of
the rigid robot in such a way that the inertia forces transmitted to the supporting
exible manipulator have a damping eect on the vibration modes. This is a true
cascade system: the control of the micromanipulator motion allows the active
damping of the macromanipulator.
chapter 5 will be extremely valuable for the design and the implementation of the
control law.
III
Multibody Dynamics
(a) (b)
(c) (d)
Figure 3.1: (a) Minimal coordinates, (b) Relative coordinates, (c) Cartesian co-
ordinates, (d) Finite Element coordinates or natural coordinates.
x1 = y1 = y6 = 0 and x6 = l4 (3.1)
x2 = x3 , y2 = y3 , x4 = x5 , y4 = y5 (3.2)
3.2 Dynamics
The Finite Element coordinates are absolute coordinates, and the total mo-
tion (rigid-body motion and elastic deformation) is directly referred to an inertial
frame. Due to the large displacements and rotations of the mechanical elements
with respect to this frame, the linear theory of elasticity is not applicable, and a
three-dimensional nonlinear theory is necessary.
An updated Lagrangian point of view for the rotation parameters has been
recommended by Géradin and Cardona [GC01]. This means that the rotations
3.2. DYNAMICS 43
L = K − V, (3.6)
δW = δqT Q (3.7)
and Φ are m scleronomic constraints. The formulation may also include non-
holonomic constraints, but we shall not insist on this point here.
The problem (3.5) can be replaced by the equivalent formulation:
Z t2
δ (L(q, q̇) + W − λT Φ) dt = 0 (3.8)
t1
where λ is the m×1 vector of Lagrange multipliers associated with the constraints.
Performing the variation and the integration by part, one obtains the Lagrange
equations:
d ∂L ∂L
− + ΦTq λ = Q
dt ∂ q̇ ∂q
(3.9)
Φ(q) = 0
where Φq denotes the constraint gradient. Assuming that K can be expressed as
a quadratic form of generalized velocities with a symmetric mass matrix M(q):
1 T
K= q̇ M(q) q̇ (3.10)
2
the equations of motion can be put in matrix form:
(
M q̈ + ΦTq λ = g(q, q̇)
(3.11)
Φ(q) = 0
44 THE FE APPROACH IN MULTIBODY DYNAMICS
where the vector of apparent forces g(q, q̇) collects external forces, internal forces
and complementary inertia forces:
∂V X X ∂Mki 1 ∂Mij
gk = Qk − − − q̇i q̇j (3.12)
∂qk i j
∂qj 2 ∂qk
This system of equations, and its linearized counterpart, can be directly used
for numerical analysis, e.g. in the time domain. Several enhanced formulations
have been proposed for more ecient numerical treatments. Among them, we
shall review the constraint elimination technique and the augmented Lagrangian
method.
q∗ = Φ∗ (θ) (3.13)
Since the m algebraic constraints Φ are nonlinear and implicit, the expression of
Φ∗ cannot be constructed analytically. However, the Jacobian of the transforma-
tion can be computed by implicit dierentiation:
" #
∗ I
∂q ∂q
= −Φ−1
q∗ Φθ ⇒ = qθ = (3.14)
∂θ ∂θ −Φ−1
q∗ Φθ
where k and p are the scaling and penalty factors. The dynamic equations follow:
(
M q̈ + ΦTq (p Φ + k λ) = g
(3.16)
k Φ(q) = 0
Since the penalty term vanishes at the solution point, it is easily observed that
the method provides the exact solution of the initial problem.
Since the time-integration procedure involves a Newton-Raphson procedure,
it is important to formulate the linearized equations for the corrections of the
displacements, velocities and accelerations:
" # " # " # " # " # " #
M 0 ∆q̈ Ct 0 ∆q̇ Kt k ΦTq ∆q
+ +
0 0 ∆λ̈ 0 0 ∆λ̇ k Φq 0 ∆λ
" #
−resq
= Φ
+ O(∆2 )
−res
(3.17)
where
- resq and resΦ denote the residual vectors of equilibrium and constraints,
with:
resq = −g + M q̈ + ΦTq (p Φ + k λ) (3.18)
3.3 Time-integration
In chapter 2, the relevance of the Newmark family of implicit algorithms for
the simulation of exible mechanisms was highlighted. This section presents the
original Newmark scheme [New59] and its generalized-α extension by Chung and
Hulbert [CH93]. Theoretical results will also be detailed, related with the issues
of stability and convergence.
where the constants β and γ are numerical parameters. It may be shown that the
optimal choice of the parameters corresponds to the average constant acceleration
formula, obtained with the particular values:
1 1
γ= and β= (3.23)
2 4
They give maximal second-order accuracy, and unconditional linear stability, as
will be demonstrated later. Numerical damping can be introduced in the New-
mark formula according to:
2
1 1 1
γ = +α and β= γ+ α>0 (3.24)
2 4 2
where α is the numerical damping parameter. This choice allows to increase the
numerical damping in the system while remaining on the stability boundary of the
3.3. TIME-INTEGRATION 47
and to satisfy the Newmark formulae (3.22), the iterative corrections should sat-
isfy:
1
4q̈n+1 = 4qn+1
β h2
γ
(3.26)
4q̇n+1 = βh
4qn+1
so that the linearized system (3.17) becomes:
" # " # " #
Sqq k ΦTq ∆q −resqk
t
= (3.27)
k Φq 0 ∆λ −resΦ
k
Time incrementation
-
t=t+h
Initial prediction
?
Evaluation of residues
?
resq (q, q̇, q̈, λ), resΦ (q)
Incrementation
?
q, q̇, q̈, λ
the Newmark formulae (3.22), whereas the residual equations are modied by
averaging the dierent contributions between both time instants:
(
∗
(1 − αm ) (M q̈)n+1 + αm (M q̈)n + (1 − αf ) gn+1 + αf gn∗ = 0
(3.29)
(1 − αf ) k Φn+1 + αf k Φn = 0
where g∗ is a notation for ΦTq (k λ + p Φ) − g, while αm and αf are numer-
ical parameters. The Hilber-Hughes-Taylor algorithm is obtained for αm = 0
and αf ∈ [0, 1/3]. The optimal parameters of the generalized-α method can be
computed from the desired spectral radius at innity ρq∞ :
2 ρq∞ − 1 ρq
αm = q and αf = q ∞ (3.30)
ρ∞ + 1 ρ∞ + 1
Dening αf m = αf − αm , the Newmark parameters are now given by:
2
1 1 1
γ = + αf m and β= γ+ (3.31)
2 4 2
It will be demonstrated later that this scheme is second-order accurate even for
αf m > 0.
The numerical solution is obtained with a predictor-corrector strategy as for
the Newmark algorithm. At the correction step, the augmented tangent matrix
becomes:
γ 1
Sqq
t = (1 − αf ) Kt + (1 − αf ) Ct + (1 − αm ) M (3.32)
βh β h2
and the linearized equations:
" # " # " #
Sqq k (1 − αf ) ΦTq ∆q −resqk
t
= (3.33)
k (1 − αf ) Φq 0 ∆λ −resΦ
k
Unconstrained system
The stability analysis of an integration scheme applied to a linear system of
equations:
M q̈ + K q = 0 (3.36)
q̈ + ω 2 q = 0 (3.37)
Complementing this equation with the Newmark formulae (3.22), we get the
matrix relation
(1 − αf )ω 2 0 (1 − αm ) q αf ω 2 0 αm q
1 0 −βh2 q̇
+
−1 −h −( 21 − β)h2
q̇ = 0
0 1 −γh q̈ 0 −1 −(1 − γ)h q̈
n+1 n
(3.39)
and we seek its eigensolutions, characterized by a synchronous behavior:
q q0
with (3.40)
q̇ = ϕn q̇0 ϕn+1 = ζ ϕn
q̈ q̈0
n
ϕ is the amplitude, [q0 q̇0 q̈0 ]T is the eigenvector, and ζ is the eigenvalue. If we
dene the polynomials,
P (ζ) = (1 − ) ζ + = {αm , αf }
Pγ (ζ) = γ ζ + (1 − γ)
(3.41)
Pβ (ζ) = β ζ + ( 12 − β)
P1 (ζ) = ζ − 1
the characteristic equation is
Pαf (ζ)ω 2 0 Pαm (ζ)
q
(ζ) = det (3.42)
Pωh P1 (ζ) −h −Pβ (ζ)h2
=0
0 P1 (ζ) −Pγ (ζ)h
which is equivalent to:
q
Pωh (ζ) = Pαf (Pγ + P1 Pβ ) (ωh)2 + P12 Pαm = 0 (3.43)
The algorithm is stable if all the roots of Pωh
q
(ζ) = 0 are inside the unit circle.
Unconditional stability requires this property to be satised whatever the value
of ωh. Therefore, an important property of the algorithm is the spectral radius
associated with the roots ζiq of Pωhq
for ωh → ∞:
ρq∞ = max |ζiq | (3.44)
i
which should be less than one for unconditional stability. This spectral radius is
computed from the characteristic equation:
q
P∞ = Pαf (Pγ + P1 Pβ ) = 0 (3.45)
52 THE FE APPROACH IN MULTIBODY DYNAMICS
If we consider the optimal parameters β and γ given by (3.31), this equation has
one simple root and one double root:
−αf 1 − αf m
ζ1q = q
ζ2,3 =− (3.46)
1 − αf 1 + αf m
Constrained system
The linearized equations of an undamped system are:
M q̈ + K q + BT λ = 0
(3.47)
Bq = 0
As pointed out by Géradin and Cardona [GC01], this system is not diagonalizable
due to the presence of the constraints. However, these authors demonstrated that
it could be transformed into a canonical form by a linear transformation:
q̈r + Ω2 qr = 0
q̈c + qλ = 0 (3.51)
qc = 0
Following the procedure detailed for the unconstrained case, the characteristic
equation is the determinant of the following matrix:
Pαf Ω2 0 0 0 0 0 Pα m I 0 0
0 0 Pα f I 0 0 0 0 Pα m I 0
0 Pα f I 0 0 0 0 0 0 0
P1 I 0 0 −h I 0 0 −Pβ h2 I 0 0
0 P1 I 0 0 −h I 0 0 −Pβ h2 I 0
0 0 P1 I 0 0 −h I 0 0 −Pβ h2 I
0 0 0 P1 I 0 0 −Pγ h I 0 0
0 0 0 0 P1 I 0 0 −Pγ h I 0
0 0 0 0 0 P1 I 0 0 −Pγ h I
(3.52)
After lengthy but straightforward algebraic manipulations, and referring to the
characteristic polynomial Pωhq
associated with the unconstrained problem, one
nds that the characteristic equation is equivalent to:
Pαf (ζ) Ω2 0 Pαm (ζ) I
q
(P∞ (ζ))2m det −h I −Pβ (ζ)h2 I = 0 (3.53)
P1 (ζ)I
0 P1 (ζ)I −Pγ (ζ)h I
54 THE FE APPROACH IN MULTIBODY DYNAMICS
The second factor can be interpreted as the characteristic equation when the
algorithm is applied to an unconstrained diagonal system, so that equation (3.53)
can be restated:
(3.54)
Y q
q
(P∞ (ζ))2m Pωi h (ζ) = 0
i
where ωi are the diagonal terms of Ω. Any root of this equation is either a
root of the polynomial Pωqi h obtained in the unconstrained case, or a root of the
same polynomial at innite frequency. We conclude that if the generalized-α
algorithm leads to a stable integration of the unconstrained problem, and if the
spectral radius at innity ρq∞ is smaller than 1, the global stability is guaranteed.
If ρq∞ = 1, it can be demonstrated that the numerical integration is weakly
unstable [GC01].
In this work, qn+1 , q̇n+1 , and q̈n+1 conventionally denote the numerical solutions,
whereas q(tn+1 ), q̇(tn+1 ), and q̈(tn+1 ) refer to the exact solutions. The numerical
solutions qn+1 , q̇n+1 satisfy the Newmark formulae (3.22), and the dierences
with the exact solutions q(tn+1 ), q̇(tn+1 ) are the local truncation errors:
with
α0 = −αm , α1 = −1 + 2αm , α2 = 1 − αm ,
µ0 = −αm , µ1 = −1 + αm ,
(3.63)
β0 = (1/2 − β) αf , β1 = 1/2 − β − αf /2 + 2 β αf , β2 = β (1 − αf ),
γ0 = (1 − γ) αf , γ1 = 1 − γ − αf + 2 γ αf , γ2 = γ (1 − αf ).
and a further elimination of the velocities in the rst equation would lead to a
three-step formula for the displacements. Equation (3.62) is a good basis for the
analysis of the local truncation error at the velocity level. Since it is a standard
multistep formula, the order 2 condition can be found in classical textbooks such
as Hairer et al. [HNW87]:
α0 + α1 + α2 = 0 (3.64)
α1 + 2 α2 = γ0 + γ1 + γ2 (3.65)
α1 + 4 α2 = 2 γ1 + 4 γ2 (3.66)
These relations come from the substitution of q̇ and f by their Taylor series
expansion in the expression of the local error. The rst two equations are auto-
matically satised by the parameters given in (3.63), whereas the last equation
yields
1
γ = + αf − αm (3.67)
2
56 THE FE APPROACH IN MULTIBODY DYNAMICS
ẋ = f (x, t) (3.69)
solved using a k-step stable method, whose local truncation error satises the
order p condition
kx(tn+1 ) − xn+1 k ≤ M hp+1 (3.70)
We also dene the (k nx ) × 1 vectors collecting the coordinate vectors of the last
k steps:
The following bound on the global error after n time-steps is given by Hairer et
al. [HNW87]
M hp nhL∗
(3.73)
∗
kX(tn ) − Xn k ≤ kX(t0 ) − X0 k enhL +
e − 1
L∗
where L∗ is a Lipschitz constant. An initial error in the initial conditions is at
most amplied by a coecient depending exponentially on the length of the time
interval nh. Therefore, it is especially important to limit this phenomenon by the
denition of consistent initial conditions.
Since the generalized-α method hides a multistep algorithm, the error prop-
agation might be characterized by a similar property, and the consistency of the
initial conditions is of practical interest to guarantee the global convergence. For
a multistep method, the initial conditions are given as the vector X0 of the solu-
tion at the k rst steps, which is necessary to initiate the integration algorithm.
However, in its pseudo-one-step formulation, the generalized-α algorithm starts
on the basis of the values q0 , q̇0 and q̈0 . For full consistency, q̈0 should satisfy:
In practice, this equation is helpless for the determination of q̈0 , and it is usually
replaced by the simplied consistency condition:
q̈0 = f0 (3.75)
This expression satises (3.74) with an O(h) error. From a dimensional analysis,
this initial error at the acceleration level will contaminate the displacements and
velocities during the rst steps with O(h3 ) and O(h2 ) errors, respectively. This
error in the initial steps will be coupled, propagated and amplied by the integra-
tion procedure, leading to a maximal global error of O(h2 ), so that second-order
convergence is still guaranteed. However, if the simplied consistency condi-
tion (3.75) is not satised, the error in the initial values may drop one order,
with disastrous consequences on the global convergence of the results.
Φ(q) = 0 (3.76)
are satised at each time-step, if they were satised by the initial conditions.
Hidden constraints should also be satised by the dynamic variables at all dier-
entiated levels, and in particular, at the velocity and acceleration levels:
Φ̇ = Φq q̇ = 0 (3.77)
Φ̈ = Φ̇q q̇ + Φq q̈ = 0 (3.78)
Those equations are, however, not considered in the algorithm, which might lead
to signicant errors.
In the special case of linear constraints, Φq is constant, and the constraints
at the acceleration level become:
Φq q̈ = 0 (3.79)
After inspection of the predictor and corrector steps of the numerical algorithm,
the constraints at displacement, velocity and acceleration levels are satised at
tn+1 if they were initially satised at tn . Therefore, the local truncation error is
of order 2, as in the unconstrained case.
Since the constraints are satised at each time grid point, the constraint vi-
olation of the numerical solution is propagated at a very high frequency which
depends on the step size. In agreement with the stability analysis, the ampli-
cation factor of this phenomenon is associated with the spectral radius of the
integration algorithm at high frequencies. Therefore, the numerical damping is
responsible for the decay of this error.
In case of nonlinear but smooth constraints, if we do not consider nonlinear
amplication eects, the violation of the hidden constraints resulting from the
nonlinearity may be seen as a disturbance, which is ltered by the algorithm in
the same manner. In this sense, those errors do not accumulate throughout the
integration process, and the convergence results demonstrated for unconstrained
systems are still relevant.
From a practical perspective, the computation of consistent initial veloci-
ties and accelerations may rely on equations (3.77) and (3.78) (or its linearized
counterpart (3.79)).
3.4. IMPLEMENTATION 59
3.4 Implementation
The concepts presented in this chapter have been implemented in the Mecano
program of the Samcef software [SAM99]. In 2005, this program has reached
industrial maturity for several years, it is well documented and it contains a huge
library of nite elements and types of analysis, which is very convenient for the
user. Mecano has been extensively exploited for problems in aerospace, robotics,
automotive engineering and machine-tools.
However, for the developments realized in this research, the Oofelie multi-
physics platform [CKG94] has been selected for several reasons. Oofelie is an
acronym for Object Oriented Finite Element Led by Interactive Executor. It is
written in C++, and modern programming concepts are favourable for an open
architecture, well-suited for new developments. Its kernel has been developed
to deal with strongly coupled multiphysics problems, with great care about the
modularity and eciency issues. Since the developments related with exible
multibody dynamics only started a couple of years ago, a restricted library of the
most signicant elements and algorithms is currently available in Oofelie. Nev-
ertheless, those capabilities are sucient to demonstrate the relevance and the
eciency of the innovative concepts developed in this research.
60 THE FE APPROACH IN MULTIBODY DYNAMICS
IV
Systems
?
Generalized-α solver
ga- wm
Mechanism
Control system
diagram language, and on the simulation of the resulting state-space model. The
simulation of the coupled system will be considered later, in section 4.2.
f s∗ (wm , x, ẋ, t) = 0
(4.1)
f o∗ (wm , x, ga , t) = 0
where the nx × 1 vector x represents the state variables, f s∗ are the nx dierential
equations of the dynamic states, whereas f o∗ are the na algebraic output equa-
tions. Control engineers are more familiar with the explicit state-space format,
which is a special case of the descriptor state-space format:
ẋ = f s (wm , x, t)
(4.2)
ga = f o (wm , x, t)
In general, the transformation from (4.1) to (4.2), is not always possible nor
trivial.
Such a global input/output black box description is compact and ecient.
However, a functional decomposition into subsystems may lead to an advanta-
geous modular approach, as illustrated in Figure 4.3. At the subsystem level,
the explicit state-space equations can be formulated more easily, and we consider
that each element e is characterized by explicit state equations, with respect to
64 CHAPTER 4. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS
ga- wm
Mechanism
Figure 4.3: Modular approach for the description of the control system.
ẋ = f s (u, x, t) (4.6)
y = f o (u, x, t) (4.7)
u = Lim wm + Lio y (4.8)
u- y
System
In this case, the implicit function theorem can be invoked to solve (4.7) and (4.8)
for the inputs and outputs:
and the numerical algorithm behaves as if it were applied to the equivalent ODE:
ẋ = a x + b u
y = cx + u (4.17)
u = y
The input/output equations hide the constraint x = 0, so that the output variable
y = u plays the role of a Lagrange multiplier in the dynamic equation. The block
diagram formalism presented here is not really appropriate for such a higher-
index DAE. In the following, Assumption 4.1 is supposed to be satised, so that
an ecient strategy can be developed with guaranteed reliability.
4.1. SIMULATION OF CONTROL SYSTEMS 67
s
(1 − δm ) ẋn+1 + δm ẋn − (1 − δf ) fn+1 − δf fns = 0 (4.20)
o
(1 − δf ) (yn+1 − fn+1 ) + δf (yn − fno ) = 0 (4.21)
ẋ0n+1 = 0
x0n+1 = xn + (1 − θ) h ẋn (4.25)
0
yn+1 = yn
The linearized form of the discretized state equations (4.20), combined with
the linearized input/output relations are:
which can be put in matrix form using (4.26), after elimination of ∆u and ∆ẋ:
" # " # !" # " #
1 I 0 −fxs −fus Lio ∆x −ressk
(1 − δm ) + (1 − δf ) =
θh 0 0 −fxo I − fuo Lio ∆y −resok
(4.30)
In chapter 3, the linear stability analysis of the generalized-α method was in-
vestigated for second-order undamped equations. We may reproduce this analysis
for rst-order equations.
4.1. SIMULATION OF CONTROL SYSTEMS 69
Time incrementation
-
t=t+h
Initial prediction
?
Evaluation of residues
?
ress (x, ẋ, y), reso (x, y)
Incrementation
?
x, ẋ, y
ϕ is the amplitude, [x∗0 ẋ∗0 ]T is the eigenvector, and ζ is the eigenvalue. The
characteristic equation is
" #
Pδf (ζ)σ Pδm (ζ)
x
Pσh = det =0 (4.40)
P1 (ζ) −Pθ (ζ)h
where ζ1x is the spurious root (equivalent to the spurious root ζ1q associated with
the simulation of a mechanical system), and ζ2x is the principal root. With the
optimal parameters for second-order accuracy θ = 1/2 + δf m , we get
1 − 2 δf m
ζ2x = − (4.45)
1 + 2 δf m
Finally, in the system (4.35), the algebraic equation associated with an output
y∗ = 0 (4.46)
The single eigenvalue of the amplication matrix is the spurious root ζ1x .
As a conclusion, |ζ1x | < 1 requires δf < 1/2, and |ζ2x | < 1 requires δf m > 0,
so that the stability is guaranteed for:
1
δm < δf < (4.48)
2
Equivalent results were obtained by Jansen et al. [JWH00]. It is easily demon-
strated that the optimal condition |ζ1x | = |ζ2x | = ρx∞ leads to
3 ρx∞ − 1 ρx∞
1
δm = and δf = (4.49)
2 ρx∞ + 1 ρx∞ + 1
These formulae are dierent from the optimal formulae (3.30) obtained for the
simulation of a purely mechanical system.
with Liq = Lim Lmq , Liq̇ = Lim Lmq̇ , and Liq̈ = Lim Lmq̈ . The generalized forces
of the actuators ga produce the mechanical virtual work
where Lqa and Lqo = Lqa Lao are the actuator and output localization matrices.
The whole set of coupled equations is then:
Equation (4.53) represents the dynamics of the mechanical system, equation (4.54),
the kinematic constraints, equation (4.55), the state dynamics, equations (4.56)
the algebraic output equations, and (4.57) the input localizations.
s s
(1 − δm ) ẋn+1 + δm ẋn − (1 − δf ) fn+1 − δf fn =0
o
(1 − δf ) (yn+1 − fn+1 ) + δf (yn − fno ) = 0
n+1 = xn + h (1 − θ) ẋn + h θ ẋn+1
x
(4.58)
with the discretized inputs:
un+1 = Liq qn+1 + Liq̇ q̇n+1 + Liq̈ q̈n+1 + Lio yn+1 (4.59)
where g∗ is a notation for ΦTq (k λ + p Φ) − g(q, q̇) − Lqo y. The rst set of
equations contains the discretized dynamic equations associated with the me-
chanical dofs and the Newmark formulae. The second contains the discretized
state equations and their time-integration formula. At this level, two remarks
may be formulated, which suggest a reformulation of those equations.
Remark 4.1 In order to uniformize the treatment of the state and displacement
variables, the dummy dynamic variables z are introduced:
Z t
z(t) = x(τ ) dτ (4.60)
0
The value of z is only meaningful at the velocity level (ż = x) and at the accelera-
tion level (z̈ = ẋ); nevertheless, an articial displacement-like Newmark formula
4.2. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS 75
Remark 4.2 According to Remark 3.1, q̈n+1 is a poor approximation for q̈(tn+1 ).
Therefore, we exclude the acceleration term from the input equation,
u = Liq q + Liq̇ q̇ + Lio y (4.62)
so that yn+1
a
is an order 2 approximation for q̈ i (tn+1 ), which can be safely con-
nected to any input port. The global output equation becomes
y = f o (u, x, t) + Loq̈ q̈ (4.65)
The prediction for qn+1 , zn+1 , q̇n+1 , żn+1 follows from the standard Newmark
formulae with the zero acceleration assumption q̈0n+1 = 0, z̈0n+1 = 0. Using the
linearized form of the input equation (4.62):
γ iq̇
∆u = iq
L + L 4q + Lio ∆y (4.67)
βh
where
γ 1
St = (1 − αf ) Kt + (1 − αf ) Ct + (1 − αm ) Mt (4.69)
βh βh2
Kt , Ct , Mt are respectively given by
Kt k ΦTq 0 −Lqo Ct 0 0 0 M 0 0 0
k Φq 0 0 0 0 0 0 0 , 0 0 0 0
,
−fus Liq 0 0 −fus Lio s iq̇
−fu L s
0 −fx 0 0 0 I 0
−fuo Liq 0 0 I − fuo Lio o iq̇
−fu L o
0 −fx 0 −Loq̈
0 0 0
(4.70)
The correction equation involves non-symmetric matrices. Thanks to the
penalization term, Kt is positive denite and it is possible to apply a direct
solver without pivoting strategy. Good performances were observed with a non-
symmetric direct solver optimized for sparse matrices. This choice is reasonable
since a mechatronic system usually involves no more than a few hundreds dofs.
For systems with much more dofs, iterative solvers could be advantageously ap-
plied.
The time-integration algorithm is described in Figure 4.6. In order to analyze
its stability and convergence properties, it is equivalent to analyze the dynamic
system after elimination of the input and output variables from the state equa-
tions. Under Assumption 4.1, the input and output variables dened by (4.62)
and (4.65) can be formulated in explicit format (4.14), (4.15), and replaced in the
4.2. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS 77
Time incrementation
-
t=t+h
Initial prediction
?
Evaluation of residues
?
resq , resΦ , ress , reso
Incrementation
?
where fes (x, q, q̇, q̈, t) = f s (feu (x, q, q̇, q̈, t), x, t).
The matrices M f , C and K include the possible contributions of the direct feed-
back control actions feq̈o , feq̇o and feqo , respectively (e.g. M
f 6= M). Hence, Mf, C
and K are not necessarily symmetric positive denite, which is an important
dierence compared to the mass, damping, and stiness matrices of a passive
mechanical system.
The velocity term C q̇ may come from the internal damping of the material,
or from a direct feedback action. It is expected that this latter contribution has a
stabilizing eect, otherwise, the whole dynamic system may get unstable as well
as the numerical simulation. Therefore, the case C = 0 can be seen as the worst
case situation, where stability in the simulation results is desirable; the damping
matrix C is thus omitted in this analysis.
The diculty may be reduced by diagonalization of matrices M f −1 K and
A:
Ω2 = T−1 q M
f −1 K Tq and Σ = T−1
x A Tx (4.75)
4.2. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS 79
q = Tq q∗ and x = Tx x∗ (4.76)
Γ = T−1 −1 −1 −1
q G Tx , Π = Tx P T q , ∆ = Tx D Tq , Υ = Tx F Tq (4.77)
where the coupling matrices Γ, Π, ∆ and Υ have no specic structure, that could
be exploited for a component-wise analysis.
Those equations are discretized according to (4.58):
(1 − αf ) Ω2 0 (1 − αm ) I (1 − αf ) Γ 0 q∗
∗
I 0 −βh2 I 0 0
q̇
∗
0 I −γh I 0 0 q̈
(1 − α ) Π (1 − α ) ∆ (1 − α ) Υ (1 − α ) Σ (1 − α ) I x∗
f f m f m
∗
0 0 0 I −γh I ẋ
n+1
αf Ω2 0 αm I αf Γ 0 q∗
∗
−I
−h I −( 21 − β)h2 I 0 0
q̇
∗
+ 0 −I −(1 − γ)h I 0 0 q̈ = 0
α Π α ∆ αm Υ αf Σ αm I x∗
f f
∗
0 0 0 I −(1 − γ)h I ẋ
n
(4.79)
where I and 0 still denote identity and null matrices of appropriate dimensions.
We seek the eigensolutions of this dierence equation, characterized by a
synchronous behavior:
q∗ q∗0
∗
q̇
q̇∗0
∗
q̈ = ϕn
∗
q̈0 with ϕn+1 = ζ ϕn (4.80)
x∗ x∗0
∗ ∗
ẋ ẋ0
n
80 CHAPTER 4. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS
In order to analyze the stability in the presence of sti dynamics in the mechanical
systems and/or in the control system, the characteristic equation is developed in
three subcases:
1. ωh → ∞:
−h −Pβ (ζ)h2 0 0
P1 (ζ) −Pγ (ζ)h 0 0
Pαf (ζ) ω det (4.83)
2
=0
Pαf (ζ) d Pαm (ζ) f Pαf (ζ) σ Pαm (ζ)
which is equivalent to P∞
q x
(ζ) Pσh (ζ) = 0.
4.2. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS 81
2. σh → ∞:
Pαf (ζ) ω 2 0 Pαm (ζ) 0
−Pβ (ζ)h2
P1 (ζ) −h 0
Pαf (ζ) σ det (4.84)
=0
0 P1 (ζ) −Pγ (ζ)h 0
0 0 0 −Pγ (ζ)h
3. ωh → ∞ and σh → ∞:
−h −Pβ (ζ)h2 0 0
P1 (ζ) −Pγ (ζ)h 0 0
Pαf (ζ) ω det
2
Pαf (ζ) d Pαm (ζ) f Pαf (ζ) σ Pαm (ζ)
would be possible if the stiness matrix K were real symmetric positive semi-
denite, and the mass matrix M f , real symmetric positive denite. Therefore, the
developments of this section assume that K and M f satisfy those conditions, as in
the passive case. This is always the case if no direct acceleration or displacement
feedback is present.
Using the notations of sections 3.3.3 and 4.2.3, we dene
Γr = TTqr G Tx , Πr = T−1 r −1 r −1
x P Tqr , ∆ = Tx D Tqr Υ = Tx F Tqr
Γc = TTqc G Tx , Πc = T−1 c −1 c −1
x P Tqc , ∆ = Tx D Tqc Υ = Tx F Tqc
(4.87)
We obtain the following set of equations:
q̈r + Ω2 qr + Γr x∗ = 0
q̈c + qλ + Γc x∗ = 0
(4.88)
qc = 0
ẋ∗ + Σ x∗ + Πr qr + Πc qc + ∆r q̇r + ∆c q̇c + Υr q̈r + Υc q̈c = 0
Assumption 4.2 The functions fes (x, q, q̇, q̈, t) and feo (x, q, q̇, q̈, t) are linear in q̈:
∗
fes = fes (ż, q, q̇, t) + feq̈s q̈ (4.94)
∗
feo = feo (ż, q, q̇, t) + feq̈o q̈ (4.95)
These equations have the same structure than the equations of a purely me-
chanical system. Hence, the same convergence properties are observed, and we
immediately conclude that the local truncation error is still of order 2.
84 CHAPTER 4. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS
Remark 4.3 Assumption 4.2 is satised for any linear control system and any
control system without acceleration measurement, and those categories cover many
practical situations.
Remark 4.4 If the functions fes and feo are not linear but only ane with respect
to q̈ ( i.e., in equations (4.94) and (4.95), feq̈s and feq̈o are not constant), the argu-
ments described in section 3.3.4 about the consequences of a non-constant mass
matrix convince us that the local truncation error might increase by one order.
For instance, such a situation arises when acceleration signals wma are obtained
from a 3-axis accelerometer xed on a moving body. The measurements are the
absolute accelerations in body axes, and they are related to the accelerations in
inertial axes q̈ by the rotation matrix of the body R:
wma = RT q̈ (4.98)
In case of very fast motion, the variations of R may not be negligible over each
time-step, leading to a small loss in accuracy.
Remark 4.5 Control systems involving non-ane acceleration feedbacks are quite
unusual for control practitioners, so that we should not be afraid about the restric-
tions involved by Assumption 4.2.
1
|ζq2,3|
0.9 |ζx2|
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 0.2 0.4 0.6 0.8 1
αfm
This explains why the optimal choices of αf and αm are dierent for a mechanical
model and a state-space model (compare formulae (3.30) and (4.49)).
According to Chung and Hulbert [CH93], the optimal parameters for the
simulation of the mechanical equations are obtained when |ζ1 | = |ζ2,3q
|. From
Figure 4.7, we then conclude
which imposes the design constraint αf < 41 . In this case, the spectral radius
associated with the mechanical variables, is always higher than the spectral radius
associated with the state variables: ρq∞ > ρx∞ .
Hence, the optimal choice of the algorithmic parameters for the coupled set
of equations is not trivial; in the numerical applications presented in this chapter,
the Hilber-Hughes-Taylor method is systematically applied.
q̈− +
n+1 6= q̈n+1 , λ− +
n+1 6= λn+1 , ẋ− +
n+1 6= ẋn+1 ,
−
yn+1 +
6= yn+1 (4.105)
q− +
n+1 = qn+1 , q̇− +
n+1 = q̇n+1 , x− +
n+1 = xn+1 (4.106)
The integration scheme developed in the previous section was able to predict
eciently the values at time t−n+1 . Since the integration formulae rely on a Taylor
4.3. SYSTEMS WITH DISCONTINUOUS DYNAMICS 87
expansion of the dynamic variables, they are not applicable over the discontinuous
transition from t−n+1 to t+n+1 , and the correction for q̈, λ, ẋ and y should be
performed independently of q, q̇ and x, which are kept constant.
This correction can be seen as an integration restart procedure, and we know
from earlier convergence results (see section 3.3.4) that the initial conditions at
time t+n+1 should satisfy the residual equations:
M q̈+ T + + +
n+1 + Φq (k λn+1 + p Φ) − g(qn+1 , q̇n+1 , t) − L
qo +
yn+1 = 0 (4.107)
k Φ(q+
n+1 ) = 0 (4.108)
ẋ+ s + + +
n+1 − f (un+1 , xn+1 , tn+1 ) = 0 (4.109)
+
yn+1 − f o (u+ + +
n+1 , xn+1 , tn+1 ) = 0 (4.110)
u+ iq + iq̇ + iq̈ + io +
n+1 − L qn+1 − L q̇n+1 − L q̈n+1 − L yn+1 = 0 (4.111)
otherwise, the order of convergence would be deeply aected. q+n+1 , q̇+n+1 and
n+1 are imposed by (4.106), and these equations can be seen as nonlinear alge-
x+
braic equations for q̈+n+1 , λ+n+1 , ẋ+n+1 , yn+1
+
. They can be solved using a Newton-
Raphson procedure with the initial values given at time t−n+1 . However, since the
constraint equation does not involve any of these unknowns, the system is ill-
dened. In order to overcome this problem, the constraints should be formulated
at the acceleration level:
k Φ̈ = k Φ̇q q̇ + k Φq q̈ = 0 (4.112)
with Φ̇q = Φ̇q (q, q̇) and Φq = Φq (q). If this hidden constraint is satised at
time t−n+1 , it implies:
k Φq (q̈+ −
n+1 − q̈n+1 ) = 0 (4.113)
which may replace equation (4.108) in order to obtain a nonsingular system1 .
The linearized equations follow:
M k ΦTq 0 −Lqo 4q̈ −resqk
4λ −resΦ
k Φq 0 0 0 k
(4.114)
=
0
0 I −fus Lio
4ẋ −ress
k
−Loq̈ 0 0 I − fuo Lio 4y −resok
1 The mechanical equations are linear in q̈, λ and y, thus, if the control system is also linear,
one Newton iteration yields convergence.
88 CHAPTER 4. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS
Time incrementation
- -
t=t+h
no
?
Update Initial prediction
yes
Discontinuity ?
x, f s , f o t− q0 , q̇0 , q̈0 , x0 , ẋ0
? ?
Evaluation of residues Evaluation of residues
resq , resΦ , ress , reso resq , resΦ , ress , reso
? ?
Check for convergence Check for convergence
yes yes
t + kresa k < a , a = {q, Φ, s, o} kresa k < a , a = {q, Φ, s, o}
no ? no ?
Evaluation of corrections Evaluation of corrections
4q̈, 4λ, 4ẋ, 4y 4q, 4λ, 4x, 4y
? ?
Incrementation Incrementation
q̈, λ, ẋ, y q, q̇, q̈, λ, x, ẋ, y
the values of q̇ should simply be updated when the impulse is detected, afterwards
the correction algorithm presented in Figure 4.8 is applicable. Therefore, the
critical point is to compute q̇+n+1 .
If ge denotes the generalized impact forces, their integral eect pe is dened
by
Z t+
pe = ge dt (4.117)
t−
where λp are Lagrange multipliers that are related to the internal impact forces
(impact reactions). The rst equation represents the conservation of momen-
tum, and the second equation states that the velocity increment has to fulll the
homogeneous velocity constraint equations.
The computation of q̇+n+1 requires the knowledge of pe , which is not always
a trivial problem. If the impulse forces result from an impact between bodies,
pe can be estimated using a physical model of their contact (e.g. an elastic or a
plastic model).
Element
4
DiscontinuousSystem ContinuousSystem
Kernel 4 4
Elements
LTISampledSystem UserSampledSystem
Figure 4.9: Element library: class diagram. LTI is an acronym for Linear Time
Invariant.
- the concepts of output and state equations, and their connection with the
Finite Element assembly procedure.
A simplied view of the library of control elements appears in Figure 4.9. The
developments realized in the kernel make the programming eort for an element
minimal, since only the functions f s and f o , and their sensitivity matrices have to
be dened (the implementation of the sensitivities can be avoided using a nite
dierence algorithm).
One may not be surprised to nd the SaturatedIntegrator among the
discontinuous elements, and the Saturation among continuous elements:
- the SaturatedIntegrator is characterized by one state variable x which
integrates the output until the saturation value xmax :
and y is continuous.
4.5 Applications
Three examples with increasing complexity are considered: a rigid four-bar
mechanism, a Scara robot, and a vehicle with semi-active suspensions.
ẋ = θ − θref (4.123)
y = −P (θ − θref ) − D θ̇ − I x (4.124)
θref - θ -
τ
- PID - Mechanism
θ̇
-
ref
xk+1 = xk + T s (θk+1 − θk+1 ) (4.126)
ref
yk+1 = −P (θk+1 − θk+1 ) − D θ̇k+1 − I xk+1 (4.127)
with the sampling period T s = 0.01 s. The discrete state xk is not treated as a dof
in the simulation, it is simply updated at every sampling period. The maximal
value of the time-step is equal to the sampling period: h(max) = T s .
In this example, the mechanism is initially at rest (q0 = qinit , q̇0 = 0),
and the time variation of the external loads and reference inputs is plotted in
Figure 4.12. We may notice that the initial values q̈ = 0 and ẋ = 0 (for the
analog control) are consistent with those initial conditions.
The concise and high level textual model denition in Oofelie is illustrated
in Figure 4.13. For the digital case, the syntax is similar, the user should simply
further specify the sampling period.
4.5. APPLICATIONS 93
g ext 6 θref 6
θtar
θinit
0 h(max) 0 h(max)
- -
t t
Figure 4.12: Four-bar mechanism with PID control: external loads and reference
input. h(max) = 0.01 s.
[Link](PID,ContinuousLTISystem_E,3,1,1);
ElemSet[PID].addInput(ThetaRefNode,GenDisp);
ElemSet[PID].addInput(HingeAngleNode,GenDisp);
ElemSet[PID].addInput(HingeAngleNode,GenVel);
ElemSet[PID].addOutput(TorqueNode);
ElemSet[PID].setStateSpaceMatrices(A,B,C,D);
Figure 4.13: Oofelie input data le - description of an analog PID controller
(3 inputs, 1 state and 1 output). PID is the element number; ThetaRefNode,
HingeAngleNode, TorqueNode are node numbers; GenDisp, GenVel mean "gen-
eralized displacement" and "generalized velocity" respectively; and A,B,C,D are
the state-space matrices associated with equations (4.123) and (4.124).
94 CHAPTER 4. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS
−3
x 10
0
1.55
−1
1.5 −2
−3
PID state
1.45
−4
θ (rad)
−5
1.4
−6
−7
1.35
Analog PID
Digital PID −8
Digital PID without correction
1.3 −9
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
t (s) t (s)
Simulation results
A small numerical damping is necessary to avoid high frequency oscillations
due to the mechanical algebraic constraints, and we have selected αf = 0.05,
αm = 0 (Hilber-Hughes-Taylor algorithm).
A rst simulation is realized with a time-step h(1) = h(max) = 0.01 s, and the
results are presented in Figure 4.14. A noticeable dierence is observed whether
the system is controlled by an analog or a digital controller. For the digital
control case, a simplied integration scheme has also been tested, where the
specic correction step after discontinuities is omitted. We observe that the
error in the simplied algorithm and the sampling eect have the same order of
magnitude, and we conclude that the correction step is essential for the accurate
representation of the sampling eect.
In order to analyze the inuence of a time-step reduction, ve dierent values
have been considered:
h(max)
h(p) = p = {1, ..., 5} (4.128)
2p−1
In Figure 4.15, the results obtained with h(1) and h(5) are very close to each other
for both the analog and the digital controller. However, for the digital controller,
the simplied algorithm without discontinuity correction leads to important er-
rors, especially for h(1) .
Assuming that the results obtained with h(5) are close to the exact solution,
4.5. APPLICATIONS 95
1.365
Analog PID / h1
1.36 Analog PID / h5
Digital PID / h1
1.355 Digital PID / h5
Digital PID without correction / h1
1.35 Digital PID without correction / h5
θ (rad)
1.345
1.34
1.335
1.33
1.325
0 5 10 15 20 25 30
t (ms)
Figure 4.15: Four-bar mechanism with PID control - inuence of the time-step
(h(1) = 10 ms, h(5) = 0.625 ms).
where qn(p) is the value of the dynamic variable predicted with a time-step h(p) at
time tn . The evolution of σ with respect to the time-step is plotted in Figure 4.16.
The expected second-order convergence is observed in all cases, excepted for the
simplied algorithm which exhibits a poor rst-order convergence. Therefore, the
simplied algorithm is denitely not acceptable for the simulation of a mecha-
tronic system with a digital controller. For the corrected algorithm, we conclude
that the order of convergence predicted by the theory is veried.
−3 −2
10 10
−3
−4
10
10
−4
10
−5
σ
σ
10
−5
10
−6
10 −6
10
The electrical current is also limited to a maximal value, and U em is the electro-
motive voltage proportional to the joint velocity:
U em = k em θ̇ (4.132)
The electromechanical equation denes the relation between the torque and the
current intensity:
T = kt i (4.133)
θ̇ -
θref -
U a- Electrical Saturated i - Electromech. T -
di
θ - PD U ref
- Limitation dt -
Figure 4.18: Scara robot - Controller and motor associated with one joint.
98 CHAPTER 4. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS
2.5 2.5
θ1
θ2
2 2 θ3
ref ref
θ1 = θ2
ref
1.5 θ3 1.5
(rad and m)
θ (rad and m)
1 1
ref
θ
0.5 0.5
0 0
−0.5 −0.5
0 0.5 1 1.5 2 0 0.5 1 1.5 2
t (s) t (s)
Figure 4.19: Scara robot, joint references (step signal) and joint coordinates.
7 7
i i1
1
6 i2 i
6 2
i i
3 3
5
5
4
4
3
i (A)
i (A)
3
2
2
1
1
0
−1 0
−2 −1
0 0.5 1 1.5 2 0 0.02 0.04 0.06 0.08 0.1
t (s) t (s)
Figure 4.20: Scara robot - electrical currents. The picture on the right is a zoom.
7
i1 0.2
6 i2 with correction
0 without correction
i3
5
−0.2
4
−0.4
3
i (A)
−0.6
i (A)
2
−0.8
1 −1
0 −1.2
−1 −1.4
−2 −1.6
0 0.5 1 1.5 2 0.98 0.99 1 1.01 1.02 1.03
t (s) t (s)
Figure 4.21: Scara robot with sampled controller - electrical current in motor 1.
The picture on the right is a zoom centered on t = 1s.
100 CHAPTER 4. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS
From this example, we conclude that our simulation method is able to deal sys-
tematically and reliably with the coupled electromechanical equations of a mecha-
tronic system. The block diagram model is modular and convenient for the dy-
namic description of the controller and the actuators. Most methods presented
by other authors were based on an adaptive time-stepping method, which is a
strong advantage over our current xed time-step implementation. It is therefore
hopeless to attempt any comparison with respect to the computational eciency,
and this weakness of our software should be addressed in future investigations.
However, the decisive advantage of our approach certainly comes from its gen-
erality, in particular its ability to deal eciently and accurately with complex
parallel topologies and mechanical deformations; those potentialities are not re-
ally highlighted in this simple example. The objective of the next section is to
demonstrate the relevance of the integrated simulation method for a mechanism
with a far more complex topology, subject to the action of a strongly nonlinear
control system.
r Controller
ab , v- Actuators ga Mechanical ab , lr , vr
iv - - -
model lr , v r - model model
Figure 4.22: Mechatronic model of the car equipped with a semi-active suspen-
sion. ab is the vector of car-body accelerations, lr is the vector of rattle extensions,
vr = l̇r is the vector of rattle velocities, iv is the vector of electrical currents, and
ga is the vector of damper forces.
relevant basis for the numerical optimization of the suspension, which will be
addressed in future research.
Mechanical model
The car is an Audi A6 (see Figure 4.23), and its rigid-body model includes
the following components, illustrated in Figure 4.24:
102 CHAPTER 4. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS
Figure 4.23: Audi A6. On the right, instrumentation with corner accelerometers.
- the car-body,
- the wheel model for the lateral force, the vertical force and the yaw torque.
The current model involves about 600 mechanical dofs, and it could be extended
to include the stiness of the suspension bushings, the exibility of the chassis,
and a longitudinal model for the wheels. However, the current model is sucient
to demonstrate the potentialities of the integrated simulation tool.
4.5. APPLICATIONS 103
Actuator model
A systematic description of the actuator model and of the controller is pre-
sented by Lauwerys et al. [LSS04].
Figure 4.25 compares a passive and a semi-active shock absorber. The passive
absorber contains two passive valves, restricting the oil ow from one chamber to
the other. Due to the motion of the rod in and out of the cylinder, the variation
of the total volume v tot available for the oil is:
where v reb and v comp are the volumes of the rebond and compression chambers.
The role of the accumulator is to compensate for the variation of v tot .
In the active shock absorber, the passive valves are replaced by check valves,
and a current controlled CVSA valve4 connects the extreme chambers. When
the rod moves up, the piston check-valve closes and the oil moves out of the
rebond chamber through the CVSA valve. When the rod moves down, the base
check-valve closes and the oil ows to the accumulator through the CVSA valve.
Therefore, both motions are aected by the controlled restriction at the CVSA
valve.
The ow variations reported in Figure 4.25 are related with a quasi-static
motion, where the pressures are in equilibrium at every instant, the uid is not
compressible, and the accumulator is ideal. Actually, the dynamic behavior is
more complex, and a nonlinear model has been calibrated by the manufacturer of
the shock-absorber. This model can be presented in nonlinear state-space format,
dening the inputs, states and output:
with lr , the rattle extension, iv , the electrical current in the CVSA valve, preb and
pcomp the pressures in the rebond and compression chambers, and g a the force
exerted by the damper. The state equations are:
Figure 4.25: Passive and semi-active damper. Quasi-static valve ows during
compression and extension sequences. The proportions on the drawings are not
realistic and, for information, typical values for the chamber diameter and the
rod diameter are 32 mm and 22 mm, respectively.
4.5. APPLICATIONS 105
vr1 vr1
vr1
iv1 iv1
gvirt1
gvirt1 vr2
vr Inverse Actuator 1
vr2 vr3
vr2
iv2 iv2
gvirt2
gvirt2 vr4
Inverse Actuator 2
vr3 ab1
vr3
iv3 iv3
gvirt3
gvirt3 ab2
gvirt Inverse Actuator 3
vr4 ab3
vr4
iv4 iv4
gvirt4
gvirt4 ab4
Inverse Actuator 4
Reconnect
Car + Dampers
gvirt Modal Virt Forces Modal Virt Forces Modal Acc Modal Acc Coupled Acc
Figure 4.26: Controller of the semi-active suspension. "vr" stands for the rattle
velocities, "ab" for the car-body accelerations, and "iv" for the valve currents.
The functions f s(damp) and f o(damp) are given by the manufacturer as C-functions,
which are linked with Oofelie and called in a user-element inheriting from the
class ContinuousSystem (see Figure 4.9). The sensitivities of those functions are
obtained using a nite dierence procedure.
Controller model
The control law of the active suspension has been developed by Lauwerys,
Swevers, and Sas [LSS04]. As illustrated in Figure 4.26, it consists of three stages:
a feedback linearization (inverse actuator models), a transformation into modal
space (coupling and decoupling operations), and a linear integral control.
The car and the dampers are represented as a black box, whose inputs are
the four CVSA electrical currents iv , and whose outputs are the rattle velocities
v r = l˙r , and the accelerations measured at the four corners of the car ab (see
106 CHAPTER 4. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS
Figure 4.23).
The feedback linearization technique seeks for virtual inputs which have the
property to inuence the outputs in a linear way. If one accepts that the non-
linearity of the mechanism is weak, the main source of nonlinearity lies in the
actuator. According to Lauwerys et al., an ecient feedback linearization law is
obtained by inversion of a simplied quasi-static model of the actuators. From
equation (4.138), it is possible to formulate the quasi-static damper force with
respect to the CVSA electrical current, and the rattle velocity:
g a = g a (iv , v r ) (4.139)
This relation can be inverted for v r 6= 0:
iv = (g a )−1 (g virt , v r ) (4.140)
where g virt is the new virtual input, which can be interpreted as a virtual damper
force. For v r = 0, the controllability is defective, and the singularity of (g a )−1 is
avoided thanks to a regularization strategy. If this inverse model is able to cancel
the nonlinearity of the actuator, the virtual input g virt is actually proportional to
the damper force. Therefore, besides the advantage of a good linearity between
g virt and the output, a force control strategy can be established on this basis.
Since the motion of the car simultaneously involves the forces applied on the
four wheels, the denition of the virtual control forces g virt for the four shock-
absorbers is a linear but coupled multi-input/multi-output problem. This prob-
lem can be simplied by a transformation into a modal space dened by the
heave (pumping), roll and pitch of the car-body. In this modal space, the sys-
tem is represented by three uncoupled single-input/single-output subsystems, for
which three independent integral controllers are designed.
The several stages of the controller are easily described using the block dia-
gram language. The inverse actuator model is implemented as a specic element,
which directly invokes the C-function implemented in the actual controller, and
the sensitivities are obtained by analytical dierentiation. All other blocks are
modeled using the element library mentioned in Figure 4.9.
Simulation results
The simulation of a lane change maneuver has been realized, and the target
trajectory, represented in Figure 4.27, corresponds to a standard qualication
4.5. APPLICATIONS 107
Figure 4.27: Lane change maneuver. The dimensions presented here correspond
to a standard qualication test.
0.05 0.5
0.04 0
slider−crank translation (m)
0.03 −0.5
0.02 −1
0 −2
−0.01 −2.5
−0.02 −3
−0.03 −3.5
−0.04 −4
−0.05 −4.5
0 1 2 3 4 5 6 0 10 20 30 40 50 60
time (s) x (m)
test. The car has a 10 m/s initial velocity, and a driver applies an open-loop
steering command, without any real-time correction (blind driver assumption).
The motor and the brakes do not produce any torque on the wheels. The time-
step for the simulation is 0.01 s, and the algorithmic parameters are αf = 0.05
and αm = 0.
Figure 4.28 illustrates the motion of the slider-crank mechanism actuated by
the driver, and the horizontal trajectory of the car.
Figure 4.29 represents the global motion of the car-body. The yaw, pitch
and roll angles of the car-body are plotted as well as the radius of gyration and
the vertical displacements. The radius of gyration becomes innite when the car
follows a straight line.
The dynamic behavior of the semi-active shock absorbers is illustrated in Fig-
ure 4.30. It is interesting to observe the pressures in the rebond and compression
chambers, depending on the sign of the rattle velocity: during the compression,
108 CHAPTER 4. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS
0.5 0.06
yaw roll
0.4 pitch
0.04
0.3
0.2 0.02
angle (rad)
−0.1 −0.02
−0.2 −0.04
−0.3
−0.06
−0.4
−0.5 −0.08
0 1 2 3 4 5 6 0 1 2 3 4 5 6
time (s) time (s)
50
0.02
45
40 0.015
vertical displacement (m)
Gyration radius (m)
35
0.01
30
25
0.005
20
15 0
10
−0.005
5
0 −0.01
0 1 2 3 4 5 6 0 1 2 3 4 5 6
time (s) time (s)
5
x 10
5
0.04
Rear right
Front right 4
0.03
Extension displacements (m)
0.02 3
Pressure (Pa)
0.01 2
0
1
−0.01
0
−0.02
−1
−0.03 Rebond chamber
Compression chamber
−0.04 −2
0 1 2 3 4 5 6 0 1 2 3 4 5 6
time (s) time (s)
1.8
400
1.6
300
1.4
Electrical current (A)
200
1.2
100
Force (N)
0.8 0
0.6 −100
0.4 −200
0.2
Rear right −300
Rear right
0 Front right Front right
−400
0 1 2 3 4 5 6 0 1 2 3 4 5 6
time (s) time (s)
Figure 4.30: Car semi-active dampers - extension, hydraulic pressure (rear right),
electrical current and damper force.
the valve between those chambers is open, and the pressure in the compression
chamber is slightly higher, but during the extension, the valve is blocked, and
the pressure in the rebond chamber can be much higher than in the compres-
sion chamber (see Figure 4.25). The saturation eect dominates the behavior of
the electrical current in the CVSA valves. The static contribution of the forces
produced by the actuators is not zero, due to non-equilibrated pressures in the
dierent chambers when the piston rod is at rest.
In this application, the integrated simulation method was able to predict the
behavior of a complex mechatronic system. The model consists of a full model
of the car and its suspension mechanisms, a nonlinear dynamic model of the hy-
draulic actuator, and the nonlinear control algorithm. The modularity of the
110 CHAPTER 4. INTEGRATED SIMULATION OF MECHATRONIC SYSTEMS
- the adaptation of the block diagram concept for the modular representation
of rst-order state equations within a Finite Element formulation,
- the extension of the generalized-α method for the strongly coupled simula-
tion of a mechatronic system represented by a mechanical Finite Element
model and a block diagram model, with stability and convergence analyses,
Multibody Dynamics
5.1 Introduction
The reduction method is an extension of the component-mode technique
established in structural dynamics, which accounts for the nonlinear kinematics
of the mechanism. Basically, the reduction procedure proceeds in two steps:
112 CHAPTER 5. NONLINEAR MODEL REDUCTION
- an order reduction of the initial Finite Element model, which can be realized
locally for any given conguration,
?
Next cfg
?
Kinematic analysis
?
Local mode synthesis
?
Model reduction
no
?
All cfg treated?
yes
?
Approximation
Figure 5.1: Simple approach for the model reduction of a exible mechanism, and
illustration with a four-bar mechanism.
114 CHAPTER 5. NONLINEAR MODEL REDUCTION
- accuracy in the bandwidth of the actuators and within the limits of the
workspace,
- low-order and free from kinematic constraints, which follows from the Global
Modal Parameterization,
M q̈ + K q = g (5.1)
q=Ψη (5.2)
n<n (5.3)
The variations of q are thus restricted to the subspace spanned by the component-
modes, and the accuracy of the reduced model depends on the ability of the
component-modes to describe the actual motion of the system.
In order to formulate the reduced equations of motion, the coordinate trans-
formation is introduced into the expressions of the kinetic and potential energies,
and of the virtual work of the external forces. Initially, we have:
1
K = q̇T M q̇ (5.4)
2
1
V = qT K q (5.5)
2
δW = gT δq (5.6)
with the reduced mass matrix, stiness matrix, and equivalent force vector:
M = ΨT M Ψ, K = ΨT K Ψ, g = ΨT g (5.10)
M η̈ + K η = g (5.11)
Φq q = 0 (5.12)
5.2. LINEAR COMPONENT-MODE SYNTHESIS 117
u = Ψuη η (5.16)
qr = θ (5.22)
This property does not hold for the constraint modes: qg 6= η γ . If the modal
masses of the internal modes are normalized, the reduced matrices and forces
have the following structure:
θθ θγ θι
0 0 0 M M M gθ gr + ΨgθT gg
γγ γθ γγ γι
, gγ =
K=
0 K 0 , M = M
M M gg
2 ιθ ιγ ιι ι
0 0 Ω M M I g 0
(5.23)
where Ω = diag(ωi ) is the diagonal matrix of internal eigenvalues.
Géradin and Rixen [GR97] interpreted a class of component-mode methods as
a truncation in the modal expansion of the mechanical impedance. In this sense,
the component-modes suggested by Hurty and Craig-Bampton are optimal.
In multibody dynamics, the linear component-mode technique is usually ex-
ploited for the compact kinematic description of an isolated exible body with
respect to a oating frame of reference, see section 2.1.2. In order to obtain a
more drastic reduction, we apply the modal parameterization to the whole mech-
anism, which is seen as a "component" of the mechatronic system. Therefore,
the concept of component-mode is replaced by the concept of local mode dened
around a conguration. The description of the conguration of a mechanism with
a suitable parameterization is addressed in the next section.
120 CHAPTER 5. NONLINEAR MODEL REDUCTION
Φ(q) = 0 (5.24)
The conguration space Ωrtq is dened as the set of kinematically admissible con-
gurations which satisfy the kinematic constraints:
Ωrt n
q = {q ∈ R | Φ(q) = 0} (5.25)
The superscript r stands for "rigid", and t for the "total" conguration space.
Assuming independent constraints, the number of kinematic modes s satises:
s=n−m (5.26)
ρ : Ωrt rt
θ → Ωq , θ 7→ q = ρ(θ) (5.27)
ϕ : Ωrt rt
q → Ωθ , q 7→ θ = ϕ(q) (5.28)
where Ωrt θ is the set of possible variations of the parameters θ . We also refer
to Ωrt
θ as the conguration space, whenever no confusion is possible. Figure 5.2
illustrates those denitions.
In general, the existence of a global parameterization with independent coor-
dinates is not possible in the total conguration space Ωrtq , and some restrictions
are necessary to dene a pragmatic solution.
We propose to dene the minimal coordinates θ as the actuated dofs, i.e. the
dofs associated with the generalized forces exerted by the actuators. For instance,
for a motorized hinge, the actuated dof is the angle between the connected links,
whereas for a linear actuator, it is the relative distance between the connected
bodies. As a consequence of this choice, the actuator dofs will appear explicitly
in the reduced model, which is extremely valuable for the design of the control
system. An implicit assumption is that the number of actuators is equal to
s. This does not imply that the reduction method is only applicable to fully
122 CHAPTER 5. NONLINEAR MODEL REDUCTION
actuated mechanisms, it only means that in any other situation, the denition of
the independent parameters is left to the user.
The actuated dofs θ are usually relative coordinates, in contrast with the
absolute coordinates q used in our Finite Element formulation. A systematic
implementation of the relation θ = ϕ(q) can be dened using mixed coordinates.
This means that the actuator coordinates θ appear explicitly among the set of n
mixed coordinates q: " #
θ
q= (5.29)
q∗
q∗ are the n − s non-actuated dofs. Since m = n − s, the number of non-actuated
dofs equals the number of kinematic constraints of the mixed formulation.
The mapping ϕ is directly characterized:
θ = ϕ(q) = [I 0] q (5.30)
The regularity of the Jacobian is not sucient to guarantee that the actua-
tor dofs θ are able to parameterize the whole conguration space. A well-dened
parameterization is globally one-to-one: for each given actuator conguration θ,
there exists only one kinematically admissible conguration q (i.e. one solution
to equations (5.31)). This condition is quite dicult to verify; however, in many
practical cases, actuator singularities separate parts of the conguration space
where ϕ is one-to-one, so that the actuator parameterization is valid in a sub-
set Ωrq ⊂ Ωrt q bounded by the singular congurations. All those concepts are
illustrated in Figure 5.3.
In this restricted part of the conguration space, the inverse map ρ is well-
dened:
where Ωrθ ⊂ Ωrt θ ⊂ R is the set of possible variations of the actuated dofs θ ,
s
away from the singularities. For usual control applications, actuator singularities
and conguration space singularities are carefully avoided, and the limitation of
the analysis to Ωrq is not restrictive.
124 CHAPTER 5. NONLINEAR MODEL REDUCTION
Figure 5.3: Rigid four bar mechanism, with an actuator at the bottom left hinge.
The conguration space is a 1-dimensional manifold, whose projection in the
plane of the coordinates ze and θ is also represented. Conguration (b) is an
actuator singularity; it separates two sub-domains of the conguration space for
which q = ϕ(θ) is one-to-one. ϕ is not one-to-one in the total conguration
space Ωrt
q , since congurations (a) and (c) are possible for the same angle θ .
5.3. PARAMETERIZATION OF THE RIGID KINEMATICS 125
Φ(q) = 0 (5.35)
The exible manifold Ωfqt is dened as the set of kinematically admissible cong-
urations which satisfy the kinematic constraints:
This problem can be analyzed numerically using the zero-strain approach. Ac-
cordingly, any equilibrium conguration achieves a minimum of the elastic poten-
tial energy, so that the kinematic problem becomes:
Given θ, nd q∗ :
min
∗
V (5.39)
q
2 The dimension here refers to the intrinsic dimension of the manifolds, as opposed to the
dimension n of the ambient space.
5.3. PARAMETERIZATION OF THE RIGID KINEMATICS 127
subject to
Φ(θ, q∗ ) = 0 (5.40)
According to the augmented Lagrangian method, the following equivalent uncon-
strained problem is considered:
min
∗
V + k λT Φ + p ΦT Φ (5.41)
q ,λ
where λ is the m × 1 vector of Lagrange multipliers, k and p are the scaling and
penalty factors, respectively. The stationarity condition leads to n + m nonlinear
equations: (
∂V
+ ΦTq∗ (k λ + p Φ) = 0
∂q∗
(5.42)
k Φ(q) = 0
Equations (5.42) contains as many equations as unknowns (q∗ ,λ), and can be
solved using standard methods.
The initial problem was to establish q (or q∗ ) from θ; it is equivalent to
construct the relation between u (or u∗ ) and θ, with:
θ " #
θ
(5.43)
u= ∗
q = u∗
λ
In the space of augmented coordinates u, the total conguration space is denoted
u . Equation (5.42) gives a condition for static equilibrium which is a necessary
Ωrt
but not sucient condition for an undeformed conguration: some solutions u
of (5.42) might not satisfy (5.38), and be out of Ωrt u . Typically, those solutions
are associated with pre-stressed congurations (λ 6= 0 or ∂V ∂θ
6= 0), which are not
considered in this research.
If the Jacobian of the system (5.42) is not singular, the implicit function
theorem guarantees the local existence of a function ρ such that u = ρ(θ). As
in the rigid case, the global validity of the actuator parameterization can only be
guaranteed in a restricted part of the conguration space Ωru ⊂ Ωrt u , where the
singularities are avoided. The parameterization mapping is dened by:
ρ : Ωrθ → Ωru , θ 7→ u = ρ(θ) (5.44)
where Ωrθ ⊂ Ωrt
θ ⊂ R is the set of authorized variations of the actuated dofs θ .
s
F(u) = 0 (5.51)
where the n reduced coordinates η are dened in an open set Ωfη ⊂ Rn , and
u = qT λT belongs to a submanifold of the exible manifold: Ωu ⊂ Ωfu . The
T f
(5.56)
f
χ : Ωfη → Ωu , η 7→ u = χ(η)
(5.57)
f
dim Ωu = n < dim Ωfu = n − m
where Ωrθ has been dened earlier and Ωfδ is a neighborhood of the origin.
Equation (5.59) can be dierentiated:
∂(Ψuδ η δ )
δu = χη δη = ρθ + δθ + Ψuδ δη δ (5.61)
∂θ
so that for an undeformed conguration (η δ = 0):
where θ are the actuator dofs, qg are the constraint dofs, and ui are the internal
dofs, including the Lagrange multipliers.
The mode shapes are selected as an optimal basis to describe the linearized
motion around an undeformed conguration, with zero velocities :
" # " # " # " # " #
M 0 ∆q̈ Kt k ΦTq ∆q 0
+ = (5.64)
0 0 ∆λ̈ k Φq 0 ∆λ 0
The approximated expression of the tangent stiness (5.48) is thus exact and can
be safely exploited.
As in the linear case, three subsets of modes are dened:
The rigid modes Ψuθ are dened from the static equilibrium (5.18); they are also
the Jacobian of the rigid kinematics: Ψuθ = ρθ . The constraint modes Ψuγ are
the solutions to the static problem (5.19), where the rigid dofs are xed, and
132 CHAPTER 5. NONLINEAR MODEL REDUCTION
Figure 5.5: Local modes of an unconstrained mechanism : 2 rigid modes and one
exible mode.
the internal modes Ψuι are dened according to the internal eigenvalue prob-
lem (5.20), where the rigid and constraint dofs are xed.
For local consistency, those modes should belong to the tangent space of the
exible manifold at point p: Tp Ωfu . Indeed, using (5.18), (5.19), (5.20) and the
structure of the tangent stiness matrix, it is easy to verify that
Φu Ψuη = 0 (5.66)
- it is globally one-to-one,
5.4. THE GLOBAL MODAL PARAMETERIZATION (GMP) 133
Figure 5.6: Local modes of an unconstrained mechanism: one rigid mode and
two exible modes.
Figure 5.7: Local modes of a constrained mechanism: one rigid mode, and one
exible mode.
134 CHAPTER 5. NONLINEAR MODEL REDUCTION
- it is continuously dierentiable,
- its rank is maximal (i.e. the Jacobian χη is of rank n),
From the previous developments, χ(η) involves two contributions:
- the kinematic mapping ρ(θ), which is well-dened in the conguration space
of interest, as demonstrated in section 5.3; in particular, it is continuous
and dierentiable,
- the exible modes Ψuδ (θ) have only received a local denition at a given
conguration; their variations in the conguration space deserve a careful
analysis.
The following analysis demonstrates that the three consistency criteria are satis-
ed by the Global Modal Parameterization.
One-to-one criterion
The exible submanifold Ωu is dened as the image of Ωfη mapped by the
f
2. the resulting nι∗ eigenmodes are sorted according to their eigenvalue, and
the rst nι < nι∗ modes are selected,
3. the modes are normalized with respect to the mass matrix so that:
with the modal mass µk = ψ Tk Mii ψ k > 0. Therefore, each eigenvalue ωk exhibits
smooth variations with respect to the coordinate θl .
From the continuity of ωk , and of the mass and stiness matrices, the rst
term in (5.71) is continuous, and so is the second one. However, the derivatives of
the eigenmodes ∂ψ k
∂θl
can take arbitrary large values in the kernel of Kii − ωk2 Mii .
If ωk has a multiplicity κ = 1, ψ k is the single vector in this kernel, which
means that the arbitrary variations are parallel to the initial vector. Those vari-
ations are forbidden due to the normalization condition (5.70), and we conclude
that ∂ψ k
∂θl
is continuous.
This property is no more guaranteed for a multiplicity κ > 1. In this case, the
kernel of Kii −ωk2 Mii can be a multi-dimensional subspace, where the eigenvectors
can take arbitrary variations.
An example is given in Figure 5.8. In the one-dimensional conguration
space, there is one conguration θ1 = θb for which ω1 = ω2 and κ = 2. For θ1 < θb
136 CHAPTER 5. NONLINEAR MODEL REDUCTION
with " #
θ − θ0
∆η = (5.74)
ηδ
138 CHAPTER 5. NONLINEAR MODEL REDUCTION
∂ 2 V qη ∂ 2 V
ΨqηT Kt Ψqη = ΨqηT Ψ = (5.80)
∂q2 ∂η 2
and we conclude:
ηη
K = ΨqηT Kt Ψqη (5.81)
Since by construction of the rigid modes Kt Ψqθ = 0, the equivalent stiness
has the structure " #
ηη 0 0
K = δδ (5.82)
0 K
with
δδ
K = ΨqδT Kt Ψqδ (5.83)
This last formula can be directly exploited in the Finite Element code to compute
the equivalent stiness. Finally, the potential energy is represented by a quadratic
form in the exible coordinates:
1 δT δδ
V(θ, η δ ) = η K (θ) η δ (5.84)
2
5.5. REDUCED-ORDER MODEL 139
The contribution ∂V
∂θ
is quadratic in the amplitude of deformation; it is neglected
under the small deformation assumption.
q̇ = χη η̇ (5.87)
we obtain:
1 T ηη
with
ηη
K= η̇ M η̇ M = χTη M χη (5.88)
2
The generalized inertia forces are dened by:
d ∂K ∂K
giner = − (5.89)
dt ∂ η̇ ∂η
A simplied expression is possible if the variations of the reduced mass matrix
M are neglected:
ηη
ηη
giner = M η̈ (5.90)
For a more thorough analysis, let us consider a partitioning of the generalized
coordinates into translation and rotation dofs:
" # " #
x χx (η)
q= = (5.91)
ψ χψ (η)
Thus, the last term in (5.96) vanishes. It is convenient to dene the third order
curvature tensor Γχx :
∂ 2 χxi
with (5.99)
x x x
Γχ
ijk = Γχ χ
ijk = Γikj
∂ηj ∂ηk
so that
(5.100)
x x
χ̇xη = Γχ χ
ij ijk η˙k = Γ . η̇ ij
where the ”.” operator is dened as a product of a third order tensor by a vector.
Hence, the translation inertia forces are concisely formulated
(5.101)
xT x
gtrans xx
χxη η̈ + χxη T Mxx Γχ . η̇ η̇
iner = χη M
- the equivalent forces associated with the curvature (or the nonlinearity) of
the coordinate reduction formula.
5.5. REDUCED-ORDER MODEL 141
α = χα (η) and α̇ = χα
η η̇ (5.106)
χΩη = T χα
η (5.107)
we have
Ω = χΩη η̇ (5.108)
and the developments realized for the translation kinetic energy can be repro-
duced:
T
∂(χΩη η̇)
grot,A
iner =χ Ωη T
Jχ Ωη
η̈+χ Ωη T
J χ̇ Ωη Ωη
η̇+ χ̇ − J χΩη η̇ (5.109)
∂η
142 CHAPTER 5. NONLINEAR MODEL REDUCTION
χΩη is not a Jacobian, and the last term does not vanish in this case. Let us
develop:
∂χα ∂ 2 χα
∂Tiq
Ωη q q
(5.110)
χ̇ ij
= + Tiq η̇k
∂ηk ∂ηj ∂ηj ∂ηk
∂χα ∂ 2 χα
∂ ∂Tiq
χΩη η̇ i =
q q
(5.111)
+ Tiq η̇k
∂η j ∂ηj ∂ηk ∂ηj ∂ηk
where
∂Tiq ∂Tiq ∂αl 1 ∂χα
= = iql l
(5.112)
∂ηk ∂αl ∂ηk 2 ∂ηk
Introducing the antisymmetric third-order tensor ΓΩ :
1 ∂χα α
q ∂χl
ΓΩ
ijk = iql with ΓΩ Ω
ijk = −Γikj (5.113)
2 ∂ηj ∂ηk
we obtain
(5.114)
α
Ωη
ΓΩ Γχ
χ̇ ij
= ijk + ijk η̇k
∂
(5.115)
χ α
χΩη η̇ i = ΓΩ
ikj + Γijk η̇k
∂ηj
(5.117)
χα
T
grot,A αT α αT
. η̇ η̇ + 2 ΓΩ . η̇ J χα
iner = χη J χη η̈ + χη J Γ η η̇
hΩ
ijk = Jjl lik (5.118)
- the inertia forces obtained when the variations of the reduced mass matrix
are quasi-static,
5.5. REDUCED-ORDER MODEL 143
- the equivalent forces associated with the curvature (or the nonlinearity) of
the coordinate reduction formula,
ggyr,A αT gyr,A
(5.120)
Ω
iner = χη h . α̇ α̇ = χαT
η giner
This generalized formulation for the inertia forces of an isolated rigid body
is also valid in the particular case where no coordinate transformation is applied:
(5.121)
α
η = α, χα
η = I, Γχ = 0
Let us demonstrate the connection of equations (5.119) with the classic Euler
equations in this case. For a rigid body described by a diagonal inertia tensor
J1 0 0
(5.122)
J=
0 J2 0
0 0 J3
hΩ
123 = −J2 , hΩ
132 = J3 ,
hΩ
213 = J1 , hΩ
231 = −J3 , (5.123)
hΩ
312 = −J1 , hΩ
321 = J2
χψT
η M
ψψ
χψ
η η̈ (5.126)
with
M = χTη M χη (5.130)
χ gyr
h=h +h (5.131)
The expression of the reduced mass matrix is classical, and two dierent third-
order tensors are responsible for quadratic forces in the modal rates.
χ
(i) The tensor h is associated with the curvature Γχ of the reduced param-
eterization:
χ
hijk = χTη M il Γχ (5.132)
ljk
gyr
(ii) The tensor h satises
gyr
. η̇) η̇ = χψT (5.133)
gyr
h . (χψ ψ ψT gyr
(h η η η̇) (χη η̇) = χη g
The mass matrix and the gyroscopic tensor are only evaluated at undeformed
congurations.
The curvature tensor is required for the construction of h. From the ex-
pression of the Jacobian (5.61), it has an interesting structure at an undeformed
conguration:
∂ 2 χi ∂ 2 ρi ∂Ψuθ
= =
ij
(5.135)
∂θj ∂θk ∂θj ∂θk ∂θk
2 uδ
∂ χi ∂Ψik
δ
= (5.136)
∂θj ∂ηk ∂θj
2
∂ χi
= 0 (5.137)
∂ηkδ ∂ηjδ
It is associated with the sensitivities of the rigid and exible modes, which are
easily estimated by nite dierence.
146 CHAPTER 5. NONLINEAR MODEL REDUCTION
(5.139)
M(θ) η̈ + h(θ). η̇ η̇ + K(θ) η = gext
An overview of the local order reduction algorithm is given in Figure 5.9. This
local analysis is combined with a conguration space approximation algorithm,
presented in the next section. Hence, a database is exploited to store and reuse the
local models. As mentioned earlier, the mode matching may fail if the reference
conguration θref is too far from θ, leading to a non-valid model.
5.5. REDUCED-ORDER MODEL 147
no
Kinematic analysis
?
u = ρ(θ)
no
yes-
Ψ uη
Valid model?
Mode matching
?
no-
MAC(Ψuι , Ψuι Non-valid model
ref )
yes
Model reduction
?
M, K, h
?
Valid model
f = f (θ) ⇒ f 'b
f (θ) ∀ θ ∈ Ωrθ (5.141)
2. run the local reduction algorithm for each input and get back f (1) ...f (r) ;
Step 1 is also referred to as the experimental design problem. Step 2 was the
subject of our previous developments. In step 3, many choices are possible, and
we restrict the study to linear approximations:
nv
(5.142)
X
fbi (θ) = Wij vj (θ)
j=1
5.6. APPROXIMATION IN THE CONFIGURATION SPACE 149
where vj (j = 1, ..., nv ) are xed basis functions, and Wij are the free weights
which allow to t the model. The linearity is with respect to the weights, and
not to the input variables. Step 4 is then a linear regression problem which can
be solved in the least-square sense using standard methods.
Radial basis functions are an interesting candidate for the basis functions:
vj = vj (kθ − θ (j) k) j = 1, ..., r (5.143)
If the number of radial functions equals the number of given congurations θ(j) ,
it is possible to interpolate exactly any set of scattered data. However, for large
data sets, the resulting function bf may become computationally inecient, and
the number of basis functions should be reduced following a non-trivial selection
procedure. Moreover, the choice of appropriate radial functions is not systematic,
and sometimes requires nonlinear optimization strategies.
In this research, we prefer low-order polynomial basis functions, which are
more simple and ecient for our problem. In order to increase the exibility of
the approximation, a piecewise strategy is adopted as proposed in [BDG06], so
that the above procedure is applied for several non-overlapping subsets of the
conguration space Ωrθ . The decomposition into subsets can be realized adap-
tively in order to satisfy a specication on the approximation error. Thus, a
general, ecient and systematic approximation procedure is developed leading to
a portable and computationally ecient nonlinear model.
The next section presents the piecewise strategy for approximation. After-
wards, the polynomial basis functions are dened and an algorithm is proposed
for automatic conguration space decomposition.
Figure 5.10: Inner and outer approximation of the conguration space using
subpavings.
Interval analysis theory [JKDW01] oers systematic methods for the approxi-
mation of a complex set using simple subsets. If [θ] ⊂ R represents a real interval,
with a lower and an upper bounds:
[θ] = [θ θ] (5.144)
an interval vector [θ] ⊂ Rs is dened by the cartesian product of s real inter-
vals [θi ]:
[θ] = [θ1 ] × [θ2 ] × . . . × [θs ] (5.145)
Roughly speaking, an interval vector is a s-dimensional box that can be repre-
sented by two corners: the lower bound θ = [θ1 θ2 . . . θs ]T and the upper bound
θ = [θ1 θ2 . . . θs ]T .
A subpaving X is a set of nb non-overlapping boxes4 :
X = {[θ]1 , . . . , [θ]nb } (5.146)
The conguration space Ωrθ can be approximated from inside and outside by two
subpavings, see Figure 5.10:
X − ⊂ Ωrθ ⊂ X + (5.147)
where vj are a low-order polynomials, dened for each boxes in local coordinates.
Hence, the construction of an approximated piecewise model in the conguration
space involves two problems:
In the following, possible choices for polynomial approximation are rst dis-
cussed; the construction of an appropriate subpaving is addressed later.
Quadratic polynomials
For a one-dimensional problem, a quadratic polynomial P1 is dened by
P1 = a + b θ + c θ 2 (5.150)
where vj(k) = vj (θ(k) ). This overdetermined set of equations can be solved in the
least square sense, using standard regression algorithms (normal equations, QR
algorithm or Singular Value Decomposition).
Several authors (e.g. Montgomery [Mon97]) argued that the data points θ(k)
should be selected according to an experimental design in order to improve the
quality of the approximation. Figure 5.11 illustrates classical denitions: the full
factorial design, the fractional factorial design, and the central composite design.
More sophisticated methods, such as the Taguchi approach, seem unnecessarily
complicated for our problem.
The central composite design is retained here, since it provides sucient
information to estimate the quadratic eects required for a second-order poly-
5 For a function f : Rs → Rt , the values of θi θj are computed once for the t output com-
ponents, and the associated computational cost is thus far less signicant than for the other
operations.
5.6. APPROXIMATION IN THE CONFIGURATION SPACE 153
nomial approximation. The axial position of the "star points" is selected at the
intersection with the border of the box.6
Due to the piecewise strategy, the approximation function exhibits a discon-
tinuous behavior at every boundary between boxes. This important drawback
can be overcome using the family of Lagrange polynomials (see [ZT89, PTVF92]
for detailed presentations), so that a C0 continuity is obtained at every boundary
between matching boxes.
Lagrange polynomials
Let us consider the approximation of a real function f : R → R, from p data:
f (k) = f (θ(k) ), k = 1, ..., p. It is well-known that an order p − 1 polynomial can
be tted exactly on those data, for instance using the p Lagrange polynomials
L(i) , which are constructed to satisfy the condition:
δik is the kronecker symbol. This condition leads to the Lagrange interpolation
formula: p p
θ − θ(k)
(5.155)
Y X
(i)
L (θ) = = Aik θk−1
θ(i) − θ(k) k=1
k=1
(k 6= i)
The last equality denes the components of the matrix A, which is implemented
in the numerical code for ecient computations. Therefore, the interpolating
polynomial of the function f is simply
p p
(5.156)
X X
P1 (θ) = f (i) (i)
L (θ) = f ∗(k) θk
i=1 k=1
with p
(5.157)
X
f ∗(k) = f (i) Aik
i=1
θ2 6
(3)
θ2 × × ×
(2)
θ2 × × ×
(1)
θ2 × × ×
-
(1) (2) (3)
θ1 θ1 θ1 θ1
with
p
∗(k)
(5.159)
X (i)
Ps−1 = Ps−1 Aik
i=1
From this analysis, Lagrange polynomials are less attractive for high-dimensional
problems. Anyway, a compromise has to be dened between:
- the accuracy,
- the continuity,
Ideally, the size of the boxes should be adapted according to the behavior of f , in
order to improve the compromise between accuracy, continuity, memory storage
and ecient construction of the approximation function. For instance, close to
singularities, a renement is suitable to track the strong variations of the dynamic
parameters. The adaptive algorithm presented here relies on the representation
of the subpaving as tree data structure and it allows a high-level denition of the
approximation error.
Before the description of the algorithm, let us dene the algorithmic object
model, at the basis of the conguration space inspection procedure. A model is
a data structure containing:
- two pointers leftModel, and rightModel to child models that might be de-
ned after the bisection of the box [θ], thus, a model is binary tree structure
which represents a subpaving,
- eval and ebi (i = 1, ..., s) are representative of the approximation error in the
box, they will be dened later,
Figure 5.14 illustrates the construction of a model associated with a box [θ].
Local models are dened for every data points of the conguration set cfgSet, and
after, the approximation function is constructed. The model is bisected if the local
analysis failed at any data point or if the approximation error is not acceptable.
In those cases, a bisection of the model leads to the recursive denition of two
child models. In the following, the error analysis, and the selection of the bisection
direction are presented with more details.
Error analysis
The relative approximation error is dened by:
f (θ) − f (θ)
b
e(θ) = >0 (5.165)
kf (θ)k
Since the components of the vector f do not have the same physical meaning
and the same order of magnitude, it is recommended to replace the norm by a
weighted norm:
√
kf kD = f T D f (5.166)
5.6. APPROXIMATION IN THE CONFIGURATION SPACE 159
Initialize
valid, ebi , eval , cfgSet, θ ref , θ
Update status
?
valid
yes
Valid model? Approximation
?
yes-
valid f
b
no
Error analysis
?
eval , ebi
no
yes-
ebi , toliw eval
< tol e
Recursive procedure
?
Build leftModel
Build rightModel
?
The validation of the approximation relies on the analysis of the error for a
set of validation congurations θval(i) (i = 1, ..., nval ). The average error of the
validation set is: n val
1 X
eval
= e θ val(i)
(5.167)
nval i=1
The criterion for the validation of the approximation is dened by
where tole is a tolerance. For critical applications, eval can be dened as the
maximal error observed in the validation set.
In order to select the optimal bisection direction, the sensitivity of the relative
error with respect to the conguration parameters θ is analyzed. In particular,
the sensitivity of e with respect to a component θi is dened by
∂e
(5.169)
∂θi
Assuming that the error vanishes at the center of the box, a rst order approx-
imation of the maximal error encountered when moving to the side of the box
along direction i is given by:
∂e θi − θi
ebi =
(5.170)
∂θi 2
Therefore, the optimal bisection direction is dened as the direction which maxi-
mizes ebi .
The sensitivities of the error are computed using a nite dierence method.
The validation set can be exploited for this purpose, for instance, by the de-
nition of validation congurations at mid-distance star-points, as illustrated in
Figure 5.15.
θ2 6
θ2
b
θ2center b × b
b
θ2
-
θ1 θ1center θ1 θ1
Figure 5.15: Validation congurations for error and sensitivity analysis in a two
dimensional conguration space. The validation points are represented by the
circles.
the box into two child boxes. The choice of the bisection direction has a critical
inuence on the eciency of the procedure.
First, it is forbidden to bisect in any direction i for which
5.8 Applications
Two applications are treated in this chapter: an academic exible four-bar
mechanism, and a rigid parallel-kinematic machine-tool, called Orthoglide.
" # " #
M 0 4q̈
M(q) q̈ + ΦTq λ − g(q, q̇, t) = 0
0 0 4λ̈
→ " # " # " #
Kt k ΦTq 4q gext
Φ(q) = 0 + =
k Φq 0 4λ 0
↓
" #
h i ∆θ
∆u = Ψuθ Ψuδ
u = ρ(θ) + Ψuδ (θ) η δ ↔ ηδ
= Ψuη ∆η
↓ ↓
1 δδ
V= 2
η δT K (θ) η δ M = ΨqηT M Ψqη
δδ
→ K = ΨqδT Kt Ψqδ
χ gyr
K= 1
2
η̇ T M(η) η̇ h=h +h
.
↓
Approximation
" #
0
M η̈ + (h.η̇) η̇ + δδ = ΨqηT gext
K ηδ
91 dofs, whereas the reduced model involves 5 modal coordinates (1 rigid mode,
1 constraint mode, 3 internal modes), which are represented for conguration (d)
in Figure 5.18. Two internal modes have an out-of-plane deformation.
For the rigid mode, Figure 5.20 illustrates the variations of the equivalent
mass and gyroscopic tensor in the conguration space. A vertical asymptote
is expected close to the extreme singular congurations. Mathematically, the
components of the rigid modes grow to innity at those points, which explains
the phenomenon. From the actuator point of view, the mechanical blocking at
the singularity is equivalent to an innite inertia. The relative error resulting
5.8. APPLICATIONS 165
Figure 5.18: Mode shapes for conguration (d). From top to bottom: rigid mode,
constraint mode (θ is xed), and 3 internal modes (θ and ze are xed).
166 CHAPTER 5. NONLINEAR MODEL REDUCTION
6 6
Frequencies (Hz)
Frequencies (Hz)
5 5
4 4
3 3
2 2
1 1
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
θ (rad) θ (rad)
700
2500
600 2000
1500
500
1000
400 500
0
300
−500
200
−1000
100 −1500
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
θ (rad) θ (rad)
0.015
0.01 0.1
0.005 0.05
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
θ (rad) θ (rad)
Figure 5.20: Regular grid with tracking strategy: variations of the equivalent
inertia associated with the rigid mode Mθθ (1, 1), and of the gyroscopic tensor
h(1, 1, 1).
6
x 10
3
9000
8000
2.5
7000
Stiffness (N/m)
Stiffness (N/m)
2 6000
5000
1.5
4000
1 3000
2000
0.5
1000
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
θ (rad) θ (rad)
Figure 5.21: Regular grid with tracking strategy: equivalent stiness associated
with the constraint mode. The picture on the right is a zoom.
168 CHAPTER 5. NONLINEAR MODEL REDUCTION
Regular grid
the adaptive grid is presented in Figure 5.22. The adaptive grid consists in 37
boxes (only 5 more than the regular grid), and it is rened close to the singu-
lar congurations, and close to conguration (c). The benets of this strategy
clearly appears in Figure 5.23, where inertia characteristics and relative errors
are compared. For the adaptive strategy, the accuracy is strongly increased near
the singularities, and small errors are tolerated in other parts of the conguration
space, allowing a coarser discretization. Figure 5.24 reveals the ability of this
approach to deal with the non-smooth behavior of the stiness of the constraint
mode around conguration (c). The regular grid leads to a bad approximation
in a large domain around the peak, which may strongly aect the quality of the
model.
Additional information can be recorded in the reduced model. For instance,
the modal contribution of the gravity forces ggrav and the amplitude of the rigid
modes for the constraint dof (Ψgθ ) are plotted in Figure 5.25.
From this example, we conclude that the reduced-order modeling technique
is able to capture accurately the nonlinear changes in the stiness and inertia
properties. The mode-tracking strategy is essential to guarantee the consistency
of the results. Close to the singular congurations, the adaptive conguration
space decomposition is valuable to reduce the approximation error and to op-
timize the computational resources. In order to demonstrate the generality of
5.8. APPLICATIONS 169
700
2500
600 2000
1500
500
1000
400 500
0
300
−500
200
−1000
100 −1500
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
θ (rad) θ (rad)
0.015
regular regular
adaptive adaptive
Rel error − gyr tensor
Rel error − inertia
0.01 0.1
0.005 0.05
0 0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
θ (rad) θ (rad)
6 5
x 10 x 10
15 10
regular regular
adaptive 9 adaptive
exact exact
8
7
Stiffness (N/m)
Stiffness (N/m)
10
6
4
5
3
0 0
0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1
θ (rad) θ (rad)
1000 6
800 4
Gravity torque (N.m)
600
2
400
0
Ψgθ
200
−2
0
−200 −4
−400
−6
−600
−8
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2
θ (rad) θ (rad)
The reduced mass matrix and gyroscopic tensor are constructed. The ap-
proximation in the conguration space is done with quadratic and Lagrange poly-
nomials. The tolerance on the relative error tole is xed to 0.5%, the minimal box
width toliw is xed to 0.05 m (i = 1, 2, 3). The discretization in the conguration
is illustrated in Figure 5.29. Quadratic polynomials require 22 boxes, whereas 8
boxes are sucient for Lagrange polynomials. Therefore, each component of the
reduced model requires the storage of 220 oating-point numbers for quadratic
polynomials, and of 216 numbers for Lagrange polynomials.
The variations of the inertia parameters along a straight line in the con-
172 CHAPTER 5. NONLINEAR MODEL REDUCTION
Figure 5.27: Model of the Orthoglide. The xation of the links on the actuator
slider is modeled as a universal joint to prevent their axial rotation.
guration space are presented in Figure 5.30. The quadratic polynomial is not
appropriate, since signicant discontinuities are observed. The results would be
improved using smaller tolerances tole and toliw , but it would lead to an exagger-
ated discretization of the conguration space. Using Lagrange polynomials, the
relative error is kept below 0.5% for the mass, and below 3% for the gyroscopic
tensor.
As a conclusion, the reduction method is applicable to a spatial mechanism
with 3 kinematic dofs. The nonlinear variations of the inertia properties are sig-
nicant but smooth, since the actuator singularities are avoided. Thus, piecewise
polynomials lead to an ecient representation of the dynamic model.
- the elastic behavior is linear, the reduced stiness only depends on the
kinematic conguration θ, and centrifugal stiening eects are neglected,
5.9. SUMMARY AND CONCLUDING REMARKS 173
Figure 5.28: Extreme congurations of the Orthoglide. In grey lines, the cong-
uration θ = [0 0 0]T m.
174 CHAPTER 5. NONLINEAR MODEL REDUCTION
0.3 0.3
0.25 0.25
0.2 0.2
0.15 0.15
θ2 (m)
θ (m)
0.1 0.1
2
0.05 0.05
0 0
−0.05 −0.05
−0.1 −0.1
−0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3
θ1 (m) θ1 (m)
5 16
Quadratic Quadratic
Lagrange 14 Lagrange
4.5
12
Gyr. tensor (kg.m2)
Inertia (kg.m 2)
4 10
8
3.5
6
3 4
2
2.5
0
2 −2
−0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3
θ1 (m) θ1 (m)
0.035 0.1
Quadratic Quadratic
Lagrange 0.09 Lagrange
0.03
0.08
Rel error − gyr tensor
Rel error − inertia
0.025 0.07
0.06
0.02
0.05
0.015
0.04
0.01 0.03
0.02
0.005
0.01
0 0
−0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3
θ1 (m) θ1 (m)
Figure 5.30: Variations of the inertia components from the extreme conguration
θ = [−0.08 − 0.08 0.26]T m to [0.26 0.26 − 0.08]T m. On the left,
the equivalent mass associated with the rst rigid mode M 11 , on the right the
θθ
component h111 .
5.9. SUMMARY AND CONCLUDING REMARKS 175
- the truncation of the modal basis (this source of error disappears when
considering a rigid mechanism),
The procedure is systematic, and the user only needs to provide high-level
information:
- the tolerance on the approximation error, and the minimum size of the
boxes.
The mode shapes are deduced automatically, even if the mechanical topology is
complex and if the distribution of physical properties in the exible bodies is non-
uniform. This is a denite advantage over assumed-mode techniques, where the
appropriate denition of the modes and boundary conditions is often intricate.
In this chapter, two examples have been successfully analyzed: a exible
four-bar mechanism, and a rigid parallel kinematic machine-tool. In chapter 6,
the reduced-order model of a exible manipulator will be exploited for the design
of a control algorithm.
176 CHAPTER 5. NONLINEAR MODEL REDUCTION
The dofs appearing in the nal model are directly associated with the actua-
tors, the sensors, and other components of the mechatronic system. The method
can be exploited for control design in (at least) three dierent ways:
Experimental Manipulator
Dierent models are elaborated and utilized for frequency domain analysis,
simulation and control design, as detailed in Table 6.1. Those models are dened
either for the mechanism, for the actuated mechanism or for the whole controlled
mechanism. An experimental validation is also conducted.
Guided by those models, a control strategy is developed and implemented for
the exible manipulator. According to a two-time-scale approach, a traditional
joint-tracking controller is complemented with a fast indirect vibration controller.
The vibration controller exploits the inertia forces of the manipulator in order to
damp the exible motion. Hence, it is an extension of inertial damping schemes
developed for the control of macro/micro-manipulators [Geo02].
Actuated Controlled
Mechanism
mechanism mechanism
Full-order
Time simulation
nonlinear
Reduced-order
Control design Time simulation
nonlinear
Reduced-order Root locus
Frequency response
linear Frequency response
Table 6.1: Dierent models and their use for dynamic analysis and control design.
The lines of the tabular are associated with the type of mechanical model, the
columns with the additional components of the mechatronic system.
120 120
110 110
100 100
90 90
α1 (deg)
α2 (deg)
80 80
70 70
60 60
50 50
40 40
30 30
0.7 0.8 0.9 1 1.1 1.2 1.3 0.7 0.8 0.9 1 1.1 1.2 1.3
θ1 (m) θ2 (m)
sure the cylinder extension. Moreover, in order to detect the vibrations of the
mechanism, two accelerometers have been placed at the tip, in two orthogonal
directions.
The control law is implemented on a real-time computer, equipped with a
National Instruments A/D board. The control law is developed in the Simulink
environment, built using Real-Time workshop, and loaded on the real-time com-
puter thanks to the XPC Target architecture.
Ampliers and power supplies are in charge of generating the electrical signals
for the hydraulic valves, as well as to amplify the accelerometer signals.
Ralf
v - - wm
Actuators + Mechanism
θiexcitation
+ vi θim
+ ?
j - Gs
θiref -
i
- Ralf -
−6
Figure 6.5: Independent SISO controller for each actuator (i = 1, 2). Low gain
values G1 = G2 = 100.
−60 −60
−70 −70
−80 0 1
−80 0 1
10 10 10 10
−80 −80
−100
Phase (deg)
Phase (deg)
−100
−120
−120
−140
20 20
0 0
−20 0 1
−20 0 1
10 10 10 10
600 0
Phase (deg)
Phase (deg)
400
−100
200
−200
0 Model Model
Experiment Experiment
−200 0 1
−300 0 1
10 10 10 10
Freq (Hz) Freq (Hz)
θ̇i = Ghydr
i vi , i = 1, 2 (6.2)
6.2. MODEL OF THE MANIPULATOR 187
where θi is the displacement of the actuator i, vi is the applied voltage, and Ghydr
i
is the hydraulic gain, whose value has been tted on the experimental data:
Ghydr
1 = 0.0236, Ghydr
2 = 0.0265 (6.3)
This ideal model does not adequately capture the dynamics of the physical
system in two ways:
- the actuator does not behave as a velocity source near the natural frequen-
cies of the system,
For these reasons, and according to the idea of Obergfell [Obe98], a more elabo-
rated model of the actuator is investigated, accounting for the limited dynamics
of an intrinsic velocity feedback.
In particular, the predicted transfer functions are shifted to the high frequencies,
which means that the model is probably too sti. Therefore, let us consider the
eect of a lumped exibility at the level of the actuators.
188 CHAPTER 6. CONTROL OF AN EXPERIMENTAL MANIPULATOR
−90 0 1
−80 0 1
10 10 10 10
0 −60
−80
Phase (deg)
Phase (deg)
−50
−100
−100
−120
−150 Model Model
−140
Experiment Experiment
−200 0 1
−160 0 1
10 10 10 10
Freq (Hz) Freq (Hz)
Input: v1 ; Output: aY Input: v2 ; Output: aY
60 40
Mag Ratio (dB)
−20 0 1
−20 0 1
10 10 10 10
300 0
200
Phase (deg)
Phase (deg)
−100
100
0
−200
−100 Model Model
Experiment Experiment
−200 0 1
−300 0 1
10 10 10 10
Freq (Hz) Freq (Hz)
Figure 6.8: Actuator model including a lumped exibility. gia denotes the gen-
erated force, kia is the nite stiness of the actuator, θi and Li the actuated dof
and the total length.
6.2. MODEL OF THE MANIPULATOR 189
Conguration 0 Conguration 1
10
0
0
−20 0 1
−10 0 1
10 10 10 10
200 200
Phase (deg)
Phase (deg)
100 100
0 0
20 20
0 0
−20 0 1
−20 0 1
10 10 10 10
200 200
Phase (deg)
Phase (deg)
100 100
0 0
40 40
20 20
0 0
−20 0 1
−20 0 1
10 10 10 10
0 0
Phase (deg)
Phase (deg)
−100 −100
−200 −200
Model Model
Experiment Experiment
−300 0 1
−300 0 1
10 10 10 10
Freq (Hz) Freq (Hz)
Conguration 2 Conguration 3
−10 0 1
−10 0 1
10 10 10 10
200 200
Phase (deg)
Phase (deg)
100 100
0 0
20 40
0 20
−20 0
−40 0 1
−20 0 1
10 10 10 10
200 300
200
Phase (deg)
Phase (deg)
100
100
0
0
−100 Model Model
−100
Experiment Experiment
−200 0 1
−200 0 1
10 10 10 10
Freq (Hz) Freq (Hz)
40 40
20 20
0 0
−20 0 1
−20 0 1
10 10 10 10
0 −50
−100
Phase (deg)
Phase (deg)
−100
−150
−200
−200
Model −250 Model
Experiment Experiment
−300 0 1
−300 0 1
10 10 10 10
Freq (Hz) Freq (Hz)
Conguration 4 Conguration 5
10
10
5
0 0 1
0 0 1
10 10 10 10
200 200
Phase (deg)
Phase (deg)
100 100
0 0
40 20
20 0
0 −20
−20 0 1
−40 0 1
10 10 10 10
200 200
Phase (deg)
Phase (deg)
100 100
0 0
40 40
20 20
0 0
−20 0 1
−20 0 1
10 10 10 10
100 −50
−100
Phase (deg)
Phase (deg)
0
−150
−100
−200
−200 Model Model
−250
Experiment Experiment
−300 0 1
−300 0 1
10 10 10 10
Freq (Hz) Freq (Hz)
m m
Input: v1 ; Output: θ1 Input: v1 ; Output: θ1
−40 −40
Mag Ratio (dB)
−100 0 1
−120 0 1
10 10 10 10
0 −50
−100
Phase (deg)
Phase (deg)
−50
−150
−100
−200
−150 Model Model
−250
Experiment Experiment
−200 0 1
−300 0 1
10 10 10 10
Freq (Hz) Freq (Hz)
m m
Input: v2 ; Output: θ2 Input: v2 ; Output: θ2
−40 −40
Mag Ratio (dB)
−50
−60
−60
−80
−70
−80 0 1
−100 0 1
10 10 10 10
−60 −50
−80
Phase (deg)
Phase (deg)
−100
−100
−150
−120
ga
v - Actuators - Mechanism - wm
Figure 6.13: Mechanical and actuator model. The model of the mechanism in-
cludes the lumped stinesses at the level of the actuators.
194 CHAPTER 6. CONTROL OF AN EXPERIMENTAL MANIPULATOR
Ω2 = diag(ωi2 ) is the diagonal matrix of the square eigenvalues. The mass matrix
and Ω2 depend nonlinearly on the conguration. This dependence is approxi-
mated by piecewise low-order polynomials. The conguration space is decom-
posed into boxes where the local polynomials are dened, in order to achieve
a relative error below the tolerance tole = 0.001. The set of boxes are repre-
sented in Figure 6.15, and the variations of the model parameters are illustrated
in Figure 6.16. The dierent orders of magnitude of the equivalent masses come
from the dierent normalization of the rigid and exible modes, as discussed in
Figure 6.14.
The relation between the end eector position xe and the modal coordinates
is:
xe = ρe (θ) + Ψeδ (θ) η δ (6.8)
Ψeθ is the 2 × 2 matrix of rigid modes at the eector. In this equation, the
eη
nonlinear contribution Ψ̇ (θ) η δ has been neglected under the assumption of
small deformations.
6.2. MODEL OF THE MANIPULATOR 195
Figure 6.14: Mode shapes, home conguration. The amplitudes of the modes
are δθ1 = 0.05 and δθ2 = 0.08, η1δ = 2.8, and η2δ = 1.45. Those dierent orders
of magnitude come from the dierent normalization of the modal cordinates:
the rigid modes are associated with unitary actuator displacements, whereas the
exible modes are normalized with respect to the mass matrix.
196 CHAPTER 6. CONTROL OF AN EXPERIMENTAL MANIPULATOR
0.3
0.2
0.1
θ2 (m)
−0.1
−0.2
6.5
5.5
ω (Hz)
5
4.5
3.5
ω1
3
ω2
2.5
0.7 0.8 0.9 1 1.1 1.2 1.3 1.4
θ1 (m)
9000 90
8000
80
7000
70
6000
M11
M13
5000 60
4000
50
3000
40
2000
1000 30
0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4
θ1 (m) θ1 (m)
Figure 6.16: Variations of the model parameters along a diagonal line in the
conguration space from θ (lower left corner) to θ (upper right corner). The
variations of the natural frequencies are in the interval [3, 7] Hz . We recall
that the modal masses of the exible modes equals 1, as a consequence of their
normalization.
θm
θ ref +
- h - vs - hv - Mechanism
+
Slow ctrl
−6 +
am
6
vf
Fast ctrl
+ j - Gs vi θim
θiref -
i
- Ralf -
−6
Inverse
θm
Mode Inertial
- cδ
η̇
des
θ̈ vf -
actuator
- -
- observer
damping
am model
Figure 6.19: Fast controller. The arrows indicate a model-based scheduling ac-
cording to the conguration changes.
measurements,
- the inertial damping control law denes the desired accelerations of the rigid
des
modes θ̈ that would produce appropriate inertial forces,
In the following, the inertial damping control law is rst presented; the ob-
server and the inverse actuator model will be described afterwards.
6.3. CONTROL DESIGN 201
Inertial damping
In order to analyze the possible damping eect of the coupling inertia forces,
the dynamic equation associated with the exible modes is extracted from the
reduced-order model (6.7):
η̈ δ + Ω2 η δ = −Mδθ θ̈ (6.13)
η̈ δ + Gf η̇ δ + Ω2 η δ = 0 (6.15)
and Gf η̇ δ has the clear interpretation of a damping term, provided that the gain
matrix Gf is positive denite. The stability of the closed-loop system is thereby
guaranteed. It is convenient to dene the gain matrix with respect to the desired
modal damping ratios ξi :
Gf = diag(2 ξi ωi ) (6.16)
As will be motivated in the root locus analysis (section 6.4.2), our design is based
on the choice
ξ1 = ξ2 = 0.1 (6.17)
In order to generate this damping eect, the subsequent problem is to control
the rigid modes so that equation (6.14) is satised. If the number of rigid modes
equals the number of exible modes, and if the inertia coupling matrix Mδθ is
not singular, this equation can be inverted, leading to the denition of the desired
des
acceleration of the rigid modes θ̈ :
des −1
θ̈ = Mδθ Gf η̇ δ (6.18)
−1
The matrix Mδθ depends on the conguration. Recorded in the reduced-
order model, it is represented by piecewise polynomial functions, which can be
−1
exported from the modeling environment to the real-time controller. Mδθ is
thus scheduled online, according to the slow conguration changes θ detected
m
6
x 10
3.5
θ2 = −0.206 m
θ2 = −0.0807 m
3
θ2 = 0.0445 m
θ2 = 0.17 m
2.5
Perf. Index
2
1.5
0.5
0
−0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4
θ1 (m)
This denition is still applicable if the coupling matrix Mδθ is not square. Since
the determinant is the product of the singular values, the higher the performance
index, the farther from inertial singularities. In Figure 6.20, we observe that no
singularity is present in the conguration space of Ralf.
Mode observer
des
In the inertial damping control law (6.18), θ̈ is computed from the modal
velocities η̇ δ , which are not directly available. The observer aims at reconstructing
η̇ δ from the measurements θ m and am . First, let us characterize the relation
between the accelerometer signals and the modal coordinates.
The accelerometers measure the absolute accelerations in body axes, and
those axes rotate during the motion. If R denotes the rotation matrix of the end
6.3. CONTROL DESIGN 203
ẍe = R am (6.20)
we obtain
ẋe ' R cma (6.23)
Inverting this expression, and introducing the modal amplitudes of the end eec-
tor displacements ẋe = Ψeη η̇ , we obtain:
The relation between the measurements and the modal coordinates is sum-
marized by:
" m # " # " #
θ̇ I 0 θ̇
= (6.25)
cma Ψaθ Ψaδ η̇ δ
Since there are as many sensors as modal coordinates, this expression can be
inverted, leading to a static estimation formula for the modal velocities:
" m #
h i θ̇
cδ = aδ −1 aδ −1 (6.26)
η̇ − Ψ Ψaθ Ψ
| {z } cma
=Ψδm
m
θ̇ -
θm Filter d
- -
dt
Ψδm - cδ
η̇
Rt cma
am - Filter -
0 dτ -
Figure 6.21: Observer. The signals of the sensors are rst ltered, then numer-
ically dierentiated/integrated with respect to time, and nally, equation (6.26)
is used to estimate the modal velocity.
In order to avoid any drift in the numerical integration, the integrator is imple-
mented as a high-pass integrator, i.e. an integrator in series with a high-pass
lter.
The construction of the inverse model from the rst actuator model has been
motivated by the direct connection that it denes between the input voltage vi and
the actuator motion θi . On the contrary, the second actuator model involves the
actuator force, so that the input/output relation between vi and θi is aected by
the whole dynamics, in particular, by the mechanical resonances. The inversion
of this model would be more complicated and sensitive to modeling error.
One may be afraid that the velocity source model is not satisfactory near the
resonances, where the fast controller is expected to be performant. However, the
root locus analysis presented in section 6.4.2 will assess the stabilization eect of
the resulting control law.
Figure 6.22: Signal-ow graph of the controlled mechanism: slow and fast control
loops.
imposing that the inertial damping is only active at high-frequencies, near the
resonances of the exible modes η δ .
In order to clarify this point, let us summarize the interactions of both con-
trollers with the exible mechanism. In Figure 6.22, the signal-ow graph illus-
trates the dynamics of the controlled mechanism; in particular, the slow and fast
control loops are clearly represented. Let us briey comment the graph. In the
second actuator model (6.4), the actuator forces ga are inuenced by the applied
voltage v and by the extension rate θ̇. According to the mechanical model (6.7),
these forces excite the rigid dynamics described by the actuator coordinates θ,
which is coupled with the exible dynamics described by the modal amplitudes η δ .
The joint-tracking controller is based on collocated measurements θm , whereas
the fast controller feeds back an estimation of the exible modal rates η̇cδ .
The interaction between both control loops is limited under the time-scale
separation assumption. Indeed, the closed-loop bandwidth of slow controller is
limited to ω s ' 1/3 ω1 , whereas the fast controller is a feedback of the exi-
ble modal coordinates, which are mainly excited near the resonances. At rst
sight, the time-scale separation seems to be satised, and good performances are
expected from the composite controller.
However, we argue that various disturbances may be responsible for unde-
206 CHAPTER 6. CONTROL OF AN EXPERIMENTAL MANIPULATOR
Simulation results
Fast control o Fast control on / η̇ δ = η̇ δ Fast control on / η̇ δ ' η̇ δ
c c
1.2
1.08 1.08
1.06 1.06 1.15
1.04 1.04
1.1
θ2 (m)
θ2 (m)
θ2 (m)
1.02 1.02
1 1 1.05
0.98 0.98
1
0.96 0.96
response response response
0.94 reference 0.94 reference 0.95 reference
0 1 2 3 4 0 1 2 3 4 0 1 2 3 4
t (s) t (s) t (s)
5 5 5
aY (m/s2)
aY (m/s2)
aY (m/s2)
0 0 0
−5 −5 −5
1 2 3 4 1 2 3 4 1 2 3 4
t (s) t (s) t (s)
sirable interactions between the slow and the fast controller. To illustrate this,
Figure 6.23 presents simulation results for a point-to-point trajectory (more de-
tails about the model will be given later). The left plot shows the slow controller
working alone. Good tracking performances are observed for θ2 , but the tip vibra-
tions decay slowly after reaching the nal conguration. In the middle plot, the
fast control law is activated, assuming a perfect estimation of the modal veloci-
ties. The improvement of the tip response is remarkable, and the joint-tracking
performances are barely aected. However, when one accounts for the non-perfect
estimation of η̇ δ due to the approximation in equation (6.23), the fast controller
is still ecient, but the joint-tracking controller is completely disturbed. This
is easily explained, since the approximation (6.23) is only justied in the fast
time-scale. In the slow time-scale, a signicant error is introduced and amplied
by the inverse actuator model, which directly aects the slow controller.
It is possible to improve the observer in order to eliminate this particular
source of error. However, in the experimental system, other disturbances may
aect the observer at low-frequencies, such as the accelerometer errors. Therefore,
there is a need to eliminate the low-frequency components induced in the fast
6.4. CLOSED-LOOP SYSTEM: SIMULATIONS AND EXPERIMENTS 207
Filters design
The fast controller is expected to be active around the rst two natural
frequencies, between 3 Hz and 7 Hz. Both the lower and higher frequency distur-
bances should be attenuated. Two Butterworth band-pass lters are thus consid-
ered, whose properties are illustrated in Figure 6.24. Filter A is a 4th order lter,
with corner frequencies at 1 and 25 Hz , whereas lter B is a 2nd order lter with
corner frequencies at 1.6 and 16 Hz . Both lters have been designed in order
to keep the phase shift below 20 deg in the interval [3, 7] Hz . Filter A has the
advantage of a larger roll-o, leading to an ecient attenuation of low-frequency
disturbances. However, its lower corner frequency is smaller than for lter B, so
that its performances at intermediate frequencies, around ω s = 1.5 Hz , are less
interesting. Both lters will be tested in simulations and experiments.
Two less critical lters can also be mentioned:
- an output high-pass lter associated with the inverse actuator model (4th
order Butterworth lter, with corner frequency at 0.1 Hz ),
- a low-pass lter for the slow controller in order to reduce the measurement
noise (2nd order Butterworh lter, with corner frequency at 100 Hz ).
- a reduced-order linear model is used for the construction of root loci, and
the analysis of frequency responses,
208 CHAPTER 6. CONTROL OF AN EXPERIMENTAL MANIPULATOR
10
Band−pass filter A
100 Band−pass filter B
0
−100
−200
−300 0 1 2
10 10 10
Freq (Hz)
150 150
100 100
Imaginary axis (rad/s)
50 50
0 0
−50 −50
−100 −100
−150 −150
−150 −100 −50 0 −150 −100 −50 0
Real axis (rad/s) Real axis (rad/s)
10 10
8 8
6 6
Imaginary axis (rad/s)
4 4
2 2
0 0
−2 −2
−4 −4
−6 −6
−8 −8
−10 −10
−10 −8 −6 −4 −2 0 −10 −8 −6 −4 −2 0
Real axis (rad/s) Real axis (rad/s)
Figure 6.24: Band-pass lters at the inputs of the fast controller. In the Bode plot,
the closed-loop transfer function of the slow controller is represented, assuming a
velocity source actuator model (actuator 1).
6.4. CLOSED-LOOP SYSTEM: SIMULATIONS AND EXPERIMENTS 209
- Actuators. The actuator model accounts for the intrinsic velocity feedback
dened by equation (6.4).
210 CHAPTER 6. CONTROL OF AN EXPERIMENTAL MANIPULATOR
- Accelerometers. In the Finite Element model, the tip accelerations ẍe are
available in inertial axes1 . In body axes, the accelerometer signals am are
computed by inversion of the formula (6.20) 2 .
- Filters and integrator. The lters are continuous linear state-space systems.
In the model, the exible modes are characterized by a small structural damping
(around 1 %). Obviously, in the closed-loop system, those poles are aected by
the behavior of the control system.
1 According to Remark 4.2 page 75, an "observer equation" is added for the accurate esti-
mation of ẍ .
e
2 Due to rather slow variations of the rotation matrix R with respect to the time-step
(h = 0.01 s), this formula is not strongly nonlinear, and does not represent a problem for
the numerical integration scheme (see Remark 4.4, page 84)
6.4. CLOSED-LOOP SYSTEM: SIMULATIONS AND EXPERIMENTS 211
Root locus
The root loci for increasing gains of the fast controller ξi (i = 1, 2) are
presented in Figure 6.26. Due to interactions with the actuators and the slow
controller, the exible modes have initially a signicant non-negligible damping.
In the absence of band-pass ltering, the vibration controller leads to signicant
damping improvements. Considering the position of the poles for ξi = 0.1, there
is room for further improvements. However, in the presence of ltering, a smaller
active damping can be reached for the rst exible mode. For large gain values,
unstable behaviors are predicted, and the value ξi = 0.1 appears as a reasonable
choice.
Frequency response
The tracking performances of the composite controller can be analyzed in
the transfer functions between the reference inputs and the dynamic responses.
Generally, a joint-tracking controller tries to optimize the responses of the actu-
ator dofs θ, but in practice one is more interested in the positioning of the end
eector. The impact of the fast controller on both types of transfer functions is
illustrated in Figure 6.27.
First, let us consider the predictions of the model when the fast controller is
turned o. In the actuator transfer function, the corner frequency is about 1.5 Hz ,
and the roll-o is -20 dB/decade. 4.4 Hz corresponds to the rst resonance of the
mechanism with clamped actuators. At this frequency, any motion of the actuator
is amplied in the whole structure and requires a huge amount of energy. This
phenomenon appears as an anti-resonance in the collocated transfer function. In
the tip transfer function, the exible modes are excited at the resonance despite
the roll-o of the joint-tracking controller.
When the fast controller is turned on, the actuator transfer function is mod-
ied around the resonance, where the inertial damping is active. The benets
mostly appear in the end eector transfer function, where the resonance peaks
are eciently attenuated. Due to phase lags, the tracking behavior is still not
ideal around the resonance, but the system is now isolated from disturbances
in this frequency range. This clearly illustrates the objective of the composite
control law: performances at low frequencies, and stabilization at the resonances.
212 CHAPTER 6. CONTROL OF AN EXPERIMENTAL MANIPULATOR
Without ltering
Increasing ξ1 Increasing ξ2
50 50
45 45
40 40
35 35
Imaginary axis (rad/s)
With lter A
(the dashed lines still refer to the unltered case)
Increasing ξ Increasing ξ
1 2
50 50
45 45
40 40
35 35
Imaginary axis (rad/s)
30 30
25 25
20 20
15 15
10 10
5 5
0 0
−30 −25 −20 −15 −10 −5 0 5 −30 −25 −20 −15 −10 −5 0 5
Real axis (rad/s) Real axis (rad/s)
Figure 6.26: Root loci. The arrows point to ξi = 0.1. In the second plot, remote
poles and loci associated with the higher corner frequency of the lters exist out
of the axes scales.
6.4. CLOSED-LOOP SYSTEM: SIMULATIONS AND EXPERIMENTS 213
Model
−50 0 1
−40 0 1
10 10 10 10
0 0
−100
Phase (deg)
Phase (deg)
−100
−200
−300
−200
−400 Fast ctrl off Fast ctrl off
Fast ctrl on Fast ctrl on
−500 0 1
−300 0 1
10 10 10 10
Freq (Hz) Freq (Hz)
Experiments
−10
−20
−30 0 1
10 10
0
Phase (deg)
−50
−100
Figure 6.28: Step force applied to the tip of Ralf (home conguration).
Simulation results
Fast control o Fast control on
15 15
10 10
5 5
aY (m/s2)
aY (m/s2)
0 0
−5 −5
−10 −10
−15 −15
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
t (s) t (s)
0.9518 0.9518
0.9516 0.9516
θ2 (m)
θ (m)
0.9514 0.9514
2
0.9512 0.9512
0.951 0.951
Experimental results
Fast control o Fast control on
15 15
10 10
5 5
aY (m/s2)
aY (m/s2)
0 0
−5 −5
−10 −10
−15 −15
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
t (s) t (s)
0.9512 0.9512
0.951 0.951
0.9508 0.9508
θ2 (m)
θ (m)
2
0.9506 0.9506
0.9504 0.9504
0.9502 0.9502
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
t (s) t (s)
Point-to-point motion
A smooth point-to-point trajectory from conguration S4 to conguration
S1 is analyzed. The trajectory θ ref (t) follows bang-bang acceleration proles, i.e.
ref ref
θ̈ is piecewise constant and the velocity θ̇ is limited. At the end conguration
(t = 0.85 s), the transition is not perfectly smooth in order to emulate possible
high-frequency disturbances.
Simulation and experimental results are compared in Figures 6.31. In the
simulation results, three subcases are studied. When the vibration controller is
turned o, the slow controller is able to track eciently the trajectory, with a
small delay. However, tip vibrations continue for a long time after the end of the
trajectory. Those vibrations are dominated by the rst exible mode. Without
ltering, the fast controller leads to an important reduction of those transient
vibrations, but the tracking controller is disturbed in an unacceptable way, as
discussed earlier. Filter A leads to a better trade-o between vibration control
and trajectory tracking; the slow oscillations at the actuator level are dominated
by the poles of the lter.
Qualitatively, the experiments agree with those observations (the non-ltered
218 CHAPTER 6. CONTROL OF AN EXPERIMENTAL MANIPULATOR
Simulation results
Fast control o Fast control on / lter o Fast control on / lter A
1.2
1.08 1.08
1.06 1.15 1.06
1.04 1.04
1.1
θ2 (m)
θ2 (m)
θ2 (m)
1.02 1.02
1 1.05 1
0.98 0.98
1
0.96 0.96
response response response
0.94 reference 0.95 reference 0.94 reference
0 1 2 3 4 0 1 2 3 4 0 1 2 3 4
t (s) t (s) t (s)
30 30 30
20 20 20
10 10 10
aY (m/s2)
aY (m/s2)
aY (m/s )
2
0 0 0
Experimental results
Fast control o Fast control on / lter A
1.08 1.08
1.06 1.06
1.04 1.04
θ2 (m)
θ2 (m)
1.02 1.02
1 1
0.98 0.98
0.96 0.96
response response
0.94 reference 0.94 reference
0 1 2 3 4 0 1 2 3 4
t (s) t (s)
30 30
20 20
10 10
aY (m/s )
aY (m/s2)
2
0 0
−10 −10
−20 −20
−30 −30
0 1 2 3 4 0 1 2 3 4
t (s) t (s)
fast controller has not been tested, fearing unpredictable motions of the mecha-
nism). However, the amplitude level of the transient vibrations after the motion
is higher than the prediction of the simulation, even when the fast controller is
disabled.
Part of this discrepancy may be attributed to unmodeled phenomena and
disturbances. For instance, clearance eects in the joints may excite the exible
motion at the end of the trajectory. In the tip-force experiment, the predictions
of the model were far closer to the experimental results, which can be explained
as follows. The response to a tip force is dominated by the second mode of
vibration, which is almost a pure bending mode of the second link, as represented
in Figure 6.14. Thus, the joints do not participate in the motion, and no clearance
phenomenon is observed. The situation is quite dierent for the point-to-point
motion, since the rst mode of vibration is then excited.
Another explanation may come from the non-perfect actuator model around
the rst natural frequency. Indeed, in Figure 6.12, it was dicult to make a deci-
sion whether the sensor measurements and the actuators were collocated or not.
For our developments, we assumed collocated measurements, but, a posteriori, an
additional simulation has been realized for the non-collocated case. The compar-
ison is given in Figure 6.32. For the noncollocated case, without active vibration
control, the vibrations detected by the linear sensors have a destabilizing eect on
the slow controller. When the vibration controller is activated, the eciency of
the active damping is overestimated by the model. We may conclude that the true
reality is between the collocated and the non-collocated situations. Experimental
data could be exploited for the optimization of an "intermediate" model. In the
following, this identication problem is left aside, and the simulation results rely
on the assumption of collocated measurements.
Figure 6.33 focuses on the behavior of the system at the end of the trajectory,
when the nal conguration has been reached. The damping performances of the
fast controller appear more clearly. Filter B turns out to be more appropriate in
this case: the slow motion of the actuator is more eciently damped.
In order to improve the computational eciency of the simulation algorithm,
the full Finite Element model of the mechanism can be replaced by the nonlinear
reduced model. The reduced model describes the dynamics in terms of modal
coordinates. The tip position is connected with the modal amplitudes by the
220 CHAPTER 6. CONTROL OF AN EXPERIMENTAL MANIPULATOR
20 20
10 10
aY (m/s2)
a (m/s2)
0 0
Y
−10 −10
−20 −20
−30 −30
0 1 2 3 4 0 1 2 3 4
t (s) t (s)
20 20
10 10
aY (m/s2)
a (m/s2)
0 0
Y
−10 −10
−20 −20
−30 −30
0 1 2 3 4 0 1 2 3 4
t (s) t (s)
Experimental results
Fast control o Fast control on / lter A
30 30
20 20
10 10
aY (m/s2)
aY (m/s2)
0 0
−10 −10
−20 −20
−30 −30
0 1 2 3 4 0 1 2 3 4
t (s) t (s)
Simulation results
Fast control o Fast control on / lter A Fast control on / lter B
1.085 1.085 1.085
θ2 (m)
θ2 (m)
1.075 1.075 1.075
aY (m/s2)
aY (m/s2)
0 0 0
−5 −5 −5
1 2 3 4 1 2 3 4 1 2 3 4
t (s) t (s) t (s)
Experimental results
Fast control o Fast control on / lter A Fast control on / lter B
1.085 1.085 1.085
response response response
reference reference reference
1.08 1.08 1.08
θ2 (m)
θ2 (m)
θ2 (m)
20 20 20
10 10 10
aY (m/s2)
aY (m/s2)
aY (m/s )
2
0 0 0
Figure 6.33: Point-to-point trajectory: zoom at the end of the trajectory, after
reaching the end point.
222 CHAPTER 6. CONTROL OF AN EXPERIMENTAL MANIPULATOR
Simulation results
Full model Nonlinear reduced model
1.08
1.06 1.05
1.04
θ (m)
θ2 (m)
1.02
1
2
1
0.98
0.96 response
response 0.95
0.94 reference reference
0 1 2 3 4 0 1 2 3 4
t (s) t (s)
30 30
20 20
10 10
aY (m/s2)
a (m/s2)
0 0
Y
−10 −10
−20 −20
−30 −30
0 1 2 3 4 0 1 2 3 4
t (s) t (s)
1.08 1.08
θ2 (m)
θ2 (m)
1.075 1.075
1.07 1.07
response response
reference reference
1.065 1.065
1 2 3 4 1 2 3 4
t (s) t (s)
5 5
aY (m/s2)
aY (m/s2)
0 0
−5 −5
1 2 3 4 1 2 3 4
t (s) t (s)
Figure 6.34: Comparison of the reduced model and the full model. Point-to-point
trajectory, normal bandwidth of the slow controller, fast controller o.
6.5. CONCLUDING REMARKS 223
Square trajectory
For a point-to-point trajectory, we are especially interested in the perfor-
mances at the nal conguration. Since the performances of the composite con-
troller might depend on the conguration, a square trajectory is analyzed in the
conguration space. The cyclic sequence S1 , S2 , S3 , S4 , S1 was illustrated in Fig-
ure 6.30. The experimental results are presented in Figures 6.35 and 6.36, and the
performances of the vibration controller are appreciable for each conguration.
- a low-order and nonlinear model of the mechanism, for the design of the
vibration controller, and for ecient simulation,
- a linearized model of the controlled mechanism, for root locus and frequency
response analyses,
224 CHAPTER 6. CONTROL OF AN EXPERIMENTAL MANIPULATOR
Experimental results
Control o Control on / lter B
1.1 1.1
response response
1.05 S1 reference 1.05 S1 reference
1 1
S2 S4 S2 S4
θ2 (m)
θ2 (m)
0.95 0.95
0.9 0.9
0.85 0.85
S3 S3
0.8 0.8
0.75 0.75
0 5 10 15 0 5 10 15
t (s) t (s)
30 30
20 20
10 10
aY (m/s )
aY (m/s2)
2
0 0
−10 −10
−20 −20
−30 −30
0 5 10 15 0 5 10 15
t (s) t (s)
Experimental results
Fast control o Fast control on / lter B
30 30
20 20
10 10
aY (m/s2)
aY (m/s )
2
S1
0 0
−10 −10
−20 −20
−30 −30
1 2 3 4 1 2 3 4
t (s) t (s)
30 30
20 20
10 10
aY (m/s2)
aY (m/s )
2
S2
0 0
−10 −10
−20 −20
−30 −30
5 6 7 8 5 6 7 8
t (s) t (s)
30 30
20 20
10 10
aY (m/s2)
a (m/s2)
S3
0 0
Y
−10 −10
−20 −20
−30 −30
9 10 11 12 9 10 11 12
t (s) t (s)
30 30
20 20
10 10
aY (m/s2)
a (m/s2)
S4
0 0
Y
−10 −10
−20 −20
−30 −30
13 14 15 16 13 14 15 16
t (s) t (s)
- a full model of the mechatronic system, for the simulation of the closed-loop
dynamics.
The predictions of those models were generally in good agreement with the exper-
imental results. However, the dynamic parameters of the actuators are critical for
the quality of the model. Their estimation was based on experimental data, but
a more systematic identication procedure would be helpful at this level. Some
discrepancies are also attributed to friction and clearance eects in the joints of
the mechanism.
A two-time-scale strategy has been selected for the design of the controller. A
slow controller is responsible for the tracking performances, and a fast vibration
controller damps the exible motion. The vibration controller is based on the
inertial damping concept, and we conclude that:
- the implementation of this control strategy is fairly simple: only two ac-
celerometers complement the hardware of the joint-tracking controller,
- a low-order model is necessary for the design and the implementation of the
control law, our reduced-order modeling technique is appropriate for this
purpose,
- one diculty is associated with the reliable estimation of the rates of the
exible coordinates, especially at low frequencies,
Conclusion
- The use of the block diagram language for the description of control systems
in the Finite Element environment. This leads to a general, systematic and
230 CHAPTER 7. CONCLUSION
- The local reduction procedure accounts for the nonlinearity of the Global
Modal Parameterization; in particular, the curvature of the parameteriza-
tion contributes to the equivalent inertia forces. Nonlinear eects caused
by large deformations are neglected.
232 CHAPTER 7. CONCLUSION
- The reduction method can be a basis for the extraction of a more structured
model (e.g., a polytopic linear model or a model based on linear fractional
transformation), for which specialized control theories are applicable.
- The inertial damping concept is applied for the vibration control of the
exible manipulator.
Various models are involved in the control design procedure, which conrms
the relevance of our original modeling and simulation tools. From the simulation
and experimental results, we conclude that:
234 CHAPTER 7. CONCLUSION
- The reliable estimation of the modal velocities is critical for the performance
of the controller. Improvements are expected if the acceleration feedback is
replaced by a strain feedback, less sensitive to low-frequency disturbances.
Acknowledgements
This work has been supported by the Belgian National Fund for Scientic Research (FNRS)
which is gratefully acknowledged. It also presents research results of the Belgian Program on
Inter-University Poles of Attraction initiated by the Belgian state, Prime Minister's oce,
Science Policy Programming (AMS IAP V/06). The scientic responsibility is assumed by its
author.
Bibliography
[Ang01] G.Z Angelis. System Analysis, Modelling and Control with Polytopic Linear
Models. PhD thesis, Technische Universiteit Eindhoven, 2001.
[Arn01] M. Arnold. Numerically stable modular time integration of multiphysical sys-
tems. In K.J. Bathe, editor,Proceedings of the First MIT Conference on Com-
putational Fluid and Solid Mechanics, pages 10621064, Cambridge, MA, June
2001. Elsevier, Amsterdam.
[BCGF04] C.L. Bottasso, A. Croce, L. Ghezzi, and P. Faure. On the solution of inverse
dynamics and trajectory optimization problems for multibody systems. Multibody
System Dynamics, 11:122, 2004.
[BCP96] K.E. Brenan, S.L. Campbell, and L.R. Petzold. Numerical Solution of Initial-
Value Problems in Dierential-Algebraic Equations. SIAM, Philadelphia, 2nd
edition, 1996.
[BDG02] O. Brüls, P. Duysinx, and J.-C. Golinval. An adaptation of the Newmark scheme
for the integrated simulation of mechatronic systems. In Proc. of the 6th Int.
236 BIBLIOGRAPHY
[BDG03] O. Brüls, P. Duysinx, and J.-C. Golinval. Computational environment for the
design of exible mechanisms with feedback control. In Proc. of the 6th National
Congress on Theoretical and Applied Mechanics, Ghent, Belgium, May 2003.
[BDG04a] O. Brüls, P. Duysinx, and J.-C. Golinval. Generation of closed-form models for
the control of exible mechanisms: a numerical approach. In Proc. of the 7th Int.
Conf. on Motion and Vibration Control (MOVIC), Saint-Louis, United States,
August 2004.
[BDG04b] O. Brüls, P. Duysinx, and J.-C. Golinval. Reliable simulation of mechatronic sys-
tems using newmark algorithms. In Colloquium of the European Mechanics Soci-
ety on Advances in Simulation Techniques for Applied Dynamics (EUROMECH
452) - Book of Abstracts, Halle, Germany, March 2004.
[BDG04c] O. Brüls, P. Duysinx, and J.-C. Golinval. A systematic model reduction method
for the control of exible multibody systems. In Proc. of the 21st Int. Congress
of Theoretical and Applied Mechanics (ICTAM), Warsaw, Poland, August 2004.
[BDG06] O. Brüls, P. Duysinx, and J.-C. Golinval. A model reduction method for the
control of rigid mechanisms. Multibody System Dynamics, in press:, 2006.
[BG03] O. Brüls and J.-C. Golinval. Simulation of an active control system in a hot-
Virtual Nonlinear
dip galvanizing line. In W. Schiehlen and M. Valásek, editors,
Multibody Systems, volume 103 of NATO Science Series II: Mathematics, Physics
and Chemistry, pages 325332. Kluwer Academic Publishers, 2003.
[Boo84] W.J. Book. Recursive Lagrangian dynamics of exible manipulator arms. Int.
J. Robotics Res., 3(3):87101, 1984.
[Boo93] W.J. Book. Controlled motion in an elastic world. ASME J. Dyn. Syst., Meas.,
and Control, 115(2B):252261, 1993.
[BP87] E. Bayo and B. Paden. On trajectory generation for exible robots. Int. J.
Robotic Systems, 4(2):229235, 1987.
[Car89] A. Cardona. An Integrated Approach to Mechanism Analysis. PhD thesis, Uni-
versité de Liège, 1989.
[CS84] H.R. Cannon and E. Schmitz. Initial experiments on the end-point control of a
exible one-link robot. Int. J. Robotics Res., 3(3):6275, 1984.
[CU02] F.M. Caswara and H. Unbehauen. A neurofuzzy approach to the control of a
exible-link manipulator. IEEE Trans. Robotics and Automation, 18(6):932944,
December 2002.
[CWM02] D. Chablat, Ph. Wenger, and J-P Merlet. Workspace analysis of the orthoglide
using interval analysis. In Proc. of the 8th Int. Symp. on Advances in Robot
Kinematics, Caldes de Malavella, Spain, June 2002. Kluwer Academic Publishers.
[CYC02] J. Cheong, Y. Youm, and W.K. Chung. Joint tracking controller for multi-link
exible robot using disturbance observer and parameter adaptation scheme. Int.
J. Robotic Systems, 19(8):401417, 2002.
[DB94] S. Devasia and E. Bayo. Inverse dynamics of articulated exible structures : Si-
multaneous trajectory tracking and vibration reduction. J. Dynamics and Con-
trol, 4(3):299309, 1994.
[DF00] P. De Fonseca. Simulation and Optimisation of the Dynamic Behaviour of
Mechatronic Systems. PhD thesis, Katholieke Universiteit Leuven, 2000.
[DJB94] J.G. De Jalón and E. Bayo. Kinematic and Dynamic Simulation of Multibody
Systems - The Real-Time Challenge. Springer Verlag, 1994.
[DLS93] A. De Luca and B. Siciliano. Inverse-based nonlinear control of robot arms with
exible links. AIAA J. Guidance, Control, and Dynamics, 16(6):11691176, 1993.
238 BIBLIOGRAPHY
[DV76] B.F. De Veubeke. The dynamics of exible bodies. Int. J. Engineering Science,
14:895913, 1976.
[EBB02] S. Erlicher, L. Bonaventura, and O.S. Bursi. The analysis of the generalized-α
method for non-linear dynamic problems. Computational Mechanics, 28:83104,
2002.
[Eck99] H. Ecker. Comparison 11 (scara-robot) with acsl, full hybrid modelling approach
- environment level. Simulation News Europe, 26:43, July 1999.
[EMO98] H. Elmqvist, S.E. Mattson, and M. Otter. Modelica - the new object-oriented
modeling language. In Proc. of the 12th European Simulation Multiconference
(ESM'98), June 1998.
[ES98] P. Eberhard and W. Schiehlen. Hierarchical modeling in multibody dynamics.
Arch. of Appl. Mech., 68:237246, 1998.
[Fav97] W. Favre. Contribution à la Représentation Bond Graph des Systèmes Mé-
caniques Multicorps. PhD thesis, INSA de Lyon, 1997.
[FLMR99] M. Fliess, J. Lévine, P. Martin, and P. Rouchon. A Lie-Bäcklund approach to
equivalence and atness of nonlinear systems. IEEE Trans. Autom. Control,
44(5):922937, May 1999.
[FLOR95] M. Fliess, J. Lévine, F. Ollivier, and P. Rouchon. Flatness and dynamic feedback
linearizability: Two approaches. In Proc. 3rd European Control Conference, 1995.
[Geo02] L.E. George. Active Vibration Control of a Flexible Base Manipulator. PhD
thesis, Georgia Institute of Technology, 2002.
[GLZ96] S.S. Ge, T.H. Lee, and G. Zhu. Energy-based robust controller design for multi-
link exible robots. Mechatronics, 6(7):779798, 1996.
[GR97] M. Géradin and D. Rixen. Mechanical Vibrations: Theory and Application to
Structural Dynamics. John Wiley & Sons, New York, 2nd edition, 1997.
[GS00] F. Ghorbel and M.W. Spong. Integral manifolds of singularly perturbed systems
with application to rigid-link exible-joint multibody systems. Int. J. of Non-
Linear Mechanics, 35:133155, 2000.
[GVVD04] K. Gallivan, A. Vandendorpe, and P. Van Dooren. Model reduction of MIMO
systems via tangential interpolation. SIAM J. Matrix Analysis and Applications,
26(2):328349, 2004.
[GW84] C.W. Gear and D.R. Wells. Multirate linear multistep methods. BIT, 24:482
502, 1984.
[Hur65] W.C. Hurty. Dynamic analysis of structural systems using component modes.
AIAA J., 3(4):678685, 1965.
[Hus91] R.L. Huston. Computer methods in exible multibody dynamics. Int. J. Num.
Meth. Eng., 32:16571668, 1991.
240 BIBLIOGRAPHY
[HW91] E. Hairer and G. Wanner. Solving Ordinary Dierential Equations II - Sti and
Dierential-Algebraic Problems. Springer-Verlag, 1991.
[Isi95] Isidori. Nonlinear Control System. Springer, 1995.
[JKDW01] L. Jaulin, M. Kieer, O. Didrit, and E. Walter. Applied Interval Analysis.
Springer Verlag, London, 2001.
[JWH00] K.E. Jansen, C.H. Whiting, and G.M. Hulbert. A generalized-a method for
integrating the ltered navier-stokes equations with a stabilized nite element
method. Comp. Meth. Appl. Mech. Eng., 190:305319, 2000.
[Kar97] D. Karnopp. Understanding multibody dynamics using bond graph representa-
tions. J. Franklin Inst., 334B(4):631642, 1997.
[KB94] D.S. Kwon and W.J. Book. A time-domain inverse dynamic tracking control
of a single-link exible manipulator. ASME J. Dyn. Syst., Meas., and Control,
116(2):193200, 1994.
[KK00] T.S. Kim and Y.Y. Kim. MAC-based mode-tracking in structural topology op-
timization. Computers and Structures, 74:375383, 2000.
[KMR00] D. Karnopp, D.L. Margolis, and R.C. Rosenberg. System Dynamics: Modeling
and Simulation of Mechatronic Systems. John Wiley & Sons Inc., 3rd edition,
2000.
[LCG04] E.V. Lens, A. Cardona, and M. Géradin. Energy preserving time integration for
constrained multibody systems. Multibody System Dynamics, 11(1):4161, 2004.
[Lee90] J.W. Lee. Dynamic Analysis and Control of Lightweight Manipulators With
Flexible Parallel Link Mechanism. PhD thesis, Georgia Institute of Technology,
1990.
[LHKB91] H.-J. Lai, E.J. Haug, S.S. Kim, and D.-S. Bae. A decoupled exible-relative
co-ordinate recursive approach for exible multibody dynamics. Int. J. Num.
Meth. Eng., 32:16691690, 1991.
BIBLIOGRAPHY 241
[LWP80] J.Y.S. Luh, M.W. Walker, and R.P.C. Paul. On-line computational scheme for
mechanical manipulators. ASME J. Dyn. Syst., Meas., and Control, 102:6976,
1980.
[MC95] L. Meirovitch and Y. Chen. Trajectory and control optimization for exible space
robots. AIAA J. Guidance, Control, and Dynamics, 18(3):493502, 1995.
[McP96] J.J. McPhee. On the use of linear graph theory in multibody system dynamics.
Nonlinear Dynamics, 9:7390, 1996.
[McP98] J.J. McPhee. Automatic generation of motion equations for planar mechani-
cal systems using the new set of "branch coordinates". Mech. Mach. Theory,
33(6):805823, 1998.
[McP03] J.J McPhee. Virtual prototyping of multibody systems with linear graph theory
and symbolic computing. In W. Schiehlen and M. Valásek, editors, Virtual Non-
linear Multibody Systems, volume 103 of NATO Science Series II : Mathematics,
Physics and Chemistry, pages 3756. Kluwer Academic Publishers, 2003.
[MEFK97] P. Maisser, O. Enge, G. Freundenberg, and G. Kielau. Electromechanical in-
teractions in multibody systems containing electromechanical drives. Multibody
System Dynamics, 1(3):281302, 1997.
[Mei90] L. Meirovitch. Dynamics and Control of Structures. John Wiley & Sons, New-
York, 1990.
[MK91] L. Meirovitch and M.-K. Kwak. Rayleigh-ritz based substructure synthesis for
exible multibody dynamics. AIAA J., 28:17091719, 1991.
[MM95] R.H. Myers and D.C. Montgomery. Response Surface Methodology. Wiley Series
in Probability and Statistics, 1995.
[Mon97] D.C. Montgomery. Design and Analysis of Experiments. John Wiley & Sons,
New York, 4th edition, 1997.
[Nel99] O. Nelles. Nonlinear System Identication with Local Linear Neuro-Fuzzy Models.
PhD thesis, TU Darmstadt, 1999.
[OC89] C.M. Oakley and R.H. Cannon. End-point control of a two-link manipulator with
a very exible forearm: Issues and experiments. In Proc. of the 1989 American
Control Conference, volume 2, pages 13811388, Pittsburgh, PA, USA, June
1989.
[OE97] M. Otter and H. Elmqvist. Energy ow modeling of mechatronic systems via
object diagrams. InProc. of the 2nd Mathmod Vienna, IMACS Symposium on
Mathematical Modelling, pages 705710, 1997.
[OL99] J.H. Oh and J.S. Lee. Control of exible joint robot system by backstepping
design approach. Automation and Soft Computing, 5(4):267278, 1999.
[OS95] B. Owren and H.H. Simonsen. Alternative integration methods for problems in
structural dynamics. Comp. Meth. Appl. Mech. Eng., 122:110, 1995.
[OS01] F. Ollivier and A. Sedoglavic. A generalization of atness to nonlinear systems
of partial dierential equations. application to the command of a exible rod. In
Proceedings of the 5th IFAC Symposium on Nonlinear Control Systems, volume 1,
pages 196200, Saint Petersburg, Russia, July 2001. Elsevier.
[PK99] F.C. Park and J.W. Kim. Singularity analysis of closed kinematic chains. ASME
J. Mech. Des., 121(1):3238, 1999.
[PLS03] J.R. Phillips, D. Luca, and L.M. Silveira. Guaranteed passive balancing trans-
formations for model order reduction. IEEE Trans. Computer-Aided Design of
Integrated Circuits and Systems, 22(8):1027 1041, August 2003.
[Pre97] A. Preumont. Vibration Control of Active Structures, An Introduction. Kluwer,
1997.
[PTVF92] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. Numerical
Recipes in C - The Art of Scientic Computing. Cambridge University Press,
2nd edition, 1992.
[Rai77] M.H. Raibert. Analytical equations vs. table look-up for manipulation: a unifying
concept. In Proc. of the IEEE Conf. Decision and Control, New Orleans, 1977.
BIBLIOGRAPHY 243
[Sco95] M.A. Scott. Time Varying Compensator Design for Recongurable Structures
using Non-Collocated Feedback. PhD thesis, University of Colorado, 1995.
[SF03] J.-C. Samin and P. Fisette. Symbolic Modeling of Multibody Systems. Kluwer
Academic Publishers, 2003.
[SH04] K. Seto and M. Hon. Modeling strategy suited to motion and vibration control
for exible structures. In Proc. of the 7th Int. Conf. on Motion and Vibration
Control (MOVIC), St-Louis, USA, August 2004.
[Sha98] A.A. Shabana. Dynamics of Multibody Systems. Cambridge University Press,
2nd edition, 1998.
[Sim85] J.C. Simo. A nite strain beam formulation. the three-dimensional dynamic
problem. Part I. Comp. Meth. Appl. Mech. Eng., 49:5570, 1985.
[SJK97] R. Sepulchre, M. Jankovic, and P. Kokotovic. Constructive Nonlinear Control.
Springer-Verlag Series in Communications and Control Engineering, 1997.
[SL91] J.-J.E. Slotine and W. Li. Applied Nonlinear Control. Prentice Hall, 1991.
[SLP03] R. Serban, D. Li, and L.R. Petzold. Adaptive algorithms for optimal control
of time-dependent partial dierential-algebraic equation systems. Int. J. Num.
Meth. Eng., 57(2):14571469, 2003.
[SM00] P. Shi and J. McPhee. Dynamics of exible multibody systems using virtual
work and linear graph theory. Multibody System Dynamics, 4:355381, 2000.
244 BIBLIOGRAPHY
[SP01] R. Serban and L.R. Petzold. COOPT - a software package for optimal control of
large-scale dierential-algebraic equation systems. Mathematics and Computers
in Simulation, 56(2):187203, 2001.
[SS90] N.C. Singer and W.P. Seering. Preshaping command inputs to reduce system
vibration. ASME J. Dyn. Syst., Meas., and Control, 112:7682, 1990.
[Str98] O. von Stryk. Optimal control of multibody systems in minimal coordinates. Z.
Angew. Math. Mech, 78:11171120, 1998.
[STW92] J.C. Simo, N. Tarnow, and K.K. Wong. Exact energy-momentum conserving
algorithms and symplectic schemes for nonlinear dynamics. Comp. Meth. Appl.
Mech. Eng., 100:63116, 1992.
[SV01] B. Siciliano and L. Villani. Two-time scale force and position control of exible
manipulators. In Proc. of the 2001 IEEE Int. Conf. on Robotics and Automation,
pages 27292734, Seoul, Korea, 2001.
[SVQ86] J.C. Simo and L. Vu-Quoc. A three-dimensional nite-strain rod model. Part II:
Computational aspects. Comp. Meth. Appl. Mech. Eng., 58:79116, 1986.
[SW83] A.A. Shabana and R.A. Wehage. A coordinate reduction technique for dynamic
analysis of spatial substructures with large angular rotations. J. Struct. Mech.,
11:401431, 1983.
[Sym04] W. Symens. Motion and Vibration Control of Mechatronic Systems with Variable
Conguration and Local Non-Linear Friction. PhD thesis, Katholieke Univer-
siteit Leuven, 2004.
[Til01] M.M. Tiller. Introduction to Physical Modeling with Modelica. Kluwer Academic
Publishers, 2001.
[TOB01] M. Thümmel, M. Otter, and J. Bals. Control of robots with elastic joints based
on automatic generation of inverse dynamics models. In Proc. of the IEEE/RSJ
Int. Conf. on Intelligent Robots and Systems, Maui, Hawaii, USA, October 2001.
[TS85] T. Takagi and M. Sugeno. Fuzzy identication of systems and its applications
to modeling and control. IEEE Trans. Syst., Man, Cybern., 15:116132, 1985.
BIBLIOGRAPHY 245
[VA99] A. Veitl and M. Arnold. Coupled simulation of multibody systems and elastic
structures. In J.A.C. Ambrosio and W.O. Schiehlen, editors, Advances in Com-
putational Multibody Dynamics, pages 635644, IDMEC/IST, Lisbon, Portugal,
September 1999.
[VBSV99] M. Valásek, P. Breedveld, Z. Sika, and T. Vampola. Software tools for mecha-
tronic vehicles: Design through modelling and simulation. Vehicle System Dy-
namics, 33:214230, 1999.
[VGVDS+ 99] A. Veitl, T. Gordon, A. Van De Sand, M. Howell, M. Valásek, O. Vaculin, and
P. Steinbauer. Methodologies for coupling simulation models and codes in mecha-
tronic system analysis and design. Vehicle System Dynamics, 33:231243, 1999.
[VNM98] M.J. Van Nieuwstadt and R.M. Murray. Real-time trajectory generation for
dierentially at systems. Int. J. Robust Nonlinear Control, 18(11):9951020,
1998.
[WC00] P. Wenger and D. Chablat. Kinematic analysis of a new parallel machine tool:
the Orthoglide. In Proc. of the 7th International Symposium on Advances in
Robot Kinematics, Slovenia, June 2000.
[WH82] R.A. Wehage and E.J. Haug. Generalized coordinate partitioning for dimension
reduction in analysis of constrained dynamic systems. ASME J. Mech. Des.,
104:247255, 1982.
[Wil86] T.R. Wilson. The design and construction of exible manipulators. Master's
thesis, Georgia Institute of Technology, 1986.
[YBS89] B.-S. Yuan, W.J. Book, and B. Siciliano. Direct adaptive control of a one-link
exible arm with tracking. Int. J. Robotic Systems, 6(6):663680, 1989.
[ZBG01] D. Zlatanov, I. Bonev, and C. Gosselin. Constraint singularities as conguration
space singularities, 2001. [Link]
[ZT89] O.C. Zienkiewicz and R.L. Taylor. The nite element method. Vol. I. Basic
formulations and linear problems. McGraw-Hill, London, 4th edition, 1989.