Soft continuum manipulators are flexible and highly deformable robots composed of soft and mostly elastic materials, and can serve as possible substitutes for rigid robots. Advantages of soft manipulator robots such as their compliance, dexterity, and adaptability to complex workspaces are driving the emergent research in this field. By contrast, rigidity of traditional rigid robots limits their use in constrained and confined environments, and reduces the possibilities for safe interaction with humans. Soft continuum manipulators have found applications in many areas, such as dexterous grasping (McMahan et al., 2006; Katzschmann et al., 2015) and assistive devices (Ansari et al., 2017), and particularly in the field of minimally invasive surgeries, such as laryngeal surgery (Simaan et al., 2004), catheter-based endovascular intervention (Grady et al., 2000; Burgner et al., 2013), and cardiovascular surgery (Kesner and Howe, 2011).

Analytical modeling of soft manipulators helps evaluate their motion and determine their workspace, in order to be used for control, motion planning, and animation purposes. Soft manipulators distinguish themselves by having an infinite number of degrees of freedom in any workspace they occupy. This characterization makes modeling complicated for soft manipulators. Several approaches have been investigated thus far in the literature. Most of the approaches consider the kinematic (i.e. static or quasi-static) modeling of the manipulators such as static analysis using virtual-work model (Xu and Simaan, 2010), Cosserat rod theory (Pai, 2002; Jones et al., 2009; Mahvash and Dupont, 2011), and α Lie group formulation (Grazioso et al., 2019). These models do not describe full dynamics of the manipulators, and they may show performance degradation when it comes to high-frequency applications or large and complex deformations. On the other hand, dynamical modeling approaches [e.g. Wen et al. (2012), Jung et al. (2014), Hyatt et al. (2019), Sadati et al. (2019), Till et al. (2019), Tariverdi et al. (2020)], contain dynamics of the manipulators and also take into account time-varying responses of manipulators, including high-frequency modes. However, the dynamic models mostly rely on traditional methods, such as finite elements and finite differences (i.e., quantitative and numerical methods), making the algorithms computationally expensive for real-time applications. In other words, to obtain sufficiently accurate solutions, methods need to deal with fine meshes, which increase memory use and computation time. Another limitation is that their solutions are discrete or not sufficiently differentiable. It is worth noting that in model-based controllers or observers, having a differentiable solution (i.e., a solution that can be evaluated continuously on the workspace) is crucial in the design process. Furthermore, when softer materials are employed for manipulator construction with more complex geometries or large deformations, modeling their behavior analytically becomes challenging. Therefore, there is a need for appropriate data-driven approaches without compromising computational bandwidths and the prediction quality.

Dynamics of soft continuum manipulators have highly nonlinear behavior and are expressed as Partial Differential Equations (PDEs). An effective approach to represent and model PDEs solutions is to use Neural Networks (NN). NN-based solutions of PDEs are infinitely differentiable by eliminating the need for interpolation. Furthermore, compared to finite elements or difference methods, solutions are represented by fewer parameters, which reduces the memory use. There are studies that use machine learning algorithms to find a solution for special types of PDEs such as (Lagaris et al., 1998; Lee and Kang, 1990; Weinan et al., 2017; Raissi et al., 2019). However, to the authors’ best knowledge, there is no study that investigates possible NN-based solutions for partial differential equations that describe the full dynamics of continuum manipulators. In this work, inspired by a time-space integration scheme and by using the Lie group variational integration method (Demoures et al., 2015), the dynamic equations for translation and rotation for each node of a soft continuum manipulator are decoupled, providing an appropriate structure aimed at developing a real-time modeling algorithm. Afterward, Recurrent Neural Networks (RNNs)-based models are employed to approximate the high-dimensional discretized equations. Additionally, external torques and forces (e.g., control inputs, friction, and gravity) are incorporated into the model in a real-time manner for control applications.

The ability of RNNs to learn and approximate large classes of nonlinear functions over sequences of inputs accurately makes them prime candidates for use in dynamic modeling of complex nonlinear systems. RNNs with Long Short-Term Memory (LSTM) layers process sequences by iterating through the sequence elements. Using an internal feedback, the network is capable of preserving long-term dependencies. Essentially, LSTM layers prevent older information from gradually vanishing. These networks also have been used for several applications in soft robotics. To name a few, Thuruthel in Thuruthel et al. (2019) proposes a model-free, real-time sensing method for soft robots perception. The authors in Thuruthel et al. (2017) uses RNNs to model and control soft robotic manipulators. Also, force and motion estimation using RNNs has been investigated in Marban et al. (2019) and Turan et al. (2018), respectively.

