First-Order Dynamic Modeling and Control of Soft Robots

Modeling of soft robots is typically performed at the static level or at a second-order fully dynamic level. Controllers developed upon these models have several advantages and disadvantages. Static controllers, based on the kinematic relations tend to be the easiest to develop, but by sacrificing accuracy, efficiency and the natural dynamics. Controllers developed using second-order dynamic models tend to be computationally expensive, but allow optimal control. Here we propose that the dynamic model of a soft robot can be reduced to first-order dynamical equation owing to their high damping and low inertial properties, as typically observed in nature, with minimal loss in accuracy. This paper investigates the validity of this assumption and the advantages it provides to the modeling and control of soft robots. Our results demonstrate that this model approximation is a powerful tool for developing closed-loop task-space dynamic controllers for soft robots by simplifying the planning and sensory feedback process with minimal effects on the controller accuracy.


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. For fluidic actuation, other models have been adopted (Marchese et al., 2016;. 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 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 .
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 . 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 .
This article investigates the viability of a first-order dynamic model for soft robots. It must be noted that unlike statespace 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 modelorder reduction provides two advantages. First, first-order dynamical systems are computationally cheaper than secondorder 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.

THEORY
Given a soft robot that can be kinematically modeled by the configuration-space q ∈ R n , the task-space variable can be obtained by the kinematic transformation: Where, x ∈ R 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: Here, M(q) represents the inertial properties, C(q,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 (whenq = 0), the first order term dominates the second order term [i.e., C(q,q)q >> M(q)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 secondorder dynamical model (Equation 2) into: 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 closedloop dynamic controller would require only the zero-order state feedback (q) for control.

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: 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: 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.
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: The operational space controller mapping now becomes: 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.

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 R 3 and its corresponding skew-symmetric matrix in so(3).

Kinematics
The relative position and orientation of a soft body i with respect to its predecessor in the chain is defined as a curve 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. 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): where ξ i (X) ∈ R 6 defines the strain state, η i (X) ∈ R 6 is the crosssection velocity twist and ad (·) ∈ R 6×6 is the adjoint operator of the Lie algebra (see Nomenclature). Going further into detail, we have with k i (X) ∈ R 3 and p i (X) ∈ R 3 the angular and linear strain; and w i (X) ∈ R 3 and v i (X) ∈ R 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: where B i ∈ R 6×n i forms a basis for the allowed motion subspace, q i ∈ R n i contains the values of the allowed strains and, ξ * i ∈ R 6 is the reference twist modeling the reference shape.
Assuming piece-wise constant strains , Equations (8) can be analytically integrated using the matrix exponential method, leading to: where Ad g i (X) ∈ R 6×6 is the Adjoint operator of SE (3), and T g i (X) = X 0 e sad ξ 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) ∈ R 6×n and its derivativeJ i (q,q, X) ∈ R 6×n (n being the total number of DOFs), which relates the generalized coordinate vector q = q T 1 q T 1 · · · q T N T ∈ R n (N being the total number of bodies) and the velocity twist η i (X), for each soft body i, as shown below.
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.
Frontiers in Robotics and AI | www.frontiersin.org where the block elements of the i th Jacobian i S (·) ∈ R 6×n (·) and its derivative iṠ ∈ R 6×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).

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 : 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);F e i (X) ∈ R 6 is the distributed external load; F a i (X) ∈ R 6 is the internal wrench due to the distributed actuation ; F i i (X) ∈ R 6 is the internal wrench due to the elasticity of the soft body; F J (·) ∈ R 6 is the wrench transmitted across joint (·) and ad * (·) (respectively Ad * (·) ) ∈ R 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: where 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: where M ∈ R n×n is the generalized mass matrix, C ∈ R n×n is the generalized Coriolis matrix, D ∈ R n×n is the block-diagonal generalized damping matrix, K ∈ R n×n is the block-diagonal generalized stiffness matrix, F ∈ R n is the vector of generalized position-dependent external forces and τ ∈ R n is the vector of applied actuation forces. Note that the dynamic Equation (13) can be written in the form required by Equation (2), with Frontiers in Robotics and AI | www.frontiersin.org C(q,q) = C q,q + D and G(q) = Kq − F q . Going further into details, the coefficient matrices take the form: 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.

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.

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 taskspace 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 secondorder model, for a fair comparison. The training parameters of the NARX network is given in Table 1.
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. 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.
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 Octopusinspired soft manipulator, which was manufactured with silicone and driven by tendons (Renda et al., 2014).
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 3 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.

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 nonredundant 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 i , x i+1 ) → τ . When testing the controller, the next task-space coordinate, x i+1 , 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.
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 0 and the two random points are p 1 and p 2 , the generated path is from p 0 → p 1 → p 2 → p 1 . 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.
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 closedloop controller. This also shows how important it is to close-theloop, even at the cost of reduced model accuracy. The tracking error of the three tests described in this scenario is shown in Figure 10.

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 nonredundant cases. Moreover, the sensory requirements for closedloop 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.