# The Logarithmic Random Bidding for the Parallel Roulette Wheel Selection
with Precise Probabilities

Koji Nakano  Graduate School of Advanced Science and Engineering, Hiroshima
University
Kagamiyama 1-4-1, Higashihiroshima, 739-8527 Japan

###### Abstract

The roulette wheel selection is a critical process in heuristic algorithms,
enabling the probabilistic choice of items based on assigned fitness values.
It selects an item with a probability proportional to its fitness value. This
technique is commonly employed in ant-colony algorithms to randomly determine
the next city to visit when solving the traveling salesman problem. Our study
focuses on parallel algorithms designed to select one of multiple processors,
each associated with fitness values, using random wheel selection. We propose
a novel approach called logarithmic random bidding, which achieves an expected
runtime logarithmic to the number of processors with non-zero fitness values,
using the CRCW-PRAM model with a shared memory of constant size. Notably, the
logarithmic random bidding technique demonstrates efficient performance,
particularly in scenarios where only a few processors are assigned non-zero
fitness values.

###### Index Terms:

roulette wheel selection, fitness proportionate selection, ant colony
optimization, traveling salesman problem, parallel heuristic algorithms

##  I Introduction

Let $f_{0},f_{1},\ldots,f_{n-1}$ denote non-negative real numbers referred to
as _fitness_. _The roulette wheel selection_ , also known as _the fitness
proportionate selection_ , is a fundamental operation that chooses one of the
$n$ indices in a way that the probability of selecting an index $i$ ($0\leq
i\leq n-1$) is directly proportional to its fitness value. Mathematically,
this probability of selecting $i$ is determined as:

$\displaystyle F_{i}={f_{i}\over f_{0}+f_{1}+\cdots+f_{n-1}}.$

The roulette wheel selection serves as a crucial component in heuristic
algorithms. For instance, in the ant colony optimization technique for solving
the traveling salesman problem (TSP) [1, 2, 3], this method is employed to
choose the next city to visit. Specifically, each edge connecting a pair of
cities is assigned a fitness value based on its suitability as a TSP route.
Subsequently, the next city to visit is randomly selected from those connected
to an unvisited city, with the selection probability being directly
proportional to the fitness value of the edge. Furthermore, roulette wheel
selection can be applied to the vertex coloring problem [4].

This study focuses on parallel roulette wheel selection, aiming to efficiently
select a processor $i$ among $n$ processors, each assigned a fitness value
$f_{i}$ corresponding to their ID from 0 to $n-1$. The objective of parallel
roulette wheel selection is to randomly select a processor $i$ with
probability $F_{i}$. For precise theoretical analysis, we assume the PRAM
(parallel random access machine model) [5], equipped with multiple processors
and the shared memory with either EREW (exclusive read exclusive write) or
CRCW (concurrent read concurrent write) operations. All processors function
synchronously, and in the EREW-PRAM model, simultaneous memory access by
multiple processors is prohibited. On the other hand, in the CRCW-PRAM model,
such simultaneous access is allowed. In this CRCW-PRAM model, if a write
conflict occurs in a memory cell of the shared memory, a randomly selected one
among the multiple memory write operations succeeds in writing the value.

Utilizing a pseudo-random number generator rand() that returns a real number
uniformly distributed in the range $[0,1)$, we can implement the roulette
wheel selection algorithm. We define $p_{i}=f_{0}+f_{1}+\cdots+f_{i}$ ($0\leq
i\leq n-1$) as the prefix-sum, with $p_{-1}=0$ for simplicity. The parallel
roulette wheel selection can be executed using the prefix-sum-based algorithm
[3] outlined below:

[Prefix-sum-based parallel roulette wheel selection]
---
1. | Compute all prefix-sums $p_{0},p_{1},\ldots,p_{n-1}$.
2. | Processor 0 computes $R\leftarrow{\rm rand()}\cdot p_{n-1}$.
3. | Each processor $i$ ($0\leq i\leq n-1$) checks if $p_{i-1}\leq R<p_{i}$.
|  | If this condition holds, processor $i$ is selected.

This algorithm selects processor $i$ with probability:

$\displaystyle{p_{i}-p_{i-1}\over p_{n-1}}$ $\displaystyle={f_{i}\over
f_{0}+f_{1}+\cdots+f_{n-1}}=F_{i}.$

Consequently, the prefix-sum-based algorithm correctly selects $i$ with the
probability $F_{i}$. We can use a parallel prefix-sum algorithm on the EREW-
PRAM [5], which operates in $O(\log n)$ time using $O(n)$ memory cells to
compute all parallel prefix-sums. Thus, the prefix-sum-based parallel roulette
wheel selection runs in $O(\log n)$ time on the EREW-PRAM with a shared memory
of size $O(n)$.

From a practical standpoint, the independent roulette selection [6] might
offer increased efficiency. The algorithm for independent roulette wheel
selection is outlined below:

[The independent roulette wheel selection]
---
1. | Each processor $i$ ($0\leq i<n$) computes $r_{i}=f_{i}\cdot rand()$.
2. | Identify the maximum $r_{i}$ among $r_{0},r_{1},\ldots,r_{n-1}$
|  | and select processor $i$.

Since $r_{i}$ is a random number in the range $[0,f_{i})$, a larger $f_{i}$
has a higher probability of becoming the maximum. Thus, a processor with
higher fitness is selected with a higher probability. However, the probability
may not equal to $F_{i}$. This discrepancy can be observed through a simple
example with $n=2$, $f_{0}=2$, and $f_{1}=1$. In this case, $r_{0}\in[0,2)$
and $r_{1}\in[0,1)$. If $r_{0}\geq[1,2)$, then processor 0 is chosen.
Otherwise, since both $r_{0}$ and $r_{1}$ are in $[0,1)$, processor 0 is
selected with a probability of $1\over 2$. Hence, processor 0 is chosen with a
probability of ${1\over 2}+{1\over 2}\cdot{1\over 2}={3\over 4}$ using the
independent roulette wheel selection, while the roulette wheel selection
requires selecting 0 with a probability of $F_{0}={2\over{2+1}}={2\over 3}$.
Thus, the independent roulette wheel selection fails to adhere to the desired
probabilities of the roulette wheel selection. For a comprehensive analysis of
selection probabilities in the independent roulette wheel selection, readers
are encouraged to refer to [7].

This paper primary focuses on presenting a novel technique for the parallel
roulette wheel selection method that ensures precise probabilities. Our new
technique called _the logarithmic random bidding_ is designed so that each $i$
($0\leq i\leq n-1$) is selected with a probability of $F_{i}$ within an
expected time complexity of $O(\log k)$ on the CRCW-PRAM with a shared memory
cell of size $O(1)$, where $k$ denotes the number of non-zero fitness values
among $f_{0},f_{1},\ldots,f_{n-1}$. Hence, the logarithmic random bidding
technique operates quite fast, especially when $k$ is small. Notably, in ant-
colony based TSP algorithms, fitness values are often set to zero for cities
that have already been visited. In such scenarios with many zero fitness
values, the logarithmic random bidding technique exhibits accelerated
performance.

This paper is organized as follows: Section II introduces our logarithmic
random bidding technique and demonstrates its ability to select processor $i$
with a probability of $F_{i}$. Section III details the implementation of a
parallel roulette wheel selection with the logarithmic random bidding
technique on the CRCW-PRAM. Finally, Section IV presents the conclusion of our
work.

##  II The parallel roulette wheel selection with the logarithmic random
bidding

This section presents _the logarithmic random bidding_ for the parallel
roulette wheel selection and provides a proof ensuring that it selects
processor $i$ ($0\leq i\leq n-1$) with a probability of $F_{i}$. This
algorithm closely resembles the independent roulette wheel selection, yet it
incorporates a crucial modification in computing $r_{i}$. Instead of the
conventional $r_{i}=f_{i}\cdot\text{rand()}$ method used in the independent
roulette wheel selection, it employs the logarithmic random bidding:
$r_{i}={\log({\rm rand()})\over f_{i}}$. The specifics are detailed below:

[Roulette wheel selection with logarithmic random bidding]
---
1. | Each processor $i$ ($0\leq i<n$) computes $r_{i}={\log({\rm rand()})\over f_{i}}$.
2. | Identify the maximum $r_{i}$ among $r_{0},r_{1},\ldots,r_{n-1}$
|  | and select processor $i$.

Next, we will demonstrate that the logarithmic random bidding technique
selects processor $i$ ($0\leq i\leq n-1$) with a probability of $F_{i}$. It is
evident that $r_{i}$ falls within the range $(-\infty,0)$. The cumulative
distribution function of $r_{i}$ is

$\displaystyle\Pr\left(r_{i}\leq x\right)$ $\displaystyle=\Pr\left({\log({\rm
rand()})\over f_{i}}\leq x\right)$ $\displaystyle=\Pr({\rm rand()}\leq
e^{xf_{i}})$ $\displaystyle=e^{xf_{i}}.$

Consequently, the probability density function of $r_{i}$ becomes:

$\displaystyle{d\Pr(r_{i}\leq x)\over dx}$ $\displaystyle=f_{i}e^{xf_{i}}$

Let us evaluate the probability that processor 0 is selected, that is, $r_{0}$
is larger than $r_{1},r_{2},\ldots,r_{n-1}$. This probability can be evaluated
as follows.

$\displaystyle\int_{-\infty}^{0}{d\Pr(r_{0}\leq x)\over dx}\cdot\Pr(r_{1}\leq
x)\cdots\Pr(r_{n-1}\leq x)dx$
$\displaystyle=\int_{-\infty}^{0}f_{0}e^{xf_{0}}e^{xf_{1}}e^{xf_{2}}\cdots
e^{xf_{n-1}}dx$ $\displaystyle=\int_{-\infty}^{0}f_{0}e^{x(f_{0}+f_{1}+\cdots
f_{n-1})}dx$ $\displaystyle=\left[{f_{0}\over
f_{0}+f_{1}+\cdots+f_{n-1}}e^{x(f_{0}+f_{1}+\cdots
f_{n-1})}\right]_{-\infty}^{0}$ $\displaystyle={f_{0}\over
f_{0}+f_{1}+\cdots+f_{n-1}}=F_{0}$

Therefore, processor 0 is selected with a probability of $F_{0}$. Similarly,
the probability of selecting any other processor $i$ can be determined in the
same manner, confirming that processor $i$ is selected with a probability of
$F_{i}$.

TABLE I: The selection probabilities of the roulette wheel selection algorithms in $10^{9}$ iterations with $f_{i}=i$ ($0\leq i\leq 9$) $i$ | $f_{i}$ | $F_{i}$ | independent | logarithmic
---|---|---|---|---
0 | 0 | 0.000000 | 0.000000 | 0.000000
1 | 1 | 0.022222 | 0.000000 | 0.022222
2 | 2 | 0.044444 | 0.000088 | 0.044446
3 | 3 | 0.066667 | 0.001708 | 0.066672
4 | 4 | 0.088889 | 0.010993 | 0.088885
5 | 5 | 0.111111 | 0.038787 | 0.111105
6 | 6 | 0.133333 | 0.094267 | 0.133340
7 | 7 | 0.155556 | 0.178238 | 0.155552
8 | 8 | 0.177778 | 0.282382 | 0.177771
9 | 9 | 0.200000 | 0.393536 | 0.200007
TABLE II: The selection probabilities of the first 10 processors of the roulette wheel selection algorithms in $10^{9}$ iterations with $f_{0}=1$ and $f_{1}=f_{2}\cdots=f_{99}=2$ $i$ | $f_{i}$ | $F_{i}$ | independent | logarithmic
---|---|---|---|---
0 | 1 | 0.005025 | 0.000000 | 0.005026
1 | 2 | 0.010050 | 0.010104 | 0.010053
2 | 2 | 0.010050 | 0.010104 | 0.010053
3 | 2 | 0.010050 | 0.010099 | 0.010048
4 | 2 | 0.010050 | 0.010101 | 0.010050
5 | 2 | 0.010050 | 0.010108 | 0.010057
6 | 2 | 0.010050 | 0.010099 | 0.010048
7 | 2 | 0.010050 | 0.010101 | 0.010050
8 | 2 | 0.010050 | 0.010100 | 0.010050
9 | 2 | 0.010050 | 0.010106 | 0.010055

To verify the probability precision of the logarithmic random bidding
technique against the conventional independent roulette wheel selection, we
conducted simulations over $10^{9}$ iterations using $f_{i}=i$ ($0\leq i\leq
9$). The Mersenne Twister random number generator [8] was employed to
implement the rand() function. The comparative analysis showcases the
inaccuracy of the independent roulette wheel selection, particularly evident
for smaller $f_{i}$ values where the selection probability significantly
deviates from $F_{i}$. In contrast, our roulette wheel selection consistently
demonstrates probabilities of selecting $f_{i}$ that closely align with
$F_{i}$.

