First-Order Dynamic Modeling and Control of Soft Robots

1. Introduction

Soft robotic technologies are becoming increasingly prevalent in the design and development of robots (Kim et al., 2013). Subsequently, there has been growing interests in the modeling and control of soft bodied systems, Unlike robots designed with rigid components, soft robotic systems present novel challenges and opportunities in developing their controllers (George Thuruthel et al., 2018).

The most common modeling and control strategy for soft robots are based on steady-state models, which, under the steady-state assumption, can be equated to the kinematic model (George Thuruthel et al., 2018). For cylindrically-shaped soft robots, this leads to the popular constant curvature model (Webster and Jones, 2010). For other shapes, geometrically exact models or Finite Element Method have been proposed (Trivedi et al., 2008; Renda et al., 2012; Duriez, 2013; Gong et al., 2018). Machine learning techniques can also be used to develop such mappings in a model-free manner (Giorelli et al., 2013; George Thuruthel et al., 2017; Jiang et al., 2017). Refer to Sadati et al. (2017) for a detailed comparison into multiple static modeling techniques. Due to their steady-state assumptions, such controllers will, however, be limited in their reachability, efficiency, and speed. Therefore, controllers developed from dynamic models are much more desirable.

A popular method for developing dynamic models for soft robots is based on the cosserat-rod mechanics. Such models have been extensively used for soft robotic manipulators driven by tendon actuation (Rucker and Webster, 2011; Renda et al., 2014, 2018). For fluidic actuation, other models have been adopted (Marchese et al., 2016; Della Santina et al., 2019). Hybrid models based on lumped mass systems also looks promising for general soft robots (Sadati et al., 2019). However, all these models will be more computationally intensive than their static counterparts. Learning-based models are a promising alternative in such cases (Thuruthel et al., 2017; Gillespie et al., 2018). Nevertheless, deriving control strategies from dynamic models, in general, introduces additional complexities in motion planning. For an alternate approach that introduces a control-oriented modeling of soft robots, the readers are suggested to look into Della Santina and Rus (2019).

Unlike static controllers, developing fully-dynamic controllers would involve a planning stage. Typically, this has to be performed using some optimization techniques irrespective of the modeling strategy. A good example is the use of trajectory optimization for the control of a soft robotic manipulator using a model-based (Marchese et al., 2016) and model-free method (Thuruthel et al., 2017). This process is time consuming and hence debilitating for closed-loop dynamic control. For fully-actuated soft robots, closed-loop dynamic controllers can be developed in the configuration space (Della Santina et al., 2018). For task-space closed loop control, model-based reinforcement learning is the only viable solution till now, however, they tend to be highly task specific and time consuming (Thuruthel et al., 2018).

This article investigates the viability of a first-order dynamic model for soft robots. It must be noted that unlike state-space dimensionality reduction methods (Thieffry et al., 2018), we are reducing the temporal dimensionality of the dynamic model. Such a model reduction is based on the hypothesis that soft robots typically have high damping and low inertial properties. This makes it possible to approximate the second-order dynamic model to a first-order dynamic model by ignoring the inertial terms (Strogatz, 2001). Even in nature, the ubiquitous muscle dynamics can be modeled as a first-order dynamical system (Zajac, 1989). This model-order reduction provides two advantages. First, first-order dynamical systems are computationally cheaper than second-order dynamical systems. Second, it opens the possibility to develop novel closed-loop control strategies using the reduced-order state feedback. Here, we show the direct learning of the operational space dynamics of the first-order dynamic model. Due to simplifying step, controllers can be easily developed using machine learning and a simple path planning algorithm. Moreover, the sensory requirements for closed-loop dynamic control is reduced because of the simplification.

We investigate the viability of this simplifying assumption using extensive simulation studies. First, we present the theoretical reasoning behind the first-order assumption and its corresponding controller. Then we briefly introduce the fully dynamic simulation model that is used to verify the learned forward models and the dynamic controller. Finally, we present details on the learning architecture together with results of the model and the closed-loop task-space controller.

2. Theory

Given a soft robot that can be kinematically modeled by the configuration-space q ∈ ℝⁿ, the task-space variable can be obtained by the kinematic transformation:

$\begin{array}{ll}x=F(q)& (1)\end{array}$

Where, x ∈ ℝ^m and m ≤ n. The task-space variable is typically the pose of the end-effector and is to be controlled. The configuration space is the set of independent variables that fully defines the state of the robot. The fully dynamic model of the soft robot can then be represented using the configuration-space variable as:

$\begin{array}{ll}M(q)\ddot{q}+C(q,\dot{q})\dot{q}+G(q)=\tau & (2)\end{array}$

Here, M(q) represents the inertial properties, $C(q,\dot{q})$ combines the coriolis, centrifugal and damping elements, G(q) represents the gravitational and stiffness effects and τ is the generalized force applied internally by the robot.

Soft robots typically have high damping values with low inertial properties. This is because they are commonly fabricated with viscoelastic materials with low material density. After the initial transient motion of a soft robot from rest (when $\dot{q}=0$), the first order term dominates the second order term [i.e., $C(q,\dot{q})\dot{q}>>M(q)\ddot{q}$]. Hence, the second-order term can be ignored without sacrificing the accuracy of the model (Zajac, 1989). The dominant modeling error will occur during the initial transient motion (Strogatz, 2001). This transforms the second-order dynamical model (Equation 2) into:

$\begin{array}{ll}C(q,\dot{q})\dot{q}+G(q)=\tau & (3)\end{array}$

