A Muscle-Driven Mechanism for Locomotion of Snake-Robots

This paper presents the concept of muscle-driven locomotion for planar snake robots, which combines the advantages of both rigid and soft robotic approaches to enhance the performance of snake robot locomotion. For this purpose, two adjacent links are connected by a pair of pneumatic artificial muscles wherein an alternate actuation of these soft actuators causes a rotational motion at the connecting joints. The muscle-based actuated linkage mechanism, as a closed six-linkage mechanism, was designed and prototyped. The linear motion and force generation of the pneumatic artificial muscle was experimentally characterized using isotonic and isometric contraction experiments. A predictive model was developed based on the experimental data to describe the relationship between the force–length–pressure of the PAMs. Additionally, the muscle-driven mechanism was kinematically and dynamically characterized based on both theoretical and experimental studies. The experimental data generally agreed with our model’s results and the generated joint angle and torque were comparable to the current snake-like robots. A skx-link planar snake robot with five joints, five pairs of antagonistic muscles, and an associated pneumatic controller was prototyped and examined for simple movements on a straight-line. We demonstrated the muscle-driven locomotion of the six-link snake robot, and the results show the feasibility of using the proposed mechanism for future explorations of snake robot locomotion.


Introduction
The demand to develop autonomous robots with the ability to adapt and operate in unstructured and unknown environments has recently emerged, with the potential for a variety of applications, such as space exploration, search-and-rescue, and agriculture. Despite the extensive application of wheeled and tracked mobile robots as unmanned ground vehicles (UGVs) for these applications, the use of these robotic systems in unstructured environments and confined spaces is a challenge. Natural organisms and animals have provided inspiration for engineers in the design and development of robotic systems with increased adaptation and versatility to address the challenges faced by conventional robots. Limbless locomotion has proven effective, based on the abundance of natural snakes and other, similar species in a variety of environments, such as rivers, forests, and deserts. Limbless-based locomotion mechanisms have several advantages over conventional mobile robots (wheeled, tracked, and legged robots) including hyper-redundancy, high adaptability, intrinsic stability, and a small body cross-section [1]. Snake-like robots have been studied over the last 50 years and their effective locomotion in confined, narrow, and irregular environments has been demonstrated [2]. These robots consist of a serial kinematic chain of rigid-links connected through joints capable of bending in one or more planes that provide many degrees of freedom, which are suitable for the creation of complex motions [2]. These advantages make them suitable for use in applications in unstructured, unknown, and dynamic environments. Most early snake robots were designed with passive wheels for locomotion across flat surfaces (ACMIII [3] and Wheeko [4]). Recent snake robots were given more complex joints for obstacle-aided locomotion on uneven and cluttered environments (Kulko [5]), as well as for climbing trees and pipes (Uncle Sam [6], OmniTread OT-4 [7]). integrate the pneumatic artificial muscles (PAMs) into the structure of a snake-like robot to generate rotational motions at the connecting joints. PAMs (a.k.a., McKibben muscles [21]) consist of an inner expandable tube/bladder inside a braided mesh sleeve, which clamp together at the ends. The inflation of the inner tube causes the radial expansion while the geometry of the outer mesh sleeve acts as a scissor linkage and translates that radial expansion into a linear contraction (for a braided angle of less than ≈57 • ). These actuators are lightweight with a high power-to weight ratio and have load-length curves that are similar to biological muscles. Typically, PAMs contract 25% of their length; however, some new designs yield 40%. These muscles are attached on each side of the connecting links of a snake-robot to provide the rotational motion through their alternating contraction and extension, as shown in Figure 1. The muscle actuators are alternately activated by applying pressurized air to move the connecting links relative to one another. The present work combines the conventional and soft robotic approaches to enhance the performance of snake robots for navigation within long range, as well as unstructured and cluttered environments. The muscle-based actuated linkage mechanism was designed as a six-bar-linkage and a prototype was manufactured using rapid additive manufacturing. The pneumatic artificial muscles' linear motion and force generation was experimentally characterized using isotonic (concentric) and isometric contraction experiments. A predictive model was developed based on the experimental data to describe the relationship between the force-length-pressure of the PAMs. Additionally, the muscle-driven mechanism was kinematically and dynamically characterized based on both theoretical and experimental studies. The experimental data generally agreed with our model's results. A six-link planar snake robot with five joints, five pairs of antagonistic muscles, and an associated pneumatic controller was prototyped and examined for simple movements on a straight-line. We demonstrated the muscle-driven locomotion of the six-link snake robot and the results show the feasibility of using the mechanism in future explorations of snake robot locomotion.
The organization of this paper is as follows. In Section 2, the kinematics of the muscledriven mechanism and overall mapping between three configuration spaces (actuator, joint, and Cartesian spaces) will be discussed. Section 3 introduces the dynamics of a single module and a N-Link snake robot based on the articulated body algorithm technique. Section 4 presents the designing, CAD modeling, and prototyping of a muscle-driven mechanism, a snake robot assembly, and a pneumatic artificial muscle fabrication. In Section 5, the experimental testing and characterization of pneumatic artificial muscles, and the kinematics and dynamics of muscle-driven mechanisms will be discussed. Finally, Section 6 discusses the results of the experimental testing, force and motion characterization, and development and validation of the predictive model to describe force, motion, and joint torque.