This paper aims to develop a real-time dynamic model for analyzing the dynamics of soft manipulators. Investigation of previous work on the modeling of the continuum manipulators suggests that existing literature focuses primarily on static or quasi-static approaches, or does not provide a real-time model. The contribution of this article is to present a scalable, parallel and real-time modeling algorithm for soft manipulators dynamics. The contributions of this paper are as follows.

• Existing approaches primarily deal with kinematic modeling methods. Nevertheless, in this study, real-time prediction of soft manipulators full spatial dynamics is considered in the proposed RNN-based algorithm by proposing multiple light-weight RNN-based models.

The remainder of this paper is organized as follows: the problem statement is given in Section 2. Section 3 describes the proposed RNN-based algorithm in details. In Section 4 and Section 5, different simulations and experimental validation are presented to demonstrate the efficacy of the proposed RNN-based method, in terms of the model performances and accurately predicting poses of manipulators. Finally, the main conclusions are stated in Section 7.

Consider a continuum manipulator with large deflections described by dynamic equations of motion [as presented in Tariverdi et al. (2020) and Demoures et al. (2015)] in the PDEs form as H ω t + ω × H ω + n × Λ − 1 ϕ x − Λ − 1 Λ x × m − m x = Λ − 1 τ M ϕ t t − Λ ( Λ − 1 Λ x × n ) − Λ n x + f c = f (1) where M = ρ × A (ρ and A are the manipulator constant mass density and its cross-section area), ω ∈ ℝ 3 is the manipulator’s angular velocity, H ∈ ℝ 3 × 3 is the manipulator’s inertia matrix, ϕ ∈ ℝ 3 is the position of the manipulator’s line of centroids in its workspace, Λ ∈ S O ( 3 ) denotes the orientation of moving cross-sections at point ϕ. Also, n ∈ ℝ 3 and m ∈ ℝ 3 are the stresses and momenta along the manipulator, f c ∈ ℝ 3 represents conservative forces (e.g. gravity). Furthermore, ( ⋅ ) x , ( ⋅ ) t , and ( ⋅ ) t t denote partial derivatives with respect to position, time, and the second partial derivative with respect to time, respectively. Finally, f ∈ ℝ 3 and τ ∈ ℝ 3 are non-conservative forces and torques (e.g., frictions and control inputs) ¹ .

Although high fidelity models given in the references can describe continuum manipulators dynamics efficiently, they suffer from limitations that are discussed in Section 1. Inspired by the structure and formulation of the dynamics based on the Lie group variational integration scheme, the aim is to propose distributed deep recurrent neural networks to capture and simulate soft manipulators dynamics in real-time to be able control them more accurately than existing models.

This section is devoted to develop a model based on the time series prediction using RNNs. To solve PDEs numerically using NNs, one approach is to utilize discrete solutions of finite element or difference methods to train an NN. A Lie group variational time integration model is employed to discretize the continuous dynamics of a soft manipulator ² . The whole manipulator is discretized into an arbitrary number of nodes where the position and orientation equations of each node are decoupled. In our study, we discretize the manipulator with equidistant nodes, but this can be changed depending on the application.

Figure 1 demonstrates a soft continuum manipulator at time t where x * is the undeformed length of Node n − 1 . The force F ( x ∗ , t ) , torque τ ( x ∗ , t ) are applied to Node n − 1 at the position ϕ ( x ∗ , t ) . Also, Λ ( x ∗ , t ) is the orientation matrix from the frame { O } to the frame { O t n − 1 } attached to the cross-section of Node n − 1 .

The discrete equations suggest an appropriate structure for the RNNs-based model. Given time-sequence inputs (as a first input layer), i.e., poses (positions and orientations) of nodes, and also forces and torques (as a second input layer) applied to each node, the RNN-based model of Node n is depicted in Figure 2A,B. For Node n, the first input layer is a time-sequence series of poses p n − 1 , p n , and p n + 1 (i.e., poses of Node n and its adjacent nodes n − 1 and n + 1 ) and the second input layer includes forces and torques of node n at time t, i.e., [ F t n , τ t n ] T which are incorporated into the model through dense layers.

The network takes specific size vectors as inputs, which are called input layers. The inputs are transformed through a series of hidden layers (LSTM, dense, or fully connected layers) to produce an output. The output vectors are called an output layer. Dense or fully connected layers perform linear operations (i.e., multiplication and summation) on their inputs. Furthermore, LSTM layers consist of LSTM units, which can process sequences of data of any length, for example, poses of nodes. An LSTM unit controls contributions of each element of the input layer in the output and keeps track of the dependencies between the elements (Hochreiter and Schmidhuber, 1997).

