# Global Control of Soft Manipulator by Unifying Cosserat and Neural Network

Haihong Li, Lingxiao Xun, Gang Zheng

► To cite this version:

Haihong Li, Lingxiao Xun, Gang Zheng. Global Control of Soft Manipulator by Unifying Cosserat and Neural Network. IEEE Transactions on Industrial Electronics, 2024, 71 (9), pp.10944-10954. 10.1109/TIE.2023.3331109. hal-04325359

HAL Id: hal-04325359

https://inria.hal.science/hal-04325359v1

Submitted on 5 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Global Control of Soft Manipulator by Unifying Cosserat and Neural Network

Haihong Li, Lingxiao Xun, and Gang Zheng

***Abstract***—Soft manipulators, a relatively novel class of robots, are increasingly prevalent in the field of robotics. Due to the material nonlinearity and the absence of links and joints which are typically used to derive their kinematics and dynamics, the exact modeling and effective control frameworks of such soft robots are much more difficult than the traditional rigid counterparts. To achieve global control of the end-effector position of the soft manipulator actuated by cables, this work proposes a hybrid control strategy by unifying the piecewise linear strain (PLS) Cosserat model and radial basis function neural network (RBFNN) to approximate the Jacobian matrix of any end-effector position in the whole workspace, and then designs a global control scheme based on the approximation of Jacobian matrix. The theoretical convergence proof and experimental validation are provided for the designed global controller. The results corroborate that the proposed control strategy has excellent accuracy and robustness in tracking desired time-varying trajectories through different experiments.

***Index Terms***—Soft manipulator, Cosserat, neural network, controller.

## I. INTRODUCTION

The biological systems in nature provide significant inspiration of soft structures which are able to bend, extend, shear and twist. Examples include elephant trunk, octopus tentacles, mammalian tongue, robotic fish, etc. Unlike the conventional rigid robots, all of such continually deformable structures have an infinite number of degrees of freedom, which provides soft robots an incredible agility to adapt to new environments and move in narrow spaces [1]. However, it remains a substantial challenge to establish a mathematical model suitable to accurately describe the kinematics and dynamics of such high-dimensional systems due to their inherent properties of compliance and flexibility. Additionally, the development of high-performance control approaches for soft robots is also a challenging research area.

To address the above-mentioned challenges, one of the most used techniques is to deduce the exact models of soft robots. In [2], the kinematic model was derived using geometric information, followed by the application of a computed torque controller to control an eel-like soft robot. The piecewise constant curvature (PCC) has proven to be quite feasible with a vast range of control applications [3], [4]. For example, their use in the inverse kinematic control [5], [6], the closed-loop dynamic feedback controllers [7], [8], [9], and feedforward dynamic control [10]. However, under the condition of external loads or contact with the environment, the validity of the CC assumption is relatively conservative and unsuitable for many

applications [11]. Other techniques have also been used to get the models. In [12], the authors employed Euler-Bernoulli beam theory to model the bending mode of the robot, and developed the position and orientation control for an eversion robot with changing structural stiffness. FEM models are suitable to handle material and geometric nonlinearities, and can provide a more accurate representation of deformations of soft robots [13], [14]. In order to provide a general method to control soft robot, an open-loop controller was designed via the FEM model [13]. In [15], the authors employed FEM and visual tracking technique to achieve the feedback control for soft robots. [16] analyzed the dynamics of soft robots, and used the reduced-order model technique to achieve fast convergence of soft robots to a desired position.

The Cosserat beam model is the geometrically nonlinear generalization of the other beam theory-based modeling methods for slender soft robots [17]. [18] overcame the expensive model computational obstacle, and adopted the Cosserat rod model to design a robust controller for continuum robots. The authors of [19] combined the Cosserat kinematic model with integrated sensing and the Jacobian matrix to provide a practically feasible closed-loop control system for application to soft robots. Previous research has also explored the use of discrete Cosserat model-based controllers. For instance, in [20], a first-order approximation assumption of piecewise constant strain (PCS) dynamics was proposed, based on which the closed-loop task-space dynamic controllers were developed. The geometric variable strain Cosserat approach [21] was used to carry out the Jacobian-based inverse kinematic control of the general soft manipulator. However, this model-based control method ignored the external forces and lacked the experimental validation [22]. The piecewise linear strain (PLS) Cosserat model proposed in [23] was used to achieve the local control of end-effector position for the cable-driven soft manipulator [24]. For those mentioned Cosserat-based modeling method where the strains were used to parameterize the manipulator's deformation, the proposed global controllers in the existing literature have only been validated through simulation, but cannot be implemented into the real prototype, as those strains are difficult to be measured in practice.

All of the aforementioned control techniques, for the tracking control of soft continuum robots, require the exact models (either kinematic, or dynamic) to be known. In contrast, model-free controllers have proposed to control the soft robots with the uncertainties and complexities of nonlinear behaviors using machine learning techniques or empirical methods. In [25], the machine learning-based method was introduced to achieve the closed loop kinematic control of
continuum robots in the task space. [26] used the kinematic and kinetic models based on artificial neural network (ANN) to control the compact bionic handling arm. For the work of [27], the authors employed the continual learning technique to design controllers able to continuously adapt to changes in robot dynamics. [28] used an ANN to approximate the input-output relation of a cable-driven soft robot, designed the robust controller, and provided the proof of the stability and robustness as well. [29] presented a machine learning based approach towards development of dynamic model and a trajectory optimization method of predictive control for the soft manipulator. In [30], a data-based control framework was introduced to solve the soft robot underwater locomotion problem using deep reinforcement learning. [31] described a Koopman-based approach and its application to model predictive control design for pneumatic soft robots. [32] employed the machine learning approach to estimate a differentiable model of the quasi-static physics of a soft robot, and then performed gradient based optimization to find the best open-loop control input. [33] presented a method to train the forward kinematic model and its Jacobian together as two neural networks to realize the real-time computation of inverse kinematics on soft robots. Despite their gratifying advances, collecting a representative and sufficient amount of data is still challenging for most learning-based control approaches for soft robots [34]. Furthermore, stability analysis and theoretical convergence proofs are difficult to provide for the learning-based controllers.

The hybrid control method combining model-based and learning-based control approaches has also been a cutting-edge control strategy for soft robots, which can be regarded as an improvement of model-free control while incorporating the strengths of physics models. In [35], a supervised learning method to solve the inverse statics of the cable-actuated soft manipulator with non-constant curvature was introduced. [36] presented a model-based policy learning algorithm for closed-loop predictive control of a soft robotic manipulator. For the work of [37], the authors put forward a hybrid modeling method which combined the first-principles model with machine learning methods to improve overall performance of the nonlinear model predictive controller. Based on FEM and ANN, [38] proposed a transfer learning scheme to minimize the effort of generating real data for neural network training, and achieved the fine control of a real soft pneumatic actuator. In [39], the authors introduced a model-based online learning and adaptive control algorithm for the wearable soft robotic glove, which enabled the soft glove to adapt to diverse hand conditions for reference tracking.

To the best of our knowledge, up to now, the discrete Cosserat static model-based control technique combined with the neural networks is scarcely explored in soft robotics. This motivates the work of this paper, where we propose a hybrid approach to globally control the end-effector position of the cable-driven soft manipulator by unifying both the PLS Cosserat static model and the neural networks to track the trajectories within its whole workspace in real time. In summary, compared to other kinematic or static model-based controller design for soft robots, the contributions of this work

<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic of the end-effector position control for the cable-driven soft manipulator.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Sub-workspaces and their associated configurations for the soft manipulator.</figcaption></figure>

include the following:

* An optimization-based approach is presented to obtain the whole workspace of the soft manipulator via PLS Cosserat, and a novel decomposition strategy of the workspace is proposed;

* A method to solve the inverse statics of the soft manipulator modeled by the PLS Cosserat approach using the minimum potential energy principle is implemented to determine the unique configuration of manipulator;

* The radial basis function neural network (RBFNN) is unified with PLS Cosserat to globally approximate the Jacobian matrix of any end-effector position of the soft manipulator in its whole workspace;

