There is a significant body of work relative to accurately modeling the kinematics and dynamics of soft robots. In Renda et al. (2012) and Thuruthel et al. (2016) the continuum joint is modeled using Cosserat-beam theory. In Kang et al. (2011) and Khalil et al. (2007) methods based on recursive Newton-Euler approaches are used, while in Tatlicioglu et al. (2007) and Godage et al. (2011) dynamic equations are derived using Lagrangian mechanics. In Zheng et al. (2012) and Giri and Walker (2011) lumped parameter models are derived by dividing the continuum joint into a number of finite length sections. The trade-off between accuracy and computational complexity in these methods can be seen by varying the number of the finite sections. The authors of Walker (2013) provide a more comprehensive review of dynamic models for soft and continuum joint robots. Notably, there has also been work to show that learned models can represent soft robot dynamics as in Thuruthel et al. (2017).

In Mochiyama and Suzuki (2002) and Mochiyama and Suzuki (2003) the authors derive the dynamic equations of a continuum arm by integrating over infinitesimal disks and using the method of Lagrange. No assumptions of constant curvature are made. These works are similar to the modeling efforts presented in this paper, the main differences being our choice of generalized coordinates and our assumption of constant curvature. These two differences allow us to derive closed-form analytical expressions for the terms in our equations of motion such as the mass and Coriolis matrices.

In Falkenhahn et al. (2014) and Falkenhahn et al. (2015) the authors derive simpler models based on the PCC assumption. However, they neglect generalized forces caused by rotational inertias. They also model the mass of each PCC section as being concentrated at a point that is fixed in some coordinate frame. In Della Santina et al. (2020b), the authors derive a similar PCC-based model (also neglecting rotational inertia) and then match it to a dynamically equivalent rigid body model. Because the mass and inertia of the joints used in our work are non-negligible, we model the mass as distributed uniformly throughout infinitesimal disks and the center of mass of each joint is calculated analytically assuming uniform density. This approach yields closed-form equations of motion for the continuum joint while more accurately representing the dynamics by including the effects of rotational inertia. This approach also illustrates to a greater extent the effect of dynamic models that include rotational inertia on the performance of highly underdamped systems when compared to the work found in Della Santina et al. (2020b).

Control strategies for soft robots vary from open-loop control such as in Shepherd et al. (2011) and Tolley et al. (2014) to Reinforcement Learning (Zhang et al., 2017) to model predictive control (Best et al., 2016). In Hyatt et al. (2019) and Hyatt and Killpack (2020) the authors demonstrate the performance of MPC on the same joints used for this work. These implementations of MPC used a learned model of the dynamics based on a less-accurate representation of the continuum joint dynamics. The model inaccuracy that resulted in less aggressive control in that work prompted the development of the more accurate model and adaptive control techniques presented in this paper.

Given a dynamic model of the correct form, the nature of soft robots is still such that certain parameters of that model may be difficult to estimate. In terms of adaptive control for soft robots, the most similar to our work is Trumić et al. (2020), where they use a similar formulation of MRAC (although with a different dynamic model and no optimal control law). Although not common in soft robotics, combining MPC and adaptive control is beginning to be an established control strategy where the strengths of MPC are combined with a variety of adaptive control schemes (see Adetola et al., 2009; Kim, 2010; Chowdhary et al., 2013; Bujarbaruah et al., 2018; Pereida and Schoellig, 2018; Abdollahi and Chowdhary, 2019; Zhang and Shi, 2020). The method developed in this paper is a unique form of adaptive MPC that borrows ideas from model reference adaptive control (MRAC) for robot manipulators (Slotine and Li, 1987). Specifically, our work can be considered an extension to the adaptive MPC presented in Terry et al. (2019). The main extensions are an adaptive law formulated specifically for robot manipulators and a regressor based on a more accurate continuum joint dynamic model. These extensions allow greater flexibility to adapt both the parameters and structure of the model.

