An Improved ZMP-Based CPG Model of Bipedal Robot Walking Searched

This paper proposed amethod to improve the walking behavior of bipedal robot with adjustable step length. Objectives of this paper are threefold. (1) Genetic Algorithm Optimized Fourier Series Formulation (GAOFSF) is modified to improve its performance. (2) Self-adaptive Differential Evolutionary Algorithm (SaDE) is applied to search feasible walking gait. (3) An efficient method is proposed for adjusting step length based on themodified central pattern generator (CPG)model.TheGAOFSF ismodified to ensure that trajectories generated are continuous in angular position, velocity, and acceleration. After formulation of the modified CPG model, SaDE is chosen to optimize walking gait (CPGmodel) due to its superior performance.Through simulation results, dynamic balance of the robot with modified CPG model is better than the original one. In this paper, four adjustable factors (Rhs,support, Rhs,swing,Rks,support, andRks,swing) are added to the joint trajectories.Through adjusting these four factors, joint trajectories are changed and hence the step length achieved by the robot. Finally, the relationship between (1) the desired step length and (2) an appropriate set of Rhs,support, Rhs,swing, Rks,support, and Rks,swing searched by SaDE is learnt by Fuzzy Inference System (FIS). Desired joint angles can be found without the aid of inverse kinematic model.


Introduction
Recently, many approaches have been adopted for generation of bipedal walking gait.Some researches [1][2][3] adopted a simplified dynamic model to generate walking gait calculated through inverse kinematic model which is complex and hence the computation load is high.
Inspired by neural science, some researchers investigated central pattern generator (CPG).The prime reason for arousing their interest is that CPG models provide several parameters for modulation of locomotion, such as step stride and rhythm, and are suitable to integrate feedback sensors.Hence, a good interaction between the robot and the environment can be achieved [4].According to Ijspeert [4], CPG becomes more and more popular in robot community.Taga et al. [5] integrated feedbacks with neural oscillators for unpredicted environment.Yang et al. [6], Shafii et al. [7], and Yazdi et al. [8] utilized TFS to formulate ZMP-based CPG model as the basic walking pattern of bipedal robot.Or [9] presented a hybrid CPG-ZMP control system for flexible spine humanoid robot.Aoi and Tsuchiya [10] proposed a locomotion control system based on CPG model for straight and curved walking.Farzenah et al. [11] noted that many researches on CPG model are designed for specific motion only and thus cannot generate arbitrary walking gait, such as changing step length and proposed 31 TS-Fuzzy systems for adjusting speed and step length.
Ijspeert [4] and Gong et al. [12] stated that stochastic population-based optimization algorithms have been chosen to optimize parameters of CPG model in many studies.Genetic Algorithm (GA) [6,13], Genetic Programming (GP) [14], Particle Swarm Optimization (PSO) [7], and Bee Algorithm [8] are adopted in searching the parameters of CPG model.Besides the above-mentioned techniques, there are still other gradient-free optimization techniques.Storn and Price [15] proposed Differential Evolution (DE) and conducted comparisons with some prominent algorithms, such as Adaptive Simulated Annealing (ASA) and the Breeder Genetic Algorithm (BGA).DE outperforms the above-mentioned prominent algorithms in terms of least number of 2 ISRN Robotics generations for finding global minimum [15].Similar results are reported in the following studies [16][17][18].Hegery et al. [16] carried out comparisons between DE and GA on N-Queen and travelling salesman problem and concluded that the performance of DE is better.Tušar and Filipič [17] carried out comparisons between DE-based variants DEMO and basic GA on multiobjective optimization problem and their result showed that DEMO outperforms basic GA.Vesterstrøm and Thomsen [18] noted that DE outperforms PSO and Evolutionary Algorithms (EAs) on majority of numerical benchmark problems.DE consists of population size (NP), scaling factor (), and crossover rate (CR) which significantly affect the performance of DE [19][20][21][22].Different problems require different parameters and strategies for effective optimization.Even in the same problem, different regions of search space may require different strategies and parameters for better performance [20].It is time consuming to search the most appropriate strategy and parameters by trial and error.Hence, Omran et al. [21] and Brest et al. [22] have proposed different methods to adjust CR and F. However, appropriate mechanism for choosing suitable strategies is not considered in [21,22].Qin et al. [20] proposed Self-adaptive Differential Evolution Algorithm (SaDE) which can adjust CR,  and choose strategy automatically during optimization.SaDE outperforms conventional DE variants and the other adaptive DE, such as SDE [21] and jDE [22], in terms of higher successful rate.
Based on the above-mentioned findings, this paper focuses on (1) CPG model for trajectory generation and (2) providing an efficient method to adjust step length.In this paper, original CPG model proposed by Yang et al. [6] is adopted since it provides a good foundation for the goal stated in this paper.The angular velocity of trajectories generated by GAOFSF [6] is usually discontinuous which has an adverse effect on ZMP.As a result, GAOFSF is modified to ensure that the trajectory generated is continuous in angular position and its first and second derivatives.After formulation of modified CPG model, parameters of CPG model are searched based on kinematic and dynamic constraints.It shows that the problem can be formulated as a multiobjectives and multiconstraints optimization.Gradient-free optimization technique is chosen since a set of parameters is searched in a highly irregular and multidimensional space which cannot be handled by standard gradient-based search method [12].SaDE is chosen as the method for optimizing the walking gait of robot in this paper because (1) its performance is superior and (2) appropriate strategies and parameters are not chosen manually.Based on [6,23], step length can be varied by simply changing several adjustable factors of GAOFSF.Look-up table proposed by Yang et al. [6] is not adopted in this paper because (1) a lot of memory is occupied if tremendous data is stored and (2) arbitrary step length within specific range cannot be commanded to the robot.To deal with this problem, four parameters ( hs,support ,  hs,swing ,  ks,support , and  ks,swing ) are added to the modified CPG model and searched by SaD.Four FIS systems are used to learn the relationship between (1) desired step lengths and  hs,support ,  hs,swing ,  ks,support , and  ks,swing .Then, desired joint angles can be found without the aid of inverse kinematics.

