Multi-Phase Joint-Angle Trajectory Generation Inspired by Dog Motion for Control of Quadruped Robot

Quadruped robots are receiving great attention as a new means of transportation for various purposes, such as military, welfare, and rehabilitation systems. The use of four legs enables a robustly stable gait; compared to the humanoid robots, the quadruped robots are particularly advantageous in improving the locomotion speed, the maximum payload, and the robustness toward disturbances. However, the more demanding conditions robots are exposed to, the more challenging the trajectory generation of robotic legs becomes. Although various trajectory generation methods (e.x., central pattern generator, finite states machine) have been developed for this purpose, these methods have limited degrees of freedom with respect to the gait transition. The conventional methods do not consider the transition of the gait phase (i.e., walk, amble, trot, canter, and gallop) or use a pre-determined fixed gait phase. Additionally, some research teams have developed locomotion algorithms that take into account the transition of the gait phase. Still, the transition of the gait phase is limited (mostly from walking to trot), and the transition according to gait speed is not considered. In this paper, a multi-phase joint-angle trajectory generation algorithm is proposed for the quadruped robot. The joint-angles of an animal are expressed as a cyclic basis function, and an input to the basis function is manipulated to realize the joint-angle trajectories in multiple gait phases as desired. To control the desired input of a cyclic basis function, a synchronization function is formulated, by which the motions of legs are designed to have proper ground contact sequences with each other. In the gait of animals, each gait phase is optimal for a certain speed, and thus transition of the gait phases is necessary for effective increase or decrease in the locomotion speed. The classification of the gait phases, however, is discrete, and thus the resultant joint-angle trajectories may be discontinuous due to the transition. For the smooth and continuous transition of gait phases, fuzzy logic is utilized in the proposed algorithm. The proposed methods are verified by simulation studies.