Kinematics
The kinematics of robots is generally described through the mapping between three configuration spaces, actuator, joint, and Cartesian/operational spaces, as shown in Figure 2. Herein, the interrelationship between these three spaces for the muscle-driven snake robot was discussed to develop kinematic models, which will be used for the dynamic modeling and control of snake robots in the future.

Actuator Space to Joint Space
The kinematics of a muscle-driven, two-link mechanism were developed and studied. The muscles are attached to the links of the snake robot at two points on the right-side and two points on the left-side, as illustrated in Figure 3. To determine the length of the muscle, the position vectors of the connecting points, p R 2,i−1 and p R 1,i , on the right side are defined with respect to the base frame, and are obtained as follows where the variables are defined as follows: • h 1 = length of each link; • h 2 = distance from the joint to the attachment point; • w = half width of each link; • l = distance from w to muscle attachment point; The relationship between the joint angle and the length of the PAMs is given by where Although Equation (4), provides two solutions, based on the feasible range of mechanism motions and a given muscle length, a unique solution will be chosen.
Assuming the muscle length changes at a constant rate, i.e.,d i = 0, and the change in α i is negligible, the joint angular velocity and acceleration terms can be determined by differentiating Equation (3) with respect to time, presented as follows;

Joint Space to Cartesian Space
To describe the relationship between joint space variables and the position and orientation of the snake robot links, the forward kinematics of the planar snake robot, as shown in Figure 4, can be derived by writing the homogeneous transformation for the consecutive links. Note that the muscle-driven mechanism is not shown in Figure 4 for the sake of clarity. Additionally, the first joint is considered a floating planar joint (two translational and one rotational DOF) and the rest of the joints are one rotational DOF. The generalized coordinates, q = [q 1 , q 2 , · · · , q N+2 ] T , describe the snake robot motion in a 2D space where q 1 and q 2 are the coordinates of the origin of the frame {b 1 } attached to the first link of the snake robot and the rest are the absolute angle of the links with respect to the inertial frame {b 0 }, as shown in Figure 4. Starting from the first link, the homogeneous transformation is defined where q 1 , q 2 , and q 3 are the generalized coordinates describing the position (i.e., x and y components) and orientation of the first link, respectively, with respect to the fixed reference frame, {b 0 }. Afterwards, for every two adjacent links, where h 1 is the length of the link and φ i is the ith joint position determined by the difference between the two connecting links' orientation, described by the generalized coordinates q i+2 and q i+3 , which is defined as follows, where b 0 b i R and 0 r i can be used to determine the orientation and position, respectively. The linear and rotational velocities of each link were defined in the form of a Plücker coordinate using spatial velocity vector notion (i.e., 6D vectors in motion space domain, M 6 , which combines the linear and rotational aspects of rigid body motion) The angular velocity of each link can be written as where the generalized angular velocity of the links (q i+2 ) are related to the joint velocityφ i , expressed in Equations (5), as followṡ The spatial velocity of any given link,υ b i , can be determined in an outward recursive form asυ R is the rotation matrix defined in Equations (7) and (8), and b i−1 r j i−1 ,j i is the position vector describing the translation from the origin of the frame {b i−1 } to {b i } defined as, and [ ] × is the operator, which converts a 3D vector to a three-by-three skew-symmetric matrix.
In Equation (13),υ J i is the velocity across a joint, which is defined aŝ where the joint variables of the snake robot are defined aṡ denotes a subspace of motion at a given joint that describes the constraint motion on the connected body with N f ≤ 6 DOF. These S i matrix/vectors are defined as follows for the snake robot Note that the first motion subspace matrix S 1 is defined for a three DOF (two translations and one rotation), while the rest of the motion subspace vectors, S i {i = 2, 3, · · · , N}, are defined for one rotational DoF joint.