The robot used for this work is composed of a continuum joint such as the one seen in Figure 1. These joints are made of four separate pressure-controlled chambers surrounding a relatively inextensible central cable. The two antagonistically placed pairs of pressure chambers allow the joint to bend about two axes. We choose to model the kinematics of this joint using arcs of constant curvature. Each arc, which traces out the path in space occupied by the inextensible spine, can be defined using three variables as described in Allen et al. (2020) (see similar derivation in Della Santina et al., 2020a). These variables are the length of the in-extensible spine (h) and two components of the axis-angle vector that describes the rotation from the bottom of the joint to the top. Because the joint cannot twist about the inextensible spine (to which the z axis is tangent) the axis-angle vector consists of only two non-zero variables which we call u and v. These values are labeled in Figure 1 and correspond to the rotation about the x and y axes, respectively. We assume that the spine is perfectly inextensible so that h in this work becomes a constant kinematic parameter.

First we note some useful kinematic relationships. Because u and v are the non-zero elements of the axis-angle vector we can write

where ϕ is the magnitude of the axis-angle vector [u, v, 0]^T, or total bend angle (see Figure 2).

Although the joint is modeled as an arc with an arc length of h, we often want to refer to a position at some intermediate point along the arc. We denote an intermediate length along the arc using the variable l where l can take on any value between 0 and h (see Figure 3). Note that a frame tangent to the arc at a length l rotates as l is increased, therefore ϕ_l, u_l, and v_l are not constant along the entire arc. However, we note that the vector ρ→ from the base of the joint to the center of curvature is the same for all points along the arc because the center of curvature does not move. At any point l along the arc this vector can be calculated as

Because the magnitude of this vector ‖ρ→‖ is the radius of curvature, we may also relate ϕ and l using the arc-length formula

We now wish to derive a means by which we can calculate u_l and v_l at any point l along the arc given only l, h, and u and v at the end of the arc. Given a point that lies at a distance l along the arc, we may say using Equation (2)

Replacing ϕ terms using Equation (3) we obtain

Differentiating with respect to time yields the relationship

In other words, the generalized coordinates u_l and v_l and their time derivatives vary linearly along the length of the arc. This becomes a very useful property of this kinematic representation when deriving equations of motion.

Using the method of Lagrange, the equations of motion for a system of rigid bodies take the form

where M(q) is the mass matrix, C(q˙,q) is the Coriolis matrix, g(q) is a vector of gravity torques, q is a vector of the generalized coordinates, and τ is a vector of the generalized torques including friction terms. These matrices are derived using partial derivatives of kinetic and potential energy terms. Since partial derivatives are easily taken using a symbolic mathematics toolbox such as Sympy (see Meurer et al., 2017), the problem of dynamic modeling is reduced to the selection of generalized coordinates and the representation of kinetic and potential energy.

In order to accurately express kinetic and potential energy we choose to model the continuum joint, as many have done before, with an infinite set of infinitesimally small disks. However, the assumption of constant curvature, the choice of generalized coordinates, and current tools in symbolic math libraries allows us to produce analytical expressions for M, C, and g, whereas previous methods have not yielded these closed-form expressions.

We can define the kinetic energy of an infinitesimally thin disc at a length l along the arc as

where μ is the linear density of the disc, dl is some infinitesimal length, p._l is the velocity of the center of the disc, ω_l is the angular velocity of the disc expressed in the disc frame, and I is the inertia of the infinitesimally thin disc expressed in the disc frame.

The linear and angular velocity of each disc (p._l and ω_l) can be found using a configuration dependent jacobian J (meaning it is a function of joint configuration variables u_l and v_l) that is defined such that

A definition of this Jacobian for the choice of u and v as generalized coordinates can be found in Allen et al. (2020).

Using this relationship, we see that we can simplify the expression for kinetic energy (Equation 8) by scaling portions of the jacobian. The new inertia-weighted jacobian is defined as

allowing us to rewrite Equation (8) for the kinetic energy of a disc as

By treating a continuum joint as a series of infinitesimal disks and integrating the kinetic energy of each disc along the length of the arc we can write the total kinetic energy of a joint as

We note here that the Jacobian can be expressed analytically at every point along the joint as a function of l and the configuration variables u_l and v_l (which are lhuh and lhvh, respectively) thanks to Equation (5). Given this analytical expression for J_weighted we can integrate this expression with respect to l over the definite bounds 0 to h to get an analytical expression for JweightedTJweighted, which we recognize as the joint space inertia matrix or mass matrix M.