* The robust control scheme for the cable-driven soft manipulator is established with complete mathematical proof. Experimental implementations on the soft prototype corroborate the feasibility and robustness of the proposed controller.

## II. PROBLEM STATEMENT

In this paper, the objective is to develop an efficient closed-loop controller based on the approximated static model for the purpose of controlling the end-effector position of the cable-driven soft manipulator within its whole workspace $W_E$.

To realize this objective, PLS Cosserat rod model [23] is adopted in our work to model the soft manipulator, as depicted in Fig. 1. The material-attached reference frame
(d₁, d₂, d₃) of the continuous Cosserat rod can be defined as $g(s) = (\mathbf{R}(s) \quad \mathbf{u}(s)) \in SE(3)$ where $\mathbf{u}(s) \in \mathbb{R}^3$ represents the position field of any cross section at $s \in [0, L]$, $\mathbf{R}(s) \in SO(3)$ stands for the material orientation in the form of an orthonormal rotation matrix w.r.t. the inertial frame $(p₁, p₂, p₃)$. In the framework of PLS, the configuration space, noted as $q \in \mathbb{R}^{6(N+1)}$, is parameterized by 6 strain components of each interpolation node of soft manipulator divided into N sections. We denote $\mathbf{u} \in \mathbb{R}^3$ as the end-effector position of the soft manipulator (workspace), and $\mathbf{T} \in \mathbb{R}^N$ as tension magnitude vector of $\overline{N}$ ($\overline{N} \ge 4$) cables (actuator space) attached at the end-effector of the manipulator, then by applying the principle of virtual work as described in [21], the PLS Cosserat static and kinematic models are derived as

$$K\mathbf{q} + \mathbf{G}(\mathbf{q}) = \mathbf{H}\mathbf{T} \quad (1a)$$

$$\mathbf{u}(\mathbf{q}) = \mathbf{P}\mathbf{g}(\mathbf{q})\mathbf{E} \quad (1b)$$

$$\boldsymbol{\eta} = \mathbf{J}(\mathbf{q})\dot{\mathbf{q}} \quad (1c)$$

where $\mathbf{K} \in \mathbb{R}^{6(N+1) \times 6(N+1)}$ is the generalized stiffness matrix, $\mathbf{G}(\mathbf{q}) \in \mathbb{R}^{6(N+1) \times 1}$ is the generalized gravitational matrix, and $\mathbf{H} \in \mathbb{R}^{6(N+1) \times \overline{N}}$ is the generalized actuation matrix; $\mathbf{P}$ and $\mathbf{E}^\top$ are the pre-defined selection matrices; $\boldsymbol{\eta} = (\boldsymbol{\Omega}^\top, \mathbf{V}^\top)^\top \in \mathbb{R}^6$ represents velocity twist where $\boldsymbol{\Omega} \in \mathbb{R}^3$ and $\mathbf{V} \in \mathbb{R}^3$ stand for the angular and linear velocity, respectively; $\mathbf{J}(\mathbf{q})$ is the geometric Jacobian matrix describing the relation between the velocity twist and the time derivative of the generalized strain of the soft manipulator. With the above PLS Cosserat models, we can obtain

$$\dot{\mathbf{u}} = \mathbf{\Gamma}(\mathbf{q})\dot{\mathbf{T}} \quad (2)$$

where $\mathbf{\Gamma}(\mathbf{q}) \in \mathbb{R}^{3 \times \overline{N}}$ is the Jacobian matrix, which determines the directional relation between the input and the output of the system. Then the classical pseudo-inverse method, based on the inverse of $\mathbf{\Gamma}(\mathbf{q})$, can be applied to design robust controller. However, it is evident that the calculation of $\mathbf{\Gamma}(\mathbf{q})$ requires the knowledge of the generalized strain vector $\mathbf{q}$, which is practically difficult to be measured.

In order to solve the mentioned issue of the PLS Cosserat model-based control design for soft robots, this work aims to globally (i.e., for all admissible $\mathbf{q}$) approximate $\mathbf{\Gamma}(\mathbf{q})$ in advance, and thus design a robust pseudo-inverse controller, based only on the measurement of end-effector position, to realize the global control of $\mathbf{u}$.

Concretely, we propose to divide the whole workspace $W_E$ into several sub-workspaces, denoted as $W_i \subseteq W_E$, $i = 1, 2, \dots, m$, where $m$ represents the number of the sub-workspaces, as shown in Fig. 2, choose a representative equilibrium point $\mathbf{u}_i \in W_i$ from the $i$th sub-workspace, and calculate the associated Jacobian matrix $\mathbf{\Gamma}_i$ in advance. Subsequently, the global approximation of $\mathbf{\Gamma}(\mathbf{q})$ will be realized via a neural network (with parameter $\alpha_i$ to be trained) as

$$\mathbf{\Gamma}(\mathbf{q}) = \hat{\mathbf{\Gamma}}(\mathbf{u}) = \sum_{i=1}^{m} \alpha_i(\mathbf{u}_i)\mathbf{\Gamma}_i$$

where $\mathbf{\Gamma}_i$ represents the $i$th representative Jacobian matrix.

In summary, to design the presented global controller by unifying PLS Cosserat static-model and neural network, the following 4 questions need to be answered:

Q1: How to estimate and divide the whole workspace $W_E$?

Q2: For each sub-workspace, how to calculate the associated $\mathbf{\Gamma}_i$?

Q3: With the pre-calculated $\mathbf{\Gamma}_i$, how to construct the neural network to approximate $\mathbf{\Gamma}(\mathbf{q})$?

Q4: With the obtained Jacobian matrix approximation $\hat{\mathbf{\Gamma}}$, how to design a robust controller to realize the global control objective of $\mathbf{u}$?

In the following, III-A will present the detailed solution to Q1, followed by Q2 addressed in III-B. After that, III-C will correspond to the solution to Q3. Finally, III-D will provide the answer to Q4.

### III. ROBUST CONTROLLER DESIGN VIA THE PLS COSSERAT STATIC MODEL AND NEURAL NETWORK

In this section, we firstly present a workspace estimation and decomposition approach for the PLS Cosserat static model, which provides the benefit for the soft robotic applications mainly related to their control. Next, the analytical solution of the matrix $\mathbf{\Gamma}(\mathbf{q})$ is provided, and an optimization approach is presented to solve the inverse PLS Cosserat static model for the purpose of computing $\mathbf{\Gamma}_i$. After that, a novel neural network based strategy for the off-line estimation of the Jacobian-based relation matrix between the input and output of the soft manipulator in the whole workspace is proposed. Finally, we design a closed-loop feedback control, and provide the theoretical convergence proof for the designed controller.

#### A. Workspace Estimation and Decomposition

The subject of workspace determination has been widely studied by the soft robotic community, which includes all its equilibrium points where the end-effector of soft robots can reach. The goal of this subsection is to estimate and decompose the feasible workspace of the soft manipulator based on the PLS Cosserat model by using the optimization-based approach. The workspace $W_E$ of the robot's end-effector can be defined as

$$W_E = \{ \mathbf{u} \in \mathbb{R}^3 | \mathbf{u} = \mathbf{P}\mathbf{g}(\mathbf{q})\mathbf{E}, \mathbf{K}\mathbf{q} + \mathbf{G}(\mathbf{q}) = \mathbf{H}\mathbf{T}, \forall \mathbf{T} \in \mathcal{T} \}$$

where $\mathbf{T} \in \mathbb{R}^N$ represents tension magnitude vector of $\overline{N}$ cables, and $\mathcal{T} = \mathcal{T}_1 \times \dots \times \mathcal{T}_{\overline{N}}$ with $\mathcal{T}_i = [T_{\min}^i, T_{\max}^i]$ being the minimal and maximal tension bounds of the $i$th cable for $i = 1, \dots, \overline{N}$.

To obtain the correct estimation of the whole workspace, an optimization-based method [40] is implemented on the soft manipulator that has been modeled by the PLS Cosserat approach. As a result, the whole workspace estimation algorithm can be formulated as