Introduction
As robotic technologies, including manufacturing, control, and system integration technologies, have been improved dramatically, many sophisticated and intelligent robotic systems have been developed in recent years. One example is a quadruped robot; the quadruped robots are regarded as a new-trend among many robotics researchers [1,2]. Due to the large base of support and high degrees-of-freedom (DoF) in the motion control, the quadruped robots show good gait stability and controllability compared to any other robotic systems. For examples, Boston Dynamics introduced Big Dog and LS3, which are able to walk on uneven terrain at the speed of about 4.4 miles per hour [3][4][5], and SpotMini, which can autonomously navigate and move in general areas, such as flat ground, rough terrain, office/construction sites, passing through steps, and avoiding human/obstacles [6][7][8]. MIT introduced the MIT Cheetah series, which are quadruped robots inspired by a cheetah, which is equipped with directly driven actuators [9]. MIT Cheetah shows a broad range of leg movement and more focus on versatility [10,11], and MIT Mini Cheetah is capable of acrobatic motions, such as a back-flip, in addition to all the capabilities of its predecessor [12,13]. Additionally, the Korea Institute of Industrial Technology developed a quadruped robot actuated by hydraulic actuators, which is able to walk at about 3.4 miles per hour with a 60 kg load [14].
For the locomotion of the quadruped robot, various trajectory generation methods have been studied. For example, a central pattern generator (CPG) is often utilized to solve this problem. The CPG consists of coupled nonlinear dynamic equations, which are usually formulated in a neural network framework. The CPG generates rhythmic control policies for the locomotion of the quadruped robot [15]. Crespi et al. utilized the CPG to generate the locomotion of a robot inspired by a snake [16]. Another approach trajectory generation method for the robotic system is a finite states machine (FSM). FSM is a computational model of the machine with an initial internal memory that defines some finite states for the machine. Liu et al. applied the FSM for the separation of the swing phase and the stance phase [17]. Hussain et al. elaborated and tested the functional structure of a control system using the logic-labeled FSM [18], and Ding et al. generated the trotting, bounding, and aperiodic motions by event-based FSM that was extended from time-based FSM through the contact detection algorithm [19].
However, these conventional trajectory generation methods (i.e., CPG, FSM, etc.) have limited degrees of freedom, because gait phases (i.e., walk, amble, trot, canter, and gallop) of the quadruped system is changed very frequently according to the gait speed. Quadruped animals effectively accelerate and decelerate the gait speed by changing their gait phases. Since quadruped animals have been evolved over millions of years, the gait phases and the associated leg motions of quadruped animals can be regarded as optimal motion for the locomotion of the quadruped animals. For example, the walk is an optimized leg motion for low-speed locomotion, and the gallop is optimized for high-speed locomotion. The efficiency of the gait phases of quadruped animals was investigated by many biologists [20,21]. Therefore, it is reasonable that the trajectory generation for the quadruped robot should realize multiple gait phases to handle the locomotion speed of a wide range. Furthermore, to efficiently change the speed of the quadruped robot, the gait phase transition should be realized smoothly.
Since the conventional trajectory generation methods have been either entirely or partly based on pre-programmed network topologies, the gait transition from low to high-speed locomotion is not possible in the conventional trajectory generation methods. Furthermore, the conventional trajectory generation methods do not consider the realization to a sufficient degree, leading to pattern generation that neither self-organizes nor is efficient in response to real-world scenarios. Although adaptive trajectory generation methods have been developed to change gait phases adaptively and to qualitatively describe the behavior of animals, the number of parameters that can be adjusted may still be limited [22]. Moreover, the result of the trajectory generation, which is the limit cycle of the coupled nonlinear dynamic functions, is not fully adaptable. Another challenge in the use of the conventional trajectory generation methods is stability; the stability of synchronization between the result of the trajectory generation is challenging to prove in general. Therefore, the conventional trajectory generation methods may not be appropriate for the quadruped robot. That is why the general quadruped robots that have been developed and announced do not consider the transition of the gait phase or use a pre-determined fixed gait phase. The SpotMini of the Boston Dynamics, which is known to be the best in this field, also had pre-determined gait phases (walk and trot) applied, and the transition of the gait phase has to be completed before the start of the locomotion [23]. Although some research teams have developed trajectory generation algorithms that take into account the transition of the gait phase, the transition of the gait phase is limited (mostly from walking to trot), and the transition according to gait speed is not considered [24,25].
In this paper, a multi-phase joint-angle trajectory generation algorithm is proposed for the quadruped robot. The main proposed ideas of the multi-phase joint-angle trajectory generation algorithm are (1) there exist the cyclic basis functions of the joint-angle trajectories of the quadruped animal; (2) there exists a unique ground contact sequence for each gait phase (i.e., walk, amble, trot, canter, transverse gallop, rotatory gallop) ac-cording to the gait speed; (3) quadruped animals effectively accelerate and decelerate the locomotion speed by smooth gait transition in their gait phases; and (4) these ideas have to be considered to generate the joint-angle trajectories for the quadruped robots. To realize these proposed ideas, the joint-angle trajectories for the control of the quadruped robot in multiple gait phases are designed by a single cyclic basis function, which is designed based on the experimental results with a dog. The input to the basis function is controlled to realize desired trajectories using a static synchronization function, and the synchronization function is formulated by which the motions of legs are designed to have proper ground contact sequences (there exists a unique ground contact sequence for each gait phase, which has been typically used to distinguish the gait phases of quadruped animals [26]) with each other. Furthermore, to realize the smooth gait transition according to the gait speed, a fuzzy logic method is applied to the proposed method. The desired ground contact sequence for each gait phase is determined according to the literature in [27,28], and they are formalized as a time-varying vector with respect to the gait speed using fuzzy logic.
This paper is organized as follows. In Section 2, studies on the gait of the quadruped animal, particularly the ground contact sequence on each gait phase and characteristics on the stance and swing periods, are explored. In Section 3, the proposed gait phase generation algorithm is introduced. In addition, the reference joint-angle trajectories inspired by the quadruped animal are introduced, and the trajectories are represented as a cyclic basis function. In Section 4, the ground contact sequences are interpolated smoothly by fuzzy logic with respect to the gait speed. Then, the ground contact sequence is transformed into a desired input of a cyclic basis function using a static synchronization function. Lastly, the desired joint-angle trajectories are decided using the cyclic basis function, and the developed algorithm is discretized for implementation to the simulation. In Section 5.2, the simulation result of the proposed algorithm is shown. Finally, Section 6 renders a summary and future works.

