Energy-shaping control of soft continuum manipulators with in-plane disturbances

Abstract

Soft continuum manipulators offer levels of compliance and inherent safety that can render them a superior alternative to conventional rigid robots for a variety of tasks, such as medical interventions or human–robot interaction. However, the ability of soft continuum manipulators to compensate for external disturbances needs to be further enhanced to meet the stringent requirements of many practical applications. In this paper, we investigate the control problem for soft continuum manipulators that consist of one inextensible segment of constant section, which bends under the effect of the internal pressure and is subject to unknown disturbances acting in the plane of bending. A rigid-link model of the manipulator with a single input pressure is employed for control purposes and an energy-shaping approach is proposed to derive the control law. A method for the adaptive estimation of disturbances is detailed and a disturbance compensation strategy is proposed. Finally, the effectiveness of the controller is demonstrated with simulations and with experiments on an inextensible soft continuum manipulator that employs pneumatic actuation.

1. Introduction

Soft continuum manipulators are devices made of low-stiffness materials, which are typically exploited to achieve large structural deformations and displacements (Rus and Tolley, 2015). Such devices offer many attractive features, including compliance, light weight, inherent safety, and relatively simple miniaturization. These features make them ideally suited for human–robot interaction and for operation in cluttered environments, such as those found in medicine (DeGreef et al., 2009; Gerboni et al., 2017). The development of soft continuum manipulators has received significant attention, particularly in recent years (Burgner-Kahrs et al., 2015; Gerboni et al., 2015; Trivedi et al., 2008; Webster and Jones, 2010; Wehner et al., 2016). Among the different actuation strategies proposed, pressurized fluids, and particularly pneumatics, is one of the most popular. This can be credited to the high power-to-weight ratio, fast response, affordability, compatibility with magnetic resonance imaging (MRI), and miniaturization potential that soft continuum manipulators with pneumatic actuation can offer. A pioneer device in this category is the flexible micro actuator (FMA) (Suzumori, 1996; Suzumori et al., 1991, 1992). Since its invention, multiple other FMA designs have been proposed and developed (e.g. Abe et al., 2007; Chen et al., 2009; Cianchetti et al., 2014; Garriga-Casanovas et al., 2018; Marchese and Rus, 2016; McMahan et al., 2006; Mosadegh et al., 2014; Wehner et al., 2016), as well as torsional soft actuators (Sanan et al., 2014). Despite all this progress on design and manufacturing of FMAs and soft continuum manipulators in general, the study of model-based control strategies specific to these systems remains a largely unexplored field (Thuruthel et al., 2018).

Control strategies for soft continuum manipulators have been proposed in the literature using model-free and model-based methods. Numerical model-free methods have the advantage of not relying on a dynamical model of the system, which can be difficult to derive analytically for general cases or can be prohibitively time consuming if computed numerically. A model-free controller based on an adaptive Kalman filter was recently proposed by Li et al. (2018), while machine learning techniques were employed by Rolf and Steil (2014) and Thuruthel et al. (2017) to construct inverse kinematic models for control purposes. A common drawback of machine learning approaches is the need for training data, while the study of stability conditions remains an open problem. In addition, low-level control of the actuators in Thuruthel et al. (2017) employs high gains. As highlighted by DellaSantina et al. (2017), high-gain feedback control imposes de facto a reduction of the compliance of the system, thus potentially defeating the purpose of physical compliance. Recent research on finite-element (FE) models (Bieze et al., 2018; Goury and Duriez, 2018) has shown promising results in reducing computation time. This has allowed using FE models for control purposes (Bieze et al., 2018; Zhang et al., 2016), but this approach is only intended for quasi-static conditions.

Model-based closed-loop control has been indicated as particularly promising for accurate positioning (Thuruthel et al., 2018). In addition, analytical control provides the tools to analyze the stability of the closed-loop system with respect to the dynamical model. Recent model-based controllers for soft continuum manipulators rely increasingly on concentrated-parameters models (DellaSantina et al., 2018; Falkenhahn et al., 2015; Godage et al., 2015; Sadati et al., 2018), which typically introduce the assumption of constant curvature (CC) or piecewise-constant curvature (PCC). The resulting discrete models are computationally more efficient than continuous models based on beam theory (Renda et al., 2014; Rucker and Webster, 2011), and thus more suitable for closed-loop control. However, the assumptions of CC and PCC might not be verified in the presence of disturbances. In addition, concentrated-parameters models only provide an approximation of the complex dynamics of soft continuum manipulators and are consequently less accurate than continuous models based on Cosserat rod theory (Alqumsan et al., 2019; Grazioso et al., 2018; Till et al., 2019). Lastly, a promising emerging trend in the modeling of soft continuum manipulators consists of employing a port-Hamiltonian approach (Ross et al., 2016), which focuses on the energy associated with the system (Moghadam et al., 2016). The main advantages of port-Hamiltonian modeling are the general applicability to different physical domains and the common formalism with energy-based control techniques.

Notable results in model-based control of soft continuum manipulators include feedback linearization techniques (Deutschmann et al., 2017; Gravagne et al., 2003) and optimal control (Falkenhahn et al., 2017). In addition, the combination of feedback and feed-forward actions was proposed by DellaSantina et al. (2017) in order to enhance robustness to uncertainties while preserving the compliance of the manipulator in closed loop, thanks to the use of small gains. More recently, sliding-mode control (SMC) algorithms, such as those reported in Slotine and Li (1991), were implemented for this class of systems by Alqumsan et al. (2019) for the case of bounded disturbances. Nevertheless, model-based dynamic controllers for soft continuum manipulators are still in their nascent stage. One of the main challenges is posed by unstructured environments due to the presence of unknown disturbances (Thuruthel et al., 2018).

Employing a concentrated-parameters model for control purposes enables a simpler controller design and a reduced computation time. Constructing a rigid-link model of a soft continuum manipulator involves associating unactuated degrees-of-freedom (DOFs) to the compliant elements, in a similar manner to flexible mechanisms (Franco et al., 2018; Yu et al., 2005). Controlling the resulting underactuated model in the presence of disturbances is a challenging engineering problem. Energy-shaping methods (Bloch et al., 2000; Ortega et al., 2002) are among the preferred choices for the control of underactuated mechanisms and they also appear particularly suited for soft continuum manipulators since they do not rely on high gains, hence they can preserve the system compliance in closed loop. In addition, robust or adaptive versions of energy-shaping algorithms have recently been developed to compensate for the effects of different types of disturbances (Donaire et al., 2017; Franco, 2019b). In particular, disturbances that only affect the actuated DOF are termed matched, while unmatched disturbances affect the unactuated DOF. One of the main difficulties associated with energy-shaping control is the need to solve a set of partial-differential equations (PDEs), which is notoriously problematic for systems with a large number of DOFs. While a variety of approaches have been investigated to circumvent this obstacle (Donaire et al., 2016; Nunna et al., 2015), they typically rely on restrictive assumptions and are difficult to generalize. In this respect, computing energy-shaping control laws in closed-form for larger classes of systems represents an open research problem with significant practical relevance, and is of immediate interest for the control of soft continuum manipulators (Katzschmann et al., 2019).

In this paper, an energy-shaping approach is investigated for the dynamic control of soft continuum manipulators consisting of one inextensible segment with a constant section, similar to Garriga-Casanovas et al. (2018), which bends due to internal pressurization resulting in a single control input, and is subject to unknown in-plane disturbances, which can be either matched or unmatched. To this end, a rigid-link model and a port-Hamiltonian formulation are employed. The main contributions of this work can be summarized as follows.

1. An adaptive algorithm for the online estimation of a class of linearly parameterized disturbances acting in the plane of bending is presented.

2. A controller design procedure is outlined and control laws are derived in closed-form for the case of an $n$ DOF rigid-link underactuated model, where $n$ can be arbitrarily large.

3. A disturbance compensation strategy is implemented and stability conditions are discussed. Differently from other approaches, prior knowledge of the disturbance bounds is not required in this case.

4. An experimental validation of the proposed controller using a soft continuum manipulator consisting of one inextensible segment with pneumatic actuation in the presence of external forces.

It should be noted that this paper focusses on a planar case similarly to other recent work on the topic (DellaSantina et al., 2017, 2018; Marchese and Rus, 2016). However, it is shown that the results can be readily applied to the control of this type of manipulator in three-dimensional (3D) space.