We use the symbolic mathematics library Sympy (see Meurer et al., 2017) to calculate JweightedTJweighted, and to integrate this expression analytically between the definite bounds 0 and h in order to obtain M(q). Once M(q) has been obtained symbolically, it is then relatively straightforward to take partial derivatives using Sympy in order to obtain an expression for the Coriolis matrix C(q˙,q) from Equation (7) using the method outlined in Bruno et al. (2010). The resulting coefficients that multiply q˙ to calculate the Coriolis matrix are commonly called the Christmases symbols of the first kind.

In order to find the gravity torques (g) we must first find the vector from the joint base to the joint center of mass (p→). By inspection we can see that a joint's center of mass must project down onto the vector ρ→ which is from the center of curvature to the base of the joint, however the vector to the center of mass must also contain some component in the z direction (orthogonal to the plane of the bottom plate of the joint). We find the components of the center of mass vector p→ by again dividing the joint into a series of infinitesimal disks of height dl.

Using the definition of the center of mass assuming the joint has uniform density along its length, the portion of p→ along the z axis is given by

Using the trigonometric relationship seen in Figure 3, namely

as well as the volume formula for an infinitesimally thin disc

we can now integrate to find z-:

Recognizing that ‖ρ→‖=hϕ,

In order to find the component of p→ that lies in the plane of u and v we follow a similar procedure. We will use x to represent the portion of p→ that lies along ρ→. Using the trigonometric relationship seen in Figure 3, namely

we can now integrate to find x-:

Recognizing that ‖ρ→‖=hϕ,

Using the derived equations for z-, x-, and the normalized version of ρ→ we obtain the vector from the base of the joint to the center of mass:

The potential energy of the joint due to gravity is simply the dot product of this vector, expressed in the inertial frame, with the gravity vector (G→) expressed in the same frame:

Having calculated the potential energy due to gravity, the gravity torques are calculated simply by taking the negative partial derivative of V with respect to q:

The method above has yielded analytical expressions for M, C, and g with the generalized coordinates u and v. Although complex, these closed-form expressions can be exported from the Sympy library into C code that can be evaluated within microseconds, allowing for real-time model-based control of these continuum joints.

In the absence of applied pressures, the joints used for this paper tend to drive themselves toward an equilibrium position at roughly u = v = 0 with slight overshoot and brief oscillation. This spring force could have been modeled as a part of the potential energy, however we choose to model the spring and damping separately from the traditional Lagrangian equations of motion. We approximate the spring forces as a linear spring term K_springq and friction as a linear viscous damping term Kdq˙. Including these terms, the final model used is,

In this section we give brief overviews of both MPC and MRAC in order to clarify notation and establish a background for the development of MRPAC. For in-depth explanations of MPC and MRAC we refer the interested reader to Hyatt et al. (2020) and Lavretsky and Wise (2013) respectively.

In the simulation portion of the experiments, a simulation is created using the model outlined in section 2.1 and this simulated system is controlled using three different controllers. The goal of each controller is to drive the system to follow a reference trajectory generated by a reference system. The three controllers implemented are MPC, MRAC, and the MRPAC algorithm detailed in section 2.2.3.

The reference system used for these experiments can be thought of as two uncoupled, critically-damped mass-spring-damper systems each modeled by the equation:

The masses (of mass m) are driven by the springs to the reference positions (r) and the damping coefficient (b) is always chosen such that the system is critically damped (b=4mk). The rise time of the reference system can be altered by varying the spring constant (k). We choose a rise time such that the system has settled to steady state within about 1 s.

As mentioned in the adaptive control literature, model parameter estimation and adaptive control schemes require sufficient “excitation” in order to converge or to adapt. We provide this excitation by changing the reference positions (r) of our system every 2 s. Reference positions are drawn from a uniform distribution bounded above and below by -π22 and π22. These bounds are chosen so that the resulting total bend angle (ϕ=u2+v2) is never greater than π2 radians.

In order to validate both simulations, we implement the same three controllers (MPC, MRAC, and MRPAC) on the soft continuum joint shown in Figure 1 and compare their performance.

The soft continuum joint used for this experiment is actuated by four plastic bellows, each of which can be controlled independently. A pressure difference in each of the bellows causes a rotation about one or both of the joint's axes. The angle about each of these axes (denoted u and v in Figure 1) is the robot's position and are the variables that we attempt to control. We expect this hardware platform to illustrate the sensitivity of each controller to both parameter mismatch and structure mismatch.