Dynamics
The dynamic model of a planar model of a muscle-driven snake-robot, as shown in Figure 4, with N links connecting N − 1 revolute joints along with the antagonist muscledriven mechanism at each joint, was derived. The dynamic model considers the snake robot without the common side-slip constraints (non-holonomic constraints). We assumed that anisotropic friction applied to each link of the snake robot with a friction larger in the normal direction than in the tangential direction along each link. Therefore, the robot has N + 2 degrees-of-freedom (DOF) and is an under-actuated dynamic system where the internal shape motion is no longer directly related to the overall displacement of the snake robot. The forward dynamics of the snake-robot, shown in Figure 4, were derived using the articulated-body algorithm [42]. The articulated-body algorithm calculates the dynamics of an N-link kinematic chain with O(N) arithmetic operations through a three-pass, outwardinward recursive procedure over the serial connection of the bodies.
Similar to the kinematics, spatial vector algebra (6D) was employed to derive the equations of motion. The spatial velocity,υ = [ω, υ C ] T ∈ M 6 and spatial force,f = [n C , f] T ∈ F 6 are combined the rotational and translational terms in one 6 × 1 vector representation where ω = [ω x ω y ω z ] T and υ C = [υ Cx υ Cy υ Cz ] T in the motion domain as well as corresponding terms, n C = [n Cx n Cy n Cz ] T and f = [ f x f y f z ] T in the force domain.

Dynamics of a Single Muscle-Driven Module
The free body diagram of two-adjacent links of the muscle-driven snake-robot is shown in Figure 5. The forces and moments acting on each body can be categorized into three groups of: (1) joint reaction forces, (2) friction forces, and (3) muscle actuation forces (on the left and right sides of each body). Based on the articulated body algorithm, we group the friction and muscle actuation forces into external forces and moments acting on a given body at the center of mass or about the center of mass.f J i andf J i+1 are the spatial forces transmitted across the proximal and distal joints associated with the link. The external friction force is applied with components in normal and tangential directions wioth respect to a given body, presented here in a spatial force vector form The friction between the snake robot and the ground was modeled as the Coulomb friction with anisotropic properties, µ t µ n .
where m i and g are the mass of the ith link and the gravitational acceleration, respectively. υ i and b i e k k = {1, 2} are the linear velocity of center of mass of each link and the unit vectors of the body frame {b i }.
As shown in Figure 5, the muscle actuation forces act on four distinct points of each link (except the tail and head links, which only have two points) and, for example, we can determine the muscle on the right bottom of the link {b i } in a spatial force form, as follows where b i e R i ∈ R 3 is the unit vector along the muscle direction and In Equation (22), f R i is the magnitude of the actuation force of the muscle on the right side of the link {b i }, which can be determined based on a predictive model that is developed and discussed in Section 6.2.1. Note that the formula for only one of out four muscle forces acting on the given body is presented here. The other forces can be determined with a similar structure, as presented in Equations (22)- (24).
Finally, the total external spatial force acting on a given link is given by The equations of motion for a single body presented in the body frame {b i } can be written asf where I i is the spatial inertial of a given body at the body's center of mass, and the spatial cross-product in force domain is defined as

