Modelling and control of a robotic arm actuated by nonlinear artificial muscles S.N. van den Brink DCT 2007.002

Master of Science Thesis Committee:

Prof. Dr. H. Nijmeijer ‡ Dr. Ir. N. van de Wouw ‡ Dr. Ir. A.G. de Jager ‡

Coaches:

Ir. R.M.G. Rijs z Ir. H.C.A. Dirkx z

‡

Technische Universiteit Eindhoven Department of Mechanical Engineering Dynamics and Control Group

z

Philips Applied Technologies Mechatronics Home robotics project

Eindhoven, 25th January 2007

Abstract Modern industrial robots are designed according to the principal design rules to construct light and stiff. These design rules are not suitable for robots to be used in a domestic environment. Industrial robots are stiff and strong and, for these reasons, dangerous to humans. To allow for robots to be safe in a domestic environment, they need to be compliant. So, in case a domestic robot comes in contact with a human, this does not lead to injuries. Several options are available to make robots compliant. One option is to use compliant actuators, for example a McKibben muscle. This is a pneumatic artificial muscle with great resemblance to a biological muscle. A McKibben muscle consists of a 15 cm piece of bicycle tube, surrounded by a nylon braid of 20 cm. The endings of the tube and the braid are connected using a fitting. Inflating the muscle yields a decrease in muscle length and an increase in diameter and stiffness. To perform research on the applicability of McKibben muscles as robotic actuators, a robotic arm is constructed at Philips Applied Technologies in Eindhoven. It has two degrees of freedom and is actuated by four McKibben muscles. The research goals are, firstly, to construct a predictive model of the robotic arm including McKibben muscles and, secondly, to construct controllers that guarantee accurate motion control of the robotic arm, using model knowledge on the robotic arm and the muscles. A predictive model of the McKibben muscle is proposed, which takes into account the change in muscle shape as well as in stiffness. This model is used in the total robotic arm model, also including the pneumatics and the rigid body dynamics of the arm. Dedicated experiments are performed to identify the parameters of the muscle model and the robotic arm model. Other experiments are used to validate whether the models give a correct prediction of the experimental robotic arm behavior. It appears that the robotic arm behavior is predicted correctly within a limited pressure range. The main model error is in the description of the muscle stiffness. The knowledge gained by the experiments and the modelling is used to define a control strategy. By driving the muscles in pairs, the robotic arm is stabilized completely using linear feedback controllers. Each muscle pair requires two controllers: one to control the robotic arm position and one to control the pressure in both muscles. The stiffness variation as a function of the muscle length of the McKibben muscles is completely compensated for by applying a notch-filter in the position controller. The controllers are implemented on both a one- and two-degree-of-freedom variant of the experimental system. Closed-loop experiments show the effectiveness of the control strategy, where the feasible bandwidth of this controlled system is 4 Hz. To improve the predictive capacities of the robotic arm model, future research on the force balance inside a muscle and the subsequent stiffness is required and advised. i

ii

Samenvatting Moderne industriële robots zijn ontworpen volgens het principe om licht en stijf construeren. Deze ontwerp principes zijn echter niet geschikt voor huiselijke robots. Industriële robots zijn stijf en sterk en zijn daarom gevaarlijk; ze staan om die rede dan ook altijd afgeschermd van mensen. Om de toepassing van robots in een woonomgeving mogelijk te maken, zullen deze compliant (minder stijf) moeten zijn zodat ze meegeven als er tegenaan wordt gestoten. Dit om verwondingen te voorkomen. Er zijn verschillende manieren om een robot compliant te maken. Eén mogelijkheid is om robots aan te drijven met actuatoren die niet stijf zijn, bijvoorbeeld een McKibben spier. Dit is een kunstmatige pneumatische spier die qua gedrag erg veel lijkt op een biologische spier. Een McKibben spier bestaat uit 15 cm fiets binnenband met daarom heen een nylon netje van 20 cm. De uiteinden van de band en het netje zijn aan elkaar bevestigd met een afsluitdop. Door de spier op te blazen wordt deze korter, dikker en stijver. Om de toepasbaarheid van McKibben spieren als robot actuatoren te onderzoeken, is bij Philips Applied Technologies in Eindhoven een robot arm gebouwd. Deze heeft twee graden van vrijheid en wordt aangedreven door vier McKibben spieren. De onderzoeksdoelen zijn, ten eerste, om een voorspellend model van de robot arm inclusief de McKibben spieren te ontwerpen en, ten tweede, om met kennis van de robot arm en van de spieren, regelaars te ontwerpen die gecontroleerde bewegingen van de robot arm garanderen. Een model voor de McKibben spieren is bedacht, dat zowel de vorm verandering als de verandering in stijfheid beschrijft. Dit model is vervolgens gebruikt in het totale robot model, inclusief pneumatiek en de dynamica van de robot arm. Met behulp van experimenten is het gedrag van de spieren en van de robot arm in kaart gebracht. Deze experimenten zijn gebruikt om de parameters in het model af te schatten en om te verifiëren of de modellen een realistische beschrijving geven van de werkelijkheid. Uit validatie experimenten blijkt dat het gedrag van de robot arm in een beperkt druk goed voorspeld wordt. De grootste model fout zit in het voorspellen van de stijfheid van de McKibben spieren. De kennis, verworven door de metingen en het model, is gebruikt om een regelstrategie te formuleren. Het blijkt dat, als de spieren in paren worden aangestuurd, het mogelijk is om met behulp van eenvoudig lineaire regelaars de robot arm volledig te stabiliseren. Voor elk paar spieren zijn twee regelaars nodig; één voor de positie van de robot arm en één voor de druk in beide spieren. Essentieel in het ontwerp van de positie regelaar is een ”notch” filter, die de variatie in stijfheid van de spieren volledig gecompenseerd. Dit is getest op de robot arm voor zowel een configuratie met één

iii

als voor een configuratie met twee graden van vrijheid. Experimenten aan het geregelde systeem tonen de effectiviteit van de voorgestelde regelstrategie aan, waarbij de maximale bandbreedte 4 Hz is. Om het model van de robot arm beter te laten voorspellen, is verder onderzoek naar de krachtbalans in de spieren en de daaruit volgende stijfheid nodig en aanbevolen.

iv

Preface

This report describes the results of my graduation project at the Dynamics and Control group of the Technische Universiteit Eindhoven. The project is carried out under authority of the Mechatronics department of Philips Applied Technologies in Eindhoven. During my assignment I received a lot of help and support by many people. Herewith I would like to show my gratitude. In the first place, I would like to thank my coaches Rob Rijs and Bart Dirkx from Philips Applied Technologies and Nathan van de Wouw and Henk Nijmeijer from Technische Universiteit Eindhoven for their support, input and critical comments during my research. I appreciate that very much. In the second place, but invaluable, I would like to thank the people and students at the Mechatronics department of Philips Applied Technologies for numerous discussions, pleasure and for showing me the true capabilities of engineering. Finally, I like to show my gratitude towards my parents Marcelle and Wim, for their great support during almost ten years of studying, and to my girlfriend Esther, for her patience and comfort.

v

vi

Contents

Abstract

i

Samenvatting

iii

Preface

v

Nomenclature

xi

1

1 4 6

Introduction 1.1 Problem statement 1.2 Report outline

2 Description of the robotic arm and the muscles 2.1 Setup of the robotic arm 2.2 Sensors and control platform 2.3 Pneumatics and the binary valves 2.4 McKibben muscles 2.4.1 Construction 2.4.2 Working principle 2.5 Contribution of the thesis

7 7 8 9 10 10 10 12

3

15 15 18 21 21 24 26 27 27 28 29 30 31 32 33

Physical modelling 3.1 Robotic arm model 3.2 Input signals and the binary valves 3.3 Pneumatics 3.3.1 Pressure and flow relations 3.3.2 Splitter 3.3.3 Formulation of the pneumatic model 3.4 McKibben muscles 3.4.1 Braid kinematics 3.4.2 Muscle volume 3.4.3 Pressure in the muscle 3.4.4 Nominal muscle length 3.4.5 Force exerted by the muscle 3.4.6 Muscle stiffness and damping 3.4.7 Formulation of the muscle model vii

3.5

3.6

Equations of motion of the robotic arm 3.5.1 Robot kinematics 3.5.2 Robot dynamics Discussion

34 34 36 37

4 Identification and model validation 4.1 Muscle identification 4.2 Validation of the equilibria 4.3 Robotic arm identification 4.4 Validation of the linearized model 4.5 Time domain model validation 4.6 Discussion

39 39 45 48 52 55 56

5

59 59 62 63 63 67 67 69 70 71 75 78

Feedback control of the robotic arm 5.1 Single joint control strategy 5.2 Single joint control 5.2.1 Angular controller design 5.2.2 Pressure controller design 5.2.3 Performance 5.3 Dual joint control 5.3.1 Dual joint control strategy 5.3.2 MIMO system identification 5.3.3 Loop shaping of the controller 5.4 Feedback control versus reinforcement learning 5.5 Discussion

6 Conclusions and recommendations

79

A Specifications of the robotic arm parts A.1 Additional photos of the robotic arm A.2 Binary valves A.3 Festo SDE1 pressure sensor A.4 Novotechnik SP2500 potentiometer A.5 Electrical scheme of the switching board

83 83 86 88 89 90

B McKibben muscles properties B.1 Advantages and disadvantages B.2 Literature on McKibben muscles B.2.1 McKibben muscle models B.2.2 Control strategies for systems with PAM’s B.3 Example of a McKibben muscle used as actuator B.4 The tube of a McKibben muscle

91 91 92 92 94 95 97

C Pulse width modulation C.1 Formulation of the PWM algorithm C.2 Saturation function

99 101 103

viii

D Considerations on the pneumatic model D.1 Pneumatic test setup D.2 Flow modelling D.3 Valve modelling D.4 Simulation of the pneumatic test setup D.5 Two DoF pneumatic model

105 105 107 108 109 113

E Considerations on the muscle model E.1 The model of Chou – Hannaford E.2 The model of Petrovi´c E.3 The model for a Shadow McKibben muscle E.4 Analytical solution for the stiffness parameters

117 117 120 120 122

F Derivation of the multibody model F.1 Kinematics of the pneumatic robot arm F.2 Mass and inertial properties F.3 Spring F.4 Muscle forces F.5 Two DoF equation of motion F.6 One DoF equation of motion F.7 Inertial parameters of both bodies F.8 Generalized muscle forces 2 to 6

123 126 128 129 131 132 135 137 139

G Derivation of the static muscle parameters

141

H Tuning of the controllers H.1 Controller part definitions H.2 Frequency response functions of the pressure loop H.3 Nichols diagram of the open-loop pressure control H.4 Tuning of the 1DoF angle controller for joint 1 H.5 Dual joint control scheme H.6 MIMO system identification H.7 Interaction between the joint angles H.8 Stability of the dual joint controller

147 147 148 149 150 153 154 156 157

I

159

Implementation of the robotic arm model

Bibliography

169

ix

x

Nomenclature Variables A A F ~ H J L Q Rs T T V W

cross-sectional area direction cosine matrix force angular momentum inertia half muscle fibre length generalized force specific gas constant temperature time interval volume work

m2 – N kg m2 s−1 kg m2 m Nm kg J s−1 K s m3 J=Nm

b c d ~e f f~ g h k k l m p q q s r t u v x y z

damping stiffness correction parameter diameter unit vector frequency force direction vector gravitational constant muscle length stiffness PWM sample length mass pressure mass flow generalized coordinate open valve time radius time plant input / control signal binary valve signal model state model output model input / flow control signal

N s m−1 – m – Hz – m s−2 m N m−1 – m kg Pa = N m−2 kg s−1 – s m s V V

xi

V

Φ Ψ Ω

power spectrum input or decoupling matrix output matrix

α β β γ ζ η θ λ ξ ̺ ρ ς τ ψ ω ω ~

pneumatic gain flow resistance factor filter parameter braid angle in Chou–Hannaford model muscle damping parameter dynamic viscosity angle of the robotic arm muscle stiffness parameter muscle stiffness parameter density mean density damping ratio valve threshold biased muscle stiffness parameter natural eigen frequency angular velocity vector

C H HOL KI Mθ Mp Pθ Pp R S

controller frequency response function (FRF) open loop FRF interaction coefficient model linearization of θ over zθ model linearization of pA B over zp experimentally identified FRF of θ over zθ experimentally identified FRF of pA B over zp robotic arm sensitivity

Subscripts A B C D AB CD S

muscle A muscle B muscle C muscle D cumulative muscles A and B cumulative muscles C and D spring

c e f in k l m n

capacity error muscle fitting and connector flow in sample low-pass muscle notch xii

– – – m2 – rad Pa s rad Pa kg m−3 kg m−3 – s rad s−1 rad s−1

– – rad V−1 bar V−1 rad V−1 bar V−1 –

out s sv th w 0

flow out splitter sensors and valves threshold working point nominal

Notation C R

set of complex numbers set of real numbers

x ¯ x ˜ x b ~x x X

equilibrium of x perturbation of x estimation of x vector x column x matrix X

Abbreviations AI BW CL CM DoF FRF MIMO ODE OL PAM PID PWM RL SISO

artificial intelligence bandwidth closed-loop center of mass degree-of-freedom frequency response function multiple input – multiple output ordinary differential equation open-loop pneumatic artificial muscle proportional differentiating integrating controller pulse width modulation reinforcement learning single input – single output

xiii

xiv

1 Introduction

This thesis considers the modelling and control of a robotic actuator to be used in a domestic environment. The commonly known robotic actuators are industrial actuators, which are designed for application in industrial robots. In general, industrial robots are unsafe for humans and not practically applicable in a domestic environment. The difference between industrial and domestic robotics is a direct consequence of their specifications. A differentiation in robotics based upon the robot’s application is given in [Kara, 2006] and shown in Figure 1.1. Industrial robots are designed to combine speed and accuracy while domestic robots must be able to be mobile, interactive1 and safe. As a result, industrial robots are stiff, heavy (except for the moving parts) and powerful; domestic robots are small and light but less accurate and less powerful. One of the issues domestic robot research focusses on, is the actuation. Compared to industrial actuators, biologic actuators (muscles) have the ideal properties for domestic robots: strong, flexible, light and efficient. One type of actuators, approaching the biological muscle properties, are pneumatic artificial muscles or PAM’s. The best known PAM is the McKibben muscle2 , a light and strong actuator with nonlinear characteristics. A disadvantage is that McKibben muscles are difficult to control. Modern industrial motion control The design of industrial robots is characterized by lightweight and stiff structures. As a result the eigenmodes are located at high frequencies, which directly affects the feasible bandwidth of a robot. The mechanics are optimized for linearity to allow for the application of well-known classical linear control techniques. In Figure 1.2, two examples of modern industrial robots are given. The ABB IRB is less advanced and general applicable, while the ASML Twinscan is extremely advanced and dedicated for one specific task. 1

2

Interactive in the context of a domestic robot has two meanings. Firstly, domestic robots are characterized by a higher level of autonomy in a loosely structured operating environment; secondly, domestic robots must be able to interact on some level (on/off button up to talking) with people. The first known application of a McKibben muscle is in a orthotic device for a polio patients, developed by J.L. McKibben, see [Nickel, 1963] for details.

Industrial Robots automatically controlled reprogrammable multipurpose manipulator in three or more directions robots for manufacturing

Professional Robots Servicing humans, their environment and equipment maintenance, repair, inspection, surgery

Robotics

Service Robots semi or fully autonomous robot to perform services useful to the well-being of humans all but manufacturing robots

Domestic Robots home care, entertainment, disabled assistance, hobby vacuum cleaners, toys, wheelchair robots

Figure 1.1: A differentiation of robotics, as defined by [Kara, 2006]. Examples of the application area of each group are indicated in italic font.

With both robots, high accuracies and short settling times can be achieved; striking are the accelerations and masses involved. One would probably expect this to be possible due to very complex controllers; in reality straightforward linear feedback controllers are used. The high performance level is achieved by putting effort into three area’s: robot design, modelling the robot behavior and, in case of the ASML Twinscan, knowledge of the robot environment. Modelling the robot dynamics and the internal disturbances enables for a controller design that is optimized for speed and accuracy and, more important, to apply feedforward control. Feedforward tremendously reduces settling times and trajectory tracking capabilities. A dry and clean environment is a minimal requirement for each robot. Often a robot is placed on (active or passive) suspensions to prevent floor vibrations from entering the robot. If accurate knowledge of the environment is required, the robot is placed in a cleanroom. In this way, the system behavior is kept constant; for example dilatation due to changes in temperature, pressure or humidity are prevented. Robots in a domestic environment In the last decades, a vast amount of research regarding domestic robotics has been performed. This research is induced by the expectation that gradually a shift from robotics from an industrial environment to a domestic environment is happening [UNECE, 2004]. The main research focus is on three areas: autonomous behavior, (human) functionality and compliant3 actuators. 3

Compliance has two interpretations. The technical definition means the inverse of stiffness: weakness. In common language, the meaning of compliance refers to accommodating behavior.

2

Introduction

(a) ABB IRB

(b) ASML Twinscan

Figure 1.2: Example of two modern industrial robots. The ABB IRB [ABB, 2006] is a six axes industrial robot, which can be used for a variety of tasks (painting, assembly, material handling, palletizing, etc). Its weight is approximately 1000 kg and its maximum payload is 60 kg; the maximum allowable acceleration is 1.5 g. The IRB does not have special requirements concerning the environment it is used in. The ASML Twinscan [ASML, 2006] is a lithographic system for semiconductor production. The mass of the reticle stage is 20 kg, it is accelerated up to 10 g and positioned with an accuracy of 3 nm with 30 ms settling time. The Twinscan can only be used within a cleanroom.

(a) QRIO

(b) Asimo

(c) Partner robot

Figure 1.3: Three modern humanoids. The QRIO [Sabe, 2005; Fujita, 2005] is a humanoid entertainment robot; it remembers people, learns, shows emotions and is capable of having simple conversations [Tanaka, 2004]. QRIO has 38DoF, is 0.58 m high and weighs 6.5 kg. The Asimo [Honda, 2003] has been introduced in 2000, the height is 1.20 m, it weighs 52 kg, has 26DoF and is capable of walking, running and climbing stairs. The Partner Robot from Toyota, is 1.20 m high and weighs 35 kg, [Toyota, 2004]. It has soft lips which enable for playing trumpet. Also a rolling variant of the partner robot exists.

3

First, ”autonomous” in the given context means that a robot is capable of performing one or several tasks independently; an example is autonomous vacuum cleaners. ”Autonomy” can also be considered in a wider perspective, by taking into account the ability to interact with humans and learning capabilities; this is a vivid research area [Bekey, 2005]. An example of autonomous humanoid is the QRIO as given in Figure 1.3a. Secondly, the functionality research focusses on possible domestic robot layouts. Well known are humanoids, robots with the overall appearance and behavior based on that of the human body. The reason to adopt human features for a robot layout is simple: our world is designed for humans. To have a robot functioning well in our world, it is desirable to resemble us. In Figure 1.3, three humanoid examples are given. These humanoids are built using stiff materials and electric actuators, i.e. according to industrial construction principles. This makes them less suitable for the domestic environment that they are built for. For safety reasons, it is better not to have industry-like robots in a human environment. Imagine a person crossing the path of an industrial robot. Since this robot is unable to detect the person and act accordingly, the robot will move on, conflicting serious injury on the person. To avoid injuries, different design rules apply for robots in a domestic environment: lots of sensors for detecting unsafe situations, less rigid structures with low moving masses and compliant actuators. This brings us to the third domestic robot research area: compliant actuators. In [Bax, 2003], several alternative compliant actuators are considered. One of the actuators rated with the best suitability for application in a domestic environment is a PAM. An extensive overview of several PAM types, is given in [Daerden, 2002]. Two well-known muscle types are the pleated pneumatic muscle [Daerden, 1999] and the McKibben muscle. In [Tondu, 2000] a small historic overview of the McKibben muscles is given. McKibben muscles are the actuators of the robotic arm investigated in this report. In Appendix B.1, the advantages and disadvantages of McKibben muscles are discussed.

1.1 Problem statement The problem statement of this thesis project focusses on the modelling and control of a robotic system, actuated by four pneumatic artificial muscles. To investigate this, a two-degree-of-freedom (2DoF) robotic arm is built at Philips Applied Technologies in Eindhoven; an impression is given in Figure 1.4. The robotic arm is deliberately not designed according to the principal design rules to construct light and stiff [Koster, 2000]; the dynamics of the joints are coupled and the moving masses are relatively large compared to the stiffness. This is done to resemble a possible layout of a domestic robot. The robotic arm consists of two revolute joints which are connected in series. Each joint is driven by a pair of McKibben muscles in an antagonist setup. McKibben muscles are known for their nonlinear behavior. Physical modelling of this behavior is considered to be difficult. Because of such complex dynamics, the control of robots actuated by McKibben muscles, has focussed on various complex control strategies. Classical feedback control has never been applied, simply because a linear approach seemed inadequate to handle the nonlinear characteristics 4

Introduction

of McKibben muscles. This thesis will focus on the classical loop-shaping approach to design a feedback controller. The problem statement of this thesis is formulated as follows: A predictive model for the two DoF robotic arm should be built, based on physical principles and dedicated experimental identification, where special attention should be given to the modelling of the artificial muscles. This model should be validated by experiments and used to construct linear feedback controllers, using a classical feedback approach, allowing for both stabilization and tracking tasks. Finally, these controllers should be implemented and evaluated on the experimental setup. The modelling part of the thesis is started by performing dedicated tests in order to understand the behavior of the system and the McKibben muscles. Based upon these experiments and the resulting increased insight, a model based on physical relations must be formulated. An experimental model must be gained by measuring frequency response functions. Based upon the models and measurement results, a control strategy is chosen and controllers are designed. The control part of the thesis starts with a robotic arm with lower complexity. By fixing the rotation in revolute joint 1, the robotic arm is a DoF system, actuated by muscles A and B ; see Figure 1.4a. If the DoF robotic arm is successfully controlled, it is expanded with the second joint as shown in Figure 1.4b. In [Maas, 2005], the results of a complex control strategy, called reinforcement learning (RL), are presented. This is an artificial intelligence technique, which attempts to stabilize and control the DoF robotic arm by predicting the robotic arm behavior using an experimental model, generated by reinforcement learning. Part of the current thesis is to compare the (black box) RL control approach with the

joint 1

¢ ¢®

¾ body 2

A

¢

body 1

B D

joint 2

­ ­ ­ À

muscle A

°££

¢

C

£

¢¢¸

muscle B (b)

(a)

Figure 1.4: Impression of the DoF and 2DoF robotic arm.

5

spring

(model based) feedback control approach. To be able to make a fair comparison, improving the robotic arm is not in the focus of the assignment.

1.2 Report outline In Chapter 2, the robotic arm and the pneumatic and electronic components are described. As McKibben muscles are not commonly known actuators, its main characteristics and working principles are explained in detail in the second chapter. The physical model for both the DoF as for the 2DoF robotic arm is presented in Chapter 3. All aspects determining the behavior of the robotic arm are modeled individually. To understand the physical principles of every model part, dedicated experiments are performed. These experiments are also used to verify the model parts and to identify the model parameters. The model has multiple purposes. Firstly, it allows for estimating the possible bandwidth and accuracy of the controlled robotic arm on a model level and it can be used as part of the control loop (for example in a feedforward or in an input – output linearization control strategy). Secondly, the influence of every parameter on the behavior can be estimated. In other words, the model can be used for predicting the robots behavior, which can support domestic robot design using artificial muscles as actuators. The model is validated in both the time- and frequency-domain in Chapter 4. The transient behavior and equilibrium points of the experimental robotic arm and the model are compared in the time domain, the linear input – output behavior is compared in the frequency domain. The input – output behavior is determined by measuring the frequency response functions of the robotic arm and by simulation of the model. In Chapter 5, the control approach and results of both the DoF and 2DoF controlled system are presented. The frequency response functions as measured in Chapter 4, are used to determine the lay-out of the controllers. By taking the second DoF into account, the system changes from SISO to MIMO. The complexity increases due to the interaction of the multiple inputs on the outputs. The scaling from a DoF to a 2DoF system, the MIMO system identification, the control approach and the control results are discussed. Finally, the model-based controller is compared with the black-box AI controller. The conclusions and recommendations of this thesis are presented in Chapter 6. Propositions are made for future work considering the motion control of robots in a domestic environment.

6

2 Description of the robotic arm and the muscles

2.1 Setup of the robotic arm The robotic arm consists of two revolute joints. Both joints are driven by two McKibben muscles which work in an antagonist setup. Each muscle is driven by two solenoid binary valves, so there are eight valves needed to drive the two joints; in Figure 2.1a, an image of the robotic arm in an outer position is given, in Figure 2.1b a schematic impression of the pneumatic robot arm in its initial position is given. The location and designation of the four muscles, an additional spring, eight valves, four pressure sensors and the two joints are indicated. In biological systems, the muscles act relatively close to the rotation point. In case of the robotic arm, the McKibben muscles act also relatively close to the rotation points. This makes it difficult to generate a large momentum. On the other hand, small elongations of the muscle are sufficient to generate large rotations of the joints. Body 1 rotates in revolute joint 1, driven by muscles C and D . Its rotation is given by the angle between body 1 and the mounting table surface. Body 1 consists of two bars: an aluminium and a steel part (blue respectively grey in Figure 2.1b). On body 1, muscles A and B , four valves and two pressure sensors are mounted. To support body 1, a spring is mounted between the fixed world (green part in Figure 2.1b) and body 1. This spring functions as a gravity compensator for the mass imbalance of body 1. Without the spring, muscle D is limited in its working area as a relatively large force is required to keep body 1 in position. Body 2 is mounted at the end of body 1, where revolute joint 2 is located. The rotation of body 2 is given by the angle of body 2 with respect to body 1. Body 2 consists of an aluminium bar with a small additional mass on top. Due to the vertical setup of body 2 and its limited mass, static balance can easily be enforced by muscles A and B ; no gravity compensation is required here. In the following section, the functionality of the control platform of the robotic arm is explained. The setup of the pneumatic system is given in Section 2.3 and the properties and behavior of McKibben muscles are explained in Section 2.4. This

joint 1

joint 2

@ @r

¢ ¢ ¢®

body 2

A B

valves

A

valves

r ¡ r r¡r¡ r joint 1 @r ¡ @ pressure @ r¡ sensors r r

B J ] J J

D C

joint 2

air supply

r rD

body 1

spring

A

r

spring Ar

r r table

!!

!!

r !! ! ! switching board

power source (a)

(b)

Figure 2.1: The pneumatic robot arm as available at Philips Applied Technologies. The four muscles are indicated by A , B , C and D . In Figure (a), muscles A and D are completely inflated while muscles muscles B and C are completely stretched. In Figure (b) the equilibrium position is given, in case the pressure in the muscles is equal in each pair, i.e. pA = pB and pC = pD .

chapter is ended by a review on the consulted literature. More information, photographs and specifications of the robotic arm and its components, can be found in Appendix A.

2.2 Sensors and control platform For the control and data acquisition of the robotic arm, a dSPACEr system is used, [dSPACE, 2000]. The controllers are defined in M ATLABr S IMULINKr , [Mathworks, 2005] and uploaded to the dSPACE system. A switching board is built to enable the communication between the dSPACE system and the robotic arm. Moreover, the power supply to the sensors and the valves is embedded on the switching board. The eight valves applied in the robotic arm are controlled by eight binary signals from the dSPACE output. On the switching board, the dSPACE output signals are amplified to 24 V using field emitting transistors or FET’s. Four pressure sensors are used to measure the pressure in the four muscles. The angle of both joints is measured using two linear high-quality potentiometers. Only four analog inputs are available on the dSPACE system, which is less then the six available analog sensor signals. Two jumpers on the switching board enable to select which sensor signals are linked to the dSPACE inputs. The specifications of the pressure sensors and the potentiometers are respectively given in Appendix A.3 and A.4. The electronic scheme of the switching board is given in Appendix A.5. 8

Description of the robotic arm and the muscles

2.3 Pneumatics and the binary valves The pneumatic system, required for driving the McKibben muscles, consists of a number of hoses, a pressure reduction valve, a capacity of 10 liter and the valves. To control the pressure in a muscle, two valves are required: one valve controls the flow towards the muscle, the other valve controls the flow from the muscle into the open. As the valves are placed in series, the flow-in valve is overruled if the flow-out valve is enabled. In conclusion, the flow towards a muscle has three possible states: flow in, flow out and no flow. In Figure 2.2, a schematic overview of the complete pneumatic system is shown. As can be seen, eight 3/2-valves are used; i.e. each valve has three connectors and two possible states. In Appendix A.2, the specifications of the valves are given; most important is the response time: 2 ms. The eight valves are driven by four flow-in signals vin and four flow-out signals vout . The eight binary control signals are given in matrix notation by: · ¸T vin,A vin,B vin,C vin,D v= . (2.1) vout,A vout,B vout,C vout,D The capacity prevents large pressure fluctuations in the air supply due to air losses in the robot pneumatics. The pressure reduction valve ensures the pressure in the capacity to be maximally 1 bar. The original supply pressure is approximately 6 bar. The air hose leaving the capacity is split into four hoses, one hose towards each muscle. reduction valve

vin,A

vout,A

sensor A muscle A

air in air out stop

vin,B

vout,B

sensor B muscle B

capacity

vin,C

vout,C

sensor C muscle C

3/2 valve:

ª v=1 ª v=0

vin,D

vout,D

sensor D muscle D

Figure 2.2: Schematic representation of the complete pneumatic system. The eight valves are depicted for low control signals, i.e. v = 0. The two valves controlling one muscle are placed in series. As a result, the flow-in valve of muscle i is overruled if vout,i = 1.

9

Figure 2.3: A deflated McKibben muscle from the Shadow Robot Company as used in the pneumatic robot arm. Visible are the braid en the two fittings with the air hose connectors.

2.4 McKibben muscles A deflated McKibben muscle as used in the robotic arm is shown in Figure 2.3; it is produced by the Shadow Robot Company [Shadow, 2006]. Besides this muscle type, several other types are available. The first commercially available McKibben muscle was the Rubbertuator from Bridgestone, shown in Figure 2.6a. It is well known for its application in a wall climbing robot [Pack, 1996]. A large number of publications refer to Rubbertuators, unfortunately they are out of production so the presented results can not be verified. Festo produces the pneumatic muscle MAS [Festo, 2006]. It has a different working area than the Shadow muscle, consequently it behaves differently. The pressure required is between 2 and 6 bar and the stroke is smaller; an image is shown in Figure 2.6b. Finally, Merlin Robotics [Merlin, 2006] produces a McKibben muscle which is almost similar to the Shadow muscle.

2.4.1 Construction A McKibben muscle consists of four parts, a rubber tube surrounded by a braid (see Figure 2.4) and two fittings. The fittings serve as gas closures, enable for the mounting and contain the air hose connectors. The braid consists of nylon threads woven into a framework of pantographs. The nylon threads have a fixed length and are considered to be stiff compared to rubber tube, this means that the nylon fibres do not change in length when the muscle is pressurized or a when external force is applied to the muscle. As the threads are woven symmetrically, the thread orientation does not change, i.e. the fittings do not rotate with respect to each other when the muscle is being pressurized or extended.

2.4.2 Working principle In Figure 2.5, the form of a McKibben muscle is shown as a function of the pressure. At 0.0 bar1 , the McKibben muscle is deflated. The rubber tube and the nylon braid do not touch; the muscle length is determined by the rubber tube. Inflating 1

Pressure relative to the environmental pressure

10

Description of the robotic arm and the muscles

Figure 2.4: Structure of a McKibben muscle: a isobutylene (rubber) tube is surrounded by a braid. The weaving of the threads is axially symmetric and the separate threads are not connected to each other. The tube and the braid are fixed together at the ends by two fittings (not shown). Image from [Shadow, 2006].

the muscle, results in an increase in muscle length, while the tube diameter remains constant; this is typical behavior for cylindrical balloons, see [Müller, 2004]. The maximum muscle length is reached when the tube touches the braid completely at the lowest possible pressure, this is at approximately 0.3 bar. If braid and tube touch, the muscle becomes stiff and is able to deliver a force. Further inflation of the muscle results in a decrease in length and an increase in muscle diameter. This is a typical behavior for all PAM types. The shortening of the muscle occurs due to the structure of the braid. Pressurizing the muscle forces the tube to increase in volume. To increase the volume, either the diameter, the length or both must increase. The braids weaving structure forces the diameter to increase. The 200

tren d lin e

185 170

fitting rubber tube

h0 [mm]

nylon braid

0.0

0.3

0.7

1.0

p [bar] Figure 2.5: The shape and length h0 of an unloaded muscle as function of the relative pressure p. Inflating the muscle causes an increase, followed by a decrease in muscle length.

11

(a)

(b)

(c)

(d)

Figure 2.6: In Figures (a) and (b), respectively, a Rubbertuator and a MAS [Festo, 2006] are shown. Figures (c) and (d) show a Shadow muscle at respectively 0.3 and 1.0 bar. Clearly visible is the increase in diameter and decrease in length at 1.0 bar. Compared to the other muscle types, the Shadow muscle is less cylindrically shaped.

braid also determines that an increasing diameter results in a shortening of the pneumatic muscle. The combination of an increasing diameter and a decreasing length, shows that the volume of a muscle hardly changes during inflation. This also means that the additional air due to inflation, mainly results in an increase in pressure. The change in shape for a Shadow muscle is shown in Figures 2.6c and 2.6d at respectively 0.3 and 1.0 bar. The maximum allowable pressure in a Shadow McKibben muscle is 3.5 bar. In case of the robotic arm, it is chosen to resemble biologic muscle behavior; as the pressure determines the stiffness of the muscles, it is decided to use a working pressure of 1.0 bar. The working area of a muscle starts as the tube touches the braid; from this point on the muscle becomes stiff. Normally, the working area starts at 0.3 bar (see Figure 2.5). Applying a pretension to the muscle (by stretching) lowers the pressure at which the working area starts; an additional advantage is that the energy dissipation due to friction in the braid is lowered. In the robotic arm, the muscles are mounted with pretension and, consequently, the working area starts at approximately 0.2 bar. A larger working area of the muscles is useful as this allows for the control of the joint angle of the robotic arm over its entire range.

2.5 Contribution of the thesis An overview of the consulted literature concerning McKibben muscles, is given in Appendix B.2. This overview is divided in two parts: the modelling of McKibben muscle behavior and the control of systems actuated by McKibben muscles. 12

Description of the robotic arm and the muscles

An example of a McKibben muscle used as robot actuator is given in Appendix B.3. In this example, the three types of muscle behavior as discussed in the literature review are illustrated: • The active behavior comprehends the change in muscle length, as a function of the muscle pressure. This type of behavior is particular important in the pressure range below 1.5 bar. Existing models do not consider the active behavior explicitly; as a result, their predictive capacities below 1.5 bar is poor. • The passive behavior describes the stiffness and damping of the muscle. In most cases, the stiffness is derived using the Chou–Hannaford model, [Chou, 1994]. No explicit relations are found on the damping of McKibben muscles and in most some cases, viscous damping is assumed. • The volumetric behavior describes the change in muscle volume. In existing models, the volumetric behavior is based on cylindrical shaped muscles. When considering the Shadow muscles (which are applied in the robotic arm), a barrel shaped muscle is a more likely description. Existing McKibben muscle models do not explicitly make a distinction between these three types of behavior. In this work, an alternative model describing the behavior of McKibben muscles will be proposed. Explicit relations for the three types of behavior will be derived. The purpose of this model is to gain insight in the working principles of a McKibben muscle. The new model must also describe the behavior below 1.5 bar and take into account a different muscle shape. The proposed model will be identified using dedicated experiments and will be subjected to a thorough validation procedure. Most references on control strategies for robots actuated by PAM’s vary from fuzzylogic, gain scheduling, backstepping control, sliding mode control, learning control to feedforward control. These control strategies are rather ”exotic”, compared to linear feedback control. In order to apply PAM’s on a larger scale, a simple but effective control strategy is required. In [Pack, 1994; Schröder, 2003], it is proven that a linear control strategy is possible, although their performance is low. No references are found on a control strategy for robots with PAM’s based on loopshaping. This is probably due to the nonlinear behavior of the muscles: a linear control approach may not seem to represent a good option as stability can not be guaranteed in a non-local sense. Also the required experimental modelling in the frequency domain might be problematic due to nonlinear behavior. How, we conjecture that a loop-shaping approach toward designing high-performance controllers for robots with PAM’s is very well possible. The contribution of this thesis is threefold. Firstly, a model is proposed which gives a better understanding of the observed McKibben muscle behavior. Secondly, it is proven that good closed-loop performance is possible by using linear controllers; and thirdly, a decoupling is proposed which maximizes the application range of the McKibben muscles.

13

14

3 Physical modelling

In this chapter, the modelling of the robotic arm is discussed. The robotic arm model is divided into four parts: the input conversion, the pneumatic system with the valves, the McKibben muscles and the robot dynamics. The purpose of the model is twofold: firstly, understanding the behavior of the McKibben muscles and of the robotic arm; secondly, to support the control design. In Section 3.1, the complete robotic model is presented; in Sections 3.2 to 3.5, the four model parts are elaborated on. This chapter is ended by a discussion.