Especially for large $n$, the probability of selecting smaller $f_{i}$ values
can tend toward near-zero values. For instance, in a scenario employing the
roulette wheel selection across 100 processors with $f_{0}=1$ and
$f_{1}=f_{2}=\cdots=f_{99}=2$, the expected selection probability for
processor 0 stands at $\frac{1}{199}\approx 0.005025$. However, the
independent roulette wheel selection yields a probability of
$({\frac{1}{2}})^{99}\cdot\frac{1}{100}\approx 1.57772\times 10^{-32}$,
essentially zero, resulting in processor 0 never being selected. Table II
shows simulations from $10^{9}$ iterations, confirming alignment with
mathematical analyses: the independent roulette wheel selection consistently
neglects selection of processor 0, while the logarithmic random bidding
technique accurately reflects the expected selection probabilities.

##  III Implementation of the parallel roulette wheel selection algorithm

This section focuses on implementations of the parallel roulette wheel
selection. In particular, we will present how we can identify the maximum
$r_{i}$ among $r_{0},r_{1},\ldots,r_{n-1}$ in parallel. It can be identified
by a parallel reduction in an obvious way as follows. Imagine a binary tree
with $n$ leaves each associated with $r_{i}$. The maximum of the two children
are computed in very internal nodes from the leaves synchronously. Clearly,
the root will store the maximum. However, this implementation takes $O(\log
n)$ time and requires a shared memory of size $O(n)$ on the EREW-PRAM.

We aim to demonstrate that the maximum value $r_{i}$ can be identified within
an expected $O(\log k)$ parallel steps. The parallel algorithm utilizes the
CRCW-PRAM model, employing shared memory variable $s$ initialized to zero,
along with the variable ${\it output}$ to store the index of selected
processor.

Each processor $i$ iteratively writes its $r_{i}$ value into $s$ until the
condition $s\geq r_{i}$ is satisfied. Consequently, write conflicts may arise,
with one of the writing accesses successfully updating the shared memory cell.
Once all processors complete this writing operation and the condition $s\geq
r_{i}$ is met for all $i$, the value stored in $s$ represents the maximum of
all $r_{i}$’s. Each processor $i$ writes its index $i$ into ${\it output}$ if
$s=r_{i}$. Clearly, ${\it output}$ retains the index corresponding to the
maximum $r_{i}$ value. Here are the detailed steps:

[Identifying the maximum $r_{i}$]
---
Each processor $i$ ($0\leq i<n-1$) performs:
1. | while $s<r_{i}$ do $s\leftarrow r_{i}$;
2. | barrier_synchronization();
3. | if $s=r_{i}$ then ${\it output}\leftarrow i$;

Next, we will estimate the number of iteration performed in the while loop. We
assume that we have $k$ non-zero $r_{i}$’s. We say that a processor $i$ is
_active_ in an iteration of the while loop if $s<r_{i}$ holds and it performs
$s\leftarrow r_{i}$. Clearly, the first iteration, $k$ processors are active
and the while loop is iterated until no active processor exists. We say that
the an iteration is success if at least a half of active processors become
inactive. The iterations can have up to $\lceil\log_{2}k\rceil$ success
iterations and every iteration is success with probability $1\over 2$ when
$k\geq 2$. The expected number of iterations is

$\displaystyle 1\cdot{1\over 2}+2\cdot{1\over 2^{2}}+3\cdot{1\over
2^{3}}+\cdots$ $\displaystyle=2.$

Hence, $2\lceil\log_{2}k\rceil$ iterations are sufficient to have
$\lceil\log_{2}k\rceil$ success iterations, and we have the following theorem:

######  Theorem 1.

The parallel roulette wheel selection with the logarithmic random bidding for
$n$ processors, each with a non-negative fitness $f_{i}$ ($0\leq i\leq n-1$),
selects a processor $i$ with a probability of $F_{i}={f_{i}\over
f_{0}+f_{1}+\cdots+f_{n-1}}$ in $O(\log k)$ expected time on the CRCW-PRAM
with a shared memory size of $O(1)$, where $k$ is the number of non-zero
fitness values among all $f_{i}$’s.

##  IV Conclusion

This paper presented the logarithmic random bidding technique for the parallel
random wheel selection with precise selection probabilities. It runs in
$O(\log k)$ time on the CRCW-PRAM with a shared memory of size $O(1)$ and the
selection probabilities follow the requirement of the roulette wheel selection
precisely.

## References

  * [1] M. Dorigo and L. M. Gambardella, “Ant colonies for the travelling salesman problem,” _Biosystems_ , vol. 43, no. 2, pp. 73–81, 1997.
  * [2] M. Pedemonte, S. Nesmachnow, and H. Cancela, “A survey on parallel ant colony optimization,” _Applied Soft Computing_ , vol. 11, no. 8, pp. 5181–5197, 2011.
  * [3] A. Uchida, Y. Ito, and K. Nakano, “An efficient GPU implementation of ant colony optimization for the traveling salesman problem,” in _Proc. of International Conference on Networking and Computing_. IEEE CS Press, Dec. 2012, pp. 94–102.
  * [4] R. Murooka, Y. Ito, and K. Nakano, “Accelerating ant colony optimization for the vertex coloring problem on the GPU,” in _Proc. of International Symposium on Computing and Networking_ , 2016, pp. 469–475.
  * [5] A. Gibbons and W. Rytter, _Efficient Parallel Algorithms_. Cambridge University Press, 1988.
  * [6] J. M. Cecilia, J. M. García, A. Nisbet, M. Amos, and M. Ujaldón, “Enhancing data parallelism for ant colony optimization on GPUs,” _Journal of Parallel and Distributed Computing_ , vol. 73, no. 1, pp. 42–51, 2013.
  * [7] H. Lloyd and M. Amos, “Analysis of independent roulette selection in parallel ant colony optimization,” in _Proc. of the Genetic and Evolutionary Computation Conference_ , 2017, p. 19–26.
  * [8] M. Matsumoto and T. Nishimura, “Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator,” _ACM Transactions on Modeling and Computer Simulation_ , vol. 8, no. 1, pp. 3–30, Jan. 1998.


Simultaneous Position-and-Stiffness Control of Underactuated Antagonistic Tendon-Driven Continuum Robots
Bowen Yi, Yeman Fan, Dikai Liu, and Jose Guadalupe Romero

This work was partially supported by the Australian Research Council (ARC) Discovery Project under Grant DP200102497, and by the Robotics Institute at the University of Technology Sydney, Australia. B. Yi and Y. Fan contributed equally to this work. (Corresponding author: Y. Fan)
B. Yi, Y. Fan and D. Liu are with Robotics Institute, Faculty of Engineering and Information Technology, University of Technology Sydney, Sydney, NSW 2006, Australia (email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
J.G. Romero is with Departamento Académico de Sistemas Digitales, ITAM, Río Hondo 1,  01080, Ciudad de México, México (email<EMAIL_ADDRESS>
Continuum robots have gained widespread popularity due to their inherent compliance and flexibility, particularly their adjustable levels of stiffness for various application scenarios. Despite efforts to dynamic modeling and control synthesis over the past decade, few studies have focused on incorporating stiffness regulation in their feedback control design; however, this is one of the initial motivations to develop continuum robots. This paper aims to address the crucial challenge of controlling both the position and stiffness of a class of highly underactuated continuum robots that are actuated by antagonistic tendons. To this end, the first step involves presenting a high-dimensional rigid-link dynamical model that can analyze the open-loop stiffening of tendon-driven continuum robots. Based on this model, we propose a novel passivity-based position-and-stiffness controller adheres to the non-negative tension constraint. To demonstrate the effectiveness of our approach, we tested the theoretical results on our continuum robot, and the experimental results show the efficacy and precise performance of the proposed methodology.

continuum robot, nonlinear control, stiffness regulation, energy shaping

§ INTRODUCTION

Continuum robots are a novel class of soft robots within the realm of robotics. Their unique properties, such as scaled dexterity and mobility, make them well facilitated important and suitable for human-robot interaction and manipulation tasks in uncertain and complex environments. For example, they can be used for manipulating objects with unknown shapes, performing search and rescue operations, and whole-arm grasping [17]. By utilizing soft materials, continuum robots represent an concrete step towards developing robots with performance levels similar to those of biological organisms.

Despite the above advantages, rigid-body robotics currently outperforms continuum robotics in tasks requiring adaptable movement and compliant interactions with environment [9]. Consequently, many efforts have been devoted to addressing the challenges for real-time control of continuum robots that facilitates fast, efficient and reliable operation. We refer the reader to the recent survey papers  [8, 36, 31] for a comprehensive understanding of this point. The existing control approaches for soft robots can be broadly classified into two categories: data driven (or machine learning) and model-based design. Initially, the data-driven approaches dominated the research because obtaining reliable models of a continuum robot was believed to be overwhelmingly complex [8]. In recent years, various state-of-the-art learning methodologies have been applied to control of soft robotics. These include the Koopman operator [6], Gaussian process temporal difference learning [11], supervised learning via recurrent neural networks [32], and feedforward neural networks [5]. However, it is well known that these data-driven approaches have some key limitations, including stringent requirements of data sets and no guarantee of stability or safety. Nonetheless, recent efforts from the control community have been made to impose stability constraints on model and learning policy [13, 37, 33]. On the other hand, in the last few years the resurgence of interest in model-based control approaches has made them particularly appealing for soft robotics. This is because they are robust even when using approximate models, and possess properties that are both interpretable and manageable [8].

Unlike rigid manipulators, elastic deformation of continuum robots theoretically leads to infinite degrees-of-freedom (DoF) motion, e.g., torsion, buckling, extension and bending. This renders them particularly suitable to be modelled by partial differential equations (PDEs) rather than conventional ordinary differential equations (ODEs) [7]. In particular, there are two prevalent categories of modeling used in this field, namely mechanics-based and geometry-based approaches. The former focuses on studying the elastic behaviour of the constitutive materials and solving the boundary conditions problem. For instance, Cosserat rod theory and Euler-Bernoulli beam theory are two methodologies for modelling commonly used in this category [34]. To enable numerical implementation of these models, they need to be solved numerically to obtain a closed formulation for each material subdomain. The finite elements method (FEM) has proven to be successful in the design and analysis of continuum robots with high accuracy [1]. The FEM approaches have extremely heavy computational burden, thus not adopted to real-time control with a notable exception of [4] for quasi-static control. In contrast, geometrical models assume that the soft body can be represented by a specific geometric shape, e.g., piecewise constant curvature (PCC). As these modelling approaches often lead to kinematic models rather than dynamical models, they enable the design of kinematic or quasi-static controllers. For example, [16] provides the first exact kinematics in closed forms of trunk sections that can extend in length and bend in two dimensions, and shows how to connect workspace coordinates to actuator inputs. However, some results have shown that such types of kinematic controllers are likely to yield poor closed-loop performance [17].

To address the challenges previously mentioned, research has been conducted in the last few years on the dynamic modelling and model-based control of continuum robots.[In this paper, we use the term “dynamic controllers” to refer to feedback laws designed from dynamical and kinematic models for robotics, rather than the terminology in control theory which refers to feedback control with dynamics extension (e.g. adaptive and observer-based control) [25].] Several dynamical models have been adopted for controller synthesis of continuum robots, including the geometrically exact dynamic Cosserat model [27], port-Hamiltonian Cosserat model [7], rigid-link models [15, 9], and reduced-order Euler-Lagrangian model [12, 10]. These works have applied various model-based constructive control approaches, such as passivity-based control (PBC), partial feedback linearisation, proportional derivative (PD) control, and immersion and invariance (I&I) adaptive control. Among these, the paper [9] provides probably the earliest solution to the design and experimental validation of dynamic feedback control for soft robots.

As illustrated above, one of the primary motivations to develop continuum robots is to enable robots exhibit more agile and adaptable movement, and compliant interactions. Consequently, there is an urgent and rapidly growing need to develop high-performance control algorithms to regulate position and stiffness simultaneously. In certain robotic applications involving human interaction or in complicated environments (e.g. search and rescue, industrial inspection, medical service, and home living care), stiffness control outperforms position control in terms of performance. For instance, in minimally invasive surgery or natural orifice transluminal endoscopic surgery, precise control of the robot's position is crucial. Additionally, the robot must dynamically adjust its stiffness to minimize side effects or transmit force. The problem of stiffness control, together with impedance control, is well established for rigid and softly-actuated (or articulated soft) robotics [3, 20, 22]. In contrast, the problem of simultaneous position-and-stiffness control for continuum robots is still an open area of research. The first stiffness controller for continuum robots in the literature may refer to [21], which can be viewed as an extension of a simple Cartesian impedance controller for continuum robots using a kinematic model. In [2], the authors tailor the classic hybrid motion/force controller for a static model of multi-backbone continuum robots, and the proposed method needs the estimation of external wrenches. In [9], Cartesian stiffness controller was proposed to achieve dynamic control of a fully-actuated soft robot, enabling interaction between the robot and its environment. However, note that these works are not applicable to underactuated dynamical models of continuum robots.

This paper aims to address the above gap by proposing a novel dynamical model and a real-time control approach that regulates both position and stiffness concurrently for a class of underactuated antagonistic tendon-driven continuum robots. The main contributions of the paper are:

C1 We propose a port-Hamiltonian dynamical model for a class of antagonistic tendon-driven continuum robots, which features a configuration-dependent input matrix that enables to interpret the underlying mechanism for open-loop stiffening.

C2 Stiffness flexibility is one of the motivations to develop continuum robots. Using the resulting highly underactuated dynamical model, we propose a novel potential energy shaping controller. Though simultaneous position-and-stiffness control has been widely studied for rigid and softly-actuated robots, to the best of the authors' knowledge, this work is the first such attempt to design a controller, which is capable of simultaneous control of an underactuated continuum robot.

C3 We analyse the assignable equilibria set of the proposed model class, which is actuated by tendons that only provide non-negative tensions, and show how to incorporate input constraints into controller design via an input transformation.

We conducted experiments in various scenarios to validate the theoretical results presented in the paper. However, due to C2 and the lack of applicable control approaches in the existing literature, we were unable to include a fair experimental comparison with previous work in this study.

The reminder of the paper is organised as follows. Section <ref> provides a brief introduction to the proposed dynamical model of continuum robots and the problem studied in the paper. In Section <ref>, we explain the underlying reason why the proposed model can interpret the open-loop stiffening. In Section <ref> we present a simultaneous position-and-stiffness controller which is further discussed in Section <ref>. It is followed by the experimental results which were tested on a robotic platform OctRobot-I in Section <ref>. Finally, the paper presents some concluding remarks and discusses future work in Section <ref>.

Notation. All functions and mappings are assumed $C^2$-continuous. $I_n$ is the $n \times n$ identity matrix, $0_{n \times s}$ is an
$n \times s$ matrix of zeros, and $\mathbf{1}_n := \col(1, \ldots, 1) \in \rea^n$. Throughout the paper, we adopt the convention of using bold font for variables denoting vectors, while scalars and matrices are represented in normal font. For $\bfx \in \rea^n$, $S \in \rea^{n \times n}$, $S=S^\top
>0$, we denote the Euclidean norm $|\bfx|^2:=\bfx^\top \bfx$, and the weighted–norm $\|\bfx\|^2_S:=\bfx^\top S \bfx$. Given a function $f:  \rea^n \to \rea$ we define the differential operators
\nabla f:=(\frac{\partial f }{ \partial x})^\top,\;\nabla_{x_i} f:=(\frac{
\partial f }{ \partial x_i})^\top,
where $x_i \in \rea^p$ is an element of the vector $\bfx$. The set $\caln$ is defined as $\caln:= \{1,\ldots,n\}$. When clear from the context, the arguments of the functions may be omitted.

§ MODEL AND PROBLEM SET

§.§ Modelling of A Class of Continuum Robots

In this section, we present a control-oriented high-dimensional rigid-link dynamical model specifically designed for a class of underactuated continuum robots driven by tendons. This model class encompasses a wide range of recently reported continuum robots in the literature, including the elephant trunk-inspired robot [39] and our own developed OctRobot-I [14], alongside other notable examples. By employing this versatile model, we aim to provide a general framework that can effectively describe and analyse various underactuated continuum robotic systems, enabling a deeper understanding of their stiffening mechanisms and facilitating control design.

In order to visualize the modelling process, we take OctRobot-I as an example to introduce the proposed dynamical model, but keep its generality in mind that the model is not limited to this specific robotic platform. This robot imitates an octopus tentacle's structure and motion mechanism, as shown in Fig. <ref>. The whole continuum manipulator consists of several sections in order to be able to deform in three-dimensional space. Each of them is made of $n$ connected spine segments that are driven by a pair of cables. More details of the continuum robot OctRobot-I are given in Section <ref>, as well as in [14].

(a) A photo of the continuum robot OctRobot-I. (b) Mechanical structure of a class of continuum robots.

In this paper, we use a rigid-link model to approximate dynamical behaviours of continuum robots due to three considerations: 1) for the class of continuum robots which we are concerned with, there are “natural” spine segments partitioning into several links; 2) the rigid-link modelling approach is simple for control-oriented tasks; and 3) it is convenient to account for external loading.