$$\begin{array}{l}
\underset{\mathbf{q}^{(i)}, \mathbf{T}^{(i)}}{\arg \min} \quad \|\mathbf{u}(\mathbf{q}^{(i)}) - \mathcal{E}^{(i)}\|_2^2 \\
\text{s.t.} \quad \begin{cases}
\mathbf{K}\mathbf{q}^{(i)} + \mathbf{G}(\mathbf{q}^{(i)}) - \mathbf{H}\mathbf{T}^{(i)} = \mathbf{0} \\
\mathbf{T}^{(i)} \in \mathcal{T}
\end{cases}
\end{array} \quad (3)$$
with

$$
\begin{align*}
\boldsymbol{u}(\boldsymbol{q}^{(i)}) &= \mathbf{P} \boldsymbol{g}(\boldsymbol{q}^{(i)}) \mathbf{E} \\
\boldsymbol{\mathcal{E}}^{(i)} &= \boldsymbol{\delta}_x^{(i)} \boldsymbol{\delta}_y^{(i)} \boldsymbol{\delta}_z^{(i)} \boldsymbol{\mathcal{E}}_0
\end{align*}
$$

where $\mathcal{E}^{(i)}$ represents the selected radiating vector outside the workspace $W_E$, $\mathcal{E}_0$ is the initial radiating vector, $\delta_{x,i,y,z}^{(i)}$ stand for the basic rotation matrices about the x-, y-, and z-axis, respectively. The detailed algorithm procedures of the optimization-based method along with the scheme to address the non-convexity problem (i.e., how to accurately determine the interior boundary of the workspace) can be found in [40]. Eventually, we can estimate the end-effector workspace for the studied soft manipulator by solving the minimization problem (3).

The properties of Jacobian matrices of all end-effector positions in a certain sub-workspace are considered to be the same, and thus we can use a nominal Jacobian matrix to replace all Jacobian matrices in this sub-workspace. With the obtained whole workspace $W_E$, the scheme of workspace decomposition will be introduced. In terms of two different points $\mathbf{u}_i \in W_E$ and $\mathbf{u}_{i+l} \in W_E$, for $l = 1, 2, \dots, S$ (any point around $\mathbf{u}_i$), the Jacobian matrices (i.e., $\mathbf{\Gamma}_i$, and $\mathbf{\Gamma}_{i+l}$) respectively corresponding to the end-effector positions (i.e., $\mathbf{u}_i$, and $\mathbf{u}_{i+l}$) can be obtained. Thus, we propose the following decomposition rule along the workspace surface to determine whether $\mathbf{u}_i$ and $\mathbf{u}_{i+l}$ belong to the same sub-workspace or not:

$$
\mathrm{Sgn}(\mathbf{\Gamma}_i) = \mathrm{Sgn}(\mathbf{\Gamma}_{i+l}) \quad (4)
$$

where the symbol Sgn(·) applies the conventional sign func-
tion to the matrix via an element-wise way. The sign for the
elements of the Jacobian matrix Γ represents the directional
relation between cables' tension and end-effector position. If
the rule formulated in (4) holds, then uᵢ and uᵢ₊₁ will be
classified into one sub-workspace Wᵢ where uᵢ is defined
as the representative equilibrium point of the sub-workspace.
Otherwise, the two distinct sub-workspaces will be generated
for uᵢ and uᵢ₊₁, respectively. Comparing the elements' sign
from the Jacobian matrices around all equilibrium points
in Wᵢ iteratively, we can finally decompose the whole
workspace of the soft manipulator via the PLS Cosserat into
m different sub-workspaces.

B. Calculation of $\Gamma_i$ for the ith Sub-Workspace

To compute $\Gamma_i$, we need to know the analytical formula of $\Gamma(\boldsymbol{q})$ and the associated strain value of $\boldsymbol{q}$ for a given equilibrium end-effector position $\boldsymbol{u}_i \in W_i$ which in fact is the well-known problem of inverse static model. Therefore, this subsection will firstly deduce the analytical formula of $\Gamma_i$, and then present how to solve the inverse PLS Cosserat static model.

1) Analytical formula of $\Gamma(\boldsymbol{q})$: With the exact relation described by (2), the matrix $\Gamma(\boldsymbol{q})$ can be formulated as

$$
\Gamma(\boldsymbol{q}) = \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{T}} = \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{q}} \frac{\partial \boldsymbol{q}}{\partial \boldsymbol{T}} \quad (5)
$$

Obviously, the calculation of the matrix $\Gamma(\boldsymbol{q})$ can be re-
garded as the multiplication of $\frac{\partial \boldsymbol{u}}{\partial \boldsymbol{q}}$ and $\frac{\partial \boldsymbol{q}}{\partial \boldsymbol{T}}$. The exhaustive

derivation of the two parts can be found in Appendix A. From
the kinematic model (1b-1c), we have

$$
\frac{\partial \boldsymbol{u}}{\partial \boldsymbol{q}} = [\mathbf{0}_3 \quad \mathbf{R}(\boldsymbol{q})] \mathbf{J}(\boldsymbol{q})
$$

Taking the derivative of (1a) with respect to T yields

$$
\frac{\partial \boldsymbol{q}}{\partial \boldsymbol{T}} = \left( \boldsymbol{K} + \frac{\partial [\boldsymbol{G}(\boldsymbol{q})]}{\partial \boldsymbol{q}} \right)^{\dagger} \boldsymbol{H}
$$

Consequently, the analytical solution of $\Gamma(\boldsymbol{q})$ is equivalent
to

$$
\Gamma(\boldsymbol{q}) = [\mathbf{0}_3 \quad \mathbf{R}(\boldsymbol{q})] \mathbf{J}(\boldsymbol{q}) \left( \boldsymbol{K} + \frac{\partial [\boldsymbol{G}(\boldsymbol{q})]}{\partial \boldsymbol{q}} \right)^{\dagger} \boldsymbol{H} \quad (6)
$$

where $(\cdot)^{\dagger}$ implies the pseudo-inverse of $(\cdot)$.

2) Inverse PLS Cosserat Static Model: With the analytical formula of $\Gamma(\boldsymbol{q})$ deduced in Section III-B1, for a given $\mathbf{u}_i \in W_i$, if we can solve the inverse statics of (1a-1b) to get the associated $\boldsymbol{q}_i$, then we can obtain the value of $\mathbf{\Gamma}_i$ for $W_i$. Therefore, the following presents how to solve this inverse problem within an optimization framework.

For a given end-effector position $\mathbf{u}_d$, the first intuitive approach to solve the inverse statics is to minimize the following cost function

$$
\begin{equation}
\begin{split}
& \arg \min_{(\boldsymbol{q}, \boldsymbol{T})} \|\boldsymbol{u}(\boldsymbol{q}) - \boldsymbol{u}_d\|_2^2 \\
& \text{s.t.} \quad \begin{cases} \boldsymbol{K}\boldsymbol{q} + \boldsymbol{G}(\boldsymbol{q}) = \boldsymbol{H}\boldsymbol{T} \\ \boldsymbol{T} \in \mathcal{T} \end{cases}
\end{split}
\tag{7}
\end{equation}
$$

with $\mathbf{u}(\boldsymbol{q}) = \mathbf{P}\mathbf{g}(\boldsymbol{q})\mathbf{E}$, where $\| \cdot \|$ denotes the standard Euclidean norm. Theoretically, multiple solutions might exist for the above optimization problem, which heavily depends on the initial guess values of $\boldsymbol{q}$ and $\boldsymbol{T}$.

In order to overcome the issue of multiple solutions, this
paper proposes to integrate the minimum potential energy
principle into the minimal distance based cost function, since
such a principle guarantees a unique configuration of any me-
chanical system in a static equilibrium. Motivated by this idea,
for the given end-effector position $\mathbf{u}_d$ with its associated strain
configuration $\boldsymbol{q}$, the corresponding elastic potential energy $\mathcal{P}_E$
and gravitational potential energy $\mathcal{P}_G$ of total length $L$ of the
soft manipulator can be written as