3.1 Robotic arm model To gain insight in the robotic arm behavior, it is sufficient to consider the DoF robotic arm (with joint 1 fixed). For control purposes, a model describing the 2DoF robotic arm is required. In this chapter, a model describing the behavior of the DoF robotic arm is presented, but the scaling to a 2DoF model is kept in mind. The structure of the robotic arm model is visualized in Figure 3.1; the four model parts are indicated by the circles. The arrows refer to the physical input and output quantities of the respective model parts. The physical quantities are given by four input signals: u=

u

£

uA

input conversion

uB

v

uC

uD

¤T

pneumatics and valves

,

(3.1)

q p

CD muscle A B model

F h

robot dynamics

Figure 3.1: Structure of the robotic arm model. The model consists of four parts. As there are four muscles present in the 2DoF robotic arm, the part describing the muscle behavior is embedded four times; in case of the DoF robotic arm model, the muscle model is implemented two times.

θ

four binary valve-in signals and four binary valve-out signals: · ¸T vin,A vin,B vin,C vin,D v= , vout,A vout,B vout,C vout,D

(3.2)

four flows in and four flows out of the muscles: · ¸T qs,A qs,B qs,C qs,D , q= qout,A qout,B qout,C qout,D

(3.3)

four pressures in the four muscles: £ ¤T p = pA pB pC pD ,

four muscles forces: £ F = FA FB

FC

four muscle lengths: £ h = hA hB hC

FD

hD

(3.4)

¤T ¤T

,

(3.5)

,

(3.6)

and two angles of the robotic arm: £ ¤T . θ = θ1 θ2

(3.7)

The inputs to the 2DoF robotic arm are the four control signals in u. The outputs are given by the two angles in θ and the four pressures p in the muscles. The state of the robotic arm model is completely described by nine states: x =

[ x1

x2

x3

x4

x5

x6

x7

x8

x9 ]T

=

[ pc

pA

pB

θ2

θ˙2

pC

pD

θ1

θ˙1 ]T

,

(3.8)

with θ˙1 and θ˙2 the angular velocities; the pressure in the capacity x1 = pc is defined in the pneumatic model part. If revolute joint 1 of the robotic arm is rigidly fixed, only revolute joint 2 remains; this yields the DoF robotic arm. In the DoF robotic arm, only muscles A and B and joint 2 have a function. When considering Figure 3.1, this means all quantities indexed by C and D and angle θ1 become redundant; i.e. the states of the DoF robotic arm are given by: x =

[ x1

x2

x3

x4

x5 ]T

=

[ pc

pA

pB

θ2

θ˙2 ]T

.

(3.9)

The dynamic model of the DoF robotic arm model is given by the following firstorder differential equation: ¢ ¡   Rs T Vc−1 qin − qs (u) ´ ´ ³ ¡ ¢  ³   Rs T qs,A (u) − qout,A (u) − x2 h˙ A ∂ V + ∂ V ∂ d2 VA−1  ∂h ∂d ∂h   2 A ´ ´ ³  ³  ¡ ¢ −1  . (3.10) ∂ V ∂ d2 ∂V ˙ x˙ =  VB  Rs T qs,B (u) − qout,B (u) − x3 hB ∂h + ∂d2 ∂h  B    x  5   ¡ ¢¡ ¢ 2 −1 QA + QB + m2 g l12 sin(x4 ) Jzz2 + m2 l12 16

Physical modelling

In the first relation of (3.10), Vc , T and Rs define the volume of the capacity, the absolute temperature and the specific gas constant of air respectively. In (3.10), we indicate the dependency of the variables on the input u. Let us now also explicate their dependencies on the states. The two flows are defined by qin = qin (x1 ) and qs = qs (x1 , x2 , x3 , u) respectively. The flows in the second relation are given by qs,A = qs,A (x1 , x2 , x3 , uA , uB ) and qout,A = qout,A (x2 , uA ) and the flows in the third relation are given by qs,B = qs,B (x1 , x2 , x3 , uA , uB ) and qout,B = qout,B (x3 , uB ). All flow relations are defined by the pneumatic model as derived in Section 3.3. In the second and third relation, d2 = d2 (x4 ) represent diameter of the subsequent muscle and VA = VA (x4 ) and VB = VB (x4 ) indicate the muscle volumes. These relations are defined in the muscle model as derived in Section 3.4. The lengths of muscles A and B are given by respectively hA = hA (x4 ) and hB = hB (x4 ) and these are defined by the kinematics of the robotic arm as derived in Section 3.5. In the fifth relation, QA = QA (x2 , x4 ) and QB = QB (x3 , x4 ) represent the generalized muscle forces depending on the pressure and the length of the subsequent muscle, these are derived in Section 3.5. The mass of body 2 is given by m2 and Jzz2 indicates the inertia of body 2 with respect to its center of mass; l12 is the distance of the center of mass of body 2 with respect to joint 2 . The first part of the robotic arm model (Figure 3.1), is the input conversion which converts the input signals u into binary signals v. The binary valves as applied in the robotic arm have two possible states, i.e. open or closed. As a result, the valves behave highly nonlinear. By applying a pulse width modulation (PWM) algorithm to the valves (which is possible due to their fast switching time), the valves approximate linear continuous behavior. For this reason, PWM is added to the robotic arm (and to the model); this is explained in Section 3.2. In Figure 3.2, a new muscle model is proposed that describes the behavior of one muscle i ∈ {A B C D }. In this figure, also the interactions with the pneumatics and the robot dynamics are indicated. The relation defining the pressure pi is considered to be part of the muscle model. From the muscle pressure pi , the stiffness ki , damping di and nominal length h0,i are derived. The muscle pressures p are also an input to the pneumatic model; together with the valve signals v. The pneumatic model is the second part of the robotic arm model; it considers the flows q towards and out of the muscles and the pressure in the capacity pc . The pressure model is discussed in Section 3.3. The third part of the robotic arm model, considers the muscle model as presented in Figure 3.2. The internal variables of muscle i ∈ {A B C D }, are given by the volume Vi , the nominal length h0,i and stiffness ki and damping di of the muscle. In the volumetric behavior, the muscle’s volume is defined; the passive behavior defines the muscle’s stiffness and damping. The active behavior gives the nominal length of the muscle. The nominal length is the length the muscle takes for a given pressure, if no force is applied to the muscle. Remark that the nominal length in general does not equal the actual muscle length hi . The state of a muscle is defined by its pressure pi and actual length hi . The muscle’s pressure dynamics is given in the second and third relation of (3.10), using the flow qi and the muscle’s volume Vi ; in Figure 3.2, this is indicated by the pressure response. The actual muscle length hi is not defined by the muscle model, but is enforced on the muscle by the dynamics of the robotic arm. The muscle model is derived in Section 3.4. 17

m us cl e

m od el

volumetric behavior

pneumatics and valves

qi pi

pressure response

Vi

pi

hi

passive behavior

hi

robot kinematics

hi pi

active behavior

ki

h0,i

θ

di

force response

Fi

robot dynamics

Figure 3.2: The structure of the muscle model and its interaction with the pneumatics and the robot dynamics, see also Figure 3.1. For clarity, the conversion from the joint angles to the muscle lengths using the robot kinematics is indicated; remark that the robot kinematics are also incorporated in the robot dynamics.

The fourth part of the robotic arm model, considers the dynamics of the robotic arm; these are defined by the equations of motion (E O M). The inputs to the E O M are the muscle forces F , which are defined in the force response of the muscle model, as indicated in Figure 3.2. The outputs of the E O M are the revolute joint angles θ. Using the kinematic relations, the actual muscle lengths h are derived from these angles. The E O M of the robotic arm are given in the fourth and fifth relation of (3.10); these are derived in Section 3.5. The robotic arm is presented in this section. The individual parts as indicated in Figures 3.1 and 3.2 are derived in Sections 3.2 to 3.5. This chapter is completed by a brief discussion on the model setup and on the adopted numerical method for simulation of the model, see Section 3.6.

3.2 Input signals and the binary valves The binary valves used in the robotic arm are fast-switching valves. In Appendix A.2, it is shown that the response time of the valves is less then 3 ms. This is very quick compared to other time constants in the robotic arm. For this reason the dynamic behavior of the valves is not taken into account, i.e. the change in valve flow resistance due to a change in valve state is considered to be instantaneous. 18

Physical modelling

The eight binary valves are used to control the four muscle flows. Each valve has two states, either open (v = 1) or closed (v = 0). The pair of valves controlling the flow to one muscle, have only three different states: flow in, flow out and no flow. The state of the valves is determined by the valve control signal v(t), see (2.1). The state of the four valve pairs is determined by the flow control signal: u(t) =

£

uA (t) uB (t)

uC (t)

uD (t)

¤T

.

(3.11)

The flow control signals are continuous and can be positive and negative (u ∈ R4 ). A positive value means flow into the muscle while a negative value means flow out of the muscle. The discrepancy between the continuous signal u(t) and the binary signal v(t) is fixed by applying a pulse width modulation (PWM) algorithm. PWM yields the binary valves to approximate linear continuous behavior (like servo valves) while in reality the binary valves behave highly nonlinear. The PWM algorithm makes use of the fast switching time of the binary valves. Besides, the PWM algorithm is used to convert a single flow signal ui into two valve signals vin,i and vout,i : £

PWM

ui −−−→

vin,i

vout,i

¤

i∈

for

©

A

B

C

D

ª

.

(3.12)

The PWM algorithm yields an output with equal power compared to the input. This is realized by defining a PWM frequency fPWM . Every PWM sample, a pulse is generated of which the duration depends on the sampled value of u(t). First, the sampling time tk and the threshold vth (a saturated1 sawtooth function) are defined: k(t) tk (t)

= ⌊t fPWM ⌋

(3.13a)

k(t) fPWM

(3.13b)

=

uth (t)

= t fPWM − k(t)

(3.13c)

vth (t)

= sat{τ,1−τ } (uth (t)) .

(3.13d)

The saturation τ , as introduced in (3.13d), functions as a dead-zone which takes the limited valve switching time into account. It is of no use to open the valve for a time smaller then its switching time, as this will not result in a flow. Besides, this saturation prevents ”nervous behavior” in an equilibrium point of the controlled robotic arm. Next, the sampled time tk is used to apply a zero-order-hold to u(t). By applying two different saturation functions, u(tk (t)) is divided into a flow in and a flow out signal:

1

uk,in (t)

=

sat{0,1} (u(tk (t)))

(3.14a)

uk,out (t)

=

sat{−1,0} (u(tk (t))) .

(3.14b)

The saturation function y = sat{a,b} (x) is defined as x being saturated by the boundary values a and b. If x ∈ [a, b] the saturation function yields y = x, if x > b, then y = b and if x < a, then y = a. The implementation of sat{a,b} (x) and some examples are given in Appendix C.2.

19

τ

vth (t)

1

ui (t) ui (tk ) 0.5

τ 0 0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

vin,i (t)

0

t [s] Figure 3.3: Example of the PWM algorithm with fPWM = 20 Hz and τ = 0.05 s.

The valve signals v are derived by comparing both sampled signals uk,in (t) and uk,out (t) to the threshold vth (t): v in (t) = uk,in (t) > vth (t)

(3.15a)

v out (t) = uk,out (t) < −vth (t) £ ¤ v in (t) v out (t) . v(t) =

(3.15b) (3.15c)

In Figure 3.3, an example of the PWM algorithm is given. In the upper part of the figure, the threshold vth (t) and the original and the sampled flow control signal ui (t) respectively ui (tk ), are explained graphically. Also the dead-zone introduced by the switching time saturation τ is indicated. In the lower part, the resulting binary valve control signal vin,i is shown. It can be seen that if ui (tk ) > vth (t), this yields vin,i (t) = 1. The values used in this example are fictitious. The complete derivation and implementation of the PWM algorithm is given in Appendix C. Two errors are introduced by the PWM. Firstly, the input signal u(t) is sampled by a zero-order-hold at times tk . This error is linear with fPWM . Secondly, the sampled input signal is saturated by a value τ . Both errors are corrected for in the control loop. Also the flow signal u(t) is saturated between −1 and 1, however this is not an error introduced by the PWM algorithm. Values of u(t) smaller than −1 or larger than 1 will result in an opened valve for the complete duration of the PWM sample. As such, this saturation is part of the pneumatic system and the saturation occurs if the valve signals are oversteering. 20

Physical modelling

3.3 Pneumatics The general setup of the pneumatic system is explained in Section 2.3; in this section a model describing the pneumatic behavior is discussed. To explain certain choices made in the modelling process, the pneumatics of the DoF robotic arm is treated first. Next, this model is extended to the 2DoF robotic arm. No temperature dependency is taken into account in the model, the temperature is considered to be constant and at 20 ◦ C = 293K. The pressure system is modelled as a series of first-order systems. No mass effects of the gas flow (like shockwaves) are taken into account as these process are too quick to play an important role in the behavior of the robotic arm. The complete pneumatic model is based on relative pressures; as a reference the atmospheric pressure patm is taken. This is a rather common approach as most pressure sensors measure relative pressures only. This choice does not affect the model as only pressure differences are important. The gas medium used in the robotic arm is pressurized air; which is considered to be an ideal gas. The equation of state describing an ideal gas is known as Boyle-Gay Lussac’s law: p V = Rs T m ,

(3.16)

with p the pressure, V the volume, Rs the specific gas constant of (in this case) air, T the absolute temperature and m the mass of air in the volume. Air is a compressible gas and, therefore, the density ̺ depends largely on the absolute pressure pabs . This is shown by rewriting (3.16) into: ̺=

m pabs = . V Rs T

(3.17)

3.3.1 Pressure and flow relations In Figure 3.4, a graphical interpretation of the DoF pneumatic model is given, the model lay-out differs from the actual layout, which is given in Figure 2.2. Four pressures and six flows can be distinguished. The out valves are placed on the other side of the muscle, this is done to simplify the implementation of the valves in the model. Another difference between the actual pneumatics and the model is that the airin valves are not shown. In Appendix D, it is shown that modelling the valves as small volumes gives numerical problems in the resulting model. The same holds for the splitter; it can not be modeled as a small volume with multiple flows to it. This is solved by combining the valve signals vin,i and the splitter into one relation defining the splitter pressure ps . The latter is elaborated on in Section 3.3.2. Pressures In Figure 3.5, a volume with an in-going flow qin and out-going flow qout is depicted. The pressure is related to the flow by the time derivative of (3.16): p

dV dp +V dt dt

= Rs T

dm . dt

(3.18) 21

vout,A

splitter vin,A vin,B

qout,A

qs,A pA , V A

pin

qin

qs vout,B

ps

pc , V c

qout,B

qs,B pB , V B

Figure 3.4: Schematic interpretation of the pneumatic model for the 1DoF model. The flow directions are indicated by the arrows. The black solid lines represent hoses and the dotted lines represent binary signals. The following parts can be distinguished: two muscle systems given in yellow and red, the external inputs given in blue, the splitter is orange and the capacity is grey.

p1

p, V qin

p2 q

qout

Figure 3.6: Flow model.

Figure 3.5: Pressure model.

lh 2rt q

2rh

v

q

Figure 3.8: Throttle valve flow.

Figure 3.7: Laminar Poiseuille flow.

p1 pin

qin

p2 q2

q1 pc

air supply

capacity

pm valve

hose 1

hose 2

Figure 3.9: Schematic representation of the pneumatic test setup.

22

muscle

Physical modelling

The variation of mass in the volume depends on the air flowing in and out: dm dt

= qin − qout .

(3.19)

Substituting (3.19) in (3.18) and rearranging terms gives a first-order differential equation for the pressure: µ ¶ dV p˙ = Rs T (qin − qout ) − p V −1 . (3.20) dt From this relation the pressures in the muscles and the capacity will be derived. The capacity has a constant volume of Vc = 10 l; this means the dV dt term in relation (3.20) is equal to zero. In Figure 3.4, the flows towards the capacity are indicated. The pressure in the capacity is now given by: ¡ ¢ p˙c = Rs T Vc−1 qin − qs . (3.21) The pressure in the splitter ps is defined in a different manner as will be explained in Section 3.3.2 and the pressures in the muscles are part of the muscle model. Flows The flow through a hose depends on the pressure difference between the two sides of the hose, as shown in Figure 3.6. Two possible relations are considered in modelling the flow; namely, a laminar Poiseuille flow (Figure 3.7) and the flow through a throttle valve (Figure 3.8). In order to determine the type of flow present in the robotic arm, a pneumatic test setup has been built. The pneumatic test setup is schematically depicted in Figure 3.9. It is used to determine the types of flow through the hoses. Modelling the pneumatic test setup and comparing the simulations with the measurements has shown that each flow through a hose can be modelled by the flow through a throttle valve. This is discussed in more detain in Appendix D. The flow through a throttle valve according to [Polytech, 1997] is given by: p q = β sign (p1 − p2 ) ρ |p1 − p2 | .

(3.22)

The throttle parameter β depends only on the throttle radius rt as indicated in Figure 3.8: √ (3.23) β = 0.59 2 πrt2 .

The density ρ of the air depends on the absolute pressure in the hose as defined by (3.17). As the pressure in the hose is not exactly known, it is estimated by the mean of the two pressures at each side of the hose: ρ

=

p1 + p2 + 2 patm . 2Rs T

(3.24)

When considering the robotic arm, the flows being defined by throttle valves can be understood. For flow qin the throttle parameter βin is defined by the reduction valve. For qs and qout the throttle parameters βs and βout are defined by the binary valves. The flows qs,A and qs,B are defined in Section 3.3.2. 23

It should be remarked that the flows do not depend on the length of the hoses as is the case in a Poiseuille flow. In the current system setup, the flow behavior is dominated by the diameter of the valves only. In case the hoses would become very large, the Poiseuille flow behavior would become dominant. Using Figure 3.4 and (3.22), the flows can be defined by: p = βin sign (pin − pc ) ρ |pin − pc | p = βs sign (pin − pc ) ρ |pc − ps | .

qin qs

(3.25) (3.26)

The flows out of the muscles depend on one pressure only: the muscle pressure. Additionally, the flow out of the muscle depends on the status of the out valves. This is taken into account by: √ = vout,A βout ρ pA √ = vout,B βout ρ pB .

qout,A qout,B

(3.27) (3.28)

3.3.2 Splitter In the DoF robotic arm, the hose leaving the capacity is split into two hoses, one to each muscle, see Figure 3.4. As a result, the flow to each muscle depends largely on the flows towards the other muscle. The most logical way to model this behavior would be to assume one flow to enter and four flows to leave a small volume. The splitting of the flow would be given by a pressure relation as given by (3.20). In Appendix D, it is shown that modelling such small volumes leads to stiff differential equations resulting in numerical problems. Instead of using small volumes (modelled by first-order dynamics), it is decided to combine the binary valves and the splitter into a zeroth-order (algebraic) relation. This means that the pressure in the splitter changes instantaneously upon a change in the valve state. When both in-valves are closed, the pressure in the spitter is equal to the pressure in the capacity. When one of the two valves is opened, the pressure in the splitter is equal the pressure in the opened muscle. If both valves are opened, the pressure in the splitter is equal to the lowest of the muscle pressures (as qs is determined by the largest pressure difference). This behavior is shown in Table 3.1; it can be fitted into a mathematical relation by using vin,A and vin,B , which are defined in (3.15): ps =

£

1 − vin,A

vin,A

¤

·

pc pA

pB min(pA , pB )

¸·

1 − vin,B vin,B

¸

.

(3.29)

To determine the flow qs,i towards each muscle, a proportionality factor which fulfills conservation of mass is required: qs =

X

qs,i .

(3.30)

i

24

Physical modelling

0

1

qs,A qs

0

pc

pA

0

0

1

1

pB

min(pA , pB )

1

0

(pc −pA ) (pc −pA )+(pc −pB )

vin,B

vin,B

vin,A

ps

vin,A 0

1

Table 3.2: The flow from the splitter to muscle A depends on the air-in signals.

Table 3.1: The pressure ps in the splitter depends on the air-in signals.

The proportionality is explained for qs,A using Table 3.2. In case of vin,A = 0, no flow towards muscle A is possible. In case vin,A = 1 while vin,B = 0, the flow qs is determined by the pressure in muscle A according to (3.29) and the flow towards muscle A is given by qs,A = qs . In case of vin,A = vin,B = 1, the muscle with the lowest pressure determines qs according to (3.26) and (3.29). The flow qs is divided into qs,A and qs,B proportional to the pressures pA and pB . The complete relation defining the flow towards muscle A is then given by:

qs,A =

vin,A (pc − pA ) qs . (pc − pA ) + vin,B (pc − pB )

(3.31)

To prevent numerical problems, in the denominator only (pc − pB ) is multiplied by vin,A . Consequently, in case vin,A = vin,B = 0 the denominator will not equal zero. It is assumed that pA will never be equal to pc which is a reasonable assumption as pc increases when both valves are closed. The same reasoning as above applies for muscle B ; i.e. the flow qs,B towards muscle B is given by:

qs,B =

par rt,in rt,s rt,out

vin,B (pc − pB ) qs . vin,A (pc − pA ) + (pc − pB )

value 0.40 mm 0.52 mm 0.72 mm

par βin βs βout

(3.32)

value −7

2

4.194 · 10 m 7.088 · 10−7 m2 1.359 · 10−6 m2

par

value

Rs T Vc

287 J kg−1 K−1 293 K 0.01 m3

Table 3.3: The parameters used in the pneumatic model.

25

3.3.3 Formulation of the pneumatic model The formulation of the complete DoF pneumatic model is given by combining (3.21) – (3.32) into: ¡ ¢ p˙c = Rs T Vc−1 qin − qs (3.33a) ps qin qs

=

£

1 − vin,A

vin,A

= βin sign (pin − pc ) = βs sign (pc − ps )

¤

·

pc pA

pB min(pA , pB )

p ρ |pin − pc |

p

ρ |pc − ps |

¸·

1 − vin,B vin,B

¸

(3.33b) (3.33c) (3.33d)

qs,A

=

vin,A (pc − pA )qs (pc − pA ) + vin,B (pc − pB )

(3.33e)

qs,B

=

vin,B (pc − pB )qs vin,A (pc − pA ) + (pc − pB )

(3.33f)

qout,A

√ = vout,A βout ρ pA

(3.33g)

qout,B

√ = vout,B βout ρ pB

(3.33h)

ρ

=

p1 + p2 + 2patm . 2Rs T

(3.33i)

It might appear that the pressures p˙A and p˙B are forgotten in the above relations. However, as these are part of the muscle model, they are stated in the muscle model in Section 3.4.7. The parameters used in the pneumatic model are defined in Table 3.3. In Appendix D.5, a pneumatic model is proposed for the 2DoF robotic arm. This model is similar to the model as given in (3.33), where two adjustments are required: • the splitter is defined in two steps to take the four muscles into account and

• this requires the flows towards the muscles qs,i to be defined in a slightly different fashion. Remarks In case vin,A = vin,B = 1, it is possible that air flows from the muscle with the highest pressure to the muscle with the lowest pressure. This behavior is taken into account in (3.31) and (3.32), however problems may occur if the pressure in one or both muscles becomes higher than the pressure in the capacity. An increase in muscle pressure is possible while all valves are closed, due to a changed in length of the muscle as will be explained in Section 3.4.2. The relation in (3.29) yields the wrong result for ps in one very specific case: if vin,A = vin,B = 1 and one or both muscle pressures are higher than the pressure in the capacity. This is caused by the min(pA , pB ) term in (3.29), it assumes pc 26

Physical modelling

always to be the largest pressure. Better would be an alternative algorithm for the case vin,A = vin,B = 1: ps ps

= pA = pB

∀ ∀

|pc − pA | > |pc − pB | |pc − pA | < |pc − pB | .

(3.34)

For the sake of simplicity (3.29) is used in the pneumatic model instead of algorithm (3.34). During simulations this did not result in any problems as this very specific case hardly occurs.

3.4 McKibben muscles The muscle model is proposed in Section 3.1, the setup of the model is visualized in Figure 3.2. The physical meaning of the inputs and outputs of each muscle model part, is indicated by the arrows; also the interaction with the pneumatics and the robot dynamics are given. The internal variables in the muscle model are the nominal length h0,i , the muscle volume Vi , the stiffness ki and the damping di ; the state of a muscle is given by the muscle pressure pi and the actual muscle length hi , for i ∈ {A B C D }. An example on the working principle of McKibben muscles and the internal variables is given in Appendix B.3. Let us now revert to the mathematical modelling. First, we construct a relation using the geometric properties of the nylon braid to link the muscle diameter to the muscle length; next the five aspects of the muscle model, as indicated in Figure 3.2, are derived.

3.4.1 Braid kinematics The braid consists of nylon threads with a constant length. These threads are considered infinitely stiff compared to the rubber tube. In Appendix B.4, it is shown that the nylon braid completely determines the shape of a muscle; the only function of the rubber tube is to keep the air from flowing out of the muscle. Therefore, only the braid’s kinematics are taken into account in the muscle model. Three muscle parameters are defined: half the nylon braid length L, the number of windings of a thread n and the diameter of the fitting d1 . Due to the symmetrical weaving, a muscle does not unwind when stretching, i.e. n is constant. Two muscle variables are introduced: the actual muscle length h and the diameter in the middle of the muscle d2 . The muscle kinematics is based on three thought experiments, visualized in Figures 3.10a to 3.10c and discussed below. • For axial stretching, the stretch–force relation is independent of the shape of the braid. Whether the braid is barrel shaped, cylindric or (cut open and) unfolded, the stretch–force relation remains the same. In Figure 3.10b, the cut open and unfolded braid of the muscle in Figure 3.10a is given. • The sides of the unfolded braid are given by the two variable lengths h and c. The second step is to note that length c is equal to c = πnd2 . This may seem odd, as the circumference of the muscle equals πd2 . However, the threads in the unfolded braid are completely stretched. As each thread makes n rotations, length c equals πnd2 . 27

c = πnd2

d2

c

kg

h

kr

h kb

L γ d1 (a)

(c)

(b)

(d)

Figure 3.10: The derivation of the muscle kinematics. The parameters of the Shadow McKibben muscles are given by L = 127.4 mm, n = 1.25 and d1 = 25 mm.

• The threads of the braid makes up a framework of parallelograms. Since the number of parallelograms in the braid is constant, the braid’s behavior is uniform with the behavior of one single parallelogram, this is shown in Figure 3.10c. The muscle’s kinematics can be seen as a bar mechanism consisting of four bars, forming a parallelogram. The length of each bar equals half the nylon thread length L. The kinematic relations of this bar mechanism are given by: p c = 4L2 − h2 (3.35) d2

h

c 1 p 2 4L − h2 = πn πn p 4L2 − (πnd2 )2 . =

=

(3.36) (3.37)

A muscle is completely stretched if its shape is cylindrical, i.e. if d2 = d1 as also shown in Figure 2.5. As a result, the maximum muscle length hmax is defined by the diameter d1 . This yields hmax ≈ 235 mm, which is in accordance with the muscle’s maximum length. It should be noted that the kinematic relations do not directly depend on the pressure in the muscle, i.e. they can be used for both the stretched and the unloaded case.

3.4.2 Muscle volume In Section 2.5, it is stated that the volume of a Shadow McKibben muscle is best described by a barrel-shaped volume. Two barrel formulations are considered, an elliptical and a parabolic shape [Weisstein, 2006]; the elliptical shape is considered to describe the muscle shape best. In Figure 3.11, the three parameters describing an elliptical barrel are given. Parameter d1 is the fitting diameter, the variable muscle length h equals the barrel length and the barrel diameter is equal to d2 . Parameter d2 is a function of the muscle 28

Physical modelling

h

d2

d1

Figure 3.11: Three barrel parameters.

length h(t), as given by (3.36). The volume of an elliptic barrel is given by: V = V (h(t) , d2 (h(t)))

V =

(3.38)

¢ π h¡ 2 3 d1 + 4 d1 d2 (h) + 8 d22 (h) . 60

(3.39)

The actual muscle length h(t) is enforced on the muscle by the robot dynamics. Relations (3.39) and (3.36) define the volumetric behavior block in the muscle model of Figure 3.2.

3.4.3 Pressure in the muscle The pressure in a muscle depends, among other things, on its volume. In Section 3.3, it is shown that the pressure in the muscle satisfies the following differential equation: µ ¶ ¡ ¢ dVi p˙i = Rs T qs,i − qout,i − pi Vi−1 for i ∈ {A B C D } . (3.40) dt

Considering the dependency of the barrel’s volume in (3.38), the pressure in a muscle as defined by (3.40) is rewritten into: µ µ ¶ ¶ ∂V ∂ V ∂ d2 dhi p˙i = Rs T (qs,i − qout,i ) − pi + (3.41) Vi−1 , ∂h ∂d2 ∂h i dt for i ∈ {A B C D }, with ¢ π¡ 2 3 d1 + 4 d1 d2 + 8 d22 60

∂V ∂h

=

∂V ∂d2

=

∂ d2 ∂h

= −

¢ π h¡ d1 + 4 d2 15

(3.42) (3.43)

h √ . π n 4L2 − h2

(3.44)

The inputs to (3.41) are the flows qs,i and qout,i defined in Section 3.3 and the muscle volume as defined by the volumetric behavior. Relation (3.41) defines the pressure response block in the muscle model in Figure 3.2. 29

