Nonlinear energy-based control of soft continuum pneumatic manipulators

This paper investigates the model-based nonlinear control of a class of soft continuum pneumatic manipulators that bend due to pressurization of their internal chambers and that operate in the presence of disturbances. A port-Hamiltonian formulation is employed to describe the closed loop system dynamics, which includes the pressure dynamics of the pneumatic actuation, and new nonlinear control laws are constructed with an energy-based approach. In particular, a multi-step design procedure is outlined for soft continuum manipulators operating on a plane and in 3D space. The resulting nonlinear control laws are combined with adaptive observers to compensate the effect of unknown disturbances and model uncertainties. Stability conditions are investigated with a Lyapunov approach, and the effect of the tuning parameters is discussed. For comparison purposes, a different control law constructed with a backstepping procedure is also presented. The effectiveness of the control strategy is demonstrated with simulations and with experiments on a prototype. To this end, a needle valve operated by a servo motor is employed instead of more sophisticated digital pressure regulators. The proposed controllers effectively regulate the tip rotation of the prototype, while preventing vibrations and compensating the effects of disturbances, and demonstrate improved performance compared to the backstepping alternative and to a PID algorithm.


Introduction
Soft continuum manipulators have been attracting increasing attention in the research community because, unlike conventional rigid mechanisms, they can interact safely with uncertain and delicate objects [34]. This property is a consequence of the structural compliance of soft actuators and makes them ideally suited for several applications, including surgery and rehabilitation [33]. Differently from conventional rigid mechanisms where actuation is concentrated at the joints, soft actuators distribute their force in a continuous fashion over their structure. An iconic example of soft actuator is the pneumatic artificial muscle (PAM), which provides actuator compliance that resembles that of physiologi-cal systems [21]. Pneumatic actuation has been one of the preferred strategies for soft continuum manipulators due to its high power density [41], ease of miniaturization, and affordability, which is highly desirable in low-income countries. The flexible micro actuator (FMA) has been one of the first designs to achieve bending on any plane due to pressurization of its internal chambers [37] so that a manipulator can be constructed with a single FMA, thus it has inspired various subsequent versions [17]. Other actuation strategies include dielectric elastomers [26,45], cables [32], and shapememory alloy [23,42]. In general, not all degrees of freedom (DOFs) in soft continuum manipulators are actuated. This attribute, together with the presence of unknown external forces, which are common in realworld applications, makes the control of soft continuum actuators and manipulators a challenging task [39].
Early attempts to control soft continuum manipulators have been relying solely on kinematic models, which is appropriate in quasi-static conditions. If external forces are negligible, the kinematic model can be further simplified by introducing the assumption of constant curvature (CC) or piece-wise constant curvature (PCC) so that the configuration of the manipulator can be uniquely defined by a reduced number of DOFs. Alternative quasi-static approaches that do not rely on CC or PCC assumptions include fast finite element methods [27], screw theory [43], and Cosserat rod theory [3]. Conversely, a dynamical model of the system is necessary for control purposes in case of fast movements [40], and model reduction techniques are required to render the control problem tractable. However, the simplifications introduced by model reduction methods should be accounted for in the controller design and in the corresponding stability analysis. Notable approaches rely on the CC and PCC assumptions [10], on extensions of beam theory [35], or employ rigid-link models [18]. An accurate continuum-based multibody modeling method for motion and shape control of soft robots was recently presented in [36], however closed-loop control was not demonstrated experimentally. Implementations of classical model-based controllers for soft continuum manipulators include model predictive control (MPC), sliding mode control (SMC) [2,4], and feedback linearization [40]. However, classical high-gain controllers such as SMC have been shown to increase the closed-loop stiffness of the system potentially reducing the benefits of employing soft actuators. This observation has motivated research on controllers that combine feedback and feedforward actions [8]. Nevertheless, the former approaches do not consider the effect of unknown external forces or of pressure dynamics and treat the system as fully actuated. Nonlinear controllers and adaptive observers have recently been investigated to compensate the effects of disturbances and model uncertainties in soft actuators, including PAMs [6,46,47], and in manipulators actuated by shape-memory alloy [24]. In our earlier work, we have introduced an energy-based approach for the control of soft continuum planar manipulators by employing a port-Hamiltonian formulation and by including either an adaptive algorithm [12] or an integral action [13] to compensate the effect of disturbances. These results were extended to soft continuum manipulators in 3D space [14], however they did not account for the pressure dynamics in the system, since digital pressure regulators were employed for the actuation. Digital pressure regulators are ideally suited to control soft actuators [12,46,47] since they automatically adjust the air outflow in order to maintain the output pressure at the prescribed value. If the dynamics of the digital pressure regulator is sufficiently fast compared to that of the soft manipulator, it is then reasonable to neglect the pressure dynamics thus greatly simplifying the controller design. Nevertheless, relying on digital pressure regulators makes the control algorithm less general and increases the cost of soft robotic systems thus potentially limiting the access to these technologies in low-income countries.
In this work, we investigate the energy-based control of soft continuum pneumatic manipulators that can bend on any plane due to internal pressurization. A rigid-link model and a port-Hamiltonian formulation are employed to describe the closed loop system dynamics, including the pressure dynamics under some realistic simplifying assumptions. The port-Hamiltonian formulation is ideally suited to model underactuated mechanical systems including soft robots [31]. Control laws are constructed for the proposed dynamical model by employing the Interconnection and damping assignment Passivity based control methodology (IDA-PBC) [28]. IDA-PBC is widely applicable to nonlinear mechanical systems and offers an interpretation of the control action in terms of mechanical energy. In addition, adaptive observers constructed with the Immersion and Invariance (I&I) methodology [1] are included in the control laws to compensate the effect of disturbances and model uncertainties.
We have chosen I&I observers since they can be integrated with IDA-PBC control laws for underactuated mechanical systems [16]. The resulting control strategy differs from other approaches, including our earlier works [12][13][14] in that it accounts for the pressure dynamics and it does not require digital pressure regulators. Hence, the proposed controllers are more general and more widely applicable. In summary, the main contributions of this work include the following points.
1. A controller design procedure that employs the IDA-PBC methodology in a multi-step fashion to account for the pressure dynamics under some realistic simplifying assumptions. To the best of our knowledge, this is the first attempt to design an energybased controller for soft continuum manipulators that explicitly accounts for the pressure dynamics of the pneumatic actuation. 2. New nonlinear control laws that include adaptive observers to compensate the effect of unknown disturbances and model uncertainties. A control law is initially constructed for a model representative of planar soft continuum manipulators, and it is subsequently extended to soft continuum manipulators in 3D. 3. Stability conditions and their relationship with the controller parameters using a Lyapunov approach. An interpretation of the control action in terms of mechanical energy of the closed-loop system is also provided. 4. The performance of the proposed controllers is assessed with simulations and with experiments on a prototype, and it is compared with two alternatives: (i) an energy-based control law constructed using a backstepping procedure; (ii) a PID algorithm. To this end, a low-cost solution consisting of a needle valve operated by a servo DC motor and a pressure sensor has been employed instead of a digital pressure regulator.
The rest of this paper is organized as follows: the dynamical model of the system is introduced in Sect. 2; the controller design for planar soft continuum pneumatic manipulators is outlined in Sect. 3; the controller design for soft continuum pneumatic manipulators in 3D is detailed in Sect. 4; the results of simulations and experiments are presented in Sect. 5; concluding remarks are discussed in Sect. 6.

