Our soft robot, “Squishy”, is an inflatable elastomeric chamber made of Smooth-On Dragon Skin 30 silicone that has a thin band of fabric embedded longitudinally to reinforce one side of the robot. The details of the manufacturing process and basic characteristics, such as the workspace volume and inflation-displacement behavior can be found in (Williamson et al., 2021). The undeformed shape of the chamber and the fabric band are 2D circular arcs. The chamber is closed by two caps at either end. One of the two caps has an inlet so that compressed air can be pumped in. When inflated pneumatically, the unreinforced side can undergo large strains while the fabric side maintains a constant length. This causes the chamber to bend and twist. Because of the fabric’s placement on the robot, its tip can trace 3D curves rather than just 2D curves as would a straight robot. Figure 1 shows the poses of the robot at some different working pressures.

The higher the robot chamber pressure is, the more it expands and the more it bends down and twists. The entire length of the robot arm could be used to perform some tasks such as scraping or pushing against objects in its environment. The advantage of this soft robot is that it can perform some tasks in 3D space with only one actuator (which is the controlled compressed air P). It should not damage objects, nor itself in the case of collisions due to the soft body. When the robot is not in use, it can be deflated, rolled, and stored in a tiny box to save space. The robot can also be included as a single module within a larger system to perform more sophisticated tasks.

To represent the motion of the robot, which is a continuum, using a finite set of variables, we use the Disc-thread model, which is briefly reviewed in this section. This was introduced in Schultz et al. (2020), and has the advantages of describing the bending and twisting of the soft robot in an easily customized fashion by a set of differential equations. In this approach, we discretize the chamber longitudinally into a sequence of N discs, each connected to its neighbors on either side by a single inextensible thread (representing the fabric reinforcement). Each disc is considered to be a rigid body and is constrained only by the thread. The frame assignments and their implications are illustrated in Figure 2. Each pose of the soft robot is derived from the relative position and orientation between pairs of adjacent discs. Each disc is parameterized by five generalized coordinates. The thread i is parameterized by two angles α _i , β _i ∈ [0; π]. The subsequent frames will involve a translation along the inextensible thread (of a fixed distance ℓ _i ), together with three other rotations γ _i , ψ _i and ϕ _i , where the twisting of the robot is accounted by γ _i . So by relying on these fictitious discs and threads as well as angles from α ₁ to ϕ _5(N−1), this approach models the soft robot as a kinematic chain of rigid bodies. However, because the joints are not prismatic or revolute this model contains a higher number of variables than a traditional kinematic chain.

With the kinematics of the soft robot fully defined by the disc-thread model, we formulate the Lagrangian and took its derivatives to find the equations of motion using the Euler-Lagrange formulation. The elasticity of the air is modeled by treating the air as an ideal gas, which exerts a normal force on the Nth disc in the z _N direction. The elasticity of the walls is modeled by connecting springs from some origin on disk i to insertion on disk i + 1. The equations of motion has the familiar form: d d t ∂ L ∂ q ̇ − ∂ L ∂ q + A q T λ = Q (1) where q ∈ R 5 ( N − 1 ) is the configuration vector containing α ₁⋯ϕ _N−1 variables, ∂ L ∂ q is the partial derivative of the Lagrangian with respect to each generalized variable, Q ∈ R 5 ( N − 1 ) is the generalized forces due to the internal pressure acting on the surface of the last disc, and A q ∈ R 5 ( N − 1 ) × 5 ( N − 1 ) is the Jacobian of the Pfaffian constraints arranged so that it can be multiplied by the vector of Lagrange multipliers λ ∈ R 5 ( N − 1 ) . In this work, we are trying to observe the robot configuration (and contact force, if applicable), whether the robot is in free space or in contact. In contact the product A q T λ will be nonzero and we replace it by F _c , where “c” denotes the effects of the contact force on the robot’s configuration space. This term will have the same units as the generalized force Q. By using the canonical momenta vector p = ∂ L ∂ q ̇ = M q ̇ with the mass matrix M ∈ R 5 ( N − 1 ) × 5 ( N − 1 ) , we can rewrite the above equation of motion in the state space form as: q ̇ p ̇ = 0 M − 1 0 0 q p + 0 ∂ L ∂ q + Q − F c (2)