3.4.4 Nominal muscle length The active behavior of the muscle model defines a relation for the nominal length of the muscle as a function of the pressure p. In Figure 3.2, this is indicated by the active behavior block. To construct such a model, a fourth thought experiment is performed. The nominal length is considered to represent an equilibrium in the face of three forces working in the muscle. These forces are indicated by three springs in Figure 3.10d: • The air in the muscle acts like a gas spring. The gas spring stiffness kg depends on the pressure, an increase in pressure yields an increase in kg (p). As the air pushes from the inside on the rubber tube and braid, it tends to elongate the muscle. • The nylon threads keep the tube from expanding; as a result it opposes elongation of the muscle. Inflating the muscle sharpens the curvature of the muscle (at the endings yielding the orientation of the nylon threads to change what is opposed by the braid; i.e. the braid stiffness kb (p) increases with the muscle pressure. • The material behavior of the rubber tube is viscoelastic; for the sake of simplicity elastic behavior with stiffness kr is assumed, which is considered to be independent of the pressure. The rubber tube opposes elongation. To model the nominal muscle behavior, it is assumed that the braid and the rubber tube always touch. This assumption enables for the following boundary condition. In case the muscle pressure p = 0 bar, the nominal diameter d20 = 0 mm, yielding a nominal muscle length h0 equal to the braid length 2L; i.e. h0 (p = 0) = 2L. A relation determining the nominal length equilibrium is proposed, which fulfils the boundary condition and the three phenomena as described above: h0

∼ = 2L −

kg (p) . kr + kb (p)

(3.45)

Extensive research to expressions for kg (p), kb (p) and kr has unfortunately not resulted in a balanced analytical solution for (3.45). It is decided to adopt (3.45) and assuming kg (p) and kb (p) to be linear with the muscle pressure and kr to be constant. The resulting relation is normalized to two parameters: s1 [m] and s2 [Pa]. The nominal length for muscle i ∈ {A B C D } is now defined by: s1 pi h0,i = 2L − . (3.46) s2 + pi The nominal behavior is experimentally identified by measuring the length and diameter of five different unloaded muscles. Each muscle is mounted on one side and unconstrained on the other side; the length and the diameter are measured for pressures between 0.0 and 1.0 bar. In Figure 3.12, the results of this experiment are given in three plots h0 (p), d20 (p) and d20 (h0 ). Based on heuristics, parameters s1 and s2 of (3.46) are fitted using the experiment shown in the left part of Figure 3.12; this yields s1 = 0.117 m and s2 = 0.232 bar. The model is also plotted in Figure 3.12. The model-based nominal diameter d20 is derived by substituting (3.46) into (3.36). Despite the variations between the several muscles and the model assumptions, the nominal muscle length gives an accurate fit on the experiments. Using the braid kinematics also gives an accurate fit for the nominal diameter. 30

Physical modelling

52

200

50

50

190

48

48

h0 [mm]

180 170

46 44

160

42

150

40

140

0

0.2

0.4

0.6

0.8

p [bar]

1

1.2

38

muscle 1 muscle 2 muscle 3 muscle 4 muscle 5 model 0

0.2

0.4

0.6

0.8

p [bar]

1

1.2

d20 [mm]

52

d20 [mm]

210

46 44 42 40 38 150

160

170

180

190

200

210

h0 [mm]

Figure 3.12: Experimental result for the nominal length h0 and nominal diameter d20 as a function of the pressure p, using five muscles. These measurements are compared to the modelled nominal length as given in (3.46), with s1 = 0.117 m and s2 = 0.232 bar. The model based nominal diameter d20 = d2 (h0 ) is derived using (3.36) with n = 1.25 and L = 0.1274 m. The measurements in the left figure show that the rubber tube touches the braid at approximate 0.3 bar. Also some hysteresis in h0 can be seen; each measurement has two lines, one for increasing pressure (upper line) and one for decreasing pressure.

In Figure 2.5, it can be seen that the assumption, that the braid and the rubber tube always touch, is only true for p > 0.3 bar. As a result, the muscle model is only defined for muscle pressures which exceed 0.3 bar. It is noted that no clear physical interpretation can be given to s1 and s2 . Relation (3.46) defines the active behavior block in the muscle model of Figure 3.2.

3.4.5 Force exerted by the muscle From the robotic arm point of view, a muscle behaves like a spring – damper system. In the dynamics of the robotic arm, this spring – damper behavior is considered to generate a force Fi , which implemented for each muscle i ∈ {A B C D }. The force response of muscle i combines the active and the passive muscle behavior and is given by: ¡ ¢ Fi = k(pi , hi ) hi (t) − h0 (pi ) + d(pi ) h˙ i (t) , (3.47)

with k(pi , hi ) the stiffness and d(pi ) the damping. The stiffness part of the muscle force is given by k(pi , hi ) and the elongation term hi (t)−h0 (pi ). The elongation is defined by the difference between the actual length and the nominal length of the muscle. Herewith, the active muscle behavior (related to the pressure as described in previous section) enters the muscle force relation. The damping part of the muscle force is given by d(pi ) and the time derivative of the muscle length by h˙ i (t). The nominal length or its time derivative plays no direct role in the damping force as only the actual motion determines the damping. This is easily understood; imagine a sudden change in pressure. This causes a change in nominal length and stiffness, resulting in a muscle force. The force results in a motion of the load and a change in muscle length, this motion is subject to dissipation. Relation (3.47) defines the force response block in the muscle model. 31

3.4.6 Muscle stiffness and damping The stiffness k is a function of the pressure pi and the actual muscle length hi . In Appendix E.3, a model describing the stiffness of a Shadow McKibben muscle is derived as: ¡ ¢¡ ¢ k(pi , hi ) = λ + pi ξ0 + ξ1 hi + ξ2 h2i . (3.48) The muscle stiffness behavior is defined by the pressure and the square of the actual muscle length. To compensate for the fact that the stiffness is only defined for p > 0.3 bar, a bias λ on the pressure is introduced. In Appendices E.1 and E.2, the two stiffness models as defined in [Chou, 1996] and [Petrovi´c, 2002a] are given. These two models are used in the derivation of (3.48). Two reason are given to redefine the muscle stiffness model: • The new model takes the muscle volume into account by the barrel description and the braid kinematics; existing models assume a cylindrical volume with fixed diameter. • The new model is adjusted for pressures below 1.5 bar by introducing a bias λ, the existing models are only predictive for pressures exceeding 1.5 bar. No further research on modelling the damping d using other models has been performed in this work. It is assumed that the damping in the muscle is determined by the damping in the rubber tube and by the friction in the nylon braid. The latter increases as the pressure increases. For this reason, a static damping is assumed with an additional pressure-dependent damping. This is given by: d(pi )

= ζ0 + ζ1 pi .

(3.49)

The minor interest in the muscle damping is motivated by the dynamic muscle test as discussed in Section 4.1 and by the work [Klute, 2002]. The order of magnitude of the stiffness and damping is: O(k) O(d)

= 104 = 101 .

(3.50)

The force generated by the muscle has a stiffness and a damping part: F = Fk + Fd . The contribution ratio of the stiffness and damping depends on the frequency ω = θ˙ of the angle θ = sin(θt) of the muscle elongation and is given by: dω cos(θt) dω Fd (ω) = ≈ . Fk k sin(θt) k

(3.51)

Remind that in this simple analysis, a linear damping and stiffness are assumed. In order for the damping to play a significant role, Fd /Fk > 0.1. Since k ≈ 1000 d the damping is only significant if ω > 20 Hz. As the bandwidth of the robotic arm is not expected to exceed 10 Hz, the modelling of the damping is not given extensive attention. Resuming, the contribution of the muscle damping to the behavior of the robotic arm is minimal. Still, some damping must be modeled to ensure energy dissipation in the robotic arm model. The easiest way would be assuming a constant viscous damping of the muscles; however, a pressure dependant viscous damping as given in (3.49) is used. Relations (3.48) and (3.49) define the passive behavior block in the muscle model of Figure 3.2. 32

Physical modelling

3.4.7 Formulation of the muscle model Combining relations in this section, gives the muscle model for i ∈ {A B C D }: q 1 4L2 − h2i (3.52a) d2,i = πn Vi p˙i µ

µ

µ

∂V ∂h

∂V ∂d2 ∂d2 ∂h

¶

i

¶

¢ π hi ¡ 2 3 d1 + 4 d1 d2,i + 8 d2,i 2 (3.52b) 60 ½ µ ¶ ¾ ¡ ¢ ∂V ∂ V ∂ d2 ˙ = Rs T qs,i − qout,i − pi hi Vi−1 + ∂h ∂d2 ∂h i (3.52c) ¡ ¢ π = 3 d21 + 4 d1 d2,i + 8 d2,i 2 (3.52d) 60

=

=

i

¶

i

¢ π hi ¡ d1 + 4 d2,i 15

= −

h0 (pi ) = k(pi , hi ) =

h p i π n 4L2 − h2i

2L − ¡

(3.52e) (3.52f)

s1 pi s2 + pi

λ + pi

¢¡

(3.52g)

ξ0 + ξ1 hi + ξ2 h2i

d(pi ) = ζ0 + ζ1 pi Fi

¢

¡ ¢ = k(pi , hi ) hi − h0 (pi ) + d(pi ) h˙ i .

(3.52h) (3.52i) (3.52j)

The inputs to this model are the two flows qs,i and qout,i (resulting from the pneumatic model as derived in Section 3.3) and the actual muscle length hi (determined by the robot dynamics which will be derived in Section 3.5). The outputs of the muscle model are the pressure pi and the force Fi . In Chapter 4, the passive muscle parameters are identified using an experimental setup. Remarks In Figure 3.12, two phenomena that are not taken into account in the muscle model can be seen. Firstly, the five muscles used in the experiment show a large variability in their nominal behavior, approximately 10 mm variation in their nominal length. The variation between muscles is not modelled. By adjusting muscle specific values for parameters like L and n, the variation might be covered. Secondly, each muscle experiences hysteresis due to friction in the braid and due to viscoelastic rubber behavior. The measurements in the left part of Figure 3.12 contain two curves; the difference between these two curves is caused by the hysteresis. As the hysteresis gives a length variation of approximately 3 mm, it is concluded that the variability between the muscles is dominant over the hysteresis. For this reason, the position error due to hysteresis is not considered either. The energy dissipation due to the hysteresis is considered as an additional viscous damping and taken into account in damping parameter ζ0 . 33

3.5 Equations of motion of the robotic arm In this section, the equations of motion of the robotic arm are derived, using a multibody model. The model of the robotic arm dynamics is shown in Figure 3.13. A muscle is given by a spring k, a damper d and an active part h0 as shown in Figure 3.14. Each muscle is modelled by the force law in (3.47). In the multibody model, the muscles are only considered by means of these forces. In Appendix F, the complete multibody model is derived for both the DoF and 2DoF robotic arm. In this section, the setup of the multibody model is discussed and the DoF multibody model is presented. In the derivation of the equations of motion, the notation as defined in [van de Wouw, 2003] is adopted.

3.5.1 Robot kinematics The dynamic multibody model describes the motion of the rigid bodies in space and time, by taking into account the forces and moments acting on these bodies. Only the x and y dimensions, as indicated in Figure 3.13, are of interest, i.e. the model is in two dimensions. The kinematic relations are required to describe the position of the rigid bodies with respect to each other. In Figure 3.15 (see Figure F.1 for an enlarged version) the kinematic relations of the robotic arm are given. The kinematics are defined in relative coordinates. This means that the position of each mass is defined relative to the position of the mass it is connected to. Three vector frames are introduced. The point O indicates the origin of the fixed frame ~e 0 . In the center of mass of body 1 (CM1), frame ~e 1 is introduced and, in the center of mass of body 2 (CM2), frame ~e 2 is located. Furthermore, three joints are introduced. The point C00 coincides with the origin O. At C01 and C11 = C12, two revolute joints are located, connecting body 1 to the fixed world and body 2 to body 1 respectively. The dimensions of the robotic arm are given in Figure F.2 and Table F.1.

body 2 body 1 A

h0 (p)

B h S

D

k (p, h)

C

d (p)

y x Figure 3.14: Muscle dynamics.

Figure 3.13: Simple robotic arm model.

34

Physical modelling

~e22 CM2 −θ2

~e12 ~b12

A C11 = C12

~e21

~e11 ~b11

CM1 C01

~b01

B

θ1 ~rCM 2 ~rCM 1

~c0

S

C

D ~e20 ~e10 O = C00

Figure 3.15: Kinematic relations.

CM2

−F~A 3 °

F~A

2 ~b2B °

CM1 6 ~b1S ° ~b1D C01

−F~S

F~D

F~S ~b0S

F~B

−F~B

−F~C

−F~D

F~C O = C00 ~b0C

~b0D

Figure 3.16: Robotic arm forces.

35

 °

C11 = C12

~b1A

4 ° ~b1C ~ b 1B 5 °

~b2A

3.5.2 Robot dynamics All forces and their subsequent body-fixed vectors are given in Figure 3.16 (see Figure F.3 for an enlarged version). The spring S connected to body 1 is a purely elastic element; the spring force F~S generated by the spring is a conservative force and is as such implemented in the multibody model. The muscle forces F~i , for i ∈ {A B C D }, describe the complete muscle behavior. Using the principle of virtual work, the muscle forces F~i are implemented as generalized forces Qi into the Lagrange equations of motion [de Kraker, 2001]. The inertial properties of a rotating body in three dimensional are defined by an inertia tensor. In the robotic arm, only two dimensions (x and y) are considered. As a result, only the rotation in the x − y plane is considered; now, the inertial properties of each body are given by one parameter. The inertia of both bodies are estimated numerically by constructing a CAD model of the robotic arm (see Figure 1.4). The inertial properties of the body 1 and 2 with respect to respectively CM1 and CM2 and with respect to respectively C01 and C12, are given in Appendix F.7. Besides these two bodies, a number of other masses are present on the robotic arm, such as the pressure sensors, the valves, the electric wiring and the hoses. The inertia of the hoses and of the muscles are not taken into account. The inertia of the wiring, sensors and valves mounted on body 1 are implemented in one single point mass sv. It should be noted that the stiffness due to the hoses is not considered. The equations of motion for the 2DoF robotic arm are given in Appendix F.5. The equation of motion for the DoF robotic arm is derived in Appendix F.6 and given by: ¡ ¢ 2 + Jzz2 θ¨2 − m2 g l12 sin(θ2 ) = m2 l12 o n¡ ¢ ¡ ¢ FA (l12 − l2A ) l11y − l1A sin(θ2 ) + l11x − lA B cos(θ2 ) + kf~A k

(3.53)

o n¡ ¢ ¡ ¢ FB (l12 − l2B ) l11y + l1B sin(θ2 ) + l11x − lA B cos(θ2 ) . kf~B k

The inputs to (3.53) are the muscle forces FA and FB as defined in (3.47); The forces Fi are in fact the magnitude of the force vector F~i . The direction of the force vector ~ is given by the vector kff~i k , with f~i the vector between the two mounting points i

of muscle i. The gravitational constant is given by g. The numerical values of all parameters as used in (3.53), are given in Table 3.4; the length parameters of the 2DoF multibody model are given in Table F.1 in Appendix F. According to Figure 3.1, the actual muscle length hi is an output of the robot dynamics, which functions as an input to the muscle model. The actual length hi of muscle i is related to the vector f~i by: © ª hi = kf~i k − lfA B for i ∈ A B (3.54) © ª hi = kf~i k − lfCD for i ∈ C D , with lfA B the length of the fittings of muscles A and B and lfCD the length of the fittings of muscles C and D ; the fittings are always in line with the muscle. 36

Physical modelling

par

value

par

value

par

value

l1A l2A l1B l2B l12

76.0 mm 142.0 mm 24.0 mm 210.0 mm 176.0 mm

l11x l11y lA B lfA B g

453.3 mm 26.3 mm 123.3 mm 140.0 mm 9.81 m s−2

msv m1 m2 Jzz1 Jzz2

0.80 kg 2.64 kg 0.533 kg 89.3 g m2 5.70 g m2

Table 3.4: Parameter values of the 1DoF equation of motion.

3.6 Discussion In Section 3.1, the complete robotic arm model is stated and in the following four sections, all relations making up the model are elaborated. It can be concluded that the robotic arm model is rather involved. In this section, some possible improvements to the model are given. In advance of the following chapter, the adopted numerical method for performing simulations with the model is discussed. In the muscle model, the nominal behavior is defined using a force balance between three phenomena. The stiffness of the McKibben muscle is defined by applying the principle of virtual work to the air pushing the rubber tube and the braid and to the elongation of the muscle. Both types of behavior are related to each other as they depend on the the same force balance in case the muscle is unloaded. The relation between these two types of behavior has not become completely clear (as discussed in Section 3.4.4); this requires for further research. To prevent the pneumatic model (and the robotic arm model) to be described by a set of stiff differential equations, the valves and the splitter are defined by analytical relations. Changes in the valve state are considered to be instantaneous. This is a reasonable assumption as the volumes of the splitter, valves and hoses are small compared to the volume of the capacity and the muscles. The adopted relations for the valves and splitter yield a complex interweaving of the valve signals v into the pneumatic model. Disadvantage of the instantaneous pressure change on a changes in valve state is that the model has discontinuous state evolutions. The interweaving of v in the model can be reduced by modelling one hose from the capacity to each muscle; this enables for the valves to be taken into account by the flow control signals u. However, as this also leads to neglecting the division of the flow in the splitter, it must be investigated whether this is a fair assumption. This brings us to the adopted numerical method for solving the robotic arm model. A consequence of the complex interweaving of the binary valve signals v are the discontinuous states. A standard solver, such as ODE45, does not handle these instant changes very well. ODE45 decreases its integration time step if states change rapidly. A dedicated solver for stiff differential equations, such as ODE23s, is no solution either; these solvers still decrease their integration time step locally. The only solution for solving the current model effective is to use a fixed-step ODE solver; these solvers handle the discrete valve signals effectively, yielding in a considerable reduction in computation time. Disadvantage is that no relative and ab37

solute integration error can be defined. The numerical convergence is verified by comparing the solution of ODE2, ODE3, ODE4 and ODE5 with the ODE45 solution; also several integration time steps are considered. based on such comparative evaluation of these numerical schemes and the results of the simulations, it is decided to use a ODE3 solver with a time step of 1 ms. The robotic arm model is implemented in M ATLAB code as well as in a S IMULINK model. The M ATLAB code implementation is given in Appendix I. The S IMULINK implementation and the M ATLAB implementation give identical results. Both make use of an ODE3 solver with an integration time step of 1 ms.

38

4 Identification and model validation

In this chapter, the identification of the McKibben muscles and of the robotic arm is discussed, together with the validation of the robotic arm model. The validation of the robotic arm model is threefold: the equilibria, frequency response functions and time behavior of test setup are compared with the equilibria, linearization and time-domain simulations of the model. The setup of this chapter might seem a little odd. The identification of the muscle model and the validation of the equilibria are discussed in Sections 4.1 and 4.2 respectively; the results of these two sections are closely related. Also the identification of the robotic arm and the validation of the FRF’s are closely related as is discussed in respectively Sections 4.3 and 4.4 The third and last part of the validation is given in Section 4.5, in which the openloop time behavior of the model and the robotic arm are compared. Section 4.6 gives a discussion on the model validation and the applicability of the model.

4.1 Muscle identification The model as proposed in the previous chapter, contains a number of parameters. The parameters of the pneumatic model and the dynamic model are already defined in respectively Tables 3.3 and 3.4. The parameters of the muscle model are yet to be identified; this is the topic of this section. In Section 3.4.6, a model for the stiffness k and damping d of a McKibben muscle is proposed: ¡ ¢¡ ¢ k(p, h) = λ + p ξ0 + ξ1 h + ξ2 h2 . (4.1) d(p)

= ζ0 + ζ1 p .

(4.2)

The stiffness parameters λ, ξ0 , ξ1 and ξ2 and the damping parameters ζ0 and ζ1 are identified using dedicated experiments. For this purpose, five McKibben muscles are mounted into a MTS810 dynamic tensile test bench. Using dynamic testing yields a better identification for the stiffness instead of using merely static tests. Identifying damping always requires dynamic or quasi-dynamic testing. Also possible frequency dependant behavior of the McKibben muscle will be determined in this way. Images of the dynamic test setup are shown in Figures 4.1 to 4.3.

Figure 4.1: The McKibben muscle with the fittings for the MTS tensile bench.The connection between the muscle and the tensile bench is rigidly fixed, so no rotations of the muscle are possible.

force sensors

Q s Q

muscle

MB B

pressure sensor

­ ­ À

manual valve

AU

¾

position actuator

Figure 4.3: Close-up of the McKibben muscle with the pressure and force sensors.

Figure 4.2: MTS810 tensile test bench with McKibben muscle muscle mounted.

40

Identification and model validation

In the dynamic muscle tests, a length is enforced on the muscle. A low amplitude vibration is applied by the testing bench and both the force and the muscle length are measured. The approach for each experiment is given in two steps: • The McKibben muscle is rigidly fixed in the MTS810. A given length is imposed on the muscle (see Figure 4.4) and the pressure in the muscle is set to a level, subsequent to one of the equilibrium gridpoints. The valves are closed in order to keep the air mass inside the muscle constant. • A noise position signal (f ∈ [0, 20] Hz) is superimposed on the reference position. The realized position and the force response of the muscle are measured. This is repeated five times in each working point for every muscle. These tests are performed in the working range with pressures from 0.5 to 1.4 bar and with lengths between 150 to 200 mm. Totally, 40 working points are used for the measurements, the measurement grid is indicated by the green circles in Figure 4.4. In this figure, also the measured nominal lengths of five muscles are shown (the same results as given in Figure 3.12). Muscle 4 is used as default testing muscle, the other four muscles are only used in a few points (indicated by the red dots in Figure 4.4) and serve to verify the muscle 4 measurements. This experiment yields 5 × 40 data sets; the force over position FRF is derived for each measurement. Using a FRF fit, the muscle’s stiffness kw and damping dw are identified in each grid point. This is elaborated on in Appendix G and the results are given in Figure 4.5. The solid lines represent the stiffness and damping for muscle 4, the symbols give the results of the four reference muscles.

h0 muscle 1 h0 muscle 2 h0 muscle 3 h0 muscle 4 h0 muscle 5 gridpoint muscle 4 gridpoint muscle 1 2 3 5

220

210

200

h [mm]

190

180

170

160

150

140 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

p [bar] Figure 4.4: The nominal length h0 of five muscles and the measuring grid for the muscle identification.

41

1.5

34

14

30

dw [Ns m−1 ]

kw [kN m−1 ]

16

12 10 8 6

26 22 18 14

4 2 150

160

170

180

190

h [mm]

10 150

200

170

180

34

dw [Ns m−1 ]

14

160

190

200

h [mm]

muscle 1 muscle 2 muscle 3 muscle 5

16

kw [kN m−1 ]

0.5 bar 0.6 bar 0.7 bar 0.8 bar 0.9 bar 1.0 bar 1.2 bar 1.4 bar

12 10 8 6

150 mm 160 mm 170 mm 180 mm 190 mm 200 mm

30 26 22 18 14

4 2 0.4

0.6

0.8

1

1.2

10 0.4

1.4

0.6

p [bar]

0.8

1

1.2

1.4

p [bar]

20

20

18

18

16

16

14

14

k [kN m−1 ]

k [kN m−1 ]

Figure 4.5: The identified stiffness kw and damping dw as function of the pressure p and the length h of muscle 4. Remark that the colors used in the upper two figures indicate the pressure while in the lower two figures the colors indicate the muscle length. The identified stiffness and damping of muscles 1, 2, 3 and 5 are indicated by a symbol for the eight working points they have been tested in. The color of the symbols refer to the line colors. In case of the stiffness, it appears that the stiffness of muscle 4 is in almost all cases below that of the other muscles. The variation is up to 20%. This figure is similar to Figure G.4.

12 10 8 6

12 10 8 6

4

4

2

2

0 1.7 1.4

0 1.7 1.4 1.1

1.1 0.8

p [bar]

0.8 0.5 0.2

0.12

0.14

0.16

0.18

0.2

0.22

0.24

p [bar]

h [m]

0.5 0.2

0.12

0.14

0.16

0.18

0.2

0.22

h [m]

Figure 4.6: In both figures, the experimentally identified stiffness is given by the solid surface. In the left figure, the analytical solution of the stiffness model is given by the mesh, while in the right figure the mesh represents the fit of the experiment on the model. The stiffness in the analytical model is structurally lower than in the experiments and the fitted model. The experimental identified stiffness is also shown in the left part of Figure G.5; the relative error between both models and the experiment is given in Figure G.6.

42

0.24

Identification and model validation

The stiffness appears to be linear with the pressure and quadratic with the length. On the other hand, the damping does not show a clear dependency on the pressure; a slight increase with the length can be seen. Note also the difference between the five muscles, especially in the damping. The next step is to relate the experimental stiffness to the stiffness model as defined in (4.1). Due to the multiplication of h and p in (4.1), it is not possible to use a least squares fit to identify the four muscle parameters directly on the experimental data kw . This is solved by using six parameters for the least squares fit, [Kreyszig, 1999]. These six parameters are derived from the original four parameters according to: ¡ ¢¡ ¢ k = λ + p ξ0 + ξ1 h + ξ2 h2 (4.3) ¢ ¡ (4.4) k = λ ξ0 + λ ξ1 h + λ ξ2 h2 + ξ0 p + ξ1 p h + ξ2 p h2 . 1

2

3

4

5

6

Unfortunately, the least squares fit yields no exact solution for the four parameters λ, ξ0 , ξ1 and ξ2 ; i.e. the six parameters in (4.4) cannot be written into four parameters. By adding two corrections parameters c0 and c1 , the model is improved: ¡ ¢ k = λ ξ0 c0 + λ ξ1 c1 h + λ ξ2 h2 + ξ0 p + ξ1 p h + ξ2 p h2 . (4.5)

In the ideal case (as is the case for the analytical solution), the two correction factors c0 and c1 equal 1. In Table 4.1, the results of the least squares fit are given. Also the analytical values of the four parameters as derived in Appendix E.4 are given. In Figure 4.6, the stiffness as experimentally identified is given in the left and the right part by the solid surface. In the left part, also the analytical solution of the parameters as derived in Appendix E.4 is given by the meshed surface. In the right part of Figure 4.6, also the fit of the experimental data on (4.5) is given. Even though the numerical values of the analytical and the fitted experimental parameters differ a lot, the analytical and the experimental model show similar results. The main difference is the stiffness – length relation. The analytical model gives only a slight increase in the stiffness for increasing length, while according to the experiments, the stiffness increases substantially with the muscle length h. The large difference of λ for the analytical and the experimental model, is a result of using a least squares fit to derive the experimental model. As can be seen in the right part of Figure 4.6, the minimum of the quadratic relation of the experimental model is located at approximately h = 160 mm (see also Figure G.6). When considering the experimental data (solid surface), a minimum around h = 100 mm seems more realistic. This difference yields the fitted λ to be higher. This also explains part of the factor 10 difference between the experimental and analytical ξ parameters. parameter

λ [bar]

ξ0 [m]

ξ1 [–]

ξ2 c0 [m−1 ] [–]

experiment analytical

3.97 0.50

0.442 0.055

−5.06 12.80 −0.47 1.97

0.82 1.00

c1 [–]

ζ0 [N s m−1 ]

ζ1 [m s]

0.85 1.00

9.55 –

1.94 –

Table 4.1: The six stiffness parameter values as experimentally identified and fitted on (4.5) are given in the upper row, the analytical solution as derived in Appendix E.4 in the lower row. Also the experimental damping parameters are given.

43

1.2

p¯A B [bar]

−1

0

4

2

1.

−30

0.8

1.

2

8

0

1

0.6

10 1.

4

−20

0.4

1.

2

−40

0.

0.

6

4

−1 0

p¯B [bar]

6

1.

1

θ¯2 [deg]

1.

0.2

10

0

0. 4 30

0.2

0.

8

20 40 0.4

1

20

30 40

0.6

0.8

1

1.2

p¯A [bar] Figure 4.7: Solution for θ¯2 = f (¯ pA , p¯B ). Also the cumulative pressure p¯A B is indicated.

2.2

p¯A [bar]

2 1

p¯B [bar]

1

0.6

1

0.6

0.6 0.5 0.4

4

0.7

0.6 0.5 0.4 0.3

0. 2 0.2

0.3

0.4 0.2 −40

3 0.

0.8

0.8

0.9

0.

p¯A B [bar]

1.2

0. 8

0.4

0.9

0.5

0.5

1.4

0.3

1.6

0.7 0.6

1.8

−30

−20

−10

0

10

20

30

θ¯2 [deg] Figure 4.8: The required pressures p¯A and p¯B for a setpoint given by θ¯2 and p¯A B .

44

40

Identification and model validation

Considering these results, the linear stiffness – pressure relation in the model is good. The quadratic stiffness – length relation as derived in the model is present in the muscles, although it is hardly present when using the analytical parameters. This might indicate that the volume – length relation in the muscle model is not correct, as the ξ parameters are derived from this relation. Concluding, relation (4.1) seems promising but the volume – length relation needs reconsideration. Based on the most accurate fit, the parameters derived by fitting the experimental data to (4.5) are used in the remainder of this work; i.e. the upper row in Table 4.1. The dynamic experiments are also used to determine the damping in the muscles, see Appendix G. The measurements are fitted to the model as given in (4.2), yielding the two ζ parameters as given in Table 4.1 (and shown in Figure G.7). It should be noted that the dynamic muscle test is not very suitable to determine the damping, the damping is not explicitly accessed by the small perturbations. A better method to measure the damping directly, is by using a quasi-dynamic test; an indirect method is to use the FRF of a load attached to the muscle. It is decided to estimate the damping using the FRF’s from the robotic arm identification; in this way, the energy dissipation in both the muscles and in the joints of the robotic arm due to friction, are considered. The current result are just an indication of the muscle damping.

4.2 Validation of the equilibria The equilibria of the robotic arm require the valves to be closed, i.e. u = 0 and v = 0. This is easy to understand by realizing that u is in fact a flow reference signal. As a result, u affects the time derivatives of the pressures. This is confirmed by the low frequency −1-slope of the magnitude of the FRF of the robot, as is presented in next section. The equilibria of the model are derived by x˙ = 0 and u = 0, with x˙ as defined in (3.10) and the stiffness according to the experimental parameters of Table 4.1. It is easy to see that the equilibrium values: h iT x ¯ = p¯c p¯A p¯B θ¯2 θ¯˙2 , (4.6) are given by: p¯c p¯A p¯B θ¯2 ¯ θ˙2

= ∈ ∈ =

pin £ ¤ 0 pin £ ¤ 0 pin f (¯ pA , p¯B )

=

0.

(4.7)

The parameter pin is the supply pressure to the capacity; in our case pin = 1.0 bar. The pressures p¯A and p¯B are both manually fixed to any value in the interval [0 , pin ]. Angle θ¯2 is defined by the two muscle pressures; i.e. any combination of the p¯A and p¯B results in a stable equilibrium position θ¯2 , where f (¯ pA , p¯B ) is the expression for θ¯2 following from the following equilibrium equation: QA (¯ pA , θ¯2 ) + QB (¯ pB , θ¯2 ) + m2 g l12 sin(θ¯2 ) = 0 . 45

(4.8)

p¯A [bar] p¯B [bar] model experiment

2.2 2 1.8

θ¯2 [deg] p¯A B [bar] model experiment

1.2

1

0.8 1.4

p¯B [bar]

p¯A B [bar]

1.6

1.2 1 0.8

0.6

0.4

0.6 0.4 0.2 −40

0.2 −20

0

20

40

0.2

0.4

θ¯2 [deg]

0.6

0.8

1

1.2

p¯A [bar]

Figure 4.9: The equilibrium points of the robotic arm, as predicted by the identified model and as experimentally verified on the robotic arm. The equilibria are set according to θ¯2 and p¯A B as shown in the left figure, the accompanying pressures p¯A and p¯B are shown in the right figure. See also the Figures 4.7 and 4.8.

p¯A [bar] p¯B [bar] model experiment

2.2 2 1.8

θ¯2 [deg] p¯A B [bar] model experiment

1.2

1

0.8 1.4

p¯B [bar]

p¯A B [bar]

1.6

1.2 1 0.8

0.6

0.4

0.6 0.4 0.2 −40

0.2 −20

0

20

40

0.2

θ¯2 [deg]

0.4

0.6

0.8

1

p¯A [bar]

Figure 4.10: The equilibrium points of the robotic arm, as predicted by the estimated model and as experimentally verified on the robotic arm. Compared to Figure 4.9, the estimated stiffness parameters are used, instead of the experimentally derived stiffness parameters. Remark that the contour lines do not resemble the contour lines in Figures 4.7, 4.8 and 4.9.

46

1.2

Identification and model validation

No analytical solution of (4.8) is presented here. Instead, a graphical interpretation of θ¯2 = f (¯ pA , p¯B ) is given in Figure 4.7; in this figure also the cumulative equilibrium pressure p¯A B = p¯A + p¯B is indicated; this will appear to be useful in Chapter 5. In Figure 4.8, the ”inverse” of Figure 4.7 is presented. It shows the required pressures p¯A and p¯B which yield a desired combination of θ¯2 and p¯A B . The equilibria, as predicted by the model, are verified on the robotic arm in twenty equilibria. To simplify these experiments, they are performed in closed-loop. The use of a controller does not affect the location of the equilibria as in an equilibrium u = 0; i.e. the controller only guarantees that both pressures are met. In the left part of Figure 4.9, the twenty equilibria (given by θ¯2 and p¯A B ) are indicated. These equilibria are realized on the robotic arm (green circle) and in the model (magenta circle). In the right part of Figure 4.9, the pressures pA and pB belonging to the equilibrium setpoints are given. As can be seen, there is a discrepancy, indicated by the black line, between the pressures in the robotic arm and as predicted by the model. The model predicts accurately in case θ¯2 ∈ [−10, 10] deg and p¯A B ∈ [0.4, 1.0] bar. In an equilibrium, the damping is not important, i.e. the muscle force is equal to the force due to the stiffness. Each equilibrium is determined by the muscle (stiffness) forces and the gravity force, i.e. these forces cause the mismatch between the model and the experiments: • If the stiffness k(p, h) is wrongly predicted in the model, this causes an error in the muscle force and in the equilibria. The stiffness, however, is experimentally determined and for this reason considered to be reliable; namely, the FRF data as discussed in Appendix G show to be reliable. The experiments to identify the muscle stiffness are performed inbetween 0.5 and 1.4 bar. In the equilibria, the muscle pressure is inbetween 0.2 and 0.9 bar. As the stiffness parameters are fitted on the data between 0.5 and 1.4, this fitted model might predict the stiffness wrong below 0.5 bar. • The stiffness force also depends on the elongation of the muscle h − h0 (p). The muscle length h is estimated correctly, using the kinematic relations of the robotic arm. A possibly wrong estimation of the nominal length h0 (p), causes errors in the muscle force. In Figure 3.12, the modeled h0 (p) seems accurate; however, the muscles show a variation of at least 10 mm in their nominal length. • The gravity force in the model is expected to be correctly estimated, as this part involves no doubtful assumptions. • Dry friction in the revolute joints my cause a deadzone, but this is unlikely to cause an angular shift of 10 degrees in just a specific part of the application range. Besides, the muscles are under pretension. In order to check the robotic arm model setup, the experimental results are also compared to the model, now using the estimated stiffness parameters of Table 4.1. The equilibria of this model are given in Figure 4.10. The match between the model and experiments is striking. Unfortunately, the estimated stiffness is considered to be too low, but according to this figure the model setup seems correct. Concluding, the mismatch between the (fitted) model and the experiments is probably caused by an error in the estimate of the nominal length (due to the variation between individual muscles) and an error in the estimate of the stiffness parameters. The latter error is caused by the use of a least squares fit over an application range that does not resemble the application range in the robotic arm. 47

Nominal muscle behavior A second important indicator to verify the equilibria is the nominal behavior of the muscles. The nominal behavior of the muscles is explained in Section 3.4.4 and shown in Figure 3.12. The match between the experimental and modeled nominal length and diameter is good, although there is much variation between the muscles. The muscle parameters L, n, s1 and s2 are very dominant in the solution of (4.8) and the nominal behavior. The discrepancy of the modeled and experimental equilibria might be improved by reconsidering the four muscle parameters, but this will be at the cost of a less accurate match of the nominal behavior.

4.3 Robotic arm identification In Chapter 5, a control strategy is proposed that stabilizes the robotic arm. Two reference variables are used in the controlled DoF robotic arm: the angle θ2 and the cumulative pressure pA B = pA + pB , which are assumed to be completely decoupled. To control these two new outputs, an angle input zθ and a cumulative pressure input zp are defined. The identification of the robotic arm is given by the FRF of angle θ2 over input zθ around an equilibrium (θ¯2 , p¯A B ). The identification is performed in 24 equilibria, which are indicated by circles in Figure 4.9. Each equilibrium is a defined by an angle θ¯2 and a cumulative pressure p¯A B . The FRF’s are measured in closed-loop to keep the robotic arm from drifting away from its reference. For the purpose of identification, noise is injected in the θ2 loop while the cumulative pressure pA B is kept constant. The general setup of the closed-loop identification procedure is given in Figure 4.11. A plant Hx is identified for output x around an equilibrium reference xref , using a sensitivity analysis on the experimental data. Using the FRF of controller C and the power spectra Φe1 and Φe2 of the error signals e1 and e2 , Hx is derived by: Hx

= −

Φe1 . C Φe2

(4.9)

The FRF of the plant around an equilibrium (θ¯2 , p¯A B ) for the angle θ2 , is given by Pθ = zθθ2 , with zθ the control signal. First, a sensitivity analysis, using a simple gain controller C is performed, this yields two angular error signals e1 and e2 . The next step is to calculate the power spectra Φe1 and Φe2 of these error signals. Finally, the FRF Pθ is derived using the power spectra and by eliminating the controller, as is shown in (4.9). n

e1 xref

+

+ −

+

e2 zx

x Hx

C

Figure 4.11: Schematic representation of a sensitivity analysis of output x over input zx , around equilibrium x ¯. A known disturbance n is introduced in the control-loop, using controller C. Both error signals e1 and e2 are measured and used to derive Hx = zxx .

48

Identification and model validation frequency response function Pθ =

θ2 zθ

power spectra Φ

0

10

|Pθ | [rad V−1 dB]

0

−10

−2

10

|Φ| [dB]

−20 −30

−4

10

−40 −50

Φxx [V2 ] Φyy [rad2 ] Φxy [V rad]

−6

−60 −70 −80

Pθ bθ P

−1

10

10

−8

0

10

10

1

−1

10

10

0

10

1

10

coherence 1 0.8

90

2 | [–] |γhF

phase(Pθ ) [deg]

180

0.6

0

0.4

−90 0.2 −180 −1

0

10

10

0

1

10

−1

10

f [Hz]

0

10

1

10

f [Hz]

Figure 4.12: Experimental sensitivity analysis around equilibrium point: θ¯2 = 0 deg and p¯A B = 1 bar. The indices x and y are defined by x = −C e2 and y = e1 . The third-order bθ (s) as in (4.10), with α = 780, ω = 6.3 Hz and ς = 0.05 is fitted transfer function P on Pθ = zθθ2 .

In Figure 4.12, the result of the sensitivity analysis is shown for the equilibrium point with θ¯2 = 0 deg and p¯A B = 1.0 bar. The experiment shows a good coherence between f = 0.2 Hz and f = 10 Hz and is considered to be reliable in this range. The FRF is characterized by a low-frequency −1-slope, an eigenfrequency at approximately 6 Hz and a high-frequent −3-slope. A standard spring-mass-damper system is characterized by a constant gain below the eigenfrequency and a −2-slope beyond the eigenfrequency. The additional −1-slope of the robotic arm is caused by the input u, which directly affects the time derivative of the pressure, thereby introducing an additional integrator. For the moment, we assume that the two muscles act like a force actuator; the dynamics of the robotic arm around an equilibrium is modelled by the following third-order transfer function: 1 α α = 3 . (4.10) s s2 + 2 ς ω s + ω 2 s + 2 ς ω s2 + ω 2 s The gain α depends on the flow (and on the numerical values of the pressure and the angle); the damping ration is given by ς. The natural eigenfrequency ω (in radians per second) depends on the sum of the stiffness of muscles A and B , i.e. kA B = kA + kB and on the moment of inertia Jzz2 of body 2 according to: q −1 kA B (l12 − l2A )2 Jzz2 , (4.11) ω = bθ (s) P

=

with l12 − l2A the distance between the muscle mounting and joint 2 ; this is correct as muscles A and B are mounted equidistantly from joint 2 , see Table 3.4. 49

Figure 4.13: The FRF of Pθ = zθθ2 for 24 equilibrium points of the 1DoF robotic arm. In the legend, the equilibrium points are given in a short handed way; i.e. θ -30 p 0.8 indicates θ¯2 = −30 deg, p¯A B = 0.8 bar.

50

0

−10

−1

10

−180

−90

0

90

180

−80 −1 10

−70

−60

−50

−40

−30

−20

|Pθ | [rad V−1 dB]

phase(Pθ ) [deg]

0

10

0

10

f [Hz]

frequency response function Pθ =

θ2 zθ

1

10

1

10

θ -30 θ -15 θ 0 θ 15 θ 30 θ -30 θ -15 θ 0 θ 15 θ 30 θ -30 θ -15 θ 0 θ 15 θ 30 θ -30 θ -15 θ 0 θ 15 θ 30 θ -15 θ 0 θ 15 θ 0

p 0.8 p 0.8 p 0.8 p 0.8 p 0.8 p 1.0 p 1.0 p 1.0 p 1.0 p 1.0 p 1.2 p 1.2 p 1.2 p 1.2 p 1.2 p 1.4 p 1.4 p 1.4 p 1.4 p 1.4 p 1.6 p 1.6 p 1.6 p 1.8

Identification and model validation

bθ according to (4.10) is fitted on the experimental FRF Pθ . The fit is The model P shown in Figure 4.12; it allows us to estimate the damping ratio ς of the robotic arm. The damping in the robotic system is very low: ς = 0.05. This yields that the damped and the natural eigenfrequency are approximately equal: ω = 6.3 Hz. In Figure 4.13, the resulting FRF’s Pθ are given for the 24 equilibria. Six conclusions can be drawn from this figure: • The FRF’s indicate a distinct different behavior in each equilibrium point; this proves that the robotic arm behaves in a nonlinear manner. • The mass line (with slope −3) is identical in each equilibrium. As a result, it can be concluded that the muscles act like a force actuator. This is in accordance with the assumption of the third-order transfer function in (4.10). If the muscles would act like a position actuator, the actuator part in relation (4.10) would have been modelled with ω 2 in the numerator, instead of 1. In that case, a change in stiffness will cause variations in the high-frequency behavior and a constant low-frequency behavior. • The stiffness line (with slope −1 due to the pneumatics) varies with a factor 3 (10 dB). This variation is a result of the variation in muscle stiffness and of the gain α. The gain α depends on the pressure. • The damping present in the robotic arm is low; the damping ratio ς varies inbetween 0.04 and 0.10. The damping ratio is an indication of the energy dissipation in the robotic arm. Energy is dissipated by damping but also by friction in the muscles and in the joints. It can be concluded that the dry friction in the muscles is low. This conclusion is supported by the good coherence as shown in the sensitivity analysis of Figure 4.12. In case of a nonlinear phenomena such as dry friction in the muscles, the coherence would not equal 1. • Increasing the cumulative pressure p¯A B , causes an increase of the eigenfrequency ω due to the increase in stiffness. As the muscle is considered to be a force actuator, this causes the magnitude of Pθ to decrease. The decrease in Pθ magnitude is amplified by the pressure. An increasing cumulative pressure yields the flow to decrease, yielding a decrease in α. • Variations in the angle θ2 from 0 degrees, hardly affect the behavior of the robotic arm when being compared to the effect of the cumulative pressure. The eigenfrequency ω slightly increases, while the magnitude of Pθ remains constant; this means that the gain α must increase. This is caused by the gravity force; in case θ¯2 6= 0 deg, the gravity force becomes part of the equilibrium, which requires to be compensated for by the muscles. This additional force yields a higher gain, without affecting the eigenfrequency. As p¯A B remains equal, the pressure does not affect the gain via the flow. The relations given in (4.10) and (4.11) can be used to estimated the stiffness from the robotic arm identification. For the FRF around (θ¯2 , p¯A B ) = (0, 1), this is shown in Figure 4.12; using the parameters as defined in Table 3.4 yields kA B = 7.8 kN m−1 or kA = kB = 3.9 kN m−1 . Each of the muscles in this equilibrium is given by pA = pB = 0.5 bar and hA = hB = 190 mm; according to the stiffness model in (4.5), this yields kA = kB = 2.9 kN m−1 . This proves that in the given equilibrium the model predicts the stiffness too low. 51

Repeating this analysis for each of the 24 FRF’s in Figure 4.13, will allow for the verification of the stiffness as identified in Section 4.1. It also gives a proper indication of the damping present in the muscles and the robotic arm and of the pneumatic gain. This should be verified in future work. At the end of Section 4.1, it is stated that we prefer to estimate the damping using the measured FRF of the robotic arm. By comparing the linearized behavior of the model with the FRF of the robotic arm, it appears that the muscle damping as defined in Table 4.1 is indeed too low to describe all energy dissipation in the robotic arm, the additional damping is probably due to dissipation in the bearings. It is to decided to take all damping into account in the damping parameter ζ0 , as introduced in the muscle model in (3.49). Based on the simulations, this yields a damping parameter: ζ0

= 40

Ns m−1 .

(4.12)

In the remainder of this report, this value will be used in the robotic arm model.

4.4 Validation of the linearized model The second part of the model validation is to compare the linearized behavior of the model with the frequency response function of the robotic arm. The model of the robotic arm is given in (3.10). The linearization around an equilibrium point (¯ x, z¯) is defined by: ¯ ¯ ∂f (x, z) ¯ ∂f (x, z) ¯ ¯ x ¯ z˜ , x ˜ + (4.13) ˜˙ = ∂x ¯x¯,¯z ∂z ¯x¯,¯z

¯+x ˜ and z = z¯ + z˜ and where x ˜ and z˜ are small perturbations. Rewith x = x mark that in an equilibrium, no flow is possible, i.e. the input signals are given by z¯ = [¯ zθ z¯p ]T = 0. In order to be able to linearize a model, a continuously differentiable vectorfield f (x, z) in (4.13) is required. This is a problem. Due to the strong interweaving of the binary valve signals v in the model of the robotic arm, linearization is not possible. The relation determining the binary valve signal consists of a saturation and binary operators, which makes it not possible to calculate a derivatives of these functions. This issue is solved by performing a sensitivity analysis of the closed-loop model in simulation. The reason to simulate in closed-loop is to avoid the model from drifting from its equilibrium point. The closed-loop simulations are performed in a similar fashion as the experimental sensitivity analysis on the robotic arm, as discussed in the previous section. The FRF of the angular behavior of the model linearization around the equilibrium (θ¯2 , p¯A B ), is defined as Mθ = zθθ2 , with zθ the angular control signal. The model is simulated in closed-loop, using a simple gain controller, while injecting noise in the control loop, see Figure 4.11. The noise consists of an increasing and decreasing sweep signal between 0.1 and 20 Hz. The simulation data is used for a sensitivity analysis, which is given in (4.9). The FRF of the angular behavior of the linearized robotic arm model Mθ can be estimated around every equilibrium defined in (4.7). 52

Identification and model validation

frequency response function H

|H| [rad V−1 dB]

0

power spectra Φ

0

10

−10

−2

10

|Φ| [dB]

−20 −30

−4

10

−40 −50 −60 −70 −80 −1 10

Mθ bθ P cθ M

Φxx [V2 ] Φyy [rad2 ] Φxy [V rad]

−6

10

−8

0

10

10

1

10

−1

10

0

10

1

10

coherence 1 0.8

90

2 | [–] |γxy

phase(H) [deg]

180

0.6

0

0.4

−90 0.2 −180 −1

10

0

10

0 −1 10

1

10

f [Hz]

0

10

1

10

f [Hz]

Figure 4.14: The simulated sensitivity analysis around equilibrium (θ¯2 , p¯A B ) = (0, 1). The indices x and y are defined by x = −C e2 and y = e1 . A third-order transfer cθ (s) with α = 1150, ω = 5.8 Hz and ς = 0.07 is fitted on Mθ . The function M bθ as shown in Figure 4.12, is also indicated. estimated experimental FRF P

In Figure 4.14, the FRF Mθ of the model linearization around equilibrium cθ on Mθ is shown; (θ¯2 , p¯A B ) = (0, 1) is shown. In this figure, also the model fit M the model is defined in (4.10). This allows us to estimate the natural eigenfrequency, the damping ratio and the gain; i.e. ω = 5.8 Hz, ς = 0.07 and α = 1150. bθ as is shown in Figure 4.12; The third FRF shown in Figure 4.14, is the fitted FRF P bθ is defined around the same equilibrium. The model shows a slightly different P bθ and FRF when being compared to the experiment; from the fitted parameters of P c Mθ , it can be seen that ω is too low and α is too high in the model:

• The gain error is probably due to model errors in the pneumatic part of the model. In case the flow towards the muscles is modeled larger than the actual flow, then the response of the muscles and the robotic arm is quicker, yielding a higher gain in the FRF. It must be noted that the flow is very sensitive to small variations in the β parameters in (3.22).

• The natural eigenfrequency depends on the inertia properties of body 2 and the stiffness of the muscles. The inertia of body 2 is considered to be correctly estimated. This means that the stiffness is estimated too low. The stiffness depends on the pressure and the muscle length. It is very likely that differences between the muscles causes the stiffness to be estimated too low. The variation in stiffness is approximately 20%, as can be seen in Figure 4.5. 53

frequency response function H

power spectra Φ

2

10

0

0

10

|Φ| [dB]

|H| [rad V−1 dB]

20

−20

−2

10

−40 −60 −80 −1 10

Mθ bθ P cθ M

Φxx [V2 ] Φyy [rad2 ] Φxy [V rad]

−4

10

−6

0

10

10

1

10

−1

10

0

10

1

10

coherence 1 0.8

90

2 | [–] |γxy

phase(H) [deg]

180

0.6

0

0.4

−90 0.2 −180 −1

10

0

10

0 −1 10

1

10

f [Hz]

0

10

1

10

f [Hz]

Figure 4.15: The simulated sensitivity analysis of the equilibrium (θ¯2 , p¯A B ) = (30, 0.8). The indices x and y are defined by x = −C e2 and y = e1 . The third-order transfer cθ , with αM = 1900, ζ M = 0.07 and ω M = 4.7 Hz, is fitted on Mθ . The function M bθ in this equilibrium is given by αP = 830, ζ P = 0.10 and experimental FRF fit P ω P = 6.4 Hz. Using (4.11) the measured cumulative stiffness is kAPB = 8.0 kN m−1 , while the modeled cumulative stiffness is much lower: kAMB = 4.3 kN m−1 .

• As the damping ratio in model and measurement are almost equal, there is no reason to assume errors here. The model and the experiment in (θ¯2 , p¯A B ) = (0, 1) show identical behavior above 6 Hz. The difference below 6 Hz is approximately 4 dB or 30%. This is acceptable when considering that the variation in muscle stiffness is 20%, see Figure 4.5. The linearization Mθ in equilibrium (θ¯2 , p¯A B ) = (30, 0.8), together with the fitted cθ and P bθ are given in Figure 4.15. The model and the robotic arm show FRF’s M a big difference. The magnitude of the model is approximately 12 dB or a factor 4 higher than the experiment and the eigenfrequency of the model is 4.7 Hz, while the experiment yields an eigenfrequency of 6.5 Hz. In Section 4.2, it is shown that the equilibria are not predicted well for |θ¯2 | > 20 deg. The dynamic behavior in (θ¯2 , p¯A B ) = (30, 0.8) also shows errors. Based on the error in the eigenfrequency and magnitude in the model, we conjecture that the stiffness of the muscles is modeled too low. This is not in accordance with the conclusion in Section 4.2, which states that the stiffness is modeled accurately. Secondly, the gain α is modeled too large. This is likely to be due to errors in relations describing the flows and muscle pressure in the model. Five possible errors are distinguished: 54

Identification and model validation

• The flow is very sensitive to variations in its β parameter; a small error in the restriction diameter has a big impact on the flow. • It is assumed (based on experiments) that the flow through the hoses is defined by the relation for a restriction. This is however not verified on the robotic arm; a hose description might give a better result for the flow description in some cases. This should be verified in future work. • The muscle volume might be estimated inaccurately by the barrel description. This has been checked by simulating the model using a cylindrical relation for the volume. It appears that the cylindrical volume description does not yield significantly different behavior compared to the barrel shaped behavior. Since this does not prove the volume is estimated correctly by the volumetric behavior in relation (3.39), this should be verified in future work. • It is assumed that the muscle volume only depends on the length of the muscle. The volumetric behavior is based on the length–volume relation of a muscle with a constant amount of air in the muscle. The pressure is taken into account in the loop: pressure → nominal length → muscle force → robot dynamics → muscle length → volumetric behavior, see Figure 3.2. This indirect pressure–volume relation is not verified. • In the relation describing the muscle flow, the influence of the tube and braid pushing the air out of the muscle, is not modeled explicitly. This behavior is considered to be taken into account in the nominal behavior, through the same loop as given in previous item. This indirect nominal length – flow relation is also not verified. For the simulation in Figure 4.15, the experimental stiffness parameters are used. In Figure 4.10, it is shown that the equilibria of the robotic arm model using the analytical stiffness parameters, yields accurate results. However, the model with the analytical stiffness parameters is unstable. Due to the fact that the analytical stiffness is lower, the gravity force becomes dominant over the muscle forces hereby inducing instability.

4.5 Time domain model validation In previous sections, it has been shown that the model has a limited pressure range in which the behavior of the robotic arm is predicted well. The third part of the validation procedure compares the open-loop time behavior of the robotic arm with that of the model. As can be expected from previous sections, these results are not very good. On the other hand, the open-loop validation confirms the conclusions in the previous sections. The open-loop experiments and simulations are performed by providing a time-varying input signal u(t) to the robotic arm, the pressures pA and pB and the angle θ2 are measured. An identical input u(t) is applied to the model. In Figure 4.16, the match between experiment and simulation is good. As is said before, the model predicts well in the working area near θ2 = 0 deg. If the pressure in the muscles drops below 0.3 bar, the predictive quality of the model worsens; in this case the model predicts a low stiffness, which causes to gravity to become slightly more dominant. As a result, the angle increases. 55

pA pB pA pB

1

p [bar]

0.8 0.6

experiment experiment simulation simulation

0.4 0.2 0

0

1

2

3

4

5

6

7

15

θ2 [deg]

10

experiment simulation

5 0 −5 −10

0

1

2

3

4

5

6

uA uB

1

u

7

0

−1

0

1

2

3

4

5

6

7

time [s] Figure 4.16: The initial angle of the robotic arm is close to zero degrees and the initial pressures are identical. The muscles are actuated in opposite directions.

The experiments shown in Figure 4.17, has an input uA = 0, i.e. the air mass in muscle A is constant. Once pressure pB starts to increase, the arm is pulled back towards θ2 = 0 deg. As a result, muscle A is elongated. The pressure increase due to the elongation (which must be caused by a decrease in the volume), is hardly seen in the simulation. Pressure pB gives a realistic result. The angle prediction of θ2 is less accurate due to the mismatch of pA . The simulation, shown in Figure 4.18, predicts well up to 4.5 seconds. This is due to the pressure difference between both muscles, which is similar for the model and the experiment. After 4.5 seconds, the pressure in the muscles is to below 0.3 bar, which is not part of the model’s predictive pressure range. In the model, the arm becomes unstable as the gravity force is stronger than the muscles forces; the arm is pulled out of position.

4.6 Discussion The major conclusion of this chapter is the limited predictive impact of the robotic arm model, especially in open-loop. Only in a specific working area, the model appears to predict the behavior correctly. This can be seen in both the equilibria as in the linearized behavior around these equilibria. The four major points of concern in this respect are: Stiffness The stiffness of the muscles is identified in a separate identification procedure. Despite the fact that the parameters of the stiffness relation are fitted to the results of these experiments and used in the robotic arm model, the stiffness is still considered to be estimated too low. 56

Identification and model validation pA pB pA pB

1

p [bar]

0.8 0.6

experiment experiment simulation simulation

0.4 0.2 0

0

1

2

3

4

5

6

7

50

experiment simulation

θ2 [deg]

40 30 20 10 0

0

1

2

3

4

5

6

uA uB

1

u

7

0

−1

0

1

2

3

4

5

6

7

time [s] Figure 4.17: The initial pressure in the experiment of muscle B is pB = 0 bar; this means that muscle B is deflated. As this behavior is not modeled, it is decided to give the model a stable initial condition; θ¯2 and p¯A are equal to the experiment, p¯B is set to an equilibrium as shown in Figure 4.9.

Pneumatics The relations determining the pressure and flow are very sensitive too small errors in the flow parameters. If the muscle pressure is not predicted accurately, this yields errors in almost every other part of the model. Model errors The muscle model is considered to describe the muscle behavior accurately. Concerns exist in the relations describing the stiffness and the flow. The damping is not physically modeled but is also less important, for this reason viscous damping is an appropriate assumption. Muscles The muscles show a large variation when being compared to each other. No two single muscles are identical, or yield identical behavior. These variations result in for example a different nominal length, which is very important in predicting the exact muscle force. The combination of these four factors, limit the predictive capacities of the robotic arm model. Still, the model as presented in Figures 3.1 and 3.2, is considered to describe the behavior correctly on a phenomenological level. This new model does give insight of the working principles of the robotic arm, but it unfortunately fails to give an accurate description on a quantitative level. The experimental identification of the robotic arm has lead to a lot of insight in the behavior of the robotic arm and of the muscles. Despite the inherent nonlinear behavior of the robotic arm, it is shown that the eigenfrequency of the robotic arm 57

pA pB pA pB

1

p [bar]

0.8 0.6

experiment experiment simulation simulation

0.4 0.2 0

0

1

2

3

4

5

6

7

5

6

7

60

experiment simulation

θ2 [deg]

40 20 0 −20

0

1

2

3

4

uA uB

u

1

0

−1

0

1

2

3

4

5

6

7

time [s] Figure 4.18: In this experiment, the initial pressures and angle are similar for the experiment and the model. Both muscles are actuated in similar ways.

are located in a rather narrow band between approximately 6 and 7 Hz. A second important conclusion is that the muscles appear to be force actuators; probably because they are applied with a certain pretension range in the robotic arm. The results of the experimental identification of the robotic arm, as discussed in Section 4.3, are accurate and reliable. It is decided to use the identification and the gained insight to construct a control strategy for the robotic arm. This is elaborated on in the following chapter.

58

5 Feedback control of the robotic arm

This chapter presents the control of the DoF and the 2DoF robotic arm. In Section 5.1, the control strategy used for the DoF robotic arm is discussed. Using this strategy and the results of the experimental identification, controllers are designed; this will be explained in Section 5.2. In case the robotic arm is expanded by the second joint, the control becomes slightly more complicated due to interference of the feedback loops. Moreover, a different control strategy is required for the 2DoF robotic arm due to hardware limitations; only a limited number of sensor signals can be redirected to the dSPACE system. The 2DoF control strategy and the controller design are discussed in Section 5.3. In Section 5.4, an evaluation of the linear control strategy proposed here and the controller based on the reinforcement learning, as discussed in [Maas, 2005], is given. This chapter is ended by a discussion on the control results.

5.1 Single joint control strategy In the introduction, it is stated that to control the robotic arm linear feedback controllers will be used. In Section 4.3, the nonlinear nature of the robotic arm is identified. The two most important nonlinear contributions in the robotic arm are the binary valves and the McKibben muscles. By applying a PWM algorithm, the binary valves approximate linear behavior. The nonlinear McKibben muscle behavior is accounted for by the control design by ensuring sufficient robustness margins in the application range of the robotic arm. The input to the PWM of the DoF robotic arm R1 , are two mass flow signals uA and uB . The state x of R1 is defined by the joint angle θ2 and the two muscle pressures pA and pB . In matrix notation, this gives: u

=

x =

£ £

uA θ2

uB pA

¤T

pB

(5.1) ¤T

.

(5.2)

θ¯2 simulation

1.2 −40 −30 −20

−10

p¯A B = p¯A + p¯B θ¯2 experiment

1

0 −30

−15 0

p¯B [bar]

0.8

0.6

15

0.4

30

10

20 30

0.2

40 0.2

0.4

0.6

0.8

1

1.2

p¯A [bar] Figure 5.1: The equilibrium points of the robotic arm, as predicted by the model and as experimentally verified on the robotic arm. See also the Figure 4.9.

−

θ2ref

eθ

+

θ2

Pθ

plant P

+ pAref B

zθ

Cθ

Cp

ep

Pp

zp

pA B

−

zθ Ψ

v

u

x R1

PWM

zp

θ2 Ω pA B

Figure 5.2: The single joint control loop. The plant P is given by two FRF Pθ = zθθ2 and Pp = pzApB , which are considered to be decoupled. The error signals eθ and ep are the input to the controllers Cθ and Cp , which yield the new plant input signals zθ and zp . In the lower figure, the decoupling of the 1DoF robotic arm R1 is depicted, using the decoupling matrix Ψ and the output matrix Ω. The flow control signals u are converted to binary valve signals v by the PWM, as defined in (3.12), which are the input to R1 .

60

Feedback control of the robotic arm

The purpose of the robotic arm is to make controlled motions, i.e. to control the joint angles. For this reason it is decided to use angle θ2 as an output of R1 . The angle θ2 is affected by the momentum generated by muscles A and B . This momentum depends on the difference between the two muscles forces. The easiest way to generate a momentum with two muscles is by inflating one muscle and deflating the other muscle for an equal period of time, i.e. uA (t) = −uB (t). To have a balanced number of inputs and outputs, a second output signal needs to be defined. It would be most convenient to have an output that is decoupled from θ2 , as this would allow for a control strategy using two SISO controllers. Output θ2 is affected by uA = −uB ; to ensure the second output to be decoupled from θ2 , it must not be affected by uA = −uB . In case the volume of the muscles is contant and the pressure in both muscles is 0.5 bar, a change in uA = −uB , does not affect the sum of the pressure in the two muscles. It is decided to use the cumulative pressure as the second output. From this elaboration, the outputs y are constructed using the state variables (and sensor signals) x as defined in (5.2) and an output matrix Ω: y Ω

= =

£

·

θ2 1 0

pA B 0 0 1 1

¤T ¸

= Ωx

(5.3)

.

(5.4)

The output given by the cumulative pressure pA B = pA + pB is affected if both muscles are inflated for an equal period of time, i.e. uA (t) = uB (t); this dependency is the opposite of the dependency given for θ2 , and for this reason it is expected that outputs θ2 and pA B are decoupled. To verify whether the decoupling is a legitimate assumption, the equilibria of R1 are considered. In Section 4.2, the equilibria are determined for both the model and the experimental setup, see also Figure 5.1. In case the outputs are decoupled, the lines of θ¯2 and p¯A B must be perpendicular in this figure; in that case a change in one output does not affect the other. For the model simulation, this is only true if θ¯ = 0 deg. For the experiments, however, it can be seen that θ¯2 and p¯A B are approximately perpendicular, within the entire application range. For this reason, the outputs y are considered to be decoupled. The next step is to define two new inputs, the angle input zθ and the cumulative pressure input zp , which address the decoupled outputs only: z

=

£

zθ

zp

¤T

.

(5.5)

Using the previous stated dependencies of the new inputs z on the flow control signals u yields the decoupling to be defined by: ) uA = zθ + zp → u = Ψz (5.6) uB = −zθ + zp Ψ =

