Robotics and Autonomous Systems 63 (2015) 136–149
Contents lists available at ScienceDirect
Robotics and Autonomous Systems journal homepage: www.elsevier.com/locate/robot
A unified system identification approach for a class of pneumatically-driven soft actuators Xiaochen Wang ∗ , Tao Geng, Yahya Elsayed, Chakravarthini Saaj, Constantina Lekakou Faculty of Engineering and Physical Sciences, University of Surrey, Guildford, GU2 7XH, United Kingdom
highlights • A unified system identification approach is proposed. • The approach is used to identify the nonlinear pressure–shape dynamic relation. • The used auxiliary kinematic setting can be implemented by gyroscopic sensors.
article
info
Article history: Received 16 December 2013 Received in revised form 10 August 2014 Accepted 29 August 2014 Available online 16 September 2014 Keywords: Continuum robot Soft actuator System identification Orthonormal basis function
abstract The class of Pneumatically-driven Low-pressure Soft Actuators (PLSA) is a popular choice potentially used in the surgical robotic applications. One fundamental problem lying in the PLSA research is the lack of a generally validated model for the complex nonlinear dynamic behaviours. In this paper, a unified identification approach for the general PLSAs is proposed. It is a parameter-independent way directly used to identify the dynamical relation between the actuating pressures and the principal degrees of freedom of a PLSA, the bending and the steering. The approach is based on a modified auxiliary kinematic setting and a newly developed identification model structure, named DIO–PWL–OBF. Following the concluded identification procedure, the implementations for the single chamber bending and the double chamber bending and steering are demonstrated separately. The results show that the proposed approach can accurately capture the nonlinear pressure–shape dynamical relation. The approach is also efficient in realtime applications. It can be further used to improve the current control design for the PLSAs in robotic applications. Crown Copyright © 2014 Published by Elsevier B.V. All rights reserved.
1. Introduction In the fields of surgery and therapy, soft robotic instruments begin to draw researchers’ attention in the past 20 years [1,2]. Compared with traditional rigid ones, the soft robots can potentially give a more viable solution when the operation is taken in human bodies. The inherent compliance offers human friendly interaction; the kinematic redundancy provides the manoeuvre competency in highly unstructured in-body organ environments. One promising choice is the class of Pneumatically-driven Lowpressure Soft Actuators (PLSA). They are compactly designed in
∗ Correspondence to: Surrey Space Centre, Faculty of Engineering and Physical Sciences, University of Surrey, Guildford, GU2 7XH, United Kingdom. Tel.: +44 7760635800. E-mail addresses:
[email protected] (X. Wang),
[email protected] (T. Geng),
[email protected] (Y. Elsayed),
[email protected] (C. Saaj),
[email protected] (C. Lekakou). http://dx.doi.org/10.1016/j.robot.2014.08.017 0921-8890/Crown Copyright © 2014 Published by Elsevier B.V. All rights reserved.
small sizes and safely operated at low pressures. An immediately foreseeable application is endoscopes [3]. It can overperform the current rigid endoscopes by giving the doctors increased adroitness meanwhile reducing damage on the patients. Fig. 1 shows three typical PLSA examples. In this paper, we specifically use the term PLSA to distinguish it from the Pneumatic Muscle Actuators (PMA) like in [4]. The latter commonly have much larger sizes and higher operation pressures. But both the PLSAs and the PMAs are in the class of the pneumatically-driven soft continuum robots, and they share many similar properties. For the general pneumatic soft continuum robots, the kinematic modelling is still the backbone in the state-of-the-art control designs [6–10], which is based on a purely geometrical description of the actuators’ shapes. Completely validated dynamic modelling for a single section of the pneumatic soft robot arms is lagging behind. The potential applicability is therefore restricted. The modelling has been analytically studied in [11,12]. The Lagrange formation under the constant curvature condition was used to derive the dynamic relation between the individual chamber lengths and the
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
137
Fig. 1. Featured PLSA prototypes: (a) Flexible Micro-Actuator by [5]. (b) Colobot by [6]. (c) An early design of the STIFF-FLOP actuator.
shape. The worked models are ineffective in handling the coupling effect between each chamber. Unlike the PMAs, this internal effect cannot be neglected in the PLSA cases, because the latter are structured in one integral rubber material. The radial expansion of each individual chamber is not constrained as the way by the individual bellow suits used in a PMA. Besides, the above modellings do not include the complex nonlinear dynamic relation between the actuating pressures and the chamber lengths, which is an independent open problem too. A recent work by [13] started to tackle the issue of lacking the actuator dynamics based on the previous work in [12]. The proposed model correctly accounts for the identified hysteretic behaviour using the Bouc–Wen restoring force method. This research is on the right track, but it is mainly investigated for the PMAs. Due to the structural difference, the involved intermediate step of the separate pressure–elongation identification for a PMA’s individual bladders cannot be applied to the internal chambers of a PLSA. And the internal coupling issue has not been formally studied. Besides, another work by [14] directly used a second-order linear transfer function to approximate the pressures–shape relation. But only the single chamber bending situation was partially studied. The model is only valid for a small degree of bending. The deformation was measured based on the deflection distance projected onto the horizontal base plane. However this kinematic variable does not naturally represent the bending deformation in a linear manner. The primary contribution of this paper is a unified identification approach for the general PLSAs. For a single segment PLSA, it can practically give reliable approximations for the nonlinear dynamical relation between the actuating pressures and the two principal Degrees Of Freedom (DOF), the bending and the steering. The identification approach has three features: 1. It is based on the measurement of the tip surface orientation of a single segment PLSA. This auxiliary kinematic setting is incorporated into the mainstream kinematic framework and can be easily implemented by on-shelf gyroscopic technologies. 2. It is based on a newly developed model structure, of which parameters are statistically defined. Hence specific mechanical parameters such as the material elasticity, the chamber placement and the chamber dimensions are not required. This makes the modelling capable of covering general designs of the PLSAs. 3. It preserves characteristics of local dynamics through linear approximation; no weighting term is used to fit the observed nonlinearities of PLSAs. The resultant model structure becomes a type of switched linear systems. Mature linear control techniques can be used in later control design. The paper is arranged in the following way: the structure and the nonlinearities of the PLSAs are discussed in Section 2. It states the need to build a new suitable identification model structure. Section 3 studies the auxiliary kinematic setting regarding Feature 1. Some aspects in the mainstream kinematic configuration are discussed in Section 3.1 firstly. And then the new setting is introduced in Sections 3.2 and 3.3. Regarding Features 2 and
3, the unified identification approach is developed in Section 4. General concerns and fundamental theories are briefly discussed and introduced in Sections 4.1 and 4.2. Section 4.3 studies the new proposed identification model structure in detail and summarises the unified identification procedure at the end. Section 5 shows the implementation and results. The conclusion and future work are presented in Section 6. 2. The PLSAs 2.1. The structure of a single segment PLSA The PLSAs commonly have the 1-, 2-, 3- and 4-chamber designs. Among them, the 3-chamber design is the most versatile type for robotic applications. It has all three DOFs compared with the 1- and the 2-chamber ones, and at the same time it uses one less chamber in contrast to the 4-chamber one. Our work will focus on the general 3-chamber PLSA design as sketched in Fig. 2(a). Optimally it is designed in a very symmetric way. Three identical pneumatic chambers are regularly disposed at 120° apart, and they are together parallel to the central longitudinal axis of the actuator. Fig. 2(b) gives a brief explanation for the ideal bending. The outer covering of a PLSA usually wear a kind of bellow suits so that the radial expansion induced by the fed-in pressure is largely reduced and can be neglected. The individual chamber elongations are the main type of the deformations structured by the bellow. Thereby they deterministically drive the deformation of the entire actuator. The related DOFs are classified as bending, steering, and elongating. As the case in Fig. 2(b), Chamber 2 takes the longest incremental length 1lc2 as it is subjected to a higher pressure than Chambers 1 and 3, which makes Chamber 2 lead the bending direction. 2.2. The observed nonlinearities Fig. 3 shows a typical example of the single chamber deflection1 —actuation relation from our experiments. The test consists of 28 local step responses, each of which starts the initial condition from its predecessor’s final condition. And both the pressure increase and the decrease situations are taken. The first nonlinearity type that can be spotted is the deadzone effect at the initial pressure range. It is due to the gas fillin phase during the initial actuation. Secondly, in the pressure increase cases, it can be seen that the steady state values of each local step response are not in a linear relation with their input ones. The similar things happen in the pressure decrease cases. By further comparing both the increase and the decrease cases, it is not hard to figure out an underlying hysteresis loop. The cause of
1 The deflection can be seen as the complementary motion of the bending. And they are equivalent under the constant curvature condition. In the rest of the paper, we will use the two words interchangeably when no confusion is caused.
138
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
a
b
Fig. 2. (a) A simplified structural drawing for the PLSAs at inactive state. (b) The illustration of the bending forced by actuation.
Fig. 3. Globally conducted local step responses.
the static nonlinearity roots in the nonlinear elongation–pressure relation of a chamber or bladder. Especially the hysteresis effect is due to the friction with the braid [2,15,16]. At last, we extracted the data sets of the increase cases, removed their offsets on both the input and output channels, and then put them together aligned to a same trigger time as shown in Fig. 4. The post-analysis showed that all those data sets can be very well linearised, but each linearised model is formed in an inconsistent model structure with others and they behave in different transient responses. They have very similar increments at the input channel. But as indicated from the output channel, the slopes under the local responses became increasingly steep as the operating pressure rose step by step. A similar result was also obtained in the decrease cases. In run time, varying the fed-in pressure alters the actuator’s stiffness, which is the main reason that causes the dynamical nonlinearity here. On the other hand, the analysed static and the dynamical nonlinearities also inherit in the cases when operating more than one chamber. Additionally, new types of nonlinearities are introduced due to the coupling effect between the active chambers.
The state-of-the-art nonlinear identification model structures might appear to be deficient to describe the combined static and dynamical nonlinearites observed in Figs. 3 and 4. We tested the Hammerstein–Wiener (HW) and the nonlinear Auto-Regression with eXogenous variables (ARX) identification approaches. And by them, an Input–Output (IO) sequence sampled from a PLSA experiment could be well fitted. But the identified models could not be fully validated when applying them to different scenarios. The nonlinear estimators in the HW and the nonlinear ARX model structures incorporate with a fixed linear part [17,18]. The structure settings of this kind might not be capable of closely capturing the varied dynamical proprieties particularly observed in Fig. 4. The intelligent methods such as the artificial neural network are adoptable, but the implementations are usually computationally costly. And their further extensions to the control design is needed to be differently reconsidered. As the fact that the local data sets have already shown strong linearity, we considered to use them to reconstruct the global
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
139
Fig. 4. The extracted local step responses from Fig. 3.
nonlinear dynamical relation directly. In this way, the nonlinear estimators with the nonlinear weighting on the regressors are not used, and modelling error caused by the inconsistent nonlinear mapping between the input and the output channels will therefore be avoided. Alternatively, the local linearised dynamics is preserved and transferred to the global sense. Besides, the local linearity idea is potentially compatible with the mature linear robust control techniques. The two concerns stated above led to the motivation behind the later proposed new model identification structure for the PLSAs in Section 4.3.
l c1
l c3
l
a
c
l2
3. Kinematic modelling 3.1. The current kinematic setting In the mainstream approaches, the kinematic modelling of the PLSAs is studied under the constant curvature condition where the actuators or robots are approximated to one or many pieces of constant curvature arc. Fig. 5 illustrates such kinematic configuration for a single section PLSA. The constant curvature condition has also been generally held in various types of the continuum robots. Regarding this, a comprehensive review is presented in [19]. In their work, a unified constant curvature kinematic setting has been generalised. It considers three spaces, actuator, configuration and task, with the connections by two mappings, robot-independent and robot-dependent. The robot-independent mapping is set between the configuration space by the curvature κ = 1r , the plane angle φ and the arc length la , and the task space by the 3-dimension position, x, y and z. It can be applied to a wide range of continuum robots. The transformation between them can be calculated via the conventional Denavit–Hartenberg approach which is commonly used in the cases of the rigid manipulators. The robot-dependent mapping bridges the actuator space to the configuration space. It is specified by the actuation mechanism and the designed structure. As the PLSAs configuration shown in Fig. 5, the actuator space parameters are lengths of the three chambers lc1,2,3 . The robot-specific mapping with respect to Fig. 5 is described as
κ=
2
2
d
√ φ = arctan la =
2
2 lc1 + lc2 + lc3 − lc1 lc2 − lc2 lc3 − lc1 lc3 3(lc2 + lc3 − 2lc1 )
lc1 + lc2 + lc3 3
3(lc3 − lc2 )
,
,
(1)
,
(2)
(3)
Fig. 5. Kinematic parameters for a single-section PLSA.
Fig. 6. The kinematic framework for the PLSA actuators.
where d is the average distance between the actuator coaxial centre and each of the three chambers coaxial centres. The overall constant curvature kinematic setting for the single segment threechamber PLSAs is illustrated in Fig. 6. Using the robot-dependent mapping, (1)–(3), this setting is adopted by nearly all of the similar PLSA work [5,6,9,11,12]. In this way, the DOFs of the bending and the steering do not come from a straightforward measurement, but estimated from the actuator space or from the task space. From the experimental exercises, we found that the bending and the steering estimations in the configuration space might be improved by adding an auxiliary setting which can bring a direct measurement to these two important types of the shape information. The new setting is based on the equivalence between the tip surface orientation and the shape of a single segment PLSA. It will be introduced in the next two subsections.
140
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
Configuration subspace 1
Configuration subspace 2
Configuration space
Task space
Configuration subspace 3
Fig. 9. The modified auxiliary kinematic framework for the PLSA actuators.
Fig. 7. The illustration for the three configuration subspaces 1, 2 and 3. They are sectioned by the three chambers’ individual bending directions.
Fig. 8. The bending angle θ b is equivalent to the deflection angle θ d in the constant curvature condition.
3.2. The auxiliary kinematic setting The auxiliary setting is based on the tip surface orientation of a PLSA. It is currently implemented by the TrackstarTM magnetic sensor in our research, but this method can be applied to general gyroscopic sensors. This paper focuses on the one-chamber and the two-chamber scenarios. Three-chamber actuation will lead to kinematic redundancy and variable stuffiness of a PLSA. This is an open and active area of research. The extension to the threechamber case will be studied together with a ready stiffening mechanism in the future. An overall view of the new auxiliary setting is illustrated in Fig. 7. The modification takes place at the robot-dependent mapping part; the part of the robot-independent mapping is kept unchanged. The actuator space is replaced by three configuration subspaces as shown in Fig. 7. In this setting, the scenario of the actuation by any two certain chambers is uniquely confined to each of the subspaces. The three configuration subspaces are defined on the same polar plane that is for the configuration space; the latter is divided by the former. Each of the three configuration subspaces is a fan-shaped area, and its two radial boundaries are defined where the two chambers’ singly-driven deflections project onto the polar plane. The three sectors are symmetric by the conceptual definition. But due to many practical factors such as mechanical designs or fabrication, they are likely to be different. The three configuration subspaces can be calibrated according to the actual situations. Therefore this reality gap wont affect the later system identification work.
Every individual configuration subspace consists of a deflection angle θi and a steering angle φisub to represent the deflection and the steering. As the equivalence to the tip surface orientation, the term arc length la is not used to estimate the bending and the steering. Compared with the configuration space, the parameters for representing the DOFs of the bending are modified in the configuration subspaces. The deflection angle is directly used. But as shown in Fig. 8, the deflection angle θ d is equal to the bending angle θ b under the constant curvature condition. Hence the former can be used to approximate the latter. For the DOF of the steering, the plane angle φ is divided and separately brought to each of the configuration subspaces as φ1sub , φ2sub and φ3sub . The only change made to them is the reduced range which now is the intersection angle of two projected beams by the single chamber deflection. Taking the three subspace plane angles φ1sub , φ2sub and φ3sub to the configuration space simply needs to add different offsets to each of them. In run time, switching between the subspaces can be implemented under the state machine framework. The problem due to the physically adjoined upper and lower mathematical bounds in the mainstream method [19] is thereby avoided. As to the new setting, Fig. 9 shows the modified kinematic framework updated from Fig. 6. The simple setting in Fig. 9 is effective to capture the DOFs of the bending and the steering. Each subspace individually links to a deformation scenario played by two certain chambers. In this way, the identification work can be taken separately according to the three subspaces. It is a practical approach when the deformations by every two chambers become greatly unsymmetrical due to different reasons such as designs, fabrication or usage. And it is supported by the current sensor technologies. On the other hand, the new auxiliary setting can well incorporate with the current elongation measurement to estimate the task space parameters. A PLSA’s length can be accurately measured by the optical fibre sensors as in [6], or be estimated by the polynomially pre-identified pressure–length relation. With an acquired tube length la and the bending and the steering angles θ b and φ obtained from this proposed auxiliary setting, one can use the summarised constant curvature approach in [19] to easily estimate the task space parameters x, y and z. 3.3. The calculation for the deflection and the steering angles The calculation for the deflection angle and the steering angle is presented in this subsection. The configuration is based on the TrackstarTM magnetic sensor which has a sensing base and a sensing tip. The base provides the global reference frame G, and it is fixed at a place during operation. The sensing tip is attached with the local reference frame B, and it is installed perpendicular to the tip surface of a PLSA. The measurement setting is illustrated in Fig. 10, which can be adopted by the gyroscopic sensors with minor adjustments.
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
Fig. 10. Measurement setting at the tip end of a single section PLSA. The local reference frame is based on the TrackstarTM magnetic sensor.
The measurement of the deflection angle θ d is illustrated in Fig. 10. The degree of the deflection is equivalent to the intersection angle between BG eZ and B ex which can be calculated as
θ = arccos(⟨ ex , d
B
B G eZ
⟩)
(4)
B G eZ
AG = Ax,γ Ay,β Az ,α
(5)
where α , β and γ are the rotations about the local z-, y- and x-axis respectively. α , β and γ are sampled by the magnetic sensor. In detail, (5) is expressed as cβcα −c γ + c α s γ sα sγ + c α sβ c γ
c β sα c α c γ + sβ sα sγ −c α sγ + sβ c γ sα
−s β c β sγ cβcγ
(6)
where s and c denote the sin and cos functions. Substituting (5) into (4), the executable form is acquired as
θ d = arccos(⟨B ex , B AG G eZ ⟩) T where B ex = 1 0 0 and G eZ = 0
OB ⊥ ABC ,
OB′ ⊥ AB′ C .
This leads to
=
OC OB
As ̸ COB + ̸ COB′ = 90°, we get cos ̸ AOB′
sin ̸ COB
= tan ̸ COB. AOB cos ̸ COB Hence, the steering angle can be expressed as cos ̸
=
(7)
,
OB′
AB AB OC AB′ and correspondingly
=
OC OB′ AB′ OC
cos ̸ AOB = cos ̸ COB cos ̸ AOC , cos ̸ AOB′ = cos ̸ COB′ cos ̸ AOC .
φ = arctan
⟨B ex , BG eY ⟩ . ⟨B ex , BG eX ⟩
(8)
Taking the transformation of (6) into (8), the steering angle can be calculated as
⟨B ex , B AG G eY ⟩ ⟨B ex , B AG G eX ⟩ T where B ex = 1 0 0 , G eX = 1 T G eY = 0 1 0 . φ = arctan
(9) 0
T
0
and
Compared with (2) for steering angle, (9) does not use the term lc3 − lc2 , i.e. the difference between two chambers’ lengths. The latter is sensitive to disturbed sensor reading for chamber lengths, while the former has no such a problem. 4. The model identification setting
T
0 1 . Compared with (1), (7) does not require the distance parameter d, which depends on different PLSA designs and can be varied under the internal actuating pressures. Bending measurement by (7) is from the outside of a PLSA. Therefore it is more robust than (1) during the shape deformation. The measurement of the steering angle φ is also illustrated in Fig. 10. It is the intersection angle between the projection of B ex onto the OXY plane and X -axis, or alternatively Y -axis. Note that because of the similarity between φisub and φ , only the symbol φ is used here. Further specification regarding the configuration subspace is neglected. The geometric relation between φ , B ex , BG eX and BG eY is equivalently presented by a pyramid in Fig. 11. From the relation, it can be immediately obtained
OB
Fig. 11. The pyramid that structures the geometric relation between φ , B ex , BG eX and B G eY .
B
where is a unit vector of Z -axis in B with respect to G, and ex is a unit vector of x-axis in B. To transform G eZ to BG eZ , the local frame yaw–pitch–roll rotation matrix is used according to the sensorspecific triple rotation. It is expressed as B
141
4.1. The experimental design The system identification approach is proposed in the bottomup fashion. The local data sets are identified at first and they will be combined together to reconstruct the global responses. The motivation comes from the fact that each local data set can be fairly linearised. Using them would lead to a globally accurate dynamical description. One can refer to Figs. 15 and 16 for the single chamber bending identification test design, and Figs. 19 and 20 for the bending and the steering one. Before the test, a number of pressure levels on the input channel are selected. They are set as the reference points for partitioning the local data sets. For various designs of the PLSAs, densely partitioned local data sets are usually unnecessary, as they react in slow dynamical responses. The bandwidth is low for the PMAs [20], and specifically for the PLSAs, it is normally below 20 Hz. During the identification test, two local step responses in an increase and decrease pair between each two neighbour reference points are separately arranged. The increase and the decrease mean that in a local interval the step response will be further individually identified regarding the mutually inverted valuechanging directions between the two reference points. This can give an effective solution to the hysteresis effect as observed in Fig. 3.
142
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
In the two-chamber bending and steering cases, the system has two inputs and two outputs. The identification test is arranged in a decoupled way, but the main idea is still similar to the single chamber cases. Each time only one of the two input channels is allowed to excite a local step response, meanwhile the responses on the two output channel are being sampled. In this way, the Single-InputMultiple-Output (SIMO) models can be separately identified with respect to each input channel. They will be used to approximate the complete Multiple-Input-Multiple-Output (MIMO) model in a local interval. 4.2. Fundamental theories A new identification model structure will be developed in the next subsection. It is based on Orthonormal Basis Function (OBF) which will be briefly introduced below. A set of orthonormal basis function Bi (q) span an n-dimensional subspace. A stable and strictly proper transfer function G(q) can be approximated in the subspace as G(q) ≈
n
ci Bi (q)
(10)
identification in [24]. Based on the OBF-based LTI system setting of (10), the LPV–OBF model structure can be expressed as yk ≈
1, Bi1 (q), Bi2 (q) = 0,
if i1 = i2 if i1 ̸= i2
(11)
Bi (q) can be constructed recursively [21] as
Bi (q) =
i −1 1 − |ξi |2 1 − ξk q
q − ξi
q − ξk
k=0
(12)
through a repetitive pole structure [22]
{ξ1 , . . . , ξnb , ξnb +1 , . . . , ξnb +nb , . . . , }.
(13)
(13) results in an optimal selection when the poles coincide with the actual system’s ones [21]. In our work, the pole set is determined from the identified local systems. The pole selection problem is formulated as in [22] nb z − ξj min max 1 − zξ ξ1 ,...,ξn z ∈Ωp b
j
j =1
(14)
where ξi s are basic poles in(13), Ωp contains the eigenvalues from x −y
all local systems and 1−xy¯ is a metric for two complex numbers x and y. The solution is provided by a recent study in [23] where a clustering algorithm named Fuzzy–Kolmogorov c-Max (FKcM) is developed.
4.3. The DIO–PWL–OBF model structure In the light of the local approach with the designed experimental procedure, a new nonlinear model structure specified for the PLSA identification is going to be developed. It is named as Delta Input–Output PieceWise Linear Orthonormal Basis Function (DIO–PWL–OBF) model structure as the indication of the three involved theoretical constituents. Along the line of reasoning, firstly the general modification on the LPV–OBF model structure will be discussed, and it will then be brought to the context of the PieceWise Linear (PWL) systems; secondly the new ingredient Delta Input–Output (DIO) will be introduced, and eventually the specific development by taking the DIO to arrive at the DIO–PWL–OBF form will be presented. The first OBF-based adaptive model structure is the LPV–OBF one which has been profoundly developed for the LPV system
(15)
i =1
where p is the scheduling signal with respect to the frozen LTI systems and is the scheduling operator. In the behavioural perspective [25], the term (ci p) can be seen as the scheduling part which on-line updates the coefficient ci of Bi (q) according to p. The linearisation type in the LPV identification is based on a floating reference point, whereas our linearisation method relies on the number of fixed equilibrium points. The discrepancy lies in the different conceptual formulation of the scheduling part, but it does not affect the role of the OBF part. The general idea of the coefficient-adaptive OBF structure is preserved. We take (15) to the category of the piecewise linear systems which is generally expressed as yk ≈ Gl (q)uk ,
l ∈ NLoc .
(16)
Thus the model structure is modified to
i=1
where θi is the coefficient for each Bi (q). The orthonormal property means
n (ci p)Bi (q)uk
yk ≈
n
cil Bi (q)uk ,
l ∈ NLoc ,
(17)
i =1
which is temporarily named as the PWL–OBF structure. In the PWL context, the OBF selection objective is redirected for suiting our application: to optimally approximate the partitioned local data sets. The underlying pole selection problem for the related LTI systems remains the same. The OBF coefficients are varied based on the localised regions. The region that the system is currently staying at decides the related OBF coefficients currently being used. There are two issues regarding applying the OBF structure to the PWL setting. They are termed as switching effect and offset effect. The brief discussions about them below are based on the ordinary difference equations (ODE) in the standard form as na
ai q−i yk =
i=0
nb
bj q−j uk .
(18)
j =1
The switching effect happens during the transition from one local region to another. If the post local model structure hardly follows the previous one in a consistent manner like in a same order or having similar frequency characters, a relatively large modelling error will be caused. As the fixed structure dominator n in each individual OBFs, the synthesised LTI i=1 ci Bi (q)uk results na −i in a fixed transience part i ai q . Due to this, the switching effect should have been small. But because of the different local dynamical n characters, the embodied Finite Impulse Response (FIR) parts j=b 0 bj q−j on the input side commonly result in inconsistent parametrisation. This can easily enlarge the modelling error during the transition. Especially, when the transience part has a pole near the unit circle, the convergence rate of the transition modelling error will be relatively slow if the switching action on the FIR part causes a relatively large deviation. On the other hand, the offset effect is caused by the initial and final conditions in each sampled local data set. In the designed experiment, the system starts and ends at two equilibrium points in each of the partitioned local regions. This will lead to two steadystate simultaneous equations in the linear regression process as
na nb a y = bj uss1 i ss 1 i =0
j =1
nb na a y = bj uss2 i ss 2 i =0
j =1
(19)
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
where uss1 and yss1 are the input and output at the initial condition; uss2 and yss2 are at the final condition. The global static nonlinearity has been observed in Fig. 3. As a result, the inequality relation by yss1
̸=
yss2
uss1
(20)
uss2
takes over the equality relation, where the latter is only held in a globally linear sense. This fact will let the coefficient sets {ai } and {bj } extend to a relatively long size to fulfil the condition made by (19). Consequently the linearised model results in a relatively highorder model structure. In order to overcome the two issues discussed in the temporary PWL–OBF model structure, a specific but simple development is further conducted. Firstly, the offsets on both input and output channels are removed so that (19) can be combined to one equation as na
ai (yss2 − yss1 ) =
na
bj (uss2 − uss1 ).
(21)
j =0
i=0
In this way, the problem caused by the global static nonlinearity are eliminated. The local data sets can be accurately identified by a relatively low-order configuration. But after this step, the model only works in the local sense. This means such a model cannot be directly applied to the PWL setting as in (16). In parallel, our consideration reverts to the superposition principle and the convolution sum which are the foundations to evaluate the dynamics of an ODE (18). Based on the input and output time series, {uk } and {yk }, the Delta Input and Output (DIO) time series are defined as
δ yk , yk − yk−1 , δ uk , uk − uk−1 ,
δ y1 , y1 δ u1 , u1 .
(22)
It can be seen that {δ uk } and {δ yk } represent the immediate change between every two neighbour sampling points in {uk } and {yk }. The relation between the former and the latter is uk =
k
δ un ,
k
yk =
n =1
δ ym .
(23)
m=1
By substituting (23) into the standard ODE form (18), it can directly achieve that na
ai q−i
i=0
k m=1
δ ym =
nb
bi q−j
j =1
k
δ un .
(24)
n =1
It is the ODE in the DIO sense. Simplify the synthesised operators k k k−i k−j q−i m=1 and q−j n=1 on both sides of (24) to m=1 and n=1 respectively, and expand them. (24) becomes a0
k
δ y m + · · · + a1
1
= b1
k−i
δ ym
new perspective induced by (26), the two calculation rules can be applied in a separate and independent way. That is, the dynamic responses of each δ uk can be individually evaluated. In the LTI systems, this feature is redundant, since it makes no difference as putting δ uk and δ yk back into uk and yk by (23) and then doing the same evaluation according to (18). But in the nonlinear systems, this feature is very appreciated because it gives an extra degree of freedom to manipulate the dynamics. Inspired by the reasoning above, the outcome (26) is going to be taken back to the temporary PWL–OBF model structure (16). In (26), preserve the equation in the first line for the latest sample δ uk and δ yk , and use (18) to replace the rest of the equations for the time series {δ uk−1 } and {δ yk−1 }, a reduced expression can obtained as Eq. (27) in Box II which is equivalently in the transfer-function form yk = G(q)q−1 uk + G(q)δ uk
1
δ un + · · · + b j
(28)
in the term of the LTI transfer function G(q). Moreover, (28) can be more generally adopted to the nonlinear systems by simply replacing G(q)q−1 uk with yk−1 , which then becomes yk = yk−1 + G(q)δ uk . (29) In run time, the one-step-before output sample yk−1 can be directly acquired from the feedback or be estimated by iterations from the previous input samples {uk−1 }. The derived (29) is particularly suitable for integrating the PWL model structure (16) with our designed local identification approach as yk ≈ yk−1 + Gl (q)δ uk ,
l ∈ NLoc ,
(30)
because the offset issue of the locally identified models Gl (q) is solved. Each identified local model Gl (q) can be used in (30) without any modification in evaluating a global dynamic response. Besides, (30) also indicates an alternative approach to handle the hysteresis effect. Between two neighbour equilibrium points, the local models are separately identified based on the increase and decrease directions. Compared with the HW and the nonlinear ARX approaches, this manner would be more advantageous for the PLSA identification, because all the local properties of the dynamics are preserved. Finally the OBF ingredient is added to (30). It now becomes yk ≈ yk−1 +
n
cil Bi (q)δ uk ,
l ∈ NLoc ,
(31)
i =1
which is illustrated in Fig. 12. By adding the OBF part, the switching can be smooth and the side effect can be diminished significantly. Unlike in the PWL–OBF model where the switching acts on {uk } and {yk }, here it only manipulates the transience of the latest sample of δ uk and δ yk . Only the related portions of δ uk and δ yk in uk and yk are processed by the adaptive mechanism at a sample point. The rest are not affected. The n-step-ahead predicator based on (31) is formulated as nB k +n
yˆ (k + n) = y(k) +
1 k−1
143
cil Bi (q)δ u(m),
l ∈ NLoc .
(32)
m=k+1 i=1 k−j
δ un ,
(25)
1
which can be further decomposed into a set of simultaneous equations as Eq. (26) in Box I. Compared with the orthodox expression (18), the derived (26) could give us a different insight into the dynamics of a system. Instead of {uk } and {yk }, the input and output channels are now viewed in the definition of {δ uk } and {δ yk }. The entire dynamics of a system can be interoperated as the composition of all the individual DIO impulse responses by per sampling point forward along the time line. The superposition principle and the convolution sum are the unaltered foundations for the evaluation of the dynamic responses. However based on the
Similarly, for the two-input-two-output cases, (31) is extended to y1,k ≈ y1,k−1 +
+
n12
n11
l c11 ,i B11,i (q)δ u1,k
i =1 l c12 ,i
B12,i (q)δ u2,k ,
i=1
y2,k ≈ y2,k−1 +
+
n22 i=1
n21
l ∈ NLoc . l c21 ,i B21,i (q)δ u1,k
i =1 l c22 ,i
B22,i (q)δ u2,k .
(33)
144
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
a0 δ yk . .. a0 δ yi+1 .. . a0 δ y2
+ + +
a1 δ yk−1
+
.. . a1 δ yi .. .
a1 δ y1
··· .. .
+
+
··· .. .
+
···
ai δ yk−i
=
+
.. . ai δ y 1 .. .
=
+
0
=
b1 δ uk−1
+
.. . b1 δ uj .. .
b 1 δ u1
bj δ uk−j
··· .. .
+
+
··· .. .
+
.. . b j δ u1 .. .
+
···
+
0
(26)
Box I.
a0 δ yk a0 yk−1
+ +
a1 δ yk−1 a1 yk−2
+ +
··· ···
+ +
ai δ yk−i ai y(k−1)−i
= =
b1 δ uk−1 b1 uk−1
+ +
··· ···
+ +
bj δ uk−j bj u(k−1)−j
(27)
Box II. Table 1 The unified identification procedure. Step 1 Step 2 Step 3 Step 4 Step 5 Step 6 Step 7
Setup the platform. The measurement and the calculation for the bending angle and the steering angle are based on the auxiliary kinematic setting. Select a number of reference pressures on each chamber. Take identification experiment according to the local partitions. Extract the local data sets. Deoffset and linearise them. Apply the FKcM pole clustering algorithm for the locally identified poles. Construct the OBF set based on the optimally selected poles. Arrange the constructed OBFs in the DIO–PWL–OBF setting.
Table 2 Poles of the local data sets. Local data set
Identified poles
Local data set
Identified poles
locicr 2
0.9288 ± 0.0828i
locdcr 2
0.9689 ± 0.0812i
locicr 3 locicr 4 locicr 5 locicr 6 locicr 7
Fig. 12. Illustration of the DIO–PWL–OBF model structure.
The development of the new DIO–PWL–OBF model structure (31) and (33) has been presented in detail. It is capable of handling the nonlinearities discussed in Section 2.2 and suitable for the identification test designed in Section 4.1. With the use of the auxiliary kinematic setting in Section 3, the unified identification procedure for the class of the PLSAs is concluded and presented in Table 1. 5. Implementation and results
5.1. Platform setup The schematic drawing in Fig. 13 illustrates the implementation setup of Fig. 14. A desktop using Matlab Real-Time Windows TargetTM serves as the central platform to interface the other parts. The system works at a sample rate of 80 Hz, and it is constantly supplied with stable 1 bar compressed air. Every two two-port solenoid valves together with one pressure sensor operate as one unit, which is used to regulate one chamber’s pressure. A test PLSA actuator used in the following identification work is made in the silicon rubber and it is 9 cm in length and 3 cm in diameter. The operating pressure of a single chamber is below 0.4 bar. The constant curvature condition is held satisfyingly for the test PLSA. The deflection angle is accurately approximated to the bending angle. Thus the elongation measurement is neglected.
0.8878 ± 0.0775i
0.9786 ± 0.1072i
locdcr 3
0.8717 ± 0.0154i 0.9883, 0.9125 ± 0.0716i 0.9908, 0.9179 ± 0.0363i 0.9774, 0.9174 ± 0.0793i
locdcr 4
0.8930
locdcr 5
0.9699
locdcr 6
0.9690
locdcr 7
0.9665
Table 3 The pole set {ξ1 . . . ξ10 } selected for OBFs construction.
ξ1,2
ξ3,7
ξ4,8
ξ5,9
ξ6,10
0.9783 ± 0.1068i
0.8770
0.8930
0.9700
0.9910
5.2. Single chamber identification The single chamber identification follows the instruction summarised in Table 1. The chamber 1 of the PLSA is used. As shown in Fig. 15, the local data sets are partitioned based on different levels of pressures. The partition interval is chosen to be 0.05 bar. Fig. 16 demonstrates the experimental data acquired according to the partition setting in Fig. 15. The entire data can be seen as a collection of individual step responses acting at different pressure levels. By deoffsetting the initial condition on both the input and output channels, these step responses are brought to the local sense. The standard Newton–Guess method is employed to take linear identification on the extracted local data sets. The Normalised Root Mean Squire (NRMSE) is used as the measurement of the goodness dcr of the fit. locicr are neglected as their bending angle 1 and loc1 ranges are below 0.5° due to the gas fill-in phase. Despite them, other local data set can be linearised with an over 90% NRMSE fit. Table 2 lists the identified poles for each local data sets. The results justify the different local dynamics observed in Fig. 4. The FKcM clustering is taken on the identified poles. Table 3 lists the calculated pole structure for the OBF construction.
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
145
Fig. 13. Schematic drawing for the platform setup. Table 4 The fitting results for the test in Fig. 18. Results
40-step-ahead 80-step-ahead Simulation
NRMSE
MSE
Full
Icr.
Dcr.
Full
Icr.
Dcr.
97.90% 96.56% 95.28%
96.94% 95.17% 94.34%
97.17% 95.21% 93.09%
1.6585 4.4438 8.3832
2.7656 6.8719 9.4654
2.2352 6.4195 13.3304
Finally, the test in the global sense is carried out. It is a full-range step response. The test and its validation results are demonstrated in Fig. 18. The numerical fitting result is presented in Table 4. Especially in Fig. 18, the switching state is plotted in parallel with the input and output channels. The 14 states are numbered in the icr dcr dcr icr dcr order: locicr 1 , . . . , loc7 , loc7 , . . . , loc1 . Note that loc1 and loc1 are virtually included, but their OBF coefficients are all set to zeros. The transience of the global step response lasts about 1.5 s. Hence the 40- and 80-step-ahead predications which respectively allows a 0.5 and 1 s forerunner are proper settings for validation or other applications. Both the predication and the simulation results show over 90% NRMSE fits. Particularly the goodness measures, NRMSE and Mean Square Error (MSE), for the individual transience parts are given in Table 4. It can be seen that the proposed identification approach is applicable in the single chamber cases. 5.3. Double chamber identification
Fig. 14. A photo of the platform.
Fig. 15. The local partition for the single chamber identification.
The contour plot in Fig. 17 shows the OBF approximation decay rates [23] of each identified poles based on the selected pole structure in Table 3. The overall decay rates are very low. Most of the identified poles are below the decay rate of 0.04, while a few are about 0.9. The results mean that the OBF set using the selected pole structure is sufficient to approximate all the local data sets. Overall, the OBF approximation for the local data sets is over 90% NRMSE.
The double chamber steering and bending identification is also conducted by following the instructions in Table 1. Chambers 1 and 2 of the PLSA are used. As illustrated in Fig. 19, the local data sets are partitioned based on different pressure levels of both the chambers. The pressure range for each chamber is set between 0.20 and 0.35 bar. This is because the bending angle is relatively small when the pressure is below 0.20 bar, and this makes the steering motion difficult to measure by our current magnetic sensor. Fig. 20 shows the experimental data obtained for the identification. It is based on the partition in Fig. 19. Similar to the single chamber situation, each local step response depends on the previous steady state. Each local data set takes the step response driven only by one chamber, while the other chamber’s pressure is kept at its initial value throughout the time range. It is a decoupled approach to identify the SIMO relation between one chamber’s pressure and the resultant two DOF movements. After removing the offsets in both the input and output channels, linear analysis is taken individually. Similar to the single chamber situation, the local sets can be accurately linearised. Each of them has over 90% NRMSE. In the double chamber cases, the overall orders of the identified models become lower. The 1-pole0-zero structure can be properly used for every decoupled SISO dynamical relation. The reason is partly because the stiffness of the
146
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
Fig. 16. The experimental data for the single chamber identification.
Fig. 17. The OBF approximation decay rate of the identified pole based on the selected pole structure in Table 3.
Fig. 18. The global test and validation for the single chamber identification.
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
147
The single chamber actuation is arranged for the validation firstly. Fig. 22 demonstrates Chamber 1’s result. Secondly, a joystick-based pressure control was applied to the platform, by which, both the chambers can be actuated. Its validation result is plotted in Fig. 23. Table 6 provides the numeric results for the figures. The local state indexing is similar to the way used in the single chamber, or one can refer to the local partition in Fig. 19. On average, the NRMSE fit results are above 80%. The system’s dynamical responses are tracked within tolerable errors. It can be seen that the proposed identification approach is applicable in the double chamber cases. 6. Conclusion
Fig. 19. The local partition for the double chamber identification.
Table 5 ij
ij
The pole set {ξ1 , . . . , ξ5 }i=1,2,j=1,2 for the OBFs construction. Decoupled IO pair (ui , yj ) u1 u2 u1 u2
→ y1 → y1 → y2 → y2
i,j
i,j
Pole structure {ξ1 , . . . , ξ5 } 0.981, 0.969, 0.964, 0.957, 0.945 0.986, 0.975, 0.960, 0.950, 0.942 0.998, 0.979, 0.949, 0.935, 0.884 0.983, 0.973, 0.961, 0.946, 0.939
PLSA is higher in the double chamber situation than in the single one. The increased internal pressure acts on stabilising the PLSA’s deformation. The FKcM clustering is taken for the identified poles of each decoupled SISO pair. The calculated poles are listed in Table 5. They are used for the OBF construction. Fig. 21 shows the OBF approximation decay rates of the identified poles in each decoupled SISO system. All the poles are lower than the level of 0.02, which indicates that the selected poles are very suitable for the OBF construction. On average, the OBF approximation for the local data sets is above 85%.
In this paper, a unified identification approach for the general PLSAs has been presented and summarised. And its implementations for the single chamber bending and the double chamber bending and steering are demonstrated separately. The identification approach is based on the newly developed DIO–PWL–OBF models structure which is capable of effectively capturing the nonlinear pressure–shape dynamical relation of the general PLSAs. The DIO–PWL–OBF model structure is consistent with the local identification manner, and it is efficient in implementation. Besides, in the kinematic setting, the new auxiliary configuration subspaces are used. We have shown that by utilising this setting the DOFs of the bending and the steering can be captured in a parameter-free way, and the auxiliary setting can be incorporated with the mainstream kinematic setting. In the next stage, our research is going to focus on using the acquired identification model to facilitate the control design. With the proposed auxiliary kinematic setting for the bending and the steering, a PLSA’s deformation subjected to a payload at the tip end or external contacts by an obstacle can be cast to the kinematic variables, θ b and φ , in a closely relevant way. This can potentially improve the compliance control performance. Currently we are investigating two effects separately: one is that the upper segment acts on the lower, and the other is that the entire arm is forced by an object. They will be investigated and presented in our ongoing work. Acknowledgements Xiaochen Wang would like to convey special thanks to his wife Zhiheng Wu for her constant encouragement to pursue his robotics career.
Fig. 20. The experimental data for the double chamber identification.
148
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
Fig. 21. The OBF approximation decay rate of the identified pole based on the selected pole structure in Table 5.
Fig. 22. The Chamber 1 single actuation test and validation.
Fig. 23. The joystick-controlled double-chamber actuation test and validation.
X. Wang et al. / Robotics and Autonomous Systems 63 (2015) 136–149
149
Table 6 The validation results for the test in Figs. 22 and 23. Cases
Results
NRMSE y1
NRMSE y2
MSE
SIMO Ch. 1
40-step-ahead prediction 80-step-ahead prediction Simulation
95.60% 93.78% 90.72%
87.61% 80.30% 73.61%
1.7701 4.9386 6.4509
MIMO joystick
40-step-ahead prediction 80-step-ahead prediction Simulation
91.83% 88.29% 87.69%
90.98% 87.97% 77.29%
3.6655 5.2224 8.0041
The work described in this paper is partially funded by STIFFFLOP project grant from the European Commission Seventh Framework Programme under agreement 287728. References [1] R. Taylor, D. Stoianovici, Medical robotics in computer-integrated surgery, IEEE Trans. Robot. Autom. 19 (5) (2003) 765–781. [2] P. Dario, B. Hannaford, A. Menciassi, Smart surgical tools and augmenting devices, IEEE Trans. Robot. Autom. 19 (5) (2003) 782–792. [3] A. Menciassi, J. Park, S. Lee, S. Gorini, P. Dario, J. Park, Robotic solutions and mechanisms for a semi-autonomous endoscope, in: IEEE/RSJ International Conference on Intelligent Robots and Systems, Vol. 2, 2002, pp. 1379–1384. [4] W. McMahan, V. Chitrakaran, M. Csencsits, D. Dawson, I. Walker, B. Jones, M. Pritts, D. Dienno, M. Grissom, C. Rahn, Field trials and testing of the octarm continuum manipulator, in: Proceedings 2006 IEEE International Conference on Robotics and Automation, 2006, pp. 2336–2341. [5] K. Suzumori, S. Iikura, H. Tanaka, Applying a flexible microactuator to robotic mechanisms, IEEE Control Syst. 12 (1) (1992) 21–27. [6] G. Chen, M. Pham, T. Redarce, Sensor-based guidance control of a continuum robot for a semi-autonomous colonoscopy, Robot. Auton. Syst. 57 (6–7) (2009) 712–722. [7] B. Jones, I. Walker, Practical kinematics for real-time implementation of continuum robots, IEEE Trans. Robot. 22 (6) (2006) 1087–1099. [8] D. Braganza, D. Dawson, I. Walker, N. Nath, A neural network controller for continuum robots, IEEE Trans. Robot. 23 (6) (2007) 1270–1277. [9] Y. Bailly, Y. Amirat, G. Fried, Modeling and control of a continuum style microrobot for endovascular surgery, IEEE Trans. Robot. 27 (5) (2011) 1024–1030. [10] M. Rolf, J. Steil, Constant curvature continuum kinematics as fast approximate model for the bionic handling assistant, in: 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, 2012, pp. 3440–3446. [11] E. Tatlicioglu, I. Walker, D. Dawson, New dynamic models for planar extensible continuum robot manipulators, in: IEEE/RSJ International Conference on Intelligent Robots and Systems, IROS 2007, 2007, pp. 1485–1490. [12] I. Godage, D. Branson, E. Guglielmino, G. Medrano Cerda, D. Caldwell, Shape function-based kinematics and dynamics for variable length continuum robotic arms, in: 2011 IEEE International Conference on Robotics and Automation, 2011, pp. 452–457. [13] I. Godage, D. Branson, E. Guglielmino, D. Caldwell, Pneumatic muscle actuated continuum arms: Modelling and experimental assessment, in: 2012 IEEE International Conference on Robotics and Automation, 2012, pp. 4980–4985. [14] K. Suzumori, S. Iikura, H. Tanaka, Flexible microactuator for miniature robots, in: Micro Electro Mechanical Systems, MEMS’91, Proceedings. An Investigation of Micro Structures, Sensors, Actuators, Machines and Robots. IEEE, 1991, pp. 204–209. [15] C. Chou, B. Hannaford, Measurement and modeling of Mckibben pneumatic artificial muscles, IEEE Trans. Robot. Autom. 12 (1) (1996) 90–102. [16] B. Tondu, P. Lopez, Modeling and control of Mckibben artificial muscle robot actuators, IEEE Control Syst. 20 (2) (2000) 15–38. [17] S. Billings, S. Fakhouri, Identification of systems containing linear dynamic and static nonlinear elements, Automatica 18 (1) (1982) 15–26. [18] I. Leontaritis, S. Billings, Input–output parametric models for non-linear systems part I: deterministic non-linear systems, Internat. J. Control 41 (2) (1985) 303–328. [19] R. Webster, B. Jones, Design and kinematic modeling of constant curvature continuum robots: a review, Int. J. Robot. Res. 29 (13) (2010) 1661–1681. [20] S. Davis, N. Tsagarakis, J. Canderle, D. Caldwell, Enhanced modelling and performance in braided pneumatic muscle actuators, Int. J. Robot. Res. 22 (3–4) (2003) 213–227. [21] B. Ninness, F. Gustafsson, A unifying construction of orthonormal bases for system identification, IEEE Trans. Automat. Control 42 (4) (1997) 515–521.
[22] P. Heuberger, P. Van den Hof, B. Wahlberg, Modelling and Identification with Rational Orthogonal Basis Functions, Springer, 2005. [23] R. Tóth, P. Heuberger, P. Van den Hof, Asymptotically optimal orthonormal basis functions for LPV system identification, Automatica 45 (6) (2009) 1359–1370. [24] R. Tóth, Modeling and Identification of Linear Parameter-Varying Systems—An Orthonormal Basis Function Approach, Springer, 2010. [25] R. Tóth, J. Willems, P. Heuberger, P. Van den Hof, The behavioral approach to linear parameter-varying systems, IEEE Trans. Automat. Control 56 (11) (2011) 2499–2514.
Xiaochen Wang received the B. Eng. degree from the University of Manchester in 2011. He is a research student at the University of Surrey. The research topic is modelling and control design for a class of soft continuum robots for surgical applications. Tao Geng received the B.Sc. degree from Shandong University, China, in 1990, the M.Sc. degree from Nanjing University of Science and Technology, China, in 1999, and the Ph.D. from Shanghai Jiao Tong University, China, in 2003. His research interests include legged robotics and motor learning. Yahya Elsayed received the M.Eng. degree in Mechanical and Medical Engineering, from the University of Hull, UK and Ph.D. in Materials Engineering, from the University of Surrey, UK. He is a Research associate in the Department of Mechanical Engineering Sciences at the University of Surrey. His research covers modelling and simulation of materials, soft robotics, and biomedical and tissue engineering.
Chakravarthini Saaj is a Senior Lecturer and the Deputy Head of the Surrey Technology for Autonomous systems and Robotics (STAR) Lab at the Surrey Space Centre, University of Surrey. Prior to joining the Surrey Space Centre in 2005, she was a Visiting Research Associate at King’s College London from 2005 to 2007 and a Research Fellow at the Technical University Hamburg-Harburg, Germany from 2002–2003. Her research interests are mainly in the area of developing innovative robotic devices for space and medical applications. Her expertise is in the area of modelling and control of rigid and flexible manipulators, bio-robotics, medical robotics, systems engineering using Model Based Systems Engineering, locomotion of wheeled and legged planetary rovers. She is a reviewer of international journals and has over fifty publications to her credit. Constantina Lekakou received the M.Eng. degree in Chemical Engineering from NTUA, Greece and Ph.D. from D.I.C. Imperial College, London. She is a Reader in the Department of Mechanical Engineering Sciences at the University of Surrey. Her research covers experimental studies and innovation, modelling and simulations of materials, materials processing and systems in the areas of polymer composites and nanocomposites, soft robotics, biomedical materials, energy materials and devices, and electronics manufacturing.