Dynamical model of soft continuum pneumatic manipulators
The class of systems considered in this work includes the FMA [37] and the soft continuum pneumatic manipulator described in [17]. Both designs have three internal chambers spaced at 120 • , and the latter design has an inextensible thread along the central axis to prevent elongation, and a second thread wound around its external cylindrical surface to prevent radial expansion. The control input corresponds to the air flow into the chambers. In the experiments, the internal chambers are supplied with a variable flow restrictor consisting of a needle valve operated by a servo motor, while a fixed flow restrictor serves as exhaust valve (see Fig. 1). As the servo motor progressively opens the needle valve, the flow rate to the internal chamber increases resulting in a pressure increment. Conversely, when the servo motor closes the needle valve, the chamber is vented through the exhaust valve and the internal pressure drops. The pressures P 1 , P 2 , P 3 in the internal chambers cause the manipulator to bend on a plane according to Suzumori et al. [37,38], such that at equilibrium and without external forces we have where θ is the tip rotation on the bending plane, ϕ is the orientation of the bending plane with respect to a fixed reference frame, and k is the structural stiffness of the manipulator. Since for the design [17] pressurization does not cause twist, this effect is not considered here. The dynamics of the soft continuum manipulator can be approximated with a rigid-link model that has n virtual elastic pin joints on the bending plane and n virtual elastic pin joints in the direction orthogonal to the bending plane. The orientation of the bending plane corresponds to joint 2n + 1 [14] (see Fig. 1a). This approach is similar to Godage et al. [18] and it is based on the pseudo-rigid-body model [44] which allows approximating the force/deflection relationship of a flexible mechanism by introducing virtual elastic joints. Indicating with q i the angle of joint i with respect to the previous link yields θ = n i=1 q i on the bending plane, γ = 2n+1 i=n+1 q i in the direction orthogonal to the bending plane, while ϕ = q 2n+1 and ω = 2n i=n+1 q i is the rotation outside the bend-ing plane. The rigid-link model has (2n + 1) DOFs and one control input for each pressurized chamber, thus it is underactuated. The dynamics of the rigid-link model can be expressed in port-Hamiltonian form as a function of its mechanical energy H 0 = 1 2q T Mq + , where M(q) = M T (q) > 0 is the inertia matrix [28]. The potential elastic energy (q) = k 2 2n i=1 q 2 i does not depend on q 2n+1 or on k since ϕ is only defined by the chamber geometry and by P 1 , P 2 , P 3 in (1). The following simplifying assumption is introduced to account for the pressure dynamics.

Assumption 1
The pressure P in each internal chamber is uniform, the gas is ideal, the conditions can be approximated as isothermal, and all orifices experience choked flow. In addition, the energy associated to the gas flow and to heat transfers is assumed to be negligible compared to the mechanical work of the gas.
Assuming uniform pressure is a reasonable approximation for small manipulators such as FMAs, and isothermal conditions represent a good approximation for pneumatic actuators [19] and PAMs [7]. Choked flow occurs if the upstream pressure P 0 is sufficiently large (i.e. P < cP 0 and P atm < cP, with critical pressure ratio for air c = 0.528, and atmospheric pressure P atm ) [30], which is realistic for this system (see Sect. 5).
The mechanical work of an ideal gas in isothermal conditions is = − where V is the final volume of the gas, V 0 is the initial volume, P is the absolute pressure (i.e. > 0 for compression). For inextensible soft manipulators, the chamber's volume can be approximated as V = V 0 + Ar n i=1 |q i |, where A is the cross-section area of the internal chamber and r is the distance from the centroid of the chamber to the centre of the section, and it can be expressed as If one chamber is pressurized such that P = P 1 while P 2 = P 3 = P atm , the mechanical work of the ideal gas relative to P atm can be approximated according to Assumption 1 as [7,19] = (P − P atm )(V 0 + Ar θ) log (V 0 /(V 0 + Ar θ)) . (2) Note here that the pneumatic circuit in Fig. 1c can only generate pressures above P atm . Computing the time derivative of the ideal gas law in isothermal conditions at temperature T yields the pressure dynamicṡ where R s is the specific gas constant, V 0 also accounts for the volume of the pipes, the flow rate u corresponds to the control input, and δ is the outflow from the exhaust valve, which is considered unknown. Including (2) in the mechanical energy of the system yields H = 1 2 p T M −1 p + + , where the system states are the position q ∈ R 2n+1 of the virtual joints, the momenta p = M(q)q, and the pressure P. The complete system dynamics can be expressed as ⎡ External forces and model uncertainties are included in the term δ ∈ R 2n+1 which is unknown and possibly time-varying, while D indicates the physical damping. Throughout the paper we indicate with I the identity matrix of appropriate dimensions and with [1 n ] and [0 n ] the column vectors with all elements equal to 1 or 0. Finally, for a scalar function f , ∇ x f represents the vector of partial derivatives in x. System (4) is further qualified by the following assumptions.