$$
\mathcal{P}_E = \frac{1}{2} \int_0^L \boldsymbol{q}^\top \boldsymbol{\Phi}^\top \mathcal{H}(s) \boldsymbol{\Phi} \boldsymbol{q} \mathrm{d}s
$$

and

$$
\mathcal{P}_G = \int_0^L \mathcal{O}^\top(s) \mathcal{M}(s) \mathcal{G} \mathrm{d}s
$$

where $\mathbf{\Phi} \in \mathbb{R}^{6 \times 6(N+1)}$ is a generalized matrix consisted of coefficient matrices of strain interpolation nodes via the PLS assumption, $\mathcal{H}(s) \in \mathbb{R}^{6 \times 6}$ stands for the stiffness matrix of each cross section, $\mathcal{O}(s) = [\mathbf{0}_{3 \times 1} \mathbf{u}^{\top}(s)]^{\top}$ represents the position vector of each cross section with respect to the inertial frame, $\mathcal{M}(s) \in \mathbb{R}^{6 \times 6}$ is the mass matrix of any cross section, and $\mathcal{G} \in \mathbb{R}^6$ is the gravitational acceleration twist.
Consequently, the inverse problem is solved in this paper by minimizing the following cost function

$$ \begin{array}{ll} \underset{(\boldsymbol{q}, \boldsymbol{T})}{\operatorname{arg\,min}} & \|\boldsymbol{u}(\boldsymbol{q}) - \boldsymbol{u}_d\|_2^2 + \beta(\mathcal{P}_E + \mathcal{P}_G) \\ \text{s.t.} & \left\{ \begin{array}{l} \boldsymbol{K}\boldsymbol{q} + \boldsymbol{G}(\boldsymbol{q}) = \boldsymbol{H}\boldsymbol{T} \\ \boldsymbol{T} \in \boldsymbol{T} \end{array} \right. \end{array} \quad (8) $$

with $\beta$ being the positive weight.

To handle the inequality constraint in (8), the nonlinear interior point method [41] is used. Specifically, adding the logarithmic barrier function $\mathcal{B}$ to the objective function in order to remove the inequality constraints, the original NLP problem (8) can be then approximately re-formulated as

$$ \begin{array}{ll} \underset{\bar{\boldsymbol{x}}=(\boldsymbol{q}, \boldsymbol{T}, \boldsymbol{\mu})}{\operatorname{arg\,min}} & \|\boldsymbol{u}(\boldsymbol{q}) - \boldsymbol{u}_d\|_2^2 + \beta(\mathcal{P}_E + \mathcal{P}_G) - \boldsymbol{\mu}\mathcal{B} \\ \text{s.t.} & \boldsymbol{K}\boldsymbol{q} + \boldsymbol{G}(\boldsymbol{q}) = \boldsymbol{H}\boldsymbol{T} \end{array} \quad (9) $$

with

$$ \boldsymbol{\mu} = [\mu_1 \quad \dots \quad \mu_{2\bar{\mathbf{N}}}], \quad \mathcal{B} = \begin{bmatrix} \log(T^1 - T_{\min}^1) \\ \vdots \\ \log(T^{\bar{\mathbf{N}}} - T_{\min}^{\bar{\mathbf{N}}}) \\ \log(T_{\max}^1 - T^1) \\ \vdots \\ \log(T_{\max}^{\bar{\mathbf{N}}} - T^{\bar{\mathbf{N}}}) \end{bmatrix}, $$

where $\boldsymbol{\mu}$ is the barrier parameter vector. The Newton-type method is then used to find the optimal solution $\bar{\boldsymbol{x}}^*$ of such an optimization problem with the nonlinear equality constraints.

### C. Global Estimation of $\mathbf{\Gamma}(\boldsymbol{q})$ via Neural Network

Inspired by the thought of Section III-B, a series of Jacobian matrices corresponding to the representative end-effector positions of the $m$ sub-workspaces are computed to construct a set of matrix basis, represented by $\bar{\mathbf{\Gamma}} = [\mathbf{\Gamma}_1^\top \ \mathbf{\Gamma}_2^\top \ \dots \ \mathbf{\Gamma}_m^\top]^\top \in \mathbb{R}^{3m \times \bar{\mathbf{N}}}$, and the weight vector $\alpha(\boldsymbol{u}) = [\alpha_1(\boldsymbol{u}) \ \alpha_2(\boldsymbol{u}) \ \dots \ \alpha_m(\boldsymbol{u})] \in \mathbb{R}^{1 \times m}$ is then defined. Finally, we propose to use the matrix basis $\bar{\mathbf{\Gamma}}$ obtained from the PLS Cosserat model to globally approximate $\mathbf{\Gamma}(\boldsymbol{q})$ corresponding to the any end-effector position $\boldsymbol{u}$ while its corresponding weight vector $\alpha(\boldsymbol{u})$ will be determined via a trained neural network. In other words, for any end-effector position $\boldsymbol{u} \in W_E$ in the whole workspace, its associated matrix $\mathbf{\Gamma}(\boldsymbol{q})$ approximation problem can be then transformed to a regression one, and it yields

$$ \hat{\mathbf{\Gamma}}(\boldsymbol{u}) = \sum_{i=1}^{m} \alpha_i(\boldsymbol{u}) \mathbf{\Gamma}_i $$

To obtain the value of $\alpha(\boldsymbol{u})$, we use the radial basis function neural network (RBFNN) which is a commonly used three-layer feedforward network. Next, I will present how to approximate the parameter $\alpha(\boldsymbol{u})$. According to the input-output relation (2), the regression problem can be then formulated as the following linear optimization one:

$$ \underset{\boldsymbol{W}^{\mathrm{T}}}{\operatorname{arg\,min}} \quad f = \sum_{j=1}^{k} \|\boldsymbol{u}_j - \boldsymbol{u}_{j-1} - (\sum_{i=1}^{m} \alpha_i(\boldsymbol{u}_j)\boldsymbol{\Gamma}_i)(\boldsymbol{T}_j - \boldsymbol{T}_{j-1})\|_2^2 \quad (10) $$

with

$$ \begin{array}{l} \alpha^{\mathrm{T}}(\boldsymbol{u}_j) = \boldsymbol{W}^{\mathrm{T}} \bar{\boldsymbol{h}}(\boldsymbol{u}_j) \\ \bar{\boldsymbol{h}}(\boldsymbol{u}_j) = [h_1(\boldsymbol{u}_j) \ h_2(\boldsymbol{u}_j) \ \dots \ h_n(\boldsymbol{u}_j)]^{\mathrm{T}} \\ h_l(\boldsymbol{u}_j) = \exp\left(-\frac{\|\boldsymbol{u}_j - \boldsymbol{c}_l\|^2}{2b_l^2}\right), \quad l = 1, 2, \dots, n. \end{array} $$

where $\boldsymbol{u}_j \in \mathbb{R}^3$ is the end-effector position vector, $k$ represents the number of the equilibrium points evenly selected from the workspace, $\boldsymbol{W}^{\mathrm{T}} \in \mathbb{R}^{m \times n}$ is the estimated weight matrix, $n$ is the number of the neurons, and $\bar{\boldsymbol{h}}(\boldsymbol{u}_j) \in \mathbb{R}^n$ with $h_l(\boldsymbol{u}_j)$ being the Gaussian kernel function of the $l$th neuron, $\boldsymbol{c}_l$ and $b_l$ are the center and width of the Gaussian kernel function, respectively.

From the objective function of (10), we try to find the pairs of $(\boldsymbol{T}_j, \boldsymbol{T}_{j-1})$ and $(\boldsymbol{u}_j, \boldsymbol{u}_{j-1})$. To generate the training data sets, we start from an equilibrium point $j$ and change a little cable's tension from $\boldsymbol{T}_j$ to $\boldsymbol{T}_{j-1}$, and then the end-effector position will be changed from $\boldsymbol{u}_j$ to $\boldsymbol{u}_{j-1}$. In this way, we can obtain a set of training data. We can select $k$ points and repeat this operation to generate multiple training data sets. Finally, the gradient descent method is used to solve the above optimization problem.

**Remark 1:** It must be pointed out that one possible solution is to approximate the Jacobian matrix $\mathbf{\Gamma}(\boldsymbol{q})$ by using for example back propagation neural network (BPNN) as presented in [42]. Nevertheless, BPNN is susceptible to getting stuck in local minima during the optimization process, and more prone to over-fitting, especially when the data is scarce. In our work, we use a hybrid strategy of both the PLS Cosserat model and RBFNN to calculate the approximation of $\mathbf{\Gamma}(\boldsymbol{q})$, which effectively ensures accurate control since it takes full advantage of the model information. On the one hand, the offline calculation method of the Jacobian matrix $\mathbf{\Gamma}(\boldsymbol{q})$ by using the RBFNN greatly improves control accuracy, as RBFNN can avoid local minima, over-fitting, and better approximate the real Jacobian matrix. On the other hand, compared to the FEM-based gain-scheduling control of a soft trunk robot [43], this method globally approximates $\mathbf{\Gamma}(\boldsymbol{q})$ of any end-effector position $\boldsymbol{u} \in W_E$, and allows the soft manipulator to track different time-varying trajectories in real time.

### D. Global Robust Controller Design

Specifically, according to the deduced input-output relationship in (2), a robust controller (valid for the whole workspace) for the soft manipulator is developed in this subsection. The feedback control law is designed as follows

$$ \dot{\boldsymbol{T}} = \hat{\boldsymbol{\Gamma}}^{\dagger}(\boldsymbol{u})(-\boldsymbol{\lambda}_1 e - \boldsymbol{\lambda}_2 \int_0^t \mathrm{eds}) \quad (11) $$

where $\boldsymbol{\lambda}_1 \in \mathbb{R}^{3 \times 3}$ and $\boldsymbol{\lambda}_2 \in \mathbb{R}^{3 \times 3}$ are constant diagonal, positive-definite matrices of feedback gains; $e = \boldsymbol{u} - \boldsymbol{u}_r$ with $\boldsymbol{u}_r$ being a constant desired end-effector position. With the introduced notations, (2) can be written as

$$ \dot{\boldsymbol{e}} = -(\boldsymbol{\lambda}_1 e + \boldsymbol{\lambda}_2 \int_0^t \mathrm{eds}) - \Delta \mathbf{\Gamma} \hat{\boldsymbol{\Gamma}}^{\dagger}(\boldsymbol{\lambda}_1 e + \boldsymbol{\lambda}_2 \int_0^t \mathrm{eds}) \quad (12) $$
with $\Gamma = \hat{\Gamma} + \Delta\Gamma$, where $\Delta\Gamma$ represents the estimation error of the Jacobian matrix.

**Assumption 1:** With enough training samples to train the RBFNN described in Section III-C, it is assumed that there exists a constant $\varepsilon > 0$ such that

$$ \|\boldsymbol{\Upsilon}\| \leq \varepsilon < 1, \forall T \in \mathcal{T} \qquad (13) $$

with $\boldsymbol{\Upsilon} = \begin{bmatrix} \Delta\Gamma\hat{\Gamma}^{\dagger} & 0_3 \\ 0_3 & \Delta\Gamma\hat{\Gamma}^{\dagger} \end{bmatrix} \in \mathbb{R}^{6 \times 6}$, where $\|\boldsymbol{\Upsilon}\|$ represents the spectral norm of matrix $\boldsymbol{\Upsilon}$, and $\varepsilon$ is a prescribed approximation accuracy error. It is worth noting that the smaller $\varepsilon$ is, the more accurate Jacobian matrix we have obtained by using RBFNN.

**Theorem 1:** In terms of the studied soft manipulator described by (2), for given $\lambda_1, \lambda_2$, there exist symmetric positive-definite matrices $\mathbf{P}$ and $\mathbf{Q}$ that fulfill the Lyapunov equation

$$ \mathbf{P}\mathbf{A} + \mathbf{A}^{\top}\mathbf{P} = -2\mathbf{Q}, \qquad (14) $$

with $\mathbf{A} = \begin{bmatrix} 0_3 & \mathbf{I}_3 \\ -\lambda_2 & -\lambda_1 \end{bmatrix}$, if the approximation accuracy error $\varepsilon$ satisfies

$$ \varepsilon \leq \frac{\lambda_{\min}(\mathbf{Q})}{\lambda_{\max}(\mathbf{P})\|\mathbf{B}\|}, \qquad (15) $$

with $\mathbf{B} = \begin{bmatrix} 0_3 & 0_3 \\ -\lambda_2 & -\lambda_1 \end{bmatrix}$, the proposed controller (11) can then drive any $\mathbf{u} \in W_E$ to the desired constant position $\mathbf{u}_r \in W_E$, i.e.,

$$ \lim_{t \to \infty} \| \mathbf{u}(t) - \mathbf{u}_r(t) \|_2 = 0 $$

**Proof 1:** To prove the stability of the system, we introduce the state variable: $\zeta = \begin{bmatrix} \int_0^t e ds \\ e \end{bmatrix}$. Then, the system (12) can be re-written as the following state-space form:

$$ \dot{\zeta} = \mathbf{A}\zeta + \mathbf{\Upsilon}\mathbf{B}\zeta $$

Note that we want to prove $\zeta(t) \to 0$ when $t \to \infty$, and this is equivalent to prove $V(\zeta) \to 0$ when $t \to \infty$, where $V(\zeta)$ is a Lyapunov function defined as

$$ V(\zeta) = \zeta^{\top}\mathbf{P}\zeta $$

Then, the derivative of $V(\zeta)$ w.r.t. time $t$ yields

$$ \begin{aligned} \frac{\partial V}{\partial t} &= \dot{V} = \zeta^{\top}\mathbf{P}\dot{\zeta} + \dot{\zeta}^{\top}\mathbf{P}\zeta \\ &= \zeta^{\top}(\mathbf{P}\mathbf{A} + \mathbf{A}^{\top}\mathbf{P})\zeta + \zeta^{\top}(\mathbf{P}\boldsymbol{\Upsilon}\mathbf{B} + \mathbf{B}^{\top}\boldsymbol{\Upsilon}^{\top}\mathbf{P})\zeta \end{aligned} $$

Logically, if both (14) and (15) hold, we can conclude that

$$ \begin{aligned} \dot{V} &= -2\zeta^{\top}\mathbf{Q}\zeta + 2\zeta^{\top}\mathbf{P}\boldsymbol{\Upsilon}\mathbf{B}\zeta \\ &\leq -2\zeta^{\top}\mathbf{Q}\zeta + 2\|\mathbf{P}\boldsymbol{\Upsilon}\mathbf{B}\|\zeta^{\top}\zeta \\ &\leq -2\zeta^{\top}\mathbf{Q}\zeta + 2\lambda_{\max}(\mathbf{P})\|\boldsymbol{\Upsilon}\mathbf{B}\|\zeta^{\top}\zeta \\ &\leq -2\zeta^{\top}\mathbf{Q}\zeta + 2\frac{\lambda_{\max}(\mathbf{P})}{\lambda_{\min}(\mathbf{Q})}\|\boldsymbol{\Upsilon}\mathbf{B}\|\zeta^{\top}\boldsymbol{\Upsilon}\zeta \\ &\leq -2\left(1 - \frac{\lambda_{\max}(\mathbf{P})}{\lambda_{\min}(\mathbf{Q})}\|\boldsymbol{\Upsilon}\|\|\mathbf{B}\|\right)\frac{\lambda_{\min}(\mathbf{Q})}{\lambda_{\max}(\mathbf{P})}\zeta^{\top}\mathbf{P}\zeta \\ &\leq -2\left(1 - \varepsilon\frac{\lambda_{\max}(\mathbf{P})}{\lambda_{\min}(\mathbf{Q})}\|\mathbf{B}\|\right)\frac{\lambda_{\min}(\mathbf{Q})}{\lambda_{\max}(\mathbf{P})}V < 0 \end{aligned} $$

<figure><img src="image_4.png" /><figcaption>Fig. 3. Soft manipulator prototype actuated by cables [23].</figcaption></figure>

which guarantees that $V$ exponentially converges to 0. Thus, we can derive that $\zeta(t)$ exponentially converges to 0. ■

The controller designed in (11) is dependent on the pseudo-inverse of the matrix $\hat{\Gamma}$, while the RBFNN in Section III-C is not trained to provide this pseudo-inverse. To get rid of the inversion during the run time of the control, $\hat{\Gamma}^{\dagger}$ can also be directly obtained by using RBFNN as long as $\Gamma^{\dagger}$ is selected as the matrix base.

## IV. EXPERIMENTAL TESTING ON THE SOFT MANIPULATOR

First of all, we describe the soft manipulator control system, parameters setting of workspace estimation, and parameter tuning of RBFNN-based optimization as well. Next, a series of experiments are conducted to demonstrate the effectiveness and robustness of the proposed controller. In the end, the experimental comparison between the control scheme which is only valid in a certain sub-workspace and the proposed robust controller suitable for the whole workspace is carried out.

### A. Soft Manipulator Prototype Control System

The conical manipulator prototype is controlled by 4 independent cables which are attached to the tip and pass through the small holes of rigid rings mounted on the isotropic silicone body. To obtain the end-effector position of the soft manipulator, a long hole inside the whole body is made, and a magnetic sensor is placed at the tip of the soft manipulator to measure the end-effector position in real time. In detail, the exact geometric parameters of the investigated soft manipulator are illustrated in Fig. 3, and the identified material parameters from [23] are: Young's modulus $E = 2.563 \times 10^5$ Pa, shear modulus $G = 8.543 \times 10^4$ Pa, and density of material $\rho = 1.41 \times 10^3$ kg/m$^3$. In addition to the robot prototype, there are micro controller, position sensor, Matlab environment on the laptop, and 4 stepping motors for control of the cable length, as shown in Fig. 4.

As for the optimization-based approach, we propose to discretize the angles with a discretization step size of 0.05 Radian (Unit), 200 successive rays with respective direction vector $\mathcal{E}^{(i)}$, $i = 1, 2, \dots, 200$ emanating at angular intervals of the angles $\gamma_x = \gamma_y = \gamma_z = 2\pi/200$ from the initial radiating point which is exterior to the workspace. In this case, the PLS Cosserat model is divided into three sections. Finally, the
<figure><img src="image_5.png" /><figcaption>Fig. 4. Experimental setup. (a) A silicone soft manipulator; (b) End-effector; (c) Stepping motor; (d) Micro controller; (e) Polhemus magnetic position sensor; (f) Matlab environment on the laptop.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. The whole workspace of the investigated soft manipulator.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. From the left, the distribution of the representative end-effector positions selected from all sub-workspaces, followed by the front view of the selected points.</figcaption></figure>

whole workspace of the cable-driven soft manipulator shown in Fig. 5 can be obtained.

According to the proposed decomposition rule (4), the whole workspace was finally divided into 8 sub-workspaces where the 8 representative end-effector position points (as illustrated in Fig. 6) that are respectively distributed in their corresponding sub-workspace can be chosen. Calculating their associated configuration variables by solving the inverse model, we can then acquire 8 Jacobian matrices to form a matrix basis for estimating (off-line) $\Gamma(q)$ corresponding to any end-effector position $u \in W_E$.

In our test, the number of neurons required for the approximation of the Jacobian matrix $\Gamma(q)$ is determined by recording the achievable control performance with a given number of neurons, and then, increasing the number of neurons until the desired tracking effects are obtained. To determine the number of neurons, we have defined the relative error between end-effector position difference and its estimation as

$$e = \frac{1}{k} \sum_{j=1}^{k} \frac{\|u_j - u_{j-1} - (\sum_{i=1}^{m} \alpha_i (u_j) \Gamma_i)(T_j - T_{j-1})\|_2}{\|u_j - u_{j-1}\|_2}$$

and compared this metric $e$ by using RBFNN with different number of neurons. The smaller the metric $e$, the more accurate approximation of $\Gamma(q)$. To determine the number of the neurons, many trials were required in the simulation environment. Initially, 6 neurons were used, and then the

number was gradually increased to 27. We noticed that when utilizing more than 24 neurons, $e$ is almost the same as that of using 24 neurons. Based on those simulations, using 24 neurons, the Jacobian matrix has been approximated with sufficient accuracy such that the estimation error of Jacobian matrix satisfies Assumption 1, and the control performance has also reached the desired level. During the training process, the estimated weight matrix in (10) was initialized to zero, and the center of each neuron was dependent on the evenly chosen end-effector position across the whole workspace.

Since we trained (off-line) the RBFNN in the simulation environment, the sufficient data sets can be generated. We tried to evenly choose 40 equilibrium points from the whole workspace of the manipulator. For each point $j$, we can then collect 4 groups of data pairs near itself by applying a small perturbation to each cable. Thus, 160 sets of data pairs can be obtained. We selected 50% of the data sets as training sets which were enough to satisfy the control performance, and the remaining 50% for validation purpose. The convergence strategy of the RBFNN involves initializing the RBF centers appropriately and iteratively adjusting the weights using optimization techniques to minimize the objective function, and finally we use the gradient descent method to solve the linear optimization problem.

## B. Experimental Testing

To evaluate the performance of the proposed controller with the obtained Jacobian matrix, different experiments including slowly time-varying trajectory tracking tests and robustness analysis are performed on the soft manipulator prototype.

1) Slowly Time-varying Trajectory Tracking Tests: In this test, we used the RBFNN-based robust controller to track circle and star-shaped slowly time-varying trajectories. The results have been demonstrated in Fig. 7 and Fig. 8, clearly showcasing the soft manipulator can rapidly track the different slowly time-varying trajectories. Basically, the proposed PLS Cosserat static model-based controller is feasible and effective for the soft manipulator to move in its whole workspace.