To obtain the dynamical model, we make the following assumptions:

- The manipulator deflection is due to bending only, and the extension is negligible.

- In this paper, we specifically concentrate on the two-dimensional case, limiting our analysis to a single section, in order to effectively illustrate the underlying mechanism.[It is promising to extend the main results to the three-dimensional case with multi-sections, and we will consider it as a valuable avenue for further exploration.]

- The actuator dynamics is negligible, i.e., the motor is operating in the torque control mode with sufficiently short transient stages.

- The sections of continuum robots have a piecewise constant curvature (PCC), conformally with respect to segments.[Each spine segment has constant curvatures but variable in time. ]

In the rigid-link dynamical model, we use a serial chain of rigid links with $n$ rotational joints to approximate one section of the continuum robot. Then, the configuration variable can be defined as

\bfq = \begmat{q_1&  \ldots&  q_n}^\top \in \calx \subset  \rea^n,

with $q_i \in \caln$ representing the approximate link angles, where $\calx$ is the feasible configuration space; see Fig. <ref> for an illustration. Practically, all angles $q_i$ are in the set  $[-{\pi \over 2}, {\pi \over 2}]$ due to physical constraints.

We model the continuum robot as a port-Hamiltonian system in the form of [23, 35]

\begin{equation}
\label{model:pH}
\begmat{ ~\dot \bfq ~\\ ~\dot \bfp~}
\begmat{0_{n\times n} & I_n \\ - I_n & -D(\bfq)} \begmat{\nabla	_{\bfq} H \\ \nabla_{\bfp} H }+ \begmat{~\mathbf{0}_n ~ \\ ~G(\bfq) \bfu + \tau_{\tt ext}~}
\end{equation}

with the generalised momenta $\bfp \in \rea^n$, the damping matrix $D(\bfq) \in \rea^{n\times n}_{\succ 0}$, $\tau_{\tt ext} \in \rea^n$ the external torque, and the input matrix $G(\bfq) \in \rea^{n\times m}$ with $m<n$. The total energy of the robotic system is given by

\begin{equation}
\label{H}
H(\bfq,\bfp) = {1\over2} \bfp^\top M^{-1}(\bfq)\bfp + U(\bfq),
\end{equation}

with the inertia matrix $M(\bfq) \succ 0$, and the potential energy $U(\bfq)$, which consists of the gravitational part $U_{\tt G}$ and the elastic one $\ue(\bfq)$, i.e.

U(\bfq) ~=~ \ug(\bfq) + \ue(\bfq).

The potential energy function has an isolated local minimum point at its open-loop equilibrium of the origin $\bfq = \mathbf{0}_n$.

An illustration of configuration variables

The variable $\bfu \in \rea^m$ represents the control input, denoting the tensions along the cables generated by actuators. In particular, for the 2-DoF case we have $m=2$ with two cables. Due to the specific structure, these two tensions are one-directional, i.e.,

\begin{equation}
\label{constraint:direction}
u_i\ge 0 \quad (i=1,2).
\end{equation}

The gravitational and elastic potential energy functions can be derived according to the geometric deformation under the uniformity assumptions of the materials. To make the paper self-contained, the details on the modelling of the potential energy functions are given in Appendix.

For the studied case, the input matrix $G(\bfq): \rea^n \to \rea^{n\times 2}$ can be conformally partitioned as

\begin{equation}
\label{Gpart}
G(\bfq) ~=~ \begmat{~ G_1(\bfq) & G_2(\bfq) ~}.
\end{equation}

In the following assumption, some key properties of the matrix $G(\bfq)$ are underlined when modelling the continuum robot.
The input matrix $G(\bfq)$ of the continuum robot model (<ref>) – or equivalently in (<ref>) – satisfies

(1a) $G(\bfq)$ is state-dependent and $C^1$-continuous.

(1b) $G_1(\textbf{0}_n) = - G_1(\textbf{0}_n)$.

(1c) $|G_1(\bfq) + G_2(\bfq)| \neq 0$ for $\bfq\in {\calx} \backslash \{ \textbf{0}_n\}$.

Among the above three items, the state-dependency of the input matrix is a key feature of the proposed model, which is instrumental to show tunability of open-loop stiffness of tendon-driven continuum robots. We will give more details about the input matrix in the sequel of the paper. The second item means that at the open-loop equilibrium (i.e. $\bfq = \textbf{0}_n$), the tensions from two cables have the same value but with opposite directions.

Let us now consider a single link in the zoomed-in subfigure in Fig. <ref>. If the forces along the cables are assumed lossless, their directions are nonparallel to the centroid of the continuum manipulator. The torques imposed on the first approximate link are given by $u_1 \mathcal{L}_1(q_1)$ and $u_2 \call_2(q_1)$ with $\mathcal{L}_1, \mathcal{L}_2$ the lever's fulcrums, which are nonlinear functions of the configuration $q_1$. From some basic geometric relations, it satisfies $\call_1(0)= \call_2(0)$. This illustrates the rationality of Assumption <ref>.

§.§ Problem Set

In this paper, we study how to design a feedback controller that is capable of regulating the continuum robot deformation and achieving the variable stiffness capability. To be precise, the closed loop complies with the input constraint (<ref>) and achieves the following aims:

A1: In the absence of the external torque (i.e. $\tau_{\tt ext} = \mathbf{0}_n$), it achieves the asymptotically accurate regulation of the position, that is

\begin{equation}
\lim_{t\to\infty} \bfq(t)  = \bfq_\star,
\end{equation}

with the desired configuration $\bfq_\star \in \calx$.

A2: We are able to control the stiffness at the closed-loop equilibrium concurrently.

§ OPEN-LOOP STIFFENING

It is widely recognised that in tendon-driven continuum robots, the antagonism mechanism is a popular method for stiffening. By arranging a pair of cables on two sides of the robot and adjusting their tensions simultaneously, the robot body compresses or expands, and generates a reaction that counteracts the tension. As a consequence, the stiffness of the robot changes.

In this section, it is shown via an intuitive example that there is a redundant degree of freedom of input in the proposed model, which provides the possibility to regulate the stiffness in the open-loop system. This property is instrumental for the controller design to regulate position and stiffness simultaneously. To this end, the following assumption is imposed.

The axial stiffness $K_{\tt A}$ at the end-effector is infinite, i.e. $K_{\tt A} \to \infty$. In other words, the continuum manipulator is axially inextensible.

Now consider the case of the open-loop system with a pair of identical constant inputs

u_1= u_2= μ.
This balance will keep the configuration variable $\bfq$ at the open-loop equilibrium $\bfq_\star=\mathbf{0}_n$ for any feasible $\mu\ge 0$. Specially, when $\mu =0$, the manipulator is in a slack state where its stiffness corresponds to its inherent properties determined by the materials and mechanical structures. Intuitively, as $\mu$ increases, the manipulator progressively transitions towards a state of higher rigidity or inflexibility. The proposed dynamical model should be capable of interpreting the above phenomenon – physical common senses tell us that a larger value of $\mu >0$ implies a larger transverse stiffness.

In order to study the stiffness at the end-effector, we need the Jacobian $J(\bfq) \in \rea^{2\times n}$ from the contact force $\bff_{\tt ext}\in \rea^2$, i.e., satisfying

\begin{equation}
\tau_{\tt ext} = J^\top(\bfq)  \bff_{\tt ext},
\end{equation}

in which $\tau_{\tt ext} \in \rea^n$ is the external torque vector acting on each link.

In the following proposition, we aim to demonstrate the property of stiffness tunability by changing the value $\mu$ in the proposed antagonistic tendon-driven model.

Consider the antagonistic tendon-driven model (<ref>) for continuum robots, with the constant inputs (<ref>) under Assumption <ref>, and assume that the stiffness matrix is diagonal. If the following hypotheses are satisfied:

H1: For $j=1,2$

\begin{equation}
\label{eq:diag}
{\partial G_{n,j} \over \partial q_n} (\bfq) < 0, ~~ {\partial G_{n,j} \over \partial q_k} (\bfq)  = 0, ~k \in \caln \backslash\{n\}
\end{equation}

in some small neighbourhood of the origin.

H2: The forward kinematics (mapping from the configuration $\bfq \in \calx$ to the end-effector Cartesian coordinate $\bfx \in \rea^2$) is a locally injective immersion.

Then, the input (<ref>) guarantees the origin $\mathbf{0}_n$ an equilibrium in the absence of the external perturbation (i.e. $\tau_{\tt ext}=0$), and a large input value $\mu>0$ implies a larger transverse stiffness $K_{\tt T}$ of the end-effector at this equilibrium.

From the second item of Assumption <ref>, the input term

G(\bfq)\bfu |_{\bfq = \mathbf{0}_n} = \mathbf{0}_n

with (<ref>) guarantees that the origin $\bfq_\star = \mathbf{0}_n$ is an equilibrium in the case of $\tau_{\tt ext} = \mathbf{0}$.

The Cartesian coordinate $\bfx \in \rea^2$ of the end-effector can be uniquely determined by the configuration $\bfq$ as

\begin{equation}
\bfx= T(\bfq)
\end{equation}