After discretizing the equation, the dynamic model can now be represented through the mapping : (q_i, τ_i) → q_i+1. Where, q_i, q_i+1 are the current and the next configuration of the soft robot, respectively. Correspondingly, this implies that a closed-loop dynamic controller would require only the zero-order state feedback (q) for control.

2.1. Controller Design

The obtained first-order dynamical Equation (3), the first-order configuration-space term can be replaced using the well known inverse kinematics mapping:

$\begin{array}{l}ẋ=J(q)\dot{q}\end{array}$

$\begin{array}{ll}\dot{q}=J{(q)}^{†}ẋ& (4)\end{array}$

Where, J(q) is Jacobian matrix and J(q)^† is any generalized inverse matrix. Note that we have ignored the null-space terms for brevity. Now the first-order dynamic equation can be reformulated as function of the task-space variables:

$\begin{array}{ll}\tau =F(q,ẋ)& (5)\end{array}$

This functional mapping can now be directly learning by a machine learning architecture. We can now transform the mapping and introduce the target task-space variable x^d as:

$\begin{array}{ll}ẋ=x^d-x^c& (6)\end{array}$

The operational space controller mapping now becomes:

$\begin{array}{ll}\tau =F(q^c,x^d)& (7)\end{array}$

Here, (q^c, x^c) are the current configuration-space coordinates and the current task-space coordinates, respectively. For the special case when the cardinality of the configuration-space coordinate is the same as the cardinality of the task-space coordinate (mapping from q → x is bijective), the operational-space controller can be further simplified to a mapping: (x^c, x^d) → τ. A simple feedforward neural network can be used to learn this mapping (Figure 2). In this paper, we restrict our studies to this condition for simplicity. This allows us to test the learned controller by providing trajectories in the task-space without the need to solve the inverse kinematics problem. For redundant task-space controllers (non-bijective mappings), additional planning stages might be required to obtain the configuration-space trajectories. In other words, the control trajectory cannot be represented only in the task-space variables. In such a case, an augmented trajectory can be defined with the task-space variable along with a optimization routine to check for kinematic constraints as shown before in Thuruthel et al. (2017). Note that the trajectories are only zero-order task space variables and the feedback required for closed-loop control is also zero-order. This greatly reduces the requirement on the sensors, the effect of noise and the complexity of the trajectory planner.

3. Simulation Model

The dynamic model is based on the Piece-wise Constant Strain (PCS) approach for soft-rigid multibody system of Renda and Seneviratne (2018) (see Figure 1). In the following, all the quantities are expressed in the local (body) coordinate frame if not specified otherwise. The superscript ′ and ˙ represent partial differentiation with respect to the space variable and time variables, respectively. The accent ^˜ represents the usual isomorphism between a vector in ℝ³ and its corresponding skew-symmetric matrix in 𝔰𝔬(3).

FIGURE 1

Figure 1. Schematic of the soft manipulator used for our studies. The manipulator is driven by three tendons arranged in the configuration shown above. Unless stated otherwise, the single section manipulator is used for the study.

3.1. Kinematics

The relative position and orientation of a soft body i with respect to its predecessor in the chain is defined as a curve g_i(·) : X ∈ [0, L_i] ↦ g_i(X) ∈ SE(3) with

$g_i(X)=\left(\begin{array}{cc}{\text{R}}_i& {\text{u}}_i\\ 0^T& 1\end{array}\right).$

The continuous models of the position, velocity and acceleration of a soft body can be derived from the Cosserat rod theory, which gives (Boyer and Renda, 2016):

$\begin{array}{ll}{g}_i^{\prime }=g_i{\hat{\xi }}_i,& (\mathrm{8a})\end{array}$

$\begin{array}{ll}{\eta }_i^{\prime }={\dot{\xi }}_i-\text{a}{\text{d}}_{{\xi }_i}{\eta }_i,& (\mathrm{8b})\end{array}$

$\begin{array}{ll}{\dot{\eta }}_i^{\prime }={\ddot{\xi }}_i-\text{a}{\text{d}}_{{\dot{\xi }}_i}{\eta }_i-\text{a}{\text{d}}_{{\xi }_i}{\dot{\eta }}_i,& (\mathrm{8c})\end{array}$

where ${\xi }_i(X)\in {ℝ}^6$ defines the strain state, ${\eta }_i(X)\in {ℝ}^6$ is the cross-section velocity twist and $\text{a}{\text{d}}_{(·)}\in {ℝ}^{6\times 6}$ is the adjoint operator of the Lie algebra (see Nomenclature). Going further into detail, we have

${\hat{\xi }}_i(X)=\left(\begin{array}{cc}{\tilde{\text{k}}}_i& {\text{p}}_i\\ 0^T& 0\end{array}\right)\in 𝔰𝔢(3),{\xi }_i(X)={\left({\text{k}}_i^T,{{\text{p}}_i}^T\right)}^T\in {ℝ}^6,$

${\hat{\eta }}_i(X)=\left(\begin{array}{cc}{\tilde{\text{w}}}_i& {\text{v}}_i\\ 0^T& 0\end{array}\right)\in 𝔰𝔢(3),{\eta }_i(X)={\left({\text{w}}_i^T,{\text{v}}_i^T\right)}^T\in {ℝ}^6,$

with ${\text{k}}_i(X)\in {ℝ}^3$ and ${\text{p}}_i(X)\in {ℝ}^3$ the angular and linear strain; and ${\text{w}}_i(X)\in {ℝ}^3$ and ${\text{v}}_i(X)\in {ℝ}^3$ the angular and linear velocity, respectively.