Experiment Settings
Success in the realization of gait motions of the quadruped robot is highly dependent on the adequacy of generated joint-angle trajectories, as well as the mechanical performance of actuators, the power capacity, and the performance of control algorithms. The jointangle trajectories of the quadruped robot are complicated and different according to gait phases, which creates a challenge in their mathematical derivation. Therefore, in this paper, joint-angle trajectories for the control of the quadruped robot are inspired by and designed from the joint-angles of an animal that runs fast exhibiting various gait phases, as shown in Figure 1. The gait motions of an animal can be regarded as optimal motions for effective locomotion of the animal. Therefore, observation of the animal's gait motions may provide fundamentals for the realization of the successful locomotion of a quadruped robot. While quadruped animals walk or run, the joint-angle trajectories exhibit certain patterns on each gait cycle. The joint-angle patterns may provide fundamentals in the control of the joints of the quadruped robot, i.e., the joint-angles of the quadruped animal may be utilized as the reference input for feedback control.
The joint-angle trajectories are obtained from experiments with a fast-running dog in this paper. The species of the dog used in the experiments is American Akita, which is regarded as one of the fastest dogs; American Akitas exhibit the multiple gait phases, as shown in Figure 1, and are able to run at the gait speed of forty kilometers per hour.  During the gait of the quadruped animals, the whole body segments (e.g., the tail, the head, and so on) are actively involved. For example, it is well-known that the tail has a function to maintain the balance of the body in high-speed locomotion. Nevertheless, the four legs are still major components among the body segments that realize the gait, and thus, the leg motions are mainly analyzed in the experiments for the design of a trajectory generation algorithm of a quadruped robot. In addition, since the leg motions during high-speed locomotion are mainly on the sagittal plane, the joint-angle trajectories on the sagittal plane are considered in this paper.
The experimental setting is shown in Figure 2. To observe the gait of the dog from experimental results, the joint-angles should be measured. For this purpose, a high-speed camera, TS3 of Fastec Co. [29], which takes a video at 98,000 frames per second, was utilized. In order to obtain the absolute positions from the projected images obtained by the high-speed camera, red markers were attached to the locations of the dog's joints, as shown in Figure 1; the locations of the red markers are also marked as in Figure 2. Since the camera takes only a set of two-dimensional images, the camera was installed far enough, and the dog was forced to run on a straight line by a trainer. From each frame of videos, the absolute positions of the red markers were extracted, and the joint-angles were calculated by arc tangent functions. Considering that the four legs are the main components that generate locomotion, the majority of gait motions can be described by the twelve joint-angles (i.e., three joint-angles per leg), which are references for controlling a quadruped robot. Since the motions of the left and right legs share the same characteristics, and those of the front legs and the hind legs show significantly different characteristics due to the mechanical configuration, the proposed algorithm generates six unique joint-angle patterns, i.e., three joint-angle trajectories for the front leg, and three joint-angles for the hind leg. The joint-angle trajectories of the left and right legs can be generated by shifting the trajectories with an appropriate time-interval. Therefore, in this paper, the characteristic patterns of the six joints of the dog are observed and analyzed for formulating reference angles. The measurements of the joint-angles of the dog follow the general rule in the anthropometry except on the hip and shoulder joints. They are measured based on the vertical line for the sake of simplicity. In Figure 2, F and E represent the directions of flexion and extension, respectively. Figure 3 shows the obtained joint-angles of a dog in multiple gait phases. Since the period of one gait cycle, i.e., the time-interval between heel-strikes for each leg, varies according to the gait speed, the joint-angles shown in the figure were normalized to a gait cycle. Moreover, the period of a gait cycle should be multiplied to the percentage of the gait cycle in order to obtain the joint-angle information in the actual time domain. It should be noted that a unique pattern is observed in each joint-angle trajectory shown in Figure 3. For example, the shoulder and hip joints exhibit pendulum motions regardless of the gait phase, and the elbow joints show major flexion during swinging and minor flexion during stance. The similar patterns are also shown in the wrist joint angles, the knee joint, and the hock joint angles. Therefore, the observed joint-angles can be formalized into a representative function. However, the amplitudes and patterns of the joint-angle trajectories still show large variations as the gait speed changes, which creates a challenge in the derivation of a representative function.

Stance and Swing Periods
The period that a foot touches the ground is called the stance, and the period that the animal moves the foot in the air is called the swing. In Figure 3, the gait cycle starts with the stance period, and the swing period starts at the location of the continuous vertical line shown in the figures. The dashed vertical lines around the continuous vertical lines represent the standard deviation for each joint angle in each gait phase. The detailed analysis results on the observed leg motions of the dog are shown in Table 1. For the generalization of the obtained data, the gait speed was normalized with the body length of the dog. Since the body length of the dog from the shoulder to the hip is 0.6 (m), the actual gait speed can be obtained by multiplying this number. As the gait speed increases, the stance period is reduced, which implies that the time that the corresponding foot is in contact with the ground is shortened. On the other hand, it should be noted that the change of the swing period is minor (see the values of the Swing period row in Table 1). This is because (1) gait stability is achieved mainly by stance legs, and (2) the swing period is minimal for the rapid recovery of the leg position. Since the standard deviation of the swing periods were 0.02∼0.09 s, it is reasonable to assume that the swing period is constant regardless of the gait speed and phase. Consequently, the average ratio of the stance period is decreased, as the gait speed increases.

Joint-Angles in a Scaled Gait Cycle Domain
For the formalization of the obtained joint-angles in multiple gait phases, the time axes of the joint-angle trajectories are scaled. It should be noted that the characteristic patterns are shown in the joint-angle trajectories of the stance and swing periods, respectively. For example, in the elbow joint-angles, although their shapes on the actual time domain are different according to the gait speeds, the major flexion is observed in every swing phase, and the minor flexion is shown in every stance phase. Therefore, the time axes of the obtained joint-angles are scaled such that [0, 50)% of the gait cycle is the stance period, and the remaining [50, 100)% of the gait cycle is the swing phase.
The joint-angle trajectories of the dog on the scaled gait cycle domain are shown in Figure 4. It should be noted that the shapes of the joint-angle trajectories, as well as their patterns, match for all the gait phases. Therefore, the joint-angle trajectories on the scaled gait cycle domain can be formulated by a single represented function, named a cyclic basis function in this paper.