2) Robustness Analysis: In an effort to demonstrate the robustness of the proposed controller, a small Teenage Mutant Ninja Turtle which can be considered as an externally permanent disturbance hanged on the soft body. Simultaneously,
<figure><img src="image_8.png" /><figcaption>Fig. 7. Circular time-varying trajectory tracking in Y-Z plane by applying the proposed global controller.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Star-shaped time-varying trajectory tracking in Y-Z plane by using the proposed global controller.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. The robustness validation of the soft manipulator from one point to another one under the combination of the temporary and permanent disturbances.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Tracking effect comparison of the slowly time-varying curve trajectory in 3D space between the local controller (green line) with a constant Jacobian matrix and the proposed controller (red line).</figcaption></figure>

the temporary disturbances are exerted on the soft manipu-
lator with hand. The manipulator reaches the desired point
**u**<sub>ref</sub> = [15.44 4.33 3.40]<sup>T</sup> from the initial equilibrium
point **u**<sub>0</sub> = [14.83 0.00 −7.53]<sup>T</sup>. It can be clearly seen
from the experimental results depicted in Fig. 9 that the
proposed controller is robust in the sense that end-effector of
the soft manipulator still returns to the desired position after
a small fluctuation.

C. Comparison Results and Discussion

It can be observed that the size of circle and star-shaped
trajectories for the global controller tracking are slightly
smaller than those for the local controller tracking introduced
in [24], which may not fully highlight the strengths of the pro-
posed control scheme despite its good tracking performances.
Therefore, to intuitively demonstrate the advantages of the
proposed controller, we have implemented a comparison of
the tracking performance for the slowly time-varying spatial
curve trajectory by using the proposed control scheme and