To model constrained rod, such as the Kirchhoff-Love case with angular strain only, the strain field is specified as:

${\xi }_i=B_i{q}_i+{\xi }_i^{*},$

where $B_i\in {ℝ}^{6\times n_i}$ forms a basis for the allowed motion subspace, $q_i\in {ℝ}^{n_i}$ contains the values of the allowed strains and, ${\xi }_i^{*}\in {ℝ}^6$ is the reference twist modeling the reference shape.

Assuming piece-wise constant strains (Renda et al., 2016), Equations (8) can be analytically integrated using the matrix exponential method, leading to:

$\begin{array}{ll}{g}_i(X)=e^{X{\hat{\xi }}_i},& (\mathrm{9a})\end{array}$

$\begin{array}{ll}{\eta }_i(X)={\text{Ad}}_{g_i}^{-1}{\eta }_h+{\text{Ad}}_{g_i}^{-1}{\text{T}}_{g_i}{B}_i{\dot{q}}_i,& (\mathrm{9b})\end{array}$

$\begin{array}{ll}{\dot{\eta }}_i(X)={\text{Ad}}_{g_i}^{-1}{\dot{\eta }}_h+{\text{Ad}}_{g_i}^{-1}{\displaystyle {\int }_0^X}\text{A}{\text{d}}_{g_i(s)}\text{a}{\text{d}}_{{\eta }_i(s)}ds{B}_i{\dot{q}}_i+\text{A}{\text{d}}_{g_i}^{-1}{\text{T}}_{g_i}{B}_i{\ddot{q}}_i,& (\mathrm{9c})\end{array}$

where $\text{A}{\text{d}}_{g_i}(X)\in {ℝ}^{6\times 6}$ is the Adjoint operator of SE(3), and ${\text{T}}_{g_i}(X)={\int }_0^X{e}^{s\text{a}{\text{d}}_{{\xi }_i}}ds$ is the tangent operator of the exponential map, of which an analytic expression, derived from (Selig, 2007), is given in the Nomenclature.

Successive applications of the kinematics (Equation 9) for all the bodies of the system, yields to the definition of the geometric Jacobian $J_i(q,X)\in {ℝ}^{6\times n}$ and its derivative ${\dot{J}}_i(q,\dot{q},X)\in {ℝ}^{6\times n}$ (n being the total number of DOFs), which relates the generalized coordinate vector $q={\left[q_1^T{q}_1^T\cdots q_N^T\right]}^T\in {ℝ}^n$ (N being the total number of bodies) and the velocity twist η_i(X), for each soft body i, as shown below.

$\begin{array}{ll}{\eta }_i(X)={\displaystyle \sum _{h=0}^i}\text{A}{\text{d}}_{g_h\cdots g_i}^{-1}{\text{T}}_{g_h}{B}_h{\dot{q}}_h={\displaystyle \sum _{h=0}^i}{}^{i}S_h{\dot{q}}_h=J_i(q,X)\dot{q},& (\mathrm{10a})\end{array}$

$\begin{array}{l}{\dot{\eta }}_i(X)={\displaystyle \sum _{h=0}^i}{}^{i}S_h{\ddot{q}}_h+\text{A}{\text{d}}_{g_h\cdots g_i}^{-1}{\displaystyle {\int }_0^X}\text{A}{\text{d}}_{g_h(s)}\text{a}{\text{d}}_{{\eta }_h(s)}ds{B}_h{\dot{q}}_h\\ ={\displaystyle \sum _{h=0}^i}{}^{i}S_h{\ddot{q}}_h{+}^i{\dot{S}}_h{\dot{q}}_h=J_i(q,X)\ddot{q}+{\dot{J}}_i(q,\dot{q},X)\dot{q},& (\mathrm{10b})\end{array}$

where the block elements of the i^th Jacobian ${}^i{S}_{(·)}\in {ℝ}^{6\times n_{(·)}}$ and its derivative ${}^i{\dot{S}}_{(·)}\in {ℝ}^{6\times n_{(·)}}$ have been defined. Note that the last three rows of Equation (10a) provide an analytical expression of the kinematics map required by Equation (4).

3.2. Dynamics

Once a Jacobian is found, the generalized dynamics of the system can be obtained by projecting the free dynamics of each soft body by virtue of the D'Alembert's principle. The free dynamic equation, with its boundary conditions, of a soft body is given by (Renda et al., 2018):

$\begin{array}{l}{{M}}_i{\dot{\eta }}_i+\text{a}{\text{d}}_{{\eta }_i}^{*}{{M}}_i{\eta }_i={\left({{F}}_{i_i}-{{F}}_{a_i}\right)}^{\prime }+\text{a}{\text{d}}_{{\xi }_i}^{*}\left({{F}}_{i_i}-{{F}}_{a_i}\right)+{\bar{{F}}}_{e_i},\\ \left({{F}}_{i_i}-{{F}}_{a_i}\right)(0)=-{{F}}_{J_i},\left({{F}}_{i_i}-{{F}}_{a_i}\right)\left(L_i\right)=-\text{A}{\text{d}}_{g_{ij}}^{*}{{F}}_{J_j}\text{;}& (11)\end{array}$

