A practicable approach for stability analysis in soft robotics: A workspace analysis for a real soft robot

Soft robotics is a novel field that has gained importance in recent years. The idea in soft robotics is to use materials with a softness comparable to human tissue. With these materials, soft robots are capable to adapt to their environment without extensive use of a sensor system. This makes them interesting, for example, in rehabilitation or as a surgery tool. Due to their softness, soft robots are typically underactuated, which means that they have a lower number of actuators than degrees of freedom. This softness leads to significant deformation of the initial soft robot configuration when external loads are applied, for example, when carrying a weight. Hence, for practical applications, it is important to determine the workspace of a soft robot especially when the robot carrying an object. For this setup, the deformation of a pneumatically actuated self‐developed soft robot, which carries a weight, was simulated. The workspace was analyzed for different air pressures in the chambers of the soft robot as well as for different applied weights. Analytical equations for the description of the workspace were developed from the results of the simulations. From these equations, the effects of air pressure and weight on the workspace can be efficiently determined. In addition, the influences on different robot configurations (e.g., robots with different lengths) can be directly compared with each other.


INTRODUCTION
Soft robots, also called continuum robots, are inspired by biological systems and stand out from conventional robots because of their intrinsic properties [1,2].These properties lie in the use and arrangement of soft materials.Soft in this context means that the stiffness of the materials used is in the range of biological materials, with a Young's modulus between 10 4 and 10 9 Pa.An almost unrestricted deformation of their flexible structure gives the soft robots their unique mobility [3,4].The areas of application for soft robots range from medical applications to safe human-machine interaction [3,5].Due to the soft material, impacts can be largely absorbed.Conventional robots, however, require complex stability and protection algorithms for their use on and with humans.These robots are therefore used in predefined environments.Soft robots can dispense with these methods thanks to their compliance.They are therefore also designed for use in unstructured F I G U R E 1 Structure of the modular soft robot at the Institute for Computational Physics in Engineering [6].
environments [5].For the practical use of soft robots, it is essential to know their workspace.The workspace describes the area in which the robot can act due to external and internal influences.External influences are, for example, the weight of an object to be carried or contact forces between the robot and an obstacle.Internal influences are, for example, the air pressure in the chambers of a pneumatically actuated soft robot.
Due to the complex and continuous deformation, based on the materials used, the description of the motion of soft robots can quickly become very elaborate and complicated.Simulations are performed to determine the deformation.With our approach, we want to use the results of simulations to implement a simplified description of the workspace in terms of analytical equations.The aim is to determine the influence and dependence of air pressure for actuation and external loads due to carrying an object on the workspace.Using the analytical equations, an efficient control of the air pressure can be realized.For a given external load and a preset air pressure, the point to be reached in space can be calculated directly.

THEORETICAL FOUNDATIONS
The section deals with the soft robot under investigation and its mathematical description.First, the invented soft robot, for which the workspace is analyzed, is briefly described.Then, the Cosserat theory on which the simulation is based is presented.The third part explains how the soft robot was implemented in the simulation environment PyElastica.

Soft robot under investigation
For the workspace analysis, the soft robot developed at our institute, as shown in Figure 1, was taken as a reference [6].
The basic module consists of three silicone chambers, which are parallel aligned and closed by end cap.Each of these chambers can be subjected to different air pressure.To suppress radial expansion of the chambers, rings as reinforcement were added.By preventing radial expansion, the chambers elongate according to the air pressure.If one or two chambers are pressurized, a bending of the module is achieved.With all three chambers under pressure, the entire module elongates.This allows the soft robot to move in all directions in space.
To ensure that the centerlines of the chambers remain parallel within a module, linking elements are used.In addition, this also keeps the chambers at a designated distance.Connector caps are used to connect chambers along the direction of expansion.The connector caps are designed so that the air pressure can be set separately for each individual chamber within a module.The end caps are used to close the soft robot [6].