the controller with the constant Jacobian matrix. The result
indicates that the controller with a fixed Jacobian matrix for
the manipulator to track a slowly time-varying curve trajectory
around the whole workspace results in the divergence of the
system, as shown by the green line in Fig. 10. However, the
tracking effect for the global controller is excellent, as the
Jacobian matrix Γ of any end-effector position within the
whole workspace can be obtained. Thus, we can conclude that
global tracking performance has been achieved by unifying
RBFNN and PLS Cosserat.

To further showcase the superiority of the global control
method, a quantitative analysis of the tracking performance
of the studied soft manipulator given by the proposed global
controller and gain-scheduling controller in terms of maximum
error (ME) and root mean square error (RMSE) in task-space
trajectory of the end-effector has been performed and can
be found in Table I. From this table, it is confirmed that
the proposed global control scheme gives minimum tracking
errors (ME and RMSE) in X, Y and Z directions of different
trajectories (circle, star, and spatial curve) as compared to
the gain-scheduling controller. The used gains of the gain-
scheduling controller are displayed in Table II.

Through the tracking effect over the period of operation of
the system, we have found that the fine motion control possible
with the proposed global controller is extremely useful for the
positioning tasks, as it enables all the capabilities of the soft
manipulator to be exploited.

V. CONCLUSION