where ${{M}}_i(X)=\mathrm{\text{diag}}(J_{x_i},J_{y_i},J_{z_i},A_i,A_i,A_i){\rho }_i\in {ℝ}^{6\times 6}$ is the screw inertia matrix of the cross-section (J_·i(X) being the second moment of area about the axis · and A_i(X) the area of the cross-section); ${\bar{{F}}}_{e_i}(X)\in {ℝ}^6$ is the distributed external load; ${{F}}_{a_i}(X)\in {ℝ}^6$ is the internal wrench due to the distributed actuation (Renda et al., 2017); ${{F}}_{i_i}(X)\in {ℝ}^6$ is the internal wrench due to the elasticity of the soft body; ${{F}}_{J_{(·)}}\in {ℝ}^6$ is the wrench transmitted across joint (·) and $\text{a}{\text{d}}_{(·)}^{*}$ (respectively $\text{A}{\text{d}}_{(·)}^{*}$) ∈ ℝ^6×6 is the co-adjoint (respectively co-Adjoint) map of the Lie algebra (respectively Lie group) defined in Nomenclature. Regarding the internal elastic force, a linear viscoelastic constitutive model is usually chosen:

$\begin{array}{ll}{{F}}_{i_i}(X)={\Sigma }_i\left({\xi }_i-{\xi }^{*}\right)+{{\rm Y}}_i{\dot{\xi }}_i={\Sigma }_i{B}_i{q}_i+{{\rm Y}}_i{B}_i\dot{q_i},& (12)\end{array}$

where

$\begin{array}{l}{\Sigma }_i(X)=\text{diag}(G_i{J}_{x_i},E_i{J}_{y_i},E_i{J}_{z_i},E_i{A}_i,G_i{A}_i,G_i{A}_i)\in {ℝ}^{6\times 6}\text{ and}\end{array}$

$\begin{array}{l}{{\rm Y}}_i(X)=\text{diag}(J_{x_i},3J_{y_i},3J_{z_i},3A_i,A_i,A_i){\nu }_i\in {ℝ}^{6\times 6}\end{array}$

are the screw stiffness and viscosity matrix (E_i being the young modulus, G_i the shear modulus and ν_i the shear viscosity).

By Jacobian projection of the free dynamics (Equation 11), we obtain the generalized dynamics in its classical configuration-space form:

$\begin{array}{ll}M(q)\ddot{q}+(C(q,\dot{q})+D)\dot{q}+Kq=\tau (q)+F(q),& (13)\end{array}$

where M ∈ ℝ^n×n is the generalized mass matrix, C ∈ ℝ^n×n is the generalized Coriolis matrix, D ∈ ℝ^n×n is the block-diagonal generalized damping matrix, K ∈ ℝ^n×n is the block-diagonal generalized stiffness matrix, F ∈ ℝⁿ is the vector of generalized position-dependent external forces and τ ∈ ℝⁿ is the vector of applied actuation forces. Note that the dynamic Equation (13) can be written in the form required by Equation (2), with $C(q,\dot{q})=C(q,\dot{q})+D$ and G(q) = Kq − F(q). Going further into details, the coefficient matrices take the form:

$\begin{array}{ll}M(q)={\displaystyle \sum _{i=1}^N}{\displaystyle {\int }_0^{L_i}}{J}_i^T{{M}}_i{J}_{i}dX,& (\mathrm{14a})\end{array}$

$\begin{array}{ll}C\left(q,\dot{q}\right)={\displaystyle \sum _{i=1}^N}{\displaystyle {\int }_0^{L_i}}{J}_i^T\left(\text{a}{\text{d}}_{J_i\dot{q}}^{*}{{M}}_i{J}_i+{{M}}_i{\dot{J}}_i\right)dX& (\mathrm{14b})\end{array}$

$\begin{array}{l}D=\text{diag}(B_1^T{\displaystyle {\int }_0^{L_1}}{{\rm Y}}_{1}dX{B}_1,B_2^T{\displaystyle {\int }_0^{L_2}}{{\rm Y}}_{2}dX{B}_2,\cdots ,\\ B_N^T{\displaystyle {\int }_0^{L_N}}{{\rm Y}}_{N}dX{B}_N),& (\mathrm{14c})\end{array}$

$\begin{array}{l}K=\text{diag}(B_1^T{\displaystyle {\int }_0^{L_1}}{\Sigma }_{1}dX{B}_1,B_2^T{\displaystyle {\int }_0^{L_2}}{\Sigma }_{2}dX{B}_2,\cdots ,\\ B_N^T{\displaystyle {\int }_0^{L_N}}{\Sigma }_{N}dX{B}_N),& (\mathrm{14d})\end{array}$

$\begin{array}{l}\tau \left(q\right)=[{\left(B_1^T{\displaystyle {\int }_0^{L_1}}{{F}}_{a_1}dX\right)}^T{\left(B_2^T{\displaystyle {\int }_0^{L_2}}{{F}}_{a_2}dX\right)}^T\cdots \\ {{\left(B_N^T{\displaystyle {\int }_0^{L_N}}{{F}}_{a_N}dX\right)}^T]}^T,& (14e)\end{array}$

$\begin{array}{ll}F\left(q\right)={\displaystyle \sum _{i=1}^N}{\displaystyle {\int }_0^{L_i}}{J}_i^T{\bar{{F}}}_{e_i}\text{.}& (\mathrm{14f})\end{array}$

It is worth noting here the different structure of the components of the generalized dynamics Equation (13). Similarly to the minimal Lagrangian models of traditional rigid robots, inertial loads are characterized by full coefficient matrices, as can be see from Equations (14a) and (14b), while damping and stiffness loads are characterized by block-diagonal coefficient matrices, as for Equations (14c) and (14d). This is in contrast with other modeling approaches that use absolute coordinates, such as Finite Elements, for which the opposite holds. Inertial coefficient matrices are block-diagonal while damping and stiffness coefficient matrices are full. Thus, neglecting inertial terms will be well suited for minimal Lagrangian models for soft robotic manipulator, such as the PCS approach.