Assumption 2
The model parameters in (4) are exactly known. In particular, the structural stiffness k is constant and uniform along the length, while the physical damping D > D 0 + D 1 |q| 2 > 0 is constant and uniform along the length, with the diagonal matrices D 0 > 0 and D 1 > 0. This assumption is motivated by the constant section of FMAs, and by the nonlinear damping of pneumatic systems [9]. External forces and model uncertainties, including the weight of the manipulator, are accounted for in δ.

Assumption 3
The angles θ and γ and their time derivatives are known at any instant and are bounded. In addition, q i q j > 0, ∀i, j, thus all sections of the soft manipulator bend on a plane in the same direction. The internal pressure P is known at any instant.
The condition q i q j > 0, ∀i, j is less stringent than CC, which requires q i = q j , ∀i, j. In particular, θ and γ are measured with a tracking system and P is measured with a pressure sensor in the experiments (see Sect. 5). where a payload is typically applied), that is δ 0 i = 0 for i > n, while σ is the time-varying bounded part, with |σ | ≤ σ 0 |q| for a known σ 0 > 0. The disturbance δ is constant and unknown: this condition corresponds to choked flow through the exhaust valve (see Assumption 1).
ComputingṖ from (4) recovers the pressure dynamics (3). In addition, computingṗ from (4) yieldṡ where G T = 1 T n 0 T n+1 . Equation (5) indicates that the contribution of the internal pressure is modulated by θ . For conciseness, the notation h(θ ) = (1+log(1+ Ar θ V 0 )) is employed in the rest of the paper. Note that, in case V 0 Ar θ , (5) reduces to the simpler case of digitally controlled pressure discussed in [12].

Controller design for planar soft continuum pneumatic manipulators
Employing the IDA-PBC methodology, the regulation problem (q,q) = (q * , 0) for a mechanical system with energy H 0 = 1 2 p T M −1 p+ is addressed by designing a control law u 0 such that the closed-loop dynamics is representative of a new mechanical system with energy The main design parameters within IDA-PBC are the inertia matrix M d and the potential energy d , which should have a strict minimizer at the prescribed equilibrium point q * to achieve the regulation goal (see [28] for further details). Thus, the control law u 0 conveniently reshapes the energy H 0 through the parameters M d and d . Since H 0 does not include the mechanical work of the gas, a new controller design procedure is outlined here for system (4). This section focuses on the regulation goal on the bending plane θ = θ * by introducing the following additional assumption, thus it is also applicable to soft manipulators with a single internal chamber. Instead, the regulation problem in 3D (θ, γ ) = (θ * , γ * ) is addressed in Sect. 4.

Assumption 5
One chamber of the soft manipulator is pressurized (i.e.P = P 1 ), while the two remaining chambers are left at atmospheric pressure (i.e. P 2 = P 3 = P atm ).

Controller design procedure
The control law is designed such that the closed-loop dynamics in port-Hamiltonian form becomes ⎡ where H d = H d + ς 2 /2 is a positive definite and radially unbounded storage function, and ς is defined as with k p a tuning parameter, and δ 0 i the disturbance estimates to be defined in Sect. 3.2. The time-varying disturbance σ is not affected by the control action thus it appears again in (6) and its effects on stability are accounted for in Sect. 3.3. The terms S ij are computed according to the following multi-step design procedure such that the open-loop dynamics (4) matches the closed-loop dynamics (6).
Step 1: This first step aims to preserve the relationship between the position q and the momenta p. Equating the first row of (4) and of (6) gives Defining S 12 = k m I and M d = k m M, and setting S 13 = 0 verifies (8) leading to the following step.
Step 2: This step aims to reshape the kinetic energy and the potential energy of the system in closed-loop by defining appropriate expressions of the potential energy d (q) and of the inertia matrix M d (q) = M T d > 0. Equating the second row of (4) and of (6) yields Substituting S 12 = k m I and M d = k m M the kinetic energy vanishes from (9), while defining S 22 = k m D the damping terms cancel out. To ensure that the remaining terms satisfy (9) for all q ∈ R 2n+1 , we introduce the full-rank left annihilator G ⊥ , such that G ⊥ G = 0 and rank G ⊥ = 2n. Pre-multiplying both sides of (9) by G ⊥ , the terms dependent on ς or on P vanish, yielding the following partial differential equations (PDEs), which are termed potential-energy PDE and disturbance matching equation where 0 = T (q − q * ) and can be interpreted as a vector of closed-loop non-conservative forces [16].
Solving (10) and (11) while ensuring that the minimizer where Note that the first term in d preserves the potential elastic energy of the system outside the bending plane. Finally, substituting (12) and (13) in (9) and multiplying it by the pseudo-inverse G † = G T G −1 G T yields The term does not appear in S 23 since it follows from (13) that G T = 0. Note that (12) and (13) correspond to the case of digitally controlled pressure [12]. The pressure dynamics (3) is accounted for in the next step.
Step 3: equating the third row of (4) and of (6) yields where S 33 is a further tuning parameter. Computing the control input u from (15), substituting (14), and defining S 33 = k i Ar h(θ) > 0 yields the explicit expression of the control law The disturbance estimates δ 0 and δ are defined in the following section and are combined with the control law according to [16]. This approach is possible in general due to the separation principle [11,25].