Both sources of error are present in hardware. Because no system identification was performed previously, the aforementioned model parameters such as h, m, r, K_spring, and K_d are not known perfectly. Additionally, the continuum joint exhibits unknown non-linear behavior near the extremes of its range of motion or in certain directions, where its stiffness or damping vary non-linearly with respect to u and v. In addition to the non-linear effects, we observe the effects of various offset forces in the plastic bellows used to actuate the joint. For example, even with equal pressures in each of the four bellows, the continuum joint remains slightly bent, indicating the presence of some constant unmodeled forces. For our simulations (see section 3.1.2) we represented this as a constant spring offset, but the actual source of this offset is unknown. In order to allow the adaptive control methods to compensate for this constant offset force we add to the regressor an identity matrix. This identity matrix means that the adaptive parameters that multiply it will be mapped directly to generalized torques in the dynamic model.

We track the orientation of a frame on top of the joint relative to a frame below the joint in order to estimate the state of the joint in real-time. We reuse the same reference trajectory from the simulation with one minor change: the command changes every 5 s instead of every two. This was adjusted in an attempt to be conservative with experimental hardware and software while still validating the performance of each controller.

As in the simulation experiments, we excite the system with the same 150 commands used in simulation (12.5 min) before evaluating each of the controllers for the last 30 commands (2.5 min). The joint trajectories for this evaluation period are shown in Figure 9.

The joint trajectories for the hardware experiments are shown in Figure 9 and the integral of the position tracking error is reported in Table 1. It is important to note that, unlike for the simulation, we cannot separate the perfect regressor and imperfect regressor cases on real hardware. Because of the nature of the continuum joint, we expect some combination of both cases to influence the controller performance results.

Generally, we see from the results that MPC struggles to eliminate steady state error. This matches the simulated behavior in Figure 5 and is expected because MPC does not have the ability to compensate for modeling errors that exist in the continuum joint. MRAC and MRPAC, on the other hand, do have the ability to compensate for modeling errors. Consequently they both track the steady state reference trajectory much closer than MPC. This indicates that the hypothesis presented in section 3.1.1 is demonstratively true at least for this hardware platform. MRAC and MRPAC certainly compensate for the modeling errors and drive the system to follow the reference trajectory. In hardware however, we see that neither controller is capable of following the reference trajectory exactly. In other words, we do not see in hardware the same performance as we see in the simulation results in Figure 5, where both trajectories deviate very little from the reference. This is because in addition to the modeling error (parameter mismatch) for which MRAC and MRPAC can compensate, there are still system dynamics for which they cannot fully compensate (structure mismatch).

The effect of structure mismatch in simulation is shown in Figure 6. Tracking error increases for all control methods as the magnitude of these modeling errors increase, but they increase dramatically for MRAC, hence its poor simulation performance exhibited in Figure 7. Importantly, this same pattern emerges in our hardware experiments. There are several instances during the evaluation period where unknown forces cause deviation from the reference trajectory. For examples of this, see the upper plot (u) of Figure 9 at 65, 100, and 135 s and the bottom plot (v) at 30, 45, and 95 s. All controllers are negatively affected, but MPC and MRPAC are more robust than MRAC. In other words, when encountering such disturbances, MRAC is forced to artificially adapt dynamic parameters in an attempt to eliminate the error. In contrast, MPC and MRPAC are better able to respond to disturbances because they re-solve the trajectory optimization over the whole time horizon, not just a single time step. These results support the hypothesis outlined in section 3.1.2 as well. MRAC and MRPAC do not track the reference trajectory perfectly because of the unknown disturbances but MPC and MRPAC are quantifiably more robust to the structure mismatch.

The results reported in Table 1 add a quantitative performance analysis in addition to the qualitative analysis from Figure 9. From the table we can see that MRPAC accumulates about half of the integrated tracking error of the other two controllers during the 1 min evaluation. It is interesting to note that MPC and MRAC have similar integrated tracking error, although qualitatively their trajectories look different. While MPC has a good transient response and large steady state error, MRAC has a poor transient response and small steady state error. This is also reflected in the statistics, since MRAC has lower mean and median error than MPC, but a higher standard deviation. According to these results, it seems that MRPAC has taken the strengths of the two approaches yielding a good transient response and smaller steady state error.