The inflatable soft robot described in Section 2.1 will operate at a range of pressures. Because the analytical expressions of M, ∂ L / ∂ q and Q are functions of generalized coordinate q and the input pressure P, the equations of motion of the soft robot are not constant but change corresponding to each pressure (operating point). Note that the expressions of these matrices are complicated because they include partial derivatives of expressions of many variables and they are varying nonlinearly with pressure. Therefore, if we evaluate M, ∂ L / ∂ q and Q at a particular pressure and obtain the numerical values at that operating point, we can retrieve an affine linear-time-invariant state space that represents the soft robot nearby that operating point. In order to describe the robot throughout its entire range of motion, we divide the entire system into K subsystems fixed to a sequence of operating points, by evaluating the constituent matrices at corresponding pressures. This action forms a switched LPV (linear parameter-varying) system. Note that input pressure is chosen as the switching condition because it is the only input that changes the robot configuration. It is also measured easily and we don’t have to know the robot pose explicitly to select the correct operating points. The LPV system consists of a set of linear systems whose state-space operational modes are driven by an underlying decision process based on the current input pressure. The LPV system is expressed as the linear state space system: x ̄ ̇ ( t ) = A v ( u ) x ̄ ( t ) + B v ( u ) u ( t ) + Θ v ( u ) y ̄ ( t ) = C x ̄ ( t ) (3) where x ̄ = q p T = x 1 x 2 T , the output y ∈ R m consists of position measurements from motion capture markers attached to the body of the robot. The output is related to the states by C = [ C 1 m × m 0 ] ∈ R m × n which results from the linearization of the nonlinear functions mapping the measured data to the world frame. The switching rule v(u) ∈ S ^K ≐{1, 2, …K} depends on input pressure P. For each j ∈ S ^K , the subsystem matrices A _j , B _j and Θ _j are evaluated at the q corresponding the pose at the pressure P _j , resulting in constant matrices and have the forms: A j = 0 M j − 1 0 0 B j = 0 Q j Θ j = 0 ∂ L ∂ q j − F c j (4) The number of operating points and their distribution are carefully selected so that the switched LPV system can adequately represent the whole continuous system. The system is considered under the following Assumption.

Assumption 1 The generalized contact forces in Eq. 4 satisfy the following conditions: F _cj has a derivative and both are upper bounded as ‖F _cj ‖ ≤ L _j , ‖ F ̇ c j ‖ ≤ L j .

Before continuing with the discussion of the observer we discuss the hysteresis behavior of the soft robot. Hysteresis is defined as a “special type of memory-based relation between the input and the output” (Macki et al., 1993) or “rate-independent memory effect” (Visintin, 1994). Although there have been variations of the definition, hysteresis is generally used to describe a dynamical system with a phase lag that depends on the input. Hysteresis systems can be found in different areas such as physics, chemistry, engineering or economics. For our robot, we conducted experiments and observed that since our robot is made of silicone and is stretchable, its behavior is described as a hysteresis loop where the body of the robot follows different curves depending on the changing direction of the input pressure. Specifically, by using the motion tracking system to inspect the tip position in the Y direction of the robot, we could see that the tip followed the red curve while the chamber was inflating and followed the blue curve while the chamber was deflating, as shown in Figure 3. Note that the curves are generated by polynomial regressions using data (the discrete points) collected from the motion tracking system. It can be seen that the deflating curve has a higher slope than that of the inflating curve. In other words, when deflating, the robot can reach a desired position faster than when it is inflating. This behavior of the inflatable robot is interesting but it increases the complexity in the system observation and control design.

Since the system configuration depends on the direction of the input pressure, from the idea of an LPV system in Eq. 3, we can express it as a generalized Duhem model as follows (Visintin, 1994): x ̄ ̇ ( t ) = A v ( u ) x ( t ) + B v ( u ) u ( t ) + Θ v ( u ) g ( u ̇ ( t ) ) y ̄ ( t ) = C x ̄ ( t ) (5) where g ( w ) = h + w if w ≥ 0 h − w if w < 0 (6) where h ⁺ and h ⁻ are different functions corresponding to the direction of the input signal.