This contribution opens, for the first time, the possibility
of controlling the soft manipulator via a hybrid approach of
both PLS Cosserat static model and RBFNN to approximate
the Jacobian matrix between the task space and the actuator
space of the soft manipulator in its whole workspace, based on
which we designed a robust closed-loop controller to control
the end-effector position of the soft manipulator. Experimental
results for the soft manipulator following the fixed point
in the presence of external disturbances or moving along
slowly time-varying trajectories are presented to demonstrate
the performance (i.e., feasibility and robustness) provided by
the proposed global control technique. Additionally, the global
controller is independent of any specific robot's cross-section,
<table><caption>TABLE I
QUANTITATIVE ANALYSIS OF END-EFFECTOR POSITION TRACKING CONTROL ALONG X, Y AND Z DIRECTIONS FOR DIFFERENT CONTROLLERS</caption><thead><tr><th rowspan="2">Control schemes</th><th rowspan="2">Trajectory types</th><th colspan="2">u<sub>x</sub> (Unit:cm)</th><th colspan="2">u<sub>y</sub> (Unit:cm)</th><th colspan="2">u<sub>z</sub> (Unit:cm)</th></tr><tr><th>ME</th><th>RMSE</th><th>ME</th><th>RMSE</th><th>ME</th><th>RMSE</th></tr></thead><tbody><tr><td rowspan="4">Global controller</td><td>Circle</td><td>0.106</td><td>0.027</td><td>0.100</td><td>0.033</td><td>0.011</td><td>0.062</td></tr><tr><td>Star</td><td>0.091</td><td>0.021</td><td>0.096</td><td>0.011</td><td>0.098</td><td>0.042</td></tr><tr><td>Curve</td><td>0.128</td><td>0.103</td><td>0.117</td><td>0.091</td><td>0.104</td><td>0.085</td></tr><tr><td>Circle</td><td>0.209</td><td>0.142</td><td>0.214</td><td>0.161</td><td>0.257</td><td>0.194</td></tr><tr><td rowspan="2">Gain-scheduling controller [43]</td><td>Star</td><td>0.255</td><td>0.149</td><td>0.269</td><td>0.157</td><td>0.230</td><td>0.135</td></tr><tr><td>Curve</td><td>0.232</td><td>0.163</td><td>0.218</td><td>0.134</td><td>0.245</td><td>0.151</td></tr></tbody></table>

<table><caption>TABLE II
PARAMETER SETUPS OF THE GAIN-SCHEDULING CONTROLLER</caption><thead><tr><th>Trajectory types</th><th>λ<sub>1</sub></th><th>λ<sub>2</sub></th></tr></thead><tbody><tr><td rowspan="3">Circle</td><td rowspan="3">[ 15 15 3.5 ]</td><td rowspan="3">[ 14 14 2.5 ]</td></tr><tr><td rowspan="2">[ 20 20 7.0 ]</td></tr><tr><td>[ 18 18 8.5 ]</td></tr></tbody></table>

and thus it can be easily adapted to controlling a wide range of soft manipulators.

In the near future, the PLS Cosserat dynamic model-based controllers will be designed to achieve the strain and position control of the soft manipulator, and experiments on the soft prototype will be carried out to validate the tracking performances of the controllers.

## APPENDIX

### A. Derivation of $\frac{\partial u}{\partial q}$ and $\frac{\partial q}{\partial T}$

1) $\frac{\partial u}{\partial q}$: According to the definition of the velocity twist $\eta$, the linear velocity $V$ can be formulated as

$$V = R^{-1}(q)\dot{u} = R^{-1}(q)\frac{\partial u}{\partial q}\dot{q} \quad (16)$$

Combining the definition of velocity twist with differential model (1c), the linear velocity $V$ can also be written as

$$V = [0 \quad I_3]\eta = [0 \quad I_3]J(q)\dot{q} \quad (17)$$

Using the expressions (16) and (17), the partial derivative of $u$ with respect to $q$ can be deduced

$$\frac{\partial u}{\partial q} = R(q) [0_3, I_3] J(q) = [0_3, R(q)] J(q)$$

2) $\frac{\partial q}{\partial T}$: Taking the derivative of (1a) with respect to $T$, we have

$$\left( K + \frac{\partial[G(q)]}{\partial q} \right) \frac{\partial q}{\partial T} = H$$

It is then supposed that the cables are installed in a manner such that the matrix $\left( K + \frac{\partial[G(q)]}{\partial q} \right)$ deduced is invertible (if it is not, we then use pseudo-inverse), then the partial derivative of $q$ with respect to $T$ can be obtained

$$\frac{\partial q}{\partial T} = \left( K + \frac{\partial[G(q)]}{\partial q} \right)^{\dagger} H$$

## REFERENCES

[1] H. Banerjee, Z. T. H. Tse, and H. Ren, “Soft robotics with compliance and adaptation for biomedical applications and forthcoming challenges,” *Int. J. Robot. Autom.*, vol. 33, no. 1, pp. 69–80, 2018.

[2] F. Boyer, M. Porez, and W. Khalil, “Macro-continuous computed torque algorithm for a three-dimensional eel-like robot,” *IEEE Trans. Robot.*, vol. 22, no. 4, pp. 763–775, 2006.

[3] R. J. Webster III and B. A. Jones, “Design and kinematic modeling of constant curvature continuum robots: A review,” *Int. J. Robot. Res.*, vol. 29, no. 13, pp. 1661–1683, 2010.

[4] I. S. Godage, G. A. Medrano-Cerda, D. T. Branson, E. Guglielmino, and D. G. Caldwell, “Dynamics for variable length multisection continuum arms,” *Int. J. Robot. Res.*, vol. 35, no. 6, pp. 695–722, 2016.

[5] B. A. Jones and I. D. Walker, “Kinematics for multisection continuum robots,” *IEEE Trans. Robot.*, vol. 22, no. 1, pp. 43–55, 2006.

[6] A. D. Marchese and D. Rus, “Design, kinematics, and control of a soft spatial fluidic elastomer manipulator,” *Int. J. Robot. Res.*, vol. 35, no. 7, pp. 840–869, 2016.

[7] R. K. Katzschmann, C. Della Santina, Y. Toshimitsu, A. Bicchi, and D. Rus, “Dynamic motion control of multi-segment soft robots using piecewise constant curvature matched with an augmented rigid body model,” in *Proc. IEEE Int. Conf. Soft Robot. (RoboSoft)*, 2019, pp. 454–461.

[8] C. Della Santina, R. K. Katzschmann, A. Bicchi, and D. Rus, “Model-based dynamic feedback control of a planar soft robot: trajectory tracking and interaction with the environment,” *Int. J. Robot. Res.*, vol. 39, no. 4, pp. 490–513, 2020.

[9] C. Della Santina, A. Bicchi, and D. Rus, “On an improved state parametrization for soft robots with piecewise constant curvature and its use in model based control,” *IEEE Robot. Autom. Lett.*, vol. 5, no. 2, pp. 1001–1008, 2020.

[10] V. Falkenhahn, A. Hildebrandt, R. Neumann, and O. Sawodny, “Model-based feedforward position control of constant curvature continuum robots using feedback linearization,” in *Proc. IEEE Int. Conf. Robot. Autom.*, 2015, pp. 762–767.