Dynamics of N-Link Muscle-Driven Snake Robot
The articulated-body algorithm calculates the dynamics of an N-link kinematic chain with O(N) arithmetic operations through a three-pass, outward-inward recursive procedure over the serial connection of the bodies. The spatial equations of motion for a kinematic chain-based structure of the snake robot can be written as follows where S ⊥ i is the orthogonal complement of S i , including the constraint forces that impose these motion constrains at the joint.
Through the first pass, from the base (tail) to the last link (head), the velocities and bias terms p i are calculatedp the second pass, from (head) to (tail) computes the articulated-body inertia I A i and bias and the third pass from tail-to-head calculates the accelerations. Solving this for the body acceleration a i yieldsθ Solving Equations (27)- (33) in the described outward and inward passes will provide the equations of motion that govern the locomotion of the snake robot on a flat surface.

Conceptual Design and CAD Models
The concept for the muscle-driven mechanism was initially implemented for a two-link mechanism, and later modified to extend it to a multiple-links snake robot. Figure 6a shows the design of an individual link. Two pneumatic artificial muscles (PAMs) are attached to either side of the connecting links through the side features. The joints on the side of the body are designed to hold one end of the PAMs. The free rotation will help to keep the muscle straight under contraction to achieve the desired angular rotation. Cavities were made on the top and bottom side of each link to reduce the weight. The dimensions were determined based on the initial kinematic modeling to achieve a minimum 15 • rotation both clockwise and anti-clockwise. There are four rotating attachment points for the PAMs on each link to connect a series of these links with side muscles, as shown in Figure 6b. This allows for unlimited link extensions on either side of the link. The links are designed with tolerances that prevent them from interfering. Figure 6c demonstrates the greatest angle difference between the attachment points that the design allows. The rotating attachment points do not touch; therefore, there is no interference. Figure 7 demonstrates two links, attached to one another via a pin, to construct a muscle-driven module that allows for angular variation.

Prototyping
The designed parts were manufactured using a Form3 3D SLA printer (Formlabs, Somerville, MA, USA) and its standard Grey resin material. Each link was printed in one piece with tolerances that allow for full maneuverability of the rotating side attachments. Pneumatic artificial muscles were fabricated, as shown in Figure 8a-e, based on the methods developed in [43]. The muscle is fabricated from a thermoplastic elastomer (TPE) sheet and a braided mesh sleeve. A cardboard template was cut and used as a placeholder for the sizing and fabrication of the muscles and their cavity, Figure 8a. A piece of the thermoplastic sheet was wrapped around the template, as shown in Figure 8b. An impulse sealer was used to seal one edge of the muscle around the template. From one end, a piece of pneumatic tube was placed inside the muscle and sealed using hot glue, Figure 8c. The template was removed and the other end of the muscle was sealed using a stopper, shown in Figure 8d. A braided mesh was placed around the muscle and zip-ties were placed on both ends of the muscle (Figure 8e).