The rest of the paper is organized as follows. In Section 2, the model of the manipulator is briefly introduced and the assumptions that define the control problem are detailed. In Section 3 the control laws are derived and the stability conditions are discussed. In Section 4, the simulation results for a rigid-link model and the experimental results for a soft continuum manipulator are presented. Finally, concluding remarks are summarized in Section 5, together with suggestions for future work.

2. Problem formulation

2.1. Model overview

A schematic of the soft continuum manipulator that serves as reference design for this work is shown in Figure 1. The manipulator consists of a tubular structure made of a hyper-elastic material (Elastosil M 4601) with a constant cross-section defining three equal internal chambers spaced at $2\pi /3$ . An inextensible nylon fiber is embedded at the center of the cross-section along the entire device to prevent longitudinal extension and to increase the capability of supporting external wrenches (Garriga-Casanovas et al., 2018). A second inextensible fiber is wound around the outer wall to prevent radial expansion while allowing longitudinal strain. This manipulator can thus bend in any direction in 3D space by applying a specific set of pressures $P_0,P_1,P_2$ in the internal chambers, while variation in length and external diameter are negligible. In the absence of disturbances, the orientation of the bending plane $\phi$ and the rotation of the tip $\theta$ on the plane of bending at equilibrium can be approximated as in (Suzumori, 1996)

$$\theta =\frac{1}{k}\sqrt{P_0^2+P_1^2+P_2^2-P_1{P}_2-P_1{P}_0-P_2{P}_0}$$

(1.a)

$$\tan \left(\phi \right)=\frac{\sqrt{3}(P_1-P_2)}{P_1+P_2-2P_0}$$

(1.b)

where $k$ is the stiffness of the structure.

For the purpose of deriving closed-loop control laws, the manipulator is modeled as a rigid-link system, similarly to Godage et al. (2015). The link lengths are defined as $l_i$ . The joint angles are defined as $q_i$ , which represent the rotation of link $i$ with respect to link $i-1$ . A point mass $m_i$ is located at the mid-point of each link, and their sum equals the total mass of the manipulator $m_T$ . The bending moment generated by the combined effect of the internal pressures has amplitude $u=\sqrt{P_0^2+P_1^2+P_2^2-P_1{P}_2-P_1{P}_0-P_2{P}_0}$ and corresponds to the control input of the dynamical model. Setting the ratio between $P_0,P_1,P_2$ allows selecting the plane of bending, as shown in Figure 1(b) and expressed in Equation (1.b).

Fig. 1.

(a) Spatial configuration of the soft continuum manipulator and corresponding rigid-link model with $n=3$ . (b) Section view of the manipulator showing the internal chambers.

The dynamics of the $n$ DOF rigid-link model in Figure 1 is expressed employing a port-Hamiltonian formulation. To this end, we define the position $q\in {\mathbb{R}}^n$ and the momenta $p=M\stackrel{·}{q}\in {\mathbb{R}}^n$ , where $M(q)$ is the positive-definite and invertible inertia matrix. The open-loop Hamiltonian $H=T(q,p)+V(q)$ corresponds to the total energy of the system, with $T(q,p)=\frac{1}{2}{p}^{\mathrm{T}}{M}^{-1}p$ representing the kinetic energy, and $V(q)$ representing the potential energy. The control input is $u>0$ and the input mapping is $G$ , with $\mathrm{rank}\left(G\right)<n$ indicating underactuation. The damping matrix is $D>0$ , while the disturbances $\delta \in {\mathbb{R}}^n$ represent the effect of model uncertainties and of generalized external forces acting in the plane of bending. The open-loop system dynamics is defined as

$$\left[\begin{array}{c}\stackrel{·}{q}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{p}}\end{array}\right]=\left[\begin{array}{cc}0\hfill & I\hfill \\ \multicolumn{1}{c}{-I}& -D\hfill \end{array}\right]\left[\begin{array}{c}{\nabla }_{q}H\hfill \\ \multicolumn{1}{c}{{\nabla }_{p}H}\end{array}\right]+\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{G}\end{array}\right]u-\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{\delta }\end{array}\right]$$

(2)

The $n\times n$ identity matrix is indicated with $I$ , the symbol ${\nabla }_q(·)$ represents the vector of partial derivatives in $q$ , and the symbol ${\nabla }_p(·)$ represents the vector of partial derivatives in $p$ .

2.2. Assumptions

The following assumptions, which define the control problem considered in this work, are introduced and briefly discussed.

Assumption 1: The states $(q,p)\in {\mathbb{R}}^{2n}$ are exactly known.

Assumption 2: The potential energy only contains elastic terms, is quadratic in the position, and is expressed as $V\left(q\right)=\frac{k}{2}\sum _{i=1}^n{q}_i^2$ , where the stiffness of the structure $k$ is constant and identical for all joints. Similarly, damping is constant and identical for all joints, thus $D=D_{0}I>0$ . The effect of gravity, which is not included in $V\left(q\right)$ , is accounted for instead within the disturbance $\delta$ .

Assumption 3: The input matrix $G\in {\mathbb{R}}^{n\times 1}$ is constant and has identical elements equal to 1.

Assumption 4: In the presence of disturbances there exists a non-empty set ${\Omega }^{*}$ of attainable open-loop equilibrium positions $q^{*}$ that satisfy the following condition, where $G^{\perp }$ is a full-rank left annihilator of $G$ such that $G^{\perp }G=0$ and such that $\mathrm{rank}\left(G^{\perp }\right)=n-1$

$$G^{\perp }({\nabla }_{q}V(q^{*})+\delta )=0$$

(3)

Assumption 5: The disturbances only act in the plane of bending and are linearly parameterized as $\delta ={\delta }_{0}f$ , where $f(q)$ is a known scalar function of the state and ${\delta }_0$ is a vector of unknown parameters, potentially time varying, with time derivative ${\stackrel{·}{\delta }}_0$ .

In particular, Assumption 1 implies that $(q,p)$ are known from measurements or from an appropriate observer, such as Venkatraman et al. (2010). Including the observer in the control formulation is beyond the scope of this paper and is part of our future work. According to Assumption 2 the effect of gravity is not accounted for in the potential energy. Instead, the effect of the weight of the manipulator and of any payload, and any un-modeled dynamics (e.g. variations in structural stiffness due to pressurization) are included in the disturbances $\delta$ for controller design purposes. Assumption 3 implies that the bending moment generated by the control input $u$ is uniformly distributed over all joints. The expression of the matrix $G$ specifies that only the dynamics in the plane of bending is considered. This is motivated by the fact that, for the type of soft continuum manipulators considered here, the orientation of the bending plane is defined by the geometry of the internal chambers and by the corresponding pressures according to (1.b). Assumption 4 implies that under a given external load $\delta$ the manipulator acquires the configuration $q^{*}$ and is typically verified if the load is commensurate with the material properties so that neither rupture nor plastic deformation occurs, while hysteretic effects are negligible. Assumption 5 limits the scope of this work to generalized forces that act in the plane of bending. Relaxing this assumption to include out-of-plane disturbance is part of our future work. For illustrative purposes, a force ${\delta }_0$ acting along the $y$ -axis at the tip of the distal link produces a bending moment $\delta ={\delta }_{0}f$ , where the component $i$ of $f$ is

$$f_i=\sum _{k=i}^n{l}_k\sin \left(\sum _{j=1}^k{q}_j\right)$$

(4)

The effect of un-modeled dynamics is accounted for by the time-varying parameter ${\delta }_0$ .

Note finally that the rigid-link model approximates the dynamics of a soft continuum manipulator, since it is underactuated and it can include an arbitrary number $n$ of DOFs, but it is not specific to pneumatic actuation. Thus, a similar approach could be employed for different actuation strategies, such as tendons or electro-active polymer actuators (Mattioni et al., 2018).

Remark 1: If $G^{\perp }\delta =0$ the disturbances are only matched and solving (3) gives $q_1^{*}=q_i^{*},\forall i\le n$ , confirming that CC is an attainable equilibrium for system (2). In this case, the model in quasi-static conditions reduces to 1 DOF. Differently from classical examples of underactuated mechanisms (Franco et al., 2018), the class of systems considered here admits a set ${\Omega }^{*}$ of attainable equilibria instead of a single position $q^{*}$ . Thus, by employing an appropriate control action, it is theoretically possible to stabilize any CC configuration if the disturbances are only matched. The latter is the case for a moment load at the tip of the manipulator. Conversely, in the presence of unmatched disturbances, the attainable equilibrium is not CC and a 1 DOF model is no longer representative even in quasi-static conditions. This scenario includes in-plane forces applied at the tip of the manipulator, the weight of the manipulator or of any payload, and uncertain non-uniform stiffness.