[11] B. Caasenbrood, A. Pogromsky, and H. Nijmeijer, “Control-oriented models for hyperelastic soft robots through differential geometry of curves,” *Soft Robot.*, 2022.
[12] A. Ataka, T. Abrar, F. Putzu, H. Godaba, and K. Althoefer, "Model-based pose control of inflatable eversion robot with variable stiffness," *IEEE Robot. Autom. Lett.*, vol. 5, no. 2, pp. 3398-3405, 2020.

[13] C. Duriez, "Control of elastic soft robots based on real-time finite element method," in *Proc. IEEE Int. Conf. Robot. Autom.*, 2013, pp. 3982-3987.

[14] O. Goury and C. Duriez, "Fast, generic, and reliable control and simulation of soft robots using model order reduction," *IEEE Trans. Robot.*, vol. 34, no. 6, pp. 1565-1576, 2018.

[15] Z. Zhang, J. Dequidt, and C. Duriez, "Vision-based sensing of external forces acting on soft robots using finite element method," *IEEE Robot. Autom. Lett.*, vol. 3, no. 3, pp. 1529-1536, 2018.

[16] M. Thieffry, A. Kruszewski, C. Duriez, and T.-M. Guerra, "Control design for soft robots based on reduced-order model," *IEEE Robot. Autom. Lett.*, vol. 4, no. 1, pp. 25-32, 2018.

[17] M. Tummers, V. Lebastard, F. Boyer, J. Troccaz, B. Rosa, and M. T. Chikhaoui, "Cosserat rod modeling of continuum robots from newtonian and lagrangian perspectives," *IEEE Trans. Robot.*, vol. 39, no. 3, pp. 2360-2378, 2023.

[18] A. A. Alqumsan, S. Khoo, and M. Norton, "Robust control of continuum robots using cosserat rod theory," *Mech. Mach. Theory.*, vol. 131, pp. 48-61, 2019.

[19] F. Campisano, S. Caló, A. A. Remirez, J. H. Chandler, K. L. Obstein, R. J. Webster III, and P. Valdastri, "Closed-loop control of soft continuum manipulators under tip follower actuation," *Int. J. Robot. Res.*, vol. 40, no. 6-7, pp. 923-938, 2021.

[20] T. George Thuruthel, F. Renda, and F. Iida, "First-order dynamic modeling and control of soft robots," *Front. Robot. AI*, vol. 7, p. 95, 2020.

[21] F. Renda, C. Armanini, V. Lebastard, F. Candelier, and F. Boyer, "A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation," *IEEE Robot. Autom. Lett.*, vol. 5, no. 3, pp. 4006-4013, 2020.

[22] F. Renda, C. Armanini, A. Mathew, and F. Boyer, "Geometrically-exact inverse kinematic control of soft manipulators with general threadlike actuators' routing," *IEEE Robot. Autom. Lett.*, 2022.

[23] H. Li, L. Xun, and G. Zheng, "Piecewise linear strain cosserat model for soft slender manipulator," *IEEE Trans. Robot.*, vol. 39, no. 3, pp. 2342-2359, 2023.

[24] H. Li, L. Xun, G. Zheng, and F. Renda, "Discrete cosserat static model-based control of soft manipulator," *IEEE Robot. Autom. Lett.*, vol. 8, no. 3, pp. 1739-1746, 2023.

[25] T. George Thuruthel, E. Falotico, M. Manti, A. Pratesi, M. Cianchetti, and C. Laschi, "Learning closed loop kinematic controllers for continuum manipulators in unstructured environments," *Soft robot.*, vol. 4, no. 3, pp. 285-296, 2017.

[26] A. Melingui, O. Lakhal, B. Daachi, J. B. Mbede, and R. Merzouki, "Adaptive neural network control of a compact bionic handling arm," *IEEE/ASME Trans. Mechatronics*, vol. 20, no. 6, pp. 2862-2875, 2015.

[27] F. Piqué, H. T. Kalidindi, L. Fruzzetti, C. Laschi, A. Menciassi, and E. Falotico, "Controlling soft robotic arms using continual learning," *IEEE Robot. Autom. Lett.*, vol. 7, no. 2, pp. 5469-5476, 2022.

[28] G. Zheng, Y. Zhou, and M. Ju, "Robust control of a silicone soft robot using neural networks," *ISA Trans.*, vol. 100, pp. 38-45, 2020.

[29] T. G. Thuruthel, E. Falotico, F. Renda, and C. Laschi, "Learning dynamic models for open loop predictive control of soft robotic manipulators," *Bioinspiration biomimetics*, vol. 12, no. 6, p. 066003, 2017.

[30] G. Li, J. Shintake, and M. Hayashibe, "Deep reinforcement learning framework for underwater locomotion of soft robot," in *Proc. IEEE Int. Conf. Robot. Autom.*, 2021, pp. 12033-12039.

[31] D. Bruder, X. Fu, R. B. Gillespie, C. D. Remy, and R. Vasudevan, "Data-driven control of soft robots using koopman operator theory," *IEEE Trans. Robot.*, vol. 37, no. 3, pp. 948-961, 2020.

[32] J. M. Bern, Y. Schnider, P. Banzet, N. Kumar, and S. Coros, "Soft robot control with a learned differentiable model," in *Proc. IEEE Int. Conf. Soft Robot.*, 2020, pp. 417-423.

[33] G. Fang, Y. Tian, Z.-X. Yang, J. M. Geraedts, and C. C. Wang, "Efficient jacobian-based inverse kinematics with sim-to-real transfer of soft robots by learning," *IEEE/ASME Trans. Mechatronics*, 2022.

[34] T. George Thuruthel, Y. Ansari, E. Falotico, and C. Laschi, "Control strategies for soft robotic manipulators: A survey," *Soft robot.*, vol. 5, no. 2, pp. 149-163, 2018.

[35] M. Giorelli, F. Renda, M. Calisti, A. Arienti, G. Ferri, and C. Laschi, "Neural network and jacobian method for solving the inverse statics of a cable-driven soft arm with nonconstant curvature," *IEEE Trans. Robot.*, vol. 31, no. 4, pp. 823-834, 2015.

[36] T. G. Thuruthel, E. Falotico, F. Renda, and C. Laschi, "Model-based reinforcement learning for closed-loop dynamic control of soft robotic manipulators," *IEEE Trans. Robot.*, vol. 35, no. 1, pp. 124-134, 2018.

[37] C. C. Johnson, T. Quackenbush, T. Sorensen, D. Wingate, and M. D. Killpack, "Using first principles for deep learning and model-based control of soft robots," *Front. Robot. AI.*, vol. 8, p. 654398, 2021.

[38] M. Wiese, G. Runge-Borchert, B.-H. Cao, and A. Raatz, "Transfer learning for accurate modeling and control of soft actuators," in *Proc. IEEE Int. Conf. Soft Robot. (RoboSoft)*, 2021, pp. 51-57.

[39] Z. Q. Tang, H. L. Heung, K. Y. Tong, and Z. Li, "Model-based online learning and adaptive control for a "human-wearable soft robot" integrated system," *Int. J. Robot. Res.*, vol. 40, no. 1, pp. 256-276, 2021.

[40] A. Walid, G. Zheng, A. Kruszewski, and F. Renda, "Discrete cosserat method for soft manipulators workspace estimation: An optimization-based approach," *J. Mechan. Robot.*, vol. 14, no. 1, p. 011012, 2021.

[41] A. El-Bakry, R. A. Tapia, T. Tsuchiya, and Y. Zhang, "On the formulation and theory of the newton interior-point method for nonlinear programming," *J. Optimiz. theory. App.*, vol. 89, no. 3, pp. 507-541, 1996.

[42] R. Hecht-Nielsen, "Theory of the backpropagation neural network," in *Neural networks for perception*, 1992, pp. 65-93.

[43] K. Wu and G. Zheng, "Fem-based gain-scheduling control of a soft trunk robot," *IEEE Robot. Autom. Lett.*, vol. 6, no. 2, pp. 3081-3088, 2021.