·

1 1 −1 1

¸

.

(5.7) 61

In Figure 5.2, the control strategy for the DoF robotic arm is presented schematically. The plant P is defined by two new inputs z and two decoupled outputs y. The behavior of P is given by the two FRF’s: Pθ

=

θ2 zθ

, Pp

=

pA B . zp

(5.8)

The decoupling allows for the application of two SISO controllers on the plant: one angular controller Cθ and one cumulative pressure controller Cp . In the lower part of Figure 5.2, the conversion of R1 to P, using the decoupling matrix Ψ and the output matrix Ω is graphically depicted.

5.2 Single joint control The controller for the DoF robotic arm consists of two linear controllers Cθ and Cp , as shown in Figure 5.2. These two controllers are designed using the respective plants Pθ and Pp as defined in (5.8). Both controllers are tuned such that the robotic arm is locally stable around each of the measured equilibria. We will show (not prove) that this control strategy also stabilizes the nonlinear system. In Appendix H.2, the FRF’s of Pp around seven equilibria with θ¯2 = 0 are given. The seven FRF’s are almost equal, i.e. the cumulative pressure p¯A B hardly affects the gain of the FRF. The magnitude of the FRF shows a clear −1-slope over the complete frequency range. This is caused by the input u which directly affects the time derivative of the pressure. In Figure 4.13, the FRF’s of Pθ in 24 equilibria are given. This figure shows that the eigenfrequency of the robotic arm varies between 6 and 7 Hz; This limits the feasible bandwidth of the controlled robotic arm; namely, the linear controllers must be able to stabilize the robotic arm in every possible equilibrium. The servo performance is hampered by two types of disturbances: external and self-induced disturbances. The latter are due to a number of reasons: interaction between the two outputs, friction in the muscles and bearings, excitation due to the switching of the binary valves and the nonlinear stiffness change of the muscles. To avoid instability, four precautions are taken: • The controllers are designed for the worst case scenario, the open-loop behavior in every equilibrium point is used to design stable controllers. • In the PWM algorithm, a dead-zone value τ is introduced, see (3.13d). This dead-zone yields that small static errors in both angle and pressure are allowed. This keeps the binary values from continuous actuation and reduces the amount of disturbances on the robotic arm. • Wide robustness margins are used. A sensitivity margin of 6 dB, a phase margin of 45 degrees and an amplitude margin of 5 dB with respect to the most critical FRF’s are employed for the open-loop behavior. • The controllers are designed such that the angle and pressure bandwidth are approximately equal. This gives outputs θ2 of pA B the same priority. 62

Feedback control of the robotic arm

5.2.1 Angular controller design In Figure 5.3, the FRF’s of 24 measurements Pθ are given, together with the controller and the open-loop frequency response. Controller Cθ consists of a lead-lag filter, a pure notch filter and a second-order low-pass filter, the parameters are given in Table 5.1. The definition of the filters is given in Appendix H.1. Essential in the tuning of the controller is the use of the notch filter. By stretching the notch filter wide and deep (−20 dB), the eigenfrequency of the robotic arm in the entire working area is attenuated by the controller. For stability reasons, the notch frequency fn is a little lower compared to the eigenfrequency. A lead-lag filter is placed on top of the notch to prevent loss of phase margin. The bandwidth is inbetween 1 and 3 Hz, as indicated by the vertical pink lines. The amplitude and phase margin are met in every case. In the lower part of Figure 5.3, it can be seen that the sensitivity is below 6 dB for almost every frequency in each equilibrium. In Figure 5.4, the Nichols diagram is given. To guarantee stability, the HOL curve must keep the origin (−180, 0) at the left if moving from low to high frequencies. The robustness criteria are indicated by the circle around the origin. It can be seen that stability and robustness are guaranteed for every equilibrium.

5.2.2 Pressure controller design In Figure 5.5, the controller Cp , the pressure FRF’s Pp and the open-loop behavior are shown. The experimental Cp data is discussed in Appendix H.2. Controller Cp consists of a gain and a first-order low-pass filter. The two parameters of this controller are given in Table 5.2. This controller is rather simple. A controller yielding a better performance is possible. However, the bandwidth of this pressure loop is now approximately 3 Hz, which is sufficient as it is equal to the bandwidth of the angular control loop. The phase and amplitude margins are both met and the sensitivity is always below 6 dB. The local stability is also guaranteed; this is proven in the Nichols diagram as shown in Appendix H.3. gain p

55

lead-lag f1 f2

1 15

notch fn β1 β2

low-pass

6.1 0.05 5

fl βl

15 0.7

Table 5.1: The Cθ controller parameters for joint 2 of the 1DoF robotic arm. The controller consists of a gain with a notch, lead-lag and a second-order low-pass filter. All frequencies f are given in Hz.

gain p

low-pass 6

fl

10

Table 5.2: The Cp controller parameters for the 1DoF robotic arm. The controller consists of a gain and a first-order low-pass filter. The low-pass frequency fl is given in Hz.

63

40 20

|H| [dB]

0 −20 −40

Cθ [V rad−1 ] HOL [–] Pθ [rad V−1 ]

−60 −80 −100

−1

10

0

10

1

10

phase(H) [deg]

180

90

0

−90

−180 −1

10

0

10

1

10

20 10

|S| [dB]

0 −10 −20 −30 −40

−1

10

0

10

1

10

f [Hz] Figure 5.3: The two upper figures show the magnitude and phase of the controller Cθ , plant Pθ and the open-loop HOL = Cθ Pθ . The Pθ measurements are the 24 FRF’s previously shown in Figure 4.13. The phase margin is indicated by the pink line at −145 dB in the middle figure. The lower figure shows the sensitivity S = (1 + Cθ Pθ )−1 of the control loop, the sensitivity margin of 6 dB is indicated by the pink line. The bandwidth is inbetween 1 and 3 Hz, as indicated by the two vertical pink lines.

64

Feedback control of the robotic arm

30

20

|HOL | [dB]

10

0

−10

−20

−30 −360

−270

−180

−90

0

phase(HOL ) [deg] Figure 5.4: The Nichols diagram of a number of FRF’s. The open-loop HOL = Cθ Pθ must cross the origin at (−180, 0), at the right if moving from low frequencies (upper right) towards high frequencies. The phase and gain margin are indicated by the circle around the origin. The loop in the curve is introduced by the notch, which is overcompensating some HOL curves. The curves in this figure correspond with the curves with the same color given in Figure 4.13.

65

90

30 20

|H| [dB]

10 0

−10 −20 −30 −40

Cp [V bar−1 ] HOL [–] Pp [bar V−1 ] 0

1

10

10

phase(H) [deg]

180

90

0

−90

−180 0

1

10

10

20 10

|S| [dB]

0 −10 −20 −30 −40

0

1

10

10

f [Hz] Figure 5.5: The two upper figures show the magnitude and phase of the controller Cp , plant Pp and the open-loop HOL = Cp Pp . The phase margin at −145 dB is indicated by the pink line in the second figure. The lower image shows the sensitivity S = (1 + Cp Pp )−1 of the control loop, the sensitivity margin of 6 dB is indicated by the pink line. The bandwidth is approximate 3 Hz, as indicated by the two vertical pink lines.

66

Feedback control of the robotic arm

output

BW [Hz]

ts [s]

tr [s]

Mp

ess

θ2 pA B

1−3 ≈3

0.69 0.30

0.25 0.06

7 deg −

0.10 deg 0.01 bar

Table 5.3: The performance indicators for the 1DoF robotic arm, using the controllers as defined in Tables 5.1 and 5.2. The angle performance is measured for a step of 30 deg, the pressure performance for a step of 0.5 bar.

5.2.3 Performance The performance is measured by the bandwidth (BW), the settling time ts , the rise time tr , the overshoot Mp and the static error ess , [Franklin, 2001]. In Table 5.3, the values of the performance indicators for θ2 and pA B are given . The bandwidth is taken from the FRF graphs, while the four time domain indicators are taken from the step response measurements as given in Figure 5.6. It is remarkable that the pressure is affected by the step in the angle, while the angle is not affected by a step in the pressure. At 9.5 seconds, an infeasible setpoint is offered, the robotic arm stabilizes at a slightly lower angle and pressure. In Figure 5.7, a third-order angle profile is tracked by the robotic arm, the pressure reference is constant. Only a small part of the measurement is shown, as it is highly repeatable. The tracking shows a delay of approximately 0.1 sec and a some overshoot can be seen around θ = 0 degrees. It must be noted, that despite a lot of model uncertainties, non-linear characteristics and a non-perfect decoupling, the motion control performance of the robotic arm is good. The robotic arm is not sensitive to external disturbances and errors are corrected accurately; this is tested by pulling the arm out of position. In Figure 5.7, also a simulation of the experiment is shown. The simulation shows some overshoot for large angles, while the experiment becomes slow. This mismatch is probably due to the friction (damping), which increases if a muscle is highly stretched. Still, the simulations shows much resemblance with the experiment. Especially the behavior of pA B is predicted accurately by the model. Also the tracking error of the angle is similar to the experiment. The controller used in the simulation is slightly different from the controller in the experiment, the notch frequency is set to fn = 4.5 Hz. Before expanding the robotic arm with joint 1, the control strategy as used for the DoF joint 2 system, is used to control the DoF joint 1 system. The outputs of this system are the cumulative pressure pCD of muscles C and D and the angle θ1 . The pressure controller Cp is similar to the controller in Table 5.2, a new angular controller Cθ is designed. This is discussed in detail in Appendix H.4. The bandwidth of the DoF joint 1 system is between 0.5 and 1 Hz.

5.3 Dual joint control Expanding the robotic arm with the second joint has two implications, which need to be considered before the designing of the controllers can be discussed. Firstly, a different control strategy must be used due to hardware limitations and, secondly, the dynamic coupling between the two joints must be considered. 67

40

θ2 [deg]

20

0

−20

−40

reference experiment 0

1

2

3

4

0

1

2

3

4

5

6

7

8

9

10

5

6

7

8

9

10

2

pA B [bar]

1.5

1

0.5

0

time [s] Figure 5.6: Step response or the 1DoF controlled robotic arm. The angle step of 30 deg and the pressure step of 0.5 bar are used to derive the performance indicators in Table 5.3. 50 40

θ2 [deg]

30 20 10 0 −10

reference experiment simulation

−20 −30 −40

0

1

2

3

0

1

2

3

4

5

6

7

4

5

6

7

1.05

pA B [bar]

1 0.95 0.9 0.85 0.8 0.75 0.7

time [s] Figure 5.7: Simulation and experiment showing the error tracking of respectively the model and the robotic arm. The angle reference is a third-order continuous profile, the pressure reference is 1 bar. The used controllers are similar, except for the notch frequency. In case of the model, fn is adjusted to the eigenfrequency of the model.

68

Feedback control of the robotic arm

5.3.1 Dual joint control strategy The state x of the 2DoF robotic arm R2 is given by two joint angles and four muscle pressures. The input to the PWM and R2 are four flow control signals u: £ ¤T uA uB uC uD (5.9) u = £ ¤T θ 1 θ 2 pA pB pC pD x = . (5.10)

Due to hardware limitations, only four of the six state variables x of R2 can be measured. As a result, the control strategy as proposed in Section 5.1 for R1 , can not be used for R2 . This makes it also impossible to decouple the outputs of R2 . To be able to control the joint angles accurately, the two angle signals are chosen as output. This leaves only two pressure signals to be used, one signal for each joint. It is decided to control the pressure in muscles A and C . The notation as introduced for R1 in Section 5.1, is also adopted for R2 . Using the state variables defined in (5.10), the outputs are defined by: £ ¤T θ 1 θ 2 pA pC y = = Ωx, (5.11) with the output matrix:  1 0  0 1 Ω =   0 0 0 0

0 0 1 0

0 0 0 0

0 0 0 1

 0 0   . 0  0

(5.12)

Using the same reasoning as given for angle θ2 in the R1 case, angles θ1 and θ2 are affected by adjusting the pressure in respectively muscles C and D and muscles A and B in an opposing manner. To control the pressure in muscles A and C , only the subsequent inputs are used. To control the outputs y, four new inputs are generated by applying the input matrix Ψ according to: £ ¤T zθ1 zθ2 zpA zpC . (5.13) z = u

=

Ψz 

0  0 Ψ =   −1 1

(5.14)

−1 1 0 0

1 0 0 0



0 0   . 1  0

(5.15)

In Appendix H.5, the conversion of R2 to the plant P 2 is visualized. To control the outputs y, four controllers are required: Cθ1 , Cθ2 , CpA and CpC respectively. The four inputs z and four outputs y of P 2 are related to each other by: Pθ1 =

θ1 zθ 1

, Pθ2 =

θ2 zθ 2

,

PpA =

pA zpA

, PpC =

pC zpC

.

(5.16)

It should be noted that the inputs pA and θ2 are not decoupled; using Figure 5.1, it can be seen that a change in pA also affects θ2 and vice-versa. The same reasoning is also valid for pC and θ1 . As a result, interference of the two control loops may cause instability. However, this is not expected, as the two control loops work opposing to each other for muscle A ; i.e. the control signals zθ2 and zpA extinguish each other. A secondary effect of this is that the settling time of θ2 may increase. The control loops of pA and pC are decoupled, the coupling between the two angular control loops for θ1 and θ2 , is discussed in the following section. 69

Pθ11 =

10

θ1 zθ1

0 −10

−20

−20

−30

−30

−40

−40

−50

−50

−60

−60

|Pθ | [rad V−1 dB]

0

−10

−70 −80

−70 −1

0

10

10

Pθ21 =

10

−80

1

10

θ1 zθ2

−20

−20

−30

−30

−40

−40

−50

−50

−60

−60

|Pθ | [rad V−1 dB]

0 −10

−70

−70 0

10

10

0

−80

1

10

10

1

10

θ2 zθ2

Pθ22 =

10

0

−1

−1

10

−10

−80

θ2 zθ1

Pθ12 =

10

p¯C = 0.4 bar p¯C = 0.6 bar p¯C = 0.8 bar −1

10

f [Hz]

0

10

1

10

f [Hz]

Figure 5.8: The four frequency response functions of P θ for different values of p¯C . The FRF’s are experimentally determined around θ¯1 = θ¯2 = 0 deg. The pressure in muscle A is kept at p¯A = 0.4 bar in closed-loop.

5.3.2 MIMO system identification It is expected that the dynamics of the robotic arm are coupled, i.e. if one angle changes, the dynamics of the robotic arm is changed yielding a new equilibrium for both angles. In case the coupled dynamics are not considered in the design of the two angular control loops, the robotic arm might become unstable due to the fact that the two control loops might interfere. For this reason, a closed-loop sensitivity identification is used to determine the frequency response functions of the angles of the 2DoF robotic arm. This is explained in Appendix H.6, using the theory as given in [Skogestad, 1996; Tousain, 2006]. The results of this analysis are given in Figure 5.8 for three different equilibria. The direct FRF’s are given by Pθ11 and Pθ22 , the coupling from the θ2 loop to the θ1 loop is given by Pθ12 and vice-versa by Pθ21 . Angle θ1 gives a response of approximate |Pθ12 | ≈ −30 dB to an excitation of θ2 . The response of |Pθ21 | ≈ −40 dB. The magnitude of Pθ11 and Pθ22 are given by respectively −10 and −20 dB. Using these approximated magnitudes, the influence of the coupling on the input-output behavior is considered for joint i, taking into account the coupling with joint j. This yields that the contribution of the coupling on the input-output behavior is negligible: ¡

Cθi Pθii

1 + Cθi Pθii − Cθi Cθj Pθji Pθij

¢¡

1 + Cθj Pθjj

¢−1 ≈

Cθi Pθii . (5.17) 1 + Cθi Pθii

The low coupling is mainly due to Pθ21 , joint 2 hardly responses to joint 1 excitation. In Appendix H.7, the minor coupling is proven, using the interaction coefficient. 70

Feedback control of the robotic arm

5.3.3 Loop shaping of the controller The 2DoF robotic arm requires four controllers, which are all based on loop shaping. In this section the designed controllers and the performance of the closed-loop system are discussed. Angular controllers The angular frequency response functions are given in Figure 5.8. Although the coupling is of not much importance, it is considered in the loop shaping of Cθ1 and Cθ2 . The structure of the angular controllers is similar to angular controller as discussed in Section 5.2.1: a gain is combined with a lead-lag, notch and secondorder low-pass filter. Compared to the DoF case, the controller for θ1 is almost similar (see Appendix H.4). The controller design for Cθ2 is different from the DoF case as is given in Section 5.2.1. The reason for this is to check whether it is possible to compensate for the differentiating behavior of the pneumatics. The lead-lag filter is applied over a larger frequency range. This requires the gain to be lower. The result of the loop shaping of the angular controllers is given in Figures 5.9 and 5.10 for respectively Cθ1 and Cθ2 . The bandwidth of the θ1 loop is 1 Hz, the bandwidth of the θ2 loop approximately 3.5 Hz1 . The controller parameters are given in respectively Tables 5.4 and 5.5. In the lower part of Figures 5.9 and 5.10, the input sensitivity S is given. In Appendix H.8, the accompanying Nichols diagrams are shown. As can be seen, the phase, amplitude and sensitivity margins are met for both controllers.

gain p

16

lead-lag f1 f2

notch

0.2 8

fn β1 β2

1.4 0.1 4

low-pass fl βl

6 0.8

Table 5.4: The Cθ1 controller parameters for joint 1 of the 2DoF robotic arm. The controller consists of a gain with a notch, lead-lag and a second-order low-pass filter. All frequencies f are given in Hz.

gain p

5.4

lead-lag f1 f2

0.2 18

notch fn β1 β2

6.1 0.06 2

low-pass fl βl

12 0.8

Table 5.5: The Cθ2 controller parameters for joint 2 of the 2DoF robotic arm. The controller consists of a gain with a notch, lead-lag and a second-order low-pass filter. All frequencies f are given in Hz. 1

It might seem that the bandwidth of joint 2 in the 2DoF case is higher then in the DoF case. This is not true, only three FRF’s around θ¯2 = 0 deg are considered. In case the FRF’s around θ¯ 6= 0 are also considered, the bandwidth of the 2DoF θ2 loop is lower. Also, the eigenfrequency of θ2 in the 2DoF case, is slightly lower compared to the DoF case.