Control System: Pneumatics and Electronics
An Arduino-based electro-pneumatic system (shown in Figure 9a  The algorithm will actuate the muscles in a sequence that creates a sinusoidal wave in the snake mechanism. The first pass starts at the head and moves from head to tail to mimic the natural motion of a snake following the lateral undulation gait. Figure 10 shows an example of actuation sequence for a muscle-driven snake robot in two states (indicated by gray and black colors) on a Serpenoid curve path. The Serpenoid curve is defined as a curve, where its curvature changes as a sinusoidal function of the arc length. The actuation of either right-or left-side muscle depends on the location of the corresponding joint along the path, as shown in Figure 10. The sign of curvature (or joint angles) determines the actuation side (e.g., in this case, a positive curvature/angle requires left-side muscle actuation and vice versa). Since generating the lateral undulation gait requires a phase shift (time delay) between the joint angles, the corresponding actuation sequences are shown in the table, which describes the required spatiotemporal distribution of the muscle actuation along the length of the snake robot. The convention of labeling the actuators on the robot is shown in Figure 10. The percentage of muscle contraction depends on the joint angle values; the greater the joint angle, the greater the contraction. Equation (3) determines the relationship between the joint angles and the muscle length. The required actuation pressure to generate that level of contraction will be determined based on our experimental characterization of the actuator in Section 5.1.

Testing and Characterization
Experimental studies were carried out to characterize the kinematics of the proposed muscle-driven mechanism by relating the linear motion of the PAM to the joint angle of the snake robot. The obtained results were compared to the theoretical model developed for validation and modification of the model. Additionally, the quasi-static relationship between the force generation of the PAM as function of input pressure and kinematics of the muscle-driven mechanism were studied.

Kinematic Testing
A two-link module of the snake robot, Figure 11, was used to characterize the kinematics of the muscle-driven mechanism. The joint connecting the two-link module was fixed to a table and allowed for full rotation. The two links were free to move, as demonstrated in Figure 11. The muscles were pressurized on alternating sides of the two-link module, with pressures ranging from 0 to 200 kPa. The movement was recorded and analyzed using a video annotation tool (Kinovea). The relationship between the actuator space and the joint space was experimentally characterized by measuring three variables: (1) muscle length (d), (2) joint angles (φ), and (3) PAM attachment point angles (α), as shown in Figure 12.

Dynamic Testing
The pneumatic artificial muscles' linear motion and force generation were experimentally characterized using isotonic (concentric) and isometric contraction experiments. To analyze the force generated by the muscle, three variables were studied, including muscle length, input pressure, and output force. The relationship between these variables were two experiments performed and used to develop a predictive model that can be used for dynamic modeling and control of the muscle-driven snake robot. Figure 13 shows the experimental setup used to characterize the relationship between the pressure and force while the muscle length was kept constant with two ends of the muscle fixed. Input pressures ranging from 0 to 200 kPa were used to activate the muscles. One end of the muscle was fixed to a load cell (Mini Low Prf Load Cell MLP-50, Transducer Techniques, Temecula, CA, USA) with a capacity of 222.4 N (50 lbf) and 0.2% resolution. The force measurements from the load cell were recorded and used for analysis of the relationship between pressure and force.

Isotonic (Concentric) Contraction
A second experiment, an isotonic (concentric) contraction, was used to analyze the relationship between the muscle length variation and the input pressure at varying applied forces. Figure 14 shows the set-up for this experiment. One end of the muscle was fixed, while the other end was attached to a pulley system. The free end of the muscle was attached to a weight via a fishing line. The variation in the length of the muscle was measured at varying input pressures (0-200 kPa) for a constant amount of weight (force). The experiment was repeated for a range of weights between 0 and 500 g with an increment of 100 g.

Kinematic Characterization
The pneumatic artificial muscle was kinematically characterized to find the relationship between the pressure and length of the muscle. The muscle length contraction causes the change in the joint angle and attachment point angles on the snake links; see Figure 3.
The results for the muscle length at incremented pressures are shown in Figure 15. The input pressure was varied from 0 to 200 kPa. The length decreased from 60.0 mm to 46.7 mm. The coinciding joint angles that resulted from each length change are shown in Figure 16. Lastly, the results of PAM attachment point angles measured at each pressure are shown in Figure 17.  To find the joint angle variation in a single module of the snake robot in terms of input pressure, the pressure was varied between 0 and 200 kPa in the left and right muscles. A total of six muscle samples were used and tested. The mean value and standard deviation of the dataset were calculated and plotted. The joint angle, φ, varied between 0 and 30 • over the range of input pressure to one-side muscle. The achieved range of joint motion, ±30 • , is adequate for the overall maneuverability of the snake robot and is comparable to the other conventional snake robots (e.g., ±45 • reported in [2], ±34 • in [31], and ±25 • in [28]). The range of motion can be adjusted by using different dimensions/designs for the muscle-driven mechanism. Figure 17 shows the attachment point angle (also referred to as α) plotted over a change in pressure ranging from 0 to 200 kPa. The first point with applied pressure starts out positive and experiences a negative linear slope until the pressure reaches 93 kPa, where the trend changed to an oscillation at a constant value of the attachment point angle α = −10 • . Despite the observed variations in the value of the attachment point angle, the changes can be considered negligible. This assumption was used in our kinematic model to calculate the joint angular velocity and acceleration.  Figure 17 shows the angle measurements of two attachment points that are attached to the same muscle. As seen in the results, the two attachment point rotation variations follow the same trend and are closely in agreement. The obtained result validates the assumption that the adjacent attachment point angles mirror each other. This assumption was used in our kinematic model of the development of the muscle-driven mechanism.
To evaluate the derived kinematic model of the actuator to joint space in Equations (3) and (4), the results obtained from the experiments (shown in Figures 15 and 16) and the kinematic model were compared, as shown in Figure 18. The experimental data of muscle length versus joint angle was determined by eliminating the pressure between the two datasets. The experimentally obtained joint angle, φ, and attachment point angle, α, was fed into Equation (3) to calculate the predicted muscle length. The results show a similar trend with a shift between the experiment and theoretical model.
The difference (shift) between the theoretical and experimental results was caused by two effects: first, in our kinematic model, the muscle was assumed to be straight during contraction and the muscle's length is ideally measured along a straight line from end-point to end-point (Figure 11b), while a slight bending in the actual muscle after contraction was not considered in our theoretical model. Second, the radial expansion of the muscle leads to a physical contact with the connecting links (Figure 11a,c), which causes an unintended resistance to the rotation between these adjacent links; this effect was not considered in our theoretical model. Therefore, the proposed theoretical model can achieve a given angle of rotation with a smaller muscle contractions compared to the actual case. The deviation can be reduced by changing the design of the muscle-driven mechanism to eliminate the unintended resistance and reduce the muscle bending.

Dynamic Characterization
Two experiments were used to characterize the output forces of the pneumatic artificial muscles. The results from the experiment that measured the change in force at varying pressures with a constant muscle length (i.e., 60 mm) are shown in Figure 19. The applied pressures varies from 0 to 200 kPa and the measured resultant forces increase from 0 to 27.88 N. Since the weight of each muscle is, on average, 1.15 g (0.0113 N), the force-toweight ratio of these muscles is 2467.3. This is very significant in comparison to common electrical motors (actuators), which have been used in conventional snake robots [28]. Nine sets of data (three samples of pneumatic artificial muscles, and the experiment was repeated three time for each sample) were collected in this case. Figure 19 show the mean value of the measured force from the datasets, accompanied by the error bar equivalent to one standard deviation of the results at each point. Figure 20 shows the change in the length of the muscle when there is an applied force. The pressure ranges from 0 to 200 kPa. The weights we used to apply the force range from 0 to 500 g, which translates to 0-4.9 Newtons. As the applied force increases, the initial change in length decreases. Figure 21 shows the same set of data as the variation of muscle length versus variation of force for an applied pressure.

Force-Pressure-Length Model
In this section, a model was derived based on the experimental results obtained in the dynamic characterization to predict the force generated by the PAM in terms of the actuation pressure and current length of the muscle. This model will be used in our dynamic model to calculate the actuation force, Equation (22), exerted to the body of the snake robot. Since we do not consider any rate of change of these variables (i.e., F, P, and L), the muscle model will be quasi-static.
The observation of the data trend obtained for pressure versus length, shown in Figure 20, points to the use of a power function to fit the data, and we employed a logarithmic linear form of that function as follows log(L) = log(m) − n log(P) (34) note from Figure 22 that log(m) and n, the bias and slope, change with respect to the applied force, F, based on the experimental results of the concentric contraction results. The variations in these two terms as a function of force were examined, where the trend suggests a linear relationship between these terms and the force applied to a muscle. Therefore, the following linear curve-fitting models were used and fitted to the experimental data, as follows where a = 0.07, b = 4.2, c = 0.0113, and d = 0.07 are the coefficients in Equations (35) and (36), respectively, and their numerical values are obtained using a linear curve-fitting procedure in MATLAB. We used MATLAB curve-fitting functions and a linear model form with least-squares fitting type. Substituting Equations (35) and (36) into Equation (34) yields which describe the relationship between force, pressure, and muscle length. To validate the proposed model, Equation (37), we used the isometric test results, which were obtained for the constant muscle length of 60 mm, as shown in Figure 19. Note that L = 60 mm is the initial length of the muscle and remained constant while the pressure, P, varied based the experimental results. The obtained force results from the proposed model, Equation (37), were plotted against the experimental results of the isometric test, which are illustrated in Figure 23. Note that the experimental errors were determined as one standard deviation and shown as the error bars.
The results show a close agreement between the force predicted by the proposed model and the force results that were obtained experimentally. The normalized mean root squared error (NRMSE) between the model and experimental results was calculated 7.3%, which is comparable with the normalized average of errors between the three sets of experimentally measured data, 6.3%, as shown in Figure 23). Additionally, the calculated NRMSE is lower than the reported level of errors (e.g., 15-20% in [15] and 60-90% in [20]) in other PAM models. This indicates the validation of the proposed model and an accepted level of error to use the model in future dynamic and control algorithm development.