3. Controller design

3.1 Adaptive estimate of the disturbances

The unknown disturbances are estimated from the open-loop dynamics (2) using the Immersion and Invariance (I&I) method (Astolfi and Ortega, 2003; Astolfi et al., 2007), as outlined in the following proposition.

Proposition 1: Consider system (2) under Assumptions 1,5 with ${\stackrel{·}{\delta }}_0=0$ . Then, the disturbance estimate $\tilde{\delta }={\tilde{\delta }}_{0}f$ converges to the correct value exponentially with the following adaptation law and with the tuning parameter $\alpha >0$

$${\tilde{\delta }}_0=-\alpha p^{\mathrm{T}}f+\alpha \int f^{\mathrm{T}}(-{\nabla }_{q}H-D{\nabla }_{p}H+\mathrm{Gu}-{\tilde{\delta }}_{0}f)\mathrm{dt}$$

(5)

Proof: We define the vector of estimation errors $z$ as follows, where the functions ${\hat{\delta }}_0$ and $\beta \left(p\right)$ are the state-independent part and the state-dependent part of the estimate ${\tilde{\delta }}_0$

$$z={\tilde{\delta }}_0-{\delta }_0={\hat{\delta }}_0+\beta (p)-{\delta }_0$$

(6)

Computing the time derivative of (6) and substituting (2) and (6) gives

$$\stackrel{·}{z}={\stackrel{·}{\hat{\delta }}}_0+{\nabla }_p{\beta }^{\mathrm{T}}(-{\nabla }_{q}H-D{\nabla }_{p}H+\mathrm{Gu}-({\hat{\delta }}_0+\beta -z)f)$$

(7)

Substituting ${\stackrel{·}{\hat{\delta }}}_0=\alpha f^{\mathrm{T}}(-{\nabla }_{q}H-D{\nabla }_{p}H+\mathrm{Gu}-{\tilde{\delta }}_{0}f)$ and $\beta =-\alpha p^{\mathrm{T}}f$ into (7), which correspond to the adaptation law (5), gives

$$\stackrel{·}{z}=-\alpha z{f}^{\mathrm{T}}f$$

(8)

Finally, choosing the Lyapunov function candidate $W=\frac{1}{2}{z}^{\mathrm{T}}z$ , computing its time derivative, and substituting (8) gives $\stackrel{·}{W}=-\alpha {\left|\mathrm{zf}\right|}^2<0$ . Consequently, the quantity $\left|\mathrm{zf}\right|$ is bounded and converges to zero exponentially, concluding the proof ▪

Remark 2: In the case in which $f$ is constant (e.g. constant bending moment) the disturbances reduce to $\delta ={\delta }_0$ and the adaptation law (5) can be simplified as

$$\tilde{\delta }=\beta \left(p\right)+\hat{\delta }=-\alpha p+\alpha \int (-{\nabla }_{q}H-D{\nabla }_{p}H+\mathrm{Gu}-\tilde{\delta })\mathrm{dt}$$

(9)

Employing (9) in the presence of variable disturbances does not ensure convergence of the estimate to the correct value. Nevertheless, the effect of the disturbances can be compensated for by an appropriate control action (see Corollary 1). Both the adaptive estimates (5) and (9) require full state feedback, which motivates the need for Assumption 1 in this work. This point becomes apparent observing that the term ${\nabla }_{q}H$ contains the partial derivatives in $q$ of the inertia matrix $M$ , which depend on the individual joint angles (the inertia matrix for $n=3$ is reported in Appendix A). In practice, Assumption 1 can be obviated with a simplified adaptation law (see Corollary 1). Conversely, the remaining Assumptions 2–4 are not required for the purpose of disturbance estimation.

3.2 Control law for CC equilibrium

The control aim considered in this work corresponds to stabilizing the desired equilibrium $(q,p)=(q_d,0)$ with $q_d\in {\Omega }^{*}$ ; hence, it is a setpoint regulation problem. Introducing the matrix $G^{†}={\left(G^{\mathrm{T}}G\right)}^{-1}{G}^{\mathrm{T}}$ , the energy-shaping control law is defined as the sum of the term $u_{\mathrm{es}}$ , which assigns the closed-loop equilibrium $q_d$ , of the term $u_{\mathrm{di}}$ , which injects damping in the system necessary to attain stability at equilibrium (Ortega et al., 2002), and of a disturbance compensation term $u^{*}$ :

$$\begin{array}{c}u=u_{\mathrm{es}}+u_{\mathrm{di}}+u^{*}\hfill \\ \begin{array}{c}{u}_{\mathrm{es}}=G^{†}({\nabla }_{q}H-M_d{M}^{-1}{\nabla }_q{H}_d+J_2{\nabla }_p{H}_d)\hfill \\ \begin{array}{c}{u}_{\mathrm{di}}=-k_v{G}^{\mathrm{T}}{\nabla }_p{H}_d\hfill \\ \multicolumn{1}{c}{u^{*}=G^{†}\tilde{\delta }}\end{array}\hfill \end{array}\hfill \end{array}$$

(10)

where the parameter $k_v=k_v^{\mathrm{T}}>0$ is a constant gain governing the damping injection, and $J_2=-J_2^{\mathrm{T}}$ is a free matrix dependent on $p$ . The closed-loop dynamics becomes then

$$\left[\begin{array}{c}\stackrel{·}{q}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{p}}\end{array}\right]=\left[\begin{array}{cc}0\hfill & M^{-1}{M}_d\hfill \\ \multicolumn{1}{c}{-M_d{M}^{-1}}& J_2-D{M}^{-1}{M}_d-G{k}_v{G}^{\mathrm{T}}\hfill \end{array}\right]\left[\begin{array}{c}{\nabla }_q{H}_d\hfill \\ \multicolumn{1}{c}{{\nabla }_p{H}_d}\end{array}\right]$$

(11)

where $H_d=\frac{1}{2}{p}^{\mathrm{T}}{M}_d^{-1}p+V_d$ is the closed-loop Hamiltonian and $q_d=\mathrm{argmin}(V_d)$ corresponds to a strict minimizer of the closed-loop potential energy $V_d$ . Equating (2) and (11) and pre-multiplying by $G^{\perp }$ we obtain the potential-energy PDE and the kinetic-energy PDE, which should be satisfied for all $(q,p)\in {\mathbb{R}}^{2n}$

$$G^{\perp }({\nabla }_{q}V-M_d{M}^{-1}{\nabla }_q{V}_d)=0$$

(12.a)

$$G^{\perp }({\nabla }_q\left(p^{\mathrm{T}}{M}^{-1}p\right)-M_d{M}^{-1}\left(p^{\mathrm{T}}{M}_d^{-1}p\right)+2J_2{M}_d^{-1}p)=0$$

(12.b)

Finding a closed-form solution of (12) is typically problematic for systems with a large number of DOFs. This is achieved here with the kinetic-energy-shaping design $M_d=k_{m}M$ , where $k_m$ is a positive tuning parameter, which solves (12.b) with $J_2=0$ . Expressing (12.a) for system (2) under Assumptions 2 and 3 results then in a system of $n-1$ linear PDEs that admits solutions in the following form

$$\begin{array}{c}{V}_d=\frac{k}{2k_m}\sum _{i=1}^n{q}_i^2-\frac{k}{2k_m}{\left(\sum _{i=1}^n{q}_i\right)}^2+\mathrm{\Phi }\left(\sum _{i=1}^n{q}_i\right)\hfill \\ \multicolumn{1}{c}{\mathrm{\Phi }=\frac{k}{k_m}\frac{n-1}{2n}{\left(\sum _{i=1}^n{q}_i\right)}^2+\frac{k_p}{2}{(\sum _{i=1}^n{q}_i-n{q}_d)}^2}\end{array}$$

(13)

Choosing $\mathrm{\Phi }$ as in (13) satisfies the following strict-minimizer conditions for all $q_d\in {\Omega }^{*}$

$${\nabla }_q{V}_d\left(q_d\right)=0$$

(14.a)

$${\nabla }_q^2{V}_d\left(q_d\right)=n(k/k_m{)}^{n-1}{k}_p>0$$

(14.b)

Substituting (13) into (10) gives $u_{\mathrm{es}}$ for a generic $n$ DOF model