71

40 20

|H| [dB]

0 −20 −40

Cθ1 [V rad−1 ] HOL [–] Pθ2 [rad V−1 ]

−60 −80

−1

0

10

10

phase(H) [deg]

180

90

0

−90

−180 −1

0

10

10

20

|S| [dB]

10 0

−10 −20 −30

−1

0

10

10

f [Hz] Figure 5.9: The two upper figures show the magnitude and phase of controller Cθ1 , plant Pθ1 and the open-loop HOL = Cθ1 Pθ1 . The phase margin at −145 dB is indicated by the pink line in the second figure. The lower image shows the sensitivity S = (1 + Cθ1 Pθ1 )−1 of the control loop, the sensitivity margin of 6 dB is indicated by the pink line. The bandwidth is approximately 1 Hz, as indicated by the vertical pink line.

72

Feedback control of the robotic arm

40

|H| [dB]

20 0

−20

Cθ2 [V bar−1 ] HOL [–] Pθ2 [bar V−1 ]

−40 −60

−1

0

10

10

1

10

phase(H) [deg]

180

90

0

−90

−180 −1

0

10

10

1

10

20

|S| [dB]

10 0

−10 −20 −30

−1

0

10

10

f [Hz] Figure 5.10: The two upper figures show the magnitude and phase of controller Cθ2 , plant Pθ2 and the open-loop HOL = Cθ2 Pθ2 . The phase margin at −145 dB is indicated by the pink line in the second figure. The lower image shows the sensitivity S = (1 + Cθ2 Pθ2 )−1 of the control loop, the sensitivity margin of 6 dB is indicated by the pink line. The bandwidth is approximately 1 Hz, as indicated by the vertical pink line.

73

1

10

Pressure controllers Both pressure controllers CpA and CpC are equal: a gain with a first-order low-pass filter. The gain is lower compared to the Cp controller as defined for the DoF robotic arm and the low-pass filter frequency is located at 4 Hz. The reason for this is simple; to lower the bandwidth of the pressure, yielding a similar bandwidth as the θ1 loop. The settings of the pressure controllers are given in Table 5.6. gain p

low-pass 5

fl

4

Table 5.6: The pressure controllers CpA and CpC are identical. The controller consists of a gain and a first-order low-pass filter. The low-pass frequency fl is given in Hz.

Performance The performance of the 2DoF system is indicated by the bandwidth (BW), settling time ts , rise time tr , overshoot Mp and static error ess , see Table 5.7. The bandwidth is taken from the FRF graphs, while the four time domain indicators are taken from the step responses measurement as given in Figures 5.11 and 5.12. In Figure 5.11, steps are given in the desired trajectory of θ1 and pC , while the setpoint for θ2 and pA are kept constant. The low coupling between the two joints is confirmed by θ2 and pA , which show no sign of such coupling. The behavior of θ1 is very undamped, which results in a long settling time. The coupling between θ1 and pC is indeed present, as was predicted in Section 5.3.1. The coupling is clearly visible at approximately 20 seconds. Angle θ1 is heavily excited by a step in pC . The coupling increases the settling time of pC . In Figure 5.12, similar steps are set to the reference of θ2 and pA , while the desired θ1 and pC are kept constant. Again, the two joints appear to be decoupled, while θ2 and pA are coupled. Most remarkable, however, is the apparently damped behavior of θ2 . Angle θ1 does not show this behavior in Figure 5.11. The only difference between both joints is the bandwidth: the bandwidths of θ1 and pC are equal, the bandwidths of θ2 and pA are not. The damped behavior is a result of pressure controller CpA and the used control strategy. If θ2 changes, both pA and pB are affected. The change of pA is opposed by CpA ; i.e. the final rotation is mainly the result of the change in pB . For this reason, it output

BW [Hz]

ts [s]

tr [s]

Mp

ess

θ1 pC θ2 pA

1 1 3 1

4.7 3.3 2.8 0.7

0.4 0.3 2.1 0.2

7 deg − − −

0.30 deg 0.01 bar 0.50 deg 0.02 bar

Table 5.7: The performance indicators for the 2DoF robotic arm as derived from the step excitations, given in Figures 5.11 and 5.12. No overshoot of θ2 can be given.

74

Feedback control of the robotic arm

is believed that the coupling between pA and θ2 will not result in stability problems. The two opposing controllers Cθ2 and CpA are counterbalanced due to the input matrix Ψ. This results in a limited energy added to the system, which appears to be damped behavior of the joint. The reason that this occurs only in joint 2 and not in joint 1 is the lower bandwidth of pA compared to θ2 . When the angle error becomes small, the motion becomes low-frequent; a low-frequent pressure change is strongly opposed by CpA . It can be concluded that by lowering the bandwidth of the pressure loop, damping is added to the robotic arm. This will only work for the 2DoF control strategy; in the DoF strategy, the angle and pressure are decoupled.

5.4 Feedback control versus reinforcement learning In [Maas, 2005], a reinforcement learning (RL) algorithm is proposed to control the robotic arm. This controller is capable of tracking a constant setpoint. An experimental model of the robotic arm is constructed in the form of a neural network. In every experiment, the neural network is updated; so the neural network learns the behavior of the robotic arm. The neural network, modelling the robotic arm, uses the three most recent states as input. Nine different actions are defined as controller output. The controller uses a second neural network to estimate the implication of each of the nine actions. Based on the nine estimate results and the endpoint of the p2p motion, one action is selected and performed. It can be concluded that the RL controller is a feedforward controller. The current states are not used to calculate the tracking error, but to make an estimate of the implication of a certain action. The advantages and disadvantages are listed below: Advantages • No dynamic system knowledge is required to design the controller; • The RL controller is robust to small and slow system changes. Disadvantages • The controller is a black box, i.e. its behavior can not be predicted; • The controller must be ”learned” on the robotic arm, which requires hundreds of trials. When comparing the results of RL control and feedback control, no standard performance indicators (i.e. bandwidth, robustness margins, settling time, overshoot, etc) can be defined. However, three conclusions are given when comparing the implementation and result of both controllers: • The current RL controller is not suitable for trajectory tracking. In case of setpoint tracking, in 90% of the cases the target is tracked. In the other 10%, the trial is aborted after 100 actions. The feedback controller tracks a trajectory correctly. • In an equilibrium point, the RL controller suffers severely from limit cycling, while the feedback controller is stable due to the dead-zone in the PWM. 75

θ1 [deg]

40 20 0

reference experiment

−20 −40

0

5

10

15

20

25

30

0

5

10

15

20

25

30

0

5

10

15

20

25

30

0

5

10

15

20

25

30

pC [bar]

1

0.5

0

θ2 [deg]

40 20 0 −20 −40

pA [bar]

1

0.5

0

time [s] Figure 5.11: Step response or the 2DoF controlled robotic arm. The angle step of 15 deg and the pressure step of 0.25 bar are used to derive the performance indicators in the two upper rows in Table 5.7. Only θ1 and pC are excited, θ2 and pA are kept constant by the controllers.

76

Feedback control of the robotic arm

θ1 [deg]

40 20 0

reference experiment

−20 −40

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

16

18

20

pC [bar]

1

0.5

0

θ2 [deg]

40 20 0 −20 −40

pA [bar]

1

0.5

0

time [s] Figure 5.12: Step response or the 2DoF controlled robotic arm. The angle step of 15 deg and the pressure step of 0.25 bar are used to derive the performance indicators in the two lower rows in Table 5.7. Only θ2 and pA are excited, θ1 and pC are kept constant by the controllers.

77

• In the feedback controller, a notch filter is essential. The RL controller must also generate behavior similar to a notch filter; which seems unlikely as the RL controller appears to function like a feedforward besides; only three time samples are considered what is insufficient for a notch filter. The advantage of a black box approach of the RL controller in fact also counts for the loopshaping approach. In depth dynamic system knowledge is not required, overall frequency-domain system identification is sufficient. This is less time consuming than the hundreds of trials required to constructed the model in the RL approach. The tuning of controllers by loopshaping is simple and effective, even in case of a nonlinear system with low bandwidth (as proven on the robotic arm). Using linear control, performance guarantees can be given and the robustness to certain system changes can be estimated (using the sensitivity and the Nyquist stability criterion). Based on the advantages and disadvantages, it can be concluded that the control strategy as proposed in this thesis is preferred.

5.5 Discussion Both the DoF and 2DoF robotic arm are stabilized using linear feedback controllers. The robot is able to follow desired trajectories, although the bandwidth is limited. In this discussion some remarks on the two control strategies and the performance are given. The proposed control strategy for the DoF robotic arm results in two decoupled outputs. The proposed control strategy for the 2DoF robotic arm gives two coupled outputs for each joint. Both strategies stabilize the system correctly. The DoF strategy results in a maximum performance of the robotic arm, while the 2DoF allows for adding damping by adjusting the pressure controller. Obviously, the DoF control strategy can also be used on the 2DoF system by using a different dSPACE system with sufficient inputs, so all sensor signals can be used. The bandwidth of the robotic arm is 3 Hz at best. Considering the system, this is good. The main performance limiting factors are the working pressure and the low eigenfrequency of the robotic arm. The nonlinearities, introduced by the McKibben muscles, are handled correctly by the linear feedback controllers and seem no limitation. The robotic arm behaves like an integrator at low-frequencies. For this reason, it is assumed that the static error is a result of the of the dead-zone τ , as introduced in the PWM. In order to reduce the static error, experiments are performed with a smaller dead-zone. This results in nervous behavior and limit-cycling, which proves that the dead-zone is essential in the control loop. In the following chapter, the final conclusions and recommendations of this work are presented.

78

6 Conclusions and recommendations To investigate the applicability of pneumatic artificial muscles as actuators for domestic robots, research is performed on the modeling and control of a robotic arm actuated by McKibben muscles (a type of pneumatic artificial muscles). To this extent, three questions are formulated in the problem statement: A predictive model for the two DoF robotic arm should be built, based on physical principles and dedicated experimental identification, where special attention should be given to the modelling of the artificial muscles. This model should be validated by experiments and used to construct linear feedback controllers, using a classical feedback approach, allowing for both stabilization and tracking tasks. Finally, these controllers should be implemented and evaluated on the experimental setup. A robotic arm model is proposed, consisting of a pneumatic model (describing the binary valves, flows and pressures), a newly developed model for the McKibben muscles and a multibody model of the robotic arm dynamics. The pneumatic model, the muscle stiffness and the active muscle behavior are identified using dedicated experiments; the experimental identification of the robotic arm dynamics is based on frequency-domain experiments. The model is validated by comparing the equilibria and the dynamic behavior in both the time- and frequency domain of the model to those of the experimental setup. A control strategy is proposed for the decoupling of two McKibben muscles in an antagonist setup. Linear controllers are designed using loop-shaping techniques and implemented on the robotic arm. This chapter gives the conclusions divided in three parts. Firstly, conclusions are given with respect to the four model parts: valves, pneumatics, muscles and robot dynamics. Secondly, conclusions on the entire robotic arm model are drawn and, finally, conclusions on the control of the robotic arm are given. This chapter is ended by recommendations for future research.

Valves, pneumatics, muscles and robot dynamics By applying pulse width modulation to the control signals, the binary behavior of the valves becomes continuous by approximation. Limit cycling of the binary valves is prevented by introducing a dead zone in the PWM. The phase shift due to the

PWM is −10 degrees at 5 Hz. As the maximum bandwidth is 4 Hz, the binary valves are no restriction for the performance of the robotic arm. The pneumatic model uses straightforward relations to describe the pressures and the flows. Experimental identification on a test setup shows that the pneumatic model is accurate. The robotic arm identification has shown that the supply pressure is the performance limiting factor in the robotic arm. Existing McKibben muscle models appear to be inadequate for describing the observed behavior, as the interaction of the muscle with its mechanical environment is generally not taken into account. This yields a wrong description of the McKibben muscle in the pressure application range below 2 bar. A new model is proposed, which takes into account the variation in stiffness, damping and volume, depending on the muscle length and pressure. The output of the model is a force, which is a result of the interaction with the environment (in terms of the flow and muscle length). Three conclusions concerning the muscle model can be given: • The structure of the muscle model allows the muscle to be interpreted as a force or a position actuator, depending on the load of the muscle. The robotic arm identification shows that the two antagonist muscles act like a force actuator; this is confirmed by the simulations of the robotic arm model. • The stiffness of a McKibben muscle is identified on a tensile bench; these dynamic experiments are used to fit the parameters of the stiffness relation of the muscle model. The outcome of this fit yields a stiffness relation which is not in agreement with the stiffness following from the frequency-domain identification experiments on the robotic arm. Both the muscle identification and the robotic arm identification are considered to be correct. The mismatch is caused by errors in the fitted stiffness relation; these errors are a result of differences between the application range of the muscles in the robotic arm (pressure: 0.3 − 1.0 bar; muscle length: 168 − 213 mm) and the range used for the muscle identification (0.5 − 1.4 bar and 150 − 200 mm). • The muscle model considers the muscle to be barrel shaped. Simulations in which a cylindrical shaped muscle is assumed have shown that the muscle shape is of minor importance on the muscle behavior; namely, the relative volume change is similar in both cases. A multibody model is used to derive the equations of motion for the DoF and 2DoF robotic arm. This model is considered to be accurate. The length, mass and inertia parameters are determined using the test setup and Unigraphics (to access the geometrical data). The friction in the robot joints is not taken into account; this is compensated for in the muscle model in the form of additional damping.

Robotic arm The DoF robotic arm model is given by a fifth-order nonlinear differential equation; the 2DoF model is represented by a ninth-order nonlinear differential equation. The modeling has given an increased insight in the functioning of the robotic arm and the McKibben muscles. Three conclusions are given on the validation of the robotic arm model: • The main nonlinear contribution to the robotic arm behavior is given by the muscle stiffness. 80

Conclusions and recommendations

• Due to errors in the fitted muscle stiffness relation (which are stated above), the robotic arm model is only predictive if the pressure difference between two antagonist muscles is below 0.5 bar. • Three types of energy dissipation are distinguished: damping in the muscles, friction in the muscle braid and friction in the robot joints. In literature, the friction in the muscle is considered to be one of the dominant factors in the muscle behavior. This is opposed by the identification experiments performed in this work, which show that the energy dissipation in the muscle is not significant for its behavior. The dimensionless damping ratio of the robotic arm is 0.05 ± 0.02 and, according to the muscle identification, approximately half the dissipation occurs in the muscles).

Control strategy A control strategy is proposed, which yields two decoupled inputs. The pressure difference between two muscles of a muscle pair is used for generating a torque, while the cumulative pressure is used as an additional input; these two inputs are decoupled. The cumulative pressure is a measure for the robotic arm stiffness. The attained decoupling is not perfect, but satisfactory. This allows for the application of SISO loop-shaping techniques to generate controllers and the following related conclusions can be drawn: • A notch filter in the angular controller is indispensable to stabilize the robotic arm. The eigenfrequencies of the robotic arm are position dependant due to the nonlinear stiffness of the muscles; the range of eigenfrequencies is compensated for with one notch filter. • A bandwidth of approximately 4 Hz is achieved. The robotic arm is experimentally proven to be stable around each equilibrium. Due to hardware restrictions, the control strategy used for the 2DoF robotic arm is slightly different. The two inputs to each joint are chosen differently, as different outputs are available; these two inputs are not decoupled. This does not lead to stability problems, but the bandwidth drops to approximately 1 Hz; this is due to the interaction of two outputs of one joint. The two joints are completely decoupled. In most literature, complex control strategies are used to control systems actuated by McKibben muscles. However, this work proves that good results are possible by using loop shaping techniques to design a linear controller. The presented control strategy has proven to give better results than an AI approach, see [Maas, 2005].

Recommendations Two recommendations concerning the robotic arm are given: • It is advisable to apply a McKibben muscles with pre-stress. This yields a larger working area as the tube touches the braid at a lower pressure. Stretching the muscle with 10 to 15 mm is sufficient. • It is recommended to use a higher working pressure, while keeping the maximum muscle pressure at approximately 1 bar. This would increase the performance while the muscles remain compliant. 81

To improve the robotic arm model; future research is desirable: • In the McKibben muscle model, the nominal length relation and the stiffness relation are derived by two different approaches. These two relations must match, as in an equilibrium the force balance, determining the nominal length; gives the stiffness of the muscle. The physical relation behind the stiffness and nominal length requires further research. • The dedicated muscle stiffness measurements and the frequency-domain identification experiments do not agree, with respect to the muscle stiffness. It is advised to repeat the muscle identification for an application range (muscle lengths and pressures), similar to the application range of the muscles in the robotic arm. • Experiments show that stretching a muscle while keeping the air mass constant results in a significant pressure increase; this is caused by a decrease in volume. The muscle model hardly shows this behavior. This indicates that the change in muscle volume not only depends on the muscle length, but also on the pressure (or the nominal muscle length). • Measurements on the muscle volume for a range of lengths and pressures would be desirable to define a better volumetric behavior. • The additional value of the muscle model, compared to the existing Chou– Hannaford model, requires additional research. Both models are predictive in a different application range; on the overlapping application range, they should give similar results. Considering the control strategy for McKibben muscle actuated robots, two recommendations for future research are given: • The decoupling between the joint angle and the cumulative pressure is not verified. An experimental MIMO analysis would allow for a value judgment of the decoupling strategy, based on the interaction coefficient. • The linear controller is performing well, but there is clearly room for improvement. The robotic arm identification can be used to generate a nonlinear model. Based on such a model, techniques such as input-output linearization, feedforward design and gain scheduling can be used. It would be interesting to perform research on the increase in performance by using more advanced control techniques. Finally, some remarks on the application of McKibben muscles as robot actuators in a domestic environment will be made. It is questionable whether McKibben muscles are desirable for this application area; namely, they suffer from wear and friction, no two muscles are identical, their dynamic behavior is complex and pressurized air is required. Electric actuators do not suffer from these disadvantages. The compliance of McKibben muscles might be an advantage in a domestic environment, but also electric actuators can be made compliant (for example by introducing a compliant element between the actuator and end-effector). Whether a McKibben muscle or electric actuated domestic robot is used, numerous sensors are required to detect possible unsafe situations. If electric actuators are applied, it is most important to use a control strategy which does not focus on performance, but on a human friendly behavior, i.e. low bandwidth and no strong feedback. In any case (electric of McKibben muscle actuated), the trajectories must be smooth, quiet and friendly. 82

A Specifications of the robotic arm parts

This appendix presents additional information on the robotic arm. First, some robotic arm images are shown. In Appendices A.2 to A.4, the specifications of the binary valves, pressure sensors and potentiometers are given.

A.1 Additional photos of the robotic arm In Figure A.1, it is shown that each McKibben muscle has two air hose connectors; one connector is used as air in and outlet while the other connector is used to measure the pressure in the muscle. The robotic arm can serve as a DoF and as a 2DoF system. By fixing revolute joint 1 (indicated in Figure 2.1), the robotic arm becomes a DoF system, consisting of revolute joint 2 and actuated by muscles A and B . The DoF system is shown in Figure A.1, the 2DoF system in Figure A.2. A close up of the valves and sensors, required for the revolute joint 2 , is shown in Figure A.3. All revolute joints and muscle hinges are mounted in bearings, joint forces out of the working plane of the robotic arm are compensated in the bearings. Each potentiometers has an internal bearing; it is mounted to the revolute joint through a Koster coupling to avoid over-fixation of the sensor, see Figure A.4. The binary valves in the robotic arm are addressed by the controller. This controller is implemented in M ATLABr S IMULINKr and uploaded to a dSPACEr system. Inbetween the dSPACE and the robotic arm, a switching board is mounted, see Figure A.5. The switching board matches the hardware interfaces of the dSPACE, the sensors and the valves; it has three additional functions: • The communication between the dSPACE board and the valves needs amplification. This is realized using eight FET’s. • An external power source is implemented to power the switching board itself, the FET’s and the pressure sensors. • The used dSPACE board has four inputs; there are six sensor signals available. Two jumpers determine which four inputs are directed towards the dSPACE. The electrical scheme of the switching board is given in Appendix A.5.

Figure A.1: The part of the robotic arm making up the DoF system. Revolute joint 1 is fixed by a beam which is visible in Figure A.3.

Figure A.2: The 2DoF robotic arm in its initial position.

Figure A.3: The potmeter, 4 valves and 2 pressure sensors of the second DoF.

84

Specifications of the robotic arm parts

Figure A.4: Both potentiometers are mounted to the axis using a so called Koster coupling. This coupling fixes DoF (rotation) while it allows motion in the other 5DoF. This is the revolute joint 1 potentiometer.

Figure A.5: The switching board as implemented inbetween the d SPACE system and the robotic arm. In Appendix A.5, the electrical scheme is given.

85

A.2 Binary valves

Festo MHE2 solenoid fast-switching valve

Electrical data Operating voltage Connection type Power consumption With fast-switching electronics Without fast-switching electronics

[V DC]

5 ±10%, 12 ±10% or 24 ±10% Plug vanes or moulded-in cable

[W] [W]

5 for 3 ms approx. (pull 1 A), then 2.88 W 2.88

Protection class to EN 60 529 With moulded-in cable With plug socket with cable KMYZ-3 With plug socket with cable KMYZ-3 and plug M8 With plug socket with cable KMYZ-4

IP65 IP65 IP65 IP40

Response times and switching frequencies With fast-switching electronics Response time on/off Switching frequency

[ms] [Hz]

2 ±10% (from 1 Hz) 150

Without fast-switching electronics Response time on/off Switching frequency

[ms] [Hz]

7 ±10% (from 1 Hz) 50

I [mA]

Current curve

1 2 3 4

t [ms] Current in the coil Current in the supply line

Capacitor charging Controlled coil current 1 A Drop to holding current Controlled holding current 0.5 A

Figure A.6: The specification of the binary valves, [Festo, 2006].

86

Specifications of the robotic arm parts

Valve response time measurements The valve response time is checked by an experiment. The flow response is determined by measuring the pressure directly behind the valve. In Figure A.6, the current curve in the supply line is given for opening the valve. A similar curve is measured, as is shown in Figure A.7. The pressure starts increasing after 2 ms, so the flow response time for opening the valve is also 2 ms. Several measurements have been performed, yielding similar results. The flow response time for closing the valve, is given in Figure A.8. The current shows a drop from 0.1 A to 0.0 A at the moment of switching. The pressure starts increasing after 3 ms, so the flow response time for closing the valve is given by approximately 3 ms.

2.5

I [A] p [0.1 bar]

2 1.5 1 0.5 0 −0.5

0

1

2

3

4

5

t [ms]

6

7

8

9

10

Figure A.7: Current and pressure behavior for switching on the valve. The current in this figure shows great resemblance with the current as shown in the valve specification. On the x-axis, both the current in Ampere and the pressure in bar is shown. The pressure in the plot is multiplied by 10, to make the pressure increase it beter visual, this is also indicated in the legend.

0.15

I [A] p [bar] 0.1

0.05

0

−0.05

0

1

2

3

4

5

t [ms]

6

7

8

9

Figure A.8: Current and pressure behavior for switching off the valve. On the x-axis, both the current in Ampere and the pressure in bar is shown. Remark, the scaling differs from the scaling in Figure A.7.

87

10

A.3 Festo SDE1 pressure sensor

Function1)

-P-

Voltage 15 … 30 V DC

-L1)

Pressure –1 … +10 bar

For example with 1 switch output PNP and 0 … 10 V analogue

-QGeneral technical data Pressure measuring range

[bar]

0 … –1

Mechanical Method of measurement Pneumatic connection Measured variable

1) 2)

–1 … +1

0…2

0…6

0 … 10

Piezoresistive pressure sensor with display R
, R, G
or QS-4 Relative pressure Differential pressure1) ±2% of the measuring range final value 0.3% Plug M8x1 or M12x1, round design to EN 60 947-5-2 Front panel mounting or on service unit, H-rail or adapter plate Any2)

Accuracy Reproducibility of switching point Electrical connection Type of mounting Mounting position Electrical Operating voltage range Max. output current Protection against short circuit Protection against polarity reversal Switch output CE symbol

Temperature range 0 … 50°C

[V DC] [mA]

15 … 30 150 Pulsed For all electrical connections PNP or NPN 89/336/EEC (EMC)

Variants with push-in fitting QS-4 The collection of condensate in the sensor should be prevented.

Dimensions – H-rail, wall or surface mounting Pneumatic connection G


1 Plug M8x1 or M12x1 to EN 60 947-5-2 2 LCD display 3 Pneumatic connection G
4 Connection socket, straight 5 Connection socket, angled 6 Adapter plate for wall mounting 7 Centre with H-rail mounting Type SDE1-…-W18-…-M8 SDE1-…-H18-…-M8 SDE1-…-W18-…-M12 SDE1-…-H18-…-M12

B1

D1

D2

H1

H2

L1

L2

L3

L4

L5

32.3

G


M8

35.2

3.5

78

70

107

89

33

32.3

G


M12

35.2

3.5

87

70

125

104

33

88

Specifications of the robotic arm parts

A.4 Novotechnik SP2500 potentiometer

89

+24V

+24V

A.5 Electrical scheme of the switching board

BYW56

R5

U2 gain

C3

B1 Ye

A0

R15 10K

BYW56

R10

1K8

Edit

OUT Festo MHE2-MS1H-3/2

Gr

BST70A

10K

1.12

Black

B3

V4

35

Black

2

2K2

3

0.83V/Bar

+24V

R1

1

+24V

R2

BYW56 Port C4

B

V1

Port A0

R6

5V1

5K6

+24V

4

0-5V

A1

Potentiometer

C4

Black

4

Black

OUT Festo MHE2-MS1H-3/2

Gr BST70A

R16 10K

?K

3

B4

V6

36

5V6

+24V

Port A1

V5 10K

5 6 7

White Brown Green

2K2

A

+24V

Festo SDE-D10-G2-W18-L-PU-M8

10K

DLP-245PB

1 2 3 4

Port C3

X3

V3

USB / PIC16877

X2 Brown White Black Blue

R11 For schematics and source: www.dlpdesign.com/dnlda

+24V

8 9 10 11

BYW56 Port C5

R7

gain

R3

5

0.83V/Bar

B2 Ye BYW56

C5

V2

Port A2

BST70A

R12

JP3 12 13 14

Spare

IN

Bottom View

Port C6

C6

X4

CN1

Festo MHE2-MS1H-3/2

2K2

10K

R13

13 15 17 40 18

+5V

19

+5V

BYW56 GND

R9

GND

10K

10

GND GND B0

GND

U1

C4

C5 C0

C7

+5V

1 24 +Vin

+

OUTPUT

A3

C3

E0

22uF

C1

A5

10 +15V Common 15

INPUT

33 34 39 6 2 8 +24V +5V

1 2 3 4 5 6 7 8 9 10

1 2 3 4 5 6 7 8 9 10

+24Va +5Va

+24Va

100nF 12 13 -Vin

32

BYW56 Port C0

10K +24Va

Black

2

Black

2K2 +24Va

1 2 3 4

10K gain 0.83V/Bar

BYW56 Port C1

Ye

10K

1.12

BYW56

3

Black

4

Black

5V1 White Brown Green

5 6 7

3K65

+24Va

0-5V

2K2

BST70A

Port A5

10K

+24Va

Potentiometer (2.1 - 4K)

Port E0

5V6

BYW56

Port C2

10K 8 9 10 11

BST70A

Black

0.83V/Bar

10K

10K

IN Festo MHE2-MS1H-3/

+24Va

BYW56 1K8

Edit

Black

6 Gr

gain

Ye

1.12

5

2K2

Brown White Black Blue

OUT Festo MHE2-MS1H-3/2

Gr

B

+24Va

A

10K

1K8

Edit

OUT Festo MHE2-MS1H-3/2

+24Va

Brown White Black Blue

1

Gr BST70A

+24Va

2AT

C2 +5V

BYW56

Port C7

A

B

5V1 10K

Port A3

Black

8

Black

IN Festo MHE2-MS1H-3/2

Gr BST70A

7

2K2

F1

C1

14 +Vout 11

+24Va

2

Festo SDE-D10-G2-W18-L-PU-M8

1 24V supply

Festo SDE-D10-G2-W18-L-PU-M8

X1

+24V

100nF 100nF

C2

- Circuit board in Eddystone box - U1 / U2 in IC-sockets - Preferably use SMD resistors on the track side of the board - X1 - X4 are screw terminals (orange) - F1 is a glass fuse

R14

+24V +

22uF

1 1

Remarks:

Gr BST70A

R19

EXTVCC

Spare

B7

V12

31

VCC-IO

10K

5V / 24V Supply

JP3 5

V11 9

GND

DC/DC Converter 305S 24FR

Microchip PIC16F877

40

2K2

10

DLP-245PB

1 PORTVCC 2 USBDM 3 USBDP 4 GND

+24V

1

FTDI FT245BM

1 MCLR/n 2 RB6 3 SWVCC 4 GND 5 RB7

IN

20

CN1

CN1 USB Conn.

+24V

R18

USB

Black

Gr BST70A

USB

Black

8

B6

V10

38

PC Standard USB Cable

7

USB 21

U2 Upper view

V9 10K

A4

D

JP3 Programming Header

BYW56

R8 7

S

Festo MHE2-MS1H-3/2

2K2 +24V

10K

5V1

Black

G

R17

B

Black

6 Gr

BST70A

R4

5

B5

V8

37

1K8

+5V

Edit

A

10K A2

10K

1.12

V7

+24V

Festo SDE-D10-G2-W18-L-PU-M8

Brown White Black Blue

10K +24Va

90

DLP-245PB Page-1 Cor van der Klooster 20-8-2006 switchingboard.vsd

B McKibben muscles properties

This appendix presents some additional information on McKibben muscles. In Section B.1, the major advantages and disadvantages of McKibben muscles are given. The results of a literature review is presented in Section B.2, followed by an example on the application of a McKibben muscle as a robotic actuator. Section B.4 discusses the contribution of the rubber tube to the muscle behavior.

B.1 Advantages and disadvantages An advantage of all types of pneumatic artificial muscles (PAM’s) is the energy to mass ratio. Compared to conventional actuators (pneumatic, as well as hydraulic and electric), PAM’s excel in this respect [Plettenburg, 2005]. In applications in which the total amount of mass is highly critical, PAM’s might be an alternative. A second advantage of pneumatic artificial muscles is that they can be applied in a very dusty, humid or wet environment. This is contradictory to conventual (electrical, mechanical or pneumatic) actuators, which must be kept clean and dry. McKibben muscles show great similarity with biological muscles. In [Klute, 1999], a comparison, based on experiments, between pneumatic and biological muscles is made. In [Klute, 2002], the McKibben muscles are compared to the Hill model, a model describing the behavior of biological muscles. The similarity is based upon the length – stiffness relation of the muscles. The length – damping relation is different for McKibben muscles and biological muscles; the damping in McKibben muscles is much lower. The compliance of the McKibben muscle is an advantage from domestic application point of view, but a disadvantage from control theory point of view. A model describing the actuator behavior is required in order to have the actuator generating a required output force, defined by the control signal. The major disadvantage of McKibben muscles, is that the relations defining its behavior are nonlinear. Both the stiffness and damping depend in a nonlinear manner on the actual muscle length and the pressure. Besides, McKibben muscles

suffer from hysteresis due to friction between the tube and the braid and between the nylon fibres of the braid. In [Chou, 1994], it is shown that the friction in the braid is dominant.

B.2 Literature on McKibben muscles In this section, the results of a small literature review are presented. The consulted literature comprehends the modelling of McKibben muscle behavior and the control of systems actuated by McKibben muscles.

B.2.1 McKibben muscle models The muscle behavior depends largely on the surrounding dynamics and the pressure supply. In every literature source, some kind of dynamic load is taken into account, in most cases a known mass is lifted. In [Chou, 1996; van Ham, 2003], the pneumatics are taken into account in the modelling. In [Colbrunn, 2001], a muscle setup is modeled, including a complex relation for the valve. In most cases however, the pressure is considered as an input to the system. Three models describing the behavior of McKibben muscles can be distinguished, see [Chou, 1994; Tondu, 1997; Colbrunn, 2001]. These three models do not use exactly the same parameters, but can be transformed into each other [Schröder, 2003]. The most referred model is the Chou–Hannaford model [Chou, 1994]: Fm =

π d20 pm 4

µ

3 1 (1 − ε)2 − 2 tan ξ0 sin2 ξ0

¶

,

(B.1)

with the muscle force Fm , the resting diameter d0 , the muscle pressure pm , the rest braid angle ξ0 and the elongation ε. In Appendix E.1, the complete derivation of the Chou–Hannaford model is given. In [Chou, 1996; Colbrunn, 2001] the relation in (B.1) is simplified (by reconsidering the individual terms) to: Fm = kg (pm − pth ) (h − hmin ) + kp (h − hr ) ;

(B.2)

with kg a supposed constant gas stiffness, kp an elastic constant for the shearing force, the actual length h, the minimal length hmin , the resting length hr and a threshold pressure pth at which the tube and braid touch. According to the definition of minimal length, the force due to the air in the muscle is zero if the actual and the minimal length are equal. Equation (B.1) has the two variables pm and ε, while (B.2) has the two variables pm and h. The pressure in a McKibben muscle is always modelled by Boyle Gay-Lussac’s law, while the generated force is defined using a force balance. In the following, the Chou–Hannaford model and contributions to this model are discussed by considering the active (elongation), passive (stiffness and damping) and volumetric (change in volume) behavior of McKibben muscles. 92

McKibben muscles properties

Active behavior The actual muscle length is defined by ε in (B.1) or by h in (B.2). No definition or analogy of the (variable) nominal length is found. In (B.2), the minimal length hmin is defined as the length where the force generated by the muscle is zero, i.e. hmin has a constant value. In Section 3.4.4, it is shown that the nominal length approximates a constant value for pressures higher then 2 bar. As a result, in the high-pressure range1 the Chou– Hannaford definition of minimal length gives the same result as the nominal length. In the low-pressure range, different values of the muscle pressure yield a different muscle length with zero muscle force. Passive behavior The passive behavior of the Chou–Hannaford model is in fact given by (B.1). In (B.2), the stiffness is made more explicit by kg and kp . In [Petrovi´c, 2002a], the Chou–Hannaford model is converted into a quadratic stiffness relation, which is more obvious then the approximation made in (B.2). In [Tondu, 1997; Chou, 1994], the (static) Chou–Hannaford model is turned a dynamic model by modelling a damper parallel to the stiffness. In [Klute, 2000], the Chou–Hannaford model is extended by a Mooney-Rivlin description for the visco-elastic behavior of the rubber tube. Compared to the friction in the braid, the contribution of the rubber tube is rather doubtful. Due to the tight bonding of the braid on the tube, the deformation of the tube is complex. An experiment with an inflating rubber tube, shows that the pressure in the muscle only depends on the braid, see Appendix B.4. Volumetric behavior The variable volume of the muscle is required to model the pressure in the muscle correctly. In most cases the muscle is considered to be cylindrically shaped. In [Chou, 1996], a correction of the Chou–Hannaford model for the concave endings is given. The approximation of cylindrical shaped muscles does not seem to be a very good approximation when looking at the Shadow McKibben muscles. The pneumatic muscles used in [Chou, 1996; Tondu, 1997; Colbrunn, 2001] are all cylindrical shaped muscle types (like Rubbertuators). In Figure 2.6, the shape of the different muscle types can be compared. Other models Besides the Chou–Hannaford model, also other models are defined. In most cases, the muscle stiffness and damping are fitted to quasi-static time domain measurements. This approach yields no insight in the working principle of McKibben muscles and for this reason we omit detailed referencing here. On exception is the frequency-domain modelling given in [Thongchai, 2001]. A stepped sine excitation is injected at 0.5, 1.0, 5.0, 10.0 and 15.0 Hz and fitted in the s-domain to obtain a transfer function from pressure to joint angle. The results are used to generate a model-based fuzzy controller. 1

The low-pressure range are pressures below 1.5 bar while pressures above 1.5 bar are in the highpressure range. The pressures used in the robotic arm are in the low-pressure range.

93

In [Bertetto, 2004], a numerical muscle model is presented. A McKibben muscle is modelled in the FEM package A NSYS. The rubber is modeled as a hyperelastic elements with eight nodes and a Mooney-Rivlin behavior; each thread of the braid is modeled by a single tension bar element with linear elastic behavior. Simulation on the force - length relation are validate by experiments. No further relations are derived from the FEM model. This approach might lead to increased insight, but is not useful in a control loop. Discussion The parameters used in (B.2) lack a clear physical interpretation. Also ξ0 in the Chou–Hannaford model is not a well defined parameter, since it is not constant over the length of the muscle and it is hard to measure. This is subscribed by [Reynolds, 2003]; it criticizes the absence of a phenomenological model, but fails in presenting a better relation. The results of the Chou–Hannaford model are good in the high-pressure range. In the low-pressure range, the muscle force is wrongly predicted. The Chou– Hannaford model is commonly used, as in most cases the pressure application range is larger then 2 bar. This is a result of the commonly used Rubbertuator, which is, just like Festo’s MAS, meant for the high-pressure range.