Dynamic Model-Based Torque Characterization
The moment of the muscle force about a given joint is determined based on Equation (22) to characterize and study the torque generated by the muscle-driven mechanism and its variation in terms of experimentally obtained joint angle, attachment point angle, and input pressure, as follows where A and B are the coefficients, which are a function of α i and d i , defined as follows where d i is the theoretical length of muscle defined in Equation (3). Additionally, F i is the muscle force, defined in Equation (37) and calculated by substituting input pressure and muscle length. The results are shown in Figure 24. The torque at the joint was calculated based on Equation (38), where the joint angle and attachment point angle are determined from the experimental results shown in Figures 16 and 17 at varying input pressures. The obtained results are shown in Figure 25. The torque generated by the muscle-driven mechanism varied in the range of 0.033-0.22 N·m.
The resulted torque appears to be significantly lower than the reported joint torques in other planar snake-like robots (e.g., 4.5 N·m, reported in [2] for Kulko and Wheeko) or a 3D snake robot with the bellow-driven 25 N·m reported for OmniTread in [28]). However, those higher reported torques are required based on a significantly heavier and more bulky design of their individual actuation module (approximately 1.0 kg and 3.4 kg for each segment of Kulko [2] and OmniTread [7,28], respectively). In comparison, each module of the muscle-driven module (i.e., two connecting links and two antagonistic muscles, shown in Figure 7) weighs about 90 g (0.09 kg), resulting in a torque-to-weight ratio of 4.8, while the torque-to-weight ratios of Kulko and OmniTread are 4.5 and 7.35, respectively. The results indicate that our proposed muscle-driven snake robot generates a comparable torque at each joint with respect to other snake-like robots. Additionally, it is known that the force generated by PAM is a function of its length and diameter [20]. Thus, in cases where we need a higher torque, it would be possible to improve the range of the force/torque by varying the size of the muscle. In future work, the actuator space and joint space variables could be explored as a part of design optimization process to determine their optimal values for improving the energy consumption efficiency and velocity of a muscle-driven snake robot locomotion. The variation in torque as a function of φ i and α i is shown in Figure 26. Note the increasing trend in the torque value at a higher pressure (>100 kPa), see Figure 25, due to a saturation of the joint angle (φ i ) and attachment point angle (α i ) around 30 • and −10 • , respectively. This saturation leads to a constant moment arm b i r R p 1,i , while the force magnitude is still increasing. This trend is illustrated in Figure 26.