$$u_{\mathrm{es}}=\frac{k}{n}\left(\sum _{i=1}^n{q}_i\right)-k_p{k}_m(\sum _{i=1}^n{q}_i-n{q}_d)$$

(15)

Alternative expressions of $\mathrm{\Phi }$ result in different control laws (see Appendix B).

Computing the term $u_{\mathrm{di}}$ gives

$$u_{\mathrm{di}}=-\frac{k_v}{k_m}\sum _{i=1}^n{\stackrel{·}{q}}_i$$

(16)

Proposition 2: Consider system (2) under Assumptions 1–5 in closed loop with controller (10), where $u_{\mathrm{es}}$ is defined as (15) and $u_{\mathrm{di}}$ is defined as (16). Assume that the disturbances are only matched and are estimated with (5). Then the equilibrium $(q,p)=(q_d,0)$ is locally stable for all $k_v>0$ and $\alpha >(1+k_v)/4k_v$ provided that $D_0{\left|\stackrel{·}{q}\right|}^2>k_m{\left|\stackrel{·}{\delta }\right|}^2$ . In addition, $q_d\in {\Omega }^{*}$ is a strict minimizer of $V_d$ if $k_p>0$ .

Proof: To streamline the proof of the first claim we initially consider the case of null disturbances. Choosing the Lyapunov function candidate $H_d=\frac{1}{2}{p}^{\mathrm{T}}{M}_d^{-1}p+V_d$ , computing its time derivative, substituting (2), (10), and recalling that $M_d=k_{m}M$ gives

$${\stackrel{·}{H}}_d={\nabla }_q{H_d}^{\mathrm{T}}\stackrel{·}{q}+{\nabla }_p{H}_d^{\mathrm{T}}\stackrel{·}{p}=-\frac{k_v}{k_m^2}{\left|G^{\mathrm{T}}\stackrel{·}{q}\right|}^2-\frac{1}{k_m}{\stackrel{·}{q}}^{\mathrm{T}}D\stackrel{·}{q}\le 0$$

(17)

Since $D>0$ , asymptotic stability of the equilibrium is concluded for all $k_v\ge 0$ .

Matched disturbances can be expressed as $G\delta$ , which substituted in (2) gives $\stackrel{·}{p}=-{\nabla }_{q}H+G(u-\delta )-D\stackrel{·}{q}$ . We define the estimation error $z=\tilde{\delta }-\delta$ (see Proposition 1) and the Lyapunov function candidate $W^{\prime }=H_d+\frac{1}{2}{z}^{\mathrm{T}}z$ . Computing the time derivative of $W^{\prime }$ along the trajectories of the closed-loop system and substituting (5) gives

$$\stackrel{·}{W}\prime =-\frac{k_v}{k_m^2}{\left|G^{\mathrm{T}}\stackrel{·}{q}\right|}^2-\frac{1}{k_m}{\stackrel{·}{q}}^{\mathrm{T}}D\stackrel{·}{q}+\frac{1}{k_m}{\stackrel{·}{q}}^{\mathrm{T}}\mathrm{Gz}-\alpha z^{\mathrm{T}}z-z^{\mathrm{T}}\stackrel{·}{\delta }$$

(18)

Introducing the Young’s inequality $\left|z^{\mathrm{T}}\stackrel{·}{\delta }\right|\le ({\left|z\right|}^2/4+{\left|\stackrel{·}{\delta }\right|}^2)$ in (18) and employing a Schur complement argument gives

$$\begin{array}{c}\stackrel{·}{W}\prime \le -\left[\begin{array}{cc}{\left({\stackrel{·}{q}}^{\mathrm{T}}G\right)}^{\mathrm{T}}\hfill & z^{\mathrm{T}}\hfill \end{array}\right]\left[\begin{array}{cc}\frac{k_v}{k_m^2}\hfill & ★\hfill \\ \multicolumn{1}{c}{-\frac{1}{2k_m}}& \alpha -\frac{1}{4}\hfill \end{array}\right]\left[\begin{array}{c}{G}^{\mathrm{T}}\stackrel{·}{q}\hfill \\ \multicolumn{1}{c}{z}\end{array}\right]\hfill \\ -\frac{1}{k_m}{D}_0{\left|\stackrel{·}{q}\right|}^2+{\left|\stackrel{·}{\delta }\right|}^2\hfill \end{array}$$

(19)

If $D_0{\left|\stackrel{·}{q}\right|}^2>k_m{\left|\stackrel{·}{\delta }\right|}^2$ , $k_v>0$ , and $\alpha >(1+k_v)/4k_v$ then $\stackrel{·}{W}\prime \le 0$ and the equilibrium is locally stable. The strict-minimizer claim follows from (14), concluding the proof ▪

Remark 3: Neglecting the disturbances in the adaptation law has detrimental effects on the stability of the equilibrium. Defining the Lyapunov function candidate $H_d$ and computing its time derivative along the trajectories of the closed-loop system gives

$${\stackrel{·}{H}}_d=-\frac{k_v}{k_m^2}{\left|G^{\mathrm{T}}\stackrel{·}{q}\right|}^2-\frac{1}{k_m}{\stackrel{·}{q}}^{\mathrm{T}}D\stackrel{·}{q}+\frac{1}{k_m}{\stackrel{·}{q}}^{\mathrm{T}}G\delta$$

(20)

Following the same procedure as Proposition 2, the corresponding stability condition is in this case $D_0{\left|\stackrel{·}{q}\right|}^2>{\stackrel{·}{q}}^{\mathrm{T}}G\delta$ , which depends on the magnitude of the disturbance rather than on its time derivative, and thus it is typically more stringent (e.g. consider the case of a large but slowly varying disturbance).

3.3 Control law for variable curvature equilibrium

Solving (3) for the case of in-plane unmatched disturbances results in a new set of attainable equilibria ${\Omega }^{{*}^{\prime }}$ , which deviate from CC

$$q_i^{*}-q_j^{*}=({\tilde{\delta }}_j-{\tilde{\delta }}_i)/k$$

(21)

where $1\le i,j\le n$ and $i\ne j$ . Since the system is underactuated, the individual joint angles cannot be regulated to arbitrary values. Nevertheless, it is still possible to achieve the setpoint regulation for the tip rotation ${\theta }_d$ (respectively for the tip coordinate $x_d$ or $y_d$ ) with an appropriate control action.

To account for in-plane unmatched disturbances in the controller (10), the gradient vector ${\nabla }_q{V}_d$ is augmented with the non-conservative generalized forces $\mathrm{\Lambda }(q)$ as in Franco (2019a), resulting in the extended vector ${\nabla }_q{V}_d^{\prime }={\nabla }_q{V}_d+\mathrm{\Lambda }$ , where $\mathrm{\Lambda }$ is computed as

$$G^{\perp }(\tilde{\delta }-M_d{M}^{-1}\mathrm{\Lambda })=0$$

(22)

This approach preserves the solution of the PDE (13) and the expression of $u_{\mathrm{es}}$ (15). Computing $\mathrm{\Lambda }\mathrm{from}$ (22) and enforcing the minimizer conditions (14) for $V_d^{\prime }$ in $q^{*}\in {\Omega }^{*\prime }$ gives

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{\Lambda }}_i=\frac{1}{n{k}_m}((n-1){\tilde{\delta }}_i-\sum _{j=1}^n{\tilde{\delta }}_{j\ne i})}\\ -k_p(\sum _{j=1}^n{q}_j^{*}-{\theta }_d)+\gamma (\theta -{\theta }^{*})\hfill \end{array}$$

(23)

where $\gamma$ is a tuning parameter. Expressing (14) for $V_d^{\prime }$ gives then

$${\nabla }_q{V}_d^{\prime }\left(q_d\right)=\mathrm{\Lambda }\left(q^{*}\right)+{\nabla }_q{V}_d\left(q^{*}\right)=0$$

(24.a)

$${\nabla }_q^2{V}_d\left(q_d\right)=n(k/k_m{)}^{n-1}{k}_p+{\nabla }_q\mathrm{\Lambda }\left(q^{*}\right)>0$$

(24.b)

Condition (24.a) is immediately verified, while (24.b) depends on the structure of the disturbances. Considering the case of constant $f$ for illustrative purposes gives

$${\nabla }_q^2{V}_d\left(q^{*}\right)+{\nabla }_q\mathrm{\Lambda }\left(q^{*}\right)=n(k/k_m{)}^{n-1}(k_p+\gamma )>0$$

(25)

The disturbance compensation term $u^{*}$ for unmatched disturbances is then