Adaptive observers
The disturbance estimate δ 0 in (16) is computed with a modification of the I&I method [1] resulting in the adaptive laẇ similarly to the case of digitally controlled pressure [12], where α is a positive constant tuning parameter.

Lemma 1
Consider system (4) with Assumptions 1-5 and with the adaptive law (17). Define the vector of estimation errors z as where β 0 ( p) = −αp and δ 0 is computed by integrating (17) in time. Then z is ultimately bounded for all α > 1/4.
The disturbance estimate δ is computed in a similar fashion with the adaptive law where α is a further positive constant tuning parameter.

Lemma 2
Consider system (4) with Assumptions 1-5 and with the adaptive law (22). Then, the disturbance estimate δ converges exponentially to the correct value for all α > 0.
Proof We define the vector of estimation errors z as follows, where the functions δ and β (P, θ) are the state-independent part and the state-dependent part of δ Computing the time derivative of (23) and substitutinġ P from the third row of (4) giveṡ Substituting (22)  Thus z is bounded and converges to zero exponentially for all α > 0 concluding the proof The result in Lemma 1 is weaker than that in Lemma 2 since it only concludes boundedness of the estimation error z. This occurs because the adaptive law (17) does not include the contribution of the kinetic energy which depends on the virtual positions q i that are not measurable, and since it does not compensate the time-varying disturbances σ . However, the inequality α > 1/4 represents a sufficient condition because the Young's inequality used in Lemma 1 is conservative. Conversely, the adaptive law (22) does not introduce any approximation but ensures exponential convergence of the estimation error z to zero provided that the disturbance δ is constant (see Assumptions 1 and 4). In this case, the inequality α > 0 represents a necessary and sufficient condition. Further conditions for asymptotic convergence of z and z to zero are provided in Proposition 1.

Stability analysis
In this section, the effects of the adaptive laws (17) and (22) and of the control law (16) on the system dynamics are investigated and the stability conditions are discussed. Proposition 1 Consider system (4) with Assumptions 1-5 in closed-loop with the control law (16) where δ 0 is computed by integrating its time derivative (17) and δ is computed with (22). Define the parameters k i , k m , α, α such that the symmetric matrix Then the equilibrium point θ,θ = (θ * , 0) is stable and θ converges to θ * asymptotically. Additionally, Proof Substituting (16) in (4) results in the closed-loop dynamics in port-Hamiltonian form ⎡ where λ = (P−P atm )A 2 r 2 Ar θ+V 0 − k n + k p and the square matrix in (26) is skew symmetric.
Defining the Lyapunov function = H d +ϒ +ϒ , which accounts for the estimation errors z and z , and computing the time derivative˙ along the trajectories of the closed-loop system (26) yieldṡ Since the elements of M depend on the cosine of the virtual positions q i [12], we have that max {M} < m T l 2 T , with m T the mass of the soft continuum manipulator, l T its total length, and max {M} the largest element of M. Thus the inequality | p| ≤ c 1 m T l 2 T |q| and the inequality ∇ q p T M −1 p ≤ c 2 m T l 2 T |q| 2 hold for some 0 < c 1 , c 2 < 1. Recalling that by hypothesis |σ | ≤ σ 0 |q| yieldṡ Refactoring common terms in (28) yieldṡ where x T = q T (q 2 ) T ς z z and is given by (25).
Consequently˙ ≤ 0, thus x is bounded and converges to zero asymptotically and the equilibrium point x = 0 is stable. Computing (26) at the point x = 0 yields ∇ q d = 0 which implies that H d has an extremum at x = 0. Computing G T ∇ q d from (12) while recalling that G T = 0 from (13), and computing the minimizer condition ∇ 2 q d yields which hold true for all k p , k m > 0. Thus θ converges to θ * asymptotically concluding the proof Remark 1 Considering that m T l 2 T 1 and that V 0 Ar θ for small soft continuum manipulators similar to FMAs, the symmetric matrix (25) can be simplified by neglecting such terms, which yields the approximated expression The necessary conditions to ensure > 0 are thus The first inequality in (32) requires the physical damping to be sufficiently large to dominate the effect of time-varying disturbances. Solving the second inequality yields In the theoretical case of constant disturbances, this inequality can be further simplified as 4D 0 αk m > 1 which suggest that the parameters α and k m contribute alongside the physical damping parameter D 0 to the stability of the prescribed equilibrium. Finally, the third inequality in (32) indicates that a larger value of k i α is required for larger soft manipulators (i.e. with larger cross-section area of the internal chambers A). The effect of the tuning parameters on the system response is discussed in detail in Sect. 5.