Power Consumption Efficiency Analysis
The power consumption efficiency is defined based on the metric found in other works [2], as follows where, W T is the power generated in the form of a translational work to carry the total weight of the snake robot ∑ N i=1 m i g with an average forward velocity ofῡ, and W R is the input power in the form of a rotational work at the joint due to the average joint torquē τ and joint velocity ω of joints during a full locomotion cycle, T. Considering the current muscle-driven prototype with a total mass of 0.27 kg, 6 links, and 5 joints, the average joint torque 0.15 N·m, and average joint velocity of 2.27 rad/s (130 • /s), the power consumption efficiency in Equation (39) becomes where, for the similar forward velocityῡ range of 10-200 mm/s obtained by other snake robots, η varies between 0.016 and 0.32 for the muscle-driven snake robot, which is 1.5 times greater than Kulko snake robot reported in [2] with the same forward velocity.

A Muscle-Driven Six-Link Snake Robot Locomotion
To show the feasibility of the proposed muscle-driven mechanism, a six-link planar snake robot with 5 joints, 5 pairs of antagonistic muscles, and an associated pneumatic controller was prototyped and examined for simple movements on a straight-line. We demonstrate the muscle-driven locomotion of the six-link snake robot and the results show the feasibility of using the mechanism for future explorations of snake robot locomotion. Figure 27 shows a simple forward movement test with a six-link model. The algorithm that controlled the movement of the snake actuates the muscles in a sequence that creates a sinusoidal wave in the snake mechanism. The first pass starts at the tail and moves from tail to head to mimic the natural motion of a snake. The snake mechanism moved forward ≈100 mm in 36 s, as shown in Figure 27.