for some smooth function $T: \rea^n \to \rea^2$, with the open-loop equilibrium $\bfx_\star:= T(\bfq_\star)$. Note that the function $T$ depends on the coordinate selection. Without loss of generality, we assume $\bfx_\star = \mathbf{0}_2$ and the local coordinate of $\bfx :=\col(x_1, x_2)$ are selected as the tangential and the axial directions of the $n$-th link.

For a non-zero external torque $\tau_{\tt ext}$, let us denote the shifted equilibrium as $(\bar\bfq, \mathbf{0}_n)$ in the presence of the external perturbation, with the corresponding end-effector coordinate
\bar\bfx = T(\bar\bfq).

In order to study the transverse stiffness, we assume that the external force $\bff_{\tt ext}$ is only acting to the $n$-th link [29]. Substituting it into the port-Hamiltonian model (<ref>), it satisfies the following equation at the shifted equilibrium

\begin{equation}
\label{eq:nabla_U_f}
\begin{aligned}
\nabla U(\bar\bfq) & ~=~ [G_1(\bar\bfq) + G_2(\bar\bfq)] \mu + J^\top(\bar\bfq)\bff_{\tt ext}
\\
\bff_{\tt ext} & ~=~ K_{\tt C}(\bfx_\star - \bar\bfx),
\end{aligned}
\end{equation}

in which $K_{\tt C} \in \rea^{2\times 2}$ is the stiffness matrix with the partition

\begin{aligned}
K_{\tt C} &:= \diag(K_{\tt T}, K_{\tt A})
\\
J(\bfq) & = \begmat{J_1(\bfq) \\ J_2(\bfq)}
\end{aligned}

with the axial stiffness $K_{\tt A}$ and the transverse stiffness $K_{\tt T}$.

For the particular coordinate selection as mentioned above, the Jacobian matrix $J_1(\bfq)$ is in the form

\begin{equation}
J_1 = \begmat{ 0 & \ldots & 0 & l_n},
\end{equation}

in which $l_n>0$ represents the distance from the contact point to the centre of the $n$-th link.

From Assumption <ref>, the continuum manipulator is inextensible along the axial direction, which implies
\bar x_2 - {x_\star}_{2} =0.
Hence, we have

\begin{equation}
\label{nabla_U:K_T}
\nabla U(\bar\bfq) = [G_1(\bar\bfq) + G_2(\bar\bfq)]\mu + K_{\tt T} J_1^\top \left[\bar x_1 - x_{\star 1}\right].
\end{equation}

For convenience of presentation and analysis, we define a function $f_\mu : \rea^n \to \rea^n$ as

\begin{equation}
f_{\mu}(\bfq) := \nabla U(\bfq) - [G_1(\bfq) + G_2(\bfq)]\mu,
\end{equation}

which is a function parameterised by $\mu \ge 0$. Invoking that $\bfq_\star = \mathbf{0}_n$ is the open-loop equilibrium, we have

\begin{equation}
f_\mu(\bfq_\star) =0,
\end{equation}

and thus