A Cyclic Basis Function
It should be noted that the gait cycle is scaled in Figure 4 in order to formulate the joint-angle trajectories at different gait speeds into a single representative function (i.e., one cyclic basis function for each joint). When the gait cycle is scaled such that the stance and the swing sections are divided equally, the shapes of the joint-angle trajectories match, which can be represented by a single cyclic basis function. Any function is applicable for the cyclic basis function, as long as it is continuous and differentiable in the gait cycle domain: where i = 1, 2, . . . , 6, where each number represents the shoulder, the elbow, the wrist, the hip, the knee, and the hock joints. φ is the index on the scaled gait cycle domain. In order to emphasize the cyclic nature of the joint-angle trajectories, trigonometric functions are utilized for the formalization of the joint-angles, i.e., where s F is a scaling factor that determines the amplitudes of joint-angle trajectories. The cyclic basis functions for the hind leg are where s H is a scaling factor. Figure 4 shows the cyclic basis functions in Equations (2)-(7). Notice that the designed functions represent the patterns and shapes of the joint-angles in multiple gait phases. Notice that the cyclic basis functions in Equations (2)-(7) are cyclic (i.e., f i (φ) = f i (φ + 2πn) for any integer number n) and differentiable for φ. Such characteristics are only valid if the scaled gait cycle is further scaled into the range of [0, 2π). Since the scaled gait cycle in Figure 4 divides the stance period and the swing period equally, [0, π) corresponds to the stance period, and [π, 2π) represents the swing period in the scaled gait cycle domain. The cyclic basis functions are to generate cyclic and continuous joint-angle trajectories as φ increases continuously. For a better representation of the cyclic nature, the cyclic basis functions can also be represented in the r-φ coordinate system, as shown in Figure 5, where the radii at an angle (φ) are the values of the joint-angles in the scaled gait phase. The gray area in Figure 5 represents the stance period.

Scaling Factor
As the gait speed increases, the shape of the joint-angle trajectories remains the same in the scaled gait cycle domain, but their magnitudes change. Notice that scaling factors are considered in the cyclic basis functions in Equations (2)-(7), where s F is a scaling factor for the three joint-angles of the front leg, and s H is one for the three joint-angles of the hind leg. As the gait speed increases, the range of joint-angles is increased, which may require increasing the scaling factors according to the gait speed. Through the experimental results, it was observed that the amplitudes of the hind leg motions are increased proportionally, as the gait speed increases. Therefore, the scaling factors, s F and s H , are formulated as where v(t) is the gait speed in the unit of meter per second, and q F and q H are constants.
In the case of the dog used in the experiments, q H = 0.3 and q F = 0.0, which means that the amplitudes of the hind leg motions change according to the gait speed, but those of the front leg motions are constant regardless of the gait speed.

Ground Contact Sequences
The quadruped locomotion is a combination of motions of the four legs. While the joint-angles of each leg can be analyzed and simulated by the corresponding cyclic basis functions, the ground contact sequences and their time-intervals should also be considered for the generation of overall locomotive motions. The ground contact sequences and time-intervals vary according to the gait phase. Namely, there exists a characteristic ground contact sequence defined for each gait phase (e.g., walk, trot, canter, gallop, etc.), as shown in Figure 6 [26], and they should be effectively considered in the joint-angle trajectory generation. Figure 6. Ground contact sequences with respect to gait phases; LF, LH, RF, and RH refer to the left-front, the left-hind, the right-front, and the right-hind, and d means the continuous ground contact time ratio of one gait cycle period.
In the ground contact sequences shown in Figure 6, the intervals between ground contact instances are determined with respect to the LF leg, and thus, the intervals for the LF leg are all zero. The intervals shown in the figure are normalized with a gait cycle period; the period of one cycle may be multiplied to the numbers shown in the figure to obtain the actual time-intervals. The normalized intervals are used as the reference intervals of the proposed gait motion generation algorithm and are denoted as R d , i.e., The second component, R lh d , means the normalized time-interval between the LF leg and the LH leg. Similarly, R r f d and R rh d mean the normalized time-intervals between the LF and the RF leg and the RH leg, respectively. The values of the components of R di , where i = 1, 2, . . . , 5 represents each gait phase, is listed in Table 2. Remember that the intervals in Table 2 are normalized values, as in Figure 6; thus, the period of one gait cycle should be multiplied to the values. For example, assuming that the period of one gait cycle is T cycle seconds, and the gait phase of a robot is in a walking phase, the locomotion of a robot is to be controlled such that the LF foot touches the ground first, and the remaining three feet touch the ground with the time-interval of 0.25T cycle seconds in the sequence of RH, RF, and LH feet.
Notice that the reference intervals, R di s, shown in Table 2, are defined for the five phases; the walk, trot, canter, transverse gallop, and rotatory gallop. The amble phase in Figure 6 is not taken into account in R di s because the normalized intervals of the amble phase are the same as the walk phase, and the period of the amble phase is too short that it is negligible in the locomotion of typical quadruped animals.