$$\begin{array}{c}\multicolumn{1}{c}{u^{*}=G^{†}(\tilde{\delta }-M_d{M}^{-1}\mathrm{\Lambda })}\\ =\frac{1}{n}\sum _{i=1}^n{\tilde{\delta }}_i+k_p{k}_m(\sum _{i=1}^n{q}_i^{*}-{\theta }_d)-k_m\gamma (\theta -{\theta }^{*})\hfill \end{array}$$

(26)

Proposition 3: Consider system (2) under Assumptions 1–5 in closed loop with controller (10), where $u_{\mathrm{es}}$ is defined as (15), $u_{\mathrm{di}}$ is defined as (16), $u^{*}$ is defined as (26), and the disturbances (matched and unmatched) are estimated with (5). Then the equilibrium $(q,p)=(q^{*},0)$ , with $q^{*}\in {\Omega }^{*\prime }$ defined in (21), is stable for all $\alpha >(1+D_0{k}_m)/\left(4D_0{k}_m\right)$ provided that $k_v{\left|G^{\mathrm{T}}\stackrel{·}{q}\right|}^2>k_m^2{\left|\stackrel{·}{\delta }\right|}^2$ . In addition, $(q^{*},0)$ is a strict minimizer of $V_d$ if (24.b) is verified.

Proof: Following closely the structure of Proposition 2 we introduce the function $H_d^{\prime }=\frac{1}{2}{p}^{\mathrm{T}}{M}_d^{-1}p+V_d^{\prime }+\mathcal{C}$ , where $V_d^{\prime }$ is defined locally as $V_d^{\prime }=V_d+\mathrm{\Lambda }{\left(q^{*}\right)}^{\mathrm{T}}(q-q^{*})$ and $\mathcal{C}>0$ is an arbitrary constant that ensures positive definitiveness of $H_d^{\prime }$ . We define the vector of estimation errors $z=\tilde{\delta }-\delta$ and the Lyapunov function candidate $W^{″}=H_d^{\prime }+\frac{1}{2}{z}^{\mathrm{T}}z$ . Computing the time derivative of $W^{″}$ along the trajectories of the closed-loop system, substituting (5) and $M_d=k_{m}M$ gives

$$\stackrel{·}{W}″=-\frac{k_v}{k_m^2}{\left|G^{\mathrm{T}}\stackrel{·}{q}\right|}^2-\frac{1}{k_m}{\stackrel{·}{q}}^{\mathrm{T}}D\stackrel{·}{q}+\frac{1}{k_m}{\stackrel{·}{q}}^{\mathrm{T}}z-\alpha z^{\mathrm{T}}z-z^{\mathrm{T}}\stackrel{·}{\delta }$$

(27)

Introducing the Young’s inequality $\left|z^{\mathrm{T}}\stackrel{·}{\delta }\right|\le ({\left|z\right|}^2/4+{\left|\stackrel{·}{\delta }\right|}^2)$ in (27) and employing a Schur complement argument gives

$$\stackrel{·}{W}″\le -\left[\begin{array}{cc}{\stackrel{·}{q}}^{\mathrm{T}}\hfill & z^{\mathrm{T}}\hfill \end{array}\right]\left[\begin{array}{cc}\frac{D_0}{k_m}\hfill & ★\hfill \\ \multicolumn{1}{c}{-\frac{1}{2k_m}}& \alpha -\frac{1}{4}\hfill \end{array}\right]\left[\begin{array}{c}\stackrel{·}{q}\hfill \\ \multicolumn{1}{c}{z}\end{array}\right]-\frac{k_v}{k_m^2}{\left|G^{\mathrm{T}}\stackrel{·}{q}\right|}^2+{\left|\stackrel{·}{\delta }\right|}^2$$

(28)

Thus, if $k_v{\left|G^{\mathrm{T}}\stackrel{·}{q}\right|}^2>k_m^2{\left|\stackrel{·}{\delta }\right|}^2$ and $\alpha >(1+D_0{k}_m)/\left(4D_0{k}_m\right)$ then $\stackrel{·}{W}\prime \le 0$ and the equilibrium is locally stable. Finally, the strict-minimizer claim follows from (24), concluding the proof ▪

Remark 4: The difference between Proposition 2 and Proposition 3 is due to the dimension of $z$ , which is a scalar in (18), while it is a vector of dimension $n$ in (27). In the latter case, the disturbances also affect the unactuated states; thus, the condition on the parameter $\alpha$ relies on the open-loop damping $D_0$ rather than on the damping injection parameter $k_v$ . The parameter $k_m$ resulting from the kinetic-energy shaping allows additional tuning freedom for $\alpha$ with a given $D_0$ . Notably, the stability conditions provided in Proposition 2 and Proposition 3 are sufficient but not necessary. This is due to the use of the Young’s inequality, which is a conservative condition.

Remark 5: If the control aim corresponds to regulating the tip rotation so that $\sum _{i=1}^n{q}_i^{*}={\theta }^{*}={\theta }_d$ and if $\gamma =0$ then $u^{*}$ in (26) coincides with $u^{*}=G^{†}\delta$ , which refers to the matched disturbances case. Otherwise, the constant term $k_p({\theta }^{*}-{\theta }_d)$ and the linear term $\gamma (\theta -{\theta }^{*})$ contribute to the control input. In addition, taking $\gamma (q)$ dependent on the states results in a nonlinear control law. In this respect, expression (26) is more general than (10). Finally, the regulation of the tip position can be achieved in the same way as the angle regulation by employing an appropriate $q^{*}\in {\Omega }^{*\prime }$ such that $x(q^{*})=x_d$ .

Remark 6: Introducing the tip rotation $\theta =G^{\mathrm{T}}\stackrel{·}{q}$ and its desired value ${\theta }_d=n{q}_d$ in (15) reveals that the components of the control law $u_{\mathrm{es}}$ and $u_{\mathrm{di}}$ are only dependent on $\theta ,\stackrel{·}{\theta }$ . However, the adaptation law (5) depends on the individual positions of the virtual joints, which need to be computed with an appropriate observer. To simplify the controller implementation, a compact version of (26) is constructed approximating ${\nabla }_{q}H\cong {\nabla }_{q}V$ , which is reasonable for slow movements of the manipulator. The complete control law becomes

$$\begin{array}{c}u=\frac{k}{n}\theta -k_p{k}_m(\theta -{\theta }_d)-\frac{k_v}{k_m}\stackrel{·}{\theta }+u^{*}\hfill \\ u^{*}=\alpha \int (-\frac{k}{n}\theta -\frac{D_0}{n}\stackrel{·}{\theta }+u-G^{†}\hat{\delta })\mathrm{dt}+k_p{k}_m({\theta }^{*}-{\theta }_d)+k_m\gamma (\theta -{\theta }^{*})\hfill \end{array}$$

(29)

The integral term in (29) corresponds to the cumulative disturbance estimate $G^{†}\hat{\delta }$ with a constant $f$ (see Remark 2), while the second and the third terms are the same as in (26). As a result, only the cumulative effect of the disturbances (i.e. $G^{†}\delta =G^{†}{\delta }_{0}f$ ) is estimated, rather than the values of the unknown parameters ${\delta }_0$ . Since (29) only depends on $\theta$ and $\stackrel{·}{\theta }$ , in this case Assumption 1 can be removed. The stability of the equilibrium (21) for system (2) in closed loop with controller (29) is discussed in the following corollary.

Corollary 1: Consider system (2) under Assumptions 2–5 in closed loop with controller (29). Then the equilibrium $(q,p)=(q^{*},0)$ is locally stable for some $\alpha >0$ such that

$$\mathcal{A}=\left[\begin{array}{cc}\frac{D_0-\alpha {\kappa }_2{m}_T}{k_m}\hfill & ★\hfill \\ \multicolumn{1}{c}{-\frac{1}{2}(\frac{1}{k_m}+{\alpha }^2{\kappa }_2{m}_T)}& \alpha -\frac{1}{2}\hfill \end{array}\right]>0$$

(30)

and provided that $\frac{k_v}{k_m^2}{\left|\stackrel{·}{\theta }\right|}^2>({\left|\stackrel{·}{\delta }\right|}^2+\alpha {\kappa }_1{m}_T{\left|\stackrel{·}{q}\right|}^4)$ , for some ${\kappa }_1,{\kappa }_2>0$ .

