Dynamics and Computed-Muscle-Force Control of a Planar Muscle-Driven Snake Robot

: This paper presents the dynamic formulation of an artiﬁcial-muscle-driven and computed-muscle–force control for the planar locomotion of a snake robot. The snake robot uses a series of antagonistic pneumatic artiﬁcial muscles, assembled at the joints, to generate the locomotion. Kinematics of the artiﬁcial-muscle-driven robot in the joint and Cartesian spaces was derived with respect to the muscles’ motion. The Lagrangian mechanics was employed for the formulation of the dynamic model of the robot and deriving the equations of motion. A model-based computed-muscle-force control was designed to track the desired paths/trajectories in Cartesian space. The feedback linearization method based on a change of coordinate was utilized to determine an equivalent linear ( input-to-state ) system. Then, a full state feedback control law was designed, which satisﬁes the stability and tracking problems. The performance of the dynamic model and the controller were successfully demonstrated in simulation studies for tracking a circle-shape path and a square-shape path with a constant linear velocity while generating the lateral undulation gait. The results indicate a low magnitude of tracking errors where the controlled muscle force are bounded to the actual pneumatic artiﬁcial muscle’s limitations.


Introduction
Since the pioneering work by Professor Shigeo Hirose [1], who developed the first snake-like robot at Tokyo Institute of Technology, biological snakes and similar limbless animals have fascinated roboticists over half a century due to their vast abundance in a variety of environments (land, sea, and even gliding in the air [2]) as well as their versatile locomotion that adapts easily to unstructured and unknown environments [3,4]. Snakelike robots, in contrast to their biological counterparts, have shown limitations for agile mobility, long range navigation, and surveying the field with a large area mainly due to the low energy-consumption-efficiency and low power-to-weight ratio [2]. Most snake robots that have been designed and developed thus far use motorized joints for generating locomotion. The recent advancements in soft robotics help to enhance the capabilities of robotic systems by designing and developing snake robots [5][6][7][8][9][10] with compliant and structurally deformable bodies that generate more complex motion while reducing the cost, weight, and complexity of the mechanical structures [11]. That leads to potential solutions for addressing energy consumption while creating more complex motion.
An extensive body of research has been conducted on studying the dynamics, control, motion planning, and design of snake robots [1,[12][13][14][15][16]. Chirikjian and Burdick introduced an integral form of continuous backbone curves for solving the inverse kinematics of hyper-redundant robot locomotion with nonholonomic constraints following desired trajectories [12]. Wright et al. designed a modular snake robot capable of climbing inside and outside of pipes [13]. Ma [14] studied the creeping locomotion of the snake robot using the Serpentine curve in a constant steady-state velocity. A mathematical framework was established in [15] for the modeling, analysis, and synthesis of serpentine locomotion with a multi-link robotic snake without wheels that facilitates adapting to the environment at the expense of power efficiency. Tanaka and Matsuno proposed a control of snake robot that lift its parts of its body and dynamically switch lifting parts during locomotion which facilitates rapid and efficient movement across a sandy surface [16]. For a comprehensive review and list of research efforts on snake robots, the reader is referred to review papers and books [2,17]. Different dynamic modeling approaches have been utilized for modeling snake robots. Creep locomotion of snake robot on inclined surfaces without non-holonomic constrains was modeled [14]. Kane's dynamic formulation method [18][19][20] and the articulated body algorithm [21][22][23] were employed for computationally efficient modeling. Guo and Mahadaven [24] developed the governing equations for the planar lateral undulation of a continuum model of a snake robot, modeled as a thin filament, integrated with internal muscles that interacts with its surrounding environment. A geometric formulation was developed for studying a kinematic connection that decouples the low-dimensional net locomotion from the high-dimensional serpentine locomotion [25]. Modeling the transversal link displacement was introduced to derive a simpler model for a snake robot locomotion instead of the commonly used joint motion-based modeling approach [2,17].
From a control point of view, many research efforts were dedicated to the pattern gait control of snake robots [2,12,17,[26][27][28][29] with some work studied on the position/heading direction while considering side-slip constraints (nonholonomic) [30]. Ma used the Serpenoid curve equations and the body shape change with respect to the time variation of the Serpenoid curve for controlling the heading direction and the locomotion speed of the snake robot [26]. Xiao et al.,in [27], proposed a locomotive reduction approach based on the average body framework, which reduces the complexity of controlling a redundant snake robot to that of navigating a differential-drive vehicle. Hatton and Choset [28] introduced two gate generator approaches, called Annealed chain fitting and Keyframe wave extraction, that maps between the 3D shape of a snake robot's backbone curve and a set of joint angles and vice versa to address the gate design problems in snake robots. Ostrowski and Burdick [29] used geometric mechanics to describe the net motion of the snake robot as a function of variations in the body shape variables. Hicks and Ito [31] has studied the determination of optimal gaits for the control of position/heading direction for the locomotion of snake robots without side-slip constraints on a flat surface; however, the velocity control was not considered. A kinematic-based controller was studied that exploits the redundancy of a 3D snake robot for controlling its position and posture [32]. Slidingmode controller were developed to modulate the internal body motion and to generate the required serpentine locomotion by attenuating the lateral slip at each link to follow desired trajectories while compensate for uncertainties regarding the model and the environment parameters such as mass and friction coefficients [18,19]. Moreover, maneuvering control of land-based [33] and underwater [34] snake robots were investigated by enforcing a virtual holonomic constraint to produce the undulatory locomotion. In related works, a combined heading-velocity control based on undulatory gait parameters was studied for autonomous navigation and obstacle avoidance of a planar snake robot [21,22]. Liljeback et al. [17] introduced a combined control with three components, a body wave component (adjust force propulsion), a shape adapt component (adjust the body shape to the environment for obstacle-aided locomotion), and a heading direction control component for steering the snake robot. Chang et al. [35] developed a combined optimal trajectory synthesis strategy and a feedback controller for collision avoidance to adjust the snake robot's body shape for passing through narrow and confined spaces. A combined shape-velocity-heading strategy based on a variation of the lateral undulation gait parameters was studied for the maneuvering control of snake robot locomotion in navigating confined spaces with varying widths [23].
In recent years, there has been an increasing interest in studying bio-hybrid robots, which integrate living cells (e.g., excised or cultured muscle) coupled to a robot's structure to generate enhanced motion and force [36,37]. This promising concept relies on this fact that biological muscles, seamlessly integrated with sensorimotor controls, are ideal actuators with superior performance metrics. On the other hand, Pneumatic Artificial Muscles (PAMs), also known as the McKibben muscle [38], have shown to be promising in improving energy consumption efficiency compared to other actuators while generating a high powerto-weight ratio, a force-to-volume ratio, and a moderate to high order of force [38][39][40][41][42]. Similarities between the biomechanics of artificial muscles and natural muscles have made these actuators a feasible option for a variety of robotics applications including biomimemtic robots [39,43], human-robot interaction [44], assistive and rehabilitation robotic exoskeletons [41,[45][46][47][48], and robotic manipulators [42]. However, most of these applications are limited to immobile platforms [39,49]. The anatomy of natural snakes and particularly their musculoskeletal structure, which is composed of 130 to 435 vertebrae along with attached antagonist axial muscle configurations [14,50,51], can be considered as a series of rigid links hinged together, whilst axial musculature can be regarded as a series of elastic elements operating laterally to the hinges and between adjacent links. During lateral undulation locomotion, snakes alternately activate these antagonistic muscles on either side of their body to active flexural movement [24,[52][53][54]. The use of similar antagonistic mechanisms has been explored and studied for snake robots such as a tendon-driven antagonist mechanism for snake locomotion [51] and a snake-like surgical robot [55]. An "integrated joint actuator", a set of pneumatic bellows attached to links between a serpentine robot (OmniTread OT-4), which generate a 2 DOF motion, was introduced by Granosik and Borenstien [56]. Parallel elastic actuators (PEA), which combine an actuator (e.g., an electrical motor at the joint) and spring element in a parallel configuration, has been exploited for energy-efficient locomotion of snake robots [57,58].
In our previous work [59], an artificial-muscle-driven snake robot was introduced, which uses a series of antagonistic PAMs, assembled lateral to the joints, where alternate activation of the antagonistic muscles generate the internal body motion and consequently the external locomotion due interaction with its surrounding environment. Based on experimental and analytical studies and their characterizations, it was demonstrated that the kinematic and dynamic performances of the snake robots are comparable to the current state-of-the-art snake-like robots, and, in some aspects, even exceed them.
This work presents the mathematical formulation for the kinematics and dynamics of a planar artificial-muscle-driven snake robot with respect to the muscles' motion and force. The Lagrangian mechanics was employed to derive the equations of motion with muscle forces as inputs to the dynamical system. A model-based computed muscle-force-controller was designed to facilitate the tracking of the desired paths with different geometrical and topological features in Cartesian space while generating the lateral undulation gate in the body joint space. Feedback linearization was utilized based on a change of coordinate to derive an equivalent linear input-to-state system. Next, a full state-feedback control law was designed, which satisfies the asymptotic stability of the snake robot in tracking desired trajectories. Simulation studies were carried in MATLAB on a model of a six-link snake robot.
The organization of this paper is as follows. Section 2 presents the kinematics of the snake robot in terms of the muscle's motion, and overall mapping between the joint space and Cartesian space configurations will be discussed. Moreover, the dynamics of a N-Link snake robot based on the Lagrangian mechanics is derived. Additionally, the design of a computed-muscle-force controller is discussed where the model-based feedback linearization and full state feedback control laws are formulated. Finally, Section 3 presents and discusses the results obtained for two simulation case studies using a six-link snake robot model for tracking a circle-shape path and a square-shape path.

Kinematics of Artificial-Muscle-Driven Snake Robot
The kinematics of robots is generally described through the mapping between three configuration spaces, actuator, joint, and Cartesian/operational spaces. Herein, the in-terrelationship between these three spaces for the snake robot was discussed to develop kinematic models, which will be used for the dynamic modeling and control of snake robots in the future.

Joint Position-Muscle Vector Relationship
The relationship between the muscle movement and joint position of a snake robot was developed and studied. As illustrated in Figure 1a, 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, where the muscle can freely rotate about those two points which keeps it straight during the rotation of the joint. To describe the motion of each 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 body frame {i} at the joint {i − 1}, and are obtained as follows where h 1 , h 2 , w, and φ i are the length of each link, distance from the joint to the attachment point, half width of each link, and the joint angle, as shown in Figure 1a,b, respectively. The vector d i connecting these two points, corresponding to the muscle and called a muscle vector, can be obtained as Equation (3) can be written in the following form, where I ∈ R 2×2 is an identity matrix,R(φ i ) ∈ O(2) is an improper rotation matrix, i.e., det(R) = −1 in terms of the joint angle, and g is a constant vector, which includes geometric information of the link given bȳ we define D : R 2 → R 2 as an operator that maps a geometric vector g to a muscle vector d. Note thatR(φ i ) is symmetric (i.e.,R T =R) in addition to all other properties of special orthogonal matrices O(2). Therefore, the operator D has the following multiplicative properties which will be used through all kinematics' and dynamics' derivation. The length of the muscle vector d i can be obtained by where the relationship between the joint angle and the length of the PAMs is given by Solving Equation (8) for the joint position, φ i , yields Although Equation (9) provides two solutions, based on the feasible range of mechanism motions and a given muscle length, a unique solution will be chosen. The linear velocity and acceleration of the muscle vector d i can be obtained bẏ where S ∈ R 2×2 is a skew-symmetric matrix of the following form, andφ andφ are the joint velocity and acceleration, respectively. From Equation (10), it can be shown thaṫ where λ g = g T g > 0. By defining the bundled ∈ R 2(N−1) and its time-derivative˙d, the vector of joint angle velocitiesφ ∈ R N−1 can be obtained by, where J (φ) ∈ R (N−1)×2(N−1) given by,

Joint Space to Cartesian Space
To describe the position and orientation of the snake robot links in terms of the joint space variables d i , the forward kinematics of the planar snake robot, as shown in Figure 2, can be derived by writing the velocities for the consecutive links of an open kinematic chain. Note that the artificial-muscle-driven mechanism is not shown in Figure 2 for the sake of clarity. The first joint is modeled as a floating planar joint (two translational and one rotational DOF) and the rest of the joints are one rotational DOF. The snake robot motion in a 2D space is described by the coordinates of the origin of the frame {1} (x and y) attached to the first link of the snake robot and the absolute angle of the links θ i with respect to the reference frame {0}, as shown in Figure 2. The position of the head of the snake robot is defined by p = [x, y] T ∈ R 2 and the orientation of the links is selected as θ = [θ 1 , θ 2 , . . . , θ N ] T ∈ R N , a vector of joint angles is defined by φ = [φ 1 , φ 2 , . . . , φ N−1 ] T ∈ R N−1 . The link angles are related to the joint angles by the following relationship, The angular and linear velocities of each link are defined as follows, Using Equations (16)- (18), the vector of joint velocities ω ∈ R N−1 can be derived as, The linear velocities of the snake robot's links and muscles are where,ĝ = h 1 0 .
Substituting Equations (20) and (21) into Equation (23) yields, presently, by defining a bundle of velocity vectors asυ C ∈ R 2N×1 , the collection of the vectors, given in Equation (24), can be written in the following compact form,υ where A ∈ R 2N×2 and B ∈ R 2N×N are given by, By comparing Equation (22) to (23), a similar relationship can be derived forυ m , as follows,υ m = Aṗ + Gω (27) where G is given as,

Dynamics
The dynamic model of a planar model of the snake robot, as shown in Figure 2, with N links connecting N − 1, revolute joints along with the antagonist artificial-muscle-driven 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 2, were derived using the Lagrangian mechanics.

Lagrangian of the Artificial-Muscle-Driven Snake Robot
The generalized coordinates q = [q 1 , q 2 , . . . , q N ] T ∈ R N are selected where q 1 , q 2 , and q 3 are the position (i.e., x and y components) and orientation of the first link (θ 1 ), respectively, with respect to the fixed reference frame {0}. The rest of the generalized coordinates are defined as the function of the muscle vector as follows, q k+3 = 2g T d k k = 1, 2, . . . , N − 1 the time-derivative of q k+3 can be related to joint velocities by substitutinġ where presently, for the vector ofq = [q 4 ,q 5 , . . . ,q N+2 ] T we have,q = Λφ where Λ = diag{λ 1 , λ 1 , . . . , λ N−1 } is a diagonal square matrix of R N−1×N−1 with non-zero entries on the main diagonal, thus it is invertible and Λ −1 exists. The Lagrangian of the dynamical system can be written as follows [60], whereM ∈ R N×N andĨ C ∈ R N×N are given by, where M p ∈ R 2×2 and M d ∈ R N−2×N−2 are symmetric and positive-definite matrices given by, the diagonal and off-diagonal terms of M d are given as follows, Given Lagrangian L(q,q) of the robot, the equations of motion can be obtained by the Lagrange-Euler equation [60], where Q ∈ R N+2 is a vector of the generalized forces acting on the snake-robot defined in Section 2.2.2.

Generalized Forces
The free body diagram of two-adjacent links of the snake robot is shown in Figure 1c. 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). A vector of friction forces acting at the center of mass of each link can be written as follows, where the friction between the snake robot and the ground was modeled as the Coulomb friction with anisotropic properties, µ 1 µ 2 .
where m i and g are the mass of the i th link and the gravitational acceleration, respectively. As shown in Figure 2c, 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 {i} in a spatial force form, as follows where e m i ∈ R 3 is the unit vector along the muscle direction. In Equation (36), f m i is the magnitude of the actuation force of the muscle on the right side of the link {i}, which can be determined based on a predictive model that is developed and discussed in [59]. 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 (36) and (37).
The generalized forces Q j ∈ Q can be obtained by where the partial velocity terms forṗ = [q 1 ,q 2 ] T are defined, and finally, for the generalized muscle vector velocitiesq k+3 = 2g Tḋ k we have,

Equations of Motion
Plugging Equation (31) into (32) gives ∂L ∂q = 0 and, d dt (39) Note that M p is a constant matrix. The time-derivative of m θ 1 and M d can be derived as follows, and, where S ∈ R 2×2 is a unit skew-symmetric matrix, H j and H k are the jth and kth rows of H, respectively. Rearranging Equation (39) and substituting the generalized force from Equation (38) and the derivative terms expressed in Equations (40)-(42), the equations of motion can be written in the general form of a multibody dynamics where M(q) is a symmetric and positive-definite mass matrix, c(q,q) is the nonlinear acceleration terms including centrifugal and Coriolis terms, D(q) is the muscle force matrix, and J(q) T is an equivalent Jacobian matrix for the external forces acting on the snake robot given by,

Computed Muscle-Force Control
Now that the equations of motion of the snake robot is achieved, the goal here is to find a diffeomorphism coordinate transformation z = T(q) ∈ R N+2 and a feedback linearizing control lawf m = Q(z,ż) + R(z,ż)u such that it transforms the nonlinear equations of motion in Equation (43) into a input-to-state linearized system expressed by, Should this goal be achieved, a linear full state feedback control law, u, can be designed to facilitate the stability and tracking capability of the snake robot.
Given z = T(q), the first and second time-derivatives of z are obtained by, where, Solving Equation (46) forq and substituting into Equation (43) gives, where F † is a generalized inverse of F . Based on Equation (47), Q(z,ż) and R(z,ż) in the proposed feedback linearizable control law can be derived as, and, R(z,ż) = (F D) −1 F MF † (49) which transform Equation (47) into Equation (44). A linear state feedback control law is designed as follows, where z d ,ż d , andz d are the desired values of the new state variables z and its first two time-derivatives. The control gain matrices K d and K p are chosen in a way to exponentially guarantee the asymptotic stability of the linearized system, Equation (44), as well as the desired dynamical behaviors of the snake robot such as tracking the desired Cartesian path. The control objectives are set as follows: • Follow the desired path in Cartesian space based on the line-of-sight concept, which adjust the heading direction of the overall movement of the snake robot; • Generate the lateral undulation gate, which means the joints' motion of the snake robot must follow the desired sinusoidal function φ i = α sin(ω j t + (i − 1)β) + γ [1]; • Achieve the desired average linear velocity in the direction of the movement while keeping the velocity zero in the perpendicular direction (virtually generating the non-holonomic constraints).
Based on these control objectives, the new coordinate z is defined in terms of the generalized coordinates q, as follows; wherep is the geometric center andθ is the overall heading direction of the snake robot, defined as,p = p T are given in Equation (9) and H r is a row vector given as, The geometric center and the average link angle are employed to control the linear velocity and heading direction of the snake robot, respectively, in following the desired Cartesian paths, while the rest of the state variables are responsible to generate the lateral undulation gate (i.e., the Serpenoid curve [1]) at the joint space.

Results and Discussion
To examine the dynamic formulation and the performance of the designed controller, a series of simulation studies were conducted using a model of a six-link artificial-muscledriven snake robot, (N = 6) with five pairs of antagonistic muscles, as shown in Figure 2. Two case studies were considered, tracking a circle-shape and a square-shape path with constant forward velocity while following the lateral undulatory gate pattern. The simulation was developed in the MATLAB environment, which ode23 was used for simulation run. The mass, length, and width of each link are set to m = 0.09 kg, h 1 = 0.1 m, h 2 = 0.05 m, and 2w = 0.05 m, respectively. The normal and tangential friction coefficients are µ 2 = 0.5 and µ 1 = 0.05, respectively. The joint parameters and control gains are set as α = 30 • , ω 0 = 2π rad/s, β = 72 • , the control gains are k p = 1.25, and k d = 0.1. The initial joint angles and velocities are set to zero. The control goal is to show that the heading direction and forward velocity tracking errors are oscillating bounded around zero, since we expect the heading and forward velocity of the robot to display oscillating behavior during undulatory locomotion [2]. On the other hand, the joint angles φ i can follow the desired sinusoidal pattern with errors exponentially approaching zero.

Circle Path
A circle-shape path with a radius of 4.2 m and a center at (1, 4.5) m was chosen in the first case. The snake robot's head is initially located at (0, 0) with the joint angles are set to zero (i.e., a straight posture) and all linear and angular velocities are zeros. Here, we use the line-of-sight (LoS) guidance law [2] to define the heading reference angleθ d = tan −1 (p ∆ ), which will be used by the controller to adjust the snake robot's orientation (i.e., the heading direction). The look ahead distance ∆ is selected as 1.5 m, which is 2.5 times of the overall length of a snake of 0.6 m, as recommended in [61]. Figure 3 shows the snapshots fo the locomotion of the snake robot following the circle path in the counterclockwise (ccw) direction. Note that the robot was initially off the path and in a stationary condition. After reaching to the desired circular path, the snake robot's controller managed to stay on the path while finishing a full circle. The simulation time was adjusted to 100 s to ensure a full cycle of rotation, as shown in Figure 3.  The variation of the muscle's length over the period of 40 s of the simulated-time are illustrated in Figure 5. The relax length of each muscle is 100 mm. The PAMs are alternately activated, which contracts due to applying pressurized air, in a cyclic pattern according to the lateral undulatory gait pattern tracking the desired path in Cartesian space. The full length of simulation results are shown in Figure A1. The range of variation of the muscles' length are presented in the Table 1. The overall results are in the feasible range of contraction length of the PAMs. The computed-muscle-force controller determined the required force to be generated by the PAMs, which then was transformed into the controlled pressure based on the forcelength-pressure predictive model developed in our previous study [59]. The control inputs (i.e., the muscle forces) for the first 40 s of the simulated-time are shown in Figure 6. The results demonstrate a low level of force required overall to follow the circle-shape path. There is an extra effort at the beginning of the simulation due to move the snake robot from a stationary position and off the desired path. The full length of simulation results are shown in Figure A3.

Square Path
Performance of the controller is examined following a square-shape path with a constant forward velocity (0.275 m/s). The expectation was that the right-angle corners of the path would be a good challenging examination for the controller and locomotion of the snake robot. The square was considered with a side size of 8.8 m located as shown in Figure 7. The snapshots of locomotion of the muscle-driven snake robot following the path are demonstrated in Figure 7. Similar LoS guidance law with the same look ahead distance of 1.5 m was employed in this case and the controller was able to keep the snake robot on track for most part of the path. There is a deviation with a maximum error of about 0.9 m around the corner due to the gradual adjustment of the reference heading directionθ d determined by the LoS and the constant desired velocity. The error could be reduced by adjusting the look ahead distance in the LoS guidance law and using a variable/adaptive linear velocity reference input rather than a constant value. The simulation time was adjusted in order for the snake robot to complete a full cycle in the CCW direction.
The results of the actual linear forward velocity of the snake robot against the desired constant velocity are shown in Figure 8. Although the overall response of the system demonstrate that the controller was able to control the snake robot's forward velocity to exponentially approach the desired value in about 5 s in transition time and after then it stays bounded, oscillating around the desired value (with maximum and minimum values of 0.292 m/s and 0.268 m/s, respectively, due to an amplitude of 0.014 m/s with an average of approximately 0.28 m/s), but there are some small deviations (i.e., the average of the oscillating velocity is decreasing with the maximum error of 0.003 m/s) and rebounding behaviors corresponding to the turning motion of the snake robot passing through the corners. Similar to the circle-shape path, the actual forward velocity is oscillating bounded about the desired value. The variation of the muscle's length over the period of simulation time are illustrated in Figure 9. The relax length of each muscle is 100 mm. The PAMs are alternately activated, contracts due to applying pressurized air, in a cyclic pattern according to the lateral undulatory gait pattern and tracking the desired path in Cartesian space. The range of variation of the muscles' length are presented in the Table 2. However, the range of the length variation is slightly larger than the circle-shape path but the results are still in the feasible range of contraction length of the PAMs. A part of the results are magnified for better readability of the results. The control inputs (i.e., the muscle forces) are shown in Figure 10 determined by the computed-muscle-force controller. The results demonstrated a low level of force (0.2 N to 0.6 N) required overall to follow the square-shape path. There are instances corresponding to the robot's motion around the corners that require forces in the order of magnitude of 10 N.

Conclusions
This work presents the mathematical formulation for the kinematics and dynamics of a planar artificial-muscle-driven snake robot with respect to the muscles' motion and force. The Lagrangian mechanics was employed to derive the equations of motion with muscle forces as inputs to the dynamical system. A model-based computed muscle-force-controller was designed to facilitate the tracking of the desired paths with different geometrical and topological features in the Cartesian space while generating a lateral undulation gate in the body joint space. Feedback linearization was utilized based on a change of coordinate to derive an equivalent linear input-to-state system. Next, a full state-feedback control law was designed, which satisfies the asymptotic stability of the snake robot in tracking desired trajectories. Simulation studies were carried in MATLAB on a model of a six-link snake robot. The performance of the dynamic model and the controller were successfully demonstrated in simulation studies for tracking a circle-shape path and a square-shape path with a constant linear velocity, while generating the lateral undulation gait. The results indicate a low magnitude of tracking errors, where the controlled muscle force are bounded to the actual pneumatic artificial muscles' limitations. In future work, we will experimentally test and validate the computed muscle-force control implemented in our prototype.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

CCW
Counterclockwise DOF Degrees of freedom LoS Line of sight PAMs Pneumatic artificial muscles PEA Parallel elastic actuators 2D Two dimensional 3D Three dimensional Figure A1. Forward velocity of the snake robot following a circle path.
The variation of the muscle's length over the period of simulation time are illustrated in Figure A2. The relax length of each muscle is 100 mm. The PAMs are alternately activated, contracts due to applying pressurized air, in a cyclic pattern according to the lateral undulatory gait pattern and tracking the desired path in Cartesian space. The computed-muscle-force controller determined the required force to be generated by the PAMs which then was transformed into the controlled pressure based on the forcelength-pressure predictive model developed in our previous study [59]. The control inputs (i.e., the muscle forces) are shown in Figure A3. The results demonstrate a low level of force required overall to follow the circle-shape path. There is an extra effort at the beginning of the simulation due to move the snake robot from a stationary position and off the desired path. Figure A3. Time evolution of muscles' force of the snake robot following a circle path.