4. Methods and Results

This section investigates two studies. First, we validate the accuracy of the learned first-order model with respect to the learned second-order model. Second, we perform simulated experiments to validate the accuracy of the proposed controller. All the tests are performed on the fully dynamic model described in section 3.

4.1. Dynamic Modeling

The learned models are derived using a kind of recurrent neural network called a nonlinear autoregressive exogenous (NARX) model (Billings, 2013). The NARX network is particularly suited for our study as it allows us to define the feedback horizon of the recurrent connection explicitly (Figure 2). In other words, we can ensure that the neural network receives only zeroth-order feedback for the first-order model and the appended first-order feedback for the second-order model. For the single section soft manipulator, the configuration-space (q) is equivalent to the task-space variable (x), which is defined as the three-dimensional position of the end-effector. We use a recurrent network with one hidden layer with a size of 40 for both the first-order and second-order model, for a fair comparison. The training parameters of the NARX network is given in Table 1.

FIGURE 2

Figure 2. Modeling of the forward dynamics and the inverse dynamics controller. Note that for the non-redundant cases, the configuration-space (q) can be interchanged with the task-space variable (x), as done in this paper.

TABLE 1

Table 1. Parameters of the learned forward dynamic model.

Random actuation of the tendons are performed (motor babbling) for 70 s to obtain the samples for learning the forward model (Figure 3). Specifically, random actuation inputs are used to drive the manipulator and the corresponding actuator inputs and end-effector position is recorded over time. The same data samples can also be used for developing the closed-loop controller. Note that the dynamic workspace is concentrated along the direction of the three actuators. This is because tendons in tension have a strong attractor behavior. Hence, it will be easier to move along this direction. When testing our controller, careful measure is taken to ensure that our trajectories pass across this workspace regions.

FIGURE 3

Figure 3. Dynamic workspace of the manipulator obtained by motor babbling. This is obtained by recording end-effector position of the manipulator when actuated by random continuous actuation signals.

Figures 4, 5 show the performance of the learned model in comparison to the original cosserat model for a step and periodic response, respectively. Note the higher errors in the first-order model in the beginning of motion for the step response. Since the inertial effects are ignored, it is also visible that oscillations caused by overshoot is not found in the first-order model. However, the steady state error, with respect to the second-order model is relatively small. For the periodic excitation case, the difference between the first-order model and the second-order model is almost non-existent in the relevant coordinates. This is as expected since our approximations are more valid when the manipulator is in a non-stationary state.

FIGURE 4

Figure 4. Step response of the learned models in comparison to the actual analytical model. The step signal is given to single actuator at time = 2 s.

FIGURE 5

Figure 5. Periodic response of the learned models in comparison to the actual analytical model.

In order to analyze the effects of the viscosity and the inertial effects on our modeling assumption, we further perform studies on the accuracy of the first-order model for varying material viscosity and density. As material density increases and the material viscosity decreases, the inertial effects become more and more dominant. Hence, one would expect the accuracy of the first-order model to decrease and the second-order model to remain constant. However, this is not necessarily the case as the training of second-order recurrent neural networks is more prone to instabilities (Pascanu et al., 2013).

Figure 6 shows how the root sum squared (RSS) error of the first-order model is affected when the material properties of the soft arm is changed in a way that weakens our main assumption. The material density is increased up to a factor of 2 and the material viscosity is reduced by a factor of 6. The motor babbling inputs and the neural network parameters are kept the same for the tests. It is clear from the prediction errors that the second-order model always performs better than the first-order model. However, the change in accuracy of both the models are not affected significantly by the change in material properties. Note that the initial parameters of the simulated cosserat model soft arm were obtained from real experiments on an Octopus-inspired soft manipulator, which was manufactured with silicone and driven by tendons (Renda et al., 2014).

FIGURE 6

Figure 6. Root sum squared error of the learned model for varying dynamic properties of the manipulator.

Increasing the length of the manipulator is another way to increase the inertial properties and weaken the first-order assumption. For this, we test the same methodology on a 4 section manipulator, however, actuated only on the first section. Each section has approximately the same length, with the total length of the manipulator adding to 418 mm. We test two designs, one with a tapered morphology and the other with a cylindrical morphology (higher inertial properties). The radius of the tapered morphology linearly reduces from 30 to 10 mm while the cylindrical morphology has a fixed diameter of 30 mm. For this test and the following controller results, the default parameters of the manipulator is used (i.e., with the material density of 1,080 kg/m³ and viscosity of 300 Pas; Renda et al., 2014). The results of our forward dynamics prediction on both the systems are shown in Figure 7. Contrary to our expectation, increasing the length of the manipulator did not affect the performance of the first-order model when compared to the second-order model. We believe this is because the inertial effects are still compensated by the medium (water) in which the manipulator is surrounded in. Our previous studies have shown that the dynamics of the manipulator becomes chaotic and hence unpredictable without a surrounding medium and when the length increases (Thuruthel et al., 2019). Based on our results, it can be deduced that this is more because of the first-order terms rather than the second-order terms. It could be because of the increased length, gravity, centripetal/centrifugal forces, higher DoFs etc. Note that this is a limitation of learning-based approaches as increased sensitivity to initial conditions decrease the stability of the training process and hence the performance of the model. The same applies to the parameter tuning process for analytical models. It could be concluded that it is a good practice to reduce the inertial effects of a soft-bodied system to improve the predictability of its dynamics, irrespective of the order of the model.

FIGURE 7

Figure 7. Root sum squared error of the learned model for the 4-section underactuated manipulator.

