This paper introduces a general Lie group framework for modeling continuum soft robots, employing Cosserat rod theory combined with cumulative parameterization on the Lie group SE(3). This novel approach addresses limitations present in current strain-based and configuration-based methods by providing geometric local control and eliminating unit quaternion constraints. The paper derives unified analytical expressions for kinematics, statics, and dynamics, including recursive Jacobian computations and an energy-conserving integrator suitable for real-time simulation and control. Additionally, the framework is extended to handle complex robotic structures, including segmented, branched, nested, and rigid-soft composite configurations, facilitating a modular and unified modeling strategy. The effectiveness, generality, and computational efficiency of the proposed methodology are demonstrated through various scenarios, including large-deformation rods, concentric tube robots, parallel robots, cable-driven robots, and articulated fingers. This work enhances modeling flexibility and numerical performance, providing an improved toolset for designing, simulating, and controlling soft robotic systems.
Figures & tables
Symbol
Unit
Definition
˙
–
Derivative with respect to time
′
–
Derivative with respect to space
–
Converts R6 in se(3)
s
m
Reference arclength coordinate
t
s
Time
p(s,t)
m
Position vector of cross-section
Table 1: Notation and definitions.
Figure 1: Schematic diagram of the soft manipulator.
Figure 2: Bspline and Hermite curve on SE(3) manifold.
Figure 3: B-spline basis functions.
Figure 4: Cumulative basis functions.
Parameterization
Original basis Bi(⋅)
Cumulative basis Bˉi(⋅)=j=i∑nBj(⋅)
Power (monomial)
Bi(s)=si
Bˉi(s)=j=i∑nsj=1−ssi−sn+1(s=1)
B - spline (degree k )
Bi(s)=Ni,k(s)
Bˉi(s)=j=i∑nNj,k(s)
Cubic Hermite † (Bézier)
Bi(s)=Ni,3(s)
Bˉi(s)=j=i∑3Nj,3(s)
Table 2: Original and cumulative basis functions for several nonlinear parameterizations
Figure 5: The variation of the kinematic Jacobian matrix concerning the arc-length coordinate s . The left figure shows the Jacobian matrix obtained using the Piecewise Linear Strain (PLS) [ 28 ] parameterization method. In contrast, the right figure shows the Jacobian matrix obtained using the proposed Lie group B-spline parameterization method.
Figure 6: Schematic of a tree-structured robot. Each ellipse represents a soft or rigid body segment (e.g., bodyi ), and the small circles represent joints that allow rotational movement. The bodyi acts as a branching point, connecting to multiple child links bodyi+1 and bodyi+2 , illustrating a typical tree structure.
Figure 7: Three types of connections. (a) Direct rigid connection; (b) Rigid body insertion; (c) Movable joint connection.
Figure 8: Relative motion modeling in concentric tube robots. (a) General case: Rod 1 and Rod 2 are connected via a relative motion gr(s)∈SE(3) , which includes both translation and rotation. (b) Concentric configuration: The rods are constrained to have no relative translation, resulting in a pure rotational relationship Rr(s)∈SO(3) along their shared centerline.
Figure 9: Schematic of a closed-loop robotic structure. Each ellipse represents a rigid body (e.g., bodyi ), and the small circles denote rotational joints. Unlike tree-structured robots, this configuration forms a closed kinematic loop: starting from bodyi , connections propagate through bodyi+1 , bodyi+3 , and bodyi+2 , and return to bodyi . Such closed-chain mechanisms are common in parallel robots and constrained articulated systems.
Constraint Type
Constraint Function [1]
Jacobiant w.r.t. g1 [2]
Jacobiant w.r.t. g2 [3]
Fixed constraint
b=g2⊟g1
∇g1b=dτb−1Adτ(b)−1
∇g2b=dτb−1
Articulated constraint
b=P1(g2⊟g1)
∇g1b=P1dτb−1Adτ(b)−1
∇g2b=P1dτb−1
Sliding constraint
b=P2(g2⊟g1)
∇g1b=P2dτb−1Adτ(b)−1
∇g2b=P2dτb−1
Table 3: Constraint functions and their Jacobians for different constraint types
Figure 10: Integration flowchat. The discrete Euler–Lagrange equation is solved to determine the next configuration based on the previous and current ones. This process is iteratively applied to compute the entire trajectory.
Figure 11: Workflow of the proposed modeling method. The process includes control point parameterization, kinematics computation, connection updates, system assembly, and static/dynamic solving.
Parameter
Value / Unit
Length
1.0 m
Front end radius
3.0 cm
Rear end radius
1.5 cm
Young’s modulus
2.0×105 Pa
Possion’s ratio
0.45
Table 4: Physical Parameters
Figure 12: Comparison of deformation results between the proposed method and the ground truth. Blue rods represent the results from the proposed method, while gray rods correspond to the ground truth. In the upper panel, an upward vertical force is applied at the right end of the rods, increasing from 0 N to 1 N in increments of 0.25 N. In the lower panel, a moment about the negative y -axis is applied at the same location, increasing from 0 Nm to 0.2 Nm in increments of 0.05 Nm.
Figure 13: Error between the simulation results and the ground truth along the length of the soft rod under different end-tip loads.
Figure 14: Convergence behavior of the static solver. The plot shows the residual norm ∥r∥ over iterations during the static equilibrium computation.
Figure 15: Convergence of the configuration error with increasing number of control points under B-spline basis functions of order 2 to 5.
Figure 16: Initial equilibrium configuration of the soft rod. The rod is fixed at its left end and subject to a 10 N upward force at the right end, resulting in a curved static equilibrium shape. This configuration is used as the initial condition for subsequent free oscillation analysis.
Parameter
Value / Unit
Length
1.0 m
Front end radius
5.0 cm
Rear end radius
3.0 cm
Young’s modulus
2.0×106 Pa
Possion’s ratio
0.45
Viscosity modulus
0 Pa ⋅ s
Table 5: Physical Parameters
Figure 17: Time evolution of the rod tip position along the z -axis. After the external force is removed, the rod undergoes sustained vertical oscillation.
Figure 18: Total energy evolution of the system over time. The total energy is computed as the sum of kinetic and potential energy. The upper plot shows simulation results obtained using the implicit Euler integrator, while the lower plot presents results using the proposed symplectic integrator.
Figure 19: Deformation of the cable-driven soft manipulator. The left panel shows the initial configuration of the soft rod, with all cables untensioned and the structure in a straight, undeformed state. The right panel illustrates the deformed configuration resulting from the actuation of the cable. Three cables (red line) run through the vertices of the triangular plates, applying equal and opposite tension between the corresponding vertices of adjacent plates.
Parameter
Value / Unit
Length
1.0 m
Front end radius
4.0 cm
Rear end radius
2.0 cm
Young’s modulus
1.0×106 Pa
Possion’s ratio
0.45
Viscosity modulus
10 Pa ⋅ s
Table 6: Physical Parameters
Figure 20: The tree structure of the cable-driven soft manipulator. The cable-driven robot consists of soft elements and rigid triangular plates connected in series. Each pair of adjacent rigid plates is linked by a soft segment, with the connections between them modeled as rigid joints.
Figure 21: Concentric tube robot .The concentric tube robot consists of three tubes sharing a common central axis, with the third tube held stationary and fixed to the world coordinate frame.
Parameter
1st tube
2nd tube
3rd tube
Length (m)
0.6
0.4
0.05
Radius (mm)
2.5
5
50
Young’s modulus (GPa)
2
2
10
Possion’s ratio
0.45
0.45
0.45
Pre-curvature (1/m)
6
4
0
Table 7: Physical Parameters
Figure 22: The concentric structure of the three-tube concentric robot. The first (innermost) tube, serving as the main tube, possesses independent degrees of freedom. The remaining tubes are subordinate, with their degrees of freedom determined by the configuration of the main tube and their relative rotations with respect to it.
Figure 23: Representative deformation stages of the concentric tube robot. β denotes the distance from the first tube’s rotation tip to the bending entry point, and l is the overlap length of the bent sections. Initially, the tubes remain stable. At t=4.16 s, the snapping event induces instability, followed by a return to stability at t=4.35 s.
Figure 24: S-curve. Comparison between the simulation results obtained using the proposed method and the analytical solution.
Figure 25: Chiral structure under compressive loading. The left image shows the initial unloaded state, in which straight vertical rods are symmetrically arranged between three horizontal circular plates. These rods are stress-free and maintain parallel alignment without deformation. The right image illustrates the deformed state under vertical compressive loading. The rods exhibit significant twist-buckling behavior, undergoing compression, bending, and torsion simultaneously.
Parameter
Value / Unit
Length
0.8 m
Radius
4.0 cm
Young’s modulus
1.0×105 Pa
Possion’s ratio
0.45
Table 8: Physical Parameters
Figure 26: Structure of the simulated parallel soft–rigid mechanism and its Jacobian sparsity pattern. The left panel illustrates the topological model of a parallel mechanism composed of rigid bodies ( R1 – R3 ) connected by soft rods ( S1 – S8 ). Nodes represent components, while edges represent geometric constraints between them. The right panel shows the sparsity pattern of the Jacobian matrix derived from the static KKT system, which includes both the Hessian and the Jacobians of all constraints. The structured sparsity reflects the modular and localized nature of the hybrid soft–rigid system.
Figure 27: Strain–force response under different initial torques. The plot illustrates the vertical compressive force versus strain relationship for chiral structures with varying initial torques (0.1 Nm, 0.5 Nm, 1 Nm, and 2 Nm). Lower initial torque results in a stiffer response with rapid force increase, while higher torques induce more pronounced twist-buckling and bending, allowing for greater strain accommodation and enhanced energy absorption performance.
Figure 28: Deformation of a cable-driven soft finger. The left panel shows the initial undeformed configuration of the soft finger, where the rigid bone segments (yellow) are aligned and the soft joint rods (blue) remain straight with untensioned cables (red line). In the right panel, tensile forces applied to the embedded cables induce bending at the joints, causing the finger to flex.
Figure 29: Geometric design of the soft joint. The shape of the soft joint is parameterized using a cubic B-spline curve defined by five control points (frame). The left panel shows the B-spline curve (blue) along with the control polygon (dashed orange) and control points (frame), while the right panel displays the corresponding 3D revolved geometry used to construct the soft joint.
Parameter
Value / Unit
Thickness
1 mm
Width
20 mm
Young’s modulus
5 GPa
Table 9: Physical Parameters
Figure 30: Trajectory of the finger’s tip with the different soft joint combination. The labels ① and ② indicate two distinct natural shapes of the soft joints. Each subfigure corresponds to a different joint combination in the robotic finger.
Modeling soft robot dynamics is challenging due to their continuum structure and typically nonlinear dynamics. Creating models based on first-order principles is typically time-demanding, and their expressiveness is limited, whereas data-driven models lack interpretability and physical consistency. This work aims to overcome these challenges by introducing a port-Hamiltonian Gaussian Process Regression framework for learning and simulating the dynamics of planar, rod-like soft robots. In detail, the proposed model integrates Cosserat rod theory and Hamiltonian physics with data-driven inference to preserve the system's energy structure while accurately learning the rod dynamics. Numerical simulations show that we can achieve accurate and energy-consistent representations of a rod-like soft robot, showing the potential for a robust and interpretable pathway for modeling complex continuum mechanics.
Mohammad Ali, Nithin Senthur Kumar, Eric J. Barth +1
Department of Computer Science, Vanderbilt University, Nashville, TN 37235, USA · Department of Mechanical Engineering, Vanderbilt University, Nashville, TN 37235, USA
In this paper, we design feedback control laws for soft robots modelled using the Cosserat rod, which is spatially discretised using the Piecewise Constant Strain (PCS) approach. The PCS approach transforms the nonlinear PDEs describing the Cosserat rod to a system of nonlinear ODEs. This simplification results in a model describing soft robots which is similar to the serial rigid-link manipulators. We design feedback control laws for the quasi-static PCS model by using external wrenches as control inputs. The control laws are designed based on state-feedback linearisation in strain and task spaces. An extensive set of numerical results demonstrates the performance of the control laws for end-effector trajectory tracking and shape control of soft robots.
Srishti Siddharth, Domenico Campolo, Ravi Banavar
Centre for Systems and Control Indian Institute of Technology Bombay Mumbai, India · School of Mechanical and Aerospace Engineering Nanyang Technological University Singapore
Cosserat rod models for soft robots usually construct sectional stiffness by summing material properties over a common cross-section. This assumption becomes inaccurate for trimmed helicoid arms, where load-bearing helix domains are separated and connected only through sparse fused crossings. This paper formulates a separated-section constitutive law that evaluates each helix domain in its local frame and pulls its constitutive response back to the backbone, yielding an effective backbone stiffness. Sparse-fusion mechanics captures the additional compliance caused by relative motion between neighboring domains and determines channel-wise reduction profiles ηc(s/L) for bending, torsion, and extension. The resulting effective sectional stiffness is strongly anisotropic: bending and extension are reduced by about one order of magnitude, whereas torsion remains close to the effective backbone stiffness. The resulting sectional law is embedded in a geometrically exact dynamic Cosserat model with GVS discretization and routed-tendon actuation. Across 103 measured configurations, the three datasets give pooled normalized position errors of 7.7%, 6.7%, and 7.8%, while each full-arm solve requires approximately 0.3s on one CPU core (Intel Xeon, Cascade Lake, 2.8GHz), enabling rapid model-based planning, state and load estimation, and morphology--control co-design for architected soft robots.
Zhihang Qin, Linxin Hou, Zeyu Zhong +10
Department of Mechanical Engineering, College of Design and Engineering, National University of Singapore, Singapore · Advanced Robotics Centre, National University of Singapore, Singapore · Department of Electrical and Computer Engineering, National University of Singapore, Singapore +1