Conclusions
In this work, we present the design, development, kinematic and dynamic studies of a muscle-driven mechanism for snake-robot locomotion. Pneumatic artificial muscles (PAMs) were integrated into each side of the connecting links of a snake-robot to provide the rotational motion through their alternating contractions and extensions. The muscle actuators are alternately activated by applying pressurized air to move the connecting links relative to one another. A two-link mechanism was designed and prototyped using 3D-printing fabrication. PAMs were designed, characterized, and fabricated to generate the desired kinematic and dynamic responses. A pneumatic control system was developed, which can modulate the air pressure and control the airflow direction applied to PAMs. The assembled mechanism of a six-link snake robot and its control system were examined. The experimental results were used to characterize the muscle, kinematics and dynamics of the muscle-driven mechanism. The overall results show the feasibility of using the muscle-driven mechanism for snake robot locomotion.  Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: https://eltnmsu-my.sharepoint.com/:f:/g/personal/mahdihj_nmsu_edu/EksKO7 1ncIZDrIyWwJ-hnVQBdlGxCZE0IEOOYGqb490E0Q?e=Thqva3 (accessed on 18 November 2021).

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

Nomenclature
The following nomenclature are used in this manuscript: transpose of a given vector a R a vector describing a quantity on the right-side of a given link (body) a L a vector describing a quantity on the left-side of a given link (body) a J a vector describing a quantity at a given joint a b i a vector describing a quantity at a given body b î a b i a 6D spatial vector describing a quantity at a given body b i a x , a y , a z x,y, and z components of a given 3D vector