Kinematic Model and Dynamic Model of Bipedal Robot
In Figure 1, it shows that the bipedal robot consists of 12 DoF.Each leg has 6 DoF.RL and LL represent right and left legs respectively, while 1, 2, 3, . . ., 6 represent joint number.Figure 1 assumes that right leg is the support leg while the left leg is the swing leg.Joint number  of support leg and swing leg corresponds to different joint shown in Table 1.Also, axis of local coordinates attached on different joints acts as rotation axis and the direction of rotation is determined by right-hand rule.Information of physical dimension of bipedal robot is measured and simplified based on a modified Kondo-3HV.The total mass of the bipedal robot is about 1.4 kg.The mass of upper trunk and lower body is about 0.4 kg and 1 kg, respectively.Since the mass of each link in lower body is almost the same, then their masses are simply obtained by 1 kg/12 = 0.083 kg.The height of each link is shown in Table 2.
Since Denavit-Hartenberg notation [24,25] and iterative Newton-Euler dynamic algorithm [26,27] are maturely developed and commonly used in many studies of bipedal robot, these two methods are adopted to formulate forward kinematic model and inverse dynamic model, respectively, in this paper.Since this paper focuses on bipedal walking on horizontal flat plane in sagittal plane (parallel to the   −   plane of global coordinate), only ZMP  shown in (1) is considered and acts as an indicator to evaluate the dynamic equilibrium of the bipedal robot.In Figure 2, a local coordinate is attached to the center of support foot to observe the variation of ZMP  with time. . (1)