Backstepping design
For comparison purposes, this section details a new alternative controller that extends our design [12] by employing the backstepping methodology [22]. This choice is motivated by the fact that backstepping controllers are well suited for variable stiffness robots [29]. While it would be possible to employ different adaptive observers such as [13,15] or [5], the same adaptive laws are employed as in (16) to highlight the effect of the proposed controller design procedure. In this case, the controller design consists of two main steps.
Step 1: Considering initially a reduced-order system by taking the absolute pressure P * as virtual control input yields the partial system dynamics where H 0 = H − = 1 2 p T M −1 p + is the mechanical energy of the system without the mechanical work of the gas. Computing the control law as in [12] yields with δ 0 computed form (17). Substituting (34) into (33) yields the closed-loop dynamics where H d = 1 2 p T M −1 d p + d with M d = k m M, and where d is given by (12) and is given by (13). Defining a new Lyapunov function candidate 1 = H d + ϒ and computing its time derivative along the trajectories of the closed-loop system (35) in a similar fashion to (27) yields in this casė Employing the inequality | p| ≤ c 1 m T l 2 T |q| and the inequality ∇ q p T M −1 p ≤ c 2 m T l 2 T |q| 2 for some positive c 1 and c 2 , and recalling that by hypothesis |σ | ≤ σ 0 |q|, yields˙ 1 ≤ −x T 1 1 x 1 with x T 1 = q T (q 2 ) T z and the symmetric matrix Thus, the system (33) in closed-loop with the virtual control law (34) has an asymptotically stable equilibrium point in θ,θ = (θ * , 0) provided that the symmetric matrix (37) is positive definite.
Step 2: Introducing the error ζ = P − P * and substituting P = P * + ζ in (33) yields the closed-loop dynamics Defining a new Lyapunov function candidate as 2 = 1 + 1/2ζ 2 and computing its time derivative while substituting (38), andṖ from (3) yieldṡ Defining the control input u in order to cancel the cross product where P * is given in (34), δ 0 is computed with (17), δ is given in (22), and k i > 0 is a tuning parameter. (4) with Assumptions 1-5 in closed-loop with the control law (40) and with the adaptive laws (17) and (22). Define the tuning parameters k i , k m , α, α such that the symmetric matrix
Proof Substituting (40) into (39) yields˙ 2 ≤ −x T 1 1 x 1 − k i ζ 2 . Defining the new Lyapunov function candidate 2 = 2 + ϒ and computing its time derivative while substituting (22) yieldṡ Substituting (37) into (42) and refactoring common terms yieldṡ (41) is positive definite by hypothesis, it follows that˙ 2 ≤ 0 thus x 2 is bounded and converges to zero asymptotically, and the equilibrium point x 2 = 0 is stable. Computing (38) at the point x 2 = 0 results in ∇ q d = 0. Computing G T ∇ q d and the minimizer condition det ∇ 2 q d yields again (30), which holds true for all k p , k m > 0. Thus θ converges to θ * asymptotically concluding the proof Remark 2 The virtual control law (34) has two parameters, namely the proportional gain k p and the parameter α in (17). The complete control law (40) has three additional parameters, that is the parameter α in (22), the parameter k m that serves as a derivative gain, and the parameter k i that accounts for the pressure dynamics. The control law (16) employs the same five parameters as (40). In particular, k p affects the closed-loop potential energy d . Large values of k p result in faster response but also increase the closedloop stiffness [8]. In this respect, values in the range k/10 < k p < k/3 represent a good compromise. The parameter k m scales the kinetic energy in closed loop [12] and should be chosen according to the stability conditions (see Remark 1). In particular, a smaller k m results in a slower response. The parameters α and α define the convergence rate of the estimation errors z and z . Finally, the parameter k i defines the convergence rate of the pressure dynamics, thus large values are permitted. For comparison purposes, a PID algorithm intended for the same regulation goal θ,θ = (θ * , 0) contains only the parameters K p , The PID represents a valid term of comparison since it is widely used in industrial practice, and it serves as basis of recent control algorithms specifically designed for soft continuum manipulators [40]. While (44) has a simple structure, it does not provide a physical interpretation of the control action and it does not account explicitly for nonlinearities or external disturbances. In practice, (44) might require higher gains or might result in less consistent performance across different operating conditions (see Section 5). Improved performance could be achieved for the PID with gain scheduling techniques. This however would require additional assumptions on the disturbances and is beyond the scope of the present work.

Remark 3
To better highlight the difference between (16) and (40) we proceed as follows: replace k i with k i = k i Ar h(θ ) in (40); assume that V 0 Ar θ ; subtract (40) from (16), which yields Equation (45) depends onθ and contains the adaptive estimate δ 0 and its time derivative˙ δ 0 , which depend on the parameter α. Thus, differences in performance between the control laws (16) and (40) can be expected in dynamic conditions and for larger values of α. A further difference between the control laws (16) and (40) becomes apparent considering the closed-loop pressure dynamics. Substituting (16) into (3) yieldṡ Instead, substituting (40) into (3) yieldṡ Compared to (46), the pressure dynamics (47) has additional terms that depend on the angular velocitẏ θ and on the disturbance estimate δ 0 . Consequently, the closed-loop dynamics with controller (40) cannot be described by using a skew-symmetric matrix as in (26). This also results in additional cross products in the Lyapunov derivative (42), thus (41) has different off-diagonal terms compared to (25).

Controller design for soft continuum pneumatic manipulators in 3D
In this section Assumption 5 is modified and the control problem for a soft continuum pneumatic manipulator operating in 3D is discussed. Assuming that two internal chambers are pressurised independently, while the third chamber is at atmospheric pressure for simplicity (i.e. P 3 = P atm ), the system dynamics can be expressed as ⎡ where we define P 1 = P 1 log 1+ 2 , while δ 1 and δ 2 are the outflows from the exhaust valves corresponding to each chamber. The mechanical energy is defined as H = 1 2 p T M −1 p + + 1 + 2 , where 1 and 2 are given by (2) for each chamber. In this case, the orientation of the tip is defined by (θ, γ ), where γ = ϕ only in the absence of out-of-plane disturbances. The control inputs u 1 and u 2 represent the flow rates to the respective internal chambers. Since the expressions of θ and ϕ given in (1) do not account for disturbances or uncertainties in the geometry of the internal chambers (e.g. due to manufacturing tolerances), it is necessary to estimate the mapping between the internal pressures P 1 , P 2 and the angle γ , as discussed in the following lemma. To this end the following assumption is introduced.
Approximating the orientation of the bending plane with a static mapping is appropriate for soft continuum manipulators similar to FMA [17] which bend on a plane according to (1) provided that the external forces (e.g payload) act predominantly on the bending plane (see Assumption 4). γ = f (η), define the observer η with the update laẇ