B.2.2 Control strategies for systems with PAM’s Numerous references can be found about the control of McKibben muscles. The controllers are based on models and / or on measurements. In this section an overview of control strategies is given and discussed. In [Pack, 1994], three control strategies (PID, PID tuned by a fuzzy supervisor and nonlinear) are compared. A model of two muscles in an antagonist setup is generated, consisting of the robot dynamics, and the Chou–Hannaford muscle model. Simulations of the joint angle show that the PID controller gives a slow response, while the best response is given by the nonlinear state feedback controller. In [Schröder, 2003], an experimental setup of two antagonist muscles is considered. Instead of the length of force, the torque of both muscles, the joint angle and the pressure difference between the muscles is considered. A double loop of two PI controllers is defined, one loop controls the torque, the other the pressure. The torque is reconstructed by a model describing both muscles; it uses the Chou– Hannaford model extended a by fitted damping relation. Rather common is the use of some sort of learning control. In [Petrovi´c, 2002b] a predictive fuzzy control approach is used to control a single muscle in simulation. The predictive fuzzy model is partially based on measurements and partially on the model given in [Petrovi´c, 2002a]. In [van der Smagt, 1996], a neural network is used in a feedback controller; simulations of a robotic arm model (using the Chou–Hannaford model) as well as experiments are used to verify the control strategy. An adaptive control strategy is proposed by [Caldwell, 1994; Caldwell, 1995]; by updating an s-domain polynomial every sample, a reconstruction of the robot is made. The model parameters are used in the controller polynomials, this strategy is experimentally verified. 94

McKibben muscles properties

In [Maas, 2005] a predictive (neural network) model, based on experiments, is used in a reinforcement learning control approach2 . The implementation of this control strategy allows only track step responses; this is validated by experiments. No tracking guarantees can not be given. Three nonlinear control strategies are compared by simulation in [Carbonell, 2001]. A robust backstepping controller, a sliding mode controller and a gain scheduling controller are applied in M ATLAB simulations on a data fitted muscle model. The sliding mode and the backstepping controller show the best error tracking; however the sliding mode controller gives a chattered response. In the following two references, simulations of very long McKibben muscles are considered; long muscles approximate the cylindrical shape best. An antagonist setup is used by [Lilly, 2005] to apply a sliding mode tracking control approach. A variation on a sliding mode controller is proposed in [Repperger, 1998] to control a single muscle. It defines a variable structure controller, which uses a sliding mode control with a stabilized switching function. For the stabilization, a Lyapunov distance function is used. Experimental data is used to generate the controller. The design and control of a 4DoF anthropomorphic arm is discussed in [Tuijthof, 2000]. The arm is driven by changes in muscle stiffness. Open loop stiffness control is applied on one joint with two muscles. Working points are defined by a fixed combination of joint angle and rotation stiffness (pressure). No experimental time responses are given, only the equilibria are experimentally verified. Finally, some remarks on the use of binary valves. In the above mentioned articles servo valves are used. In the robotic arm, eight binary valves are used. In [van Ham, 2003], experiments with a robotic arm actuated by pleated artificial muscles and binary valves is discussed. A simple feedback (no controller) is defined. Two data conversions to drive the binary valves are considered: a bang – bang algorithm (using the sign of the error as valve signal) and a PWM algorithm. The PWM algorithm yields the best results.

B.3 Example of a McKibben muscle used as actuator Commonly used actuators are either position or force actuators. A force actuator applies a known force on a load, independent of the position; a position actuator imposes a displacement on a load, independent of the required force. The compliance of a McKibben muscle indicates that it is neither a position, nor a force actuator. The McKibben muscle is considered to be a non-stiff, force generating actuator. In Figure B.1, an example is given of a McKibben muscle used as actuator. Two length definitions are introduced. The actual length h is defined by the length of a muscle mounted in a robotic setup, i.e. the distance between the two mounting fittings. The nominal length h0 is defined as the length the muscle takes as a function of the muscle pressure pm , if no external force is applied. 2

The reinforcement learning control approach in [Maas, 2005] is validated on the same robotic arm as discussed in this report. Part of this thesis aims to compare the RL controller with a linear control strategy; this is discussed in Section 5.4.

95

m

h(θ)

θ

k(pm , h) h0 (pm )

Fm

d(pm ) m

pm , Vm (h)

qout

qin pin

vin

vout

pout

Figure B.1: A McKibben muscle is mounted between the fixed world and a load m. The pressure system with inlet pressure pin and outlet pressure pout is connected to the McKibben muscle; the flows qin and qout towards and from the muscle are controlled by respectively the valves signals vin and vout . The muscle pressure is given by pm and the volumetric behavior is given by Vm (h). In the upper part, a mechanical interpretation of the McKibben muscle is given. The active behavior is indicated by nominal length h0 and the passive behavior is given by stiffness k and damping d. The actual muscle length h(θ) depends on angle θ of load m. The generated muscle force Fm (k, d, h, h0 ) induces the acceleration of load m.

Imagine this system to be in a stable equilibrium with pm = 0.4 bar and θ = 0 rad. Enabling vin yields a flow qin towards the muscle; as a result the pressure pm in the muscle increases. The increase in pm yields three effects: • the flow qin decreases;

• the nominal length h0 decreases;

• the muscle’s stiffness k and damping d increase.

As the actual length is not instantly altered due to the inertia of the load m, the actual length is larger than the nominal length, i.e. the muscle is stretched by h−h0 . The stretched muscle in combination with the stiffness k and damping d, results in a force Fm . This muscle force works on the load m, resulting in a motion of m to the left; as a result, the actual muscle length decreases. The decrease in muscle length h yields three effects: • the muscle force decreases due to a decrease of h − h0 ; 96

McKibben muscles properties

• the muscle volume Vm depends on h; as h decreases, Vm increases;

• the stiffness k decreases as h decreases.

The change in stiffness k and actual length h, affect the muscle force Fm ; the change in muscle volume Vm affects the pressure pm . These changes, in turn, affect the nominal length h0 , the flow qin , etc. . . With this exposition, we aim to clarify the complex behavior of the McKibben muscle.

B.4 The tube of a McKibben muscle The rubber tube of a Shadow McKibben muscle is in fact a piece of a bicycle inner tire. The rubber type is isobutylene, which has an excellent airtightness property. The stress at fracture is up to 700% [Stiomak, 1997]. Relations describing the behavior of spherical and cylindrical balloons are derived in [Müller, 2004]. In Figure B.2a, the typical behavior of inflating a spherical balloon is shown. A lot of pressure is required, before the starts to inflate. From this moment on, less pressure is required for further expansion. Similar behavior is found for cylindrical balloons, as shown in Figure B.2b. The major difference is that the pressure has a true maximum here. Besides, a relatively high pressure is required before the diameter of the balloon starts to increase. In Figure B.2c, the variation in length is plotted against the variation in radius. Inflating the balloon will first cause an expansion in both diameter and length; in the end, the length will reach a maximum and only the diameter will increase. An experiment is performed with an adapted version of the McKibben muscle. The braid is stripped from the muscle; only the rubber tube and the fittings remain. In Figure B.3, two photos of the experiment are shown. In both cases, the balloon is at approximately 0.3 bar, the (big) difference is the volume of both muscles. Most interesting is that the pressure in a cylindrical balloon is bounded according to Figure B.2b. More air into the balloon will only lead to an increasing volume, and not to an increase in pressure. Similar behavior is found during the experiment. The expansion of the rubber tube starts at approximately 0.4 bar; afterwards, the pressure remains inbetween 0.3 and 0.4 bar, independent of the amount of air. As air is an ideal gas; the pressure, volume and air mass in the rubber tube are related to each other by Boyle-Gay Lussac’s law, see also (3.16). The experiment shows that the pressure remains constant and that the increased air mass yields an increase in volume. A McKibben muscle shows different behavior; an increase in air mass yields a minor change in volume and an increase in pressure. For this reason, it is believed that the ribber tube does not contribute to the pressure relation of the McKibben muscle and that the volume of the muscle is completely determined by the nylon braid. The rubber tube, however, is still likely to contribute to the passive and active behavior of the McKibben muscle. This is also shown in [Klute, 2000]; the Chou– Hannaford model is extended by visco-elastic relations, to take the behavior of the rubber tube into account.

97

(a) p(µ) spherical

(b) p(µ) cylindrical

(c) λ(µ) cylindrical

Figure B.2: A balloon’s length stretch ratio is defined by µ = ll0 , with l the balloon’s length and l0 the balloon’s resting length. A balloon’s radius stretch ratio is defined by λ = rr0 , with r the balloon’s radius and r0 the balloon’s resting radius. The pressure as function of strain for a spherical balloon, Figure (a) and a cylindrical balloon, Figure (b). The length - radius ratio for a cylindrical balloon is shown in Figure (c). The values along the axis are of minor interest, the shown behavior is important. Images taken from [Müller, 2004].

6 6

170 mm

360 mm

?

?

(a)

(b)

Figure B.3: The rubber tube of a McKibben muscle, before (a) and after (b) inflating, the pressure in both cases is approximately 0.3 bar. The length of the tube is indicated in the figures.

98

C Pulse width modulation Pulse width modulation is applied to convert a continuous input signal into a discrete output signal with equal power. The output signal consists of pulses with value 1 and a variable pulse width. In the robotic arm setup, a distinction is made between a positive or a negative value of the input signal. A positive input signal refers to the flow-in valve while a negative input value refers to the flow-out valve. The PWM algorithm uses a sample frequency fPWM . On every sample the value of the input signal is taken as a reference for the pulse width. The PWM frequency is based on a comparative assessment of four properties: • A high fPWM enlarges the bandwidth of the PWM algorithm, a rule of thumb 1 fPWM . is BW = 10 • The number of valve actions every second is equal to fPWM , so a high frequency might make the robotic arm act nervously. • The valves have a limited switching frequency (in this case fv = 300 Hz). If fPWM approaches fv , there is not much variation possible in pulse width. • The use of hoses in the pneumatic system introduces low-pass behavior. A sufficiently high fPWM smoothens the pulsating air flow. A test setup, as given in Figure C.1, is used to determine the influence of fPWM on the time behavior of the pressure. In Figures C.2 to C.4, three experiments with different fPWM are shown. As can be seen the pressure in the volume pout smoothens rather well if fPWM = 50 Hz; at fPWM = 100 Hz, pout hardly shows any ripple at all. Based on the value of fPWM and pout , it is decided to set fPWM = 50 Hz. A higher fPWM is considered to be too close to fv . Notice that fPWM does not affect the pressure rise time. Due to the choice for fPWM , the bandwidth of the PWM algorithm is 5 Hz. This is in the same range as the expect bandwidth of the robotic arm.

pin

pout pvalve

Figure C.1: Three pressures are measured to determine the influence of fPWM .

1

p [bar]

0.8

0.6

0.4

pin pvalve pout

0.2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

t [s] Figure C.2: Inflating of a constant volume through a binary valve. The valve is controlled by a block signal of 25 Hz and the pulse width is 50%.

1

p [bar]

0.8

0.6

0.4

pin pvalve pout

0.2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

t [s] Figure C.3: Inflating of a constant volume through a valve. The valve is controlled by a block signal of 50 Hz and the pulse width is 50%.

1

p [bar]

0.8

0.6

0.4

pin pvalve pout

0.2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

t [s] Figure C.4: Inflating of a constant volume through a valve. The valve is controlled by a block signal of 100 Hz and the pulse width is 50%.

100

2

Pulse width modulation

C.1 Formulation of the PWM algorithm A PWM algorithm samples the input signal by a frequency fPWM . This means that the pulse width is calculated for the value of the input signal at the sample time. It is decided to use a zero-order hold (ZOH) as this introduces the least phase delay. For the ZOH, a discrete PWM time tk (t) is introduced, the samples are given by k(t): k(t)

= ⌊t fPWM ⌋

tk (t) =

(C.1)

k(t) . fPWM

(C.2)

The first step in the algorithm is to apply the saturation, sign splitting and sampling on the input signal u(t), this is shown in (C.3). The sampling is performed by taking u(tk ) of u(t). The sign splitting is taken into account by applying two saturations with different boundaries on u(tk ). uk,in (t) = uk,out (t) =

sat{0,1} (u(tk )) sat{−1,0} (u(tk )) .

(C.3)

The parameters describing a pulse are shown in Figure C.5, the pulse width is given 1 . by s, where s ∈ [0, TPWM ], with TPWM = fPWM 1

vin

0 s TPWM Figure C.5: Pulse parameters.

The basic property of the PWM algorithm is to keep the power of the input and the output signal equal. This property is used for the inputs uk,in (t) and uk,out (t), but not for the original input u(t). The algorithm is explained for uk,in (t). Solving the integral of uk,in (t) and vin (t) over one sample time, determines the value of s: T ZPWM

T ZPWM

vin (t) dt

(C.4)

uk,in (t) TPWM

= 1 · s + 0 · (TPWM − s)

(C.5)

uk,in (t) TPWM

= s.

(C.6)

uk,in dt

0

=

0

101

To ensure the valves to open at t = tk and to close t = tk + s, a threshold is defined. The threshold uth ∈ [0, 1] is a sawtooth function which runs from 0 to 1 every TPWM , as given by: uth (t) = t fPWM − k(t) .

(C.7)

The valves have a limited switching time. It is of no use to demand the valves to open for a time less than the switching time; namely, this will result in errors and possible limit cycling of the control loop. To prevent limit cycling, a saturation τ is introduced on uth : vth (t) = sat{τ,1−τ } (uth (t)) ,

(C.8)

with vth ∈ [τ, 1 − τ ]. A rule of thumb for τ is τ=

2 , fv

(C.9)

where fv is the maximum valve switching frequency. The PWM output signal is given by comparing the threshold vth (t) to the sampled input uk,in (t). Finally, the valve signals are derived by comparing the flow signal uk,in (t) with the threshold vth (t). This is shown in (C.10) for both vin (t) and vout (t); remark the additional minus in the vout (t) relation. vin (t)

= uk,in (t)

> vth (t)

vout (t)

= uk,out (t)

< −vth (t)

(C.10)

In Figures 3.3 and C.6, the PWM results are shown for two different fPWM . In case of the robotic arm, four input signals u(t) are transformed into eight discrete output signals v(t), defined by: £ ¤T uA (t) uB (t) uC (t) uD (t) u(t) = (C.11) · ¸T vin,A (t) vin,B (t) vin,C (t) vin,D (t) v(t) = (C.12) vout,A (t) vout,B (t) vout,C (t) vout,D (t) Remarks The main advantage of this algorithm is the use of the threshold vth (t) as it allows for real time implementation in dSPACE. Errors due to the zero-order hold are not compensated for. As the PWM algorithm is part of the control loop, errors are corrected directly. Besides, a static error is introduced by the dead-zone due to the saturation τ .

102

Pulse width modulation

τ

vth (t)

1

ui (t) ui (tk ) 0.5

τ

0 0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

vin,i (t)

0

t [s] Figure C.6: Example of the PWM algorithm with fPWM = 50 Hz and τ = 0.01 s.

C.2 Saturation function It is assumed that a standard saturation function is defined with the boundaries −1 and 1 as given by: y

=

SAT (x)

y y y

= −1 = x = 1

∀ ∀ ∀

x < −1 −1≤x≤1 x>1.

(C.13)

In this thesis, an alternative notation is introduced which enables for a saturation with variable boundaries: y

=

sat{a,b} (x)

y y y

= a = x = b

∀ ∀ ∀

xb,

(C.14)

with a the lower boundary and b the upper boundary. The implementation of (C.14) is given as a function of (C.13) according to: ¶ µ x−α +α, (C.15) sat{a,b} (x) = β SAT β with a+b (C.16) 2 b−a . (C.17) β = 2 In Figure C.7, four examples of the saturation function as defined in (C.14) are given. In the upper right figure is shown that sat{−1,1} (x) yields the same result as SAT (x). α

=

103

y

y = SAT (x) = sat{−1,1} (x)

y = sat{−1,0} (x)

2

2

1

1

0

0

−1

−1

−2 −2

−1

0

1

−2 −2

2

y

y = sat{0.1,0.9} (x) 2

1

1

0

0

−1

−1

−1

0

x

0

1

2

y = sat{−1.5,0.3} (x)

2

−2 −2

−1

1

−2 −2

2

−1

Figure C.7: Example of the used saturation function.

104

0

x

1

2

D Considerations on the pneumatic model

In this appendix, the building and modelling of a pneumatic test setup is discussed. Based on these results, the model of the pneumatics in the robotic arm are defined. Two types of behavior are investigated with the pneumatic test setup; the type of flow through the hoses and how to implement the valves in the model. In Appendix D.1, the setup of the pneumatic test setup discussed. In Appendix D.2, two flow descriptions are defined and, in Appendix D.3, the modelling of the valves is discussed. Finally, in Appendix D.4, simulations will show which flow and which valve model is preferred.

D.1 Pneumatic test setup To model the pneumatic system, a simple test setup has been built which resembles the pneumatics for one muscle in the robotic arm, see Figure D.1 and Figure D.2. The test setup consists of a capacity connected to the pressure from the wall. A smaller capacity (with a fixed volume) is used instead of a muscle. A muscle has a variable volume which is not desirable in the modelling of the pneumatic system as the change of volume is not exactly known. A similar solenoid valve as in the robotic arm is used mounted inbetween the two capacities. Pressure sensors are mounted on three places: on the two capacities and behind the valve. A current sensor for the valve is made of several resistances in order to be able to determine the exact moments of activating the valves. The manual switch is mounted to release the air from the small capacity. A signal generator and an amplifier are used to direct the solenoid valve. As the valves in the robotic arm are directed by binary signals, similar signals are used in the test setup. Two types of signals are used: • The response of the system for instant opening of the valve is measured, • The response of the system using a pulsating valve signal. The frequency and the duty-cycle of the pulse can be varied. Further, a dSPACE system is used for data acquisition. The measurements are imported into M ATLAB for comparison with the simulations.

signal generator pressure in

amplifier capacity current sensor sensor 3

sensor 1 hose 1

muscle

manual switch sensor 2

valve

hose 2

Figure D.1: Pneumatic test setup.

p1 pin

qin

p2

q1

q2

pc air supply

pm

capacity

valve

hose 1

hose 2

muscle

Figure D.2: Schematic representation of the pneumatic test setup; the flows are indicated by q and the pressures by p.

p1

p2 q Figure D.3: Flow model.

106

Considerations on the pneumatic model

D.2 Flow modelling Every pneumatic system must fulfill mass conservation. For this reason flowing air is described by mass flow, while in general flow velocity is measured. Velocity v and mass flow q are related to each other by the density ρ and the flow area A: Z q = ρ v dA . (D.1) A

The flow through a hose depends on the pressure difference between the two sides of the hose, as shown in Figure D.3. In this section, two relations for q(p1 , p2 ) are considered. Flow behavior dominated by the hose geometry and flow behavior dominated by the geometry of throttle, which is present in the hose. Laminair flow through a hose In this case, it is assumed that the flow through the hoses is laminar and that the hose itself determines the flow, see Figure D.4. A laminar flow through a round hose is described by the Hagen-Poiseuille law, [Polytech, 1997]. This is the solution for (D.1) in case of the quadratic velocity profile v: ¡ ¢ q = γ ρ p1 − p 2 , (D.2) with the flow resistance γ=

πrh4 , 8ηlh

(D.3)

the radius of the hose rh , the length of the hose lh and the dynamic viscosity of air η. The density of the air is a function of the pressure. As the pressure changes with the length, the mean of the two pressures p1 and p2 as defined in Figure D.3, is used to derive the density: ρ=

p1 + p2 + 2 patm . 2Rs T

(D.4)

Flow through a throttle In the hoses, valves are present. These valves act as throttles in case the valve radius is small compared to the hose radius, see Figure D.5. In this case, the flow

lh 2rt q

v

2rh

Figure D.4: Laminar Poiseuille flow.

q

Figure D.5: Throttle valve flow.

107

behavior is determined by the valves. The relation describing the flow through a throttle valve [Polytech, 1997] is given by: q=β with

p

ρ (p1 − p2 ) ,

β = 0.59

√

(D.5)

2 πrt2

(D.6)

and rt the throttle radius. The density is defined by (D.4).

D.3 Valve modelling The valves used in the robotic arm and in the pneumatic test setup are fast-switching solenoid valves. In Section 3.2, it is shown that the dynamic behavior of the valves is not important for modelling the robotic behavior. For this reason, a change in valve state (open or closed) is instantaneous. In that section, two methods for modelling the valves (taking its status in account) are considered. In order to explain the two valve implementations, the general setup of the pneumatic test setup model is discussed first. A schematic overview of the pneumatic test setup is already shown in Figure D.2. The pressure in both volumes satisfy (3.20) and the flow in the hoses is given by either (D.2) or (D.5). Second-order valve description One way to model the solenoid valve in Figure D.2, is to consider it as two small volumes which are interconnected by a short hose. Depending on the valve status, the short hose is open or closed. This is visualized in Figure D.6. The small volumes are defined by the volume in solenoid valve and of the hoses. A consequence of this approach is that each volume adds one order to the model as it is described by the differential equation given in (3.20). When considering Figure 2.2, it can be noticed that in case of the robotic arm, five orders would be added to the model; one volume for the splitter and four volumes for the backside of the valves. Zeroth-order valve description Another way to model the solenoid valve is given in Figure D.7. The valves are not modeled at all, only the function of the valves is implemented into the hose between the capacity and the muscle volume. The status v of the solenoid valve determines whether flow is possible or not. Using this approach, no additional orders are added to the model as no volumes are introduced. This might be an advantage in solving the model. In case of the robotic arm, multiple McKibben muscles and valves are connected to the capacity, by a single hose which is split into four hoses, see Figure 2.2. Similar as for the valves, the splitter can not be modeled as a small volume, a different relation is derived for the implementation of the splitter behavior. This is presented in Section 3.3.2. 108

Considerations on the pneumatic model

p1 , V 1

p2 , V 2

capacity pin

qin

muscle

pc , V c

q1

qv , v

q2

pm , V m

Figure D.6: Model of the pneumatic test setup with a second-order valve description.

capacity pin

qin

pc , V c

muscle v

qs

pm , V m

Figure D.7: Model of the pneumatic test setup with a zeroth-order valve description.

D.4 Simulation of the pneumatic test setup In this section, the two flow relations of Appendix D.2 and the two valve relations of Appendix D.3 are implemented into a model of the pneumatic test setup. The first model is called PTS1. It combines the second-order valve model with the Poiseuille flow. Applying (D.2) for the hoses and (3.20) for the four volumes, gives: ¡ ¢ p˙c = Rs T Vc−1 qin − q1 ¡ ¢ p˙1 = Rs T V1−1 q1 − qv ¡ ¢ (D.7) p˙2 = Rs T V2−1 qv − q2 p˙m

= Rs T Vm−1 q2 ,

with the Poiseuille flows given by: ¡ ¢ qin = γin ρin pin − pc ¡ ¢ q1 = γ1 ρ1 pc − p1 ¡ ¢ qv = v γv ρv p1 − p2 ¡ ¢ q2 = γ2 ρ2 p2 − pm .

(D.8)

In the second model (called PTS2), the second-order valve model is combined with the throttle flows of (D.5). The pressures are also given by (D.7), but the flows are given by: p qin = βin ρin (pin − pc ) p q1 = β1 ρ1 (pc − p1 ) p (D.9) qv = v βv ρv (p1 − p2 ) p q2 = β2 ρ2 (p2 − pm ) .

All flows and pressures are indicated in Figure D.6. The flow resistance γ, the density ρ and the throttle parameter β are defined by respectively (D.3), (D.4) and (D.6). 109

1

0.9

0.8

0.7

p [bar]

0.6

0.5

0.4

0.3

simulation pc simulation p1 simulation p2 simulation pm experiment pc /p1 experiment p2 experiment pm

0.2

0.1

0

0

0.5

1

1.5

2

2.5

3

3.5

4

t [s] Figure D.8: Simulation of PTS1 compared to the measurement for instant activation of v. An ideal fit would be if the red (pm ) and green (p2 ) dashed lines would match the red and green solid lines. The solid blue lines must be inbetween the dashed blue (pc ) and black (p1 ) line. 1 0.9 0.8 0.7

p [bar]

0.6 0.5 0.4 0.3

simulation pc simulation p1 simulation p2 simulation pm experiment pc /p1 experiment p2 experiment pm

0.2 0.1 0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

t [s] Figure D.9: Simulation of PTS2 compared to the measurement for instant activation of v.

110

1.6

Considerations on the pneumatic model

The valve status v is a binary signal; the gas temperature T and the specific gas constant Rs are constants. In the pneumatic test setup, v is generated by a signal generator. Both models are compared with an experiment. In PTS1, the γ parameters are defined using the hose dimensions. Unfortunately, this yields PTS1 to be unstable. To check whether PTS1 will ever be able to yield a good fit, the γ parameters are tuned to the best possible fit between the simulation and the experiment, This is shown in Figure D.8. The dashed lines represent the model, while the solid lines represent the measurements. As can be seen, the Poiseuille flow it is not likely to describe the time behavior of the pressure accurately. It can be concluded that the flow is not dominated by Poiseuille behavior and that the discrepancy between the simulation experiment is not due to errors in the parameters. Using the the diameters of the reduction valve and solenoid valve, the parameters β in PTS2 are defined. in Figure D.9 the simulation is compared to the experiment. In this case, the model shows a good resemblance with the experiment. Also simulations with combinations of (D.8) and (D.9) are performed, although this yielded no improvement of the simulation as given in Figure D.9. It can be concluded that the flow in the pneumatic test setup and in the robotic arm, are determined by the geometry of the reduction valve and the solenoid valves. In Figure D.10, a second example of PTS2 is shown. The response of PTS2 and a measurement are compared to each other when v is a pulsating signal at 25 Hz with 50% duty cycle is given. The match between experiment and simulation is good. In Figure D.9, it can be seen that PTS2 becomes slightly unstable after simulating 1.2 seconds. From this point on, the pressures are almost equal to each other. As a result, numerical errors occur in the simulation. The reason for this is twofold: cancelation (lost of significant numbers due to subtraction) combined with stiff differential equations [Heath, 2002]. Using the second-order valve model, two small volumes are introduced. These small volumes have a small time constant compared to the relatively large volumes of the capacity and the muscle yielding a set of stiff differential equations. In Figure D.11, the flow belonging to the simulation in Figure D.9 is given. As can be seen, large fluctuations occur due to cancelation in q1 and q2 . To prevent numerical errors, it is decided to model the valves using the zeroth-order model. The third model (called PTS3) combines the zeroth-order valve model with the throttle flow relation as given by: ¡ ¢ p˙c = Rs T Vc−1 qin − qs p˙m

qin qs

= Rs T Vm−1 qs

(D.10)

p ρin (pin − pc ) p = v βs ρs (pc − pm ) . = βin

The flows and pressures are indicated in Figure D.7. This model has proven to work in simulations of the robotic arm model. In Figure D.12, a simulation of PTS3 is given, which proves that this model is accurate and stable. 111

1 0.9 0.8 0.7

p [bar]

0.6 0.5 0.4 0.3

simulation pc simulation p1 simulation p2 simulation pm experiment pc /p1 experiment p2 experiment pm

0.2 0.1 0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

t [s] Figure D.10: Simulation of PTS2 compared to the measurement. The valve is excitated by a block signal at 25 Hz with a 50% duty cycle.

−4

10

x 10

qin q1 q2

8

q [kg s−1 ]

6

4

2

0

−2

−4

0

0.2

0.4

0.6

0.8

1

1.2

1.4

t [s] Figure D.11: The flow belonging to the simulations of PTS2 as given in Figure D.9. No experimental flow is displayed as the flow is not measured. The numerical errors occur as the large fluctuation in flow starting at approximate 1.25 seconds.

112

1.6

Considerations on the pneumatic model

1 0.9 0.8 0.7

p [bar]

0.6 0.5 0.4 0.3 0.2

simulation pc simulation pm experiment pc /p1 experiment p2 experiment pm

0.1 0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

t [s] Figure D.12: Simulation of PTS3 compared to the measurement for instant activation of v, this experiment is similar to the experiment given in Figure D.9. The simulation of PTS3 returns only two pressures, see Figure D.7, which are stable from 1.25 s on.

D.5 Two DoF pneumatic model Given the DoF pneumatic model in (3.33), the 2DoF pneumatic model can be derived. In Figure D.13, a schematic representation of the 2DoF pneumatic model is shown. The major difference is the implementation of the splitter. The relations of the other parts are derived similarly as in Section 3.3. The virtual pressure in a splitter with two outputs is described by (3.29). A splitter with four outputs can be considered as two splitters with two outputs, resulting in two virtual pressures: ps,A B and ps,CD . The virtual pressure ps is now defined by the lowest of the two virtual pressures: ps = min(ps,A B , ps,CD ) .

(D.11)

The relation describing the flow out of the splitter is slightly changed. For the case of flow qs,A , the proportionality factor must take the other three possible flows into account, according to:

qs,A =

vin,A (pc − pA ) qs . (D.12) (pc − pA ) + vin,B (pc − pB ) + vin,C (pc − pC ) + vin,D (pc − pD )

Similar relations hold for qs,B , qs,C and qs,D . 113

vout,A qout,A

qs,A pA , V A

vout,B qout,B

qs,B

vin,A vin,B

pB , V B pin

qin

qs pc , V c

vout,C

ps

qout,C

qs,C pC , V C vin,C vin,D vout,D

qout,D

qs,D pD , V D Figure D.13: Schematic representation of the 2DoF pneumatic model.

114

Considerations on the pneumatic model

The complete pneumatic model for the 2DoF robotic arm is now given by: ¡ ¢ p˙c = Rs T Vc−1 qin − qs (D.13) ps,A B

=

ps,CD

=

ps

=

£

£

1 − vin,A

1 − vin,C

qs,C

=

qs,D

=

ρ

=

¸·

1 − vin,B vin,B

¸

(D.14)

·

pc pC

pD min(pC , pD )

¸·

1 − vin,D vin,D

¸

(D.15)

p ρ (pc − ps )

= vout,i βout

=

pB min(pA , pB )

(D.17)

= βs

qs,B

pc pA

p ρ (pin − pc )

qs

=

¤

·

(D.16)

= βin

qs,A

vin,C

¤

min(ps,A B , ps,CD )

qin

qout,i

vin,A

p

(D.18)

ρ (pi − patm )

for

i=

£

A

B

C

D

¤

(D.19)

vin,A (pc − pA )qs (pc − pA ) + vin,B (pc − pB ) + vin,C (pc − pC ) + vin,D (pc − pD ) (D.20) vin,B (pc − pB )qs vin,A (pc − pA ) + (pc − pB ) + vin,C (pc − pC ) + vin,D (pc − pD ) (D.21) vin,C (pc − pC )qs vin,A (pc − pA ) + vin,B (pc − pB ) + (pc − pC ) + vin,D (pc − pD ) (D.22) vin,D (pc − pD )qs vin,A (pc − pA ) + vin,B (pc − pB ) + vin,C (pc − pC ) + (pc − pD ) (D.23) p1 + p2 + 2patm . 2Rs T

(D.24)

115

116

E Considerations on the muscle model

In the literature, three models are stated to describe the stiffness of a McKibben muscle. These models are given in [Chou, 1996; Tondu, 1997; Carbonell, 2001]. In these references, the same approach in obtaining the models is adopted; based on the principle of virtual work, a relation between the length, pressure and force is determined. The difference between the three models is that a different set of parameters are used in each reference. However, in [Schröder, 2003], it is stated that the three models can be transformed into each other. In Section E.1, the Chou– Hannaford model as defined in [Chou, 1996] is explained. A minor adjustment to the model Chou–Hannaford is given in [Petrovi´c, 2002a]. The Petrovi´c model is treated in Section E.2. In Section E.3, a model is stated, that is similar to the Petrovi´c model, but is adapted to the Shadow muscles. This stiffness model is used in the model of the McKibben muscles, as presented in Section 3.4.

E.1 The model of Chou – Hannaford The variables introduced in this section, differ from the variables in as defined in [Chou, 1996], but are chosen according to the nomenclature on page xi. The muscle is considered as a cylinder with variable length and diameter. The length of the cylinder is given by h and the diameter by D. Since the diameter of the cylinder does neither correspond to d1 nor to d2 , a new diameter D is introduced. The length 2L of each shell thread is considered to be constant. A variable γ is defined which is a measure for the angle characterizing the orientation of the thread, see Figure E.1. This angle is considered to be constant (a result of the assumption that the muscle is cylindrical), while in reality it depends on the position along the muscle. The number of windings of one thread is given by n. From these parameters, it can be seen that the thickness of the braid and the tube

γ

1 thread

h

unfold

2L

nπD Figure E.1: Muscle parameters in the Chou–Hannaford model.

are not taken into account. From Figure E.1, the following relations can be derived: h =

2L cos γ

(E.1a)

D

=

2L sin γ nπ

(E.1b)

V

=

π 2 2L3 D h = sin2 γ cos γ 4 πn2

(E.1c)

In [Carbonell, 2001], it is noticed that the braid angle γ is difficult to measure. For this reason, a different set of parameters is introduced. The static physical relations are derived by considering the pneumatic to mechanical energy conversion. This conversion is realized by transferring the pressure on the inner surface of the bladder into a shortening tension. The virtual input work δWin is done when gas pushes the bladder surface: Z δWin = p δ I~ · ~si = p δV , (E.2) Si

with p the relative pressure in the muscle, δ I~ the inner surface virtual displacement, ~si the area vector, Si the total inner surface and δV the virtual volume change. The virtual output work δWout is done when the actuator shortens with the volume change: δWout

= −F δh ,

(E.3)

with F the force performed by the muscle and δh the virtual length change. Application of the principle of virtual work for a system in a static equilibrium gives: δWout

= δWin

(E.4)

−F δh

= p δV

(E.5) 118

Considerations on the muscle model

This involves two assumptions: no energy is stored in the muscle and there are no other system losses. Using (E.1), we can write: ¢ 2L3 ¡ 2 sin γ cos2 γ − sin3 γ δγ 2 πn δh = −2L sin γ δγ .

δV

=

(E.6)

Substituting of (E.6) in (E.5) and realizing that it should hold for any δγ gives the Chou–Hannaford model as: F

¢ L2 ¡ δV = p 2 cos2 γ − sin2 γ 2 δh πn ¢ L2 ¡ 3 cos2 γ − 1 . = p 2 πn

= −p

(E.7)

This completes the static muscle relation of the Chou–Hannaford model, in [Chou, 1994] almost the same theory is discussed. In [Tondu, 1997] a similar relation for the muscle force F is derived, but as a function of the muscle contraction ratio ε. ¶ µ 3 L2 1 2 . (E.8) F = p (1 − ε) − πn2 tan2 γ0 sin2 γ0

The newly introduced parameter γ0 = 20◦ is the reference braid angle. In this reference, it is also noted that (E.8) is similar to the Chou–Hannaford model. In order to simplify the model, the tension is considered as a function of the pressure and the length. This means that the muscle is like a variable stiffness element. Its stiffness k is proportional to the pressure and the stiffness per unit pressure kg approximates a constant. The stiffness k and the stiffness per unit pressure kg are defined according to: k

=

kg

=

dF dh dk . dp

(E.9)

As a result, the force can be linearized as: F

= kg p (h − hmin ) ,

(E.10)

with hmin the theoretical minimum muscle length (when F = 0). In order to taken into account the stored energy in the bladder and shell, a threshold pressure pth is defined for the pressure to overcome the radial elasticity of the bladder for expanding. This gives F

= kg (p − pth ) (h − hmin )

if p > pth ,

F

= kp (h − hmin )

if p ≤ pth .

(E.11)

This model is compared with measurements as described in [Chou, 1996]. The tests are quasi-static tests: a slow elongation of the muscle under a given pressure. The model fits the measurements rather well for pressures above 2.0 bar. Also dynamic tests are performed and some related conclusions are drawn in the reference: 119

• The tension - length behavior is velocity independent at low velocities; • The hysteresis is history dependant; • The velocity independent hysteresis is most likely due to Coulomb friction, which dominates the total friction of the actuator.

E.2 The model of Petrovi´c The model in [Petrovi´c, 2002a] uses the Chou–Hannaford model. Petrovi´c avoids the simplification as given by (E.10); it is considered for the following reason. The force F is a linear function of the pressure p and quadratic with the length h, which h as given in (E.1a) into (E.7) and reorganizing can be seen by substituting cos γ = 2L terms: F

= p

S(q)

3h2 − 4L2 = p S(q) 4πn2

= s0 + s1 q + s2 q 2 ,

(E.12a) (E.12b)

with q s0

= h − hr ¡ 2 ¢ 3hr − 4L2 = 4πn2

(E.13a) (E.13b)

s1

=

3hr 2πn2

(E.13c)

s2

=

3 , 4πn2

(E.13d)

and hr a constant reference length. This simple quadratic relation for S(q) is an alternative for the gas spring constant kg (h − hmin ) as derived in (E.10).

E.3 The model for a Shadow McKibben muscle The derivation of the Chou–Hannaford model using virtual work, seems very elegant from the authors point of view. For this reason, the simplification as proposed in (E.10) seems somewhat odd, especially when considering the relation as derived in (E.12). On the other hand, the model in (E.10), appears to be a good predictor for the muscle behavior for pressures of 2.0 bar and higher. In this section, a relation is derived using the muscle properties for the volume V as derived in Section 3.4 and the the Petrovi´c model, based on virtual work. The volume of a Shadow McKibben muscle is defined in (3.39) as: V (h)