Proof: Computing the time derivative of the estimation error $z=\hat{\delta }+\beta \left(p\right)-\delta$ and substituting $\stackrel{·}{\hat{\delta }}=\alpha (-{\nabla }_{q}V+\mathrm{Gu}-D{\nabla }_{p}H-{\hat{\delta }}_0)$ results in the following error dynamics

$$\stackrel{·}{z}=-\alpha z+\alpha ({\nabla }_q\left(p^{\mathrm{T}}{M}^{-1}p\right)-\alpha p)-\stackrel{·}{\delta }$$

(31)

Substituting $\tilde{\delta }=\hat{\delta }$ in (26) recovers (29). Computing the time derivative of $W^{″}$ as in (27) and substituting (29) and (31) gives

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{W}″=-\frac{k_v}{k_m^2}{\left|G^{\mathrm{T}}\stackrel{·}{q}\right|}^2-\frac{1}{k_m}{\stackrel{·}{q}}^{\mathrm{T}}D\stackrel{·}{q}+\frac{1}{k_m}{\stackrel{·}{q}}^{\mathrm{T}}(z+\alpha p)}\\ -\alpha z^{\mathrm{T}}z-z^{\mathrm{T}}\stackrel{·}{\delta }+z^{\mathrm{T}}\alpha ({\nabla }_q\left(p^{\mathrm{T}}{M}^{-1}p\right)-\alpha p)\hfill \end{array}$$

(32)

We recall that ${\nabla }_p{H}^{\mathrm{T}}=\stackrel{·}{q}$ and $G^{\mathrm{T}}{\nabla }_{p}H=\stackrel{·}{\theta }$ . In addition, motivated by the structure of $M$ (see Appendix A), we can write ${\nabla }_q\left(p^{\mathrm{T}}{M}^{-1}p\right)\le {\kappa }_1{m}_T{\left|\stackrel{·}{q}\right|}^2$ and $\left|p\right|\le {\kappa }_2{m}_T\left|\stackrel{·}{q}\right|$ for some positive constants ${\kappa }_1,{\kappa }_2$ . Rewriting (32) and regrouping the common terms gives

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{W}\prime \prime \le -\frac{k_v}{k_m^2}{\left|\stackrel{·}{\theta }\right|}^2-\frac{1}{k_m}{D}_0{\stackrel{·}{q}}^{\mathrm{T}}\stackrel{·}{q}+\frac{1}{k_m}{\stackrel{·}{q}}^{\mathrm{T}}\left(\alpha {\kappa }_2{m}_T\stackrel{·}{q}\right)+\frac{1}{k_m}{\stackrel{·}{q}}^{\mathrm{T}}z}\\ -\alpha z^{\mathrm{T}}z+\left|z^{\mathrm{T}}\stackrel{·}{\delta }\right|+z^{\mathrm{T}}\alpha ({\kappa }_1{m}_T{\left|\stackrel{·}{q}\right|}^2+\alpha {\kappa }_2{m}_T\stackrel{·}{q})\hfill \end{array}$$

(33)

Introducing the Young’s inequalities $\left|z^{\mathrm{T}}\stackrel{·}{\delta }\right|\le ({\left|z\right|}^2/4+{\left|\stackrel{·}{\delta }\right|}^2)$ and $\left|z\right|{\left|\stackrel{·}{q}\right|}^2\le ({\left|z\right|}^2/4+{\left|\stackrel{·}{q}\right|}^4)$ in (33) and employing a Schur complement argument gives

$$\stackrel{·}{W}\prime \prime \le -\left[\begin{array}{cc}{\stackrel{·}{q}}^{\mathrm{T}}\hfill & z^{\mathrm{T}}\hfill \end{array}\right]\mathcal{A}\left[\begin{array}{c}\stackrel{·}{q}\hfill \\ \multicolumn{1}{c}{z}\end{array}\right]-\frac{k_v}{k_m^2}{\left|\stackrel{·}{\theta }\right|}^2+{\left|\stackrel{·}{\delta }\right|}^2+\alpha {\kappa }_1{m}_T{\left|\stackrel{·}{q}\right|}^4$$

(34)

where $\mathcal{A}$ is given in (30). Thus, $\stackrel{·}{W}″\le 0$ and the equilibrium is locally stable if (30) is verified and provided that $\frac{k_v}{k_m^2}{\left|\stackrel{·}{\theta }\right|}^2>{\left|\stackrel{·}{\delta }\right|}^2+\alpha {\kappa }_1{m}_T{\left|\stackrel{·}{q}\right|}^4$ , which concludes the proof ▪

Note finally that, if $m_T<<1$ , inequality (30) can be approximated neglecting the terms $m_T^2$ , which yields the condition

$$(2D_0-{\mathrm{\Delta }}_1)/{\mathrm{\Delta }}_2+1/6<\alpha <(2D_0+{\mathrm{\Delta }}_1)/{\mathrm{\Delta }}_2+1/6$$

(35)

where ${\mathrm{\Delta }}_1=2\sqrt{D_0^2-{\mathrm{\Delta }}_2(4D_0+3/k_m)/12}$ and ${\mathrm{\Delta }}_2=\left(6{\kappa }_2{m}_T\right)$ . Inequality (35) admits real solutions provided that $D_0>{\kappa }_2{m}_T+3\sqrt{{\kappa }_2{m}_T/k_m}$ .

4. Results of simulations and experiments

The proposed controllers are implemented for the soft continuum manipulator described by Garriga-Casanovas et al. (2018), by employing a $3$ DOF rigid-link model with link lengths $l_0=l_3=0.125l_T$ and $l_1=l_2=0.375l_T$ , where $l_T$ is the total length of the device. These values were chosen as a result of an optimization routine to minimize the difference between the potential energy of the rigid-link model and that of the computer-aided design (CAD) model computed with FE simulations. A point mass $m_i$ is located at the mid-point of each link, and is proportional to its length so that the sum of the masses is equal to the mass of the manipulator $m_T$ . It must be noted that the proposed controllers do not depend on this particular choice.

4.1 Simulations

Simulations were conducted in MATLAB, by employing the solver $\mathrm{ode}23$ with the following parameters representative of an ideal manipulator: $k=5,l_T=0.1,m_T=1.5,D_0=0.01,k_v=4,k_p=0.01,k_m=20,\alpha =10$ , and alternatively $\alpha =20$ , chosen according to Propositions 2 and 3. The value of desired tip rotation is ${\theta }_d=\pi /2$ , which corresponds to equal joint angles $q_d=\pi /6$ in the case of CC equilibrium. The control input was computed with (10), (26) and the disturbances were estimated with (5) using full state feedback. Three types of disturbance were employed for illustrative purposes: matched disturbances consisting of a constant bending moment $\delta =-1\mathrm{Nmm}$ (see Figure 2); matched and unmatched disturbances consisting of a tip force $P=-10\mathrm{N}$ in the $y$ direction (see Figure 3); and the weight of the manipulator amplified by a factor 10 (see Figure 4). Notably, the force $P$ corresponds to the weight of a mass attached at the tip of the manipulator, as might occur during a manipulation task (see Section 4.2). A baseline proportional–integral–derivative (PID) controller defined as $u=K_p({\theta }_d-\theta )-K_v\stackrel{·}{\theta }+K_i\int ({\theta }_d-\theta )\mathrm{dt}$ was implemented for comparison purposes. Note that, because of the close similarity between the structure of the PID and of the controller (10), (26), the following equalities hold: $K_p=k_p{k}_m,K_v=k_v/k_m,$ while $K_i$ corresponds to $\alpha$ . Consequently, the PID parameters are chosen as $K_p=0.2,K_v=0.2,K_i=10$ , or $K_i=20$ . Thus, the corresponding tuning parameters are identical in both controllers. Low values of $k_p$ and corresponding values of $K_p$ are chosen to avoid reducing the compliance of the manipulator in closed loop (DellaSantina et al., 2017).

Fig. 2.

Simulation results for tip moment $\delta =-1$ : (a) time history of the tip rotation with controller (10) and (26) using $k_v=4,k_p=0.01,k_m=20$ ; (b) corresponding control input; (e) disturbance estimate ${\tilde{\delta }}_0$ ; and (b) time history of the tip rotation with the baseline proportional–integral–derivative and equivalent tuning parameters $K_p=0.2,K_v=0.2$ ; (d) corresponding control input; (f) contribution of the integral action.

Fig. 3.

Simulation results for tip force $P=-10$ : (a) time history of the tip rotation with controller (10), (26) using $k_v=4,k_p=0.01,k_m=20$ ; (c) corresponding control input; (e) disturbance estimate ${\tilde{\delta }}_0$ ; and (b) time history of the tip rotation with the baseline proportional–integral–derivative and $K_p=0.2,K_v=0.2$ ; (d) corresponding control input; (f) contribution of the integral action.