Controller design procedure
The controller design follows similar steps to those outlined in Sect. 3.1. Considering for conciseness an axial symmetric manipulator with equal cross-section areas of the internal chambers A 1 = A 2 = A, the term ς is defined as The closed-loop dynamics in port-Hamiltonian form is where H d = H d + ς 2 /2, ς is given in (52), and H d = Step 1: Equating the first row of (48) and of (53) gives Defining S 12 = k m I , M d = k m M, and S 13 = S 14 = 0 verifies (54) leading to the following step.
Step 2: Equating the second row of (48) and of (53) while recalling the definition of G yields Substituting S 12 = k m I and M d = k m M the kinetic energy vanishes from (55), while defining S 22 = k m D the damping terms cancel out. Pre-multiplying both sides of (55) by the left annihilator G ⊥ , the terms dependent on ς and on P 1 , P 2 vanish, yielding the potential-energy PDE (10) and the disturbance matching equation (11). Solving (10) while ensuring that the potential energy d has a strict minimizer in θ = θ * yields (12), while solving (11) for an unknown constant disturbance δ 0 yields (13). Finally, substituting (12) and (13) in (55) and multiplying it by G † yields where is a free term to be defined in the next step.
Step 3: Equating the third row of (48) and of (53) yields Similarly, equating the fourth row of (48) and of (53) yields . In order to define , we compute the time derivative of η = P 2 /(P 1 + P 2 ) substituting (49), which yieldṡ Substituting (59) into (60) yields Computing from (61) yields finally θ Computing u 1 and u 2 from (57) and (58) while substituting (62) and assuming k i = k i for conciseness yields the complete control law where δ 1 and δ 2 are given by (22), ς is given by (52) and contains the adaptive estimate δ 0 computed with (17). The general case A 1 = A 2 is conceptually similar and is omitted for brevity.

Stability analysis
The stability of the closed-loop system (53) is discussed in this section using a similar procedure to Proposition 1.
Proof Define the Lyapunov function candidate = H d + ϒ +ϒ 1 +ϒ 2 , where ϒ 1 and ϒ 2 refer to the estimates of the outflows δ 1 and δ 2 in the pressurized chambers and correspond to the estimation errors z 1 , z 2 . Computing the time derivative of along the trajectories of the closed-loop system (53) yieldṡ Substituting in (65) the inequalities | p| ≤ c 1 m T l 2 T |q| and ∇ q p T M −1 p ≤ c 2 m T l 2 T |q| 2 with the constants 0 < c 1 , c 2 < 1, recalling that by hypothesis |σ | ≤ σ 0 |q|, and refactoring common terms yieldṡ which can be rewritten in compact form aṡ where x T = q T (q 2 ) T ς z z 1 z 2 and is given by (64). Consequently˙ ≤ 0, thus x is bounded and converges to zero asymptotically, and the equilibrium point x = 0 is stable. Computing (53) for x = 0 results in ∇ q d = 0 and computing the minimizer conditions yields again (30) which hold true for all k p , k m > 0. Finally, it follows from Lemma 3 that η = P 2 /(P 1 +P 2 ) converges to η * for all k I > 1 , where γ * = f (η * ). Thus (θ, γ ) converge to (θ * , γ * ) asymptotically concluding the proof

Simulations
Simulations have been conducted in MATLAB for the planar system (4) and for the 3D system (48) with n = 3 employing the model parameters k = 2, A = 10, r = 0.003, and D 0 = 0.03, D 1 = 0.01 for illustrative purposes. The model consists of n = 3 virtual links of equal length with mass concentrated at their midpoint for simplicity, such that the total mass and the total length of the rigid links are m T = 1.5 and l T = 0.1. These values refer to an ideal soft continuum pneumatic manipulator and have been chosen to better highlight the differences between the control algorithms. The links at rest are collinear (i.e. q = 0), the effect of gravity is neglected, and an unmodelled force f 0 = 2 parallel to the neutral axis at rest acts on the tip of the manipulator resulting in a bending moment against the direction of positive θ . The remaining parameters are the dead volume V 0 = 10 Al T , the gas constant for air R s = 287, the temperature in Kelvin T = 293, and the supply pressure in bar P 0 = 6. The tuning parameters for controller (16) have been set as k p = 0.5; k m = 2.5; k i = 30; α = 4; α = 2.5. These values ensure > 0 in (31) for all σ 0 ≤ D 0 /15. The parameters of controller (40) have been set as k p = 0.5; k m = 2.5; k i = 5; α = 4; α = 2.5 to achieve a similar response to that of controller (16). Similarly, the parameters of the PID algorithm (44) have been chosen empirically as K p = 50; K d = 10; K i = 10 to obtain a comparable response. The simulation results in Fig. 2 show the time histories of the tip rotation θ for system (4) highlighting the effect of the physical damping D 0 , which is instrumental to stability (see Remark 1). All controllers correctly achieve the regulation goal. In particular, controller (16) results in a smooth transient for both values of D 0 . Controller (40) results in comparable performance to (16) when D 0 = 0.03, but it leads to vibrations when D 0 = 0.01. The PID (44) results in a larger control input during the transient, and it also leads to vibrations in case of lower D 0 . In summary, the simulation results indicate that the control law (16) is less sensitive to physical damping. It must be noted that the improved performance of controller (16) is due to the proposed controller design procedure, which applies the IDA-PBC methodology in a multi-step fashion. The adaptive observers (17) and (22), although important to achieve the prescribed regulation goal, are not responsible for the difference in performance since they are common to the backstepping design (40). This observation complements the findings of our previous works [12,13,15], which have shown that similar performance can be achieved for this class of systems by combining different adaptive observers with the same energy-based control law. The effect of the tuning parameters in the control law (16) is shown in Fig. 3. A larger k m , a larger k p , and a larger k i result in a faster transient at the cost of a higher control action but they do not trigger vibrations. In particular, a larger k p results in a steeper gradient of the potential energy d , while a larger k m results in a larger closed-loop inertia. Instead, a larger k i increases the convergence rate for the term ς (see Remark 2). The parameters α and α have a less noticeable effect on performance, provided they meet the stability conditions outlined in Proposition 1 and in Remark 1. This result is in line with our previous work [12,13] and it highlights the benefits of the I&I methodology [1], which allows designing the observer independently from the control law. Figure 4 shows the time histories of the in-plane rotation θ and of the out-of-plane rotation γ for system (48). Also in this case, the unmodelled tip force f 0 = 2 generates a moment on the bending plane against the direction of positive θ . The tuning parameters in (63) have been set similarly to (16), that is k p = 0.5; k m = 2.5; k i = k i = 20; k I = 20; α = 4; α = 2.5. In this case the PID algorithm (44) has been employed to compute u 1 and u 2 , where the former depends on θ and the latter on γ , with the parameters K p = 50, K d = 10, and K i = 10. Both controllers achieve the regulation goal in a smooth fashion with the chosen parameters: the control law (63) results in a faster convergence of θ and γ ; conversely, the PID leads to slower response and also requires a larger control input during the transient. Increasing the proportional gain in the PID (e.g. K p = 100) leads to faster convergence but also to larger vibrations in θ , which are undesirable. Since the simulations indicate that controller (16) is preferable to controller (40), the latter has not been extended to the regulation in 3D. Figure 5 shows the comparison between controller (63), the PID (44), and the SMC algorithm [2] reported below in two different operating conditions corresponding to f 0 = 2 and f 0 = 10. The SMC algorithm has been employed to compute u 1 and u 2 , where the former depends on θ and the latter on γ , while the function tanh(·) has been used to reduce chattering. The tuning parameters have been chosen empirically as K e = 5, K i = 30, μ = 10, β 0 = 5 to achieve a similar response to controller (63) with the tip force f 0 = 2.
The results indicate that the controller (63) produces a consistent response with the tip force f 0 = 10. Instead, the PID and SMC algorithms result in a slower convergence of θ and in a noticeable overshoot on γ with f 0 = 10. In addition, both PID and SMC produce a higher and more oscillatory control input compared to controller (63).