More specifically, this hysteretic system can be writen as an SPV system as: x ̄ ̇ ( t ) ≜ x ̇ ( u ) = A v ( u ) + x ( t ) + B v ( u ) + u ( t ) + Θ v ( u ) + if u ̇ ≥ 0 A v ( u ) − x ( t ) + B v ( u ) − u ( t ) + Θ v ( u ) − if u ̇ < 0 y ̄ ( t ) ≜ y ( u ) = C x ̄ ( t ) (7) This SPV system contains twice the number of subsystems in Eq. 3, including a half subsystem to describe the LPV system while the robot is inflating and the other half to represent the system in the deflating direction of the input. The switching action takes place when the control input changes its direction.

In the current formulation, we have to do backward-in-time identification for the hysteresis curve and the input pressure, which are not suitable when designing conventional observers and controllers. To avoid such kind of identification and checking the derivative of the input, we applied a transformation by introducing a monotonically increasing independent variable, u ̄ ∈ [ u m i n , 2 u m a x − u m i n ] . The transformation is inspired by the work of Oh and Bernstein (2005). This helps us to rewrite the model as a new semilinear Duhem model as follows: x ̇ ( u ̄ ) = A v ( u ̄ ) + x ( u ̄ ) + B v ( u ̄ ) + u ̄ + Θ v ( u ̄ ) + if u m i n ≤ u ̄ ≤ u m a x A v ( u ̄ ) − x ( u ̄ ) + B v ( u ̄ ) − u ̄ + Θ v ( u ̄ ) − if u m a x < u ̄ ≤ 2 u m a x − u m i n y ( u ̄ ) = C x ( u ̄ ) (8) where y ̄ ( u ̄ ) ≜ y ̄ + ( u ̄ ) if u m i n ≤ u ̄ ≤ u m a x y ̄ − ( u m a x + u m i n − u ̄ ) if u m a x ≤ u ̄ < 2 u m a x − u m i n (9) The new hysteresis curves describing the tip of the robot are displayed in Figure 4 where the inflation curve is the red line and the deflation curve is the blue line. Compared to the version in Figure 3, the deflating curve is now flipped over so that the two curves form a continuous function with respect to the new monotonically increasing input u ̄ . A system identification process has been performed to define each matrix in the SPV system in Eq. 8 at each working point to complete all of its subsystems. Then the SPV model can be used for observer and controller design.

Due to the deformable body of the soft robot, the numerical values of q, p and F _c in Eq. 2 can not be obtained directly from measurements from the sensors we are using. The model of the soft robot in this study includes a large number of state variables, and we need an observer to estimate both states and contact forces. For these reasons, we selected the high-order sliding mode observer, because it has a simple rule to select the observer gains, no matter the size of the state space in observation; and based on our knowledge, it is the best observer that works well for both states and disturbance observation. The other alternatives are more challenging to be applied to our system. The following Lemma is re-stated to summarize the result of finite-time stability of the dynamical system in the study of Levant and Livne (Levant, 2003), which is applied to develop the observer:

Lemma 1 Consider the following system: ε ̇ 0 = − λ 0 | ε 0 | n / n + 1 s i g n ( ε 0 ) − η 0 ε 0 + ε 1 , ε ̇ 1 = − λ 1 | ε 1 − ε ̇ 0 | ( n − 1 ) / n s i g n ( ε 1 − ε ̇ 0 ) − η 1 ( ε 1 − ε ̇ 0 ) + ε 2 , … ε ̇ n − 1 = − λ n − 1 | ε n − 1 − ε ̇ n − 2 | 1 / 2 s i g n ( ε n − 1 − ε ̇ n − 2 ) − η n − 1 ( ε n − 1 − ε ̇ n − 2 ) + ε n , ε ̇ n = − λ n s i g n ( ε n − ε ̇ n − 1 ) − η n ( ε n − ε ̇ n − 1 ) − 1 L 0 f ( t ) (10) where n is the relative degree, ɛ ₀, …, ɛ _n are the state variables, λ ₀, …, λ _n and η ₀, …, η _n are appropriate positive scalar constants and the perturbation f(t) satisfies the condition |f(t)| ≤ L ₀ with L ₀ a proper positive constant. Then the system converges to the origin in finite time.