Fig. 4.

Simulation results with manipulator weight amplified by a factor of 10: (a) time history of the tip rotation with controller (10), (26) using $k_v=4,k_p=0.01,k_m=20,\alpha =20$ ; (b) time history of the tip rotation with the baseline proportional–integral–derivative and $K_p=0.2,K_v=0.2,K_i=2$ .

In the case of matched disturbances both controllers correctly achieve the regulation goal, which corresponds to CC equilibrium (see Figure 2). Their smooth and slow response is due to the presence of open-loop damping and to the use of small gains. In particular, the adaptive controller (10), (26) results in a slower and smoother response compared to the PID baseline with the same tuning. Employing $\alpha =20$ , the disturbance estimate reaches the correct value more quickly without drastically changing the closed-loop response. Conversely, the PID with $K_i=20$ results in a larger overshoot (8.7% with $K_i=10$ , and 22.5% with $K_i=20$ ). This is due to the contribution of the integral term, which shows a similar trend and also converges to a different value compared to the case with $K_i=10$ . Employing larger values of $k_p$ and $K_p$ (e.g. $k_p=1,k_m=2$ , and $K_p=k_p{k}_m=2$ ) results in a faster response for both controllers (see Figure C1 in Appendix C). Also in that case, however, the baseline PID shows a noticeable overshoot with $K_i=20$ . Finally, setting $\alpha =K_i=0$ results in large steady-state errors because the disturbances are not compensated.

The simulation results for the case of unmatched disturbances are shown in Figure 3. Also in this case, the regulation goal is correctly achieved by both controllers. The system response with the adaptive controller (10), (26) is similar for different values of $\alpha$ and is comparable to that of Figure 2 even though the disturbance has a different magnitude. With small values of $K_i$ (e.g. $K_i=2$ ) the performance of the PID becomes similar to that of controller (10), (26). However, the system response changes more drastically with larger $K_i$ , resulting in noticeable overshoot (11% with $K_i=10$ , and 25.4% with $K_i=20$ ). This suggests that the performance of the proposed controller is less sensitive to the gain governing the integral action compared to the PID. In practice, this can be advantageous since it could avoid time consuming tuning procedures. Finally, the adaptation law (5) correctly estimates the force $P$ for different values of $\alpha$ , which could serve the purpose of estimating interactions with the environment. Instead, the contribution of the integral action in the PID changes depending on $K_i$ .

The simulation results for a different type of unmatched disturbances representative of the weight of the manipulator amplified by a factor 10 are shown in Figure 4 considering two opposite mounting positions (i.e. gravity vector pointing up or down with respect to the robot frame). In both cases the manipulator at rest is horizontal while the weight of each segment acts in the vertical direction. The same tuning is employed for the parameters $k_v,k_p,k_m,K_p,K_v$ as in Figure 3. The values $\alpha =20$ and $K_i=2$ are used in an attempt to match the response of the controller (10), (26) with that of the PID (see Figure 3). The results show that the energy-shaping controller leads to a similar response for different directions of the load. Conversely, the performance of the PID changes more drastically. Increasing $\alpha$ results in a similar response, while increasing $K_i$ leads to a faster convergence with a noticeable overshoot for the “down” orientation, since gravity acts in the direction of motion.

The time history of the position $q$ for the previous operating conditions is shown in Figure 5. In the presence of matched disturbances, the position deviates from the CC configuration during the transient. This confirms that, although the regulation of the tip rotation might appear to be a 1 DOF problem in proximity of the equilibrium, a 1 DOF model is not appropriate outside of quasi-static operating conditions. Instead, in the presence of unmatched disturbances the equilibrium is no longer a CC configuration. However, it is still possible to stabilize the equilibrium $\theta ={\theta }_d$ or alternatively $x=x_d$ employing an appropriate value of ${\theta }^{*}$ in (26). The value of ${\theta }^{*}$ is typically different from ${\theta }_d$ and can be either computed analytically from a kinematic model or numerically based on the time history of $\theta$ and $x$ .

Fig. 5.

Simulation results with controller (10), (26) and $k_v=4,k_p=0.01,k_m=20,\alpha =20$ : (a) time history of the position q with $\delta =-1$ ; (c) corresponding terms of the control input. And (b) time history of the position with $P=-10$ ; (d) corresponding terms of the control input.

The system response with controller (29) employing the same tuning parameters is shown in Figure 6 for both disturbances and for different values of $n$ . If the parameter $n$ in (29) matches the number of DOFs of the model, the resulting response is very close to that of the controller (10), (26). Using a different $n$ in (29) results in noticeable but limited differences that do not compromise performance (e.g. overshoot). In addition, the regulation goal is still achieved in a similar way for disturbances of different magnitude. A notable difference is that with the controller (10), (26) and with the adaptation law (5) the unknown parameter ${\delta }_0$ is correctly estimated (see Figures 2 and 3). Instead, controller (29) compensates for the cumulative effect of the disturbances but does not provide an estimate of the parameter ${\delta }_0$ .

Fig. 6.

Simulation results: (a) time history of the tip rotation with $\delta =-1$ and controller (29) with $k_v=4$ , $k_p=0.01$ , $k_m=20$ , $\alpha =20$ ; (b) corresponding results with $P=-10$ .

4.2 Experiments

4.2.1. Experimental setup

The controller (29) was tested on the soft continuum manipulator described by Garriga-Casanovas et al. (2018), shown in Figures 7(a) and (b), by employing the setup illustrated in Figure 7(c). The controller (29) was employed for the experiments since it does not require a kinematic observer to compute the states of the virtual joints (see Remark 6). The rotation $\theta ,\phi$ and the position coordinates $x,y,z$ of a 5 DOF sensor mounted at the tip of the manipulator are measured with an electromagnetic tracking system (Aurora, NDI, Canada, 0.70 RMS) at 40 Hz sampling rate. A MATLAB script records the sensor readings and computes the control input $u$ . The corresponding pressure values are then computed from (1) and are communicated to a digital microcontroller (mbed NXP LPC1768) via a serial interface. The latter communicates control voltages in the range of 0–3.3 V to proportional pressure regulators (Tecno Basic, Hoerbiger, Germany), which supply the chambers of the manipulator. The prototype employed in the experiments measures 30 mm in length and 6 mm in diameter, and its mass is approximately 1.5 g. The stiffness of the structure ( $k=4\mathrm{Nmm}/\mathrm{rad}$ ) was determined from FE simulations and experiments. The latter confirmed that the tip rotation at equilibrium varies almost linearly with the internal pressure up to 2.5 bar, and thus Equation (1.b) suggests that employing a constant parameter $k$ is a reasonable approximation in this case. The damping $D_0$ was estimated experimentally by applying a step command of 1 bar in different chambers and recording the values of $\theta ,\stackrel{·}{\theta }$ . Computing $D_0$ from (2) evaluated at the maximum velocity and averaging the values for different chambers resulted in $D_0=0.03\mathrm{Nms}$ . The following numerical values of the controller parameters were thus used for the experiments: $k=4,k_v=1,k_p=0.01$ or $k_p=0.025,k_m=20,\alpha =10$ , and alternatively $\alpha =20$ . Both values of $\alpha$ fulfil the conditions of Corollary 1 for all ${\kappa }_2<0.6$ (e.g. $1<\alpha <21$ if ${\kappa }_2=0.6$ ). In particular, a small value of $k_p$ is employed to avoid high feedback gains in spite of using larger values of $k_m$ . In turn, employing a large value of $k_m$ allows a greater variation of $\alpha$ according to condition (35). Given a desired value of the tip rotation ${\theta }_d$ , different orientations of the bending plane can be achieved by setting $P_1=\mu P_2$ with appropriate values of $\mu$ . Considering the small size of the internal chambers of the manipulator ( $<<1\mathrm{ml}$ ), and the comparatively large flow rate that the digital pressure regulators can supply ( $>>10\mathrm{l}/\mathrm{s}$ ) as well as their fast response ( $\cong 10\mathrm{ms}$ ), the pressure dynamics of the system was neglected in this work. The motivation of this choice is that in this case the pressure wave propagates considerably faster than the reaction time of the pressure regulator (0.3 ms propagation time for a 0.1 m long pipe). The pressure dynamics would not be negligible with longer supply pipes (e.g. > $10\mathrm{m}$ ), which are required in specific applications, such as MRI-compatible robotics (Franco et al., 2016). This specific scenario will be investigated as part of our future work. The disturbances acting on the manipulator during the experiments include the following: the weight of the manipulator itself and of the sensor, which is mounted at the tip; uncertainties in stiffness and damping parameters, which might not be uniform due to manufacturing inaccuracies and might be affected by pressurization; uncertainty in the bending moment generated by the applied pressure. Finally, an additional disturbance was included by attaching a small mass at the tip of the manipulator.