Gait Transition by Fuzzy Logic
Although the definitions of gait phases are discrete, as in Table 2, the gait speed of actual animals continuously varies over time, and their leg motions are smooth and continuous. In a robotic system, the smooth and continuous leg motions are also desired for natural and effective locomotion. Otherwise, the abrupt change of joint-angle trajectories may cause actuator saturations and instability of the gait.
In the proposed algorithm, the reference (i.e., desired) ground contact sequences for each gait phase is R di . The discrete values should be effectively interpolated with each other, such that the value of every component of a resultant reference ground contact sequence, R d f (v), is smooth and continuous (i.e., differentiable for the time). For this purpose, a Mamdani type fuzzy logic system is applied to the proposed algorithm [30][31][32].
The characteristic speeds that distinguish the desired gait phases may vary according to the species of quadruped animals. For example, the locomotion of the dog used in the experiments was definitely in the walk phase at the speed of 1.19 (m/s) (i.e., 1.98 bodylength per second), the trot phase at the speed of 2.69 (m/s), and the gallop phase at the speed of 4.34 (m/s), as shown in Table 1. Namely, if the gait speed is the same as one of the characteristic speeds, then the desired gait phase can be determined deterministically. However, the transition of the gait phases may occur in between the characteristic gait speeds. In this paper, the transition of gait phases is realized by the linear interpolation of two gait phases adjacent to the current gait speed. Following this logic, the likelihoods of the gait phases are modeled as membership functions (i.e., µ i , where i = 1, 2, . . . , 5 means each gait phase: walk, trot, canter, transverse gallop, and rotary gallop), as shown in Figure 7a. v iS is the starting point of the likelihood changing at each gait phase, and v iE is the ending point of the likelihood changing at each gait phase.
where r representatives the changing ratio of the likelihood. In this condition, the membership functions of each gait phase are defined as The magnitudes of the membership functions represent the likelihood of each gait phase. For example, at the gait speed of 0.5 bodylength per second, µ 1 = 1 and µ 2,3,4,5 = 0, which implies that the desired gait phase is definitely the walk phase. On the other hand, if the gait speed is nearly by 3 bodylength per second, µ 1,2 = 0.5 and µ 3,4,5 = 0, it means that the likelihoods of both the walk and trot phases are the same. Notice in the figure that the characteristic speeds are represented as ranges, i.e., µ i = 1 and  The membership functions for modeling the likelihoods of gait phases, however, are not smooth. Since the proposed algorithm changes the gait phases (i.e., gait transition) by changing the membership functions based on the fuzzy logic, the non-smooth characteristic of the membership functions can cause a drastic change in the ground contact sequence. This is the main problem of gait stability, such as dynamics change, control performance, and so on. Therefore, the membership functions have to generate as smoothly as possible. For further smoothing of the membership functions, new membership functions are defined as functions of µ i (v)'s, i.e., where i = 1, 2, . . . , 5 is the gait phase index, which is the same as in µ i 's, and s ∈ R + is a sensitivity factor that changes the slope of the membership functions in transition. The shapes of M i (v)s are represented in Figure 7b. If s is too small, the deterministic section (i.e., the section where M i (v) = 1) does not appear. On the other hand, if s is too big, the phase transition becomes too drastic. Thus, an appropriate s should be used; an appropriate value for s was empirically found to be 6.0. It should be noted in Equation (14) that M i (v) is normalized by the sum of all M i (v)s, which is to make the sum of all M i (v)s one in any situation. In the proposed algorithm, M i (v) is regarded as the likelihood of ith gait phase.
Since the total likelihood of gait phases (i.e., ∑ 5 i=1 M i (v)) is always one, M i (v)s can be used as weighting factors for the calculation of the desired ground contact sequence, i.e.
Since M i (v) returns a value close to one in the characteristic speed range, i.e., the gait phase is not deterministic, M i (v)s smoothly interpolate the values of R di s adjacent to the given gait speed.