For the training process, data-sets contain time-sequence inputs and forces and torques applied to each node. Also, for each node, the poses of the node and its neighbors are considered features, as shown in Figure 2A,B. The first and second input layers proceed through LSTM layers and dense layers as hidden layers, respectively. Finally, output layers have resulted from fully connected layers.

By augmenting the given models for all nodes (see Figure 2B) as a series, the proposed RNN-based models of the whole continuum manipulator with N nodes with non-conservative forces and torques are depicted in Figure 2C. Output of every node is updated at each time step by using a history (at least two previous time steps) of neighboring node outputs. Therefore, the proposed architecture suggests a suitable framework to construct a parallel modeling algorithm.

In this section, we consider different examples and evaluate the performance of the proposed RNN-based model in Figure 2C. It is worth mentioning that data-sets play a crucial role in efficiency and accuracy in machine learning-based algorithms. The data acquisition process from a robot in real-world environments is both time and cost-consuming (implementation of multiple sensors, data filtering, and fusion, etc.). As an alternative approach, the required data can be acquired through simulations of high fidelity models. The obtained data can thus be transferred to train the algorithms to be implemented in real-world scenarios. In this section and for the presented examples, required data for the training of the proposed RNN-based model are acquired through simulations of the algorithm presented in Tariverdi et al. (2020), Sec. 2). For clarity, this model is henceforth referred to as the analytical dynamic model. In addition, since thin rods are considered in the examples, orientations of cross-sections are not of any concern. Also, it should be noted that orientations, except the twisting angle, can be reconstructed from manipulators’ configuration. Therefore, to obtain a computationally light model, the focus of attention is only on the prediction of positions.

This section is devoted to the experimental validation of the presented model. To that end, we fabricated a soft manipulator on which magnetic fields are used to produce necessary forces and torques. Compared to the simulations in which positions are predicted, time-sequence input is composed of orientations of nodes in the experiment. Furthermore, to show the performance of the algorithm, results from the presented method and a Cosserat rod-based theoretical model are compared to show the efficiency of the proposed RNN-based model. The Cosserat rod model of the soft manipulator is detailed in Appendix.

The distal and proximal node rotations are predicted both by Cosserat rod model and the proposed model, and the results are shown in Figure 12. Also, the maximum and mean absolute errors are stated in an ordered pair in Table 1.

The computation time required to find a solution of the manipulator statics from a Boundary Value Problem (BVP) with Cosserat rod theory depends on the quality of the initial solution guess, i.e., n ( s ) and m ( s ) at s = 0 , the tolerable error ( E ∈ ℝ ), and the number of nodes ( N ∈ ℕ ) used to discretize the manipulator.

A tolerable error describes the error between the distal internal forces and moments obtained from forward integration which are called n f d and m f d , and distal boundary condition, i.e., n b d and m b d . The tolerable error can be written as [ n f d − n b d , m f d − m b d ] 2 ≤ E .

Decreasing the tolerable error increases the solution accuracy, but potentially requires more time to solve convex optimizations for the BVP. Increasing the number of nodes is necessary to describe complex manipulator geometries, but should be chosen to minimize the required steps during forward integration.

To visualize how the required computation time changes with the number of nodes and the tolerable error, multiple simulations were performed by assigning known torques τ m and forces f m for m = 1,2 , to the manipulator, and finding a valid solution from solving the BVP. Changes to tolerable errors (E) and number of nodes (N) were made manually. For example, an error of 2 % in the initial solution guess was obtained by multiplying the valid solution with 0.98. After each change the BVP was solved again fifty times. The obtained mean and standard deviation of the computation times are shown in Figure 13. By taking into account all the aforementioned variables, i.e., number of nodes and the tolerable error, the Cosserat rod model is capable of achieving real-time performances of the bandwidths between 8.33 and 50 Hz.

Figure 13A shows how the computation time required for solving a solution to the BVP changes with decreasing tolerable error (E) and increasing percentage errors from a valid solution (at 0 % ), for a constant number of nodes. Also, Figure 13B shows how the computation time required for solving a solution to the BVP changes with an increasing number of nodes (N) and increasing percentage errors from a valid solution, for a constant tolerable error. However, It should be mentioned that the proposed RNN-based model shows a real-time performance with a bandwidth of 60.75 Hz on average for the given architecture in Figure 7, number of epochs = 25, and batch size = 1. In addition, Figure 14A demonstrates the computation bandwidth required for the prediction of the next step using the trained model with a different number of LSTM units and a different size of time history horizons in Data-set IV. The figure suggests that computation bandwidths are fairly unchanged with the number of LSTM units; however, increasing the length of time history reduces bandwidth. The optimal region maximizing the bandwidth is approximately with time history size in ( 2,20 ) and LSTM unit size in ( 5,15 ) . Figure 14B suggests that RMS error of the prediction decreases by increasing the number of LSTM units and the optimal area minimizing RMS errors is approximately with time history size in ( 2,20 ) and LSTM unit size in ( 20,25 ) .

