Mimicking highly dexterous and deformable biological bodies has been a trending topic of multi-disciplinary research, called soft robotics, using intrinsically soft materials in the form of continuum robotic platforms (Rus and Tolley 2015). Performing delicate tasks (Cianchetti et al. 2014), high maneuverability in unstructured and confined environments (Burgner-Kahrs et al. 2015;Cianchetti and Menciassi 2017;Walker et al. 2016), dexterous grasping (Katzschmann et al. 2015), mimicking biological tissue and organs (He et al. 2018), bioinspired dynamic locomotion (Wehner et al. 2016) such as crawling (Rich et al. 2018), terrestrial (Godage et al. 2012) or submerged locomotion (Cianchetti et al. 2015) are among the promises made by the researchers in the field. Soft robots are appealing to investigate new design and theoretical concepts such as variable stiffness structures (McEvoy and Correll 2018), morphological computation (Nakajima et al. 2018) and embodied intelligence (Nakajima et al. 2015), to simplify the control and sensing tasks through robot embodiment (Füchslin et al. 2012;Thuruthel et al. 2018a).
However, compliance has disadvantages, such as uncertain deformations, limited control feedback, reduced control bandwidth, stability issues, underdamped modes, and lack of precision in tasks involving working against external loads (Blanc et al. 2017;Cianchetti et al. 2013). These usually result in modeling and/or control challenges for such designs. Besides, there is an urgent need for unified frameworks to transfer our well-established knowledge of dynamic system analysis, path planning and control design for rigid-body robots to soft robotics research (Kapadia et al. 2010;Renda and Seneviratne 2018;Renda et al. 2018;Della Santina et al. 2018b) and to model hybrid rigid-soft body systems (Sadati et al. 2018c;Paternò et al. 2018). Such frameworks should be as simple as possible and easy to use, to be widely accepted by the ever-growing soft robotics research community that gathers researchers from different disciplines and backgrounds. It should provide fast computational performance, to be suitable for control and design problem of soft systems with large state spaces. Also, it needs to be integrable with standard software platforms, e.g. C/C++ language, Matlab software, ROS (Robotic Operation System, see ros.org), etc., widely used in the community.
Here, we introduce two new modeling approaches for continuum rods and actuators, a general reduced-order model (ROM), and a discretized model with absolute states and Euler-Bernoulli beam segments (EBA). These models enable us to perform more accurate simulation of continuum rod manipulators as well as extending the solution to modeling 2D and 3D continuum geometries, which is missing in similar recent research (Renda et al. 2018). Besides, a new formulation is presented for a recently introduced discretized model by (Renda et al. 2018;Shiva et al. 2018) which is based on Euler-Bernoulli beam theory and relative states (EBR). These models are implemented in a Matlab software package, that we name T M T Dyn, to establish a new modeling and simulation tool for hybrid rigid-continuum body systems. The package is improved with a new High-Level Language (HLL) text-based interface, an elementwise representation of deriving Lagrange EOM in a vector formalism (called TMT method (Wisse and Linde 2007)), a CAD-file import module, automatic formation of the system EOM for different modeling and control tasks, implementing C-mex functionality for improved performance, and other modules for static and linear modal analysis of a hybrid system. Our main goal is to make the tasks of deriving EOM of hybrid rigid-continuum body robots, performing dynamic system analysis, state observation, and control system design more accessible to the interdisciplinary soft robotics research community and people with limited expertise in dynamic system modeling.
In this paper, first, in a brief state-of-the-art review, we discuss High-Level Languages in robotic system modeling, system dynamics modeling, and mechanics of continuum structures with a focus on continuum rods as the different elements of this research. Then, we discuss hybrid system kinematics, based on unit and non-unit quaternions, and how rigid-body kinematics framework can be extended to the variable curvature and discretized continuum kinematics cases, with relative and absolute states. A new general yet efficient reduced-order solution for the rod backbone is discussed based on truncated series. Then, an elementwise form for the TMT method is presented to derive a separate set of equations for each body/element in a dynamic system. The Principle of Virtual Work (PVW) is adapted to derive the system linear and nonlinear compliance actions. Using this framework, hybrid system dynamics is discussed where a lumped-system representation of the Cosserat rod theory and a new discretization method based on absolute (independent) states is presented. As a result, the relations for a one-dimensional continuum element can be generalized to model two-dimensional (membrane or fabric) and three-dimensional (tissue block) geometries. A simple example is employed to showcase the derivation steps and the use of T M T Dyn package text-based user interface.
High-Level and Domain-Specific Languages Domain-Specific Languages (DSLs) and Model-driven Engineering (MDE) are interesting emerging areas in the robotics research community, e.g. distributed robotics, system control, and vision, with significant potential in facilitating the programming of future robots. A DSL is a dedicated programming language for a particular problem domain, offering specific abstractions and notations, to decrease the coding complexity and increase programmer productivity. DSLs have been used for programming complex systems, such as robots, control systems, etc, for which traditional general-purpose languages do not provide a good correlation between the implementation requirements and language features. To address this, DSLs are powerful and systematic ways to provide two main features;
1. quick and precise adaptation by domain experts, who are not familiar with general-purpose programming languages; 2. hiding the architecture complexity by software engineers to facilitate complex configuration and design architectures before transferring to domain experts (Fowler and Parsons 2011). From this point of view, DSLs are often High-Level Languages (HLL) introducing a higher level of abstraction suitable for the intended application.
ROS, Orocos (Open Robot Control Software, see orocos.org), SmartSoft (see smart-robotics.sourceforge.net), OpenRTM (see openrtm.org), Matlab Robotics System Toolbox (see mathworks.com/products/robotics), and Robotics Toolbox (see petercorke.com/wordpress/toolboxes/roboticstoolbox) are some robotic software platforms developed to make robotics programming and configuration as accessible as possible to experts from different application domains. Robotics Toolbox by Corke (2019) is probably the most comprehensive text-base toolbox developed in Matlab software environment. The Toolbox, currently in its 10 th release, provides many functions for kinematics, dynamics, and trajectory generation of classical arm-type robotics based on methods in (Corke 2017).
However, none of the aforementioned robotics software platforms are fully compatible with soft structure robots, which need modules for handling highly articulated geometries with repeated elements, discretization or employing reduced-order model assumptions to simplify the modeling and control state space for a continuum geometry. Moreover, it is not always straightforward to access the derived Equations of Motion (EOM) which is necessary for off-the-shelf dynamic system analysis, design optimization, and control system design tasks. On the other hand, most recent efforts for publishing modeling packages for soft robots (Coevoet et al. 2017;Gazzola et al. 2018;Renda and Seneviratne 2018;Hu et al. 2018) were more focused on the feasibility of their modeling and control approach, and cannot be extended easily due to not providing the classical EOM Prepared using sagej.cls for the system (mostly because of using differential equations and Finite Element Method (FEM) solvers, e.g. SOFA (Simulation Open Framework Architecture, see www.sofaframework.org) package (Coevoet et al. 2017) and Chain-Queen (Hu et al. 2018)). Almost none of them provide an easy-to-use HLL to reach researchers with less expertise in the mechanics of continuum structures.
As a part of this research, we address these challenges by developing a Matlab-based package for modeling and control of hybrid robotic systems with minimal modeling and control states. To this end, the key challenge is to integrate continuum mechanics with traditional rigid body dynamics, which is widely being used in robot control research. To address this issue, we start with a brief review of the modeling methods for mechanical systems and, later, continuum elements.
Deriving EOM for modeling hybrid rigid-continuum body mechanisms has been a challenge in soft robotic research (Boyer 2014;Rus and Tolley 2015;Burgner-Kahrs et al. 2015;Sadati et al. 2017b). The Newton-Euler method (Jung et al. 2011), Lagrange dynamics (Mustaza et al. 2019;Godage et al. 2016), the Principle of Virtual Work (PVW) (Trivedi et al. 2008;Sadati et al. 2016Sadati et al. , 2017a)), TMT methods for deriving the Lagrange dynamics (Sadati et al. 2018b) and Bond Graph approaches (Sutar and Pathak 2017) have been used to derive EOM of such systems. Here, we disregard the Cosserat rod dynamics (Till and Rucker 2017), FEM (Coevoet et al. 2017) and the Moving Least Squares Material Point Method (MLS-MPM) (Hu et al. 2018), methods that result in a system of differential equations which are not inherently compatible with control design methods for rigidbody systems.
Most commercially available dynamical system modeling software, e.g. MSC. ADAMS (Negrut and Dyer 2004) and Robotics Toolbox (Corke 2019), utilize the Lagrange dynamics formulation. Lagrangian method results in a hardto-interpret final set of equations. Hence, usually, an extra step is needed to collect the final EOM in a closed-form vector formalism (see Appendix 2).
Recently, Wisse and Linde (2007) have employed a vectorform of Lagrange dynamics to model a class of passive biped walkers in which the coefficient matrices of the EOM vector formalism are derived directly and mostly separate from each other. According to Wisse and Linde (2007), "this method is based on parts of Lagrange's investigations on variational calculus and analytical mechanics, dealing with generalized coordinates, virtual work and inertial forces (presented earlier by D'Alembert in 1743), which was published in 1788 and before his well-known Lagrange method". Wisse and Linde (2007) have named it "TMT", not as an acronym but because of the easy formation of the system inertial matrix (M ) final form in the generalized coordinates space (T M T ). Here, T is the Jacobian transformation matrix between the Cartesian and generalized coordinates spaces and superscript ( ) is the matrix transpose operator.
Although it is not a unique feature of the TMT method (many standard approaches for dynamic system analysis result in the same term, e.g. D'Alembert, Lagrangian, Recursive-Newton Euler, Gauss-Principle, Jourdain's Principle, etc. method (Lynch and Park 2017)), we continue with this choice of name for simplicity and consistency with previous research. The main advantage of this method compared to the standard Lagrange method is that a final step for collecting coefficients of the system states and arranging them in a vector-form is not needed which tends to be complex and time consuming for large systems. Furthermore, the derived terms are independent, hence suitable for parallel analytical derivation and numerical simulation. As a result, it is faster to derive a system EOM in a closed vector form using commercially available symbolic mathematical toolboxes, e.g. Matlab Symbolic Toolbox, Maple, etc. The derivation steps of TMT method are compared with the Lagrange method in Appendix 2.
We have recently developed a Matlab software package, called AutoT M T Dyn to derive TMT EOM of rigid-body mechanisms (Sadati et al. 2018d(Sadati et al. , 2015)). AutoT M T Dyn * was originally developed for deriving the TMT EOM of rigid-body systems (Sadati et al. 2015) and used for analyzing free-fall righting maneuvers of a robot cat (Sadati and Meghdari 2017), lumped system modeling of continuum appendage (Sadati et al. 2017b), and dynamic analysis of a spider web structure (Sadati et al. 2018a).
In this paper, we introduce a new version of AutoT M T Dyn, now called T M T Dyn † , which is equipped with a new HLL text-based interface, CAD-file import module, automatic formation of the system EOM for different modeling and control tasks, implementing Matlab C-mex functionality for improved performance, and modules for static and linear modal analysis of a hybrid system. We employ the TMT method to derive EOM of continuum bodies based on discretized and reduced-order solutions. As a result, a unified software package is provided for deriving EOM, control design, and numerical simulation of hybrid rigid-continuum body systems. To this end, a brief review is provided on different modeling elements and assumptions for continuum rods.
If taking a theoretical approach, two key stages can be identified which determine the modeling strategy of a soft robot (Burgner-Kahrs et al. 2015;Gazzola et al. 2018;Sadati et al. 2017b). (1) Modeling assumptions for (1-I) the system kinematics, (1-II) system conservation law (system mechanics), and (1-III) material constitutional law (material mechanics). This stage results in a system of differential equations for the mechanics of a continuum medium. (2) The method to solve the resulting system, which can be based on (2-I) direct or (2-II) indirect methods. Please note that by a theoretical approach, we mean in contrast to pure learning (Braganza et al. 2007;Thuruthel et al. 2018b), combined reduced-order solution and learning (Thieffry et al. 2018a) approaches, and beyond the distinction between 2D versus 3D and static versus dynamic models. Table 1 has summarized the different methods used for modeling a continuum beam in soft robotic literature. Here, we focus on modeling methods for 1-dimensional (1D) continuum elements, (continuum rods), as the most studied continuum structure in soft robotics research. Later we explain how to generalize models for continuum rods to 2-dimensional (2D) (membrane) and 3-dimensional (3D) continuum structures.
Two methods are widely used to describe continuum rod kinematics; 1-I-i) continuous and 1-I-ii) discretized Kinematics. 1-I-i-a) Variable Curvature (VC-finite strains assumption), 1-I-i-b) truncated series shape functions, and 1-I-i-c) Constant Curvature (CC-as a subset of general shape function approach), which is probably the most simple and widely used assumption for soft manipulator modeling report instances of employing continuous kinematics.
Instances of using a discretized representation of such system kinematics are as follows. 1-I-ii-a) Employing series rigid-body kinematics, by simplifying a continuum rod as a hyper-redundant mechanism with finite but large enough number of segments, based on transformation matrices for consecutive but distinct rotational and translational joints, or methods based on 1-I-ii-b) Screw Theory and 1-I-ii-c) forward discretization of VC differential equations, where a skew-symmetric matrix of local curvatures/torsion vector is used to describe the local relative rotations along the backbone, are instances of using discretized representation of such system kinematics.
1-II-i) The Cosserat Rod method, 1-II-ii) Principle of Virtual Work (PVW), 1-II-iii) Beam Theory, and 1-II-iv) Lagrange Dynamics are used to derive the system governing equation (conservational law). The material constitutional law (material mechanics) is usually derived based on 1-IIIi) linear elasticity theory (Hooke's law), 1-III-ii) finite strain theory (considering large strain in hyperelastic materials, such as Neo-Hookean, Mooney-Rivlin, Gent, etc), or 1-IIIiii) by considering hyper-viscoelastic properties.
Any combination of the above choices results in a system of Ordinary (ODE) or Partial Differential Equations (PDE) to be solved numerically based on the system initial and boundary conditions. Using shape functions or discretized kinematics results in PDEs with decoupled spatial and time domains where direct solutions based on 2-I-i) analytical, if possible, or 2-I-ii) numerical forward integration steps in spatial and time domain can be used to solve the resulting initial value problems. If static solutions are sought, such systems turn into a Boundary Value Problem (BVP) where forward integration is valid when distributed loads, e.g.body weight, are neglected (Shiva et al. 2018).
Alternatively, indirect solutions can be sought. 2-II-i) Optimization-based methods, i.e. single shooting, multiple shooting, and concatenation methods, are suitable for BVPs resulting from static models with general loads, or for learning the coefficients of an approximate series solution or gains in a Neural Network (NN) model. 2-II-ii) Finite Element Methods (FEM) or similar segmentation methods are suitable if spatial (kinematics) discretization methods are used where, instead of a forward integration over the spatial domain, a system of nonlinear equations is formed with a large but sparse coefficient matrix. The system equilibrium point in static cases or at every time step of the dynamic simulation is found by calculating the pseudo-inverse of the coefficient matrix while satisfying all the geometrical, dynamical and optimal control constraints (Cianchetti and Menciassi 2017). While considering truncated series solutions as the system kinematics, 2-II-iii) reduced-order or Ritz method for solving a PDE problem, different choices of weighting functions can be used to improve the accuracy of the solution, e.g. in the case of Ritz-Galerkin methods. Finally, 2-II-iv) a combination of the above methods can be used, usually for solving PDEs resulting from complex geometries.
As an example of general practice in many commercially available FEM solvers, Tunay (2013) used a discrete Galerkin method, where weighted governing equations are used to construct the FEM solution for pneumatic actuators with general deformation. Sadati et al. (2017a) and Sadati et al. (2018b) used forward integration on the spatial domain for the PDEs resulting from employing reduced-order solutions for continuum manipulator kinematics, and then combined that with single shooting optimization method to find the system static solution under excessive external tip loads. Till et al. (2019) implicitly discretized the PDE problem time domain and solved the resulting BVP in arc length for each time step. Bieze et al. (2018) combined FEM and optimization methods to solve the closed-loop control problem of continuum manipulators. Gazzola et al. (2018) combined FEM with forward integration on time-domain in dynamic simulations. Thieffry et al. (2018b) constructed a reduced-order model based on dominant deformation modes that are found from multiple FEM based simulations of a system under different loading conditions. The coefficients of such solution were then optimized to solve for general cases. Cianchetti and Menciassi (2017), Bieze et al. (2018), andThieffry et al. (2018b) used a SOFA FEM modeling package for real-time dynamic simulation of soft structures.
In a comparative study with experimental results with a single module STIFF-FLOP appendage (Sadati et al. 2017b), we have recently shown the advantage of:
for dynamic analysis and traditional control design CC and modified CC (1-
for considering structural complexity and design parameter study
for accuracy in general cases -reduced-order series solutions (1-I-i-b & 2-II-i) for real-time performance, all based on 1-III-i if needed. We show that combining reduced-order kinematics, the Cosserat rod mechanics, numerical integration on spatial domain and optimization-based solution (1-
) produces most of the aforementioned advantages, i.e. accuracy, simple control design, real-time performance, considering structural complexity for a single STIFF-FLOP appendage in planar motion with excessive external load at the tip (Sadati et al. 2018b).
In this paper, we generalize our solution for multisegment arms in general 3D dynamic motion and compare the accuracy and numerical performance of the results with models with other assumptions. Additionally, the discretization method presented by Renda et al. (2018), which is based on Screw Theory and transformation matrices, is modified to use absolute (independent) states to achieve discretized models for multi-dimensional continuum geometries with a large number of states and significantly improved numerical efficiency.
A dynamic system with inertial, compliant and constraining elements can be expressed as a set of lumped (point) masses, usually assumed at the system elements' Center of Mass (COM) locations, with moments of inertia which are connected with spring/dampers and joints to the adjacent lumped masses. For a continuum system, where usually a system of differential equations describes the system mechanics, a differential format of the lumped-system approach can be employed. To this end, first, the free body diagram of the load balance in a single differential element is drawn, then the lumped-system equivalence of the system is assumed where the parameters are differential terms.
The main assumptions made in this paper about modeling a hybrid rigid-continuum system are as follows.
-The dynamic motion of a multi-link system is derived where external/input loads, geometrical constraints, rope elements, and soft impacts can be modeled.
Each element in the system can be assumed as a combination of separable inertial, linear elastic, viscous damping with power law, and external load elements, each with 3D elements.
--The system may have finite (only for the ROM elements) or infinite number of elements but must have a finite number of states, forming an ODE to be integrated numerically over time.
-External/input loads, elastic, damping, geometrical constraint, soft contact, and directional elements (such as string and membrane).
-1D continuum elements can be modeled as a finite number of interconnected Euler-Bernoulli elastic beam elements (discretization), or as continuous beam elements with predefined polynomial deformation shape functions (ROM).
-2D & 3D continuum elements can be modeled as wire meshes in which edges are 1D Euler-Bernoulli beams and connections are point masses.
-Hyperelasticity is not captured directly but can be added by updating an element stiffness matrix in an intermediate step during the numerical simulation. As a result, our modeling assumption has the following limitations that should be considered before putting into use.
-Material nonlinearity due to hyperelasticity is not directly addressed. Coulomb friction is not supported.
--Geometrical inputs, e.g. displacement and velocity, need to be transformed to geometrical constraints before being able to deal with.
-Hard contacts, i.e. impulse model, are not addressed although the provided methods can be adapted for such cases.
-Continuum systems with an infinite number of states, i.e. with PDE governing equations, are not addressed, hence the accuracy of the results is limited.
-2D & 3D continuum systems cannot be modeled as planar or 3D meshes where the use of shear rates instead of bending/twist improves the numerical performance of simulations.
-numerical stability and performance of the presented methods are noted but not thoroughly investigated.
-Large systems, especially with multi-material layers, are hard to implement and computationally expensive to simulate.
-Solutions for the system kinematic singular points, the sensitivity of the dynamic system to highly nonlinear material, or sudden transition in the system states are not looked into resulting in instability of numerical simulations in such cases. In the following sections, first, the continuum-body system kinematics using quaternions and the reduced-order method of using truncated polynomial series is described. Then, the TMT dynamics of such systems are discussed. Finally, we explain how to derive the lumped-system equivalence of Cosserat rod and reduced-order methods with relative and absolute (independent) states. The rigid-body kinematics based on quaternions are reviewed in Appendix 1 to provide a unified framework for the less experts in dynamic system modeling.
As an example, the different presented derivations are showcased to derive the governing equations for the planar motion (2-DOF: a bending about y-axis θ 1 and an elongation along the backbone l) of a beam with a single element (n d = 1 for the SRL, EBR, and EBA) or expressed with a first-order polynomial (n r = 1 for the ROM). We refer to this example as "our simple example" in the rest of this paper. Fig. 1 presents the parameters of our simple model. Here, x -z i and Q i are the local frame directors and quaternion representation of rotation with 0 subscript for the reference frame. ρ 1 is the position vector for the 1 st point mass. l 1 is the axial location of the 1 st point mass, l 0 is the beam initial length, 0 is the axial deformation for this section, and ξ 0 = l 1 -1 is the section axial strain. θ 1 is the bending angle of the 1 st frame, θ 0 is the frame initial bending angle, α 0 is the change in the bending angle for this section, and ζ 0 = θ 0 is the section bending. f l is a point load at axial location l and g is the gravity acceleration vector.
In the rest of this paper, parts of the derivations that are presented for the first time and considered novel are
Parameters for planar motion of a simple beam (in light gray) as our simple example in this paper. This system is used to showcase the application of the presented models.
More precisely, we have θ1, α0 2 , α0 = [0, α0 2 , 0], ξ0 = [0, 0, l 1 -1], ζ0 = [0, θ 1 , 0], Q0 = [1, 0, 0, 0], and Q1 = [cos(θ1/2), 0, sin(θ1/2), 0]. Here l0 is the initial length of the beam and θ0 = 0. placed in a box frame. Parts of the derivations are from the literature but presented here using our choice of parameters to provide a unified modeling framework. These parts are either adopted as an intermediate step of our final derivations, to be implemented as a part of our software package, or used only for comparison purpose. The rest of the derivations are related to the presented simple example system to showcase the implementations of our methods. The latter two cases are indicated in the text.
System kinematics describes the geometric relations between the system elements in terms of rotation and translation. As a result, equations describing the position vector and orientation of each point of the system (usually only COMs and joint axis) are derived.
The general transformation of a rigid body in 3D space can be described with a combination of a rotational and a translational transformation. The rotations/orientations can be described in terms of rotation matrices R [3×3] , quaternions Q [1×4] or exponential coordinates based on screw theory and Lie group notations. The translation is a 3D vector ρ t . Rotation matrices and Screw Theory result in the same final set of equations, despite their differences in notation. Non-unit quaternion and axis-angle representation of rotations define a rigid body orientation with four states, where unit quaternions, rotation matrices, and exponential coordinates (screw theory) requires three states. A non-unit quaternions representation of rotation, implemented in the T M T Dyn package, does not suffer from inherent singular points (known as Gimbal lock orientation where the axes of two of the three rotations are driven into a parallel configuration) associated with the transformation matrices and screw transformations. It also does not suffer from the singularity at π [rad] rotations associated with a unit quaternion. Also, Matlab is shown to be faster in optimizing the derived equation when quaternions are used. However, quaternion arithmetic requires more numerical operations, as described in Appendix 1. In the rest of this paper, the Jacobian of the transformation map to the system state space (T m ) is
where T ρm = ρ m,q is the Jacobian of the transformation that maps the link COM position vector between the Cartesian and the system state spaces, ρ m is the COM position vector with respect to reference frame, q is the system state (DOF) vector, T Qj = (2(Q -1 j ) × Q j,q ) [2:4,1:nq] is the Jacobian of the transformation that maps the link COM orientation quaternion to the system state space, and n q is the number of system states (DOFs). Here, × represents quaternion multiplication as in Rucker (2018), Q -1 represents quaternion inverse (equivalent to conjugate for unit quaternions), superscripts ( ˙) and (¨) are for the first and second temporal derivatives, and subscript ( , ) represents partial derivatives as X ,x = ∂X/∂x. All the vectors in this paper are with respect to and expressed in the system reference frame unless stated otherwise. See Appendix 1 for a complete explanation and derivation of rigid-body system kinematics based on quaternion representation of rotations.
Continuum Rod VC Kinematics-We use 1-dimensional (1D) continuum elements, i.e. continuum rods, as the basis of modeling continuum geometries in this work. A higher-dimensional geometry (e.g. a mesh geometry) can be constructed based on lumped masses (e.g. nodes of a mesh) interconnected by these 1D elements (e.g. edges of the mesh).
Using the Cosserat rod method, which considers all six translational (strains-ξ) and rotational (curvatures/torsionζ) differential states, is beneficial for 1D continuum elements. It is not the most efficient method for higher dimensional geometries, where the strain field is sufficient to
Discrete: calculate the distortions too; however, it is well suited for our main purpose here, which is deriving easy to interpret EOM of hybrid systems suitable for system dynamics analysis and controller design. Such models are not aimed at exact geometrical analysis and design, for which FEM methods are more suited. They rather are built upon simplifying assumptions for improved performance to investigate the large state space of such systems to clarify their underlying dynamics and/or control opportunities/challenges. We have recently used Variable Curvature (VC), based on rotation matrices, and Beam theory to investigate continuum manipulator mechanics in static and dynamic motions (Sadati et al. 2017b(Sadati et al. ,a, 2018b) ) and for a hybrid rigidcontinuum body system proposing of a highly-articulated inter-locking interface for stiffness control of a continuum appendage (Sadati et al. 2018c,a). We showed that employing a beam theory approach simplifies solving the BVP for static analysis of a continuum rod with numerical optimizationbased or reduced-order model methods. However, it is not a good candidate as a part of a unified framework for modeling hybrid systems and especially in a lumped system approach framework. VC kinematics and the Cosserat rod theory are used to model the static mechanics of continuum rods, based on rotation matrices (Gazzola et al. 2018) and quaternion (Trivedi et al. 2008;Burgner-Kahrs et al. 2015) representation of rotations.
We start with the variable curvature kinematics as in (Rucker 2018;Trivedi et al. 2008) The rod backbone curve spatial configuration (ρ s ) and 1 × 4 rotation quaternion unit vector (Q s ), expressed in inertial Cartesian coordinates ([ î, ĵ, k]), are derived according to VC
where * is the operator for quaternion rotation of a vector as in Rucker (2018). Here, the equations are derived in a local physical curvilinear coordinates [ d1 , d2 , d3 ], where s is the variable along the backbone, d3 is tangent to the backbone, and at the rod base we have [ d1 , d2 , d3 ](s = 0) = [ î, ĵ, k] (Fig. 2.a). The simple implementation of quaternion rotation is more computationally expensive than using rotation matrices R s as R s = Rζ × s and ρ s = R(ξ s + [0, 0, 1]), where × for a 3-element vector denotes the standard mapping from R 3 to so(3) (Burgner-Kahrs et al. 2015). However, quaternions are reported to be better in terms of numerical integration accuracy and preserving frame orthogonality and vector length (Trivedi et al. 2008;Rucker 2018).
For our simple example we have, Q 1s = [cos(θ 1s /2), 0, sin(θ 1s /2), 0],
ζ 0s = [θ 1s , 0, 0], and ξ 1s = [0, 0, l 1s -1], where subscript show that the variables are a function of curve unit length s. Trivedi et al. (2008). Then Eq. 2 reduces to Eq. 3. System states in this equation are continuous. To incorporate the VC kinematic of Eq. 2 in a hybrid modeling framework, three different discretization (the SRL, EBR, EBA) and a general reduced-order (ROM) methods are discussed in the following sections. Fig. 3 presents a simple graph abut the assumptions made for each method.
Discretized VC Kinematics-Discretized versions of Eq. 2 with rotation matrices are discussed in (Takano et al. 2017;Renda and Seneviratne 2018;Renda et al. 2018;Shiva et al. 2018). Shiva et al. used
) (i is the element numerator,), which is probably the simplest assumption, failing to conserve the principal properties of a rotation matrix. Renda et al. used the same method in the context of screw theory as R i+1 = R i e ζi∆s , where e ζ∆s = ζ∆s + I [3×3] (Renda and Seneviratne 2018;Renda et al. 2018). Takano et al. used
in the case of using quaternions representation of rotation. The above relation is used for the comparison purpose in our study. A similar representation is discussed in (Shiva et al. 2018), Appendix section, arguing that the order of the rotations is not important as long as small enough elements are considered along with the backbone (infinitesimal curvatures/torsion). They showed that using any of the above methods does not affect the accuracy of modeling a short appendage with beam theory, even for large deformations. The same can be done by quaternion transformations too.
In the case of our simple example, we have Q 0 = [1, 0, 0, 0], Q θ = [cos(θ 1 /2), 0, sin(θ 1 /2), 0, 0], and hence for the transformation pair (Ξ 1 = {Q 1 , ρ t1 }) we have
We can define the system states as q SRL = [ 01 , 03 , θ 1 ].
Here we use the first representation used by Shiva et al. and Renda et al., since it is easy to interpret its inverse as
∆s, which is necessary for modeling a continuum rod with absolute (independent) modeling states (a new contribution of this paper). Using quaternions and their properties, the final form of the discretized equations are
, for their inverse. Here, α = ζ∆s is the local bending/twist angle vector, and = ∆ρ = (ξ + [0, 0, 1])∆s is the deformed local position/translation vector. Notice that Q i is the absolute orientation quaternion at each point, but α describes the relative orientation of consecutive elements. Note that [1, α i /2] in Eq. 6 is a non-unit quaternion by definition but we assume it is a unit one as long as α i ≈ 0 holds, i.e. small bending/torsion angles.
In the case of our simple example for Eq. 6 we get (see Fig. 1)
where 01 = x 1 , 03 = ∆l = l 1 -l 0 and α 0 = ∆θ = θ 1θ 0 . We can define q EBR = [ 01 , 03 , α 02 ] as the system states if using this equation set. This results in a set of relative states with respect to the local (segment fixed) reference frame to be used with the SRL or EBR model. The same set of states is used for the SRL case too, but the kinematic relations that
define the segments' relative transformation as in Eq. 4. Note that Q 1 from Eq. 8 is an approximation for this term in Eq. 5 for small angles θ ≈ 0.
For their inverse from Eq. 7 we get
It is clear form Eq. 9 that this assumption is only valid if cos(θ 1 /2) -1 ≈ 0, e.g. small bending angles (θ 1 ≈ 0). Then we have α 02 ≈ θ 1 . More importantly, Eq. 8 does not predict any lateral movement of the tip. Also, any lateral movement of the beam tip in the case of Eq. 9 is solely resulted from shear ( 01 = x 0 ) and not the beam elongation ( 03 = z 0l 0 ). These cases are similar to an Euler-Bernoulli beam where the lateral and longitudinal deformations are decoupled. As a result, any realistic model for continuum rods should have more than one segment along the backbone. Despite this, we continue with our simple example that only serves to present the implementation of the presented methods in a simple yet generalizable example. From Eq. 9 we can define q EBA = [x 1 , z 1 , θ 1 ] as the system states. This results in a set of absolute states with respect to the reference frame to be used with the EBA model. Eq. 6 & 7 show that the deformation of a discretized element can be modeled as a quaternion transformation pair Ξ = {Q α , }, which is a 3D translational joint with state space and initial value 0 , followed by a 3D rotational joint with state α, initial value α 0 and quaternion representation of
Having an element's initial bending/twist angle (α 0 ) and a local translation vector ( 0 ), the local deformation of the discretized geometry (∆α = α -α 0 and ∆ = -0 ) can be calculated for deriving the element viscoelastic mechanical action due to system deformation.
Kinematics-Reduced-Order Models (ROM) for continuum rod kinematics is discussed by Godage et al. (Godage et al. 2011, 2016) based on the pressure chambers' length for a pneumatic soft manipulator. The presented solution is hard to interpret, results in complicated dynamic derivations, and the mechanical coupling between the actuation chambers' input pressure and length are not considered. The need to learn a large number of coefficients that should be learned through experimental trials is another drawback of such a method. In our previous work, we showed the advantageous numerical performance and accuracy of using a truncated Lagrange polynomial series passing through some arbitrary points along the backbone (Sadati et al. 2018b). The proposed solution is easy to interpret for shape estimation and controller design, since the used polynomial is constructed using Cartesian coordinates of physical points, and has a small number of states (6 for a short appendage consisting of Cartesian coordinates of 2 points at the appendage tip and mid-length). Both methods solve the singularity problem of using Constant Curvature and rotation matrix representations.
However, we used the CC assumption to compensate for the imaginary torsion of a Frenet-Serret frame and to find the physical torsion of the appendage cross-section based on the input chambers' pressure. Besides, the cross-section shear was neglected and a mean axial strain is assumed along the backbone. The kinematics was combined with Beam theory for static modeling and PVW for dynamic modeling, using Ritz and Ritz-Galerkin solutions. Duriez, Bieze, and Thieffry have recently generalized the same concept to modeling and control of complex continuum geometries by extracting the dominant deformation maps using the SOFA FEM package (Cianchetti and Menciassi 2017;Bieze et al. 2018;Thieffry et al. 2018b).
Here, we drop using Frenet-Serret frames and present a new general ROM approach to account for the cross-section local strains, as well as dealing with curvatures/torsion without any secondary assumptions, e.g. CC. Additionally, a simple polynomial is used instead of a Lagrange polynomial that results in simpler and faster derivation of the system kinematics. An inverse linear problem is solved to find the initial value of the polynomial coefficients based on the position and orientation of some nodes along the rod backbone. The final solution is more suitable for Cosserat rod and PVW methods. We assume that the manipulator geometry is defined by 6 truncated polynomial series of order n r + 1 for the position and n r for the orientation map, as
where Q ρ is the vector part of the quaternion representation
is the polynomial coefficient matrix which is considered as the system modeling states, S i = [s, s, s, 1, 1, 1]s i and S 0 = [0, 0, s, 0, 0, 0] are the shape function matrices that satisfies the rod base boundary conditions, i.e. being perpendicular to the base. Alternatively, the orientation can be expressed with a nonunit quaternion ( Q) to avoid singularities for curvatures with rotations larger than π. Then, we have [ρ, Q] in Eq. 10, and
Defining S = S 1:nr , we can rewrite the above equation in a vector-form as
(11) The state (coefficient matrix) initial values (C r0 ) are found based on position (ρ 0 ) and orientation (Q ρ0 ) of a few points along the manipulator backbone (s r ) by solving the following inverse problem,
The above inverse problem can be solved efficiently using Matlab inv function. ρ 0 and Q ρ0 can be simply measured from experimental observations using magnetic or visual trackers.
For the local strain ξ and curvatures/torsion ζ, from Eq. 2 and similar to the inverse map in Eq. 7, we obtain
ξ and ζ are used to calculate the mechanical action of the rod structural compliance. Alternatively, the system geometry can be described using only four polynomials, three for ρ and one for twist angle (θ 1), where the bending angles are found by compensating the Frenet-Serret frames non-physical orientation using θ 1. Such a solution results in a system with a smaller number of states, but more complex equations to handle, and inherent singular points.
In the case of our simple example, if we choose a unit quaternion representation of rotation and consider a 1 st order polynomial(n r = 1) for each of the three states ([x 1s , z 1s , Q ρ1 2s ], where
where q ROM = [C r1 , C r2 , C r3 ] is the system states. Knowing the initial configuration of the system tip position (s = l 0 ) to be [0, l 0 , 0], from Eq. 12 we get q 0 = C r 0 = [0, 0, 0]. For Q ρ1 s we get
Then, Eq. 16 can be derived from Eq. 13.
Note that ξ 01 has a value, meaning that this method can capture the material shear deformation along the local frame x-axis.
Dynamic modeling and controller design for hybrid systems are practiced, by having a kinematic model with a finite number of states. The discussed discretization method and reduced-order model are the key elements in modeling such systems. We use the TMT method as in Wisse and Linde (2007); Sadati et al. (2015) for deriving vector formalism of inertial terms in Lagrange EOM. The TMT method benefits from a smaller number of steps and simpler procedures compared to deriving EOM based on a system Lagrangian.
In the following section, the TMT method is adapted to derive a separate set of equations for each body/element in the system that speeds up the derivation and optimization of derived equations, as well as providing the possibility of implementing parallel numerical simulation methods for large systems. Then, the Principle of Virtual Work is used to derive the actions for the system springs, dampers, and input forces. In the rest of this paper, the link (rigid-body) inertial matrix (M ) is
where m is the body mass, I m is the body second moment of inertia [3 × 3] matrix, and I [n×n] is an [n × n] unity matrix. M for our simple example in the SRL, EBR, and EBA is
if the beam is assumed to be a uniform cylinder with radius r 1 , density σ 1 .
Element-wise TMT Method-From Newton's second laws of motion for the inertial terms (related to M ) and gravitational forces (related to g) of a rigid body EOM in general free motion, we have M χ = f χ , where χ = [ρ, Q ρ ], M as in Eq. 17, f χ are the forces in Cartesian coordinates, e.g. body weight f χ = f g = M [g, 0, 0, 0], and g = [0, 0, 9.81] [m/s 2 ] is the gravity vector. From mapping the system states from Cartesian coordinates to state space coordinates using the Jacobian transformation matrix, we have T m = χ ,q , where q is the vector of system states), and hence χ = T m q, χ = T m q + (T m q) ,q q.
Combining the above relations and transforming the equation to the state space (using T m ) we obtain, T m M T m q = T m (-M D m q + f g ), where D m = (T m q) ,q . While the above relation is valid for a large dynamic system, deriving the necessary vectors (M, T m , D m , f g ) for individual bodies in a system provides flexibility in dealing with bodies of a different type. So we have
where
, i is a general numerator for the body number in a system with n m bodies, q [nq×1] = q 1:nm is the vector of all the system states, and n q = Σ nm i=1 n qi is the number of states in the system. Note that the Jacobian matrices are calculated with respect to q (all the system states) and not q i . Eq. 19 can be solved for q to form an ODE problem, and then integrated over time using a numerical integration method, e.g. 4 th -order Runge-Kutta method implemented in Matlab software ode15s or ode113 function.
There is only one element in our simple example in the case of the SRL, EBR, and EBA. For T m as in Eq. 1 we have
Springs, Dampers, External & Input Loads-Following PVW, the mechanical action of such elements in the system state space (q) is wδq = f T f δq, where f is the force vector exerted by the element, and T f = χ f,q is the Jacobian that transforms the loads exerting location/orientation from Cartesian space to the system state space. The term w = T f f will be added to the right side of Eq. 19 as
where
represents other mechanical actions in the system due to viscoelastic elements, external, internal, and body forces, etc., and n fi is the dimension of vector f i .
For external loads (f l ) with action point ρ l , we have
If the elements are acting directly on the system states, i.e. elements parallel to the DOFs, T q = q is for the acting point, f lq is the DOF direct input,
is the force of a parallel spring, q 0 k is the resting value of the spring, and
is the force of a parallel viscous damper with power ν.
For simple axial (1D) elements, connecting two points of the system (ρ f1 &ρ f2 ), we have
where ρf = ∆ρ f = ρ f1 -ρ f2 is the line of action,
is the external force with value |f l | along the element unit direction vector where ρ = ρf / ρf ρ f ,
is the spring force vector, l 0 k is the spring resting length, and
is the viscous damping force vector. Note that the order of ρ a|b is not important in calculating ρf , as long as it remains consistent.
Finally, by monitoring the sign of the deformation of the element ∆l k = ρf ρ f -l 0 k , tension only (∆l k > 0, e.g. rope) and compression only (∆l k < 0, e.g. impact) elements can be modeled. The same procedure can be adapted for continuum elements, e.g. an Euler-Bernoulli beam or reduced-order model of a deforming beam, which is discussed in the following sections. The derivations for our simple example is provided in the next sections.
Continuum Body Dynamics-Eq. 21 can be easily adapted for a discretized continuum rod using the lumped mass method. The differential form of TMT terms are
where l i is the length of the i th the ROM element. For each of the ROM element in the system, the above spatial integrals can be handled with a numerical forward integration method, e.g. trapezoidal rule implemented in Matlab software trapz function, in each integrating time step of Eq. 19. The Cosserat rod method is used to derive the TMT method differential terms in Eq. 29 for the ROM elements. To this end, we start with the differential form and then differential lumped mass representation for Cosserat rod model.
For T m as in Eq. 1 we have T m SRL|EBR = T mEBA = I [2×2] . For our simple example in the case of the ROM, we have dm = m 1 /l 0 ds, dD m = 0, f g = [m 1 g/l 0 , 0]ds, and hence
Cosserat Rod Mechanics-The Cosserat rod theory presents the conservation law to balance the material local internal loads (f s , τ s ) due to local external (f l , τ l -external loads at the rod tip or body), internal (f u , τ u -due to internal actuation pressure or tendon tension) and body loads (f g , τ gdue to body weight or other uniform loads, e.g. a magnetic field). Different methods of deriving Cosserat rod mechanics are presented in the literature, based on distributed load balance in an infinitesimal element (Trivedi et al. 2008), and differentiation of the shear force balance on a long segment (applying variational calculus) (Rucker et al. 2010;Burgner-Kahrs et al. 2015;Renda et al. 2018). Governing equations for the above approaches can be derived using the Principle of Virtual Work too (Grazioso et al. 2018). Here, we follow (Trivedi et al. 2008) by considering the distributed load balance for free body diagram of a single differential element along the rod backbone (Fig. 2.a). For the load balance, expressed in the reference frame, we have
where ρ s and Q s are as in Eq. 2, I ms is the cross-section second moment of inertia, f g = σa s g, σ is the material density, and a s is the rod cross-section area. Note that the effect of f u and τ u are considered in the force ((Q s * f u ) ) and moment ((Q s * τ u ) ) balance equations after being transformed to the reference frame.
External local loads f ls , τ ls are local shear loads considered as local boundary conditions
where ∆ represents the variation between the integration boundaries. Knowing all the boundary conditions (ξ, ζ, f s , τ s ) at a point and all the external contact forces along the backbone f l , τ l , Eq. 33 can be integrated between the external contact points and free ends, using Eq. 32 to update f s , τ s at each external contact point (between two consecutive spatial integration steps). In a discretized or reduced-order modeling framework, Eq. 33 and 32 are handled alongside after direct (ROM) or indirect (discretization methods) spacial integration of Eq. 33. Hooke's law (linear stress-strain relation), expressed in reference frame, is the usual choice for the system constitutional law as
where k ξ and k ζ are diagonal stiffness matrices based on the rod material.
In the case of our simple example, we assume f u = [0, 0, f u3 ] and τ u = [0, τ u2 , 0] to be the force and moment from internal actuation, e.g. due to internal pressurized chambers or tendons.
is the tip external force vector which should be considered as a boundary condition of the form Eq. 32.
The above formulation results in a BVP problem which is hard to integrate into a unified modeling framework for hybrid systems. In the next sections, three methods for dynamic modeling of a continuum rod are discussed, using the discretized (Eq. 6 & 7) and reduced-order model (Eq. 10 & 13) kinematics, discussed earlier, and using TMT differential terms (Eq. 29), and PVW for compliance elements and loads explained above.
Discretized Continuum Dynamics with Relative States-Discretizing Eq. 33, a highly articulated system with length l, n d elements, and relative states (q = [ , α]) is formed with the kinematic relation expressed in Eq. 6. M and T m are found by substituting ρ i and Q i from Eq. 6 into Eq. 1 & 17, to find the TMT inertial terms as in Eq. 19. This forms the right hand of Eq. 33. The external loads are handled based on their exerting point, found from Eq. 6. In such systems, beam elasticity and damping, and the internal pressure/tendon tension acts parallel to the states q, so we set k|µ q = [k|µ , k|µ α ] and f lq = [f u , τ u ], and follow the relevant procedure for compliance elements and loads explained above. These form the left-hand side of Eq. 33.
Finally, the above terms are used alongside other terms in Eq. 21. The proposed procedure is easy to implement; however, the derived equations tend to be complex after having less than ten elements, which results in long segments, resulting in inaccurate results, slow derivation and simulation (Sadati et al. 2018d;Takano et al. 2017). The method is not suitable for large system models.
In the case of our simple example, from Eq. 8 and Hooke's Law, we have q SRL|EBR = [ 01 , 03 , α 02 ] = [x 1 , l 1 -l 0 , θ 1 ], for θ 0 = 0. Then, the stiffness coefficient matrices are
where J 1 = πr 4 1 /4. For the damping coefficient matrix, we assume µ , µ α to be constant values. The action for the compliant elements, i.e. continuum body, is directly related to the system states as w k|v = w k|vq , where
This is valid for both the SRL and EBA cases, while their difference is in the kinematic relations that define the relative transformation between the segments. Notice that, usually there is no direct input force on the direction of shear deformation (x-axis) when tendons or pressurized chambers parallel to the rod backbone are used as the actuation method.
Discretized Continuum Dynamics with Absolute States-To avoid complex derivations for a high number of elements, we may assume the discretized system states to be the lumped-masses' Cartesian positions and take the vector part of their unit quaternion orientation as q = [ρ, Q ρ ], or all the four elements of a non-unit quaternion as q = [ρ, Q]. Here we continue with the unit quaternion where Q 0 is derived based on Q ρ to form a unit quaternion.
The system kinematics is the same as q and increasing the number of elements does not increase the complexity of the derivations. M and T m are found by substituting ρ and Q into Eq. 1 & 17, to find the TMT inertial terms as in Eq. 29. The external loads become loads directly acting on system states f lq = f l . The inverse map presented in Eq. 7 is used to derive i and α i , based on which the beam elasticity (w k ) and damping (w v ), and the internal pressure/tendon tension (w u ) actions are calculated as
where χ b0 is the beam initial position vector and bending/twist angle that can be fund based on the system states' initial condition q 0 as χ b0 = χ b (q 0 ). The above terms are used alongside other terms in Eq. 21. The proposed method allows for handling a large number of elements. To the best of our knowledge, this is the first time that such a discretization method is used for modeling an actuated continuum manipulator, as well as its integration to a unified hybrid system modeling framework.
However, similar discretization methods are widely being used to solve hyperbolic PDEs, e.g. Eq. 33, numerically. Meeting Courant-Friedrichs-Lewy condition is necessary for converging the solution which usually results in systems with a very large number of elements and hence slow performance (Skeel and Berzins 1990). We do not analyze the convergence criteria in this paper, but a comparison with experimental results and the other presented methods in this paper are provided later.
For our simple example and from Eq. 9, we have
and
. k and µ are similar to the SRL and EBR models. Hence, Eq. 36 becomes
where w k|vq = 0.
Reduced-Order Model Dynamics-In the case of reducedorder model kinematics, the system spatial and temporal domains are decoupled. So we keep the differential form of Eq. 33, and perform a forward numerical integration over the ROM terms in each time step of the final system EOM numerical temporal integration.
Here the states are the elements of C r in Eq. 10 which gives 6 × n r states. The system kinematics is presented in Eq. 10 as χ s = [ρ s , Q ρs ], and T m is found by substituting ρ s and Q ρs into Eq. 1. for dM we have, dm = σads and I ms is found based on the second moment of inertia for planar objects with the shape of the rod cross-section. TMT differential terms are found from Eq. 29. The contact point kinematics of an external load at location s l along the backbone is found by substituting s = s l in Eq. 10.
Using the inverse map in Eq. 13 to find ξ, ζ, the differential form of Eq. 36 is used to find the action derivatives for viscoelastic structure and internal pressures/tendon tensions. This method does not suffer from discretization inaccuracy; however, the modeling accuracy depends on the order of the polynomial, while a higher number of terms does not necessarily improve the accuracy. Initial bent configurations, rods with initial arbitrary geometries are easy to handle, by choosing appropriate values for χ s0 .
For our simple example, the beam shear deformation along the local frame x-axis is derived in Eq. 16 too. Similar to the SLR, EBR, and EBA cases but in a differential form we get Eq. 38, where ξ, ζ, q ROM , and C r0 are as in Eq. 16.
Step for Higher Order Nonlinear Terms-High order nonlinear terms due to soft structure hyper viscoelastic behavior can be handled by considering nonlinear stiffness and damping coefficients. Such assumptions do not change the method of deriving such elements' actions as for compliance elements and loads explained before. However, the value of the nonlinear coefficients should be updated during the numerical simulation.
An intermediate numerical integration step is introduced to update the nonlinear coefficients based on the system current states. We have used this approach in our earlier work to account for the braid constraint of pneumatic actuators in continuum manipulators and the material hyperelastic deformation (Sadati et al. 2018b;Shiva et al. 2018). This intermediate step can be used for event handling and saturation constraints. For example, in the case of our simple example, this intermediate step calculates and updates the values for I|J 1 when the beam cross-section alters. Further details are provided in the later sections and in (Sadati et al. 2017a(Sadati et al. ,b, 2018b)).
Linear Modal Analysis-The concepts of exploiting natural dynamics of a continuum system for design, optimization and control of soft robots (Giorgio et al. 2017;Della Santina et al. 2018b;Sayahkarajy 2018) and continuum sensors (Trinh et al. 2018) have been investigated recently. Although it is not the main focus of this paper, the following section is dedicated to linear modal analysis of a hybrid system based on which a module is designed in the T M T Dyn software package. The theory presented below provides the necessary tool for such analysis on a hybrid system.
A linearized version of Eq. 21, without time-varying external forces, at a given point q 0 is im
where M = M,q(q0) , V = w v,q(q0) , and K = w k,q(q0) are linearized matrix coefficients (derivatives with respect to q evaluated at q 0 ) for inertia, stiffness, and viscous damping respectively, and A 0 is a constant term due to linearizion of d, gravitational forces, springs resting values, etc.. A 0 does not contribute to the system modal analysis.
Eq. 39 can be used for linear modal analysis of a system without damping or a proportionally damped system, subject to proper choice of damping coefficients ( V = A 1 M + A 2 K) and ν = 1, where A i is a general constant. The eigenvalue problem for Eq. 39 can be solved with the Matlab eig function to find the system undamped natural frequencies (ω) and matrix of mode shapes (Φ) as [Φ ω , ω] = eig(-M -1 K). Then the system modal damping ratio is η ω = (ωM ω ) -1 V ω , where M ω = Φ ω M Φ ω and V ω = Φ ω V Φ ω are the linearized inertial and viscous damping matrices in modal space.
The derived models for our simple example are already linear for the SRL, EBR, and EBA cases. From Eq. 18, 35 & 37 fro Eq. 39 we get
In the case of the ROM, the above terms are found by disregarding higher order terms with respect to the system states q ROM and then rearranging in a vector-form.
Constraints, Controller & Observation Design-Constraints can be modeled as soft constraints, e.g. elastic connections, and hard ones, e.g. geometric constraints. The former can be handled by adding viscoelastic elements to the system, while for the former the following closed form expression can be adapted from Eq. 21
where T c = χ c,q is the Jacobian transformation matrix for the constraint, χ c = 0 is the geometrical constraint relation, λ c is Lagrange multiplier, d c = (T c q) ,q q and u c is a control term that can be used to set a desired acceleration for the constraint term as u c = χc + C c . u c can be used to design a Jacobian based nonlinear controller by setting it to the desired acceleration control input for the constraint geometry (χ c ). In this case, T c is substituted with -T u and λ c with f u . C c can be a constant or a PID term, e.g. C c = C c P χ c + C c I t 0 χ c dt + C c D χc , to compensate for the numerical integration and/or inversion errors for a constraint, or to cancel out the control tracking error in a controller design. As a result, f u is derived alongside the rest of system states q in a way that satisfies the desired behavior that is expressed by the constraint terms T c , d c , and χc .
As an example, Eq. 41 can be adapted to derive an inverse-Jacobian based controller to control the beam tip configuration (position and orientation). Assuming the EBA model (q EBA ), for the tip configuration we have
as the desired constraint, where z c , θ c are the desired tip position and orientation angle. Hence
Note that we do not have any control on the tip lateral deformation in any of the SRL, EBR, or EBA modeling methods, since it is solely governed by the external load f l1 and not the system inputs. A model with at least two segments is needed to capture the effect of f u3 , τ u2 on this system. Adding the system inputs f u3 , τ u2 to the system states, and by substituting Eq. 20 and 37 in 41, we get to Eq. 43 for the control system EOM. Note that the system is not fully defined and a pseudo matrix inversion method is needed to solve the resulted system of equations. Eq. 43 is equivalent to a constrained system EOM but with an PD error compensation terms C c .
The discussed modeling frameworks for soft robots have enabled us to incorporate their states into the traditional modeling framework of rigid body dynamics.
In the following section, the T M T Dyn Matlab software package modules are explained to automate the derivation, simulation, and visualization of the EOM of a hybrid system. The T M T Dyn package flowchart are presented in Fig. 4. The package structure is presented in Table 10 and 11 of Appendix 3. The package consists of 5 main stages and 10 sub-stages, implemented in 13 modules each in a separate Matlab file. Some modules have sub-functions, embedded in the same file with local only access. The connection between the stages/modules are provided by three means: (i) passing a Matlab structure variable, named par and defined in the system.m file, between the modules, (ii) generating the necessary functions (Matlab, C or C-mex format) after completion of some stages and calling them in the next stages, and (iii) calling the main modules in a single file (system.m) that returns the resulting parameters to the Matlab workspace environment for later use.
The system geometry, in the form of a set of bodies that are interconnected with compliant elements and constraints, is described in the User Interface module (system.m) using an HLL. Next, the check.m module is called to check the defined parameters and assign default values if needed. The system is passed to the tmt eom derive.m module to derive the hybrid system TMT EOM. Depending on the method of choice for optimization and storing the derived EOM, Matlab, C-mex or C++ functions are generated by calling the save f unc.m and save mex.m modules.
The intermediate steps in the numerical simulation should be defined in int mid step.m by the user. Then, int eom mex.m is called to generate the code that is used for numerical static and dynamics simulations. equil.m, modal.m, and dyn sim.m are called to perform static equilibrium analysis, linear modal analysis, and numerical dynamic simulation. Finally, the results are passed to anim.m, providing a simple animation, and to post proc.m (provided by the user) for any intended post-processing of the results. The functionality and implementation of the HLL interface elements are explained in Appendix 3. Pseudocodes for each file are presented in Appendix 4.
Two cases of hybrid rigid-continuum systems are modeled and the results are verified in comparison to experimental results. The studied cases are as follows.
1. Dynamic motion of a single STIFF-FLOP continuum appendage in the presence of external loading (E1). 2. Dynamic deformation of a fabric sleeve worn on an elbow-like rigid-link pendulum (E2).
Different modeling assumptions for continuum rods are tested and compared for E1 in terms of modeling complexity, EOM derivation time, dynamic and static simulation time, and accuracy of the results in comparison with experiments. As a result, the most efficient and accurate method for modeling 1D continuum elements, i.e. continuum rods, is identified. The challenges in modeling a 2D (fabric) continuum medium are investigated in E2. The case studies show how the proposed unified modeling framework and software package help in deriving simple Lagrange EOM for hybrid systems. Further discussion is provided on the advantages of such models for controller and observer design tasks.
Where applicable, unit quaternions are used to model the systems since a rotation larger than 180 [deg] is not observed in our test cases. The mean (M) and Mean Standard Deviation (MSTD) values for the experimentally measured values are compared with the numerical simulation results to evaluate absolute (Err) and percentage error (%Err) for each case study. A Lenovo Yoga 3 Pro 1370 laptop computer (Intel Core M-5Y71 CPU, 2 × 1.2 -2.9GHz cores, 8 GB of RAM) was used to perform the analysis in this paper.
STIFF-FLOP Continuum Appendage-A STIFF-FLOP (STIFFness controllable Flexible and Learn-able manipulator for surgical OPerations) module (Fras et al. 2015) is a pneumatic continuum appendage (Fig. 5). Earlier tests show their high repeatability and negligible performance change due to aging and fatigue (Shiva et al. 2018).
The manipulator is made of silicone elastomer (Ecoflex 50 from Techsil) and selective actuation of the three braided pneumatic chamber pairs (6 chambers in total) via separate electronic proportional micro pressure regulators (Camozzi K8P) provides 3-DOF (one axial elongation and two sidebendings) of the appendage tip (Appendix 6 Fig. 14). The first chamber pair is placed along with the manipulator +yaxis with 120 [deg] offset from the other pairs (Fig. 5 The pressure regulators are connected to a compressor (BAMBI MD Range Model 150/500) and controlled via a data acquisition (DAQ) board (National Instruments Inc. NI-DAQmx USB-6411) and its control software. A control diagram for this study is presented in Appendix 6 Fig. 13. An ATI Mini40 6-axis force sensor is mounted at the appendage base and an ATI Nano17 6-axis force sensor is connected at the manipulator tip to measure the external loads at the appendage tip (f lt ). i=1 p i a p and the vector of bending/torsion is τ u = Σ 6 i=1 p i a p r Oi × [0, 0, 1], where r p and a p = πr 2 p are the pneumatic chamber inner radius and area.
.
Links are considered as EB beams with linear elasticity
).diag(J) which are diagonal stiffness matrices associated with strains ( ) and curvatures/torsion (α) respectively, in the di frame. Here, E and G ≈ E/3 are the material elasticity and shear modulus, J = π/4(r 4 a2 -r 4 a1 -6r 4 p ).[1, 1, 2] -a p diag(r O .r O ) is a 1 × 3 vector consisting of the cross-section second moments of areas. r O is a matrix of which rows are position vectors of the chambers in the manipulator cross-section plane as
where r o is the radial offset of the chambers from the center.
The finite length or infinitesimal length masses are assumed to be cylinders which are attached to the point of interests, e.g. joint in discretized methods, at their COM. As a result, we have a symmetric mass distribution in the body-fixed local frames. For the ROM the model we assume a hollow disk with differential moment of inertia
, where dm 2 = m 2 /l m2 ds. For discretized methods we have
). The effects of actuation chambers are neglected here.
The four modeling assumptions presented in Fig. 3 are implemented using the T M T Dyn package as in Tables 3 &4. The properties of these models and the abbreviation used to describe each are presented in Table 5.
Numerical Integration Mid-step Function-An intermediate step is applied in the numerical integration to account for the material hyper-elastic deformation (Gazzola et 2018), except for r p which remains constant due to the dense braid constraint (γ ≈ π/2) and in calculating f g which is independent of the cross-section deformation. | 3 | is the rod mean axial stretch which is l s=0 ξ 3 ds for the ROM model, where ξ i3 = sprdmp(i).dl. For the SRL and EBR models we have
, where is a part of system total states (Z). The same relation is valid for the EBA model, but i3 = sprdmp(i).dl is the linking EB beam lengths.
In addition, the force and moment due to input pressure f up , τ up needs to be updated during the simulation. The code to implement these considerations is provided in int mid step.m module. We introduce some symbolic variables in par.sym for | 3 |, f up , τ up and update their corresponding values in par.var in int mid step.m module.
Experimental Results & Discussion-Two sets of experiments were carried out with and without external load at the manipulator tip, and for different input pressures. Each experiment takes about 55 [s] and dynamic data for the actuator inputs, the manipulator tip position, orientation, and force were recorded. Sample experimental steps and recordings from the two experiments in comparison to simulation results from the EBR model with n d = 4 are presented in Fig. 6 &7.
Matlab ode15s was used to speed up the numerical integration in our simulation in this section. Considering hyper-elasticity effects, if possible, improves the simulation results accuracy by up to %6. This is in accordance with our previous observations in (Sadati et al. 2018b). We set µ ξ = 0.1, µ ζ = 1e -5 for the ROM and SRL, µ ξ = 100, µ ζ = 0.5 for the EBR, and µ ξ = 100, µ ζ = 0.01 for the EBA model. EBR model required higher damping coefficients to compensate the system sensitivity to rapid changes in the states in dynamic simulations.
Fig. 6-8 presents a comparison between computational performance (Fig. 6) and accuracy of these models in predicting experimental results with (Fig. 7) and without (Fig. 8) external loads. The ROM consumes the least memory and computer CPU time to derive the EOM, but the EBA is the best in terms of CPU time for equation optimization. The SRL is the worst in this regard. It takes hours to optimize the EOM for a system with more than three consecutive links. The change in system links, either to improve accuracy of a single rod model as presented here or in a system with multiple links, affects the EBA model CPU usage the least. For planar geometries (results are not provided here), a system with three times n r or n d of those presented here consumes the same memory and CPU time. As a result, the EBA is the best model for systems with a large number of bodies. The EBR and EBA have the best simulation time, as well as static and dynamic, performance. However, the EBR showed to be very sensitive to sudden changes in the input pressure and external force values. High viscous damping values were considered to prevent exponentially growing errors (numerical analysis diverge) in this case. As a result, the EBR simulation outputs were not reliable for fast dynamic motions. All the models show real-time performance (CPU time < 1 [s]) except the ROM which has the highest CPU time demand in simulations.
The ROM presents significantly lower errors, even for n r = 1, compared to other models, with 6 [mm] absolute error (Abs. Err.) & 9% normalized error (Norm. Err.) for static motion and 3.5 [mm] & 5.5% for dynamic motion in experiments without external loads (Fig. 7). These values are 12.5 [mm]& 19% for static and 9 [mm]& 14% for dynamic motion in presence of external tip loads (Fig. 8). The accuracy of other methods increases rapidly with a higher number of states and even slightly surpass the accuracy of the ROM for experimental cases without external tip load (Fig. 7). This accuracy increase is less noticeable and even reverses in some cases (static motion with external load, Fig. 8) for the ROM.
The ROM remains the most accurate model for experimental cases with external load. The EBA and the EBR methods show similar performance for static analysis of the experimental case without an external load while the EBA performs better for static analysis of cases with external load. However, the EBA method accuracy is much higher than the EBR method in dynamic simulations. The SRL method shows closer results to the EBA, but we did not report the results for a model with more than three segments since it takes hours to derive and optimize their EOM. Results for the EBR model were not reliable for dynamic simulations since high viscous damping values filter parts of the dynamic motion.
Our dynamic simulation results are more accurate than the static ones. However, we did not try to optimize the manipulator parameters or numerical analysis properties to find the best results, since this is not our main purpose in this study. It is possible to improve the accuracy by tuning these parameters or finding a way to effectively implement the hyper-elasticity assumption for all the modeling cases.
Table 8 presents the calibrated damping coefficient for each method and the cases where the cross-section deformation is neglected. Where the models are highly sensitive to the nonlinearities due to the change of crosssection, this effect is neglected to avoid large errors. This is observed for the SRL, EBR, and EBA methods in the experimental case with an external load where a sudden change in the external load is troublesome. The sensitivity increases as the number of system states increases. The EBA method is less sensitive to this increase. Furthermore, the dynamic simulations are slightly more robust and the modification is less needed, especially in the case of the EBA method. The ROM method showed less sensitivity to such nonlinearity and did not need any modification. On the other hand, the EBR model required higher damping coefficients to compensate the system sensitivity to rapid changes in the states in dynamic simulations. Considering these, the EBR method was more prone to be sensitive to numerical instability.
Overall, we observed better accuracy and less sensitivity to nonlinearities for the ROM, especially for a system with a small number of states and cases with external tip loads. The EBA method provides better numerical performance among all the methods, especially for large systems, and less sensitivity to the nonlinearities compared to the SRL and EBR methods. In the next section, the EBA is adapted to model a hybrid system with a 2D continuum membrane.
Pendulum with Fabric Sleeve Setup-A fabric sleeve, made of Jersey fabric, was cut and clamped on a rigid-link pendulum, made of ABS clear plastic (Fig. 9). The pendulum resembled a standard-sized human arm and is used to simulate the effect of clothing movement given wearer motion. The pendulum was fixed with a 1 DOF joint at the top and passively swings. The model was intended to capture the fabric dynamics due to the pendulum free motion. Three magnetic trackers were used to measure the link COM motion, and deformation of two points on the fabric (s 1 , s 2 ).
Capturing the dynamics of soft fabrics can be useful for research on wearable sensors and textile-embedded human motion analysis such as those used for computer animation and rehabilitation feedback Michael andHoward (2017, 2018). The fabric parameters are not known and identified to present the best correlation with the experiments. Table 10 presents the fabric and setup dimensions and simulation parameters.
The fabric can be modeled as a membrane which is a 2D tension-only continuum geometry that does not withstand bending or compression. This can be done by assuming the fabric as a net of equally distributed masses with connecting linear springs. Setting dir = 1 in the package express that the springs are tension-only elements. We have used a similar method to model a spider web with T M T Dyn recently (Sadati and Williams 2018). A drawback of such assumptions is that the model does not capture the crumbling Table 3. T M T Dyn package input for different models of a continuum rod, based on a pneumatically actuated STIFF-FLOP continuum appendage. Model labes are as in Fig. 3. Continues in Table 4... world.g = [0, 0, -g]; gravity body(1).m = m 2 /lm 2 ; ROM body body(1).I = I 2 ; joint(1).second = [1, lm 2 ] ; ROM joint joint(1).rom.order = nr; joint(1).tr.trans = [inf, inf, inf ]; joint(1).tr.rot = [0, inf, inf, inf ]; joint(1).dof(3).init = (1 : nr)lm 2 ; DOF properties joint(1).spring.coeff = [diag(Kv) , diag(Ku) ]; EB beam stiffness joint(1).damp.visc = [µ ξ , µ ζ ]; EB beam viscous damping joint(1).damp.pow = ν; viscous damping power law index joint(1).input = [fu p , τu p ]; pressure inputs body(2).m = m 3 ; tip force sensor body(2).I = I 3 ; body(2).l com = [0, 0, lm 3 /2]; joint(2).first = [1, lm 2 ] ; tip force sensor joint joint(2).second = 2; mesh = [ ]; no mesh body exload(1).exbody = 2; tip load exload(1).ftau = [f l t , τ l t ]; exload(1).tr(1).trans = [0, 0, lm 3 ]; ; joint(1).dof(i h ).damp.pow = ν; joint(1).dof(i h ).input = f χ(i h ) ; end body(2).m = m 3 ; tip force sensor body(2).I = I 3 ; body(2).l com = [0, 0, lm 3 /2]; joint(2).first = [1, n d ] ; tip force sensor joint joint(2).second = 2; mesh = [ ]; no mesh body exload(1).exbody = 2; tip load exload(1).ftau = [f l t , τ l t ]; exload(1).tr(1).trans = [0, 0, lm 3 ];
Discrete:
of the fabric between the mesh nodes. Also, the membrane assumption may not be accurate for thick fabrics, such as the Jersey fabric used in these experiments. This is more important when the fabric takes shapes such as hollow columns, e.g. a fabric sleeve, that is more resistive against buckling and bending. The fabric behaves like a shell, a 2D continuum geometry that withstands bending and compression too, in these scenarios. We used EB beam elements to resolve this issue without the need to increase the number of nodes or introducing diagonal connections, e.g. forming a tetrahedron mesh. This is a simplifying assumption that may be accurate enough to capture the underlying physics of a system with a thin membrane or shell geometries.
To model the system, we focus on the fabric model and import the link motion in the form of a constraint that follows an already recorded path (χ c = ρ com ) based on the experimental recording of ρ com . We use Eq. 41 to model the proposed constrained system. One geometric constraint is enough to fully define the link 1 DOF motion. Due to the constrained motion of the link, the values for m 1 & I m1 are not important for our analysis. Index 1 is used to define the link parameters and index 2 for the fabric ones. The fabric deforms when clamped on the link. The overall geometry of the clamped fabric is modeled with FreeCAD software as a wireframe sketch with a 3 × 5 grid of n d = 15 nodes and 22 edges as in Fig. 9.c. The CAD model is stored in IGES format to be imported into the model later. Here, nodes 12 & 13 are equivalent to s 1 & s 2 . The fabric is clamped at nodes 14 &
Discrete: The fabric mesh is modeled with lumped masses at the CAD-file wireframe nodes that are interconnected with EB beams. The system states are described with absolute Table 6. Experiments on free and forced deformations of a STIFF-FLOP continuum appendage. An external load is exerted by a force sensor at the tip of an actuated module while the appendage tip deformation is measured by a magnetic tracker. states of the masses. This is similar to the EBA model for continuum rods that was discussed in the previous section. The nodes are rigid lumped masses with an equally distributed mass of m 2 /n d . The relation for thin plates is used to derive their second moment of inertia as
Links are considered as EB ribbons with linear elasticity
).diag(J) as in the case of continuum rods in previous section. Here, J = [l 3 m2y l m2 b , l m2y l 3 m2 b , l 3 m2y l m2 b + l m2y l 3 m2 b ] is a 1 × 3 vector consisting of the cross-section second moments of area, where l m2 b = (l m2x + l m2z )/2 is the mean width of the ribbons in the x, z-axis directions.
The only parameter that needs updating during simulation is u c , the desired acceleration of the pendulum COM position, either in the xor z-axis direction. To map the motion of the nodes to the beams deformation map, one way is to define the xaxis vector that defines beam specific frames. This results in better accuracy but cumbersome calculations to derive and simulate the system EOM. An alternative assumption is to have the beams in the local body frames and initially deformed to reach the second connecting body. This is possible by setting init = nan which adopts the beam initial geometry ( 0 , α 0 ) to the system initial condition. We use the second method where 2D continuum geometries are modeled. This results in an almost 10 times decrease in the size of the file that stores the system EOM. The inputs for the T M T Dyn package to model this setup is as in Table 9.
Experimental Results & Discussion-Two sets of experiments were carried out to record the fabric deformation at points s 1 and s 2 . The static solution for the fabric deformation equilibrium point is challenging to find. The final result is sensitive to the initial guess for the system states that is passed to Matlab's f solve solver and unrealistic solutions may be found if the guessed states are not close enough to the equilibrium states. It is more accurate if a dynamic simulation is performed, by knowing the system initial conditions, to find the deformed state of the system. To test this, the pendulum was set free from an initial angle θ h0 ≈ 85 [deg] to swing under gravity and its joint internal viscous damping (µ h ). The fabric and pendulum 3D motions were recorded using NDI Aurora six DOF sub-millimeter accuracy magnetic trackers.
The pendulum Swung passively under gravity. The recorded data for the pendulum was imported as the desired rigid link trajectory in the model (Fig. 11.left), and the simulation results were compared with the experimental data for fabric deformation (Fig. 11.right). The simulation results show good accuracy in predicting node motion in the x-axis direction. The overall absolute error is ≈ 10 [cm] (≈ 40% normalized error based on l zs ).
Modeling such complex hybrid structures can have other benefits other than capturing the real dynamics with minimum error. For example, it is interesting to observe that the oscillation frequency of the node motion along the z-axis direction is twice the frequency of their oscillation in the xaxis direction. Also, the fabric is slightly deformed toward gravity (the decreasing mean value of sensor readings along the z-axis in Fig. 11) during the pendulum swing in both the simulation and experimental results. A simple model for the hybrid system helps in capturing such basic physics which are useful for system design, control, and observation.
Simple, reliable, and easy to use models for hybrid rigidcontinuum systems which are suitable for controller and observer design are highly sought in the field of soft robotics research. Current models are not suitable for control tasks Table 9. T M T Dyn package input for the sleeve fabric models clamped to a rigid-link pendulum. Model labes are as in Fig. 9.
world.g = [0, 0, -g]; gravity body(1).m = m 2 ; pendulum rigid link body(1).l com= [0, 0, -lm 1 ]; pendulum COM joint(1).second = 1; joint(1).tr.rot = [2, inf ]; 1 DOF rotation around y-axis joint(1).dof.init = θ h 0 ; pendulum initial angle joint(1).dof.damp.visc = µ h ; pendulum joint viscous damping Import mesh geometry: mesh.file name= cad.iges ; CAD-file name mesh.tol= 1e -3; geometry import tolerance mesh.tr.trans= [0, 0, -lc z ]; mesh geometry initial position/orientation mesh.tr.rot= [2, θ h 0 ]; mesh.body.m = m 2 /n d ; equally distributed fabric mess over the nodes mesh.body.I = I 2 ; Describing the mesh absolute DOF with mesh.joint(1): mesh.joint(1).tr.trans = [inf, inf, inf ]; masses absolute state as system DOFs mesh.joint(1).tr.rot = [0, inf, inf, inf ]; unit quaternion representation of orientation Describing the mesh EB beam connections with mesh.joint(2): mesh.joint(2).spring.coeff = [diag(K ), diag(Kα)]; linear elasticity of beams mesh.joint(2).spring.init = nan.ones(1, 6); beam initial state from system geometry mesh.joint(2).damp.visc = [µ , µα]; linear viscous damping Pendulum control constraint: joint(2).second = 1; joint(2).tr2nd.trans = [0, 0, -lcom 1 ]; joint(2).fixed = 1; constraining x-axis joint(2).control = uc; Fabric clamps: joint(3).first = 1; joint(3).second = 16; clamp at node 16 based on mesh file plot joint(3).tr.trans = [lc x , 0, -lc z ]; joint(3).fixed = ones(1, 6); fully constrained joint joint(4).first = 1; joint(4).second = 14; clamp at node 14 based on mesh file plot joint(4).tr.trans = [-lc x , 0, -lc z ]; joint(4).fixed = ones(1, 6);
fabric mesh (e.g. FEM) due to a large state space, reliable in general cases (e.g. machine learning methods), adequately accurate (e.g. lumped system method), or compatible with energy and Lagrangian based modeling and controller design techniques (the Cosserat beam method). The performance and accuracy of the newly introduced discretized (Renda et al. 2018) and reduced-order models (Sadati et al. 2018b;Thieffry et al. 2018b) that have promised to overcome these challenges are not yet thoroughly investigated. Additionally, there is no unified framework to implement these methods for hybrid systems with a combination of multi-dimensional rigid and continuum elements.
In this paper, we developed two new models for continuum rods and actuators: a general reduced-order model (ROM), and a discretized model with absolute states and Euler-Bernoulli beam segments (EBA). These models enable us to perform more accurate simulation of continuum manipulators, as well as modeling 2D and 3D continuum geometries, which has so far been missing in similar investigations (Renda et al. 2018). Furthermore, a new formulation is presented for a recently introduced discretized model by (Renda et al. 2018;Shiva et al. 2018) which is based on Euler-Bernoulli beam theory and relative states (EBR). These models are built in a new Matlab software package, called T M T Dyn (Sadati et al. 2015), to develop a modeling and simulation tool for hybrid rigid-continuum systems. The package performance is boosted using a new High-Level Language (HLL) text-based interface, a CADfile import module, automatic formation of the system EOM for different modeling and control tasks, implementing Cmex functionality for improved performance, and other modules for static and linear modal analysis of a hybrid system.
The package is used herein to compare and validate the aforementioned modeling methods in comparison to experimental results on the general motion of a STIFF-FLOP continuum appendage under external loads. We observed higher simulation accuracy (with as little as 8-14% normalized error) and numerical robustness (enabling consideration of material hyper-elasticity) of the ROM model, while the EBA was less computationally cumbersome to derive and simulate with near real-time performance. The EBR showed high sensitivity to sudden changes in the system actuation inputs and external loads and relatively higher computational cost both in derivation and simulation. The lumped system approach for modeling continuum rods as a hyper-redundant series-rigid-link system (SRL) is investigated as well. The SRL had the highest computational cost to derive and optimize the system EOM. This showed the importance of our investigation on discretized and reduced-order methods for modeling hybrid systems.
Finally, the package was successfully tested for modeling a system with 2D continuum geometries. A CAD-file import module was introduced to ease importing simple 2D and 3D continuum geometries in the form of wireframe objects stored in IGES CAD-file format. A fabric sleeve was fixed to a rigid link pendulum to study cloth dynamics. We observed reliable consistency between our simulation and experimental results. To the best of our knowledge, this study Sym.
Value Met. Sym. Value Met.
is the first that investigated low-cost analytical models for deformation of fabrics as a part of a hybrid system. Such models are important for wearable robotics research, and service robots which may deal with fabrics in their tasks.
As a part of this case study, we showcased how the package can be used to easily design a nonlinear controller for a hybrid system based on the system nonlinear map (system Jacobian).
As mentioned earlier, a fluent user interface needs to be as close as possible to human language. On the software package development side, we plan to develop a dynamic High-Level Language for our package to make it more accessible for less-expert researchers. Some preliminary results on this ongoing effort are presented in (Sadati et al. 2019b). Our long term goal is enabling robotic agents to flexibly produce abstract models of other robots and mechanisms in their environment. On the hybridsystem modeling side, we plan to test our package for observation and control of multi-arm continuum appendages and manipulators for medical applications, e.g. narrow port surgery and soft tissue probing. As a result, the package will be benchmarked for a hybrid system with 3D continuum geometry (the tissue sample).
A similar study to what presented here for 1D continuum geometries (rods) on the numerical performance and accuracy of different modeling assumptions is needed to compare different modeling methods for 2D and 3D continuum geometries in a hybrid system. The current package develops multiple copies of mostly similar codes in the case of the EBA for which most of the terms can be considered sparse vectors/matrices. We are looking into more efficient ways of employing the computational capacity and memory of the computing unit by only storing the none zero elements and collocating the system EOM in the run-time stage, using Matlab sparse function.
One of the interesting observations in our comparative study was the lower sensitivity of the EBA to nonlinearities introduced by considering the material hyperelasticity. We plan to look into why the same is not observed for the EBR, giving the same nature of both techniques. We also plan to further investigate the numerical stability, convergence, and computational efficiency of the introduced methods and the T M T Dyn software package in the near future. This is possible with performing standard numerical and modeling tests in comparison to analytical solutions or results from high fidelity modeling methods such as FEM. Our final goal is to use this package for hybrid force and position control and geometry, force and stiffness estimation in soft robotics research. Appendices
Unit quaternion representation of a rotation φ around a unit vector ρ
where Q 0 = cos(φ/2) and Q ρ = sin(φ/2)ρ, hence, basic rotation around frame x-, y-, or z-axis is simply setting ρ to î, ĵ, or k respectively. Alternatively, if a general 3D rotation/orientation is described with the three elements of Q ρ , we have
The method of defining a rotation with only three parameters does not suffer from singularity around 0 [deg] (which is associated with the singularity of CC and series rigid-link kinematics of continuum rods at straight configuration). Note that
Following any of the above formulations, by definition, can never produce a quaternion with both positive and negative scalar part and employ the whole space of quaternions, rather only the ones with positive or scalar parts. As a result, this method will suffer from inherent representational singularities, just like Euler angles do but for rotations approaching ±180 [deg] about any axis. One should use such method only for rotations ∈ (-π, π), not inclusive of ±π. An alternative is defining a non-unit quaternion rotation with all the four elements Q where (Rucker 2018) where the normalization step is implemented in the
is also possible to modify the existing operators, e.g. multiplication, to incorporate a normalizing term to deal with a non-unit quaternion, as in (Tunay 2013;Till and Rucker 2017;Rucker 2018). We use the former in this paper and implemented in the software package. While most of the derivations are based on a unit quaternion, non-unit representations are mentioned throughout the text where applicable.
Consecutive local rotations are handled by simple righthand multiplication R 1:n = Π n i=1 R i , or matrix multiplication of quaternions as
where Q × ( × for a quaternion vector) is equivalent to a matrix product as
and 1 : n means considering all the instances with index 1 to n. For rotating/transforming a vector ρ to ρ r with quaternions, we have ρ r = Rρ with rotation matrices and [0,
for a unit quaternion. We present the above quaternion manipulation by ρ r = Q * ρ in this paper for simplicity.
A translation (ρ t ) followed by a rotation (R or Q) can be represented by a 4 × 4 transformation matrix Γ = {R, ρ t } as Γ(1 : 3, 1 : 4) = [R, ρ t ], Γ(4, 1 : 4) = [0, 0, 0, 1], where consecutive transformations are handled by simple multiplication of Γ as Γ = Π n i=1 Γ i , and a vector transformation by ρ r = Γρ. We use a similar quaternion transformation pair as Ξ = {Q, ρ t }, where consecutive transformations are handled by Ξ 1:n = {Q 1:n , ρ t1:n }, in which ρ t1:n = ρ t1:n-1 + Q 1:n * ρ tn and Q 1:n = Q t 1:n-1 imesQ n . We present the above quaternion pair transformation as Ξ 1:n = Ξ 1:n-1 Ξ n in this paper for simplicity. All the vectors in this paper are with respect to and expressed in the system reference frame unless stated otherwise.
Each joint can be described as a set of transformations (Γ j or Ξ j ) and the system kinematics represents the joint axis transformations with respect to and expressed in the system reference (inertial) frame. Then, the linear velocity of the body COM, to which the joint (j) is immediately connected in the system main kinematic chain, is ρm = T ρm q, where T ρm = ρ m,q is the Jacobian of the transformation that maps the link COM position vector between the Cartesian and the system state spaces, ρ m = Ξ j * ρ com is the COM position vector, q is the system state (DOF) vector, ρ com is the COM local vector (with respect to and expressed in local frame), superscripts ( ˙) and (¨) are for the first and second temporal derivatives, and subscript ( , ) represents partial derivatives as X ,x = ∂X/∂x. The COM rotational velocity Ω m is derived with respect to the reference frame but expressed in the link local frame, since it is easier to calculate the link second moment of inertia with respect to the link local frame. We have, Ω m = T Qj q, where T Qj = (2(Q -1 j ) × Q j,q ) [2:4,1:nq] is the Jacobian of the transformation that maps the link COM orientation quaternion to the system state space, and n q is the number of system states (DOFs). Finally, the Jacobian of the transformation map to the system state space (T m ) become
Derivation steps of the TMT are compared with those of the Lagrange method for an unconstrained system with conservative forces is presented in Table 12. The fewer steps associated with the TMT method is the reason for implementing it in our TMTDyn package.
Potential H and kinematic K energies and forming Lagrangian L L,q, L, q dL q dt
Lagrange EOM terms Collecting coeff.s of q and form M
M , T = χ,q, and f for inertial, loads, and compliant elements
EOM closed vector-form M q = ...
Derivation steps of the TMT vs. the Lagrange method for an unconstrained system with conservative forces. Both methods start with the same system kinematics and result in the same closed form for the EOM. Terms in the same box can be handled in parallel. Lagrange method needs two extra steps. The step involving the collection of system states coefficients and forming M is very tedious and computationally time-consuming when using commercially available symbolic mathematical packages, e.g. provided by Matlab or Maple software.
The TMTDyn package elements are explained further below.
The interface of our package, or more broadly a HLL or DSL, need be fluent, i.e. it has a human languagelike flow. Properly named elements and utilizing a human language-like structure are the necessary considerations to this end. Use of method completion, where suggestions are provided to complete a syntax, or a dynamic language environment, where an agent can interpret the intention of the user from a less strict syntax structure, are the more advanced steps toward a fluent user interface Fowler and Parsons (2011). While method completion and dynamic language are not implemented in the current step of our package development, we have tried to use proper names and develop a sensible structure based on standard interpretation of kinematics of a multi-link mechanism, e.g. similar to DH (Denavit-Hartenberg) parameters Craig (2009). As a result, basic knowledge about the kinematics and dynamics of multi-link mechanisms helps to understand the implemented user-interface. The package control parameters, and the system properties (geometrical, inertial and stiffness) and inputs (actuator inputs, external loads) are described in system.m file, using a text-based High-Level Language. The HLL inputs are Matlab structure variables that store the package control parameters (par), properties of the simulation environment (world), and a dynamic system consisting of bodies (body), joints (joint), mesh elements from a CAD (Computer-Aided Design) file (mesh), and external loads (exload). We continue with our simple example to show how different elements of the package user-interface can be used. The fields, types, input options, default values, and suggested unit for these structure variables are explained in detail in Appendix 4 Table 13-16. Then different modules are called respectively based on the package control parameters par. These modules derive the EOM (tmt eom derive(...)), solve static equilibrium (equil(...)), perform linear modal analysis (modal(...)), simulate the system dynamics (dyn sim(...)), present and record simple animation of the system geometry and motion (anim(...)), and perform user-provided post processes on the results (post proc(...)). Pseudocode for system.m is presented in Appendix 5 Algorithm 1.
par contains the package control (par.anim, par.mov, par.derive, par.f un, par.mex, par.equil, par.modal, par.dyn, par.nint) and modeling (par.sym, par.var) parameters. The modeling parameters, e.g. system dimensions, can be set to numeric values or symbolic variable, stored in par.sym. It is possible to vary the symbolic values, stored in par.sym, after deriving the system EOM to provide flexibility. The numeric values remain constant speeding up the derivation, function optimization/generation and simulation processes. par.sym is considered as an input for the generated functions from derived EOM. par.var contains the numeric values for par.sym to be set after EOM derivation, e.g. for numerical simulation. For example, if we want to derive the EOM for our simple example which gives us flexibility to change the structural parameters (l 0 , m 1 , I 1 ), input (f u3 , τ u2 ), or control (u cz , u c θ ) values later, we can use the following syntax in Matlab language
Then par.var stores the numeric values for par.sym that can be changed before or during numerical simulation. For detailed explanation of package control parameters see Appendix 4 Table 13141516world.g stores the gravity vector as [g x , g y , g z ]. It is world.g = [0, 0, -9.81] in the case of our simple example. body, joint, and exload store the system structural parameters. Assume a system with n m different bodies or body sets (body), n j joints or joint sets (joint), and n l external loads (exload). A single joint or exload my describe a set of n d elements with similar properties, which we call a mesh. Series of links with relative DOFs, nodes with absolute DOFs, and any arbitrary interconnections between the nodes in a mesh are possible. A joint with rom field defines a continuum rod modeled with the ROM. A mesh geometry cannot be defined with a ROM joint. Our simple example has a body (n m = 1), a joint (n j = 1) in the case of the ROM, SRL, and EBR, or two joints (n j = 2) in the case of the EBA. The two joints for the EBA model is to define the body absolute DOF as well as the Euler-Bernoulli beam connecting it to the ground. We can assume an external load exload at the tip too (n l = 1).
The first joint that is connected to a body is a member of the Main Kinematic chain (MK). Any other joint connected to a body is a compliant element (if spring, damp, or input fields are defined) or a geometric constraint (if f ixed field is defined). A body that is connected to a mesh joint defines the inertial properties of the mesh elements. An alternative way to define a mesh body is importing a mesh CAD file. The mesh structure stores the CAD file name, a body field to be assigned to the mesh nodes, and two joint fields to define the bodies absolute DOFs (with respect to reference frame) and to assign to the mesh edges. To generate a mesh geometry with a ROM joint, a mesh with such joints can be defined to assign to a CAD file imported geometry. Brief descriptions for the fields in each structure variable are provided below. We use SI units throughout this paper.
body has the following fields. m is the body mass, I is the 3 × 3 inertia matrix, .l com is the COM position vector, and tip is the link tip position vector, both in the local frame. If a body describes a mesh with n d elements, each of the above fields can have different values for each mesh elements. A n d × 1 column vector is assigned to each parameter that each row stores the value for a mesh element with the corresponding row number. We follow the convention throughout this paper, that individual values for mesh elements are stored in a different column of a vector or matrix, except for I that becomes a 3 × 3 × n d cube variable. where eye(3) is a Matlab function that generates a 3 × 3 unit matrix. In case of the ROM continuum elements, all the fields present differential values with "per unit length" unit. n m is the total number of bodies defined, i m is their unique numerator, and n md = Σ nm i=1 n di is the total number of masses defined in the system. joint describes the system geometric connections and defines one of the following: i) system DOFs if being a member of MK, ii) a continuum ROM element if it contains rom filed and is in MK, iii) a compliant connection, e.g. spring, viscous damper, or actuator input, if not being a member of MK and contains any of spring, damp, or input fields, iv) a geometric constraint if the f ixed field is defined and joint is not in MK. It can define a mesh too if it is connected to multiple instances of a body, in any of the above cases, except for a ROM joint.
joint has the following fields. rom defines if the joint is a ROM link. rom has one filed order setting the ROM polynomial order (n r ) that should be defined (SBD). If the ROM with a 1 st -order polynomial is considered to model our simple example, we write joint(1).rom.order = 1;
(50) Otherwise, we do not need to define it. f irst and second are 1 × 2 or 1 × n d + 1 row vectors defining the body number and instance/axial location that the joint ends are connected to. The first element of the row defines the body number. Either of the ends should be connected to a single body. In the case of a 1 × 2 vector, i) the second element defines the ROM continuum beam length, if joint is a ROM (has rom field) in MK, ii) it is axial location of contact on the continuum beam, if joint is not in MK and the target body is a ROM, or iii) it is the instance of the body that the joint is connected to if joint is not in MK and the target body is not a ROM.
There is a single joint in our simple example that connects the beam to the ground. As in Appendix 4, we can leave joint(1).f irst undefined since its default value is zero, i.e. the joint first end is connected to the ground. For the joint other end, depending on the chosen modeling approach, we write joint(1).second = [1, l 0 ] ; (ROM) (51) joint(1).second = 1; (SRL, EBR, EBA)
In the case of a 1 × n d + 1 vector, i) the joint creates a mesh with n d elements if joint is in MK, and ii) it connects to instances of the body defined in 2 : n d + 1 elements if joint is not in MK. joint cannot be a ROM when defining a mesh. Leaving these fields completely empty (fields' default value) means connecting to the reference frame (ground), and if a scalar is assigned, the first instance of the body is used. f irst and second must have the same number of elements if both are defined and their assigned vectors have three or more elements.
Each joint can define n t number of consecutive transformation between its ends that are defined in the tr field. Each tr has a translation followed by a rotation pair defined in trans and rot fields with zero default values. trans is a 1 × 3 position vector in local frame. rot can be i) a principal axis and rotation set defined as elements of a 1 × 2 vector respectively (eg. [2, θ] defines a rotation of θ around local frame y-axis), ii) a 1 × 3 vector defining the axis of a unit quaternion (Q r ), iii) a non-unit quaternion 1 × 4 vector if the first element is inf , or iv) an angle-axis 1 × 4 vector, if the first element is any other number, with the first element as the angle (cannot handle DOF).
trans and rot can define a fixed transformation, if their elements are set to numeric values, or a free DOF if set to inf . The properties for each DOF are defined in a dof field with zero default value. n h is the number of DOF definitions and n q = Σ nj j=1 n hj n dj is the total number of states (generalized coordinates) in a system.
In the case of the SRL, EBR, and EBA there are two transformations (n t = 2) for the joint in our simple model, one the axial transformation of l 0 and one the shears and bending DOFs. However, the way the DOFs are defied is the key element that distinguishes between the four presented kinematic models. For the SRL case, we have transformation 0 along local z-axis, followed by a rotation α 02 around reference frame y-axis, joint(1).tr(1).trans = [0, 0, l 0 ];
(52) joint(1).tr(2).trans = [inf, 0, inf ]; joint(1).tr(2).rot = [2, inf ];
The EBR case is similar but with a translational strain 0 along local z-axis, followed by a bending strain α 02 around reference frame y-axis, joint(1).tr(1).trans = [0, 0, l 0 ];
(53) joint(1).tr(2).trans = [inf, 0, inf ]; joint(1).tr(2).rot = [0, inf, 0];
The EBA defines the general translation and rotation (expressed with quaternions) of a free mass in 2D space, joint(1).tr(1).trans = [0, 0, l 0 ];
(54) joint(1).tr(2).trans = [inf, 0, inf ]; joint(1).tr(2).rot = [0, 0, inf, 0]; Note that a unit quaternion is used, resulting in a system with three DOFs [x 1 , z 1 , Q r2 ]. Alternatively, a non-unit quaternion can be used as joint(1).tr(1).rot = [inf, 0, inf, 0]; resulting in a system with four DOFs
In the ROM case, join have only one transformation (n t = 1) that defines the free directions of the curvilinear frame joint(1).tr(1).trans = [inf, 0, inf ];
(55) joint(1).tr(1).rot = [0, 0, inf, 0];
Here, unit or non-unit quaternions should be used to define the curve bending/torsional DOF. Similar to the EBA case, the system may have three or four DOFs depending on our choice of quaternion representation of rotation.
dof has the following fields mostly with zero default values: initial value (init), initial values axial location for a ROM Joint (init s ), kinematic constraint between DOFs (equal2), elastic properties in spring (spring coefficient(spring.coef f ), spring initial value (spring.init), and initial compassion ratio (spring.compr)), active direction (dir), damping properties in damp (viscous damping coefficient (damp.visc) and damping power (damp.pow)), and direct actuator input (input).
init and init s are for the ROM elements, where the element base is considered fixed at the base of the joint local frame. Hence, there is no need to define an initial value at s = 0.
The spring/damping/input elements act in parallel to the DOF (having the same displacement). If the spring.init is set to nan, this value is assigned automatically based on the system initial configuration (system DOFs' initial value).
dir sets the active direction of these elements, not the DOF itself, and accepts a scalar value in [0, ±1].
equal2 handles kinematic constraints between different DOFs. Instead of defining a new state, the new DOF is set to equal with the referenced one in the equal2 value. However, a new set of spring, damping, and input elements are defined with their own dir property.
In the case of our simple example, all the models have three DOFs with an initial value of 0. It is also bidirectional without a compression ratio. So we can leave most of the dof elements undefined. The beam stiffness acts in parallel to these states for the SRL and EBR models, where for the elements in q with the order that they are appeared in tr(1), we have The above syntax changes slightly for the EBA model where there is no physical stiffness and damping that act parallel to the system states. So we leave all the elements of joint(1).dof (1) undefined.
For a ROM joint, init is a matrix with 3 rows, for the 3-value Cartesian location of a point at an axial location defined by init s, and an unlimited number of columns, the number of points along the axis which is equal to the number of init s elements. An inverse problem is solved in the check.m file to find proper initial values for the coefficients of an n r -order polynomial that passes through all these points. An initially curved beam can be defined as a ROM body by setting proper values for the init and init s fields.
In the case of our simple example with the ROM, we leave all the elements of joint(1).dof (1) undefined except for the initial condition in the case of the ROM. Any number of points along the beam can be used to define the initial shape of the polynomial curve. As a general rule, we may use (1 : n r )l m for initially straight curve as joint(1).dof (2).init = l 0 ;
(57)
where the 1 st and 3 rd DOFs remain zero along the curve initial shape. The code automatically adds the base point ρ (s=0) = [0, 0, 0] to the above set of points.
A joint defines a mesh by assigning matrices and vectors with n d rows to its fields, except for a ROM joint. All instances of a mesh joint have the same tr but can have different DOF properties. This is possible by setting the dof sub-fields to matrices with n d rows.
A joint in MK does not need any definition for a transformation regarding the connection to the second body frame. In the case of MK, the second body local frame is assumed to be attached to the joint itself. tr2nd defines such transformation with respect to the second body local frame if joint is not in MK. We name these "linking" joints. They have the same fields as tr but cannot have inf elements, i.e. define any new DOFs. A linking joint can have dir, spring, damp, and input fields as stated above. A linking joint can define an elastic Euler-Bernoulli beam by assigning sixcolumn vectors/matrices to these fields. The six columns correspond to [ξ, ζ] states of an Euler-Bernoulli beam as in Eq. 7.
dir field, here, can be a vector of values in [0, ±1] to set the active direction for each of the local frame directors individually, or be a scalar and assign the same feature to all the directions. The former feature can be used to model soft contacts with a surface where the contact point moves along the surface, e.g. bouncing of a ball on a surface while it can freely move parallel to the surface.
spring.init values can be set to nan for a linking EB beam. By default, such a beam is defined along the z-axis of the frame defined by tr. Alternatively, a new right-hand orthogonal frame is calculated using the aforementioned zaxis and the director axis defined in xaxis field. The matrices assigned to these fields can have n d rows, to define and/or assign different values to elements of a mesh joint. If a row vector is assigned for a mesh joint, same values are assigned to all the mesh elements.
For our simple example, in the case of the EBA, an Euler-Bernoulli beam connects m 1 (joint(2).second = 1) to the ground (there is no need to define joint(2).f irst element.). This is defined using a linking joint(2) as joint(2).second = 1;
(58)
joint(2).spring.init = [0, 0, nan, 0, 0, 0]; joint(2).damp.visc = [µ , 0, µ , 0, µα, 0]; joint(2).damp.pow = 2; joint(2).input = [0, 0, fu 3 , 0, τu 2 , 0];
Table 10. TMTDyn package structure. Underlined file names are provided or edited by the user (continued in Table 11).
Note that, although the model is a planar model, we defined the linking beam joint as a 3D link.
A ROM joint is similar to a linking joint in having the above fields and definition of an Euler-Bernoulli beam with a few differences. The assigned values for a ROM joint should have the same number of columns as the order of the ROM polynomial. Also, spring.init values should not be set to nan, neither a mesh can be defined for a ROM joint. In the ROM case for our simple example, the same fields as in Eq. 58 are defined for joint(1).
Finally, ref body, f ixed, and control define the reference frame in which the constraint is defined, constrained directions, and their desired control value. ref body is assigned similar to f irst or second fields of a body element. The constraint is defined in the local frame of the body which is defined by ref body. If ref body is not defined, the constrained are defined in the reference frame. When a ROM body is set as ref body, the axial position of the local frame should be provided or the tip location is assumed to be the default position of the body local frame (see Appendix 4).
f ixed accepts a boolean vector and a geometric constraint is defined for each direction that is set to 1. A 1 × 3 vector can be used to constrain local Cartesian directions ([x, y, z]) or a 1 × 6 vector to constrain the different states of an EB beam ([α, ]). This does not override the spring, damp, or input fields for the constrained directions. The dynamic actions associated with the compliant elements parallel to the constrained directions should result in zero. A spring and damper in parallel to these directions act as PD control terms to minimize the numerical simulation errors in satisfying the constraints.
control can be set to a symbolic variable and then updated during the simulations with the desired acceleration (since EOM is of 2 nd differential order) of the constrained geometry. Setting it to zero (its default value) fixes the constraint. It is useful for designing a Jacobian nonlinear controller. The dir field is ignored for constraints. Also f ixed field is ignored for a ROM joint.
For the constrained dynamic model presented in Eq. 43, we need to define a constraint joint (joint(3) for the ROM, or joint(2) for the SRL, EBR, and ROM) that defines a constraint between the beam tip (body(1)) and the ground and controls the tip vertical position z 1 and angle θ 1 as joint(3).second = 1;
(59) joint(3).f ixed = [0, 0, 1, 0, 1, 0];
where u c is the control term, e.g. PID term, to be updated in the int mid step.m file for each time step during the numerical simulation.
An alternative way to define a complex geometry is importing a CAD mesh file with extension "IGES" or "STL". Matlab igesToolBox provided by Ellenberg (2019) is used to import the files. The mesh structure has a f ile name to store the CAD file name, tol to set the points' overlap tolerance of the CAD file, tr to set initial position and orientation of the imported body, a body and two joint substructures with the same fields as stated before. A body is assigned to each node of the CAD file that at least two edges are attached to it. Nodes that are only attached to a single edge are assumed fixed connecting points with the ground. Each "body" has absolute DOFs, described in mesh.joint(1). A linking joint is assigned to each edge of the CAD file based on mesh.joint(2). The CAD-file coordinates transform based on tr filed and then used to define the initial condition of the imported elements. The links are numbered based on the order of extracted lines from the CAD file and overlapping of the nodes. Hence, it is hard to predict the assigned numbers to the links and masses. A plot is shown at the end of the importing process with labels showing the assigned mass numbers.
The imported bodies will be added after all the previously defined body instances, hence their labels may start with numbers greater than 1. The imported joints will be added before all the previously defined joints to satisfy our definition of MK joints. The import process is handled in the mesh import.m file. This module returns structure vectors for created body and joint instances. Only one mesh element can be defined (one CAD-file can be imported) at the moment. Pseudocode for this file is presented in Appendix 5 Algorithm 2.
External loads on the system are defined in exload. n l is the number of defined external loads. exload has exbody, ref body, and f tau fields. exbody defines the exerting body number and the body instance, in case of a mesh, or axial location of the external load, in case of a ROM body. exload can define a set of loads if a matrix with multiple rows is assigned to exbody, the same as defining a mesh with a joint that is explained earlier. The exerting point location in the local frame is defined by the trans and rot fields. The external force is defined as a 1 × 6 vector, in the form of [f, τ ] l , in the f tau field. f tau elements can be set to symbolic parameters and later updated in the int mid step.m file during the simulation steps. f tau is measured in the local frame of the body which is defined by ref body (similar to a constraint joint). If ref body is not defined, f tau is then measured in the reference frame. When a ROM body is set as ref body, the axial position of the local frame should be provided or the tip location is assumed to be the default position of the body local frame (see Appendix 4).
The external tip load in our simple example are exerted at body(1) but is defined in the reference frame. We may define it as exload(1).exbody = 1;
(60) exload(1).f tau = [f l1 , 0, f l2 , 0, 0, 0]; Most of these fields can be left undefined or empty. The check.m module assigns the default values, formats the input matrices, and performs some simple complication checks for the parameters. For a complete explanation of the above fields, and their acceptance and default values see Appendix 4. Pseudocode for this file is presented in Appendix 5 Algorithm 3. Examples of different hybrid systems are provided later to clarify the application of the HLL and text-based user interface. The inputs are then passed to the tmt eom derive.m module to derive the system TMT EOM. After calling the mesh import.m and check.m modules, the code proceeds to derive the system TMT EOM and generate optimized Matlab, C-mex, and/or C functions.
The system TMT EOM is derived in a set of Matlab structure variables. First the code iterates through all the joint vector elements (n j ), their instances (n d ) and the ROM order (n r ). This step identifies inf elements in trans and rot subfields, and generates the system states (q, q) and collects the elastic f kq , viscous damping (f vq ) and input (f uq ) actions in parallel to each DOF for the whole system. This step calculates each joint instance rotation (joint(i j ).Q(i d ).loc) and transformation (joint(i j ).T Q(i d ).loc) in the local frame too.
Then the code iterates through all body vector elements (n m ) and finds their MK joint. The joint type (the ROM or not) and the number of instances (n d ) define the type and number of mass elements in the system. The MK joints rotation (joint(i j ).Q(i d ).abs) and transformation (joint(i j ).T Q(i d ).abs) in the reference frame are calculated. The terms for each inertial element (M for Mi , T for T mi , Dd for D mi , and f g for f gm i ) are stored as fields of a mass structure vector. These terms are collected for the whole system in single variables (M for M , T for T m , Dd for D m , and f gv for f gm ). r jtips [ n md × 6] matrix stores the base and tip position vector of all the mass vector elements. rom.mass stores zeros for rigid links and length of the beam for the ROM rods.
Then the code iterates through all the joint vector elements (n j ) and their instances (n d ) to find the linking joints and constraints. The joint type (axial element, an EB beam or a ROM beam) defines the type of compliant elements in the system. The terms for each linking joint (T t for T j , kx to calculate f kj , vd to calculate f vj , dl for ρj , in for f uj , k mat for K j , v mat for V j , and dir for joint active direction) are stored as fields of the sprdmp structure vector. These terms are collected for the whole system in single variables (f j k for f kj , f j vd for f vj , f j in for f uj , f j sdi for f kj + f vj + f uj , f j k mat for K j , and f j vd mat for V j ). The r ks [ n m d × 6] matrix stores the base and tip position vector of all the sprdmp vector elements. rom.sprdmp stores zeros for rigid links and length of the beam for the ROM rods. The code checks if any of the link directions are fixed. lambda row vector and cnst structure vector store the terms for constrain (lambda for λ ci , T for T ci , D for D ci , and in for u c ). T cn for T c , Dcn for D c , and Ccn for u c are single variables that collect all the constraint terms in the system.
Finally, the code iterates through all the exload vector elements (n j ) and their instances (n d ). The terms for each external load term (T t for T li and f tau for f li ) are stored as fields of the loads structure vector. These terms are collected for the whole system in single variables (T tef for T l and f tau ef for f l ).
The structure vectors are used in a framework as in Eq. 21 while the single variables that collect the terms for the whole system form the EOM closed vector-form for the whole system. We use the structure vectors to have more flexibility in handling different link types and for faster derivation and optimization of the derived parameters. For a complete definition of the variable names refer to Appendix 5 Algorithm 4-6 and source codes available at (Sadati 2017).
Derived parameters are passed to save f unc.m and save mex.m modules for converting to Matlab, C, and C-mex functions and are stored as separate files in a folder named "code" in the active Matlab folder. Matlab matlabF unction, ccode, and codegen functions are used to perform this task. These functions optimize the derived equations by searching for and collecting repeated terms in the code. C-mex functions are only available for structure vectors. The generated Matlab functions are used to debug the model and code, by setting par.equil, par.dyn, and par.modal to 1. Once the model is debugged, C-mex functions can be used for increased numeric performance, by setting par.equil, par.dyn, and par.modal to 2. Pseudocode for these modules is presented in Appendix 5 Algorithm 7.
The system EOM are constructed using the generated functions in the previous step. The save eom mex.m module generates the EOM for dynamic and static analysis. The system states are named Z consisting of the geometric states q and Lagrange multiplier λ c vector and their temporal derivatives.
This module generates the EOM code as a string variable (string all) and stores it as a Matlab function file, using the Matlab f printf function. The module calls f j k, f j vd, and f j in to account for the compliance and inputs in parallel to the system states. Then it iterates through all the mass vector elements and generates a piece of code that calculates and sums up the summation in Eq. 19. The code performs a spatial integration if the mass is a ROM element based on its corresponding element in rom.mass. Then it generates the necessary code for the linking joints in sprdmp, external loads in loads, and constraints in the cnst structure vectors. Finally, the code to calculate Eq. 21 is added to string all. The differential terms are omitted for static equilibrium analysis. string all is written in EOM.m file for dynamic simulation and in EOM eq.m for static analysis.
A C-mex function is generated for each of the above Matlab functions. Pseudocode for this module is presented in Appendix 5 Algorithm 8. A similar module (save modal mex.m) generates necessary equations for linear modal analysis based on the system linearized EOM as in Eq. 39. The generated functions are stored in the EOM modal.m file. par.equil, par.modal, and par.dyn values control if an analysis is performed and which generated function is used.
Static equilibrium analysis is used to solve for the static model of a system or initial equilibrium point before a dynamic simulation. If par.equil is not zero, the equil.m module is used to solve the system static equilibrium. Matlab's f solve function is used to solve the static equilibrium problem. EOM eq.m is used if par.equil = 1, otherwise C-mex function is used. The dyn sim.m module is used to perform a dynamic simulation, if par.dyn is not zero. Matlab's ode15s or ode113 functions are used to perform a numerical temporal integration, using the EOM.m function if par.dyn = 1 or using the generated Cmex function otherwise.
The results (analysis time t and states Z) are passed to the anim.m module to generate and record a simple animation of the results. The simple plots show rigid links with continuous lines, the ROM continuum links with continuous line curves, and compliant joints (axial elements or beams) with dashed lines. par.anim and par.mov control if an animation is generated and/or recorded to a video format file. The post proc.m module is considered to perform userspecified post-processes on the generated data. This file does not have any code by default.
In "Case Studies & Validation" section, different examples of modeling multi-dimensional continuum and hybrid systems are provided. T M T Dyn enables us to provide comparisons between numerical performance and accuracy of using different modeling assumptions and complexity levels. This helps with deciding the most appropriate method for similar hybrid-system modeling tasks. : rotation expressed with a quaternion or Euler axis-angle rotation, set any of the 1 st -4 th element to inf form a DOF,... set the 1 st element to 0 form a unit quaternion (Q). If so, Q0 is set to QρQ ρ ,... set the 1 st element to inf form a non-unit quaternion ( Q). If so, Q is set to Q/ Q Q ,...
Prepared using sagej.cls
Journal Title XX(X)