Front. Robot. AIFrontiers in Robotics and AIFront. Robot. AI2296-9144Frontiers Media S.A.55783010.3389/frobt.2021.557830Robotics and AIOriginal ResearchModeling and Reconstruction of State Variables for Low-Level Control of Soft Pneumatic ActuatorsIbrahim et al.Soft Robotics ModelingIbrahimSerhat*KrauseJan ChristophOlbrichAlexanderRaatzAnnikaInstitute of Assembly Technology, Leibniz Universität Hannover, Hannover, Germany
Edited by:Concepción A. Monje, Universidad Carlos III de Madrid, Spain
Reviewed by:Kean C Aw, The University of Auckland, New Zealand
Jorge Muñoz, Universidad Carlos III de Madrid, Spain
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
To further advance closed-loop control for soft robotics, suitable sensor and modeling strategies have to be investigated. Although there are many flexible and soft sensors available, the integration into the actuator and the use in a control loop is still challenging. Therefore, a state-space model for closed-loop low-level control of a fiber-reinforced actuator using pressure and orientation measurement is investigated. To do so, the integration of an inertial measurement unit and geometric modeling of actuator is presented. The piecewise constant curvature approach is used to describe the actuator’s shape and deformation variables. For low-level control, the chamber’s lengths are reconstructed from bending angles with a geometrical model and the identified material characteristics. For parameter identification and model validation, data from a camera tracking system is analyzed. Then, a closed-loop control of pressure and chambers’ length of the actuator is investigated. It will be shown, that the reconstruction model is suitable for estimating the state variables of the actuator. In addition, the use of the inertial measurement unit will demonstrate a cost-effective and compact sensor for soft pneumatic actuators.
soft roboitcssensorization of soft robotsmodeling of soft robotscontrol of soft robotsnonlinear controltest bench design1 Introduction
Soft robots, with flexible shape and infinite configuration possibilities, offer completely new capabilities compared to conventional industrial robots [Trivedi et al. (2008) and Marchese et al. (2014)]. Due to their compliance, soft robots adapt to their environment. This makes them suitable for grippers handling objects with undefined shapes. Since there is no risk of damage in the event of a collision, they are also suitable for human-robot collaboration. A decisive factor determining the movement of soft pneumatic actuators is their design. New actuator designs and mechanisms have been developed for this field of research [Runge and Raatz (2017), Galloway et al. (2013) and Garcia et al. (2020)]. The soft and flexible structures with mostly nonlinear material properties and hyperelasticity present a challenge for modeling, sensing and control. Especially the use of suitable sensors for state detection of the actuator needs to be researched. Due to the low force and high deformability of the actuator, conventional strain gauges cannot be used for this purpose. One option is the use of contact-free camera tracking systems Runge and Raatz (2017). The disadvantage, however, besides the high costs, is the use in confined spaces (high space requirement of the cameras) to avoid covering in cluttered scenes. For use in confined spaces, sensors, which are integrated into the actuator, are more suitable [Szelitzky et al. (2014)].
Table 1 shows different methods for measuring the bending of soft actuators with integrated sensors. Roduit et al. (1998) and Gibbs and Asada (2005) use resistance measurements to determine a bending angle. Roduit et al. (1998) use the difference in position of two parallel cables and Gibbs and Asada (2005) use conductive fibers. Felt et al. (2016) present an inductive measuring method. Wire ties are cast around the fins of a soft pneumatic actuator and its inductance is measured. Here, an inductance change of up to 19% is observed for bends up to 190°. Besides mapping quantified elongation into a geometric shape, another method is the estimation by covering change in pose with an inertial measuring unit. In this case, the pose of an object is observed based on acceleration and rotation rates as well as magnetometer data. Best et al. (2015) and Seel et al. (2014) use this method to measure the bending between rigid links. Seel et al. (2014) achieve an accuracy of 3° at a frequency of 60Hz.
Overview of different measurement methods for the determination of actuator deformation compiled from the current literature.
References
Measurement
Uncertainty
Frequency
DoF
Al Jaber and Althoefer (2018)
Optical
—
—
2
Best et al. (2015)
IMU
—
—
3
Donno et al. (2008)
Optical
0.01∘
1kHz
1
Felt et al. (2016)
Inductive
2∘
—
1
Gerboni et al. (2017)
Conductive
1.08∘
40Hz
1
Gibbs and Asada (2005)
Resistor
2.4∘
2Hz
1
Roduit et al. (1998)
Resistor
2∘
—
2
Seel et al. (2014)
IMU
3.3∘
60Hz
3
Visentin and Fiorini (2018)
Impedance
—
—
2
Yuen et al. (2018)
Capacive
—
10Hz
1
Gerboni et al. (2017) use a commercial flex bend sensor based on conductivity measurements for a soft pneumatic actuator with one degree of freedom (DoF). In the experiment with a closed-loop control, an accuracy of 1.08° is achieved at a clock rate of 40Hz. Yuen et al. (2018) describe the manufacture of strain sensors, which are directly integrated into several film layers in a soft pneumatic actuator. The capacitive based sensor consists of multiple layers with silicone-based conductive electrodes and silicone elastomers as the dielectric. For the measurement using electrical impedance tomography, shredded carbon fibers are arranged as electrodes in the actuator, as described in Visentin and Fiorini (2018). The change in electrical conductivity is measured to reconstruct the bending. The optical bending sensor presented by Donno et al. (2008) is a very accurate measuring method. Non-polarized laser light is polarized by a filter and sent through an optical fiber. If this optical fiber is bent, its polarization changes. Then the change in angle can be recorded via a photo electrode with a second polarizing filter. The accuracy for measurements with up to 1kHz is specified as 0.01°. The use of alloys that are liquid at room temperature should also be mentioned here. EGaIn sensors can also be used for bending measurement [Mengüç et al. (2013)]. Their support fixtures are based on similar or same material as the actuator to avoid inflecting the behaviour of the actuator. However, the production of such sensors is proving to be difficult, and for this purpose separate system components must be developed.
In addition to the sensors, models or neural networks are also used to estimate the state parameters for closed-loop control [Runge and Raatz (2017), Tan et al. (2019) and Katzschmann et al. (2019)] from, for example, pressure measurements. Katzschmann et al. (2019) have published an approach for closed-loop control, where a reduced order finite element model is used for the feedback.
The research presented here aims to enable a low-level control for a three DoF fiber-reinforced actuator (FRA) using orientation measurement of the actuator’s tip. For this purpose, an inertial measurement unit (IMU) is studied. The low-level system description is done at chamber level, where the chamber’s pressure and length are considered. A reconstruction model is developed to observe the state variables, which are relevant for the control. In particular, the observed states include the lengths of the individual actuator chambers, which cannot be measured directly. To build the measurement model, an actuator segment is assumed to have a shape with a piece-wise constant curvature. The parameters are identified using particle swarm optimization and the validation of the measurement model is performed using a camera tracking system. The developed models are used for chamber length control and pressure control. Kinematic relationships between actuators chambers are not modeled. They are included in this concept as unknown disturbances. Compared to Katzschmann et al. (2019), we focus on the closed-loop control of individual segments at chamber level. For this purpose, we consider all components from the effector to the actuator chambers. The lumped second order dynamic model from Skorina et al. (2015) is on low level as well. In contrast to our work the effector system with a pneumatic valve is neglected for modeling.
For the test bench, a PC with Simulink Real-Time as operating system is used. It communicates with the Beckhoff IO-devices over EtherCAT bus. Three Enfield LS-V05 5/3 proportional directional valves are connected to regulate the airflow to the three FRA chambers. To reduce measurement noise, a peripheral EK1100 EtherCAT bus coupler with analog inputs connects five pressure sensors by First Sensors to measure pressure in all chambers, as well as supply and atmospheric pressure. To detect the orientation, the IMU is connected via a microcontroller with an EtherCAT shield. All components are commercially available.
2 MODELING OF THE SYSTEM
In the following, the system components are modeled for use in a closed-loop control (Figure 1). First, the behavior of the valves that regulates the airflow is described. The valve model is needed for the development of the sliding mode control (Section 5.1). Then the connection tubes between valves and actuator chambers are considered. Afterwards the actuator is modeled. For this purpose, the individual dynamic modeling of the chambers are combined to form a complete description of the entire actuator’s geometry. Finally, a state-space representation of the soft robot system is set up.
The closed-loop control uses the difference between observed chamber length and the ones from PCC configuration as feedback.
2.1 Model of Valve
The valve model is based on the work of Ben-Dov and Salcudean (1995) and Richer and Hurmuzlu (2000a), Richer and Hurmuzlu (2000b). For a detailed description of the valve modeling, we refer to our preliminary work in Ibrahim et al. (2019). The air mass flow m˙v through an orifice A of the valve is described withm˙v=cfApuTΨ(pd,pu).
This mass flow depends on the upstream pressure pu and downstream pressure pd as well as the temperature T. Here, the temperature is assumed to be homogeneous throughout the system. The flow coefficient describes the ratio of real and ideal volume flow withcf=V˙realV˙ideal.
The flow function can be calculated withΨ(pu,pd)={κR(2κ+1)κ+1/κ−1pdpu≤pcrit2κR(κ−1)(pdpu)1/κ1−(pdpu)κ−1/κpdpu>pcrit.
The flow Ψ depends on the critical pressure pcrit, which is calculated withpcrit=(2κ+1)κ/κ+1,with κ as the heat capacity ratio. Here, Ψ is constant for a pressure ratio pd/pu, which is smaller than the critical value. Differently, it is a nonlinear function, which depends on the upper- and downstream pressure.
The orifice A of the valve depends on the spool position xs. Assuming that a rectangular slider with an edge length b covers a circular opening with a radius r, the effective area is calculated as a circle segment. With the coordinatexe=xs−b2+r,the area A can be calculated withA(xe)={0xe<0(xe−r)2rxe−xe2+r2arccos(r−xer)0≤xe≤2rπr2xe>2r.
The dynamic of the spool displacement is described with the second order differential equationmx¨s=−Ff+FS−2kxs−dx¨s.Ff stands for the frictional force that occurs during the movement. The spool of the valve has a damping d and a stiffness 2k. Its force is calculated withFS=KSiS=Kτu,with the spool current iS and the motor constant KS, as well as input voltage u, gain K and time constant τ.
2.2 Model of Connecting Tubes
The tubes, which connect the valves with the actuator chambers, affect the air mass flow. The friction in the tube leads to a loss of flow, which causes a time delay, which is based on the sonic speed csonic. The incoming mass flow m˙in can be calculated using the tube diameter d and the tube length l withm˙in(t)=ϕm˙v(t−lcsonic).
Here, ϕ is the attenuation coefficient and it is calculated withϕ=e−RtRT2pendlcsonic.
The pressure pend is measured at the end of the tube. The friction resistance is described in Ibrahim et al. (2019) asRt=0.158d2(BT3/2T+S)1/4(4m˙vdπ)3/4.with the model of Sutherland and its constant B=1.4747×10−6Pas/k and a substance-specific temperature S. By using short and wide tubes, the friction and time delay are minimal and can be neglected. With airtight connectors and tanks, the leakage is also minimal and therefore neglected as well.
2.3 Model of Soft Pneumatic Actuator
In the following actuator modeling is presented using the example of a FRA made of Dragonskin 10 silicon [Polygerinos et al. (2015)]. Figure 2A shows the actuator segment. The Deformation of the actuator is due to expansion of the chambers, which are located along the actuator length. First the entire actuator is considered and a geometric model is created. Then the dynamics of an individual chambers of the actuator are considered and modeled.
Geometry of soft pneumatic actuator with PCC parameters ϕ,θ,l, top view (A), side view (B) and coordinate frames used for the reconstruction model (C).
2.3.1 Geometric Modeling
Based on work from Webster and Jones (2010) the actuator’s shape is approximated with a piece-wise constant curvature (PCC). The configuration is described with the arc length la, the bending direction ϕ and the bending angle θ=(1/r)la with bending radius r. As seen in Figure 2, a segment with length la consists of three symmetrical arranged chambers with a length of li. An elongation of at least one of these chambers leads to a bending and extension of the segment. For a mapping between the PCC parameters (ϕ,θ,la) and the task space coordinates the homogeneous transformation matrixT01PCC=[R01PCCt(1)01]=[cϕ2cθ+sϕ2sϕcϕ(cθ−1)cϕsθcϕ(1−cθ)laθsϕcϕ(cθ−1)sϕ2cθ+cϕ2sϕsθsϕ(1−cθ)laθ−cϕsθ−sϕsθcθsθlaθ0001],from Webster and Jones (2010) is used. This describes the transformation between the segment’s base (CF)0 and the end of the PCC part (CF)1 (Figure 2). For static parts of the segment a linear displacement dE, that results in a transformation T1E with R1E=I and thus only a translation vector t(E)=(0,0,dE,1)T, is added.
The mapping between the PCC parameters and the chamber length is also given in Webster and Jones (2010). For each chamber i the arc length isli=la−dicos(σi−ϕ)θ.
The chamber position is specified by the distance di to the central axis and the angle σi to the x(0)-axis.
2.3.2 Chamber Modeling
The basis for chamber dynamics is the low-level model from Ibrahim et al. (2019). This describes the pressure curve depending on incoming and outgoing mass flow m˙in and m˙out respectively, as well as changes in volume V˙. Considering a chamber with a volume V, the pressure change is described withp˙=RTV(αinm˙in−αoutm˙out)−αpiVV˙.
The chamber’s volume is affected by the difference between internal pressure and atmospheric pressure. Since the FRA only expands the chambers in one axial direction, the spherical approach from Ibrahim et al. (2019) is not suitable. For this reason, a cylinder model is constructed in the following. The volume of a cylinder isV=πrc2l,and it is described by the circular base with the chamber radius rc and the chamber length l. At idle state, the pressure in the actuator chamber is patm and the volume is V=πrc2l0 with the initial length l0.
Analogous to the sphere model, the force due to the pressure difference pdiff=p−patm isFp=pdiffπrc2and the force of the material tension isFσ(l)=σc(ϵ)2πrcwc,with wc as the chamber wall thickness. The elongation depends on the length change ϵ=l−l0/l0 of the actuator chamber. Based on Ibrahim et al. (2019) and the cylindrical shape the chamber’s dynamic is modeled with a second order nonlinear differential equationMl¨+Dl˙=Fp−Fσ(l).
The coefficient M is the chamber’s mass and the coefficient D describes the chamber’s damping.
2.4 State Space Representation
Using the equations above, a state-space representation of the soft robot system is set up. For this purpose, a segment with three chambers is considered. In Figure 1 the control loop is shown. For each chamber a valve is used to regulate the in- and outgoing mass flow. This flow causes a pressure change in the actuators chambers and as a result the chambers in- or deflate. The pressure in each chamber is measured, as well as the actuator’s orientation at a certain point along the arc.
2.4.1 System Dynamics
The system dynamicsx˙=f(x,u)=(p˙1,p˙2,p˙3,l¨1,l¨2,l¨3,l˙1,l˙2,l˙3)T,describes the change of the state variables x=(p,l˙,l)T. These are the chamber pressures p=p(1,p2,p3)T, the chamber lengths l=(l1,l2,l3)T and its derivative l˙=(l˙1,l˙2,l˙3)T. The product of the valve’s opening cross-section Ai and the flow coefficient cf,i is selected as system input u=(ui,u2,u3)T=(cf,1A1,cf,2A2,cf,3A3)T. Introduced in Ibrahim et al. (2019), the use of fast switching valves allows to neglect the spool dynamic Eq. 7.
The states are described, based on Eq. 14, asp˙i=RTaπrc2li(αin,im˙in,i(ui,pi)−αout,im˙out,i(ui,pi))−αipilil˙i,and based on Eq. 18 asl¨i=rc,i2Mi((pi−patm)π−σc,i(ϵ)2π−Dirc,i2l˙i),with i=1,2,3. For mass flow equation please refer to (Section Model of Valve).
2.4.2 Measurement Model
The measurement model is used for mapping between the state space x and the measurement output y=(p1,p2,p3,γx,γy)T with bending angles γx and γy. The pressures p=(p1,p2,p3)T are both state and measurement quantities. The orientation at any point of the arc is represented with Euler angles in RPY notation. The corresponding rotation is described with the rotation matrixRRPY=[cγzcγycγzsγysγx−sγzcγxcγzsγycγx+sγzsγxsγzcγysγzsγysγx+cγzcγxsγzsγycγx−cγzsγx−sγycγysγxcγycγx].
From the states l=(l1,l2,l3)T, the PCC parameters (ϕ,θ,la) are calculated first and then the bending angles γx and γy at arc position lm are determined by comparison of the entries of the rotation matrices Eqs. 12, 22.
If the quotientl2−l1l3−l2=d1cos(σ1−ϕ)−d2cos(σ2−ϕ)d2cos(σ2−ϕ)−d3cos(σ3−ϕ),is formed from Eq. 13 and the addition theorem cos(α−β)=cosαcosβ+sinαsinβ is applied, the equationtanϕ=l2−l1l3−l2(d2cosσ2−d3cosσ3)−(d1cosσ1−d2cosσ2)(d1sinσ1−d2sinσ2)−l2−l1l3−l2(d2sinσ2−d3sinσ3),results. Withl2−l1=θ(d2cos(σ2−ϕ)−d1cos(σ1−ϕ)).the angleθ=l2−l1d2cos(σ2−ϕ)−d1cos(σ1−ϕ),is determined. Transposing Eq. 13, the arc length isla=l1+d1cos(σ1−ϕ)θ.
Getting from PCC parameters to RPY angles, we first determine the rotation matrix at the measuring position. Based on the PCC parameters and the measurement position lm we determine the rotation matrix R0m(lm). For a measurement position at any point l=lm on the central arc, the transformation matrix T0m(lm) is based on T01PCC from Eq. 12 with bending angleθm=θlmla,and arc length lm.
From the comparison of the rotation matrices Eqs. 12, 22 followsγy=arctan2(−r31,r112+r212).
If cos(γy)=0, the angle γx=0. For other casesγx=arctan2(r32cosγy,r33cosγy),applies. From these equations the measurement model can be set up withy=g(x)=[pγxγy].
The state space representation consists of x˙=f(x,u) and y=g(x) is used to analyze and simulate the system. In practical applications not only the description of the system behavior is relevant. Furthermore, a consideration of the states x during operation is essential.
In summary, the measurement model is based on the correspondence of the rotation matrices, which was established on the one hand by the sensor values in RPY coordinates and on the other hand by the approximation of the actuator shape by the PCC parameters. The offset between the actuator’s tip and the measurement position is also included by shifting the position with the PCC parameters.
3 Reconstruction
In the control loop shown in Figure 1, the controlled variable is the length of the chambers. Since the lengths are not directly measurable, a reconstruction is necessary. In the following a model for reconstruction is described, which determines the state of the actuator from sensor measurements.
3.1 Reconstruction With a Static Inverse Measurement Model
The pressures p and the orientation at a certain point on the arc are available as measured system outputs. The relationship between the measured and state variables is determined by the measurement model from Eq. 31. The pressure is directly mapped from state to system output. The inverse function of the measurement model is not sufficient to determine the actuators shape, because the mapping is not bijective. With the measured orientation, only a relative chamber length is captured. Therefore the overall actuator length is unknown and a reconstruction of the length is performed. Mapping the system output to the state variables, the inverse measurement model[pl1l2l3]=g−1(pγxγy),is formed with the measurement equations that have already been established as well as the system dynamics. In a first step the orientations measurement is used to calculate the PCC parameters ϕm and θm at a measurement position lm. The arc length la cannot be determined, because it does not affect the orientation as seen in Eq. 12.
Unlike the two measurement variables γx and γy the rotation γz around the z-axis is unknown. Checking matrix Eq. 12, it becomes apparent that the entries r12 and r21 are identical. To match the rotation matrices, this must also apply to R10RPY on Eq. 22. Thus followsr12,RPY=r21,RPY,sinγzcosγy=cosγzsinγysinγx−sinγzcosγx,and therefore for the angleγz=arctan2(sinγxsinγy,cosγy+cosγx).
If the rotation matrices Eqs. 12, 22 are compared with each other, the bending direction can be found intanϕm=sinϕmcosϕm=r3,2r3,1.
This Results in the Following Anglesϕm=arctan2(cγysγx,−sγy)+π,
To get the bending angle θ we needtanθm=sinθmcosθm=r1,3cosϕ+r2,3sinϕr3,3,and so it isθm=|arctan2[(cγzsγycγx+sγzsγx)cϕ+(sγzsγycγx−cγzsγx)sϕ,cγycγx]|.at the measurement position lm. Since the arctan definition range is [−π,π], Eq. 37 is shifted by π, so ϕm is in the [0,2π] PCC definition range. The angle θm is positively defined, hence no full consideration of all quadrants of the inverse angle function arctan is necessary for the PCC parameter.
In contrast to ϕm and θm, the arc length la cannot be reconstructed from the orientation measurement. A reconstruction based on the actuator’s model is necessary. Considering a static case the force equilibrium isFp=Fσ.
With Eqs. 16, 17, strain, based on pressure, isϵi=σc−1(pdiff,irc2wc).
Through the defined strain the chamber lengthsli=(1+ϵi)l0,i,can be determined. With Eqs. 24, 26 the arc length la is known from Eq. 27.
From the reconstructed arc length la the bending angle at segment endθ=θmlalmcan be derived with Eq. 28. In a last step the chamber lengths li are determined with Eq. 13.
Reconstruction of the state variables was also performed by comparing the rotation matrices. With the orientation measurement, the PCC parameters can be determined at the measurement position. Since the actuator length cannot be found using the orientation measurement, it was necessary to look at the actuator forces. For static case the state variables can be reconstructed now.
3.2 Measurement Devices for Shape Sensing
In this research, a camera tracking system and an IMU are used to capture actuator’s shape. The camera tracking system is used for identification and validation experiments and the IMU is used for orientation measurement of the actuator’s tip (Figure 2A).
3.2.1 Camera Tracking System
An OptiTrack Flex three camera system is installed to track the FRA’s segment tip. The system is infrared based, therefore reflecting markers are attached to the end of the FRA. With a resolution of 100 frames per second, 2D images of six cameras are reconstructed into a 3D representation, thus calculating the tip’s position. As a result, it is possible to record the position of the actuator with the cameras in a cycle of 100 Hz.
3.2.2 Inertial Measurement Unit
The IMU Waveshare12476 has an ICM20948 chip, which includes a compass, a gyroscope and an accelerometer. The rotations γxIMU, γyIMU and γzIMU in the frame of the IMU (CF)IMU can be estimated. The IMU uses the earth’s gravitational force (direction of the z-axis) and the earth’s magnetic field (direction of the y-axis) for the orientation of the basic coordinate system. Furthermore, the IMU includes a processor for motion processing algorithms, which forwards the data via the I2C bus to the host processor. In this setup, a clock rate of 40 Hz is achieved.
4 Identification
In the previous sections model equations, which depend on various parameters, have been derived. Therefore the parameters have to be determined. Some parameters are based on literature, others can be found in CAD models or can be measured directly. However, a few parameters cannot be determined directly and thus they must be identified. In the following, parameters to be determined are highlighted and their identification procedures are described.
4.1 Parameter of Valve Model
The function of the valves is described with the mass flow Eq. 1. The following parameters have to be defined:
• The ideal gas constant R and the isentropic exponent κ,
• discharge coefficient cf,
• as well as the mapping between the valve’s orifice A(u) and the input voltage u.
A detailed description of parameter choice and identification can be found in Krause et al. (2019).
4.2 Parameter of Actuator Model
Regarding the actuator, a distinction is made between chamber modeling and geometric modeling of the entire actuator. First, the chamber dynamics is considered. For Eqs. 14, 18, the parameters needed are
• the coefficients αin and αout based on the occurring heat transfer,
• the stress-strain curve σ(ϵ),
• the chamber radius rc and wall thickness wc and the
• chamber’s mass M and damping D.
The identification process of these parameters is also mentioned in Krause et al. (2019). In addition to the procedure mentioned above, an identification of the actuator geometry is carried out. Also, a more precise volume description for identifying the stress strain curve is possible. The actuator’s geometric model is parameterized with
• the chamber positions, that consist of the angle σi and the offset di to the central axis,
• the offset from the end of the PCC segment to the end effector dE,
• as well as the length of the unstressed PCC segment l0.
The chamber position is based on the design of the actuator’s mold. If the three chambers are arranged as in Figure 2, their position is specified with the angleσi=2i+13π.
Assuming a symmetric design, the offsets are equal with di=dc. The chamber displacement dc and the linear distance to the end effector dE are based on the actuator’s CAD data. The initial actuator length l0 of the PCC segment needs to be identified. For this purpose, the pressure control from Ibrahim et al. (2019) is used to deflect the actuator in different bending directions and angles. The true position is recorded with the camera tracking system described in Section Camera tracking system. A marker is attached at the end effector with a displacement dmk. The camera tracking system records the marker position r(cam)mk in the camera frame (CF)cam. This is calibrated with a ground plane to match its frame orientation Rcam0 and origin to the actuators base (CF)0. Only a displacement dtop at the top of the actuator mount is left. Thus, the transformation iscamT0=[camR0(cam)rtop01],with Rcam0=I and t(cam)top=[0,0,dtop,1]T. With the homogeneous transformation matrix T01 from Webster and Jones (2010), the marker position r(1)mk=[0,0,dmk,1]T is mapped to the base coordinate system (CF)0. Hence, the estimated marker’s position isr^(cam)mk=Tcam0T01PCCr(1)mk.
To estimate the marker position, the transformation matrix T01PCC and therefore the PCC parameters are needed.
For measurement, an IMU is used. As described in Section Reconstruction with a static inverse measurement model, an estimation of the actuator elongation is necessary. A special actuator design leads to constraints for the arc length la. If there is a construction with high stiffness in the longitudinal side, one can set la=l0=constant. This is the case for the 3D printed PneuNet actuator from Garcia et al. (2020). While the chambers lengthen and shorten during the bending of the actuators, the central axis does not change in length.
If only positive elongation of the chambers is possible, the arc length la is not fixed and approximated as a function of the bending angle θ. Assuming a linear relationship, we estimate the change of arc length to bela=l0+dnθ.
This is equivalent to a displacement of a neutral axis in bending direction, where there is no strain. This behavior is typical for the fiber-reinforced actuator with at least one relaxed chamber. For identification with length approximation and camera tracking system the.
• displacement dtop of the camera base frame,
• the displacement of the neutral axis dn
• and the marker offset dmk
are also needed.
4.2.1 Orientation of the Inertial Measurement Unit
The orientation of the IMU is recorded at lm=la, the tip of the segment. For IMU measurement the transformation RIMU1 is unknown. The rotation matrixRIMU1=R(ϕz,ϕy,ϕx),is built with RPY angles (ϕz,ϕy,ϕx), which must also be identified. To determine the PCC parameter, the rotation matricesR01PCC=R0IMURIMU1must be equal. For the identification routine, we first determine the yaw angle θz. This is accomplished similar to Eq. 35. By inspecting the PCC rotation matrix Eq. 12, notably the entries r1,2 and r21 are identical. From this relation and from Eq. 49 the yaw angle is determined withtanθz=(sϕycθx−cϕysϕzsθx)sθy−cϕycϕzcθy−sϕxcϕysθx+(cϕxcϕz+sϕxsϕysϕz)cθx(sϕxcϕycθx+(sϕxsϕysϕz+cϕxcϕz)sθx)sθy+(sϕxsϕycϕz−cϕxsϕz)cθy+sϕysθx+cϕysϕzcθx.
With Eqs. 36, 38 the bending direction ϕ and the bending angle θ can be derived from the rotation matrix of the IMU. The arc length la is approximated with Eq. 47.
With these PCC parameters the homogeneous transformation matrix T01PCC is built and the marker position r^(cam)mk can be estimated from IMU measurements with Eq. 46.
4.2.2 Optimization
With particle swarm optimization, the parameter values are optimized to fit the calculated positions r^(cam)mk from the sensor to the true data from the camera tracking system. The cost function is built with the Euclidean distance of the marker positions. In order to consider the measurements in the deflected state more intensely, the costs are increased with the Euclidean distance in the (xy)-plane. Consequently, the cost function for N measurements isc=∑i=1N||r^(cam)mk,i−r(cam)mk,i||(x^(cam)mk,i−x(cam)mk,i)2+(y^(cam)mk,i−y(cam)mk,i)2N.
The identification measurement is recorded with pressure steps in each chamber separately and in pairs of two. This movement covers many operation points. After the oscillation has subsided, the measurement data of each stage i is recorded and averaged for noise reduction. This provides the identification data set. The identification of the IMU sensor results in a mean error ||r^(cam)mk−r(cam)mk|| of 1.6mm and a standard deviation of 0.9mm. For validation, sine pressure curves with different phase shifts are recorded. Here, again a path with different operation points is selected. The path in Cartesian x, y and z-direction is shown in Figure 3. Furthermore, the deviations of the individual coordinates between IMU and the camera tracking system are shown in Figure 3. The validation results in a mean error of 4.1mm with a standard deviation of 0.9mm. The largest deviations occur at changes of the moving direction.
Validation for IMU identification.
For length reconstruction, based on strain from Eq. 41, the relation σ(ϵ) is needed. At steady state, the pressure pi and the chambers’ lengths li are recorded. If only one chamber is actuated, there is a bending dependent extension of the arc length. First, the PCC parameters from Eqs. 36, 38 and the length la=l0+dnθ as well as the real bending angle Eq. 43 are determined. With this configuration the chambers’ lengths can be calculated by Eq. 13. To prevent falsification due to wrong identification of chamber radius rc and wall thickness wc, the augmented stiffnessS(ϵ)=σ(ϵ)wcrc=pdiff2,is identified. With different steady states, a look-up table for the stress is filled. The results for all three chambers are shown in Figure 4. The values of the individual chambers differ due to manufacturing tolerances. It should be noted here that the elongation of an individual chambers refers to the length l0 of the central axis of the entire PCC segment. Therefore an elongation ϵ≠0 is possible although the material is not under tension.
Identified stiffness parameters of the actuator chambers.
5 Closed-Loop Control Using the Reconstruction Model
In this section the previously described reconstruction model is used for low-level closed-loop control of pressure and chambers’ length of the FRA. For the pressure control, a sliding mode control (SMC) is used. In our previous research [Ibrahim et al. (2019)], it was shown, that a SMC was worse than a PI controller, due to the lack of information about the volume of the chamber of the actuator. With the information of the chamber length from the reconstruction model in this research and the known radius of the chambers, their volume can be calculated. This is used to design the SMC and the results are compared with a PI controller.
In addition, a closed-loop control for the chambers’ length using a PID controller is implemented and evaluated. Here, a path is also traced and the PCC and Cartesian coordinates are considered. In Figure 1, the layout of the control system is shown with w as reference input and y as feedback.
5.1 Closed-Loop Control of the Pressure With a Sliding Mode Control
The control law for sliding mode control isu=ueq−ξsat(sζ),with its parameters ξ, the maximum gain, and ζ, which depends on a feasible error p˜max and control frequency fc byζ=2πfc5p˜max.
The tracking error is p˜=pa(t)−pd(t) and leads tos(p˜)=(ddt+λ)n−1p˜(t).
Its first order n=1 becomess(p˜)=p˜=pa−pdand its time derivative with Eq. 14 iss˙(p˜)=RTV(αinm˙in−αoutm˙out)−αpaV˙V−p˙d.
The condition for equivalent controls˙(p˜)=0,is converted to u. If pd≥pa the air has to flow into the actuator. It results inm˙out=0,andm˙in=VαinRTa(p˙d+αpaV˙V).
Inserting Eq. 1 in Eq. 60 gives the equivalent controlueq=cfA=VRTa(p˙d+αpaV˙V)psΨ(ps,pa)ϕ(u(t−T),ps,pa)αin.
For calculating the attenuation coefficient from Eq. 10 with tube resistance Eq. 11, the previous mass flow and therefore the previous input u(t−τ) is used. The calculation of the equivalent control for pd<pa is determined analogously.
The sliding mode controller was compared with a PI controller (Figure 5). Here, two different operation (1.8×10−5Pa and 2.5×10−5Pa) points were approached in one jump and one stair function. The evaluation of the control quality for the steps is shown in Table 2. Here, the overshoot, the rising time, the settling time (5%) and the control deviation are considered. The SMC has a lower overshoot at all steps compared to the PID controller. The greater the height of the step, the greater the difference between SMC and PID overshoot (comparison 5s and 15s). Since the PID controller is set dynamically, the rising time is shorter than the time of the SMC. The SMC performs better than the PID controller in terms of settling time. The control deviation shows a weakness of the SMC. While with rising steps (5s, 15s, 25s and 30s) the control deviation between SMC and PID is comparable, SMC shows a clear deviation for falling steps (10s, 20s and 35s). This problem can be solved by optimizing the controller parameters of the SMC [Ibrahim et al. (2019)]. It can be seen that the SMC in combination with the reconstruction model and the IMU provides a better performance than a PID controller for pressure control. Especially with different operation points, the advantages of the SMC become clear.
Comparison between SMC and PID Controller for pressure control.
Performance of SMC and PID controller for pressure control.
Time [s]
Controller
Overshoot [Pa]
Rising time [ms]
Settling time [ms]
Control deviation [Pa]
5
SMC
5,218
92
275
117
PID
8,943
78
698
124
10
SMC
263
191
4,778
1939
PID
7,838
92
672
241
15
SMC
3,315
182
275
427
PID
10,511
78
723
476
20
SMC
368
250
349
1884
PID
11,920
80
545
219
25
SMC
2,714
90
164
117
PID
8,185
77
713
131
30
SMC
2,832
136
497
398
PID
7,187
102
598
475
35
SMC
5,449
114
316
123
PID
10,078
68
546
129
5.2 Closed-Loop Control of the Chambers’ Lengths
Beside the pressure control, a closed-loop control with the previously described state variables li (chambers’ length) as feedback is considered. The reference variables are the bending direction ϕ(t), the bending angle θ(t) and the segment length la. With Eq. 13 each chamber length is calculated and is used as control variable. As a controller, a PID controller designed with Ziegler-Nichols’ method is used [Ziegler and Nichols (1942)]. During controller design, it was found, that the chambers behave differently, which can be attributed to manufacturing tolerances. Thus, separate controllers are designed for each chamber of the actuator. The step response of the three chambers for two different operation points (0.144m and 0.151m) is shown in Figure 6. Furthermore, the deviation of the individual chamber lengths can be taken from Figure 6. It is shown, that apart from the steps, the measured chambers lengths follow the desired chambers lengths. The evaluation of controller performance for the first two steps is shown in Table 3. For the control deviation, a good value is reached with <0.05×10−3m for all chambers. Since the controller is set dynamically, the overshoot is large (<<6×10−3m) but the rising time is small. Within the chambers, chambers one and three show stronger overshoots than chamber 2 with similar rising time. These differences can be explained by the manufacturing tolerances. The overshoots from Table 3 and Figure 6F do not match, because the Δ of the chambers lengths is shown in the Figure. At the moment of the step there is a dead time, so the delta is greater than the overshoot.
Step response of PID controllers of chambers’ lengths.
Performance of the PID controller for closed-loop control of the chambers’ lengths.
Time [s]
Chamber
Overshoot [m]
Rising time [ms]
Settling time [ms]
Control deviation [m]
20
1
3.88×10−3
136
1870
0.02×10−3
2
2.22×10−3
132
2,613
0.01×10−3
3
3.71×10−3
124
1,431
0.04×10−3
40
1
−0.81×10−3
193
761
0.03×10−3
2
−0.43×10−3
173
940
0.02×10−3
3
−0.63×10−3
149
1854
0.03×10−3
60
1
5.99×10−3
131
2,605
0.04×10−3
2
3.28×10−3
160
1,442
0.02×10−3
3
5.71×10−3
142
1,444
0.01×10−3
80
1
−0.84×10−3
164
930
0.02×10−3
2
−0.51×10−3
218
1,098
0.03×10−3
3
−0.71×10−3
170
2,436
0.05×10−3
For testing the controller performance, different PCC parameters are specified. With bending angles θ=15∘ and θ=25∘ and actuator length of la=0.15m, multiple bending directions ϕ=k30∘ with k=0,1,…,11 are used to compute the reference in Eq. 13. Figure 7 shows the lengths of the chambers during the movement. The overshoots are clearly visible in the steps, whereby these increase with increasing step height. This is clearly shown in the Δ of the chamber lengths in Figure 7. First l1 has the largest overshoots, then l2 and finally l3. This is due to the dynamic setting of the PID controller. It also becomes, clear that the desired value is not achieved with small chamber lengths. One reason for this could be the stretching of the chamber during the previous actuation. Since no negative pressure is generated, the desired length cannot be achieved. Figure 8 shows the movement in the PCC parameters ϕ and θ of the actuator. The initial position of the actuator is not defined for the PCC parameters (singularity). For a better view the measurement is set to ϕ=0. Also the reference of ϕ=0∘ lead to results in a neighborhood of ϕ=360∘. Thus ϕ oscillates at the beginning of the experiment in Figure 8. The steps of the desired angle ϕ are well achieved. At the angle θ the larger steps are not quite reached.
Desired and measured values of the chambers’ lengths during validation of PID Controller.
Desired and measured PCC parameters ϕ and θ during validation of PID Controller.
The controller performance is validated with the camera tracking system. For this purpose, the desired marker position is determined based on work from Section Reconstruction with a static inverse measurement model and Parameter of Actuator Model. As shown in Figure 9 there is a mean deviation of 4.3×10−3m between the desired path and the reconstructed position. The reconstruction differs from the validation data from the camera tracking system with a mean of 1.9×10−3m and a standard deviation of 2.2×10−3m. The overshoots in x and y are similar in size and the overshoots in z are smaller by a factor of 3. The reason for this is that the influence of the chamber length on x and y is greater than on z.
Desired and measured Cartesian parameters during validation of PID Controller.
In addition to a set point stabilization that is done with the steps in the validation above, a control for path tracking is considered, too. For that a circle with a radius of about 31×10−3m is constructed with the PCC parameters θ(t)=25∘, ϕ(t)=360t/T∘ and la=0.1402m is given as reference path. The results for times of circulation T=200s, T=100s and T=5s are shown in Figures 10–12. It can be seen that there is still a maximum deviation between 4×10−3m and 14×10−3m. The error increases with decreasing of path time for the circle path. The reason for this is the feedback frequency of the IMU. This shows, that a reconstruction model an IMU can be used for suitable low level closed-loop control of a soft pneumatic actuator.
Movement in a circular path in T=200s.
Movement in a circular path in T=100s.
Movement in a circular path in T=5s.
6 Conclusion
In this research, a model for reconstruction of state variables of a soft pneumatic actuator with an inertial measurement unit was demonstrated. A fiber-reinforced soft pneumatic actuator was chosen for the investigation. With the PCC approach, the shape and the deformation variables of the actuator were described and a geometrical model was developed. Then the dynamics of the actuator chambers were modeled using a nonlinear second order differential equation. A state space representation of the soft robotic system was set up with the air pressure, the chambers’ length and the first and second time derivation of this as state variables. A measurement model was set up to map between the state variables and the measurement data of the IMU. With the geometric model and data of the pressure and orientation measurement, a reconstruction model for the deformation angles was set up, concerning the specific material properties of the actuator. The reconstruction model was used to determine the volume for a sliding mode controller of pressure. Furthermore, the control of the chambers’ lengths of the actuator was investigated.
In the validation of the reconstruction model, a mean error of 4.1mm with a standard deviation of 0.9mm results from the camera data for a sinusoidal signal. Also, no large deviations between the reconstruction model and the camera data were detected, during the test of the controller. These have a mean error of 2.2mm with a standard deviation of 1.9mm. When designing the PID controller using the reconstruction model for the closed-loop control of the chambers’ length, a good control quality were evaluated with settling time <2605ms and an control deviation <0.05×10−3m. The controller was set dynamically so that overshoots were present in the step response. For the pressure control, a SMC using the information of the chambers’ lengths was designed. The evaluation shows a better performance of the SMC compared to the PI controller, especially with different operation points.
To increase the performance of the controller, it is necessary to increase the feedback frequencies of the IMU. A filtering of the measurement signals can also be considered. Due to the fact that the pressure dynamic differs from the actuator dynamic, the reconstruction of the chambers’ length with pressure measurement is insufficient. Therefore, an observer with known model dynamic is necessary. In further work, Kalman-filtering approach for state estimation is recommended. In this approach, different measurement rates and noises from sensors are concerned.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.
Author Contributions
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.
Funding
Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)-405032969.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
ReferencesAl JaberF.AlthoeferK. (2018). “Towards creating a flexible shape senor for soft robots,” in IEEE International Conference on Soft Robotics (Robosoft), Livorno, Italy, 114–119. April, 2018. 10.1109/ROBOSOFT.2018.8404906Ben-DovD.SalcudeanS. E. (1995). A force-controlled pneumatic actuator. BestC. M.WilsonJ. P.KillpackM. D. (2015). “Control of a pneumatically actuated, fully inflatable, fabric-based, humanoid robot,” in IEEE-RAS 15th international conference on humanoid robots (humanoids). New Jersey, NJ: IEEE, 1133–1140. 10.1109/HUMANOIDS.2015.7363495DonnoM.PalangeE.Di NicolaF.BucciG.CiancettaF. (2008). A new flexible optical fiber goniometer for dynamic angular measurements: application to human joint movement monitoring. FeltW.SuenM.RemyC. D. (2016). “Sensing the motion of bellows through changes in mutual inductance,” in International Conference on Intelligent Robots and Systems IROS 2016. New Jersey, NJ: IEEE, 5252–5257. 10.1109/IROS.2016.7759772GallowayK. C.PolygerinosP.WalshC. J.WoodR. J. (2013). “Mechanically programmable bend radius for fiber-reinforced soft actuators,” in 16th international conference on advanced robotics (ICAR), November 2013, Montevideo, Uruguay. 10.1109/ICAR.2013.6766586GarciaM., D. S.IbrahimS.CaoB.-H.RaatzA. (2020). Cluj-Napoca, Romania EuCoMeS 2020: New Trends in Mechanism and Machine Science.488–495. 10.1007/978-3-030-55061-5_55GerboniG.DiodatoA.CiutiG.CianchettiM.MenciassiA. (2017). Feedback control of soft robot actuators via commercial flex bend sensors. GibbsP. T.AsadaH. H. (2005). Wearable conductive fiber sensors for multi-Axis human joint angle measurements. IbrahimS.KrauseJ. C.RaatzA. (2019). Linear and nonlinear low level control of a soft pneumatic actuator. In IEEE International Conference on Soft Robotics (RoboSoft), February 2019, 434–440. 10.1109/ROBOSOFT.2019.8722737KatzschmannR. K.SantinaC. D.ToshimitsuY.BicchiA.RusD. (2019). “Dynamic motion control of multi-segment soft robots using piecewise constant curvature matched with an augmented rigid body model.” in IEEE International Conference on Soft Robotics (RoboSoft), February 2019, Seoul, South Korea, 454–461. 10.1109/ROBOSOFT.2019.8722799KrauseJ. C.IbrahimS.RaatzA. (2019). “Evaluation environment for control design of soft pneumatic actuators,” in MarcheseA. D.KomorowskiK.OnalC. D.RusD. (2014). “Design and control of a soft and continuously deformable 2D robotic manipulation system,” in IEEE international Conference on Robotics and automation (ICRA). 2189–2196, September, 2014, Hong Kong, China. 10.1109/ICRA.2014.6907161MengüçY.ParkY.-L.Martinez-VillalpandoE.AubinP.ZisookM.StirlingL. (2013). “Soft wearable motion sensing suit for lower limb biomechanics measurements.” in IEEE Int.ernational Conf.erence on Robotics and Automation, May 2013. 5309–5316. 10.1109/ICRA.2013.6631337PolygerinosP.WangZ.OverveldeJ. T. B.GallowayK. C.WoodR. J.BertoldiK. (2015). Modeling of soft fiber-reinforced bending actuators. RicherE.HurmuzluY. (2000b). A high performance pneumatic force actuator system: Part II-nonlinear controller design. RicherE.HurmuzluY. (2000a). A high performance pneumatic force actuator system: Part i—nonlinear mathematical model. RoduitR.BesseP.-A.MicallefJ.-P. (1998). Flexible angular sensor and biomechanical application. RungeG.RaatzA. (2017). A framework for the automated design and modelling of soft robotic systems. SeelT.RaischJ.SchauerT. (2014). IMU-based joint angle measurement for gait analysis. SkorinaE. H.LuoM.OzelS.ChenF.TaoW.OnalC. D. (2015). “Feedforward augmented sliding mode motion control of antagonistic soft pneumatic actuators,” in IEEE international conference on robotics and automation (ICRA), May 2015, Seattle, United States, 2544. 10.1109/ICRA.2015.7139540SzelitzkyE.KuklyteJ.MandruD.O’ConnorN. (2014). Low cost angular displacement sensors for biomechanical applications - a review. TanC. P. E.NurzamanS.LooJ. (2019). “Non-linear system identification and state estimation in a pneumatic based soft continuum robot,” in 2019 3rd IEEE Conference on control technology (CCTA), August 2019. Hong Kong, China, 10.1109/CCTA.2019.8920693TrivediD.RahnC. D.KierW. M.WalkerI. D. (2008). Soft robotics: biological inspiration, state of the art, and future research. VisentinF.FioriniP. (2018). “A flexible sensor for soft-bodied robots based on electrical impedance tomography,“in IEEE international conference on soft robotics (RoboSoft), Livorno, Italy.158–163. 10.1109/ROBOSOFT.2018.8404913WebsterR. J.JonesB. A. (2010). Design and kinematic modeling of constant curvature continuum robots: a review. YuenM. C.Kramer-BottiglioR.PaikJ. (2018). “Strain sensor-embedded soft pneumatic actuators for extension and bending feedback,” in IEEE international conference on soft robotics (RoboSoft), April 2018. 10.1109/ROBOSOFT.2018.8404920ZieglerJ. G.NicholsN. B. (1942). Optimum settings for automatic controllers. Glossarycϕ