\begin{equation}
\label{eq:f_mu'_star}
f_\mu(\bar\bfq) - f_\mu(\bfq_\star) = K_{\tt T} J_1^\top(\bar\bfq)  [\bar x_1 - x_{\star 1}].
\end{equation}

Invoking the local injectivity hypothesis H2, there exists a left inverse function $T^{\tt L}$ – which is defined locally – of the function $T$ such that

\begin{equation}
\bfq = T^{\tt L} (T(\bfq))
\end{equation}

in a small neighborhood of $\bfq_\star$.

From the equation (<ref>), the transverse stiffness $K_{\tt T}$ at the open-loop equilibrium $\bfq_\star$ is defined by taking $\bar\bfq \to \bfq_\star$, i.e.

\begin{equation}
\label{eq:K_T}
\begin{aligned}
K_{\tt T} &~=~   {J_1 (\bfq_\star) \over |J_1(\bfq_\star)|^2 } \lim_{\bar x_1 \to x_{\star 1}} {  f_\mu(\bar\bfq) - f_\mu(\bfq_\star) \over \bar x_1 - x_{\star1}}
\\
&~=~ {J_1 (\bfq_\star) \over |J_1(\bfq_\star)|^2 } [\nabla f_\mu(\bfq_\star)]^\top \nabla_{x_1} T^{\tt L}(\bfx_\star)
\\
&~=~ {J_1 \over |J_1|^2 }\big[\nabla^2 U - \mu\big( \nabla G_1  + \nabla G_2 \big)^\top\big]
\nabla_{x_1} T^{\tt L} \bigg|_{\bfq_\star}.
\end{aligned}
\end{equation}

The hypothesis H1 guarantees that the $(n,n)$-element of $\nabla G_j(\bfq_\star)$ for $j=1,2$ is negative. On the other hand, from the local coordinate selection, the variation of $x_1$ implies that $q_n$ will also change accordingly, as consequence the last element of $\nabla_{x_1} T^{\tt L}$ is non-zero (indeed positive). It is straightforward to see that $K_{\tt T}$ is increasing by selecting a larger $\mu>0$.

The above calculation shows the underlying mechanism on tunability of stiffness in the open loop. Later on, we will illustrate that the tendon difference $(u_1 - u_2)$ provides another degree to regulate the robot configuration.

The hypothesis H1 is a technical assumption, which means $G_{n,j}$ only depends on the state $q_n$ rather than other configuration variables. It is used to simplify the presentation and analysis. Indeed, it is unnecessary, and with only $J_1(\nabla G_1 + \nabla G_2) \nabla_{x_1} T^{\tt L} \neq 0$, we are able to show the ability to tune the stiffness. The assumption H2 means that we are able to find a unique inverse kinematics in a small neighborhood of a given configuration $\bfq_\star$. Though it is generally not true to guarantee the existence of a global inverse, it is possible to achieve it locally.

§ CONTROL DESIGN

In this section, we will study how to design a state feedback law to regulate position and stiffness simultaneously using the proposed dynamical model in Section <ref>.

To facilitate the controller design, we additionally assume the following for the input matrix $G(\bfq)$ in terms of the geometric constraints.

The matrix $G(\bfq)$ given by (<ref>) is reparameterised as

\begin{equation}
\begin{aligned}
G_1(\bfq) & ~=~  \bfg_1(\bfq) + \bfg_0 \\
G_2(\bfq) & ~=~  \bfg_1(\bfq) - \bfg_0
\end{aligned}
\end{equation}

with a constant vector $\bfg_0 \in \rea^n$ and a $C^1$-continuous function $\bfg_1:\rea^n \to \rea^n$ satisfying the following:

(3a) $\bfg_1(\bfq)$ is a smooth odd function;

(3b) $\bfg_1(\bfq)$ is full column rank for $\bfq\in \calx \slash \{\mathbf{0}_n\}$;

The constant vector $\bfg_0$ and the vector field $\bfg_1(\bfq)$ can be re-parameterised as

\begin{equation}
\label{eq:g0gq}
\bfg_0 = g_0 \mathbf{1}_n , \quad \bfg_1(\bfq) = g_1(\bfq) \mathbf{1}_n,
\end{equation}

and $g_0+ g_1(\bfq) \neq 0$ for all $\bfq$.

Clearly, the above is compatible with Assumption <ref>, and the function $\bfg_1(\bfq)$ is related to the open-loop stiffness tunability outlined in Proposition <ref>.

§.§ Assignable Equilibria

For underactuated mechanical systems, it is crucial to identify the set of assignable equilibria, also known as achievable or feasible equilibria. Although extensively explored, tendon-driven robots face a significant obstacle in the form of the one-directional input constraint (<ref>).

In order to facilitate the control design, we make the following input transformation:

\begin{equation}
\label{trans:input}
\tau = T_{\tt u} \bfu, \quad  T_{\tt u}:=\begmat{1 & -1 \\ 0 & 1}
\end{equation}

with new input control $\tau= \col(\tau_1,\; \tau_2) \in \rea^2$. For $\tau_1$, there is no sign constraint for its value; the other input channel verifies $\tau_2 \ge 0$, and we define the admissible input set as

\begin{equation}
\label{E_tau}
\cale_\tau := \{\tau \in \rea^2: \tau_1 \in \rea, ~\tau_2 \ge 0\}.
\end{equation}

Invoking the intuitive idea in Section <ref>, we may use these two inputs $\tau_1,\tau_2$ to regulate the position and stiffness concurrently. For convenience, we define the new input matrix as

\begin{equation}
G_\tau(\bfq) = \begmat{\rho_1(\bfq) & \rho_2(\bfq)} := G(\bfq)T_{\tt u}^{-1}
\end{equation}

\begin{aligned}
\rho_1(\bfq) &  ~=~ \bfg_0 + \bfg_1(\bfq)
\\
\rho_2(\bfq) &  ~=~ 2\bfg_1(\bfq),
\end{aligned}

and then the controlled model becomes

\begin{equation}
\label{model:pH2}
\begmat{ ~\dot \bfq ~\\ ~\dot \bfp~}
\begmat{0_{n\times n} & I_n \\ - I_n & -D(\bfq)} \begmat{\nabla	_{\bfq} H \\ \nabla_{\bfp} H }+ \begmat{~\mathbf{0}_{n} ~ \\ ~G_\tau(\bfq) \tau ~}
\end{equation}

in the absence of the external perturbation $\tau_{\tt ext}$.

Indeed, the real constraint for $\tau_1$ should be $\tau_1 \ge -\tau_2$ rather than $\tau_1\in \rea$ in order to guarantee the constraint $u_1,u_2\ge0$. Since we are able to set the value of $\tau_2$ arbitrarily, we consider the admissible input set $\cale_\tau$ defined above for convenience in the subsequent analysis.

According to [23] and invoking the full-rankness of $T_u$, if there were not input constraints, the assignable equilibrium set  would be given by
\{ \bfq\in \rea^n : G(\bfq)^\bot \nabla U (\bfq) =0 \}.
Clearly, this does not hold true for our case, since the feasible solution cannot be guaranteed to live within the set $\cale_\tau$ rather than $\tau\in \rea^2$. To address this point, in the following proposition we present the assignable equilibria set for the studied case with constrained inputs.

(Assignable Equilibria)
Consider the model (<ref>) and the input transformation (<ref>) with the input constraint $\cale_\tau$ in (<ref>). All the assignable equilibria are given by the set $\cale_q \cap \calx$, with the definition

\begin{equation*}
\cale_q := \left\{ \bfq\in \rea^n \Bigg|
\begin{aligned}
(\bfg_1^\bot \bfg_0)^\bot \bfg_1^\bot \nabla U =0
\\
\bfg_1^\top \nabla U - \bfg_1^\top(\bfg_0+\bfg_1) (\bfg_1^\bot \bfg_0)^\dagger \bfg_1^\bot \nabla U \ge 0
\end{aligned}
\right\}.
\end{equation*}

In terms of Assumption (<ref>b), we can always find a left annihilator $\bfg_1^\bot(\bfq) \in \rea^{(n-1)\times n}$, which is full rank for all $\bfq \in \calx\slash\{\mathbf{0}_n\}$. For an equilibrium $(\bfq,{\bf 0})$, there should exist $\tau_1$ and $\tau_2$ satisfying

\begin{equation}
\label{id:eq1}
\begin{aligned}
\nabla U (\bfq) & ~= ~ G_\tau(\bfq)\tau
\\
& ~= ~
[\bfg_0 + \bfg_1(\bfq)]\tau_1 + 2 \bfg_1(\bfq)\tau_2.
\end{aligned}
\end{equation}

Considering the full-rankness of the square matrix

\begmat{~\bfg_1^\bot (\bfq) ~\\ ~\bfg_1^\top (\bfq)~} \in \rea^{n\times n}

for $\bfq \in  \calx\slash\{\mathbf{0}_n\}$, we have

\begin{equation}
\eqref{id:eq1} ~\iff ~ \begmat{\bfg_1^\bot  \\ \bfg_1^\top } \nabla U
\begmat{ \bfg_1^\bot  \bfg_0 \tau_1 \\ \bfg_1^\top(\bfg_0 + \bfg_1) \tau + 2|\bfg_1|^2\tau_2 }.
\end{equation}

Its solvability relies on finding all the points $\bfq\in \calx \subset \rea^n$ satisfying

\begin{align}
\bfg_1^\bot \nabla U & = \bfg_1^\bot \bfg_0 \tau_1 \label{c1}
\\
\bfg_1^\top \nabla U & = \bfg_1^\top(\bfg_0 + \bfg_1) \tau_1 + 2|\bfg_1|^2\tau_2  \label{c2}
\end{align}

at the same time under the constraint $\tau \in \cale_\tau$. Clearly, all the feasible equilibria satisfying (<ref>) live in the set $\{\bfq\in \rea^n:(\bfg_1^\bot(\bfq) \bfg_0)^\bot \bfg_1^\bot(\bfq) \nabla U(\bfq) =0  \}$, and the corresponding input $\tau_1$ is given by

\begin{equation}
\label{solution:tau1}
\tau_1  = (\bfg_1^\bot \bfg_0)^\dagger \bfg_1^\bot \nabla U.
\end{equation}

On the other hand, Assumption (<ref>b) imposes the constraint $|\bfg_1|^2>0$, thus (<ref>) admits a positive solution to $\tau_2 \ge 0$ if and only if

\begin{equation}
\bfg_1^\top \nabla U  -  \bfg_1^\top(\bfg_0+\bfg_1) \tau_1 \ge 0.
\end{equation}

By inserting (<ref>) into the above equation, we complete the proof.

After imposing Assumption (<ref>c) to the input matrix $G(\bfq)$, we are interested in a class of particular equilibria, and we call them the homogeneous equilibria that are characterised by the set

\begin{equation}
\cale_\theta :=\{ \bfq\in \rea^n :  \bfq_\star :=  \col(\theta, \ldots, \theta) \}
\end{equation}

for some constant $\theta$. This definition is tailored for the proposed continuum robot model under the assumptions in the paper.

In the following, we show all homogeneous equilibria belong to the assignable equilibrium set $\cale_q$ as defined in Proposition <ref>.

Consider the model (<ref>) of the continuum robot under Assumptions <ref>-<ref>. Then, all homogeneous equilibria are assignable, i.e. $\cale_\theta \subset \cale_q$.

For the case with $\theta = 0$, since $g_1(\mathbf{0}_{n}) =0$, the equilibrium $\bfq_\star= \theta \mathbf{1}_n$ makes the equation (<ref>) solvable with $\tau_1 = 0$ and any $\tau_2  \ge 0$.

For the case with $\theta \neq 0$ and any fixed $\tau_1\ge 0$, the determination of the set $\cale_q$ is equivalent to solving (<ref>), which can be written as

\begin{equation}\label{pde:tau_n}
\begin{aligned}
\nabla U(\bfq) & =  [\bfg_0 + \bfg_1(\bfq)]\tau_1 + 2 \bfg_1(\bfq)\tau_2
\\
& = G_{\tt N}\tau_{\tt N}
\end{aligned}
\end{equation}

with the definitions

\begin{equation}
\begin{aligned}
G_{\tt N} & ~:=~ \mathbf{1}_n
\\
\tau_{\tt N}(\tau) & ~:=~ [g_0+g_1(\bfq)]\tau_1 + 2g_1(\bfq) \tau_2.
\end{aligned}
\end{equation}

For any fixed $\tau_2\ge 0$, invoking (<ref>) from Assumption <ref>, the mapping $\tau_1 \to \tau_{\tt N}$ is a diffeomorphism from $\rea \to \rea$. It implies that there is no constraint for $\tau_{\tt N}$. As a consequence, the PDE (<ref>) becomes

\begin{equation}
\label{GN_nabla_U}
G_{\tt N}^\bot \nabla U(\bfq_\star) =0
\end{equation}

at the desired equilibrium $\bfq_\star$.

A feasible full-rank annihilator of $G_{\tt N}$ is given by

\begin{equation}
G_{\tt N}^\bot = \begmat{1& -1 & 0& \ldots & 0 \\ 0 & 1& -1 & \ldots &0 \\  && \ddots&\ddots  \\ 0 & \ldots& 0& 1& -1} \in \rea^{(n-1)\times n},
\end{equation}

and the Jacobian $\nabla U $ at the desired equilibrium $\bfq_\star$ is in the form

\begin{equation}
\nabla U (\bfq_\star) = \alpha_1 \underbrace{\begmat{~\sin(q_\Sigma)~ \\ \vdots\\ ~\sin(q_\Sigma)~}}_{ \sin(n\theta) \mathbf{1}_n }
+ \alpha_2 \underbrace{\begmat{~\bfq_{\star,1} ~\\ ~\vdots~ \\ ~\bfq_{\star,n} ~}}_{\theta \mathbf{1}_n}
\end{equation}

q_\Sigma := \sum_{i\in \caln} \bfq_i.

It is straightforward to verify that (<ref>) holds true for any $\tau_2\ge 0$ with a homogeneous equilibrium $\bfq_\star$. Since there is no constraint for the input variable $\tau_1$, the equilibrium for this case is also assignable under the constraint (<ref>). We complete the proof.

In the sequel of the paper, our focus will be on control design aimed at regulating certain homogeneous equilibria that have been demonstrated to be assignable within the proposed class of models for continuum robots.

§.§ Simultaneous Position-and-Stiffness Control

We now aim at stabilising an arbitrary homogeneous equilibrium $\bfq_\star$ in the subset of $\cale_\theta$ with a tunable stiffness of the closed loop.

Towards the end, we will employ the passivity-based control (PBC) method since it has a clear energy interpretation and simplifies both modelling and controller design. This makes it suitable for continuum robotics to preserve the system compliance [15].

Our basic idea is to fix $\tau_2$ at some constant value $\tau_2^{\star}\ge 0$. We then utilise the input $\tau_1$ to achieve potential energy shaping for the regulation task. Compared to the more general approach of interconnection and damping assignment (IDA) PBC [24], on one hand, potential energy shaping may provide a simpler controller form, and on the other hand, as pointed out in [18] changing the inertia is prone to fail in practice – albeit being theoretically sound with additional degrees of freedom.

For a given input $\tau_2 = \tau_2^\star \ge 0$, the actuation into the dynamics is given by

\begin{equation}
\label{tau_N}
\begin{aligned}
G_\tau(\bfq) \tau & ~= ~\rho_1(\bfq) \tau_1 + \rho_2(\bfq) \tau_2^\star
\\
&:=~ G_{\tt N} \tau_{\tt N}(\col(\tau_1, \tau_2^\star)),
\end{aligned}
\end{equation}

with a new function $\tau_{\tt N}: \rea^2 \to \rea$. From Assumption <ref>, the vector field $\rho_1(\bfq) \neq 0$ for all $\bfq \in \calx$. Now the design target becomes using the control input $\tau_1$ (with a fixed $\tau_2^\star$) to shape the potential energy function $U(\bfq)$ into a new shaped potential energy function $U_{\tt d}(\bfq)$. To this end, we need to solve the PDE [23]

\begin{equation}
\label{pde:shaping}
G_{\tt N}^\bot \big[ \nabla U(\bfq) - \nabla U_{\tt d}(\bfq) \big] = 0.
\end{equation}

Note that the solution must adhere to the constraints

\begin{align}
\nabla U_{\tt d}(\bfq_\star) & ~=~ 0
\\
\nabla^2 U_{\tt d}(\bfq_\star) & ~\succ ~ 0,
\end{align}

in order to make the desired configuration $\bfq_\star$ an asymptotically stable equilibrium.

We are now in the position to propose the controller for simultaneous control of position and stiffness.

Consider the continuum robotic model (<ref>), (<ref>) with the constraint (<ref>) satisfying Assumptions <ref>-<ref>. The controller

\begin{equation}
\label{u:control}
u = T_u^{-1} \tau
\end{equation}

with the transformed input

\begin{equation}
\label{tau+}
\tau =  \tau_{\tt es}+ \tau_{\tt da}  + \tau_{\tt st}
\end{equation}

and the terms

\begin{equation}
\label{tau:3}
\begin{aligned}
\tau_{\tt st} & = \begmat{
-{ 2g_1(\bfq) \over g_0+g_1(\bfq)}  \\  1
\\
\tau_{\tt es} & = \begmat{
{1\over g_0+g_1(\bfq)} G_{\tt N}^\dagger (\nabla U_{\tt d} - \nabla U)  \\  0
\\
\tau_{\tt da} & = \begmat{ - {1\over g_0+g_1(\bfq)} G_{\tt N}^\top K_{\tt d} M^{-1}(\bfq)\bfp  \\ 0},
\end{aligned}
\end{equation}

where $G_{\tt N} = \mathbf{1}_n$, $\tau_2^\star >0$, $K_{\tt d} \succ 0$ is gain matrix, and the desired potential energy function is given by

\begin{equation}
\label{Ud}
\begin{aligned}
U_{\tt d}(\bfq) &~=~ - \gamma \cos(q_\Sigma - q_{\Sigma}^\star) + {\alpha_2 \over 2}|\bfq - \bfq_\star|^2
\\
q_{\Sigma}^\star &~=~ \sum_{i \in \caln} \bfq_{\star,i}
\end{aligned}
\end{equation}

and the gain $\gamma>0$, achieves the following closed-loop properties:

P1: (Position regulation in free motion) If the external force $\tau_{\tt ext} =0$ and $\gamma <\alpha_2 $, then the desired equilibrium point $\bfq_\star$ is globally asymptotically stable (GAS) with

\begin{equation}
\lim_{t\to + \infty} \bfq(t) = \bfq_\star.
\end{equation}

P2: (Compliant behavior) The overall closed-loop stiffness (i.e., from the external torque $\tau_{\tt ext} \in \rea^n$ to the configuration $\bfq \in \rea^n$) is

\begin{equation}
\label{K_O}
K_{\tt O} = \gamma \mathbf{1}_{n\times n} + \alpha_2 I_n,
\end{equation}

where $\mathbf{1}_{n\times n}$ is an $n\times n$ matrix of ones.

First, it is straightforward to verify that the vector $\tau_{\tt st}$ in (<ref>) is in the null space of $\tau_{\tt N}(\tau)$ for any $\tau_2^\star$, i.e.,

\begin{equation}
\label{tntst=0}
\tau_{\tt N}(\tau_{\tt st}) = 0, \quad \forall \tau_2^\star \in \rea_{\ge 0}.
\end{equation}

It means that the first term $\tau_{\tt st}$ does not change the closed-loop dynamics.

Now, let us study the effect of the potential energy shaping term $\tau_{\tt es}$. The Jacobian of the desired potential energy function $U_{\tt d}$ is given by

\begin{equation}
\nabla U_{\tt d}(\bfq) = \gamma \sin(q_\Sigma - q_\Sigma^\star) \mathbf{1}_n + \alpha_2(\bfq - \bfq_\star).
\end{equation}

It satisfies the following:

\begin{align}
\label{Jacobian:Ud}
\nabla U_{\tt d} (\bfq_\star) & ~=~0
\\
\nabla^2 U_{\tt d} (\bfq) & ~=~  \gamma \cos(q_\Sigma - q_\Sigma^\star)\mathbf{1}_{n\times n}
+ \alpha_2 I_n \succ 0,
\end{align}

in which the  second equality holds true for all $\bfq\in \rea^n$ by noting that the eigenvalues of the symmetric matrix $\nabla^2 U_{\tt d}$ are given by

\{\underbrace{\alpha_2, \ldots, \alpha_2}_{n-1}, \alpha_2 + \gamma  |\cos(q_\Sigma - q_\Sigma^\star)|\}

with all elements positive from the condition $\gamma < \alpha_2$ in P1. This implies that the desired potential energy function $U_{\tt d}$ is convex and achieves its global minimum at $\bfq_\star$.

For the function $U_{\tt d}$, we have

\begin{equation*}%\small
\begin{aligned}
& G_{\tt N}^\bot [\nabla U - \nabla U_{\tt d}]
\\
= & G_{\tt N}^\bot [ \alpha_1 \sin(q_\Sigma)\mathbf{1}_n + \alpha_2 \bfq - \gamma \sin(\bfq_\Sigma - \bfq_\Sigma^\star) \mathbf{1}_n - \alpha_2 (\bfq -
\bfq_\star ) ]
\\
= & G_{\tt N}^\bot [ \alpha_1 \sin(q_\Sigma)\mathbf{1}_n  - \gamma \sin(q_\Sigma - q_\Sigma^\star) \mathbf{1}_n + \alpha_2 \bfq_\star  ]
\\
=& 0,
\end{aligned}
\end{equation*}

where in the last equation we have used the fact $\bfq_\star \in \cale_\theta$, so that  the PDE (<ref>) is verified. Together with (<ref>), the controller (<ref>) makes the closed-loop dynamics have the form

\begin{equation}
\label{cl_dyn:ext}
\begmat{ ~\dot \bfq ~\\ ~\dot \bfp~}
\begmat{0_{n \times n} & I_n \\ - I_n & -\mathsf{D}} \begmat{ \nabla_\bfq H_{\tt d} \\  \nabla_\bfp H_{\tt d}} + \begmat{~\mathbf{0}_n ~ \\ ~ \tau_{\tt ext}~},
\end{equation}

\begin{equation}
\begin{aligned}
H_{\tt d}(\bfq,\bfp) & := {1\over2} \bfp^\top M^{-1}(\bfq) \bfp + U_{\tt d}(\bfq)
\\
\mathsf{D}(\bfq) &:= D(\bfq) + G_{\tt N} K_{\tt d} G_{\tt N}^\top \succ 0.
\end{aligned}
\end{equation}

For free motion (i.e. $\tau_{\tt ext} = 0$), following the standard Lyapunov analysis we have

\begin{equation}
\dot {H}_{\tt d} \le -(\nabla H_{\tt d})^\top  \mathsf{D}(\bfq) \nabla H_{\tt d} \le 0.
\end{equation}

For the system (<ref>), the set $\{(\bfq, \bfp): (\nabla H_{\tt d})^\top  \mathsf{D}(\bfq) \nabla H_{\tt d} =0\}$ contains only a single point $(\bfq_\star, \mathbf{0}_n)$. From the LaSalle's invariance principle <cit.>, we are able to show the global asymptotic stability of the desired equilibrium $(\bfq_\star, \mathbf{0}_n)$. Hence, we have proven the property P1.

The next step is to verify the stiffness property P2 in a small neighborhood $B_\varepsilon(\bfq_\star)$ of $\bfq_\star$ with a sufficiently small $\varepsilon>0$. For a constant external force $\tau_{\tt ext}$, the shifted equilibrium point $(\bar\bfq,\mathbf{0}_n)$ should satisfy

\begin{equation}
- \nabla_\bfq U_{\tt d} + \tau_{\tt ext} =0,
\end{equation}

or equivalently

\begin{equation}
\label{phi:q}
\phi(\bar\bfq):=  \gamma \sin(\bar q_\Sigma - q_\Sigma^\star) \mathbf{1}_n + \alpha_2 (\bar\bfq - \bfq_\star) = \tau_{\tt ext},
\end{equation}

with $\bar q_\Sigma:=\sum_{i \in \caln} \bar \bfq_i$. Note that $\nabla \phi = \nabla^2 U_{\tt d} \succ 0$, which means that
\phi: \bfq\mapsto \tau_{\tt ext}

is a (locally) injective immersion. Hence, in the small neighborhood $B_\varepsilon(\bfq_\star)$ of $\bfq_\star$, there is a unique solution $\bar\bfq$ to (<ref>) for a given $\tau_{\tt ext}$.

We show that the shifted equilibrium $(\bar\bfq, \mathbf{0}_n )$ is asymptotically stable by considering the Lyapunov function

\begin{equation}
V(\bfq,\bfp) = H_{\tt d}(\bfq,\bfp) - \bfq^\top \tau_{\tt ext}.
\end{equation}

From the above analysis, it is clear that

\begin{equation}
\begin{aligned}
\nabla V(\bar\bfq,\mathbf{0}) & = \begmat{\phi(\bar\bfq) - \tau_{\tt ext} \\ \mathbf{0}} = \mathbf{0}
\\
\nabla^2 V(\bar\bfq,\mathbf{0}) &\succ 0.
\end{aligned}
\end{equation}

Hence, $V$ qualifies as a Lyapunov function. Its time derivative along the system trajectory is given by

\begin{equation}
\begin{aligned}
\dot V &~=~ - \|\nabla_\bfp H_{\tt d}\|_{\mathsf{D}}^2 - (\nabla_\bfp H_{\tt d})^\top \tau_{\tt ext} - \dot \bfq^\top \tau_{\tt ext}
\\
&~= ~ - \|\nabla_\bfp H_{\tt d}\|_{\mathsf{D}}^2.
\end{aligned}
\end{equation}

It yields the Lyapunov stability of the closed-loop dynamics (<ref>) in the presence of a constant external torque $\tau_{\tt ext}$, with all the system states bounded. On the other hand, the set

\cale_{u}:= \{(\bfq,\bfp):\|\nabla_\bfp H_{\tt d}(\bfq)\| =0 \}

only contains a single isolated equilibrium $(\bar\bfq,\mathbf{0})$. According to LaSalle's invariance principle, $(\bar\bfq,\mathbf{0})$ is an asymptotically stable equilibrium, in which $\bar\bfq$ depends on the (arbitrary) constant torque $\tau_{\tt ext}$ – in terms of the unique solution of the algebraic equation (<ref>).

The overall stiffness $K_{\tt O}$ is defined by $\tau_{\tt ext} =  K_{\tt O}(\bar\bfq-\bfq_\star)$. Substituting it into (<ref>), we have

\begin{equation}
\gamma \sin(\mathbf{1}_n^\top \bar\bfq - q_\Sigma^\star) \mathbf{1}_n + \alpha_2 (\bar\bfq - \bfq_\star) =  K_{\tt O}(\bar\bfq-\bfq_\star).
\end{equation}

The stiffenss at the desired equilibrium $\bfq_\star$ is calculated from any direction of the limit $\bar\bfq\to \bfq_\star$, thus obtaining

\begin{equation}
\begin{aligned}
K_{\tt O}\Big|_{\bfq_\star} & ~=~ {\partial \phi\over \partial \bfq}(\bfq_\star)
\\
& ~=~ \lim_{\bar\bfq\to \bfq_\star} \gamma \cos(\mathbf{1}_n^\top \bar\bfq - q_\Sigma^\star) \mathbf{1}_{n\times n} + \alpha_2 I_n
\\
& ~=~ \gamma \mathbf{1}_{n\times n} + \alpha_2 I_n.
\end{aligned}
\end{equation}

This verifies the property P2, and we complete the proof.

The above shows that the proposed controller (<ref>)-(<ref>) can achieve the end-effector position regulation with the closed-loop stiffness $K_{\tt O}$ given by (<ref>). It implies our ability to set a prescribed stiffness by selecting the control gain $\gamma>0$ properly. More discussions about the control law will be provided in the next section.

§ DISCUSSIONS

The following remarks about the proposed controller are given in order.

In the property P2, we study the overall stiffness – from the external torque vector $\tau_{\tt ext}$ to the configuration $\bfq \in \calx$, rather than the transverse stiffness at the end-effector. Consider the external force $f_{\tt ext}\in \rea$ at the end-effector along the transverse direction of the $n$-link with the Jacobian $J = [0,\ldots, 0, \ell]$. With a small force $f_{\tt ext}$, the coordinate of the end-effector would shift from

\left(
\ell\left(\sum_{k \in \caln}\sin(k\theta_\star)\right),
\ell  \left(\sum_{k \in \caln}\cos(k\theta_\star)\right)
\right)

(F_x( f_{\tt ext}), F_y(f_{\tt ext}))

\begin{pmatrix}  ~F_x ~\\ ~F_y~ \end{pmatrix}
\begin{pmatrix}
\ell \sin(\beta) + \ell\sum_{k \in \caln\backslash \{n\}}\sin(k\theta_\star + \gamma\ell f_{\tt ext})
\\
\ell \cos(\beta) + \ell\sum_{k \in \caln\backslash \{n\}}\cos(k\theta_\star + \gamma\ell f_{\tt ext})
\end{pmatrix}

and $\beta := n\theta_\star + (\gamma + \alpha_2) \ell f_{\tt ext}$. Hence, the transverse stiffness is given by

\begin{equation*}
\begin{aligned}
K_{\tt T} = \lim_{ f_{\tt ext} \to 0 }
{ \sqrt{[F_x(f_{\tt ext}) - F_x(0)]^2 + [F_y(f_{\tt ext}) - F_y(0)]^2} \over  f_{\tt ext}} .
\end{aligned}
\end{equation*}

As a result, we have

\begin{equation}
\label{K_T:prop}
\boxed{
~ K_{\tt T} ~\propto~  \kappa_1\gamma + \kappa_2 ~
\end{equation}

with some non-zero constants $\kappa_1, \kappa_2$ for $\theta_\star \neq 0$. It means that, for a given desired equilibrium $\bfq_\star \in \calx\backslash\{\mathbf{0}\}$, the transverse stiffness is affine in the gain $\gamma$, thus providing a way to tune the closed-loop stiffness linearly.

The proposed controller can be roughly viewed as a nonlinear PD controller. The first term $\tau_{\tt st}$ is used to compensate the “anisotropy” in the input matrix $G(\bfq)$ due to its state-dependency property; the potential energy shaping term $\tau_{\tt es}$ and the damping injection term $\tau_{\tt da}$, indeed, play the role of nonlinear PD control. To be precise, the term $\tau_{\tt es}$ is the error between the nonlinear functions of the position $\bfq$ and its desired value $\bfq_\star$; and the term $\tau_{\tt da}$ can be viewed as the negative feedback of velocity errors. This is not surprising, since the original idea of energy shaping has its roots in the pioneering work of Takegaki and Arimoto in robot manipulator control [30], in which they proposed a very well-known “PD + gravity compensation” feedback [25].

To ensure that $\nabla^2 U_{\tt d}(\bfq_\star) \succ 0$, it is necessary to impose the condition $\gamma< \alpha_2$ on the control gains. However, this condition may restrict the range of closed-loop stiffness values within an interval. If this condition is not imposed, it is only possible to guarantee the positive definiteness of $\nabla U_{\tt d}$ in the vicinity of $\bfq_\star$, which would result in local asymptotic stability. We provide some experimental evidence regarding this point in the next section.

Let us now look at the proposed controller (<ref>)-(<ref>). Note that the term
$M^{-1} (\bfq)\bfp$ corresponds to the generalised velocity of $\bfp$. Thus, the controller
depends on only three plant parameters ($\alpha_1, \alpha_2$ and $g_0$) and a nonlinear function $g_1$ – which need to be identified beforehand – and two adaptation gains (i.e. $K_{\tt d}$ and $\gamma$). This means that it is unnecessary to identify all parameters and functions in the plant model. This makes the resulting controller robust vis-à-vis different types of uncertainties.

§ EXPERIMENTAL RESULTS

§.§ Experimental setup

Photo of the entire experimental platform

In this section, the proposed control approach was tested using the OctRobot-I, a continuum robot developed in our lab at the University of Technology Sydney [14]. We considered the planar case with six segments (i.e. $n=6$) with an overall length of 252 mm, and a diameter of approximately 50 mm, which meets the critical assumptions outlined in the paper. Notably, the OctRobot-I has a jamming sheath that can provide an extra degree of freedom for stiffening, though the present paper does not delve into this feature's stiffening capabilities. Additional details of the OctRobot-I can be found in [14].

As shown in Fig. <ref>, the test platform used in the experiments consists of the one-section robot (OctRobot-I), two servo motors (XM430-W350, DYNAMIXEL) with customized aluminium spools, three force sensors (JLBS-M2-10kg), a linear actuator, and an electromagnetic tracking systems (Aurora V3, NDI). The control experiments and data collections were conducted using the software MATLAB™.

The force sensors used in the platform have a response delay and provide noisy measurements when used for real-time control. On the other hand, the servo motors in the platform can provide accurate position information with high accuracy, making it easier to control cable lengths between the servo motors and the actuator unit. Using Hooke's law, it is possible to approximately consider the cable length proportional to the force for each cable, with some coefficients to be identified off-line using some collected data sets. To verify the linear relationship between the cable length and the tension force, as well as to obtain the coefficients, we conducted a group of experiments with different configurations, and recorded the cable lengths and the corresponding forces. Each configuration was repeatedly conducted three times under the identical condition, and all the data were utilized for identification. In Fig. <ref>, we plot the relation between the right cable length $L_2$ and the corresponding force, and the one between the length difference $\Delta L:= L_1 - L_2$ and the force difference $\tau_1 := u_1 - u_2$ of these two cables. The correlation coefficients are 0.9977 and 0.9987, which imply the strong linearity between cable lengths and forces. Thus, it is reasonable to use the cable lengths – driven by motors – as the “real” input signals.

The linearity between forces and lengths: The length $L_1$ vs the force $u_1$ of the right cable; and the length difference $\Delta L:= L_1 - L_2$ vs the force difference $\tau_1 := u_1 - u_2$ (“$\times$” represents test data, and the dash lines are the fitted functions.)

§.§ Open-loop stiffening experiments

In this subsection, we aim to verify the results regarding open-loop stiffening that were presented in Section <ref>. For this purpose, we utilized a linear actuator positioned at the end-effector to generate a small displacement $\delta x > 0$, as illustrated in Fig. <ref>(a). The actuator was attached to the force sensors to measure the external force $f_{\tt ext}$ related to the displacement. By dividing the measured force by the applied displacement, i.e., ${f_{\tt ext} \over \delta x}$, we were able to estimate the transverse stiffness, as long as $\delta x$ was sufficiently small. This procedure allowed us to validate the findings related to open-loop stiffening presented in Section <ref>.

[Experimental setup]
[Experimental results of transverse stiffness]fig/plot_fig1.tex
Experiments for open-loop stiffening

We measured the stiffness values under different open-loop tendon forces $\mu>0$ in the interval $[0,45]$ N. Each experiment was repeatedly conducted three times under the same condition in order to improve reliability. The experimental results are shown in Fig. <ref>(b), where “$\times$” represents the mean values of the calculated stiffness for all $\mu$, and the error bars are $\pm 1$ standard deviation. This clearly verifies the theoretical results in Proposition <ref>. The correlation coefficient between $\mu$ and the stiffness is 0.989, which illustrates the strong linearity – exactly coinciding with the equation (<ref>).

§.§ Closed-loop experiments

In order to apply the proposed real-time control algorithm to the experimental platform, we first conducted the identification procedure to estimate the parameters outlined in the fourth discussing point in Section <ref>. It relies on the fact that at any static configuration $\bfq$ (i.e. with $\bfp =0$) the identity $\nabla_\bfq U(\bfq) = G(\bfq) \bfu$  holds true, and thus $J(\theta) =0$ with the cost function

\begin{align}
& J(\theta,u_1,u_2) \\
& \hspace{.3cm} := \Big|\alpha_1 \sin(n\theta) + \alpha_2\theta - [g_0 +g_1(\theta)] u_1 + [g_0 - g_1(\theta)]u_2\Big|^2,\nonumber
\end{align}

that contains all the quantities to be identified. Since the function $g_1(\cdot)$ is infinite-dimensional, we simply parameterised as $g_0 + g_1(\theta) = c_1 + c_2 \sin(\theta)$ with two constants $\alpha_1>0$ and $\alpha_2<0$, complying with the assumptions on the input matrix. We regulated the continuum robot to different equilibria $\theta^j$ ($j=1,\ldots, w$ with some $w\in \mathbb{N}_+$) by driving the cables, and recorded the corresponding forces $(u_1^j, u_2^j)$.

The identification procedure boils down to solving the optimisation problem

\begin{equation}
\underset{c_1, \alpha_1,\alpha_2>0, c_2<0}{\arg\min} ~\sum_{j\in \{1,\ldots, w\}}  J(\theta^j, u_1^j, u_2^j) .
\end{equation}

We ran the identification experiments to collect data at 15 equilibria points (i.e. $w = 15$), and repeated for six times. Using this data set, the identified parameters were
c_1= 1.2143,  c_2 = -2.9015,  \alpha_1 = 8.6114
and $\alpha_2 = 0.001$.

To evaluate the performance of position control, we first consider a desired configuration $\mathbf{q}_\star = [\theta_\star, \ldots, \theta_\star]^\top$ with $\theta_\star = 5$ deg for the proposed control scheme. We conducted experiments with various values of the gains $\gamma$ and $K_{\tt d}$, as shown in Figs. <ref>-<ref>, respectively. The second row of these figures depicts the configuration variable at the steady-state stage during $[2,8]$ s. It is worth noting that the control inputs $u_i$ ($i=1,2$) are mapped to the cable length $L_i$, as explained in Section <ref>. In all these scenarios, the transient stage lasted for less than 1.5 seconds, and the configuration variable quickly converged to small neigborhoods of the desired angle, demonstrating the high accuracy of the proposed control approach. There were no apparent overshootings in configuration variables. Our results indicate that selecting either a sufficiently small or large $K_{\tt d}$ can negatively affect the control performance during the transient stage. On the other hand, setting a large $\gamma>0$ may lead to chattering due to measurement noise at the steady-state stage, which is well understood as the deleterious effect of high-gain design in the control literature [25].

Photo sequence of the position control with $\theta_\star = 10$ deg and $\gamma=1, K_{\tt d}=0.1$

[$K_{\tt d} = 10, ~ \gamma =1$]FigPlot/Inputs_5_deg_a
[$K_{\tt d} = 5, ~ \gamma =1$]FigPlot/Inputs_5_deg_b
[$K_{\tt d} = 0.1, ~ \gamma =1$]FigPlot/Inputs_5_deg_c
[$K_{\tt d} = 0.01, ~ \gamma =1$]FigPlot/Inputs_5_deg_d

Position control performance for the desired configuration $\theta_\star = 5 ~{\rm deg}$ with different $K_{\tt d}$

[$K_{\tt d} = 1, ~ \gamma =0.01$]FigPlot/Inputs_5_deg_i
[$K_{\tt d} = 1, ~ \gamma =0.1$]FigPlot/Inputs_5_deg_j
[$K_{\tt d} = 1, ~ \gamma =1$]FigPlot/Inputs_5_deg_k
[$K_{\tt d} = 1, ~ \gamma =5$]FigPlot/Inputs_5_deg_l

Position control performance for the desired configuration $\theta_\star = 5 ~{\rm deg}$ with different $\gamma$.

We conducted additional experiments to test the proposed control approach in different scenarios, including the desired configurations of $\theta_\star = 10$ deg and $15$ deg, which are shown in Fig. <ref>. These results demonstrate that the algorithm is capable of achieving high accuracy and performance for position control. To quantify the steady performance, we study the configuration trajectories during $\cali_{\tt s}:=[4,8]$ s, since for all these scenarios the system states arrive at the steady-state stage. For these two desired equilibria, the proposed design achieved high accuracy, verifying the property P1 in Proposition <ref>.

We summarise the accuracy achieved in these experiments with different equilibria ($5$ deg, $10$ deg and $15$ deg) and gains of $\gamma$ and $K_{\tt d}$ in Table <ref>, where $[\theta_{\tt min}, \theta_{\tt max}]$ represents the minimal and the maximal values during the interval $\cali_{\tt s}$, and $\bar \theta$ and $\pm\sigma_\theta$ are the average value and the $\pm 1$ standard deviation. For $\theta_\star = 5$ deg, it achieved the highest accuracy among the three equilibria, for which the selection of $\gamma$ as $1$ and $5$ degraded the steady-state accuracy a little bit. In Fig. <ref>, we present a photo sequence of one of the scenarios with the desired configuration $\theta_\star = 10$ deg, and the gains $K_{\tt d} =0.1$ and $\gamma=1$. This sequence serves as an intuitive illustration of the dynamic behaviour of the closed loop.

In addition, we report the result with a large $\gamma=10$. However, as explained in the discussion point 3) in Section <ref>, a large $\gamma>0$ may make the desired potential energy function $U_{\tt d}$ non-convex, resulting in instability. This is consistent with the experimental results, as we observe the neutral stability with oscillating behaviours at the steady-state stage for $\gamma=10$; see Fig. <ref>.

[$\theta_\star=10$ deg]FigPlot/Inputs_10_deg
[$\theta_\star=15$ deg]FigPlot/Inputs_15_deg
Position control performance for the desired configurations $\theta_\star = 10 ~{\rm deg}$ and $15$ deg with $K_{\tt d} = 1$ and $\gamma =0.1$

The steady-state errors in the time interval $[4,8]$ s of different scenarios (Unit: deg)

2l                     2|c|$K_{\tt d}=10$                                                 2c|$K_{\tt d}=1$                                                                                2c$K_{\tt d}=0.1$
2l|        $\gamma$                   [$\theta_{\tt min}, \theta_{\tt max}$]     $\bar \theta \pm \sigma_\theta$     [$\theta_{\tt min}, \theta_{\tt max}$]     $\bar \theta \pm \sigma_\theta$     [$\theta_{\tt min}, \theta_{\tt max}$]     $\bar \theta \pm \sigma_\theta$

3*$5$ deg     $0.01$      [4.8876, 4.9403]     4.9336 $\pm$ 0.0097
[4.9492, 4.9532]      4.9506 $\pm$ 0.0009     [4.9236, 4.9515]     4.9413 $\pm$ 0.0093

$0.1$       [4.9538, 4.9687]      4.9605 $\pm$ 0.0037
[4.9387, 4.9556]      4.9505 $\pm$ 0.0043
[4.9480, 4.9528]     4.9503 $\pm$ 0.0012

$1$      [4.9148, 4.9265]     4.9219 $\pm$ 0.0024
[4.9148, 4.9265]    4.9219 $\pm$ 0.0024      [4.9156, 4.9195]     4.9179 $\pm$ 0.0008

$5$      [4.9325, 4.9367]     4.9351 $\pm$ 0.0010
[4.9333, 4.9425]    4.9383 $\pm$ 0.0034     [4.9316, 4.9473]     4.9408 $\pm$ 0.0041

3*$ 10$ deg     $0.01$      [9.5437, 10.2709]     9.8678 $\pm$ 0.1055
[9.8328, 9.9597]     9.8824 $\pm$ 0.0276      [9.8458, 9.9596]     9.8902 $\pm$ 0.0021

$0.1$       [9.8320, 9.9048]      9.8803 $\pm$ 0.0200
[9.8328, 9.9597]    9.8824 $\pm$ 0.0276
[9.7990, 9.9007]    9.8827 $\pm$ 0.0232

$ 1$      [9.9074, 9.9295]     9.9216 $\pm$ 0.0060      [9.7814, 9.9044]    9.8454 $\pm$ 0.0309
[9.8204, 9.8457]    9.9324 $\pm$ 0.0081

$5$     [9.8706, 9.8857]     9.8756 $\pm$ 0.0040
[9.8955, 9.8701]     9.8884 $\pm$ 0.0071
[9.8942, 9.8442]     9.8812 $\pm$ 0.0153

3*$ 15$ deg     $0.01$     [12.0252, 14.9968]     14.5734 $\pm$ 0.6592
[14.7067, 14.7753]     14.7650 $\pm$ 0.0139     [14.7333, 14.7798]        14.7711 $\pm$ 0.0075

$0.1$      [14.7071, 14.7555]     14.7375 $\pm$ 0.0147     [14.5142, 14.7350]     14.6810 $\pm$ 0.0711
[14.7174, 17.7939]     14.7710 $\pm$ 0.0228

$ 1$      [14.6276, 14.9968]     14.8058 $\pm$ 0.0673    [14.7499, 14.7794]    14.7632 $\pm$ 0.0089
[14.7228, 14.7597]     14.7481 $\pm$ 0.0116

$5$      [14.6863, 14.7590]    14.7208 $\pm$ 0.0279     [14.7061, 14.7307]     14.7236 $\pm$ 0.0048
[14.7051, 14.7301]    14.7186 $\pm$ 0.0069

A large $\gamma =10$ leads to unstable performance

Finally, we present experimental results demonstrating closed-loop stiffness regulation around the desired equilibria, which is related to P2 in Proposition <ref>. While measuring the overall stiffness is generally not manageable, we can test the transverse stiffness as outlined in Item 1) of Section <ref>. To this end, we equipped a linear actuator perpendicularly to the tangential direction of the continuum robot at the end-effector, as shown in Fig. <ref>. We repeated the experiments for two different desired equilibria, namely 8 deg and 10 deg. We collected stiffness data using different gains $\gamma$ and plotted the results in Figs. <ref> and <ref>. The results match the equation (<ref>) in Section <ref> that the closed-loop stiffness is affine in the control gain $\gamma$. This implied that we were able to identify the parameters $\kappa_1$ and $\kappa_2$, and use them to tune the controller for a prescribed stiffness around the desired configuration.

It is important to note that although various control approaches have been proposed for continuum robotics, their suitability for achieving simultaneous control of position and stiffness in underactuated robots is limited, particularly considering the variations in actuation mechanisms across different continuum robotic platforms. Given the absence of applicable control strategies in the existing literature, this study did not provide experimental comparisons to previous works. However, our aim is to lay the groundwork for future exploration and development of experimental studies in this area.

Stiffness regulation with the desired configuration $\theta_\star =$ 8 deg ($\kappa_1 = 0.0288$ and $\kappa_2 = 0.2332$)

Stiffness regulation with the desired configuration $\theta_\star =$ 10 deg ($\kappa_1 = 0.0805$ and $\kappa_2 = 0.4305$)

§ CONCLUSION

In this paper, we studied the modelling and control of a class of underactuated antagonistic tendon-driven continuum robots. The proposed model possesses a configuration-dependent input matrix, which effectively captures the mechanism for open-loop stiffening through cable tension regulation. We have thoroughly analysed the assignable equilibria set and devised a potential shaping feedback controller that enables simultaneous position-and-stiffness regulation while adhering the non-negative input constraint. To the best of the authors' knowledge, this is the first design for such a problem. The experimental results on the robotic platform OctRobot-I demonstrate the effectivenss and reliability of the proposed approach. Since the proposed approach relies on only a few intrinsic parameters of the model, rather than fully utilising the exact dynamical model, lending it remarking robustness to modelling errors.

Along the research line, the following problems are considered as the potential future works:

1) As per Proposition <ref>, we impose the condition $\gamma<\alpha_2 $ to guarantee the convexity of the desired potential energy function $U_{\tt d}$. It should be noted that the parameter $\alpha_2$ is a intrinsic characteristic of the continuum robot, and consequently the range of choices for the adaptation gain $\gamma$ is limited, which restricts our ability to control the stiffness in a relatively narrow interval. Our experimental results support this assertion. To enlarge the closed-loop stiffness range, a potential way is to make full use of jamming in the continuum robot via changing the compression level of jamming flaps [14].

2) Exploring alternative desired potential energy functions may offer a promising approach to enhance closed-loop performance. In addition, applying state-of-the-art energy shaping methodologies, such as those demonstrated in e.g. [38, 28], could prove valuable for solving more complex tasks for continuum robots, e.g., path following and robust simultaneous position-and-stiffness control.