Cosserat beam theory
Beam theories are useful for modeling slender structures like continuum robots.Thereby the three-dimensional robot is reduced to a one-dimensional curve  of length .Here, we use the Cosserat beam theory, presented in [7] to calculate the dynamics of the one-dimensional curve.This curve represents the centerline of the entire soft robot with where  is the arc length (0 ≤  ≤ ) and  the time.According to Figure 2A, second coordinate system is set up.This coordinate system is defined by the directors  1 ,  2 ,  3 .The directions of  1 and  2 coincide with the direction of the principal axes of inertia of the beam cross section and The rotation of the directors with respect to  1 ,  2 ,  3 can be realized by three successive rotations, summarized in cos  cos  + sin  sin  sin  sin  cos  − cos  sin  + sin  sin  cos  − sin  cos  + cos  sin  sin  cos  cos  sin  sin  + cos  sin  cos  cos  sin  − sin   cos Here  is the angle between  1 and  1 ,  is the angle between  2 and  2 and  is the angle between  3 and  3 .The dynamics of the beam can be described by an equilibrium of forces and moments ()()   =   (, ) +  (, ), (, ) =   (, ) + (, ) × (, ) + (, ), where () is the density, () is the cross sectional area and (, ) describes the angular momentum.(, ) is a vector of internal shear forces and normal force and (, ) is a vector of bending moments and twisting moment. (, ) and (, ) are external forces and moments per length.An external force can be, for example, the weight of an object to be carried.  and   give derivatives with respect to time, and   the derivative with respect to arc length.Furthermore, (, ) represents the tangent of the curve (, ) The quantities ,  and  can be calculated by the constitutive relations considering a linear elastic material behavior with Here,  is the Young's modulus,  is the shear modulus,   is the second area moment of inertia and  is a numerical factor depending on the shape of the cross section.The constitutive relations introduce the rotations  and the angular velocities  of the directors.The following relation is valid for the description of  and      (, ) = (, ) ×   (, ), (, ) = (, ) ×   (, ). (11)

Implementation in PyElastica
PyElastica is a Python implementation of Elastica, an open-source software package for simulating dynamic behavior of Cosserat rods [8,9].The beam theory presented in Section 2.2 is implemented in this simulation environment and is solved numerically.However, it is not possible to perform a static simulation in PyElastica directly.To find a static solution, it is necessary to add damping. Figure 3 illustrates the vertical movement of the end point of a horizontally clamped beam under gravitational load.Considering damping, the system can remain in the static rest position.
The following simulations for the workspace analysis are performed with respect to Figure 2 in the  1 - 3 plane.The vertical direction is the z-axis and the horizontal direction is the x-axis.Within the simulation, one end of the beam is fixed at a point on the z-axis.The propagation of the beam is in negative z-direction.To model the bending of the module due to the air pressure in the chambers, described in Section 2.1, the following relation between the bending moment , occurring due to air pressure in a chamber, and curvature  of the robot is used.The bending moment due to the air pressure is calculated from the force   generated by the chambers and the distance  between the chambers to the centerline of the soft robot with The force   is determined by the applied pressure and the cross-sectional area of the chamber.With the help of the curvature , it is possible to determine the shape of the centerline of the beam.
To explain the workflow of the simulation, the following describes the basic steps.First, the shape of the beam is predetermined by a given force   .Second, an external force   is applied to the free end of the beam in the negative z-direction.Finally, the deformation of the beam is simulated for the defined forces   and   .This procedure is performed for all possible combinations of   and   from a preset interval.

RESULTS AND DISCUSSION
To describe the reachable end positions of the soft robot with respect to its actuation pressure and a load, the simulation results are transformed into analytical equations.With these analytical equations, the dependencies of   and   on the achievable end coordinates of the beam can be determined.In Section 3.1, the steps to reach the analytical equations are presented.These equations are illustrated and described in Section 3.2.The simulations and results discussed in the following refer to the mechanical and geometrical properties of the soft robot of our institute (see Section 2.1).