Experiments
The proposed control laws have been tested on a soft continuum pneumatic manipulator prototype that measures 6 mm in diameter, 30 mm in length, and that has an approximate mass m T = 1.5 grams. The tip rotations θ and γ have been measured with an electromagnetic (EM) tracking system (Aurora, NDI, Canada) and an EM sensor (Aurora 5 DOFs sensor with 0.5-mm diameter and 8-mm length, part number 610061, Northern Digital Inc) that has a root mean square (RMS) accuracy of 0.2 • . The flow rate to the internal chambers is provided by a needle valve (part number 7770 06 00, Legris) operated by a servo motor (Servomotor RC 6V, Parallax Inc) as shown in Fig. 1c. The pressure relative to atmosphere is measured with a pressure sensor (MPX2200GP, NXP Semiconductors), as shown in Fig. 1b. A MATLAB script records the pressure measurement, the tip rotations θ and γ , and sends the control signal to the servo motor using an embedded microcontroller (mbed NXP LPC1768, NXP Semiconductors) via serial link (baud rate 921600). A manual pressure regulator is employed to set the supply pressure to a constant value P 0 = 4 bar. It must be noted that, differently from the simulations, the control input u in the experiment is a normalized signal between 0 and 1 that corresponds to the position of the servo motor between 0 and 180 • . In isothermal and choked flow conditions (see Assumption 1) the mass flow rate from an orifice can be expressed as Q = C P 0 ρ 0 , where C is the conductance of the orifice and ρ 0 the density of the gas [30]. For simplicity it is assumed that the conductance varies with the control input u in a linear fashion, that is C = C 0 u where C 0 = 1. To avoid damaging the prototype, the exhaust valves have been set such that the absolute pressure in the internal chambers is smaller than 3 bar when the position of the servo motor is 180 • . According to Assumption 5, two chambers are left at atmospheric pressure for the experiments in 2D, that is P 2 = P 3 = P atm , while only one chamber is at atmospheric pressure for the experiments in 3D, that is P 3 = P atm . In the latter case, two needle valves with servo motors and two correspond-ing exhaust valves have been employed to pressurize two chambers of the prototype. The control laws have been implemented with the model parameters m T = 1.5, l T = 0.03, k = 1, D 0 = 0.03, D 1 = 0.015, which have been estimated experimentally [14]. The remaining model parameters are the same as in the simulations. The tuning parameters in the control laws (16) and (40) have been set as in the simulations, that is k p = 0.5; k m = 2.5; α = 4; α = 2.5 with k i = 30 for (16) and k i = 5 for (40). The parameters of the PID (44) have been set empirically as K p = 50, K d = 5, K i = 70 to obtain a comparable response. For the experiments in 3D, the tuning parameters of controller (63) have been set in a similar way to the simulations, that is k p = 0.5; k m = 2.5; k i = k i = 30; k I = 30; α = 4; α = 2.5. In this case the PID algorithm (44) has been employed to compute both u 1 and u 2 , where the former depends on θ and the latter on γ , using the parameters K p = 100; K d = 5; K i = 100. These values have been chosen empirically in an attempt to achieve a comparable response to that of controller (63). The SMC algorithm [2] has not been employed in the experiments since the simulations indicate that it yields a similar performance to the PID.
The system response on the bending plane with the control laws (16) and (40) and with the PID (44) is shown in Fig. 6. With the chosen tuning parameters, the response is similar for all controllers. Nevertheless, the PID results in small vibrations in θ that also appear in the control input and that affect the pressure dynamics. Controllers (16) and (40) lead to a similar response due to the large damping of the system, as anticipated by the simulation results. In the presence of an unmodelled external load due to a mass m 0 = 3.5 g attached at the tip of the manipulator, the angle θ requires a longer time to reach the prescribed value with all controllers. However, the difference is less noticeable with the proposed controller (16). Further experimental results that refer to a different setpoint θ * are shown in Fig. 7. In this case, the transient response for controllers (16) and (40) is similar. Instead, the response with the PID is noticeably different from Fig. 6 and it is characterized by a more oscillatory nature during the transient and in proximity of the equilibrium. In summary, the results suggest that the control law (16) leads to a more consistent response for different setpoints θ * and in the presence of disturbances, while the PID might have to be re-tuned depending on the operating condition.
The experimental results for the prototype in 3D are shown in Fig. 8. The controller (63) and the corresponding PID correctly achieve the regulation goal in the presence of the external load due to the tip mass m 0 = 3.5 grams, which affects the system dynamics predominantly on the bending plane. With the tuning employed, the PID leads to noticeable vibrations in θ , both during the initial transient and around the equilibrium. Instead, the controller (63) results in a smooth time history of the tip angles and of the pressures in the internal chambers. In particular, the transient response of θ is similar to that in Fig. 7. Differently from the simulations, a small overshoot occurs on γ and it is ascribed to the differences between the two needle valves supplying the two chambers (e.g. different friction forces and different conductance of the corresponding orifices). Nevertheless, this does not trigger unwanted vibrations.
For completeness, Fig. 9 shows that the proposed controller (16) yields similar performance to our previous implementation which relies on digital pressure regulators [12] where δ 0 is computed from (17). In both cases the tip rotation θ reaches the prescribed value without overshoot: the settling time is approximately 6.5 seconds for controller (16), while it is approximately 5 seconds for controller [12] with the tuning parameters employed (i.e. k p = 0.3; k m = 20; k v = 1; α = 10). This difference is due to the faster response of the digital pressure regulator (≈ 10 ms) compared to the servo motor with the needle valve. It must be noted that the control input represents the position of the servo motor for controller (16), while in [12] it corresponds to the output pressure of the digital pressure regulator, thus it is not directly comparable.
Compared to the simulations, the experimental results show a slower response with all controllers. This is due to the following factors: the model uncertainties and the disturbances affecting the system in the experiments also include the weight of the prototype, the structural stiffness, and the unknown conductance of the needle valve; a lower supply pressure P 0 has been employed to avoid damaging the prototype; the response of the servo motor is not instantaneous. In  (16); d corresponding control input; g pressure P 1 − P atm ; b tip rotation θ with controller (40); e corresponding control input; h pressure P 1 − P atm ; c tip rotation θ with PID (44); f corresponding control input; i pressure P 1 − P atm . The tip mass is m 0 = 3.5 grams particular, employing a lower P 0 limits the pressure in the internal chambers thus reducing the responsiveness. In addition, while the effect of the control input is instantaneous in the simulations, the servo motor has a nominal maximum speed of 140 rpm at no-load, and the experimental measurements indicate that a 180 • rotation takes approximately 0.7 seconds when the motor shaft is attached to the needle valve. Finally, the valve conductance might not change linearly with the position of the servo motor, which has also a non-negligible dead-band. In this respect, employing a faster servo motor with a specially designed needle valve could fur-ther improve the performance with all controllers and shall be investigated in our future work.