3) Similar to the recent works [15, 9], the approach in the paper is developed for the planar case. It is underway to extend the main results to multi sections in a three-dimensional space.

4) Our proposed approach does not consider the actuation dynamics, opting instead to utilise a high-gain design to enforce time-scale separation and disregard these dynamics. It would be advantageous to take the actuation dynamics in to the controller synthesis by incorporating advanced robustification techniques [26].

§ A. MODELLING OF POTENTIAL ENERGY FUNCTIONS

In this section, we provide additional details on the model, in particular the potential energy functions of the continuum robotic platform.

§.§.§ Gravitational energy
In order to approximate the gravitational potential energy $U_{\tt G}$, we make the hypothesis that the mass is lumped at the centre of each link with the link lengths $l_i>0$ and the masses $m_{i}>0$ ($i\in \caln$). From some basic geometric relations, we have

\begin{equation}
U_{\tt G}(\bfq) = \sum_{i \in \caln} {l_i m_i \over 2} \left[ \cos(q_0 +\ldots+  q_{i-1})  - \cos(q_0+ \ldots +q_i) \right]
\end{equation}

with the parameter $q_0 =0$, which satisfies $U_{\tt G}(\mathbf{0}_n) = 0$.

Besides, we make the following assumption on the mass and length.

The continuum robot satisfies the uniformity assumptions:

(4a) The masses verify $  m_i = m_j, ~ \forall i,j \in \caln.$

(4b) The lengths satisfy the relation: $l_0 =\ell$ and $l_i = 2\ell$ ($i \in \caln$). The radius of the beam is $r$.

With the above assumption, the potential energy becomes

\begin{equation}
\begin{aligned}
U_{\tt G}(\bfq) &~ = ~\alpha_1 (1-\cos(q_\Sigma))
\\
q_\Sigma &~ :=~ \sum_{i\in \caln} \bfq_i,
\end{aligned}
\end{equation}

with some coefficient $\alpha_1>0$.

§.§.§ Elastic energy
In the designed continuum robot, each spine segment contains two pair of helical compression springs. Since we limit ourselves to the 2-dimensional case, we only consider a pair of springs as illustrated in Fig. <ref>, and make the assumption below.

The deformable part of the continuum manipulator is composed of fixed number of segments with constant curvature with differentiable curves everywhere [9].

Illustration of the elastic energy in springs

In terms of the above assumptions, the boundary lengths in the $i$-th segment are given by

h_{i,1} = q_i\Big[ \ell \cot\Big({q_i \over 2}\Big) + r \Big],
\quad
h_{i,2} = q_i\Big[ \ell \cot\Big({q_i \over 2}\Big) - r \Big].