Fig. 7.

(a) Side view of the prototype of soft continuum manipulator from (Garriga-Casanovas et al., 2018); (b) detail view of the internal chambers; (c) Test setup.

4.2.2. Experimental results

The controller (29) correctly achieved the regulation goal ${\theta }^{*}={\theta }_d=\pi /6$ consistently for different orientations of the bending plane (Figure 8). This result confirms that the control of the manipulator in three dimensions is a direct extension of the two-dimensional (2D) problem in the absence of out-of-plane disturbances. The system response is slower compared to the simulations because of the higher damping of the prototype. The effect of the parameter $\alpha$ and of the parameter $k_p$ is illustrated in Figure 9, which shows a close similarity to the simulation results. In particular, larger values of $\alpha$ and of $k_p$ result in a faster response.

Fig. 8.

Experimental results: (a) tip rotations $\theta$ and $\phi$ with controller (29) using $\alpha =10$ and $P_1=P_2$ ; (c) corresponding control input and disturbance compensation $u^{*}=G^{†}\hat{\delta }$ . And (b) tip rotations $\theta$ and $\phi$ with controller (29) using $\alpha =10$ and $P_1=2P_2$ ; (d) corresponding control input and disturbance compensation.

Fig. 9.

Experimental results: (a) tip rotations $\theta$ with controller (29) and $k_p=0.01$ for different values of $\alpha$ ; (c) corresponding control input. And (b) tip rotation $\theta$ with $\alpha =10$ for different values of $k_p$ ; (d) corresponding control input.

The repeatability of the closed-loop system response was assessed for the two regulation goals ${\theta }_d=\pi /6$ and ${\theta }_d=\pi /8$ , each repeated five times. The standard deviation of the tip rotation $\theta$ was evaluated at time $t=1$ s and at $t=2$ s and remained below 0.01. Finally, the controller (29) was tested with $u^{*}=0$ , resulting in a large error on the tip rotation $\theta$ . In this case, incrementing $k_p$ in an attempt to reduce the error increases the closed-loop stiffness of the system and might lead to instability since the disturbances are not accounted for (see Remark 3).

A mass of 3 g, which is approximately twice the mass of the manipulator, was attached at its tip to introduce additional disturbances due to gravity (Figure 10). The results show that the controller (29) with $\alpha =20$ correctly achieves the regulation goal ${\theta }^{*}={\theta }_d=\frac{\pi }{6}$ . Similar to the simulation study, the baseline PID was employed for comparison purposes using the tuning parameters $K_p=k_p{k}_m=0.2,K_v=k_v/k_m=0.05,$ which correspond to those of the energy-shaping controller (Figure 10). In this case, the value $K_i=10$ was employed for the PID since the system becomes unstable with $K_i=\alpha =20$ . Although the PID results in a faster response, it also leads to high overshoot and oscillations, which increase in amplitude with larger $K_i$ . Employing lower values of $K_i$ the system response becomes more similar for both controllers (see Figure C2 in Appendix C, which shows results for a different setpoint and different tip forces). This confirms the validity of linear controllers in quasi-static operating conditions (Bieze et al., 2018; Zhang et al., 2016). Comparing the control law (29) with the PID highlights a close similarity and suggests that their different performance could be due to the following reasons: (a) the PID lacks a feed-forward term that accounts for the structural stiffness of the manipulator; (b) the disturbance compensation in (29) accounts for the damping and for the stiffness of the manipulator within the adaptive law, which are ignored in the PID.

Fig. 10.

Experimental results with 3 g mass attached to the tip of the manipulator from the start of the test: (a) tip rotation $\theta$ for controller (29) with $\alpha =20$ and for PID with $K_i=10$ ; (c) corresponding control input; (d) corresponding disturbance estimate. And (b) schematic of the loading condition.

A further set of experiments was conducted attaching the 3 g mass to the tip of the manipulator during the test (Figure 11). Also in this case, controller (29) results in a smooth response and the increment in control input after the mass is attached to the manipulator corresponds to the increase in the cumulative disturbance estimate $G^{†}\hat{\delta }$ . Comparing Figure 10(a) and Figure 11(a) confirms that the system response with the controller (29) is similar during the first 10 seconds even though the disturbances are different and even if different values of $\alpha$ are employed. Instead, with the tuning employed the PID results in a faster response and higher overshoot before the mass is attached to the tip of the manipulator, since the actuator encounters less resistance. The onset of the disturbances is indicated on Figure 11 by vertical dashed lines. In particular, the mass is attached manually to the tip of the manipulator and this operation requires a few seconds. This variability contributes to the different deviation between $\theta$ and ${\theta }_d$ after the disturbance onset. Afterwards, the controller (29) brings the tip rotation $\theta$ back to the desired value while the PID results in small oscillations around the setpoint. A faster response is achieved by employing a larger value of $k_p$ in (29) and a corresponding value of $K_p$ in the baseline PID (see Figure C3 in Appendix C). In the latter case, the PID results in higher overshoot and oscillations after the disturbance onset, while controller (29) still achieves a smooth response.

Fig. 11.

Experimental results with 3 g mass attached to the tip of the manipulator during the test: (a) tip rotation $\theta$ for controller (29) with $\alpha =10$ , $k_p=0.01$ and for PID with $K_i=\alpha =10$ , $K_p=k_p{k}_m=0.2,$ $K_v=k_v/k_m=0.05$ ; (c) corresponding control input; (d) disturbance estimate. And (b) picture of the loading condition. The disturbance onset is indicated with vertical dashed lines.

Note finally that a very similar response of the controller (29) and of the PID can be achieved by reducing $K_i$ (e.g. $K_i=5$ ) and by increasing $n$ (e.g. $n=30$ ), as shown in Figure 12. In this case, the feed-forward term that accounts for the structural stiffness of the manipulator in the control law (29) is reduced. The system response is similar for both loading conditions considered (i.e. tip mass introduced either at the beginning of the test or when the rotation reaches the setpoint) and the PID only shows a small overshoot (approximately 4%). The performance of both controllers could be further improved with an optimized parameter tuning, which we intend to explore as part of future work.

Fig. 12.

Experimental results with different loading conditions: (a) tip rotation $\theta$ for controller (29) with $\alpha =20$ and $n=30$ (b) tip rotation $\theta$ for PID with $K_i=5$ , $K_p=k_p{k}_m=0.5$ , $K_v=k_v/k_m=0.05$ . The constant load refers to the same test condition as Figure 10, while the variable load denotes the same test condition as Figure 11.

5. Conclusions

An adaptive energy-shaping controller for a class of soft continuum manipulators subject to in-plane disturbances was designed by employing a rigid-link model and a port-Hamiltonian formulation. The rigid-link model approximates the dynamics of the soft continuum manipulator and simplifies the controller design. The effects of the discrepancies between the model and the real system are treated as disturbances and are compensated adaptively. Closed-form expressions of the control law for an $n$ DOF model, where $n$ can be arbitrarily large, were presented and stability conditions were discussed. The effects of external disturbances on the equilibrium were highlighted, confirming that CC is not achievable with unmatched disturbances. The effectiveness of the controller was demonstrated with simulations and with experiments on a soft continuum and inextensible manipulator with pneumatic actuation. A comparison with a baseline PID highlighted that the energy-shaping controller results in a smoother but slower response when employing the same tuning parameters. The results indicate that, while similar performance can be achieved with both controllers for a given loading condition, the energy-shaping controller appears less sensitive to the tuning of the integral action.

Although the controller presented in this work was implemented for a specific soft continuum manipulator consisting of one inextensible segment, we believe that the proposed approach is of wider applicability. Future work aims to extend the results to out-of-plane disturbances and to investigate the control of manipulators consisting of multiple segments arranged in series. In addition, we intend to include the pressure dynamics in the port-Hamiltonian model and to investigate the control of soft manipulators with long supply pipes within the context of a medical application. Finally, we intend to apply the proposed control approach to soft continuum manipulators that employ different types of actuation.