Given position data from the marker set, we built a set of observers based on the third order sliding mode approach of Levant (2003), which has fast convergence and high robustness. The number of observers is the same as the number of subsystems in Eq. 8 The observer for subsystem j is of the form: x ̂ ̇ 1 j = M j − 1 x ̂ 2 j − k 1 j L j 1 / 3 ‖ x ̂ 1 j − x 1 j ‖ 2 / 3 s i g n ( x ̂ 1 j − x 1 j ) − k 2 j ( x ̂ 1 j − x 1 j ) x ̂ ̇ 2 j = − k 3 j L j 1 / 2 ‖ M j − 1 x ̂ 2 j − x ̂ ̇ 1 j ‖ 1 / 2 s i g n ( M j − 1 x ̂ 2 j − x ̂ ̇ 1 j ) − k 4 j ( M j − 1 x ̂ 2 j − x ̂ ̇ 1 j ) + ∂ L ∂ q j + Q j + F ̂ c j F ̂ ̇ c j = − k 5 j L j s i g n ( M j − 1 x ̂ 2 j − x ̂ ̇ 1 j ) − k 6 j L j ( M j − 1 x ̂ 2 j − x ̂ ̇ 1 j ) (11) where x ̂ 1 j , x ̂ 2 j are vectors of estimated states, k _1j to k _6j are the observer gains to be designed and F ̂ c j is the vector of estimated generalized contact forces. Note that x _1j can be obtained from the multiplication of the inverse of C ₁ and y. Since we have more outputs than the number of states (one marker provides three outputs including the 3D coordinate of that marker in X, Y, Z direction), we can form a non-singular square matrix C ₁ so that it is invertible.

Theorem 1 For system Eq. 8, if the observer set is designed as in Eq. 11 and the observer gains are selected properly, then the estimated states and contact forces will converge to the true values in finite time, which follows from Lemma 1 (Levant, 2003).

Proof. The estimation error variables are defined by the following formulation: ϵ 1 j = x ̂ 1 j − x 1 j L j ϵ 2 j = M j − 1 ( x ̂ 2 j − x 2 ) L j = M j − 1 x ̂ 2 j − x ̂ ̇ 1 j L j = M j − 1 ( x ̂ 2 j − M x ̂ ̇ 1 j ) L j ϵ 3 j = F ̂ c j − F c j L j (12) the dynamics of the estimation errors are then obtained as: ϵ ̇ 1 j = 1 L j − k 1 j L j 1 / 3 ‖ L j ϵ 1 j ‖ 2 / 3 s i g n ( ϵ 1 j ) − k 2 j ( L j ϵ 1 j ) + M j − 1 x ̂ 2 j − M j − 1 x 2 = − k 1 j ‖ ϵ 1 j ‖ 2 / 3 s i g n ( ϵ 1 j ) − k 2 j ( ϵ 1 j ) + ϵ 2 j , ϵ ̇ 2 j = 1 L j − k 3 j L j 1 / 2 ‖ L j ϵ 2 j ‖ 1 / 2 s i g n ( ϵ 2 j ) − k 4 j ( L j ϵ 2 j ) + F ̂ c j − F c j , = − k 3 j ‖ ϵ 2 j ‖ 1 / 2 s i g n ( ϵ 2 j ) − k 4 j ( ϵ 2 j ) + ϵ 3 j ϵ ̇ 3 j = − k 5 j s i g n ( ϵ 2 j ) − k 6 j ϵ 2 j − 1 L j F ̇ c j (13) If the conditions in Assumption 1 are satisfied, it follows from Lemma 1 that system Eq. 13 is finite-time stable, which implies that the estimation variables converge to true variables in finite time. According to the rule recommended by Levant (2005), the convergence can be guaranteed by defining the gains as k _1j ≥ 3, k _3j ≥ 1.5, k _5j ≥ 1.1 for this third-order system. Note that the larger the gains, the faster the convergence and the higher sensitivity to input noises and the sampling step. k _2j , k _4j , k _6j can be selected through trials in the simulation.

The block diagram in Figure 5 summarizes the design process for the robot systems and the observer. The observer block stands for the sub-observer j serving the corresponding subsystem j working at a certain pressure value.