=

¢ π h¡ 2 3 d1 + 4 d1 d2 (h) + 8 d22 (h) . 60 120

(E.14)

Considerations on the muscle model

In (E.5), the result of the virtual work is given by: F

= −p

δV δh

(E.15)

Assuming the muscle to be in an stable equilibrium; δV is rewritten to: δV

dV δh . dh

=

(E.16)

Substituting (E.16) into (E.15) yields: F

= −p

dV . dh

(E.17)

Now, the muscle’s stiffness k can be defined as: k

=

d2 V dF = −p . dh dh2

(E.18)

In this relation, the pressure is considered to be an input; as such, it is not a func2 tion of the length h. The derivative ddhV2 represents the stiffening of the muscle due to elongation of the muscle. The stiffness is linear with the relative pressure p, however one aspect is missing. The stiffness of the muscle is biased. For p = 0 bar, relation (E.18) returns k = 0 Nm−1 . This is not correct. Below 0.3 bar, the muscle behaves differently as the braid and the rubber tube do not touch, see also Figure 2.5. For this reason, the stiffness of the muscle is biased. This bias is implemented by introducing an additional parameter λ on the pressure. The change in stiffness as a function of the pressure also appears if p ≈ 0; for this reason, the pressure bias λ is applied directly on the actual muscle pressure p according to: k

= − (λ + p)

d2 V . dh2

(E.19)

2

As the stiffening, given by ddhV2 , is rather involved, this function is approximated by a truncated Taylor series. Simulations with several orders for the Taylor series have shown that a second-order Taylor series yields satisfying results. The simulations are compared to experimental stiffness data. These experiments are discussed in Section 4.1. The application of the Taylor series yields: −

d2 V dh2

≈ ξ0 + ξ1 h + ξ2 h2 ,

(E.20)

Combining (E.19) with (E.20), yields the final relation for the muscle stiffness: k

=

¡

λ+p

¢¡

¢ ξ0 + ξ1 h + ξ2 h2 . 121

(E.21)

E.4 Analytical solution for the stiffness parameters In the stiffness relation in (E.21), four parameters are present. The three ξ param2 eters are derived from the Taylor series of ddhV2 . The exact expression is: µ ¶2 ∂ V ∂ 2 d2 ∂2V ∂ 2 V ∂ d2 d2 V ∂ 2 V ∂ d2 + = + 2 . (E.22) + 2 dh2 ∂d2 ∂h2 ∂h ∂d2 ∂h ∂d2 2 ∂h |∂h {z } =0

This relation is estimated by a second-order Taylor series around the reference length hr = 0.19 m. This value is given by the muscle length h in case the angle of joint 1 of the robotic arm is equal to θ = 0 deg. The Taylor series for (E.22) is: ¯ ¯ ¯ d3 V ¯¯ 1 d4 V ¯¯ d2 V d2 V ¯¯ + (h − hr ) + (h − hr )2 . ≈ dh2 dh2 ¯h=hr dh3 ¯h=hr 2 dh4 ¯h=hr (E.23)

From this relation, the three Taylor constants s are given by: ¯ d2 V ¯¯ ψ0 = = 0.0373 m dh2 ¯h=hr ¯ d3 V ¯¯ ψ1 = = 0.2825 dh3 ¯h=hr ¯ 1 d4 V ¯¯ = 1.9691 m−1 . ψ2 = 2 dh4 ¯h=hr

(E.24)

The numerical values of ψ0 , ψ1 and ψ2 are determined in M ATLAB. Making use of 2 these three parameters, the relation for ddhV2 would be defined as:

d2 V ≈ ψ0 + ψ1 (h − hr ) + ψ2 (h − hr )2 , dh2 while in (E.20), it is stated that:

(E.25)

d2 V ≈ ξ0 + ξ1 h + ξ2 h2 . dh2 This discrepancy is solved by rewriting the ψ parameters to ξ parameters: ξ0

= ψ0 − ψ1 hr + ψ2 h2r

=

ξ1

= ψ1 − 2ψ2 hr

= −0.4658

= ψ2

=

ξ2

0.0547

1.9691

(E.26)

m (E.27) m

−1

.

The parameter λ is approximated using a simple experiment. A mass is attached to the rubber tube and the static elongation is measured, this yields a stiffness of 1.5 kN m−1 . Substituting (E.27) in (E.21) for p = 0 bar and h = hmin = 160 mm, yields: ¡ ¢¡ ¢ k = λ + 0 ξ0 + 0.16 ξ1 + 0.162 ξ2 = 1500 N m−1 λ

=

0.5 bar.

(E.28)

This numerical value is a coarse estimate, but it is a good order estimate of the order λ. An interpretation of the four muscle stiffness parameter is difficult, as no physical meaning can be given to the parameters and its values. 122

F Derivation of the multibody model The equations of motion of the robotic arm, will be derived by making use of the multibody theory. The notation as defined in [van de Wouw, 2003] is adopted in the following derivation. In Figure F.1, the setup of the multibody model is given; in Figure F.2 and in the upper half of Table F.1 the accompanying dimensions of the robotic arm are given. In the lower half of Table F.1, three additional lengths are defined which will be declared in the remainder of this appendix. The length lS refers tot the unloaded spring length. The length of the muscle fittings for respectively muscles A and B is given by lfA B and the length of the fittings for muscles C and D by lfCD . The kinematics are defined in relative coordinates. This means that the position of each mass is defined relative to the mass it is connected to. Three vector frames are introduced. Point O indicates the origin of the fixed frame ~e 0 . In the center of mass of body 1 (CM1), the frame ~e 1 is introduced an in the center of mass of body 2 (CM2), the frame ~e 2 is located. Furthermore, three joints are introduced. C00 is rigidly fixed to the origin O. C01 and C12 are revolute joints connecting body 1 to the fixed world respectively body 2 to body 1. Next, the vectors indicated in Figure F.1 will be explained using the dimensions as given in Figure F.2. dim

length [mm]

dim

length [mm]

dim

length [mm]

l0 l01x l01y l11x l11y l1x l1y

375.0 66.7 13.7 453.3 26.3 520.0 40.0

l1A l2A l1B l2B lA B l12 lsvx

76.0 142.0 24.0 210.0 123.3 176.0 110.0

l0C l1C l0D l1D l0S l1S lsvy

20.0 35.0 20.0 105.0 170.0 185.0 26.0

75.0

lf A B

140.0

lfCD

180.0

lS

Table F.1: Values of all dimensions as used in the multibody model.

~e22 CM2 −θ2

~e12

~b12

A C11 = C12 ~e21 ~e11 ~b11

CM1 ~b01

C01

B

θ1 ~rCM 2 ~rCM 1

~c0

S

C

D

~e20 ~e10 O = C00 Figure F.1: Kinematics of the multibody system.

124

Derivation of the multibody model

sv lsvx lsvy

~bsv C01

CM1 ~rsv

~c0

CM2

l1x l11x l01x

l2A

l1A

l12

C11 = C12

lA B l2B l1B

l1C

l1S

l1D C01

l11y l01y

l0

CM1

l1y = l01y + l11y

l0S l0D

l0C

O = C00 Figure F.2: Dimensions of the multibody system and the point mass sv parameters.

125

F.1 Kinematics of the pneumatic robot arm The body fixed vectors indicating the position of the revolute joints with respect to the center of mass are given by: £ ¤ ~b01 = −l01x −l01y ~e 1 £ ¤ ~b11 = (F.1) l11x l11y ~e 1 £ ¤ 2 ~b12 = 0 −l12 ~e . The generalized coordinates are given by: £ ¤T x1 y1 θ1 x2 y2 θ2 q = .

The variables x1 , y1 , x2 and y2 indicate the constraint vectors: £ ¤ x1 y1 ~e 0 ~c0 = £ ¤ x2 y2 ~e 1 , ~c1 =

(F.2)

(F.3)

where ~c0 is indicated in Figure F.1 (~c0 = l0~e20 ) and ~c1 is the relative position vector of C12 with respect to C11 (i.e. ~c1 = ~0). Consequently, the constraint equations related to the revolute joints are given by: q1 q2 q4 q5

= x1 = y1 = x2 = y2

= = = =

0 l0 0 0.

(F.4)

Applying the constraint equations to (F.2) leaves only two unconstrained generalized coordinates: £ ¤T (F.5) θ1 θ2 . q =

The direction cosine matrices, linking the three frames to each other are given by: ¸ · cos(θ1 ) sin(θ1 ) A10 = − sin(θ1 ) cos(θ1 ) · ¸ cos(θ2 ) sin(θ2 ) 21 A (F.6) = − sin(θ2 ) cos(θ2 ) · ¸ cos(θ1 + θ2 ) sin(θ1 + θ2 ) , A20 = A21 A10 = − sin(θ1 + θ2 ) cos(θ1 + θ2 ) where ~e 1

= A10 ~e 0

~e 2

= A20 ~e 0

= A21 A10 ~e 0 .

(F.7)

Now, the position of the center of mass for body 1 (CM1) is given by: ~rCM 1 r0CM 1

T = ~c0 − ~b01 = r0CM 1 ~e 0 £ ¤T £ ¤ 0 l0 = + l01x l01y A10 · ¸ l01x cos(θ1 ) − l01y sin(θ1 ) = . l0 + l01x sin(θ1 ) + l01y cos(θ1 )

126

(F.8)

Derivation of the multibody model

The position of the center of mass for body 2 (CM2) is given by: T = ~rCM 1 + ~b11 + ~c1 − ~b12 = r0CM 2 ~e 0 £ ¤ £ ¤ = r0CM 1 + l11x l11y A10 − 0 −l12 A20 · ¸ l1x cos(θ1 ) − l1y sin(θ1 ) − l12 sin(θ1 + θ2 ) = . l0 + l1x sin(θ1 ) + l1y cos(θ1 ) + l12 cos(θ1 + θ2 )

~rCM 2 r0CM 2

(F.9)

In order to define the velocity of a rigid body in space, the time derivative of the direction cosine matrices in (F.6) are required. 10 A˙ 21 A˙ 20

A˙

¸ −θ˙1 sin(θ1 ) θ˙1 cos(θ1 ) −θ˙1 cos(θ1 ) −θ˙1 sin(θ1 ) · ¸ −θ˙2 sin(θ2 ) θ˙2 cos(θ2 ) = −θ˙2 cos(θ2 ) −θ˙2 sin(θ2 ) =

·

21

(F.10)

10

= A˙ A + A A˙ · −(θ˙1 + θ˙2 ) sin(θ1 + θ2 ) = −(θ˙1 + θ˙2 ) cos(θ1 + θ2 ) 10

21

(θ˙1 + θ˙2 ) cos(θ1 + θ2 ) −(θ˙1 + θ˙2 ) sin(θ1 + θ2 )

¸

.

The angular velocity of a frame is given by the angular velocity vector ω ~ according 1 to ~e˙ = ω ~ × ~e 1 . The angular velocity vectors are given by: 10 21 20

ω ~ ω ~ ω ~

= = =

£

£

³

0 0

θ˙1

¤

¤

θ˙2 ´ θ˙1 + θ˙2 ~e30 , 0 0

~e 0 = θ˙1 ~e30 £ ¤ ~e 1 = 0 0 θ˙2 ~e 0 = θ˙2 ~e30

(F.11)

where ij ω ~ is the angular velocity vector of frame i with respect to frame j and ~e30 = ~e10 × ~e20 . The angular velocity vector can also be used to derive the time derivative of the direction cosine matrix by using the Poisson equation; this yields the same result as (F.10). The velocity of the center of mass CM1 of body 1 in space is derived by calculating the time derivative of the position vector: ~r˙CM 1

r˙ 0CM 1

³ ´ d ~ = − dt b01 =

d dt

³

b01 1

T

´

T 1 ~e 1 − b01 1 ~e˙

T T 10 = −b01 1 A˙ ~e 0 = r˙ 0CM 1 ~e 0 · ¸ −l01x θ˙1 sin(θ1 ) − l01y θ˙1 cos(θ1 ) = . l01x θ˙1 cos(θ1 ) − l01y θ˙1 sin(θ1 )

(F.12)

d (~c0 ) = 0 is valid as ~e 0 is a fixed reference frame and ~c0 is The simplification dt ¡ 1¢ d b01 = 0 as ~b01 is a body fixed a body fixed vector. The same can be said for dt vector of body 1.

127

The same approach goes for mass 2: ³ ´ d ~r˙CM 2 = dt ~rCM 1 + ~b11 − ~b12

T T 1 2 = ~r˙CM 1 + b11 1 ~e˙ − b12 2 ~e˙

(F.13) T T T 10 10 20 = −b01 1 A˙ ~e 0 + b11 1 A˙ ~e 0 − b12 2 A˙ ~e 0 ´ ³³ ´ 10 T T T T 20 ~e 0 = r˙ 0CM 2 ~e 0 = b11 1 − b01 1 A˙ − b12 2 A˙ · ¸ −l1x θ˙1 sin(θ1 ) − l1y θ˙1 cos(θ1 ) − l12 (θ˙1 + θ˙2 ) cos(θ1 + θ2 ) = . l1x θ˙1 cos(θ1 ) − l1y θ˙1 sin(θ1 ) − l12 (θ˙1 + θ˙2 ) sin(θ1 + θ2 )

r˙ 0CM 2

F.2 Mass and inertial properties Both bodies have a a steel and a aluminium part; this is taken into account in the calculation of the inertial properties. The inertia tensors of body 1 and body 2 are determined by modelling the robotic arm in the CAD package Unigraphics; in Appendix F.7, two output files are given. Each output file shows two types of inertia parameters: inertia in the revolute joint of the body, indicated by work; and inertia around the center of mass of the body, indicated by centroidal. The coordinates indicated by centroid give the location of the center of mass of each body, these parameters are used in Table F.1. The inertial parameters for body 1 and 2 are given by respectively Ixx and Iyy in their subsequent output file. Also the principal axes of inertia and the moments of inertia for the principal axes are calculated. The principal axes of body 1 do not coincide with the axes of ~e 1 . This is confirmed by the non-zero products of inertia. In case of body 2 , the principal axes of inertia do coincide with the axis of ~e 2 . This yields the inertia matrix to be a diagonal matrix: the three moments of inertia equal zero. The robotic arm allows only for rotation around the z axis (i.e. in the x − y plane, see Figure F.1). As a result only the Jzz inertial parameter is required. The angular momentum of both bodies is given by: ~ CM 1 H

= Jzz1 θ˙1 ~e30

~ CM 2 H

= Jzz2 (θ˙1 + θ˙2 ) ~e30 .

(F.14)

Besides the two bodies, a point mass sv is defined representing the four valves, two pressure sensors and the wiring which are mounted on body 1. In Figure F.2, point mass sv is indicated. The position and velocity vectors of sv are defined by: ~rsv ~bsv r0sv

r˙ 0sv

= ~c0 − ~bsv = r0sv ~e 0 £ ¤ −lsvx −lsvy ~e 1 = £ ¤T £ ¤ 0 l0 = + lsvx lsvy A10 · ¸ lsvx cos(θ1 ) − lsvy sin(θ1 ) = l0 + lsvx sin(θ1 ) + lsvy cos(θ1 ) · ¸ −lsvx θ˙1 sin(θ1 ) − lsvy θ˙1 cos(θ1 ) = . lsvx θ˙1 cos(θ1 ) − lsvy θ˙1 sin(θ1 ) 128

(F.15)

Derivation of the multibody model

All inertias and masses used in the multibody model are given by: msv m1 m2 Jzz1 Jzz2

= 0.80 kg = 2.64 kg = 0.533 kg = 89.3 · 10−3 kg m2 = 5.70 · 10−3 kg m2 .

(F.16)

The mass of the valves, sensors and wiring are determined experimentally using a balance. The mass of body 1 and 2 are derived from Appendix F.7. It should be noted that a number of masses are neglected. The mountings of the muscles consist of an aluminium part (36 g) and a plastic part (the muscle fitting, 10 g). The muscles are also neglected.

F.3 Spring On body 1, a spring is mounted to compensate for the momentum due to the effective mass unbalance around C01. In the initial position of the pneumatic robot arm (θ1 = θ2 = 0), the unloaded spring length lS must fulfill lS < l0 S . The vectors indicating the acting points of the spring are visualized in Figure F.3 and given by: ~b0S

=

~b1S

=

£

£

¤

~e 0

−l0S

0

−l1S

−l01y

¤

(F.17)

~e 1 .

The spring force vector F~S can be defined by the spring direction vector f~S , the unloaded spring length lS and the spring stiffness kS : F~S

³ ´ f~ S = kS kf~S k − lS , ~ kfS k

(F.18)

with f~S f 0S

T = ~rCM 1 + ~b1S − ~b0S = ~c0 − ~b0S − ~b01 + ~b1S = f 0S ~e 0 ¢T ¡£ ¤ £ ¤ l0S l0 + l01x − l1S 0 A10 = · ¸ l0S + (l01x − l1S ) cos(θ1 ) = . l0 + (l01x − l1S ) sin(θ1 )

The spring stiffness is experimentally determined on kS = 1 kN m−1 .

129

(F.19)

CM2

~b2A −F~A

~b1A

~b1S

° 6

~b1D

5 °

C01

F~A

3 °

CM1 ~b1C

 ° C11 = C12 2 ~b2B °

4 °

~b1B

−F~B F~B

−F~C −F~S

−F~D

F~C

F~D

F~S

~b0S

O = C00 ~b0D

~b0C

Figure F.3: Force vectors and related acting points of the multibody system.

130

Derivation of the multibody model

F.4 Muscle forces In Section 3.4, the behavior of the muscles is discussed. In the robot, four muscles are present. The body fixed vectors indicating the muscle mountings are visualized in Figure F.3. Using Figure F.2, the body fixed vectors are defined as: £ ¤ ~b0C = l 0 ~e 0 £ 0C ¤ ~b0D = −l0D 0 ~e 0 £ ¤ 1 ~b1A = l l ~e £ A B 1A ¤ 1 ~b1B = lA B −l1B ~e £ ¤ (F.20) ~b1C = −l1C −l01y ~e 1 £ ¤ 1 ~b1D = −l1D −l01y ~e £ ¤ ~b2A = 0 −l2A ~e 2 £ ¤ 2 ~b2B = 0 −l2B ~e .

The muscle forces acting in these points are given by the force vectors F~i , as indicated in Figure F.3: © ª ~ for i ∈ A B C D . F~i = Fi kff~i k (F.21) i

Herein, the muscle forces Fi are defined by (3.47), the force direction vectors f~i are given by: ³ ´ T f~A = ~rCM 2 + ~b2A − ~rCM 1 + ~b1A = ~b11 − ~b1A − ~b12 + ~b2A = f 0A ~e 0 f 0A

=

=

f~B f 0B

f 0C

f~D f 0D

·

l11x − lA B

l11y − l1A

¤

A10 +

£

0

l12 − l2A

¤

A20

¢T

(F.22) ¸ (l11x − lA B ) cos(θ1 ) − (l11y − l1A ) sin(θ1 ) − (l12 − l2A ) sin(θ1 + θ2 ) (l11x − lA B ) sin(θ1 ) + (l11y − l1A ) cos(θ1 ) + (l12 − l2A ) cos(θ1 + θ2 )

³ ´ T = ~rCM 2 + ~b2B − ~rCM 1 + ~b1B = ~b11 − ~b1B − ~b12 + ~b2B = f 0B ~e 0 =

=

f~C

¡£

¡£

·

l11x − lA B

l11y + l1B

¤

A10 +

£

0

l12 − l2B

¤

A20

¢T

(F.23) ¸ (l11x − lA B ) cos(θ1 ) − (l11y + l1B ) sin(θ1 ) − (l12 − l2B ) sin(θ1 + θ2 ) (l11x − lA B ) sin(θ1 ) + (l11y + l1B ) cos(θ1 ) + (l12 − l2B ) cos(θ1 + θ2 )

T = ~rCM 1 + ~b1C − ~b0C = ~c0 − ~b0C − ~b01 + ~b1C = f 0C ~e 0 ¢T ¡£ ¤ £ ¤ −l0C l0 + l01x − l1C A10 = · ¸ −l0C + (l01x − l1C ) cos(θ1 ) = l0 + (l01x − l1C ) sin(θ1 ) T = ~rCM 1 + ~b1D − ~b0D = ~c0 − ~b0D − ~b01 + ~b1D = f 0D ~e 0 ¢T ¡£ ¤ £ ¤ l0D l0 + l01x − l1D A10 = · ¸ l0D + (l01x − l1D ) cos(θ1 ) = . l0 + (l01x − l1D ) sin(θ1 )

131

(F.24)

(F.25)

An output from the robot dynamics to the muscle model is the actual muscle length hi (t) for i ∈ {A B C D }. This length is derived from the force direction vector f~i (t) as defined above. The force direction vector is the vector between the two muscle mounting points. This vector replaces the muscle and the fittings. The actual muscle length is now given by the absolute value of this vector, minus the fitting lengths: © ª hi (t) = kf~i (t)k − lfA B for i ∈ A B (F.26) © ª hi (t) = kf~i (t)k − lfCD for i ∈ C D .

F.5 Two DoF equation of motion

The necessary inputs to state the equations of motion for the pneumatic robot arm are stated. First, the energy equations will be formulated, thereafter the equations of motion will be derived using the Lagrangian approach. Kinetic energy The kinetic energy of the robot arm is given by: ³ ´ T = 21 msv ~r˙sv · ~r˙sv + m1 ~r˙CM 1 · ~r˙CM 1 + m2 ~r˙CM 2 · ~r˙CM 2 ³ ´ ~ CM 1 · 10 ω ~ CM 2 · 20 ω ~ + H ~ . + 21 H

Substituting (F.11) to (F.15) in (F.27) and reorganizing terms gives: ³ ´ ¢ ¢ ¡2 ¡2 2 2 T = 12 msv lsvx + Jzz1 θ˙12 + lsvy + m1 l01x + l01y ³ ´³ ´2 ¡ ¢ + 1 m2 l2 + l2 θ˙2 + 1 m2 l2 + Jzz2 θ˙1 + θ˙2 2

1x

1y

³

1

− m2 l12 θ˙1 θ˙1 + θ˙2

2

´³

12

(F.27)

(F.28)

´ l1x sin(θ2 ) − l1y cos(θ2 )

Potential energy The total potential energy equals the sum of the potential energy of the conservative forces and the internal energy. V

= V ex + U in .

(F.29)

The potential energy due to the gravity force is given by V ex

= −msv ~g · ~rsv − m1 ~g · ~rCM 1 − m2 ~g · ~rCM 2 ,

with the gravity vector defined by £ ¤ 0 −g ~e 0 . ~g =

(F.30) (F.31)

Application of (F.15), (F.8), (F.9) and (F.31) to (F.30) gives the gravitational energy: ³ ´ V ex = msv g l0 + lsvx sin(θ1 ) + lsvy cos(θ1 ) ³ ´ + m1 g l0 + l01x sin(θ1 ) + l01y cos(θ1 ) (F.32) ³ ´ + m2 g l0 + l1x sin(θ1 ) + l1y cos(θ1 ) + l12 cos(θ1 + θ2 ) . 132

Derivation of the multibody model

The internal energy available in the spring is given by ´2 1 ³ ~ kS kfS k − lS . U in = 2

(F.33)

Applying (F.19) to (F.33) gives the internal energy. ³³© ª2 l0S + (l01x − l1S ) cos(θ1 ) U in = 12 kS ´2 © ª2 ´ 21 − lS . + l0 − (l01x − l1S ) sin(θ1 )

(F.34)

Generalized muscle forces The remaining forces acting on the robotic arm are the muscle forces, as defined in (F.21). Each muscle consists of a pressure dependent length, stiffness and damping. The generalized muscle forces Q are defined by: ¶T nF µ X ∂~rj

Q =

j=1

∂q

· F~j .

(F.35)

In Figure F.3, the six circled numbers indicate the positions at which the muscle forces act, this means nF = 6 in (F.35). The circled numbers correspond with the number of the generalized muscle forces. The relation between the generalized muscle forces Qj for j ∈ {1 2 3 4 5 6} as defined in (F.35) and the muscles forces F~i for i ∈ {A B C D } as defined in (F.21), is also given by: Q1 Q2 Q3 Q4 Q5 Q6

→ → →

−F~A −F~B F~A F~B

→

−F~C → −F~D →

(F.36)

.

The generalized muscle force for j = 1 is given by: £ ¤ ~r1 = ~rCM 2 + ~b2A = ~rCM 2 + 0 −l2A A21 A10~e 0 · ¸ l1x cos(θ1 ) − l1y sin(θ1 ) − a1 0 r1 = l0 + l1x sin(θ1 ) + l1y cos(θ1 ) + b1 a1

=

b1

=

∂r01 ∂q Q1

(F.37) (F.38)

(l12 − l2A ) sin(θ1 + θ2 )

(l12 − l2A ) cos(θ1 + θ2 ) · −l1x sin(θ1 ) − l1y cos(θ1 ) − b1 = l1x cos(θ1 ) − l1y sin(θ1 ) − a1 = −

FA kf~A k

µ

∂r01 ∂q

¶T

−b1 −a1

¸

f 0A .

(F.39)

(F.40)

The derivation of Qi for i = 2 . . . 6 is given in Appendix F.8. 133

Equations of motion F X T d ³ ´ T,q˙ − T,q + V,q = Qnc i dt i=1

n

(F.41)

As all individual terms of the equations of motion are rather involved, the complete equations will not be stated. Each term is given below.  ³

       T,q˙ =        



          d ³ ´  T,q˙ =   dt           

´ ¡2 ¢ ¡2 ¢ 2 2 msv lsvx + lsvy + m1 l01x + l01y + Jzz1 θ˙1 ´³ ´ ³ ¡2 ¢ 2 2 + m2 l1x + l1y θ˙1 + m2 l12 + Jzz2 θ˙1 + θ˙2 ³ ´³ ´ − m2 l12 2θ˙1 + θ˙2 l1x sin(θ2 ) − l1y cos(θ2 ) ³

´³ ´ 2 m2 l12 + Jzz2 θ˙1 + θ˙2 ³ ´ − m2 l12 θ˙1 l1x sin(θ2 ) − l1y cos(θ2 )

T               

³

´ ¢ ¡2 ¢ ¡2 2 2 + Jzz1 θ¨1 msv lsvx + lsvy + m1 l01x + l01y ³ ´³ ´ ¡2 ¢ 2 2 + m2 l1x + l1y + Jzz2 θ¨1 + θ¨2 θ¨1 + m2 l12 ³ ´³ ´ − m2 l12 2θ¨1 + θ¨2 l1x sin(θ2 ) − l1y cos(θ2 ) ´³ ´ ³ − m2 l12 2θ˙1 + θ˙2 l1x θ˙2 cos(θ2 ) + l1y θ˙2 sin(θ2 ) ³

´³ ´ 2 m2 l12 + Jzz2 θ¨1 + θ¨2 ³ ´ − m2 l12 θ¨1 l1x sin(θ2 ) − l1y cos(θ2 ) ³ ´ − m2 l12 θ˙1 l1x θ˙2 cos(θ2 ) + l1y θ˙2 sin(θ2 ) 0

³ ´³ ´  − m2 l12 θ˙1 θ˙1 + θ˙2 l1x cos(θ2 ) + l1y sin(θ2 )



³ ´ msv g lsvx cos(θ1 ) − lsvy sin(θ1 ) ³ ´ + m1 g l01x cos(θ1 ) − l01y sin(θ1 ) ³ ´ + m2 g l1x cos(θ1 ) − l1y sin(θ1 ) − l12 sin(θ1 + θ2 ) ³ ´¡ ¢¡ ¢ − kS 1 − kfl~S k l01x − l1S l0S sin(θ1 ) + l0 cos(θ1 ) S

−m2 g l12 sin(θ1 + θ2 ) 134

T                       

T

T,q =        V,q =      

(F.42)

(F.43)

(F.44)

T            

(F.45)

Derivation of the multibody model

F.6 One DoF equation of motion In case the first joint at C01 is fixed at θ1 = 0, muscle C and D and the spring become redundant. The position vectors given in (F.8) and (F.9) can be simplified to the rather obvious relations:

~rCM 1

·

=

l01x l0 + l01y

·

¸T

~e 0

(F.46)

l1x − l12 sin(θ2 ) = l0 + l1y + l12 cos(θ2 ) · ¸ −l12 θ˙2 cos(θ2 ) = . −l12 θ˙2 sin(θ2 )

~rCM 2 r˙ 0CM 2

¸T

~e 0

(F.47) (F.48)

The muscle force vectors (F.22) and (F.23) can be rewritten to respectively:

f~A

=

f~B

=

·

·

l11x − lA B − (l12 − l2A ) sin(θ2 ) l11y − l1A + (l12 − l2A ) cos(θ2 ) l11x − lA B − (l12 − l2B ) sin(θ2 ) l11y + l1B + (l12 − l2B ) cos(θ2 )

¸T

¸T

~e 0

(F.49)

~e 0 .

(F.50)

Only two generalized muscle forces are present in the DoF system:

f 0A r01 ∂r01 ∂q Q1

=

·

l11x − lA B − (l12 − l2A ) sin(θ2 ) l11y − l1A + (l12 − l2A ) cos(θ2 )

=

·

l1x − (l12 − l2A ) sin(θ2 ) l0 + l1y + (l12 − l2A ) cos(θ2 )

=

·

−(l12 − l2A ) cos(θ2 ) −(l12 − l2A ) sin(θ2 )

= − =

FA kf~A k

µ

∂r01 ∂q

¶T

¸

f 0A

¸

¸

(F.51)

(F.52)

(F.53)

(F.54)

n¡ o ¢ ¡ ¢ FA (l12 − l2A ) l11y − l1A sin(θ2 ) + l11x − lA B cos(θ2 ) , kf~A k 135

·

l11x − lA B − (l12 − l2B ) sin(θ2 ) l11y + l1B + (l12 − l2B ) cos(θ2 )

f 0B

=

r02

=

·

l1x − (l12 − l2B ) sin(θ2 ) l0 + l1y + (l12 − l2B ) cos(θ2 )

∂r02 ∂q

=

·

−(l12 − l2B ) cos(θ2 ) −(l12 − l2B ) sin(θ2 )

Q2

FB = − kf~B k =

µ

∂r02 ∂q

¶T

¸

¸

¸

f 0B

(F.55)

(F.56)

(F.57)

(F.58)

o n¡ ¢ ¡ ¢ FB (l12 − l2B ) l11y + l1B sin(θ2 ) + l11x − lA B cos(θ2 ) . kf~B k

The terms in the Lagrange equation related to the kinematic and potential energy are simplified to: ¢ 1 ¡ 2 T = m2 l12 + Jzz2 θ˙22 (F.59) 2 ¡ ¢ 2 + Jzz2 θ˙2 (F.60) T,q˙ = m2 l12 d ³ ´ T,q˙ dt T,q

V ex U in V,q

=

=

¡

¢ 2 m2 l12 + Jzz2 θ¨2

0

¢ = m2 g l0 + l1y + l12 cos(θ2 )

=

0

¡

= −m2 g l12 sin(θ2 ) .

(F.61)

(F.62) (F.63) (F.64) (F.65)

Finally, the equation of motion for the DoF system is given by: 2 X d ³ ´ T,q˙ − T,q + V,q = Qj T dt j=1

¡

¢ 2 m2 l12 + Jzz2 θ¨2 − m2 g l12 sin(θ2 ) =

(F.66)

o n¡ ¢ ¡ ¢ FA (l12 − l2A ) l11y − l1A sin(θ2 ) + l11x − lA B cos(θ2 ) + (F.67) kf~A k o n¡ ¢ ¡ ¢ FB (l12 − l2B ) l11y + l1B sin(θ2 ) + l11x − lA B cos(θ2 ) . kf~B k 136

Derivation of the multibody model

F.7 Inertial parameters of both bodies Using Unigraphics, the mass, center of mass position and the inertial parameters are estimated. Two output files for respectively body 1 and 2 are given. The values indicated with Work give the inertial parameters with respect to the working point of the body. The working point is chosen according to the revolute joint the body rotates in. The values indicated with Centroidal give the inertial parameters with respect to the center of mass of the body. The inertial parameter of body 1 is given by Ixx = 89307997... ============================================================ Measurement Mass Properties Displayed Mass Property Volume = Area = Mass = Radius of Gyration = Centroid =

Values 555797.462344 mm^3 99441.694269 mm^2 2.638858 kg 184.085874 mm -0.000000, 66.694511, 13.749857 mm

============================================================ Detailed Mass Properties Analysis calculated using accuracy of 0.990000 Information Units Grams - Millimeters Density Volume Area Mass

= = = =

0.00474788 555797.46234368 99441.69426889 2638.85946627

First Moments Mxc, Myc, Mzc

= -0.00000000, 175997.44173566, 36283.94065550

Center of Mass Xcbar, Ycbar, Zcbar

= -0.00000000, 66.69451101, 13.74985713

Moments of Inertia (Work) Ixxw, Iyyw, Izzw = 101544959.62521537, 1909101.10982269, 99869140.25131625 Moments of Inertia (Centroidal) Ixx, Iyy, Izz = 89307997.30881111, 1410202.10958365, 88131076.93515104 Moments of Inertia (Spherical) I = 89424638.17677291 Products of Inertia (Work) Pyzw, Pxzw, Pxyw = 9492917.38440122, 0.02202123, 0.00000000 Products of Inertia (Centroidal) Pyz, Pxz, Pxy = 7072977.70474444, 0.02202123, 0.00000000 Radii of Gyration (Work) Rxgw, Rygw, Rzgw = 196.16478856, 26.89715427, 194.53937866 Radii of Gyration (Centroidal) Rxg, Ryg, Rzg = 183.96577916, 23.11705727, 182.74958837 Radii of Gyration (Spherical) R = 184.08587437 Principal Axes Xp(X), Xp(Y), Xp(Z) Yp(X), Yp(Y), Yp(Z) Zp(X), Zp(Y), Zp(Z)

= 1.00000000, 0.00000000, -0.00000004 = 0.00000004, -0.08076016, 0.99673356 = -0.00000000, -0.99673356, -0.08076016

Principal Moments Ixxp, Iyyp, Izzp

= 89307997.30881111, 88704163.69514164, 837115.34959304

============================================================

137

The inertial parameter of body 2 is given by Iyy = 5697561... ============================================================ Measurement Mass Properties Displayed Mass Property Volume = Area = Mass = Radius of Gyration = Centroid =

Values 118335.121191 mm^3 37374.670265 mm^2 0.533309 kg 104.107195 mm 0.000000, -0.000000, 175.950872 mm

============================================================ Detailed Mass Properties Analysis calculated using accuracy of 0.990000 Information Units Grams - Millimeters Density Volume Area Mass

= = = =

0.00450677 118335.12119085 37374.67026454 533.30934225

First Moments Mxc, Myc, Mzc

= 0.00000000, -0.00000000, 93836.24369690

Center of Mass Xcbar, Ycbar, Zcbar

= 0.00000000, -0.00000000, 175.95087178

Moments of Inertia (Work) Ixxw, Iyyw, Izzw = 22251745.53414509, 22208130.50291576, 121603.69416256 Moments of Inertia (Centroidal) Ixx, Iyy, Izz = 5741176.65130608, 5697561.62007675, 121603.69416256 Moments of Inertia (Spherical) I = 5780170.98277269 Products of Inertia (Work) Pyzw, Pxzw, Pxyw = -0.00000000, 0.00000000, -0.00000000 Products of Inertia (Centroidal) Pyz, Pxz, Pxy = 0.00000000, 0.00000000, -0.00000000 Radii of Gyration (Work) Rxgw, Rygw, Rzgw = 204.26428898, 204.06400440, 15.10023786 Radii of Gyration (Centroidal) Rxg, Ryg, Rzg = 103.75543588, 103.36057572, 15.10023786 Radii of Gyration (Spherical) R = 104.10719541 Principal Axes Xp(X), Xp(Y), Xp(Z) Yp(X), Yp(Y), Yp(Z) Zp(X), Zp(Y), Zp(Z)

= 1.00000000, 0.00000000, 0.00000000 = 0.00000000, 1.00000000, 0.00000000 = 0.00000000, 0.00000000, 1.00000000

Principal Moments Ixxp, Iyyp, Izzp

= 5741176.65130608, 5697561.62007675, 121603.69416256

============================================================

138

Derivation of the multibody model

F.8 Generalized muscle forces 2 to 6 Muscle force 2: ~r2

= ~rCM 2 + ~b2B = ~rCM 2 + ·

=

a2 b2

= (l12 − l2B ) sin(θ1 + θ2 ) = (l12 − l2B ) cos(θ1 + θ2 )

Q2

=

·

0

−l2B

l1x cos(θ1 ) − l1y sin(θ1 ) − a2 l0 + l1x sin(θ1 ) + l1y cos(θ1 ) + b2

r02

∂r02 ∂q

£

−l1x sin(θ1 ) − l1y cos(θ1 ) − b2 l1x cos(θ1 ) − l1y sin(θ1 ) − a2