Synchronization of Leg Motions
The inputs of the cyclic basis functions in Equations (2)-(7) are controlled to generate appropriate joint-angle trajectories with desired ground contact sequences obtained in the previous section. Three different cyclic basis functions are applied for the three joints of a leg (i.e., Equations (2)-(4) for the front leg, and Equations (5)-(7) for the hind leg), while the inputs to the three cyclic basis functions are identical for each leg. The input of a cyclic basis function is an index on the scaled gait cycle domain, and thus, the locomotion of four legs is realized by determining four indexes on the scaled gait cycle domain, i.e., where each component is the index on the scaled gait cycle domain for each leg, i.e., the input of the cyclic basis functions. The components of Φ(t) should be effectively updated considering the desired ground contact sequences and the gait speeds. Namely, the adequacy of the joint-angle trajectories generated by the proposed algorithm is dependent on how the components of Φ(t) are effectively updated. In order to reflect the desired ground contact sequences in the calculation of Φ(t) in real-time, the desired time-intervals of ground contact instances (i.e., the moments of heel-strike) should be obtained first. Remember that the desired ground contact sequence, i.e., R d f (v) in Equation (15), includes the desired intervals normalized to one gait cycle. Therefore, the actual time-intervals are obtained by multiplying the total period of one gait where T cycle (v) is the total period of one gait cycle. Since the first component of R d f (v) is always zero, the time-interval for the LF leg, τ l f (v), is zero. Remember that the total period of one gait cycle varies according to the gait speed, as shown in Table 1. The total period of one gait cycle consists of periods of the stance (T st (v)) and the swing (T sw ), i.e., As discussed in Section 2.3, the period of the swing is constant regardless of gait phases and speeds. In the case of the dog used in the experiments of this paper, T sw = 0.34 s. On the other hand, the period of the stance is reduced as the gait speed increases. Noting that the body moves forward during the stance period, as shown in Figure 8, T st (v) can be calculated as where l leg is the length of a leg (i.e., the distances between the wrist and shoulder joints and the hock and hip joints), and ∆θ e is an effective angle change during the stance period.
Since the angle variation of the shoulder and hip joints are dominant in the stance period, ∆θ e can be regarded as ∆θ shoulder and ∆θ hip for the sake of simplicity. Once the total period of one gait cycle is obtained as in Equation (18), the input of the cyclic basis functions can be calculated. Remember that the cyclic basis functions are on the domain of the scaled gait cycle. Therefore, the rate of the scaled gait cycle index in the stance period (i.e., [0, π) on the scaled gait cycle domain) may be different from that in the swing period (i.e., [π, 2π) on the scaled gait cycle domain). The two rates are defined as where ω st (v) means the rate of φ(t) in the stance period, and ω sw is the rate in the swing period. Figure 9 depicts how the index of the scaled gait cycle domain is updated in the actual time domain. It can also be formalized as where t 0 = ω −1 st π. Since the gait motions are cyclic, g(t) should also be cyclic with the period of T cycle (v), i.e., g(t + nT cycle (v)) = g(t) (23) where n is any integer number. Figure 9. Rates of the input of cyclic basis functions; t 0 = ω −1 st π, ω sw and ω st are the rates in the swing and stance periods, and T cycle means the total period of one gait cycle.
The four indexes in the scaled gait cycle domain are updated by Equation (22) with the desired time-intervals in Equation (17). Therefore, the inputs of the cyclic basis functions are calculated as

Discretization for Implementation
The overall procedure of the proposed joint-angle trajectory generation algorithm is shown in Figure 10. Firstly, for a given (i.e., current) gait speed, R d f (v) is obtained by interpolating the desired ground contact sequences defined for each gait phase, R di . A fuzzy logic method is applied for a smooth and continuous transition of the gait phases. Then, the rates of the scaled gait cycle indexes are calculated for each leg, as in Equation (22). With the calculated rates, the inputs of the cyclic basis functions are updated with the desired time-intervals. Finally, the joint-angle trajectories are obtained by the cyclic basis functions for each joint of each leg.
Since the proposed algorithm is to be implemented in a digital computer, it should be discretized. Notice that the majority of the variables in the proposed algorithm, such , and so on, are all functions of the gait speed, except the inputs of the cyclic basis functions in Equation (24). Therefore, the overall algorithm can be implemented and simulated in real-time once g(t) is discretized into g(k), i.e., where T is the sampling period. As g(t) is cyclic, g(k) is also cyclic, i.e., where n is any integer number.

Simulation Environments
The proposed joint-angle trajectory generation algorithm for quadruped robots is verified by simulation studies in this section. In order to obtain reasonable information for verifying the proposed algorithm accurately in the simulation, the dynamic model of the simulation is necessary. In this simulation, Newtonian mechanics were applied to the dynamic model [33][34][35]. Newtonian mechanics describe the motion of rigid bodies. Additionally, Newtonian mechanics provide very accurate results as long as the model of a mechatronic system is an accurate interpretation of the actual physics. The detailed parameters of the Newtonian mechanics in this simulation are shown in Table 3. Furthermore, indirect inverse dynamics based on a feedback control approach were used [36][37][38]. In the indirect inverse dynamics method, the input (i.e., the actuation force) is determined by the feedback controller such that an error between the desired motion by the proposed algorithm and the simulated output is small. In this simulation, the PD (proportional-derivative) feedback controller was applied as where u is the input, θ d is the desired joint-angle trajectory by the proposed algorithm, θ m is the simulated output by Newtonian mechanics model, and ω is the angular velocity. Furthermore, K P and K D were 15,000 and 2500, respectively, in this simulation. The simulation results were carried out using Matlab with the sampling period of 1 millisecond.