The soft robot and position of the fabric can be seen in Figure 1 with the detailed dimensions (in mm) of the undeformed shape are given in Figure 6. To apply the Disc-Thread modeling into the simulations and experiments, our soft robot is discretized into four discs and three threads. The first and last discs are at the base and the tip of the robot, correspondingly. The second disc is considered to go through two points O ₂ and O 2 ′ and the third disc is consider to go through O ₃ and O 3 ′ (see Figure 7). The thread lengths are set at ℓ ₁ = 7 cm, ℓ ₂ = 9 cm and ℓ ₃ = 7 cm. Note that the robot can be modeled with a higher number of discs but it will result in longer expressions without much improvement in model accuracy. For a robot of this size and aspect ratio, this four-disc model can adequately describe its behavior. We selected four operating points at different increasing pressures at (1.5, 2.5, 3, 3.5) psi (see Figure 1) and evaluate M, ∂ L / ∂ q and Q at these points to generate four corresponding subsystems that represent the entire working range of the soft robot during its inflation. The set of these matrices are provided in the supplemental data. Note that these four subsystems are the same for the LPV system and the inflating branch of the SPV system. Considering that working with the whole SPV system is too long for this paper, and it does not improve the illustration of the observer performance, we only use the inflating branch of the SPV system for our simulation and experiment. A set of sub-observers is designed using the form of Eq. 11. The observer gains are chosen as k ₁ = 3, k ₂ = 4, k ₃ = 1.5, k ₄ = 3, k ₅ = 1.1, k ₆ = 2 for all sub-observers. The disturbance is assumed to be bounded by L = 20. At the run-time of the robot, when a new subsystem is switched on due to the changing pressure, the new corresponding sub-observer with updated system matrices is switched on as well. The switching rule is illustrated in Figure 8.

The proposed observer was validated by simulations and experiments. The simulations were performed in Matlab Simulink using S-functions for the dynamics of the plant and observer. The experiments were conducted in our laboratory environment. They are described in detail as follows:

Simulation 1: We examined the ability of the observer to estimate the space variables q _i at four operating points during quasistatic motion in free space for the whole continuous system.

Experiment 1: We conducted this experiment to gather data from a real contact event between the soft robot and its surrounding environment. The data will be used to validate the observer performance after it estimates the generalized contact force. In order to measure the real contact forces, we built a system that includes a load cell connected to an acrylic plate on one end and connected to a frame on the other (see Figure 7). The load cell is an Interface 6A27 which can measure forces simultaneously in three mutually perpendicular axes and three simultaneous torques about those same axes. Initially, the robot was kept in a stable position and the plate was located at 3.38 cm above the robot. We inflated and deflated the robot with the input pressure history shown in Figure 9. We observed that at first the robot raised its head and then bent down and twisted (as can be seen in Figure 1) when inflated, and returned to its original pose when deflated. The input pressure and the initial gap between the plate and the robot were chosen so that the robot was in contact with the plate surface at the location corresponding to disc three when the chamber pressure was in the range [1.5–3.3] psi. During the contact event, the contact forces in each direction were recorded by the load cell. Then a set of 15 generalized contact forces (which have units of torque) were calculated using F c = J b c T F (where F is the three components of the contact forces measured by the load cell and J _bc is the relative Jacobian between the base and the contact point). These values were the ground-truth generalized contact forces. To simulate the observer estimating generalized contact forces, the measured forces were considered as disturbances to the robot model. When we ran the Simulink model of the plant and the observer, the observer could detect the disturbances (generalized contact forces) and produced the close estimates of these same measured forces.

Experiment 2: In this experiment, the robot was moving in free space without any contact with the environment. We inflated and deflated the robot with the same input pressure history shown in Figure 9. Since the disc-thread model represents a discretization of a continuum soft robot, the discs themselves are fictitious and do not exist in the physical system. Thus, there is no ground-truth value against which to directly compare the state estimates (α ₁ to ϕ ₃). However, we can compute what the marker locations should be based on the estimated states and compare them to the measured marker locations. The poses assumed by the robot during its movement were recorded by seven markers of the Polhemus electromagnetic motion tracking system. The 3D location of the markers on the body of the robot are considered to be the measured output of this robot system. We then feed the data measured from the markers into the observer into the Simulink model to estimate the soft robot poses.