Conclusion
This paper presented a new energy-based control strategy for a class of soft continuum pneumatic manipulators that can bend on any plane. The controller design procedure employs the IDA-PBC methodology in a multi-step fashion to account for the pressure dynamics of the pneumatic actuation. Nevertheless, we believe that the proposed approach has a more general value, and it could be applied to other types of actuation (e.g.  (16); d corresponding control input; g pressure P 1 − P atm ; b tip rotation θ with controller (40); e corresponding control input; h pressure P 1 − P atm ; c tip rotation θ with PID (44); f corresponding control input; i pressure P 1 − P atm . The tip mass is m 0 = 3.5 grams hydraulic). The simulation results indicate that the proposed controllers correctly achieve the regulation goal without triggering vibrations. In comparison, an alternative control law that has been constructed with a backstepping method can lead to unwanted vibrations for systems with smaller physical damping. A comparative simulation study with a PID and an SMC algorithm suggests that the proposed controller achieves the regulation goal with a smoother control action, and that the transient performance is less affected by external disturbances. The validity of the proposed approach has been confirmed with experiments on a prototype in the presence of unmodelled external forces acting on the bending plane and due to a tip mass. The comparison with a PID algorithm indicates that the proposed controllers lead to a more consistent performance across different operating conditions, without changing the tuning parameters. The experimental results also confirm that employing a needle valve operated by a small servo motor allows achieving the regulation goal, thus it could represent an alternative to digital pressure regulators. This approach could help reducing the cost of soft robotic systems thus promoting their adoption in low-income countries.  Fig. 9 Experimental results showing the regulation on the bending plane without tip mass for controller (16) indicated as 'Needle valve' and for our previous implementation [12] indicated as 'Pressure regulator': a time history of the tip rotation θ; b control input some of the simplifying assumptions introduced in this paper, such as those of isothermal conditions and choked flow, to integrate proprioceptive sensing in our prototypes, and to develop a bespoke needle valve. We also plan to extend the control formulation to tracking tasks and to apply it to more complex soft robotic systems consisting of multiple manipulators connected in series and in parallel as well as to soft manipulators that can produce twisting moments. Finally, we aim to compare the proposed controllers with a wider range of competitor solutions as part of a larger experimental study that considers different types of disturbances.