Simulation Results
The gait speed was set to smoothly increase from zero to 7.5 bodylength per second, which is 4.5 (m/s) for the bodylength of 0.6 (m), as shown in Figure 11. The first step of the proposed algorithm is to identify the likelihoods of gait phases for the current gait speed (i.e., M i (v)) and to determine the desired ground contact sequences (i.e., R d f (v)). This process was introduced in Section 4.1. Figure 12 shows the four components of R d f (v) with respect to the simulation time. Since the first component of R d f (v) was used as the reference, it was zero for the entire time range, as shown in the figure. The remaining three components of R d f (v) showed certain values with the desired intervals from the first component. It should be noted that the values of R d f (v) were continuous and smooth for the entire time range, which implies that the smooth and continuous gait phase transition was achieved successfully.  In addition to the identification of the likelihoods of gait phases and the desired ground contact sequences, the total period of one gait cycle is necessary for generating an input of cyclic basis functions, i.e., the index on the scaled gait cycle domain. Figure 13a shows the total period of one gait cycle calculated as in Equation (18). Notice that T cycle (v) was reduced as the gait speed increased, which implies that the cyclic locomotive motions occur more frequently at high speeds. In addition, a ratio between the swing and the stance periods was also changed according to the gait speed; the higher the gait speed, the lesser the ratio of the stance period, as shown in Figure 13b. Consequently, the ground contact time (i.e., the stance period) was shortened significantly as the gait speed increased, which may affect gait stability at high speeds. Since the total period of one gait cycle and the desired ground contact sequences were obtained, as shown in Figures 12 and 13b, respectively, the desired time-intervals in the actual time domain can be calculated, as in Equation (17). Figure 14a,b show the desired time-intervals (i.e., τ l f (v), . . . , τ rh (v)) calculated in the simulation. Notice that the time-intervals changed as the gait phases switched, and their magnitudes were reduced as the total period of one gait cycle was reduced. Nevertheless, the values of the time-intervals were smooth and continuous in the entire time range and successfully represented the desired ground contact sequences.
(a) From walk to trot  Based on the calculated desired time-intervals, the scaled gait phase indexes can be obtained for the calculation of joint-angle trajectories, as shown in Figure 15. Remember that the scaled gait phase indexes are the input variables of the cyclic basis functions that represent the shapes of joint-angle trajectories. The scaled gait phase indexes continuously increased in the object's continued locomotion; the values shown in Figure 15 are wrapped into the range of [0, 2π) for better representation of the figure. The wrapping is automatically achieved in the proposed algorithm because g(t) in Equation (22) (or, g(k) in Equation (25)) is cyclic over [0, 2π). The patterns shown in Figure 15a,b are different because the desired time-intervals and the rates of the indexes were different due to the gait speeds. The line shape of φ i in one gait cycle is not linear, and the lines are bent in the middle of the gait cycle. It is because the kernel function makes Φ(k) from T(k). The kernel function composed of two linear functions that have different slopes: ω st and ω sw . Thus, each φ i has the same shape of the kernel function. Angle (rad) (a) From walk to trot Angle (rad)  Figures 16-18 show the generated joint-angle trajectories in transitions of gait phases. It should be noted that the transition of the gait phases (i.e., the changes of the joint-angle amplitudes and intervals) took place smoothly in the generated joint-angle trajectories. The time axis in Figure 16 corresponds to the range of v 1E and v 3S ; remember that the gait speed was set to monotonically increase, as in Figure 11. Therefore, Figure 16 shows the gait phase transition from the walk phase to the canter phase, while passing through the trot phase. On the other hand, the time axis in Figure 17 corresponds to the range of v 3E and v 4S , and thus, Figure 17 shows the gait phase transition from the canter phase to the transverse gallop phase. The time axis in Figure 18 shows the range of v 4E and v 5S , which means that Figure 18 verifies the gait phase transition from the transverse gallop phase to the rotatory gallop phase. The sub-figures labeled as (a) and (b) represent the generated joint-angle trajectories for the front leg and the hind leg, respectively. The simulation results in Figures 16-18 verify that the gait transition was realized smoothly and naturally without abrupt changes in the generated joint-angle trajectories.     In order to observe the ground contact sequences in the simulation results, the stance and swing phases were analyzed, as shown in Figure 19. In the figure, the value of 1 means the stance phase, and the value of 0 represents the swing phase. Notice in the figure that the width of the stance phase was deceased as the gait speed increased, while that of the swing phase remained the same regardless of the gait speed. For each of the five characteristic speeds (i.e., v 1,2,...,5 ), one gait cycle is highlighted by gray boxes, as shown in Figure 19. In each gray box, the ground contact sequences can be observed as marked by numbers 1 , 2 , 3 , and 4 in the figure. It should be noted that the numbers match the ground contact sequences shown in Table 2. The generated joint-angle trajectories were applied to the three-dimensional model of a quadruped robot, as shown in Figure 20. Notice that the simulation model successfully realized locomotive motions in multiple gait phases. The proposed algorithm generated the joint-angle trajectories in real-time without changing any parameters in the simulation study.