Formulation of Modified ZMP-Based CPG Model
Hip and knee trajectories of GAOFSF [6] consist of two different sections.For hip joint, two different Truncated Fourier Series (TFS) are used to formulate the upper portion ( ℎ + ) and lower portion ( ℎ − ) of trajectory.For knee joint, TFS and lock phase are joined together to formulate the whole trajectory.Based on Figure 3, the angular position is observed to be continuous.However, angular velocity of searched trajectory generated by GAOFSF is usually discontinuous (1) at the transition between  ℎ + and  ℎ − ( 3 and  6 or  0 ) and (2) at the beginning ( 1 or  4 ) and the end of lock phase ( 2 or  5 ).Abrupt change in angular velocity has an adverse effect on ZMP and hence the dynamic equilibrium of bipedal robot.In this paper, the GAOFSF is modified to ensure that the trajectories generated are continuous in angular position and its first derivative and second derivative.[ and  ℎ − are different, two different TFS are required to formulate the hip trajectory.Yang et al. [6] have proposed 5th order TFS and showed that the amplitudes of 4th and 5th orders are too small which can be neglected.3rd order TFS is adopted.In ( 2)-( 3), angular velocity of θ − ℎ and θ + ℎ at  3 and  6 is set as equal to ensure angular velocity is continuous.Thus, two more orders are added to  ℎ − to satisfy this constraint and are calculated by ( 4)- (5).Consider

Knee Trajectory in Sagittal
Plane.In the knee trajectory, Yang et al. [6] proposed a lock phase ([ 1 ,  2 ] and [ 4 ,  5 ]) in knee trajectory which assumes constant joint angle and zero angular velocity and acceleration.This introduces abrupt change in angular velocity.The main goals to be achieved in the modified knee trajectory are (1) continuity in angular position, velocity, and acceleration and (2) the advantage proposed in GAOFSF that can be maintained.Then, a new formulation is proposed in the following equation: In (11), lock phase is canceled to ensure continuous angular velocity.By adjusting  [6], GAOFSF can achieve energy efficient and stable "bent knee" walking gait.This advantage is still remained in the proposed modified model.The general shape of the hip and knee trajectories generated by modified CPG model is shown in Section 6.

Ankle Trajectory in Sagittal
Plane.Right and left ankle joint trajectories are simply formulated as (12) to ensure that the trunk is upright and swing foot is parallel to the horizontal flat plane.Consider

Optimization of Basic Walking Pattern by SaDE
This section focuses on how to search the basic walking pattern of bipedal robot by SaDE [19,20].To simplify the whole process in this section,  hs,support ,  hs,swing ,  ks,support , and  ks,swing are set as 1.The following are the main objectives of basic walking pattern to be achieved through using SaDE.
(1) The desired step length is set as 0.05 m.
(2) Upper bound and lower bound of swing height are set as 0.02 m and 0.01 m, respectively.(3) Premature landing should not occur throughout the walking cycle.(4) ZMP is within the area of support polygon throughout the walking cycle.
The procedure of SaDE takes the following steps.Also, a flow chart of SaDE is shown in Figure 30 of the Appendix to facilitate the understanding of readers.Details of the procedure are discussed in Sections 4.1-4.7.
(2) Code the parameters to be searched to form target vector ( , ).(3) Initialize the first generation ( 1 ) of population ().( 4 where  is number of data.Strike velocity during landing should be as small as possible since impact between landing foot and the ground can cause mechanical wear of parts and unstable walking. 2 is stated as follows [6]: If ZMP  is within the area of support polygon [0.061, −0.061] throughout the walking cycle, dynamic equilibrium of the robot is satisfactory.To deal with the discrepancies between simulation model and physical test bed, a safety factor (0.01 m) is added to ensure good performance during experiment. 1 is formulated as follows: where   is length (0.121 m) of foot sole/2.The robot is assigned to walk forward along positive axis.To ensure correct direction,  2 and  3 are designed as follows: 2 = max (− trunk, , 0) ,  3 = max (− swing foot, , 0) , (16) where  trunk, is mean trunk velocity in -direction and  swing foot, is mean swing foot velocity in -direction.

ISRN Robotics
To achieve natural human-like walking gait, 6 is designed to ensure that the robot can achieve desired step length: To ensure reasonable swing height, S 7 is designed to ensure that the swing height is within the [0.01, 0.02] m.
≥  lower && max ( swing foot ) ≤  upper , S 7 = 0. Else desired −  swing Foot,      . End. ( 8 is designed to prevent the swing foot from penetrating through the ground and hence reduce the chance of premature landing while  9 is designed to ensure swing foot lands on the ground ( ground = 0 m) at landing time which is  2 and  5 : The total scores used for evaluating each target vector ( , ) are formulated as follows: These weightings are set by trial and error to ensure their importance is almost the same and walking gait with satisfactory performance can be obtained: = [2500, 285, 265, 100, 2200, 12000, 300, 960, 25000] , where  , is weighting of fitness functions and  , is weighting of constraints.

Initialization of Parameters of Target Vectors (𝑋
Similar approach is applied to adjust CR  and   .During learning period (LP = 5 generations), for strategy , CR  , and   , corresponding to the trial vectors ( , ) that successfully enters next generation are stored in database.Once learning period is completed, mean value (CR  and   ) of the stored value corresponding to strategy  is calculated and is used to generate a new set of CR   , and   , .Then, storage data of CR  , and   , are removed for the next learning period to eliminate the effects of past data.The coefficients of  joint,support/swing are solved by the following constraints:

Fuzzy Inference System (FIS).
Four FISs are adopted to learn the relationship between (1)  hs,support ,  hs,swing ,  ks,support , and  ks,swing and (2) desired step length (5 cm, 4 cm, 3 cm, 2 cm, 1 cm, and 0 cm).The architecture of FIS proposed by [28] is shown in Figure 4.The desired step length (()) is the input of FIS while  joint,support/swing is the output of FIS which is obtained through the summation of   (()).Based on [24],   (()) is a function of desired step and can be represented as a first order polynomial stated in (40).
In this section, the membership function (  ) is set as Gaussian function as follows: where    is degree of membership function of   ,   is parameter that affects the width of Gaussian function of   ,   is parameter that affects the center of Gaussian function of   ,   is firing strength of rule , and   is normalized firing strength.This section focuses on adjusting the consequent parameters ( ,1 ,  ,2 ) of   () by least square estimation (LSE) [28] stated in (41).Gradient descent algorithm is not adopted to adjust premise parameters since the size of training data ( = 6) is relatively small.Hence, premise parameters are fixed during the training process.In this paper,  1−5 are set as 0.2, 0.4, 0.6, 0.8, and 1, respectively, while  1−5 are set as 0.3.After 50 training cycles, an appropriate set of consequent parameters is found.Details of result in this section are stated in Section 6. Consider where  is vector of desired step length and  is vector of desired output.shows that step achieved by bipedal robot is the same as the desired step length (0.05 m).In Figure 8, maximum height of swing foot is satisfactory based on Table 3. Swing foot lands on the ground at desired landing time ( 2 = 0.7 s and  5 = 1.7 s).The joint trajectories of modified CPG model are shown in Figure 9. Based on Figures 10 and 11, angular position and velocity of each joint trajectory is observed to be continuous and smooth.Because of the heavy trunk mass, abrupt change in  trunk, should be prevented.From Figure 12, it shows that  trunk, is continuous and smooth except the landing time.

Dynamic Aspects of Modified CPG Model.
Figure 13 shows that dynamic equilibrium of the robot is satisfactory since ZMP  is within the area of support polygon [−0.061, 0.061] m.The length of foot sole is 0.122 m while ZMP  margin is defined to be (length of foot sole/2) − max(|ZMP  |).Based on Figure 13, ZMP  margin is observed to be at least 0.03 m which is large enough to allow larger step.ZMP  is smooth and the only abrupt change happens at the landing time ( 2 = 0.7 s and  5 = 1.7 s) since ZMP  shifts to the support foot of the next step at these two moments.that mutation strategy 2 and strategy 3 do not favor for searching feasible walking gait since they are suppressed completely after 60 generations (Figure 14). Figure 5 shows that scores of target vectors tend to converge in the range of generation 100 to 200.Strategy 1 can help trial vector converge faster since it involves the best trial vector in the mutation operation.Hence, from generation 100 to 200, the probability for choosing strategy 1 is always higher than that of strategy 4.

Results of Original CPG Model.
Except constraint 1, the original CPG model [6] is searched by SaDE under the same setting.Safety factor is set as 0 m since it becomes difficult to search a feasible walking gait for original CPG model if safety factor is set as 0.01 m.Based on Figure 15, the scores of the best trial vector is −1105.9.Also, the searched parameters are [0.05190.0012 0.0070 −0.5237 −0.0053 −0.0068 0.3579 0.1459 0.0025 0.0891 0.1736].The values of fitness functions and constraints are [0.0694ms −1 , 0.0484 ms −1 ] and [0.0002 m, 0 ms −1 , 0 ms −1 , 0 rad, 0 rad, 0 m, 0 m, 0 m, 0 m], respectively.Compared with the performance of modified CPG model, the original CPG model is less satisfactory.The walking gait of original CPG model searched by SaDE is visualized in Figure 16.Also, in Figure 17, desired step length (0.05 m) is successfully achieved.In Figure 18, maximum height of swing foot is satisfactory based on Table 3.Also, swing foot lands on the ground at desired landing time ( 2 = 0.7 s and  5 = 1.7 s).
The joint trajectories of modified CPG model are shown in Figure 19.Based on Figures 20 and 21, angular velocity of each joint trajectory is observed to be discontinuous.The standard deviation ( original = 0.0694 ms −1 ) of  Trunk, is larger than that ( modified = 0.018 ms −1 ) of modified CPG model.The prime reason is that, in Figure 22,  Trunk, is observed to be discontinuous because of discontinuity in angular velocity.In Figure 23, obvious abrupt change in ZMP  is observed.Also, ZMP  exceeds the area of support polygon ( 1 = 0.0002 m).This shows that dynamic equilibrium of bipedal robot is affected by discontinuous angular velocity.Although the area of foot sole can be made larger to tolerate the abrupt change of ZMP  , the agility of the robot is sacrificed.For the original CPG model, a larger step is not allowed since ZMP  has exceeded the support polygon when the robot walks with 5 cm step length.4. In Figure 24, it shows that the bipedal robot can achieve desired step length (4 cm, 3 cm, 2 cm, and 1 cm).Also, in Figure 25, the maximum swing height corresponding to different desired step length reaches reasonable swing height based on Table 3 and no premature landing occurs during walking.In Figure 26, ZMP  is always within the area of support polygon [0.061, −0.061].Hence, the dynamic equilibrium of    the robot for different desired step length is satisfactory.Then, an attempt is made to show that when desired step length is adjusted in the next step, ZMP  can still be maintained within the area of support polygon [−0.061, 0.061].The bipedal robot is at rest initially and is commanded to change its step length at  2 = 0.7 s.In Figure 27, it shows that if change in step length is larger, then change in ZMP  is larger.It also shows that, in the case of the largest change in step length (0 cm to 5 cm), ZMP  is still within area of support polygon.Hence, robot can keep its dynamic balance well when desired step length is changed in the next step.5. Arbitrary desired step lengths (0.041 cm, 0.033 cm, 0.025 cm, and 0.017 cm) within a specific range (0 cm-5 cm) are commanded to the robot.In Figure 28, it shows that the robot can achieve arbitrary desired step length successfully.In Figure 29, premature landing does not occur throughout the whole walking cycle.Hence, the relationship between (1)  hs,support ,  ks,support ,  hs,swing and  ks,swing , and (2) desired step length is learnt successfully.

Figure 2 :
Figure 2: Local coordinate system attached on the foot sole (top view).

Figure 8 :Figure 9 :Figure 10 Figure 11
Figure 8: Variation of height of swing foot with time (modified CPG model).

Figure 14 :
Figure 14: Probability of mutation strategy  to be chosen (modified CPG model).

Figure 15 :
Figure 15: Average scores and scores of best target vector of original CPG model.

Figure 16 :
Figure 16: Walking gait of original CPG model searched by SaDE.

Figure 18 :Figure 19 :
Figure 18: Variation of height of swing foot with time (original CPG model).

Figure 20 :
Figure 20: θ hs versus  hs of original CPG model.

Figure 25 :Figure 26 :0
Figure 25: Variation of height of swing foot with time (different desired step length).

Figure 27 :Figure 28 :
Figure 27: Variation of ZMP  corresponding to different desired step length change.

Figure 29 :
Figure 29: Height of swing foot (arbitrary desired step length command).

6. 7 .
Relationship between  hs,support ,  ks,support ,  hs,swing , and  ks,swing and Desired Step Length Learnt by FIS.The parameters of FLS tuned by LSE are shown in Table

Table 1 :
Joint number of support leg and swing leg corresponds to different joints.

Table 2 :
Height of each link in support leg and swing leg.

Table 3 :
Lower bound of desired swing height corresponding to different step lengths. ).Then, ST  and FT  are reset to zero for the next learning period to eliminate effect of the past data:

hs,support , 𝑅 hs,swing , 𝑅 ks,support , and 𝑅 ks,swing ) for Step Length Adjustment
The objective functions and constraints mentioned in Section 4 are adopted to search appropriate value of  hs,support ,  hs,swing ,  ks,support , and  ks,swing final corresponding to different desired step length (0.04 m, 0.03 m, 0.02 m, and 0.01 m).
5.1.Objective Functions and Constraints of SaDE.joint,support is changed to  joint,swing since the support leg is changed to swing leg in the next step.In order to have a smooth transition, fifth order polynomial shown in equation (34) is utilized.The period of transition time (  −  0 ) is set as 0.4 s which is 40% of time duration (1 s) of one step to balance between rapid transition time and good dynamic equilibrium.

Table 5 :
Parameters of FLS tuned by LSE.