4.2. First-Order Dynamic Controller

The closed-loop task-space controller is derived by learning the operational-space dynamics mapping as described in section 2.1. As the mapping is not recursive, it can be learned using a simple feedforward neural network, as shown in Figure 2. For our non-redundant case, we can replace the configuration-space variable, q, with the task-space variable, x. The samples for learning the mapping is obtained through the same motor babbling process as described in the previous section. For training the controller the mapping is defined as: (xⁱ, xⁱ⁺¹) → τ. When testing the controller, the next task-space coordinate, xⁱ⁺¹, is replaced by the desired task-space variable x^d. The parameters of the neural network used for learning the first-order inverse dynamics mapping is given in Table 2.

TABLE 2

Table 2. Parameters of the first-order inverse dynamics controller.

Path planning is usually a complex problem in inverse dynamics based controllers, but as our inverse model is developed with only zeroth order state feedback, the development of the task-space trajectory is greatly simplified. Acceptable paths can easily be generated using the data points obtained from the workspace of the manipulator, which is obtained during the motor babbling phase. The desired paths can be generated by picking reachable points from the workspace and routing a path through them, ensuing that there is sufficient time for the manipulator to reach adjacent points. This can be easily done by fixing a cap on the maximum distance between adjacent task-space variables.

To test our controller, we generate randomized linear paths for the end-effector of the manipulator to follow. This is done by picking two random points from the robot workspace and linearly interpolating a trajectory between them and from the initial position of the end-effector. If the initial position of the end-effector is p⁰ and the two random points are p¹ and p², the generated path is from p⁰ → p¹ → p² → p¹. Note that the intermediate points are not necessarily reachable by this naive approach. Accurate trajectories can be generated by searching for adjacent points in the workspace or by projecting the trajectory onto the workspace surface.

The results of the trajectory tracking is shown in Figure 8. The performance of the controller is excellent considering the fact that it is myopic with no step-ahead planning and the task-space trajectory being non-optimal. Due to the low computational cost in running the inverse-dynamic controller, we are able to run closed-loop task-space controller at a very high control frequency of 100 Hz. This will also allow the controller to compensate for any modeling errors incurred by the approximation.

FIGURE 8

Figure 8. Controller results for a continuous path. The motion of the end-effector is shown here.

The quality of the task-space trajectory can be analyzed by running the controller in open-loop. This can be done by assuming the manipulator is able to reach each trajectory point perfectly and obtaining the best control action at each time step. The results of such a scenario is shown in Figure 9, along with the corresponding control inputs. As expected, the open-loop controller performs worse than the closed-loop controller even though there are no external disturbances in the simulation. The same scenario is repeated with the controller inactive for the first 0.2 s in Figure 9. As the desired targets are now farther from the current position of the manipulator, it is not necessary that the controller is able to follow the trajectory accurately. However, as evident from the results, the controller is able to recover and remarkably converge to the same solution as the original closed-loop controller. This also shows how important it is to close-the-loop, even at the cost of reduced model accuracy. The tracking error of the three tests described in this scenario is shown in Figure 10.

FIGURE 9

Figure 9. Testing the first-order controller in open-loop, closed-loop, and with a late start.

FIGURE 10

Figure 10. Tracking error for the scenario shown in Figure 9.

5. Conclusion

This paper presents and verifies a model simplifying assumption for soft robots. The core idea of the assumption is that soft robots, by definition, tend to have low inertial and high viscoelastic properties. This leads to dynamic behaviors which are well approximated by a first-order system, as typically observed in nature. We verify this assumption using a simulated fully dynamic model of an Octopus-like manipulator and a type of recurrent neural network called NARX network. Finally, we develop easy-to-develop closed-loop task-space dynamic controllers based on this assumption. Our results indicate that controllers developed on this assumption can compensate the errors in modeling accuracy with the increased control frequency. Our method makes path planning simpler for non-redundant cases. Moreover, the sensory requirements for closed-loop dynamic control is reduced because of the simplification. This is because the state feedback required for the controller is only the zeroth-order component. Our work also indicates that the additional modeling complexity that soft elements introduce can, to some extent, be reduced by designing low inertial highly visco-elastic soft robot designs.

Although we use machine learning tools to test our modeling assumption and develop our controller, the approach is equally suited for analytical approaches. In fact, ignoring the inertial elements would greatly simply the modeling and parameter estimation process involved in model-based control of soft robots. Not only can we reduce the states of the dynamical system, but we also avoid the problem of estimating and inverting the full mass matrix (see section 3.2). However, it must be kept in mind that first-order systems present numerical challenges in their implementation. Such problems are not found in a learning based approach and hence desirable in that respect. Typically, the first-order model leads to a stiff differential equation and requires specialized techniques for solving them. Interested readers are suggested to refer to Strogatz (2001) for further information. Future work involve extending the work to redundant systems.

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/s.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Funding

This publication was partially supported by Khalifa University of Science and Technology under Award Nos. FSU-2018-08, RC1-2018-KUCARS and Research Excellence (AARE) 2018-05 and by the SHERO project, a Future and Emerging Technologies (FET) programme of the European Commission (grant agreement ID 828818).

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.

References

Billings, S. A. (2013). Nonlinear System Identification: NARMAX Methods in the Time, Frequency, and Spatio-Temporal Domains. Chichester: John Wiley & Sons.

Google Scholar

Boyer, F., and Renda, F. (2016). Poincare's equations for cosserat media: application to shells. J. Nonlin. Sci. 27, 1–44. doi: 10.1007/s00332-016-9324-7