FB = − kf~B k

µ

∂r02 ∂q

¶T

¤

A21 A10~e 0

(F.68a)

¸

−b2 −a2

(F.68b)

¸

(F.68c)

f 0B .

(F.68d)

Muscle force 3: £

¤

~r3

= ~rCM 1 + ~b1A = ~rCM 1 +

r03

=

·

(l01x + lA B ) cos(θ1 ) − (l01y + l1A ) sin(θ1 ) l0 + (l01x + lA B ) sin(θ1 ) + (l01y + l1A ) cos(θ1 )

∂r03 ∂q

=

·

−(l01x + lA B ) sin(θ1 ) − (l01y + l1A ) cos(θ1 ) 0 (l01x + lA B ) cos(θ1 ) − (l01y + l1A ) sin(θ1 ) 0

Q3

=

FA kf~A k

µ

∂r03 ∂q

¶T

lA B

l1A

A10~e 0

(F.69a) ¸ ¸

f 0A .

(F.69b)

(F.69c)

(F.69d)

Muscle force 4: ~r4

= ~rCM 1 + ~b1B = ~rCM 1 +

£

lA B

−l1B

¤

A10~e 0

=

·

(l01x + lA B ) cos(θ1 ) − (l01y − l1B ) sin(θ1 ) l0 + (l01x + lA B ) sin(θ1 ) + (l01y − l1B ) cos(θ1 )

∂r04 ∂q

=

·

−(l01x + lA B ) sin(θ1 ) − (l01y − l1B ) cos(θ1 ) 0 (l01x + lA B ) cos(θ1 ) − (l01y − l1B ) sin(θ1 ) 0

Q4

=

FB kf~B k

r04

µ

∂r04 ∂q

¶T

f 0B .

(F.70a) ¸ ¸

(F.70b)

(F.70c)

(F.70d)

139

Muscle force 5: £

~r5

= ~rCM 1 + ~b1C = ~rCM 1 +

r05

=

·

(l01x − l1C ) cos(θ1 ) l0 + (l01x − l1C ) sin(θ1 )

=

·

− (l01x − l1C ) sin(θ1 ) (l01x − l1C ) cos(θ1 )

∂r05 ∂q Q5

FC = − kf~C k = −

FC kf~C k

−l1C ¸

0 0

−l01y

¤

A10~e 0

(F.71b)

¸

(F.71c)

µ

∂r05 ∂q

·

¡ ¢ ¸ (l01x − l1C ) l0 cos(θ1 ) + l0C sin(θ1 ) . 0

¶T

(F.71a)

f 0C

(F.71d)

Muscle force 6: ~r6 r06 ∂r06 ∂q Q6

= ~rCM 1 + ~b1D = ~rCM 1 +

£

−l1D

=

·

(l01x − l1D ) cos(θ1 ) l0 + (l01x − l1D ) sin(θ1 )

=

·

− (l01x − l1D ) sin(θ1 ) (l01x − l1D ) cos(θ1 )

= −

FD kf~D k

FD = − kf~D k

¶T

0 0

¸

−l01y

¤

A10~e 0

¸

µ

∂r06 ∂q

·

¡ ¢ ¸ (l01x − l1D ) l0 cos(θ1 ) − l0D sin(θ1 ) . 0

(F.72a) (F.72b)

(F.72c)

f 0D

140

(F.72d)

G Derivation of the static muscle parameters

The identification of the muscles is performed in 40 working points, see Figure 4.4. A working point w = (hw , pw ) is defined by a length h which is imposed by the tensile bench and a pressure p which is imposed manually using a valve. This means that a working point is not necessarily an equilibrium of the muscle, i.e. the muscle might generate a force in a working point. Before each measurement, a working point w is imposed on the muscle. During ˜ ) is imposed on the length hw . The perturbathe measurement a perturbation h(f ˜ ) is a noise signal, containing all frequencies f between 0 and 20 Hz. In tion h(f order to access the low-frequency behavior of the muscle, each measurement takes 60 seconds. In each working point five measurements are performed; this yields 5 × 40 data sets for muscle 4. The five FRF’s in each working point are used to estimate the force over length FRF: Hw (f )

Fw (f ) , hw (f )

=

(G.1)

bw , which is derived using with Fw the muscle force. The FRF Hw is estimated by H the power spectra of the sensor signals, [de Kraker, 2000]: H1,w (f )

=

ΦhF,w (f ) Φhh,w (f )

H2,w (f )

=

ΦF F,w (f ) , ΦhF,w (f )

(G.2)

with ΦhF,w (f ) the cross-power spectrum of h and F , Φhh,w (f ) the auto-power spectrum of h and ΦF F,w (f ) the auto-power spectrum of F . By applying the Welch estimating technique on H1,w (f ) and H2,w (f ), Hw is estimated. bw , also the coherence γ 2 (f ) of each measurement is calculated using Besides H hF the power spectra. The coherence is a measure for the linearity of the input–output behavior in the measurement. The coherence can be assessed by: 2 γhF,w (f )

=

H1,w (f ) . H2,w (f )

(G.3)

bw frequency response function H

15

0

10

|Φ| [dB]

10

5 0

−2

10

−5

−4

−10

power spectra Φ

2

10

bw | [N m−1 dB] |H

20

0

10

1

10

10

2

Φhh [m2 ] ΦF F [N2 ] ΦhF [N m] 0

10

10

1

2

10

10

coherence 180

bw ) [deg] phase(H

1 0.8

2 γhF [–]

90

0

0.6 0.4

−90

0.2

−180 0

10

1

10

0

2

10

0

10

f [Hz]

1

2

10

10

f [Hz]

b w in w = (017, 0.8). The upper and lower left picture give the Figure G.1: FRF H b w . The power spectra are given in the upper right figure. The magnitude and phase of H coherence in the lower right figure shows that data up to 18 Hz is reliable.

bw is given in the working point w = (0.17, 0.8). It In Figure G.1, the estimate H bw (f ). shows three power spectra, the coherence and the phase and magnitude of H This analysis is repeated for the 40 working points of muscle 4, the resulting FRF’s are given in Figure G.2. Each analysis shows an excellent coherence up to f = 18 Hz. Remarkable is that the magnitude of the FRF increases slightly with the frequency. This phenomena, the so-called stiffening, is probably caused by friction and by the rubber behavior. Also note the non-zero phase at f = 0 Hz. It is assumed that the muscle behavior in a working point is given by a stiffness and viscous damping; the next step is to fit the stiffness kw and damping dw : Hw (f )

= kw + j 2πf dw ,

(G.4)

bw . The results of this fit for five working points with on te experimental data H pw = 0.9 bar,is given in Figure G.3. The stiffness is easy to plot and gives a reliable result. The damping is rather difficult to identify, as it is hardly present in the measurements (see also Figure G.1). The non-zero phase at f ≈ 0 Hz and the slight increase in |Hw | due to the stiffening complicate the damping fit further. Nevertheless, (G.4) is used to derive kw and dw and is believed to give an accurate and meaningful value, especially on the stiffness. Finally, the stiffness kw and damping dw are given as a function of the working points w in Figure G.4 and as a surface plot in Figure G.5. Besides the experiments on muscle 4, also eight reference measurements are performed on four other muscles, see Figure 4.4. The results of these measurements are given in Figure 4.5. 142

Derivation of the static muscle parameters

Comparing the stiffness and damping models with the measurements In Section 4.1, the identified stiffness is fitted on the muscle stiffness model given in (4.1). The relative error between the analytical stiffness model and the identified stiffness, is given in the left part of Figure G.6. In the right part, the relative error between the experimental stiffness model and the identified stiffness is given. The relative error of the experimental model has a maximum of 5%, while the relative error of the analytical model is up to 60%. In Figure G.7, the fit and the relative error of the experimental damping and the experimental damping model is given. The relative error has a maximum of 30%. From this figure, it can be concluded the experiments do not give an clear view on the damping and that the damping in the muscles is almost constant and low. For this reason, the model is considered to be satisfactory.

15

10 5 0

0

5

10

15

160 mm

15 10 5 0

0

5

10

15

20

10

10

5 0

5

10

15

20

180 mm 190 mm

0

5

10

15

20

0

5

10

15

20

0

5

10

15

20

0

5

10

15

20

0

5

10

15

20

0

5

0

5

10

15

20

0

15

15

10

10

5

5 0

5

10

15

20

0

15

15

10

10

5

5

0

20

10

5

0

15

15

10

0

10

5

15

200 mm

0 15

0

5

5

15

bw | [N m−1 dB] |H

170 mm

0.5 bar 10 0.6 bar 0.7 bar 5 0.8 bar 0 200.9 bar 0 1.0 bar 15 1.2 bar 10 1.4 bar

bw ) [deg] phase(H

150 mm

15

0

5

10

15

20

f [Hz]

0

f [Hz]

Figure G.2: The magnitude and phase of the 40 measured FRF’s from muscle 4. Each b w , for one muscle length. The figures on the right figure on the left gives the magnitude of H give the phase for one muscle length. The color of each line, refers to a pressure level. As b w | increases with the pressure and with the muscle length. can be seen, the magnitude |H

143

150 mm

6.5

20 d = 15.7 Ns/m

6 10 5.5 k = 5.69 kN/m 5

0

5

10

15

0

20

7

160 mm

6

15

20

10

15

20

10

15

20

10

15

20

10

15

20

10 k = 5.93 kN/m 5

10

15

170 mm

6 k = 6.45 kN/m 0

5

10

15

0

20

20

8.5

phase(Hw ) [deg]

0

7

5

180 mm

10

20

8

0

5

20 d = 14.1 Ns/m 10

0

0

5

20 d = 14.7 Ns/m

8 10 7.5 k = 7.71 kN/m 7

0

5

10

15

20

11

190 mm

5

d = 14.8 Ns/m

5

|Hw | [N m−1 dB]

0

0

0

5

20 d = 16.6 Ns/m

10.5 10 10 k = 10.2 kN/m 9.5

0

5

10

15

20

f [Hz]

0

0

5

f [Hz]

b w in five working Figure G.3: The magenta lines are the experimentally estimated FRF H points w = (hw , 0.9) with hw ∈ {0.15, . . . , 0.19} m. The black lines represent Hw with the fitted stiffness kw and damping dw as defined in (G.4). The left figures give the magnitude and the right figures the subsequent phase. The fitted stiffness kw is given in the magnitude plot and the damping dw in the phase plot. As the bandwidth of the robotic arm is not likely to exceed 10 Hz, only the experimental data up to 10 Hz is considered to fit the stiffness and damping. This hardly affects the value of kw , but the damping is not very well estimated using this kind of experiment.

144

16

34

14

30

dw [Ns m−1 ]

kw [kN m−1 ]

Derivation of the static muscle parameters

12 10 8 6

0.5 bar 0.6 bar 0.7 bar 0.8 bar 0.9 bar 1.0 bar 1.2 bar 1.4 bar

26 22 18 14

4 2 150

160

170

180

190

10 150

200

160

170

180

190

200

hw [mm]

16

34

14

30

dw [Ns m−1 ]

kw [kN m−1 ]

hw [mm]

12 10 8 6

150 mm 160 mm 170 mm 180 mm 190 mm 200 mm

26 22 18 14

4 2 0.4

0.6

0.8

1

1.2

10 0.4

1.4

0.6

0.8

pw [bar]

1

1.2

1.4

pw [bar]

Figure G.4: The stiffness kw and damping dw as function of the pressure pw and the length hw of muscle 4. Remark that the use of colors in the two upper figures indicate a pressure while in the two lower figures the use of colors indicates the muscle length.

20

24

dw [N s m−1 ]

kw [kN m−1 ]

22 15

10

5

0 1.7

20 18 16 14 12 1.7

1.4

0.24 0.22

1.1 0.8

pw [bar]

0.14 0.2

0.12

0.24 0.22 0.2

0.8

0.18 0.16

0.5

1.4 1.1

0.2

pw [bar]

hw [m]

0.18 0.16

0.5

0.14 0.2

0.12

hw [m]

Figure G.5: Surface plot of the results given in Figure G.4. It can be seen that the stiffness shows a quadratic dependency on the length, while it shows a linear dependency on the pressure. This is in accordance with the stiffness model given in (4.1). The damping does not show a clear dependency on either the pressure or the length.

145

1.4

1.3

1.3

1.2

1.2

1.1

1.1

p [bar]

p [bar]

1.4

1 0.9

1 0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5 0.15

0.16

0.17

0.18

0.19

0.5 0.15

0.2

0.16

0.17

h [m] −0.6

−0.55

−0.5

−0.45

0.18

0.19

0.2

h [m]

−0.4

−0.35

−0.3

−0.04

−0.02

0

0.02

0.04

0.06

Figure G.6: The relative error between the analytical stiffness model and the identified stiffness (left figure) and between the experimental stiffness model and the identified stiffness (right figure). The match between the models and the measurements is given in Figure 4.6.

1.4 1.3 24

1.2 1.1 20

p [bar]

dw [N s m−1 ]

22

18 16

1 0.9 0.8

14 12 1.7

0.7 1.4

0.24 0.22

1.1

pw [bar]

0.6

0.2

0.8

0.18 0.16

0.5

0.14 0.2

0.12

0.5 0.15

hw [m]

0.16

0.17

0.18

0.19

0.2

0.1

0.2

h [m] −0.3

−0.2

−0.1

0

Figure G.7: The damping as derived from the muscle experiments, is given by the solid surface in the left figure. The experiments do not give a clear picture of the damping, as function of the pressure and muscle length. This is mainly due to the type of experiment. The damping model as defined in (3.49), is fitted on the data and shown by the mesh surface in the left figure. In the right figure, the relative error between the identified and the experimental damping model is given.

146

H Tuning of the controllers

This appendix presents additional information on the design of controllers. First, the implementation of the controllers using filters is given. In appendix H.2, several frequency response functions of the pressure behavior Pp are given. In Appendix H.1, the stability of the DoF pressure controller Cp is proven, using the Nichols diagram. In the remaining parts of this appendix, the identification and stability of the 2DoF robotic arm are discussed.

H.1 Controller part definitions The linear controllers consists of a gain with one or more additional filters. Below, the general purpose and implementation of the five used filters is given. See also [Franklin, 2001]. Gain The gain of the feedback loop, with s ∈ C: P (s) = p .

(H.1)

Lead-lag filter A lead-lag filter is used to gain phase lead within a certain frequency range. A lead filter is also called a tame d-action: 1 s+1 2πf1 L(s) = . 1 s+1 2πf2

(H.2)

Notch filter A notch filter is used to eliminate certain frequencies for the contoller. A pure notch filter has equal frequency parameters fn1 = fn2 , while a tilted notch filter has unequal frequency parameters fn1 6= fn2 . The depth of the notch filter is defined by ββ12 ; the width is related to the relative damping of the filter (which is defined by β2 ): 1 2β1 s+1 s2 + (2πfn1 )2 2πfn1 . N (s) = 2β2 1 2 s + s + 1 (2πfn2 )2 2πfn2

(H.3)

First order low-pass filter Eliminate frequencies above fl : F1 (s) =

1 . 1 s+1 2πfl

(H.4)

Second order low-pass filter Eliminate frequencies above fl . Depending on the parameter βl , the filter is more or less sharp: F2 (s) =

1 . 2βl 1 2 s + s+1 (2πfl )2 2πfl

(H.5)

The controllers as defined in this work, all consist of a gain P (s), combined with one or more filters. If for example a controller consisting of a gain P (s) with leadlag filter L(s), notch filter N (s), the implementation is given by defining the filters according to (H.1), (H.2) and (H.3); the filters are transformed to the z-domain and implementing in S IMULINK blocks. These blocks are applied to the error signal one by one in a S IMULINK model, the model is build and implemented in the dSPACE.

H.2 Frequency response functions of the pressure loop frequency response function Pp =

|Pp | [bar V−1 dB]

20

pA B zp

10 0

−10 −20 −30

θ 0 p 0.8 θ 0 p 0.8 θ 0 p 1.0 θ 0 p 1.0 θ 0 p 1.2 θ 0 p 1.2 θ 0 p 1.4

−40

0

1

10

10

phase(Pp ) [deg]

180

90

0

−90

−180 0

1

10

10

f [Hz] pA B zp

Figure H.1: The FRF of Pp = in four equilibria with varying p¯A B and θ2 = 0 deg. In the legend, the equilibria are given short handed. As can be seen, additional experiments are performed to identify the behavior below 1 Hz. All FRF’s are approximately identical.

148

Tuning of the controllers

H.3 Nichols diagram of the open-loop pressure control 30

20

|HOL | [dB]

10

0

−10

−20

−30 −360

−270

−180

−90

0

90

phase(HOL ) [deg] Figure H.2: The Nichols diagram of the open-loop HOL = Cp Pp , with Pp = pzApB . The phase and gain margin are indicated by the circle around the origin (−180, 0). As can be seen, the controller yields a stable closed loop system.

149

H.4 Tuning of the 1DoF angle controller for joint 1 Before the 2DoF robotic arm is controlled, we aim to control joint 1 in a similar way as presented for joint 2 . In this section, the identification of Pθ around joint 1 and the design of Cθ for joint 1 is presented. During the identification and control of θ1 , muscles A and B are constantly fully inflated. This yields the angle of joint 2 , θ2 , to be 0 deg, with a maximum stiffness of muscles A and B . This way, the dynamic coupling between revolute joints 1 and 2 is minimized. The frequency response functions of Pθ = zθθ1 around several working points are indicated in green in Figure H.3. Using these FRF’s a controller Cθ is designed consisting of a gain with a lead-lag, notch and second-order low-pass filter. The parameters of this controller are given in Table H.1. The bandwidth of the controlled system is approximately 0.5 and 1 Hz. This is indicated by the vertical pink lines in Figure H.3. In this figure, also the sensitivity of the closed-loop system is given. In Figure H.4, the stability of the closed loop system is proven in a Nichols diagram. Using these figures, it can be seen that all robustness margins are met. The pressure controller Cp is designed using a gain with a first-order low-pass filter, such that the bandwidth is approximately 1 Hz, similar to the bandwidth of θ1 .

gain p

lead-lag 6

f1 f2

0.1 6

notch fn β1 β2

1.4 0.1 3.5

low-pass fl βl

6 0.8

Table H.1: The Cθ controller parameters for joint 1 of the 1DoF robotic arm. The controller consists of a gain with a notch, lead-lag and a second-order low-pass filter. All frequencies f are given in Hz.

150

Tuning of the controllers

40 20

|H| [dB]

0 −20 −40

Cθ [V bar−1 ] HOL [–] Pθ [bar V−1 ]

−60 −80

−1

0

10

10

phase(H) [deg]

180

90

0

−90

−180 −1

0

10

10

20 10

|S| [dB]

0 −10 −20 −30 −40

−1

0

10

10

f [Hz] Figure H.3: The two upper figures show the magnitude and phase of the controller Cθ , plant Pθ = zθθ1 and the open-loop HOL = Cθ Pθ for joint 1. The lower image shows the sensitivity S = (1 + Cθ Pθ )−1 . The sensitivity margin of 6 dB is indicated by the pink line. The bandwidth is inbetween 0.5 and 1 Hz, as indicated by the two vertical pink lines.

151

30

20

|HOL | [dB]

10

0

−10

−20

−30 −360

−270

−180

−90

0

phase(HOL ) [deg] Figure H.4: The Nichols diagram of the open-loop HOL = Cθ Pθ for joint 1, with Pθ = zθθ1 . The phase and gain margin are indicated by the circle around the origin (−180, 0). As can be seen, the controller yields a stable closed-loop system.

152

90

Tuning of the controllers

H.5 Dual joint control scheme The control strategy for the 2DoF robotic arm R2 is shown in Figure H.5. The implementation shows much resemblance with the control strategy for R1 as presented in Section 5.1 and Figure 5.2. The figure below is not complete, the coupling between the different inputs and outputs are not shown. It is believed that pA and pC are not coupled. As a result, it is also assumed that the coupling between θ1 and p¯C and between θ1 and p¯C is absent. The coupling between θ1 and θ2 can be neglected as is shown in Appendix H.7. Finally, the coupling between θ1 and p¯A and between θ2 and p¯C is taken for granted as it is believed that this will not result in stability problems.

plant P 2

+

eθ1

θ1ref

+

eθ2

θ2ref

+

epA

pCref

+

epC

-

z

Ψ

x R2

PWM

pC

PpC

v

u

pA

PpA

zpC

CpC

θ2

Pθ2

zpA

CpA

θ1

Pθ1

zθ 2

Cθ2

pAref

zθ 1

Cθ1

y Ω

Figure H.5: The embedding of the 2DoF robotic arm R2 in the control loop, yielding the plant P 2 . The input matrix Ψ and the output matrix Ω are defined in Section 5.3.1.

153

H.6 MIMO system identification In a MIMO system identification, the plant is identified by determining the input and process sensitivity experimentally. The theory of MIMO system analysis and control is treated in [Tousain, 2006; Skogestad, 1996]. In Figure H.6, the input-output relations of a 2DoF MIMO system are given. Also the three types of signals are indicated, the outputs θ, the noise n and the error signals e. The controllers, frequency response functions and signals are stored in matrices: ¸ · Cθ1 0 (H.6) Cθ = 0 Cθ2 · ¸ Pθ11 Pθ12 Pθ = (H.7) Pθ21 Pθ22 n = e = θ

=

£ £

£

n1

n2

e1

e2

θ1

θ2

¤T

¤T

¤T

(H.8) (H.9) .

(H.10)

The input sensitivity SI is derived from e and n by: e

=

(I + C θ P θ )

−1

n

SI

=

(I + C θ P θ )

−1

,

(H.11)

while the process sensitivity SP is derived using θ and n according to: θ SP

= P θ (I + C θ P θ )

−1

n

−1

= P θ (I + C θ P θ )

.

(H.12)

Multiplying the process and the inverse input sensitivity yields the plant P θ : Pθ

= SP S−1 I .

(H.13)

Essential in the experimental identification is that the noise signals are not injected simultaneous. First, noise in n1 is injected and e(1) and θ(1) are measured. Using the power spectra Φ of n1 , e(1) and θ(1) , the input and process sensitivity terms depending on n1 are calculated for each frequency ωi of the power spectra:   Φe(1) n1 (ωi ) 1 . . .  Φn (ωi )  1   (H.14) SI (ωi ) =    Φ (1) (ωi )  e2 n1 ... Φn1 (ωi )

SP (ωi )



Φθ(1) n1 (ωi ) 1  Φn (ωi ) 1  =   Φ (1) (ωi ) θ2 n1 Φn1 (ωi )

... ... 154



   . 

(H.15)

Tuning of the controllers

−

+

Cθ1

Pθ11

+

+

θ1 +

e1

n1

Pθ12 Pθ Pθ21 n2 +

Cθ2

−

e2 Pθ22

+

+

θ2

+

Figure H.6: The input-output relations as present in a coupled 2DoF MIMO system. The identification is performed in closed-loop, using the controllers Cθ1 and Cθ2 The direct frequency response functions are given by Pθ11 and Pθ22 , the frequency response functions determining the coupling are given by Pθ12 and Pθ21 . Noise is injected in either n1 or n2 ; the error signals e1 and e2 and the output signals θ1 and θ2 are measured.

Next, noise is injected in n2 and e(2) and θ(2) are measured. The parts of the input and process sensitivity depending on n2 are calculated, using the power spectra Φ of n2 , e(2) and θ(2) :   Φe(1) n1 (ωi ) Φe(2) n2 (ωi ) 1 1  Φn (ωi ) Φn2 (ωi )  1   SI (ωi ) =  (H.16)   Φ (1) (ωi ) Φ (2) (ωi )  e2 n1 e2 n2 Φn1 (ωi ) Φn2 (ωi )

SP (ωi )



Φθ(1) n1 (ωi ) 1  Φn (ωi ) 1  =   Φ (1) (ωi ) θ2 n1 Φn1 (ωi )

 Φθ(2) n2 (ωi ) 1 Φn2 (ωi )    . Φθ(2) n2 (ωi )  2 Φn2 (ωi )

(H.17)

Finally, the plant P θ is estimated by a frequency-wise matrix manipulation of the process and input sensitivity: P θ (ωi ) = SP (ωi ) S−1 I (ωi ) .

(H.18)

155

H.7 Interaction between the joint angles A measure for the interaction between two coupled outputs, is the interaction coefficient KI which is defined by: KI (ω) =

H12 H21 , H11 H22

(H.19)

with H = Pθ , as defined in Figure 5.8. In Figure H.7, the interaction coefficient for coupling between the two joint angles is given. As can be seen, is the interaction below −20 dB up to 4 Hz. So it is legitimate to consider the two joints as decoupled within the bandwidth of the controlled robotic arm.

20 10

p¯C = 0.4 bar p¯C = 0.6 bar p¯C = 0.8 bar

0 −10

|KI | [dB]

−20 −30 −40 −50 −60 −70 −80

−1

10

0

10

1

10

f [Hz] Figure H.7: The four frequency response functions of P θ for different values of p¯C . The FRF’s are experimentally determined around θ¯1 = θ¯2 = 0 deg. The pressure in muscle A is kept at p¯A = 0.4 bar in closed loop.

156

Tuning of the controllers

H.8 Stability of the dual joint controller 30

20

|HOL | [dB]

10

0

−10

−20 −270

−180

−90

0

phase(HOL ) [deg] Figure H.8: The Nichols diagram of the open-loop HOL = Cθ1 Pθ1 for joint 1, with Pθ1 = zθθ1 . The phase and gain margin are indicated by the circle around the origin 1

(−180, 0). As can be seen, Cθ1 yields a stable closed-loop system.

157

30

20

|HOL | [dB]

10

0

−10

−20 −270

−180

−90

phase(HOL ) [deg] Figure H.9: The Nichols diagram of the open-loop HOL = Cθ2 Pθ2 for joint 2 , with Pθ2 = zθθ2 . The phase and gain margin are indicated by the circle around the origin 2

(−180, 0). As can be seen, Cθ2 yields a stable closed-loop system.

158

0

I Implementation of the robotic arm model

Main file for simulating the robotic arm model % Model of pneumatics and one muscle % Master’s project Stef van den Brink % Philips AppTech %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Set simulation parameters fpwm = 50 fs = 1000 fvalve = 300 tend = 20 Tpwm = 1/fpwm Ts = 1/fs

;% ;% ;% ;% ;% ;%

% Set simulation parameters [air,beta] = presparameters par = nominalbehav Pin = [1e5 1]

;% load pneumatic parameters ;% load robot and muscleparameters ;% [Pa] air supply and outer pressure

sample frequency PWM controller frequency Valve frequency limit (acceptable error band) end time of simulink simulation pwm sample time controller sample time

% Initial condition and excitation x0 = [1e5 .5e5 NaN 0 0] ;% [Pa Pa Pa deg rad/s] Starting conditions x0(4) = x0(4)*(pi/180) ;% x0(3) = equilibria(x0(4),x0(2),par) ;% %------------------------------------------------------------------% Solve ODE for pressure t=0:Ts:tend; x=ode3(@PneuTwoMuscleModelDV,t,x0,air,beta,par,fpwm,Pin,exc); % Recap system values n=length(t); Q=zeros(n,6); HH=Q; vv=zeros(n,5); KK=zeros(n,4); FF=KK; U=zeros(n,2); DD=U; VV=U; PS=zeros(n,1); for i=1:n [xdot,q,u,v,Ps,K,F,H,D,V]=PneuTwoMuscleModelDV(t(i),x(i,:),air,beta,par,fpwm,Pin,exc); % [hA,hB,dhA,dhB]=kinematics(x(i,4),x(i,5)); Q(i,:)=q; U(i,:)=u; PS(i,:)=Ps; HH(i,:)=H; vv(i,:)=v; DD(i,:)=D; VV(i,:)=V; KK(i,:)=K; FF(i,:)=F; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Setting the pneumatic parameters function [air,beta]=presparameters; % Pneumatic model parameters % Master’s project Stef van den Brink % Philips AppTech %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Constants & Air properties (Polytechnisch zakboekje table 1.23 on page E1/24) air.frac = [78.09 20.95 0.93 0.03] ;% [-] fractions N2 O2 Ar CO2 air.M_frac = [28.013 31.999 39.948 44.01] ;% [kg/kmol] fractions N2 O2 Ar CO2 air.M = air.frac*air.M_frac’/sum(air.frac) ;% [kg/kmol] molecular mass air air.R = 8314.51 ;% [J/(kmol K)] universal gas constant air.Rs = air.R/air.M ;% [J/(kg K)] specific gas constant air air.T = 293 ;% [K] temperature air.eta = 17.1e-6 ;% [Pa s] dynamic viscosity air %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Robot pneumatics parameters R.vonce = 7.2e-4 R.vtwice = 5.2e-4 R.reduc = 4e-4

;% [m] radius of one valve ;% [m] radius of two serial valves ;% [m] radius of reduction valve

beta.reduc = 0.59*sqrt(2)*pi*R.reduc^2 beta.vsingle = 0.59*sqrt(2)*pi*R.vonce^2 beta.vdouble = 0.59*sqrt(2)*pi*R.vtwice^2

;% [m^2] restriction gain reduction valve ;% [m^2] restriction gain valves ;% [m^2] restriction gain valves

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Setting the robotic arm parameters function par=nominalbehav % Determine nominal behavior of muscle % Master’s project Stef van den Brink % Philips AppTech % % par=nominalbehav equil = 1 comp = 1 obs = 0

;% [0 1] calculate equilibrium points ;% [0 1] compare equilibria with experiments ;% [0 1] observability (comp=1 required)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Set parameters par.Vc par.m par.J par.L12 par.LfA par.LfB par.d1 par.n par.L par.s1 par.s2

= = = = = = = = = = =

10e-3 0.533 5.7e-3 0.176 0.14 par.LfA 0.025 1.25 0.1274 0.117 0.232e5

;% ;% ;% ;% ;% ;% ;% ;% ;% ;% ;%

[m^3] capacity volume (10 l) [kg] 0.533 [kg m^2] 5.697561e-3 [m] 0.17595 [m] length of muscle A fitting [m] length of muscle B fitting [m] fitting diameter [-] number of windings of one thread [m] thread length [m] muscle parameter 1 [Pa] muscle parameter 2

%------------------------------------------------------------------% Stiffness and damping [kd,mdl] = FitMuscleModel %kd = mdl

;% load passive parameters ;% use modeled stiffness

par.k = [kd.p0*kd.k0*kd.c0 kd.p0*kd.k1*kd.c1 kd.p0*kd.k2... kd.p1*kd.k0 kd.p1*kd.k1 kd.p1*kd.k2]’; par.d = [40 kd.d1]’; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

160

Implementation of the robotic arm model

Derive the passive muscle parameters from the experimental identification function [kd,mdl]=FitMuscleModel % The function FitMuscleModel fits the stiffness % model on experimentally identified muscle stiffness. % Master’s project Stef van den Brink % Philips AppTech % % Usage: [kd,mdl]=FitMuscleModel %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % load measured data if ~exist(’S’,’var’); load([base,’\MuscleFit\S10.mat’]); end %-------------------------------------------------------------------------% set figure options if plot is required if nargout==0; dis=1 ; clc else dis=0; end %-------------------------------------------------------------------------for z = 1:3 % select data for least squares fit and store in columns %-------------------------------------------------------------------------if z3; chck=varargin{4}; else; chck=0; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Determine equilibrium %%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

162

Implementation of the robotic arm model

syms pB [hA hB dhA dhB A B] = kinematics(phi,0,par); KA1 KB1 KA2 KB2 KA KB Fz h0A h0B

= = = = = = = = =

[1 hA hA^2] * par.k(1:3) [1 hB hB^2] * par.k(1:3) [pA pA*hA pA*hA^2] * par.k(4:6)*1e-5 [pB pB*hB pB*hB^2] * par.k(4:6)*1e-5 KA1+KA2 KB1+KB2 par.m*9.81*par.L12*sin(phi) (2*par.L-(par.s1*pA)/(par.s2+pA)) (2*par.L-(par.s1*pB)/(par.s2+pB))

;% ;% ;% ;% ;% ;% ;% ;% ;%

[N/m] [N/m] [N/m] [N/m] [N/m] [N/m] [N] [m] [m]

Z1 = A*KA*(hA-h0A) + B*KB*(hB-h0B) + Fz; Y=eval(solve(Z1,pB)); pB=Y(2); if imag(pB)~=0 | pB1e5; ButtonName=questdlg([’Pressure in muscle B = ’,n2s(pB*1e-5,3),’ bar, continue?’], ... ’Unrealistic output’, ... ’Yes’,’No’,’No’); switch ButtonName, case ’Yes’; disp([’Pressure in muscle B = ’,n2s(pB*1e-5,3),’ bar’]) case ’No’; error(’Simulation aborted by user.’) end end end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Damping and stiffness of the muscles %%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if nargout uth xoutA = u(1) < -uth

;% filter very small values ;% filter values of almost 1

;% [kg/s] control signal from zero order hold ;% u(1) block function ; u(2) sinus function ;% [N] no external excitation ;% [kg/s] no valve signal ;% [N] external force excitation

;% muscle A air in signal ;% muscle A air out signal

165

xinB = u(2) > uth xoutB = u(2) < -uth

;% muscle B air in signal ;% muscle B air out signal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Volume change %%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [hA hB dhdtA dhdtB A B] = kinematics(x(4),x(5),par) ;% apply kinematic relations to get muscle variables h0A d2A VA dVdhA

= = = =

2*par.L - par.s1*x(2)/(par.s2+x(2)) ;% [m] nominal length of muscle A sqrt(4*par.L^2-hA^2)/(par.n*pi) ;% [m] actual diameter d2 for muscle A pi/60*hA*(3*par.d1^2+4*par.d1*d2A+8*d2A^2) ;% [m^3] actual volume of muscle A pi/20*par.d1^2 + ... par.d1/15/par.n*(4*par.L^2-hA^2)^(1/2) + ... 8/15/par.n^2/pi*par.L^2 - 2/5/par.n^2/pi*hA^2 - ... par.d1/15*hA^2*par.d1/par.n/(4*par.L^2-hA^2)^(1/2) ;% [m^2] volume change per length change

h0B d2B VB dVdhB

= = = =

2*par.L - par.s1*x(3)/(par.s2+x(3)) ;% [m] nominal length of muscle B sqrt(4*par.L^2-hB^2)/(par.n*pi) ;% [m] actual diameter d2 for muscle B pi/60*hB*(3*par.d1^2+4*par.d1*d2B+8*d2B^2) ;% [m^3] actual volume of muscle B pi/20*par.d1^2 + ... par.d1/15/par.n*(4*par.L^2-hB^2)^(1/2) + ... 8/15/par.n^2/pi*par.L^2 - 2/5/par.n^2/pi*hB^2 - ... par.d1/15*hB^2*par.d1/par.n/(4*par.L^2-hB^2)^(1/2) ;% [m^2] volume change per length change

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Pneumatic model for two muscles and capacity %%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Pressure in supply line Ps = [1-xinA xinA]*[x(1) x(3) ; x(2) min(x(2:3))]*[1-xinB ; xinB]; % Air density rho = ([P0(1)+x(1) ; x(1)+Ps ; x(2)+P0(2) ; x(3)+P0(2)]/2 + 1.01325)/(air.Rs*air.T); % Pressure change in volumes dP = [(P0(1) - x(1)) (x(1) - Ps) (xoutA*(x(2) - P0(2))) (xoutB*(x(3) - P0(2))) ]; % Flow in hoses qin = beta.reduc * qs = beta.vdouble qsA = qs * (xinA * qsB = qs * (xinB * qoutA = beta.vsingle qoutB = beta.vsingle

sqrt(rho(1)*abs(dP(1))) * sign(dP(1)) * sqrt(rho(2)*abs(dP(2))) * sign(dP(2)) (x(1)-x(2))) / (x(1) - x(2) + xinB*(x(1)-x(3))) (x(1)-x(3))) / (xinA*(x(1)-x(2)) + x(1) - x(3)) * sqrt(rho(3)*abs(dP(3))) * sign(dP(3)) * sqrt(rho(4)*abs(dP(4))) * sign(dP(4))

;% ;% ;% ;% ;% ;%

to capacity to valves to muscle A to muscle B out muslce A out muscle B

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Damping and stiffness of the muscles %%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% pA = x(2)*1e-5; pB = x(3)*1e-5; kA dA kB dB

= = = =

[1 [1 [1 [1

hA hA^2 pA pA] * par.d hB hB^2 pB pB] * par.d

pA*hA pA*hA^2] * par.k pB*hB pB*hB^2] * par.k

;% ;% ;% ;%

[N/m] [Ns/m] [N/m] [Ns/m]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Forces acting on the robot dynamics %%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FA = kA*(hA-h0A) + dA*dhdtA FB = kB*(hB-h0B) + dB*dhdtB

;% [N] force generated by muscle A ;% [N] force generated by muscle B

QA = FA*A QB = FB*B Fz = par.m*9.81*par.L12*sin(x(4))

;% [N] Non-conservative force of muscle A ;% [N] Non-conservative force of muscle B ;% [N] gravity force on the beam

166

Implementation of the robotic arm model

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Differential equation for capacity and muscle pressure %%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if nargout