Note that the above functions are well-posed when $q_i \to 0$, i.e.,

\lim_{q_i \to 0} h_{i,1} =2\ell, \quad \lim_{q_i \to 0} h_{i,2} =2\ell.

Hence, the elastic energy can be modelled as

\begin{equation}
\label{UE:real}
U_{\tt E} = \sum_{i \in \caln} {k_i} \left[ q_i^2\Big(\ell^2 \cos\Big({q_i \over 2}\Big)^2 + r^2\Big) - \ell^2 \right] + k_i' q_i^2,
\end{equation}

in which $k_i>0$ and $k_i'>0$ are some elastic coefficients to characterise the elastic energies caused by elongation and bending of springs.

In the high-dimensional rigid-link model, each configuration variable $q_i$ would generally be small, i.e., $q_i\in [- {\pi\over 12}, {\pi\over 12}]$, for which the term $\cos({q_i \over 2})^2$ takes values within $[0.983,1]$. Then, it is reasonable to make the following quadratic assumption to approximate the highly nonlinear function in (<ref>).

The elastic energy function $U_{\tt E}$ has the form

\begin{equation}
U_{\tt E}(\bfq) = \bfq^\top \Lambda \bfq + U_0
\end{equation}

with a constant coefficient $U_0$ and a diagonalisable matrix $\Lambda:=\diag(\alpha_2, \ldots \alpha_2) \succ 0$.

[1]
C. Armanini, F. Boyer, A. T. Mathew, C. Duriez, and F. Renda, “Soft robots
modeling: A structured overview,” IEEE Trans. Robot., pp. 1–21,
2023, (early access).

[2]
A. Bajo and N. Simaan, “Hybrid motion/force control of multi-backbone
continuum robots,” Int. J. Robot. Res., vol. 35, no. 4, pp. 422–434,

[3]
C. M. Best, L. Rupert, and M. D. Killpack, “Comparing model-based control
methods for simultaneous stiffness and position control of inflatable soft
robots,” Int. J. Robot. Res., vol. 40, no. 1, pp. 470–493, 2021.

[4]
T. M. Bieze, F. Largilliere, A. Kruszewski, Z. Zhang, R. Merzouki, and
C. Duriez, “Finite element method-based kinematics and closed-loop control
of soft, continuum manipulators,” Soft Robot., vol. 5, no. 3, pp.
348–364, 2018.

[5]
D. Braganza, D. M. Dawson, I. D. Walker, and N. Nath, “A neural network
controller for continuum robots,” IEEE Trans. Robot., vol. 23, no. 6,
pp. 1270–1277, 2007.

[6]
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.

[7]
B. Caasenbrood, A. Pogromsky, and H. Nijmeijer, “Energy-shaping controllers
for soft robot manipulators through port-Hamiltonian Cosserat models,”
SN Comput. Sci., vol. 3, no. 6, p. 494, 2022.

[8]
C. Della Santina, C. Duriez, and D. Rus, “Model-based control of soft robots:
A survey of the state of the art and open challenges,” IEEE Control
Syst. Mag., vol. 43, no. 3, pp. 30–65, 2023.

[9]
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.

[10]
B. Deutschmann, A. Dietrich, and C. Ott, “Position control of an underactuated
continuum mechanism using a reduced nonlinear model,” in Proc. IEEE
Conf. Decis. Control.1em plus 0.5em minus 0.4emIEEE, 2017, pp.

[11]
Y. Engel, P. Szabo, and D. Volkinshtein, “Learning to control an octopus arm
with Gaussian process temporal difference methods,” in Adv. Neural.
Inf. Process. Syst., vol. 18, 2005.

[12]
V. Falkenhahn, A. Hildebrandt, R. Neumann, and O. Sawodny, “Dynamic control of
the bionic handling assistant,” IEEE/ASME Trans. Mechatron., vol. 22,
no. 1, pp. 6–17, 2016.

[13]
F. Fan, B. Yi, D. Rye, G. Shi, and I. R. Manchester, “Learning stable
Koopman embeddings,” in Proc. Am. Control Conf., 2022, pp.

[14]
Y. Fan, D. Liu, and L. Ye, “A novel continuum robot with stiffness variation
capability using layer jamming: Design, modeling, and validation,”
IEEE Access, vol. 10, pp. 130 253–130 263, 2022.

[15]
E. Franco and A. Garriga-Casanovas, “Energy-shaping control of soft continuum
manipulators with in-plane disturbances,” Int. J. Robot. Res.,
vol. 40, no. 1, pp. 236–255, 2021.

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

[17]
A. D. Kapadia, K. E. Fry, and I. D. Walker, “Empirical investigation of
closed-loop control of extensible continuum manipulators,” in Proc.
IEEE/RSJ Int. Conf. Intell. Robots Syst.1em plus 0.5em minus
0.4emIEEE, 2014, pp. 329–335.

[18]
M. Keppler, D. Lakatos, C. Ott, and A. Albu-Schäffer, “Elastic structure
preserving (ESP) control for compliantly actuated robots,” IEEE
Trans. Robot., vol. 34, no. 2, pp. 317–335, 2018.

[19]
H. K. Khalil, Nonlinear Systems, 3rd ed.1em plus 0.5em minus
0.4emPatience Hall, 2002.

[20]
M. J. Kim, A. Werner, F. Loeffl, and C. Ott, “Passive impedance control of
robots with viscoelastic joints via inner-loop torque control,” IEEE
Trans. Robot., vol. 38, no. 1, pp. 584–598, 2021.

[21]
M. Mahvash and P. E. Dupont, “Stiffness control of surgical continuum
manipulators,” IEEE Trans. Robot., vol. 27, no. 2, pp. 334–345,

[22]
R. Mengacci, F. Angelini, M. G. Catalano, G. Grioli, A. Bicchi, and
M. Garabini, “On the motion/stiffness decoupling property of articulated
soft robots with application to model-free torque iterative learning
control,” Int. J. Robot. Res., vol. 40, no. 1, pp. 348–374, 2021.

[23]
R. Ortega, M. W. Spong, F. Gómez-Estern, and G. Blankenstein,
“Stabilization of a class of underactuated mechanical systems via
interconnection and damping assignment,” IEEE Trans. Autom. Control,
vol. 47, no. 8, pp. 1218–1233, 2002.

[24]
R. Ortega, A. van der Schaft, B. Maschke, and G. Escobar, “Interconnection and
damping assignment passivity-based control of port-controlled hamiltonian
systems,” Automatica, vol. 38, no. 4, pp. 585–596, 2002.

[25]
R. Ortega, A. J. Van Der Schaft, I. Mareels, and B. Maschke, “Putting energy
back in control,” IEEE Control Syst. Mag., vol. 21, no. 2, pp.
18–33, 2001.

[26]
R. Ortega, B. Yi, and J. G. Romero, “Robustification of nonlinear control
systems vis-à-vis actuator dynamics: An immersion and invariance
approach,” Syst. Control Lett., vol. 146, 2020, Art. no. 104811.

[27]
F. Renda, M. Giorelli, M. Calisti, M. Cianchetti, and C. Laschi, “Dynamic
model of a multibending soft robot arm driven by cables,” IEEE Trans.
Robot., vol. 30, no. 5, pp. 1109–1122, 2014.

[28]
J. G. Romero, A. Donaire, and R. Ortega, “Robust energy shaping control of
mechanical systems,” Syst. Control Lett., vol. 62, no. 9, pp.
770–780, 2013.

[29]
J. K. Salisbury, “Active stiffness control of a manipulator in Cartesian
coordinates,” in Proc. IEEE Conf. Decis. Control.1em plus
0.5em minus 0.4emIEEE, 1980, pp. 95–100.

[30]
M. Takegaki and S. Arimoto, “A new feedback method for dynamic control of
manipulators,” ASME J. Dyn. Syst. Meas. Control, vol. 102, pp.
119–125, 1981.

[31]
T. G. 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.

[32]
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,

[33]
H. Tsukamoto, S.-J. Chung, and J.-J. E. Slotine, “Contraction theory for
nonlinear stability analysis and learning-based control: A tutorial
overview,” Annu. Rev. Control, vol. 52, pp. 135–169, 2021.

[34]
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., 2023, (early access).

[35]
A. van der Schaft, $L_2$-Gain and Passivity Techniques in Nonlinear
Control.1em plus 0.5em minus 0.4emSpringer, 2000.

[36]
J. Wang and A. Chortos, “Control strategies for soft robot systems,”
Adv. Intell. Syst, vol. 4, no. 5, 2022, Art. no. 2100165.

[37]
B. Yi and I. R. Manchester, “On the equivalence of contraction and Koopman
approaches for nonlinear stability and control,” IEEE Trans. Autom.
Control, pp. 1–16, 2023, (conditionally accepted).

[38]
B. Yi, R. Ortega, D. Wu, and W. Zhang, “Orbital stabilization of nonlinear
systems via mexican sombrero energy shaping and pumping-and-damping
injection,” Automatica, vol. 112, 2020, Art. no. 108661.

[39]
J. Zhang, Y. Li, Z. Kan, Q. Yuan, H. Rajabi, Z. Wu, H. Peng, and W. Jianing,
“A preprogrammable continuum robot inspired by elephant trunk for dexterous
manipulation,” Soft Robot., no. 1, 2023, (early access).

[]Bowen Yi obtained his Ph.D. degree in Control Engineering from Shanghai Jiao Tong University, China in 2019. From 2017 to 2019 he was a visiting student at Laboratoire des Signaux et Systèmes, CNRS-CentraleSupélec, Gif-sur-Yvette, France.

He has hold postdoctoral positions in Australian Centre for Field Robotics, The University of Sydney, NSW, Australia (2019 - 2022), and the Robotics Institute, University of Technology Sydney, NSW, Australia (Sept. 2022 - now). His research interests involve nonlinear systems and robotics. He received the 2019 CCTA Best Student Paper Award from the IEEE Control Systems Society for his contribution on sensorless observer design.

[]Yeman Fan (Student Member, IEEE) received the B.E. degree in Machine Designing, Manufacturing and Automation, and the M.E. degree in Agricultural Electrification and Automation from Northwest A&F University, Yangling, China, in 2016 and 2019, respectively. He is currently pursuing the Ph.D. degree with the Robotics Institute, University of Technology Sydney, Sydney, NSW, Australia.

His research interests include continuum robots and manipulators, robot control systems, and jamming technology for robotics.

[]Dikai Liu (Senior Member, IEEE) received the Ph.D. degree in Dynamics and Control from the Wuhan University of Technology, Wuhan, China, in 1997.

He is currently a Professor in Mechanical and Mechatronic Engineering with the Robotics Institute, University of Technology Sydney, Sydney, NSW, Australia. His main research interests include robotics, including robot perception, planning and control of mobile manipulators operating in complex environments, human-robot collaboration, multi-robot coordination, and bioinspired robotics.

[]José Guadalupe Romero (Member, IEEE) obtained the Ph.D. degree in Control Theory from the University of Paris-Sud XI, France in 2013.

Currently, he is a full time Professor at the Instituto Tecnológico Autónomo de México (ITAM), Mexico City, Mexico, and since 2019 he has been the Director of undergraduate mechatronics engineering program. His research interests are focused on nonlinear and adaptive control, stability analysis and the state estimation problem, with application to mechanical systems, aerial vehicles, mobile robots and multi-agent systems.