CrossRef Full Text | Google Scholar

Della Santina, C., Katzschmann, R. K., Bicchi, A., and Rus, D. (2019). Model-based dynamic feedback control of a planar soft robot: trajectory tracking and interaction with the environment. Int. J. Robot. Res. 39, 490–513. doi: 10.1177/0278364919897292

CrossRef Full Text | Google Scholar

Della Santina, C., Katzschmann, R. K., Biechi, A., and Rus, D. (2018). “Dynamic control of soft robots interacting with the environment,” in 2018 IEEE International Conference on Soft Robotics (Livorno: RoboSoft), 46–53.

Google Scholar

Della Santina, C., and Rus, D. (2019). Control oriented modeling of soft robots: the polynomial curvature case. IEEE Robot. Autom. Lett. 5, 290–298. doi: 10.1109/LRA.2019.2955936

CrossRef Full Text | Google Scholar

Duriez, C. (2013). “Control of elastic soft robots based on real-time finite element method,” in 2013 IEEE International Conference on Robotics and Automation (Karlsruhe), 3982–3987.

Google Scholar

George Thuruthel, T., Ansari, Y., Falotico, E., and Laschi, C. (2018). Control strategies for soft robotic manipulators: a survey. Soft Robot. 5, 149–163. doi: 10.1089/soro.2017.0007

PubMed Abstract | CrossRef Full Text | Google Scholar

George Thuruthel, T., Falotico, E., Manti, M., Pratesi, A., Cianchetti, M., and Laschi, C. (2017). Learning closed loop kinematic controllers for continuum manipulators in unstructured environments. Soft Robot. 4, 285–296. doi: 10.1089/soro.2016.0051

PubMed Abstract | CrossRef Full Text | Google Scholar

Gillespie, M. T., Best, C. M., Townsend, E. C., Wingate, D., and Killpack, M. D. (2018). “Learning nonlinear dynamic models of soft robots for model predictive control with neural networks,” in 2018 IEEE International Conference on Soft Robotics (Livorno: RoboSoft), 39–45.

Google Scholar

Giorelli, M., Renda, F., Ferri, G., and Laschi, C. (2013). “A feed-forward neural network learning the inverse kinetics of a soft cable-driven manipulator moving in three-dimensional space,” in 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems (Tokyo), 5033–5039.

Google Scholar

Gong, Z., Cheng, J., Hu, K., Wang, T., and Wen, L. (2018). “An inverse kinematics method of a soft robotic arm with three-dimensional locomotion for underwater manipulation,” in 2018 IEEE International Conference on Soft Robotics (Livorno: RoboSoft), 516–521.

Google Scholar

Jiang, H., Wang, Z., Liu, X., Chen, X., Jin, Y., You, X., et al. (2017). “A two-level approach for solving the inverse kinematics of an extensible soft arm considering viscoelastic behavior,” in 2017 IEEE International Conference on Robotics and Automation (ICRA) (Singapore), 6127–6133.

Google Scholar

Kim, S., Laschi, C., and Trimmer, B. (2013). Soft robotics: a bioinspired evolution in robotics. Trends Biotechnol. 31, 287–294. doi: 10.1016/j.tibtech.2013.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Marchese, A. D., Tedrake, R., and Rus, D. (2016). Dynamics and trajectory optimization for a soft spatial fluidic elastomer manipulator. Int. J. Robot. Res. 35, 1000–1019. doi: 10.1109/ICRA.2015.7139538

CrossRef Full Text | Google Scholar

Pascanu, R., Mikolov, T., and Bengio, Y. (2013). “On the difficulty of training recurrent neural networks,” in International Conference on Machine Learning (Atlanta, GA), 1310–1318.

Google Scholar

Renda, F., Boyer, F., Dias, J., and Seneviratne, L. (2018). Discrete cosserat approach for multisection soft manipulator dynamics. IEEE Trans. Robot. 34, 1518–1533. doi: 10.1109/TRO.2018.2868815

CrossRef Full Text | Google Scholar

Renda, F., Cacucciolo, V., Dias, J., and Seneviratne, L. (2016). “Discrete cosserat approach for soft robot dynamics: a new piece-wise constant strain model with torsion and shears,” in 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) (Daejeon), 5495–5502.

Google Scholar

Renda, F., Cianchetti, M., Abidi, H., Dias, J., and Seneviratne, L. (2017). Screw-Based Modeling of Soft Manipulators With Tendon and Fluidic Actuation. J. Mechanisms Robotics. 9(4):041012. doi: 10.1115/1.4036579

CrossRef Full Text | Google Scholar

Renda, F., Cianchetti, M., Giorelli, M., Arienti, A., and Laschi, C. (2012). A 3D steady-state model of a tendon-driven continuum soft manipulator inspired by the octopus arm. Bioinspir. Biomimet. 7:025006. doi: 10.1088/1748-3182/7/2/025006

PubMed Abstract | CrossRef Full Text | Google Scholar

Renda, F., Giorelli, M., Calisti, M., Cianchetti, M., and Laschi, C. (2014). Dynamic model of a multibending soft robot arm driven by cables. IEEE Trans. Robot. 30, 1109–1122.

Google Scholar

Renda, F., and Seneviratne, L. (2018). “A geometric and unified approach for modeling soft-rigid multi-body systems with lumped and distributed degrees of freedom,” in 2018 IEEE International Conference on Robotics and Automation (ICRA) (Brisbane, QLD), 1567–1574. doi: 10.1109/ICRA.2018.846118

CrossRef Full Text | Google Scholar