Simulation process
In Section 2.3, it was described that for predefined forces   and external force   , the deformation was simulated.
Figure 4A shows simulation results of the beam deformations at a constant force   with variable forces   .From the results of these simulations, the coordinates of the end point of each beam were extracted.These coordinates are marked with a red circle in Figure 4A.For all combinations of   and   , this procedure was repeated.The description of the workspace is done by the reachable coordinates of the endpoint of the soft robot.The coordinates were summarized in an overview plot in Figure 4B.In this plot, each line represents a constant external force   .The filled circles indicate the force   in the chamber according to Figure 4A, which increases on the line to the right.Figure 4B was created by a studied interval from 0 to 10 N for   and from 0 to 6 N for   (Δ  = Δ  = 0.5 N).
To describe the workspace with analytical equations, the x-coordinates as well as the z-coordinates are considered separately.As a function of the forces   and   , the corresponding values are plotted in Figure 5.The blue points in Figures 5A  and 5B represent the simulated data.Using these data points, a polynomial is fitted.The order of this polynomial can be chosen arbitrarily.In Figure 5, a third-order polynomial was selected, which is composed as (  ,   ) =  0 +  1   +  2   +  3 When fitting this polynomial to the data points, the constants  0 to  9 are adjusted as best as possible with respect to the data points.The representation of the polynomial for the z-coordinate is analogous to Equation ( 14).

F I G U R E 5 (A) Extracted
x-coordinates from the simulations as blue points, the surface shows the fitted polynomial.(B) Extracted z-coordinates from the simulations as blue points, the surface shows the fitted polynomial.The order of the polynomial was chosen to be of third order.

TA B L E 1
Numerical values of the constants for the third-order polynomial of the x-coordinate from Figure 5A.

Workspace description
For the fitted polynomial of the x-coordinates shown in Equation ( 14), the coefficients  0 to  9 are given in Table 1.Table 2 contains the coefficients for the fitted polynomial of the z-coordinates.By inserting the numerical values into the equations,   and   remain as unknowns.For a given external force   , all reachable  −  coordinates can be calculated as a function of the force   .
From the numerical values of the constants, it can be concluded that   has mainly a linear influence on the xcoordinate.In comparison,   has a much stronger quadratic influence on the z-coordinate.  has a strong quadratic influence on the x-coordinate.On the z-coordinate, the influence predominates in a linear way.Analogous analysis can be performed for other soft robot configurations.To do this, the mechanical and/or geometrical properties of the soft robot must be changed in the simulation, as desired.Simulations for shorter soft robots revealed strong linear dependencies of   and   on the x-coordinate as well as on the z-coordinate.For longer soft robots, the influence shifts to higher order terms.

CONCLUSION
In this paper, an approach is proposed to determine the workspace of a soft robot.With this approach, it is also possible to identify the dependencies of the two input variables, air pressure inside a chamber and an external load, on the deformation of a soft robot.
In order to fully understand the approach, the soft robot used as an example was described at the beginning (see Section 2.1).Such slender structures are modeled by the Cosserat beam theory, which was explained in Section 2.2.The application of the simulation environment PyElastica was highlighted in Section 2.3.A description of the approach was provided in Section 3.
Due to its simple structure, the developed approach offers a fast evaluation of the deformation behavior of a soft robot.Using the analytical equations available from the approach, two properties can be determined with reduced effort.On the one hand, these are the reachable coordinates of the end point of the robot, which here describe the workspace.On the other hand, the dependencies of the forces   and   on the deformation can be obtained.The force   is directly related to the air pressure inside the chambers (see Section 2.3).  describes the external load on the soft robot.It is also possible to compare different soft robot configurations using the equations.Thus, it could be observed that for shorter soft robots, the dependence of   and   shifts toward linear terms.For longer soft robots, the dependencies shift to higher order terms.
A drawback of this approach is the dependence on simulations.If a property of the sot robot changes, such as the length, new simulations must be performed.It has also been found that extrapolations quickly exhibit high errors.Therefore, the intervals, in this case the intervals for   and   , in which the soft robot operates, must be known in advance.However, interpolations can be performed with a small error.
The presented approach offers a fast and simple solution to compare different soft robot configurations.Influences of the air pressure in the chambers as well as external loads can be read from analytical equations.Likewise, the analytical equations open an efficient option to control the air pressure for a given load.

F I G U R E 3
Difference between the beam movement without and with damping.(A) Without and (B) with damping.

F
I G U R E 4 (A) Deformation of a beam with constant external force and variable force of the chambers   and (B) endpoints of all beams from each combination of   and   .
Numerical values of the constants for the third-order polynomial of the z-coordinate from Figure5B.