Conclusions
In this paper, an algorithm that generates joint-angle trajectories for the control of a quadruped robot was proposed. The proposed algorithm was inspired by a quadruped animal and designed to realize locomotive motions by cyclic basis functions that represent the characteristic natures of the joint-angle trajectories of the animal. The proposed algorithm followed a hierarchical procedure: in the lowest level, it determined the most appropriate gait phase for a given gait speed. When the appropriate gait phase was not deterministic (e.g., the transition of gait phases), the fuzzy logic smoothly interpolated the likelihoods of two gait phases adjacent to the current gait speed. The estimated likelihoods of gait phases were utilized as weighting factors in the calculation of a resultant desired ground contact sequence. At the same time, the proposed method also calculated the stance and swing periods for the given speed, and the total period of one gait cycle was multiplied to the desired ground contact sequence in order to obtain the desired time-intervals in the actual time domain. At the middle level of the proposed algorithm, the indexes on the scaled gait cycle domain were calculated based on the desired gait phase (i.e., the desired time-intervals) and the gait speed. The indexes were the inputs of the cyclic basis functions. At the highest level, the joint-angle trajectories were calculated by the cyclic basis functions. The proposed algorithm was verified by the simulation results. It generated the joint-angle trajectories smoothly and continuously, even when the gait phases changed frequently. Since there was no abrupt change in the generated joint-angle trajectories, the proposed algorithm provided quadruped robots with effective reference joint-angle trajectories.
However, the proposed algorithm generates joint-angle trajectories considered in the sagittal plane only. In the real robot system, joint-angle trajectories have to be generated with considering three-dimensional space, including the frontal and transverse planes that are closely related to gait stability, balance, and locomotion performance (i.e., curve running, pass through the tilted and rough terrain). To generate the joint-angle trajectories in the three-dimensional space, the degree-of-freedoms should be added to the robot system, such as the abduction and adduction of the hip joint, rotation and flexion of the spine, movement of the pelvis and scapula, and the joint-angle trajectories of additional degree-of-freedoms should be obtained and formulated. Furthermore, the joint-angle trajectories of additional degree-of-freedoms will change according to the gait speed similarly to the degree-offreedoms in the sagittal plane. Therefore, this paper has a meaningful contribution to the starting line of the research about obtaining and formulating the joint-angle trajectories of a robot system by taking account of the gait transition based on natural biology.
The feasibility of the proposed ideas was verified through the simulation results, but experimentation of the real robot system was not dealt with in this paper. The following research will introduce the proposed method's experimental results with a high-speed quadruped robot. To progress the experiment of the high-speed quadruped robot, a feedback control method that robustly controls actuators to achieve the desired motions and guarantees the gait stability by reacting to disturbances has to be considered. In the real world, the more demanding conditions robots are exposed to, the more challenging the control of robotic legs becomes. As the gait speed increases, the frequency bandwidth of reference signals (i.e., desired motions) is also significantly enlarged. Thus, the significant enlargement of the closed-loop bandwidth of a control system is necessary. Consequently, high gain control with a high-speed data acquisition system is unavoidable, which also necessitates the accurate measurement of physical states without noise. Moreover, the dynamics of a robotic leg are highly uncertain and time-varying due to the ground contact and the ground condition, as well as the unknown payload. Therefore, a control algorithm that can achieve both great control performance and stability robustness is essential for the effective control of a quadruped robot.
Furthermore, an external disturbance (e.x., ground reaction force, impact, change of the load condition. etc.) has to be considered in terms of the control method. A disturbance observer (DOB) is a good choice for dealing with external disturbance. The DOB estimates external disturbances by comparing an actual measurement with a simulated output. By filtering the output discrepancy with an inverse of the system model, the estimation of the exogenous disturbance is achieved. The DOB can also be applied as a feedback controller, in which case the disturbance is rejected, and the overall system is controlled to follow the nominal model. Since it attenuates the model variation and makes the closed-loop system robust to external disturbances, the DOB is appropriate for the control of the quadruped robot.