Rucker, D. C., and Webster, R. J. III. (2011). Statics and dynamics of continuum robots with general tendon routing and external loading. IEEE Trans. Robot. 27, 1033–1044.

Google Scholar

Sadati, S. H., Naghibi, S. E., Shiva, A., Michael, B., Renson, L., Howard, M., et al. (2019). TMTDyn: A matlab package for modeling and control of hybrid rigid–continuum robots based on discretized lumped systems and reduced-order models. Int. J. Robot. Res. doi: 10.1177/0278364919881685. [Epub ahead of print].

CrossRef Full Text | Google Scholar

Sadati, S. H., Naghibi, S. E., Shiva, A., Walker, I. D., Althoefer, K., and Nanayakkara, T. (2017). “Mechanics of continuum manipulators, a comparative study of five methods with experiments,” in Annual Conference Towards Autonomous Robotic Systems (Guildford, UK: Springer), 686–702.

Google Scholar

Selig, J. (2007). Geometric Fundamentals of Robotics. Monographs in Computer Science. New York, NY: Springer.

Google Scholar

Strogatz, S. (2001). Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (Studies in Nonlinearity). Boca Raton, FL: CRC press.

Thieffry, M., Kruszewski, A., Duriez, C., and Guerra, T.-M. (2018). Control design for soft robots based on reduced-order model. IEEE Robot. Autom. Lett. 4, 25–32. doi: 10.1109/LRA.2018.2876734

CrossRef Full Text | Google Scholar

Thuruthel, T. G., Falotico, E., Renda, F., Flash, T., and Laschi, C. (2019). Emergence of behavior through morphology: a case study on an octopus inspired manipulator. Bioinspir. Biomimet. 14:034001. doi: 10.1088/1748-3190/ab1621

PubMed Abstract | CrossRef Full Text | Google Scholar

Thuruthel, T. G., Falotico, E., Renda, F., and Laschi, C. (2017). Learning dynamic models for open loop predictive control of soft robotic manipulators. Bioinspir. Biomimet. 12:066003. doi: 10.1088/1748-3190/aa839f

PubMed Abstract | CrossRef Full Text | Google Scholar

Thuruthel, T. G., Falotico, E., Renda, F., and Laschi, C. (2018). Model-based reinforcement learning for closed-loop dynamic control of soft robotic manipulators. IEEE Trans. Robot. 35, 124–134. doi: 10.1109/TRO.2018.2878318

CrossRef Full Text | Google Scholar

Trivedi, D., Lotfi, A., and Rahn, C. D. (2008). Geometrically exact models for soft robotic manipulators. IEEE Trans. Robot. 24, 773–780. doi: 10.1109/TRO.2008.924923

CrossRef Full Text | Google Scholar

Webster, R. J. III, and Jones, B. A. (2010). Design and kinematic modeling of constant curvature continuum robots: a review. Int. J. Robot. Res. 29, 1661–1683. doi: 10.1177/0278364910368147

CrossRef Full Text | Google Scholar

Zajac, F. E. (1989). Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Crit. Rev. Biomed. Eng. 17, 359–411.

PubMed Abstract | Google Scholar

Nomenclature

[ad_{ξ, η}] $(X)=\left(\begin{array}{cc}\tilde{\text{k}},\tilde{\text{w}}& 0_{3\times 3}\\ \tilde{\text{q}},\tilde{\text{v}}& \tilde{\text{k}},\tilde{\text{w}}\end{array}\right)\in {ℝ}^{6\times 6}$. adjoint representation of the strain twist vector in SE(3).

[$\text{a}{\text{d}}_{\xi ,\eta }^{*}$] $(X)=\left(\begin{array}{cc}\tilde{\text{k}},\tilde{\text{w}}& \tilde{\text{q}},\tilde{\text{v}}\\ 0_{3\times 3}& \tilde{\text{k}},\tilde{\text{w}}\end{array}\right)\in {ℝ}^{6\times 6}$. coadjoint representation of the strain twist vector in SE(3).

[Ad_g] $(X)=\left(\begin{array}{cc}R& 0_{3\times 3}\\ \tilde{\text{u}}R& R\end{array}\right)\in {ℝ}^{6\times 6}$. Adjoint representation of SE(3).

[$\text{A}{\text{d}}_g^{*}$] $(X)=\left(\begin{array}{cc}R& \tilde{\text{u}}R\\ 0_{3\times 3}& R\end{array}\right)\in {ℝ}^{6\times 6}$. coAdjoint representation of SE(3).

[$e^{X\hat{\xi }}$] $=I_4+X\hat{\xi }+\frac{1}{{\theta }^2}(1-cos(X\theta )){\hat{\xi }}^2+\frac{1}{{\theta }^3}(X\theta -sin(X\theta )){\hat{\xi }}^3$. Exponential map in SE(3) with θ² = k^Tk.

[T_g] $(X)=X{I}_6+\frac{1}{2{\theta }^2}(4-4cos(X\theta )-X\theta sin(X\theta ))\text{ a}{\text{d}}_{\xi }+\frac{1}{2{\theta }^3}(4X\theta -5sin(X\theta )+X\theta cos(X\theta ))\text{ a}{\text{d}}_{\xi }^2+\frac{1}{2{\theta }^4}(2-2cos(X\theta )-X\theta sin(X\theta ))\text{ a}{\text{d}}_{\xi }^3+\frac{1}{2{\theta }^5}(2X\theta -3sin(X\theta )+X\theta cos(X\theta ))\text{ a}{\text{d}}_{\xi }^4$. Analytic expression of the tangent operator of the exponential map in SE(3).