To sum up, this experiment demonstrates that not only can the presented RNN-based model outperform classical modeling approaches such as the Cosserat rod model, but also it shows possibilities to use the model in practice for closed-loop control applications.

This work suggests a distributed architecture for modeling complex dynamical systems by using multiple light-weight RNN-based models. As a result, the architecture would be easier to design and debug, and also benefits from faster convergence compared to one large network. Furthermore, large networks may take longer times to be trained, and they may not show an acceptable performance and readjusting (hyper-)parameters and restarting the training process might be necessary.

Increasing the size of history horizons in training stages may reduce the error to some extent, but on the other hand, it makes the model slower. Based on conventional dynamical models, the length of the history size should be at least 2. To reach a state-of-the-art performance, i.e., having less error and faster model simultaneously, one may prefer varied batch sizes in the training and run-time phases. As a suggestion, we can use different batch sizes for training and run-time stages. A model can be trained with appropriate batch sizes such that the model performance suits the given criteria. Afterward, one can create a new network with the pre-trained weights compiled with a batch size of 1.

The performance, i.e., the convergence and stability, of the presented algorithm in this paper, unlike conventional algorithms, is independent of the number of nodes considered for the whole manipulator. To be specific, in the analytical model, there might be a need for several discretization nodes to achieve a convergent solution with a specific tolerable error; however, in the RNN-based model, only specific points/points of interests (e.g., two actuation points in the experiment) are considered. In other words, in the experiment, 13 nodes (4 for each flexible subsection and two for each magnet, and 1 for the base) were chosen for solving the Cosserat rod model, but two nodes were selected for the RNN-based model. However, the complexity of dynamical systems (i.e., PDEs) affects the complexity of the architecture used in the RNN-based model, i.e., the number of layers and LSTM units and generally how deep the model is. Nevertheless, the suggested model suits parallel implementation and can benefit from a high bandwidth for closed-loop control applications. Furthermore, the architectures of the proposed RNN-based model can be optimized by reducing the number of layers and trainable parameters to maximize the achievable bandwidths.

The evaluations showed that incorporating poses of adjacent nodes and also wrenches as a separated input might help to have, to some extent, a generalizable model rather than just purely learning the structure of data. However, supervised learning methods likely tend to preserve structure of data, and these models might not entirely respect underlying physics (conservation laws). In other words, these methods might not be wholly physics-aware and applicable for untrained/unprecedented dynamics or geometries without any adjustment, re-training, or using techniques such as transfer learning, etc. One possible and interesting solution (Psichogios and Ungar, 1992; Lagaris et al., 1998; Raissi et al., 2019) to overcome this problem and move toward fully physics-aware neural networks is revisiting lost functions for the training process. To be specific, it is mentioned in Problem Statement Section that the idea is finding solutions for PDEs given in Eq. 1, i.e., Λ ( x , t ) and ϕ ( x , t ) for sufficiently large number (e.g., N f ) of pair ( x i , t i ) ∈ ( 0 , L ) × [ 0 , T ] in which L is the unreformed length of the manipulator and parameter T is a user-defined time. Considering Eq. 1, a neural network can be learned by minimizing the mean squared error loss 1 N f ( ∑ i = 0 N f ‖ J ω t i + ω × J ω + n × Λ − 1 ϕ x i − Λ − 1 Λ x i × m − m x i − Λ − 1 τ ‖ { x i , t i } 2 + ∑ i = 0 N f ‖ M ϕ t i t i − Λ ( Λ − 1 Λ x i × n ) − Λ n x i + f n c − f ‖ { x i , t i } 2 )

This modified loss function enforces the structure imposed by Eq. 1 for large number (e.g., N f ) of pair ( x i , t i ) ∈ ( 0 , L ) × [ 0 , T ] and the trained neural network will be aware of governing PDEs.

This paper describes an approach for the real-time prediction of dynamics for general continuum soft manipulators, based on machine learning techniques and Lie group variational integration methods. Poses of a soft, polymer-based manipulator, in the presence of conservative and non-conservative wrenches, are predicted and validated experimentally. The comparison results of the proposed model and a well-known model for continuum manipulators, i.e., Cosserat rod theory, are also provided, revealing the practical effectiveness of the proposed model. The presented method can be extended to different soft robots with different shapes and materials. In addition, training of physics-aware neural networks for solving PDEs and the procedure of a model-based controller design are topics of research to be studied as future work.