Trajectory Planning of Flexible Walking for Biped Robots Using Linear Inverted Pendulum Model and Linear Pendulum Model

Linear inverted pendulum model (LIPM) is an effective and widely used simplified model for biped robots. However, LIPM includes only the single support phase (SSP) and ignores the double support phase (DSP). In this situation, the acceleration of the center of mass (CoM) is discontinuous at the moment of leg exchange, leading to a negative impact on walking stability. If the DSP is added to the walking cycle, the acceleration of the CoM will be smoother and the walking stability of the biped will be improved. In this paper, a linear pendulum model (LPM) for the DSP is proposed, which is similar to LIPM for the SSP. LPM has similar characteristics to LIPM. The dynamic equation of LPM is also linear, and its analytical solution can be obtained. This study also proposes different trajectory-planning methods for different situations, such as periodic walking, adjusting walking speed, disturbed state recovery, and walking terrain-blind. These methods have less computation and can plan trajectory in real time. Simulation results verify the effectiveness of proposed methods and that the biped robot can walk stably and flexibly when combining LIPM and LPM.


Introduction
Compared with other types of robots, humanoid robots have good adaptability to the environment, stronger obstacle avoidance ability, and a smaller moving blind area, which has attracted the attention and in-depth research of scholars [1][2][3][4][5][6][7][8][9]. At present, biped robots are still quite far away from the real sense of anthropomorphism, and there are many problems to be solved in this field. For example, due to the inherent instability of biped walking, walking stability analysis is still an important issue for biped robots. In addition, the biped robot is a high-order and strong coupling nonlinear system, which makes the trajectory planning and control difficult. The realization of stable walking is the primary task in the research of humanoid robots.
There are many methods for gait planning of biped robots. These methods could be divided into two classes. The first uses the accurate information of dynamical parameters to generate walking patterns. Joint angle trajectories [10][11][12][13][14] or trajectories of some key parts [15][16][17][18], e.g., hip and/or feet, are usually fitted by spline or polynomial functions, then the coefficients of spline or polynomial functions are determined by parameter optimization technique. However, these gait-planning methods need a lot of computation and cannot meet the requirement of trajectory planning in real time. The more degrees of freedom and the higher the order of the polynomials, the more computation time is needed for solving the optimization problem.
The other class is based on a simplified model to generate walking patterns. Inverted pendulum [19,20] is widely used because of its simplicity. A biped robot is usually regarded as a concentrated mass and massless leg. The trajectory of the center of mass (CoM) is planned with a simplified model, and then the angles of other joints are solved by inverse kinematics. One of the widely used methods is the linear inverted pendulum model (LIPM) [21,22]. The advantage of LIPM is that the trajectory of the CoM has an analytical solution. Moreover, its forward and lateral motions are decoupled. Another model is the inverted pendulum model (IPM) with constant leg length [23][24][25]. In this model, the CoM moves along an arc. Although the dynamic equation of IPM is simple, there is no analytical solution due to its nonlinearity. The disadvantages of LIPM and IPM are that they can only generate the trajectory of the single support phase (SSP), but cannot generate the trajectory of the double support phase (DSP).
From an application perspective, when the biped robot is walking outdoors, due to the unstructured ground environment, the robot is required to have the ability of real-time gait generation according to the current environment. However, the more accurate the model is, the more computation is needed. Hence real-time gait planning may become very difficult. Therefore, the simplified model is a feasible and very useful method for real-time gait planning.
On the other end of spectrum, there is little attention on the DSP. Many gait-planning methods consider only the SSP and ignore the DSP, or the DSP is assumed to be instantaneous. In this situation, the center of pressure (CoP) or zero-moment point (ZMP) [2] needs to transfer from the trailing foot to the leading foot instantaneously when the support leg is switched. This requires an impulsive force between the rear foot and the ground. The emerge of impulsive force could lead to some adverse factors: (1) it has a negative effect on the walking stability analysis; (2) generating impulse force needs a sufficiently large joint torque that the joint driving motors may not provide; (3) it may damage the hardware of the robot. The introduction of the DSP can reduce the impact between the foot and the ground, make a smooth ZMP transition from the trailing foot to the leading foot, and improve walking stability. In addition, the support polygon area of the DSP is larger than that of the SSP, so the ground can provide greater external torque to the robot during the DSP. Therefore, the robot has stronger state-adjustment ability during the DSP. During the SSP, because the robot's foot is small, the ground cannot provide a large enough external torque to avoid the robot falling down; as a result, the robot needs more adjustment of the internal state. Kajita and Tani reported that adding the DSP to the LIPM reduced the loss of the CoM's velocity when the support leg exchanges [26].
To overcome shortcomings of models without the DSP, some scholars introduced the DSP in gait planning. Kajita et al. [27] planned the CoM's trajectory of the DSP as a fourth-order polynomial function. The coefficients of the polynomial are determined by the boundary condition and the specified duration of the DSP. Motoi et al. [28] designed the CoM's trajectory in the DSP as a fifth-order polynomial function. The disadvantage of their method is that walking stability was not taken into consideration during the DSP. With the increase of the order of the polynomial, unexpected oscillation of CoM may occur; as a result, unexpected oscillation of the ZMP may occur during the DSP. In addition, the displacement of the CoM during the DSP is not intuitive and cannot be perfectly integrated with LIPM.
Shibuya et al. [29] proposed the linear pendulum model (LPM) to plan the trajectory of the CoM in the DSP, and determine the appropriate suspension point, which can ensure that the acceleration of the CoM is continuous at the moment of the switch between the SSP and the DSP. However, they only plan the cyclic gait of the robot on the horizontal ground, and do not give the gait-planning method when the robot faces a more complex environment. Shibuya et al. [30] extended the results of [29] to generate DSP trajectories in two situations. One is to land the swing leg earlier than planned, and the other is trajectory planning to stop walking in the DSP. However, they still did not put forward the method in more situations.
In our previous work, the two-point-foot walking model was proposed. A planar walking pattern was designed in [31]. Then this method is extended to 3D biped walking in [32]. Furthermore, it was implemented to the speed control for dynamic walking [33]. None of these works have fully considered the role of the DSP. In this study, we extend our works and add a DSP in gait planning.
In order to meet the requirement of real-time trajectory generation in complex environments, the gait-planning method should satisfy either planning simplicity or walking stability. In this paper, LIPM and LPM are used to plan the trajectories of the SSP and the DSP, respectively. The dynamic equations of LIPM and LPM are linear, so they have analytic solutions. Trajectory planning only needs a small amount of computation. Through dynamic analysis of two pendulum models and their ZMP, the stability of gait can be guaranteed. Moreover, LPM is well-compatible with LIPM. Not only does the trajectory of CoM have an analytical solution, but the displacement of CoM in the SSP and DSP also has a very intuitive geometric representation.
The novelty of the paper is to propose a trajectory-planning method using LIPM and LPM, which can achieve flexible walking for biped robots. According to our proposed method, biped robots can generate trajectory online for several cases, such as periodic walking, changing walking speed, disturbance recovery, and walking on uneven ground, etc. The main difference between this paper and other papers is the trajectory-planning method in the double support phase. In some papers, the trajectory of the CoM is planned by LIPM in the SSP and is designed as a polynomial in the DSP. However, the trajectory of the CoM is planned by LIPM and LPM in the SSP and the DSP, respectively. According to the characteristics of LIPM and LPM, some geometric relations can ensure that acceleration of the CoM is continuous at the switch of SSP and DSP. Moreover, according to the regulation of the DSP, the biped robot can walk flexibly. The main work of this paper is as follows: Firstly, LIPM and LPM are introduced, their dynamic equations are analyzed, and their analytical solutions are given. Secondly, the walking stability of the robot during the SSP and DSP is analyzed. Then, several trajectory-planning methods are proposed under different situations. At last, simulations are carried out and simulation results show that the biped robot can walk stably and flexibly based on our proposed gait-planning method. This validates the effectiveness of the proposed method.

LIPM and LPM
To simplify the biped model, we consider the biped robot as a concentrated mass with massless legs. In this paper, we only consider the motion of the sagittal plane. LIPM and LPM are very simple and their dynamic equations have analytical solutions. We use LIPM to plan the CoM's trajectory for the SSP and LPM for the DSP. In this section, the dynamic equations of LIPM and LPM and their analytical solutions are discussed.

LIPM
When the CoM of the robot moves along a horizontal straight line under the force of the massless leg, we call it a linear inverted pendulum, as shown in Figure 1. Readers interested can refer to Kajita et al. [27].

Dynamic Equation with Constant Height of CoM
When the CoM moves along a horizontal straight line, its resultant force in the vertical direction is zero. Therefore, the dynamic differential equation of the CoM in the horizontal direction could be obtained by: ..
where (x, z s ) is the position of the CoM. Because z s remains constant, the solution of the ordinary differential equation can be easily obtained by: .

Dynamic Equation with Constant Height of CoM
When the CoM moves along a horizontal straight line, its resultant force in the vertical direction is zero. Therefore, the dynamic differential equation of the CoM in the horizontal direction could be obtained by: where ( , ) is the position of the CoM. Because remains constant, the solution of the ordinary differential equation can be easily obtained by: is the time constant, which is related to the height of the CoM and the acceleration of gravity. (0), (0) is the initial state.

Orbital Energy
Multiplying on both side of the dynamic Equation (1), and integrating it, we can obtain: The quantity is called orbital energy of LIPM. Equation (4) shows the conservation of orbital energy of LIPM.

Orbital Energy
Multiplying .
x on both side of the dynamic Equation (1), and integrating it, we can obtain: The quantity E s is called orbital energy of LIPM. Equation (4) shows the conservation of orbital energy of LIPM.

Transfer Time
In lots of situations, the time in which the CoM moves from one point to another is required. By doing some algebraic operations on Equations (2) and (3), we can get the following formula for calculating the transfer time: or: The results of Equations (5) and (6) are the same, unless one of them is singular by the numerator or denominator being zero.

The CoM Moving along a Constrained Line
If the CoM moves along an oblique line, as shown in Figure 2, the equation of the oblique line is: where k s is the slope of the line, z s is the intersection of the line and the Z-axis. This line is called the constraint line. If the CoM moves along an oblique line, as shown in Figure 2, the equation of the oblique line is: where is the slope of the line, is the intersection of the line and the Z-axis. This line is called the constraint line. The dynamics of the CoM for LIPM in X-direction can be obtained by: It is found that the equation of motion (8) is independent of the slope of the constraint line, and when CoM moves along the oblique line under the push force , the dynamic Equation (8) is the same as the dynamic Equation (1) in X-direction. In other words, when the initial conditions in X-direction, ( , ), are the same and the intersection points of constraint lines and Z-axis, , are also the same, the three motions shown in Figure 3 are the same in X-direction. The dynamics of the CoM for LIPM in X-direction can be obtained by: ..
It is found that the equation of motion (8) is independent of the slope of the constraint line, and when CoM moves along the oblique line under the push force f , the dynamic Equation (8) is the same as the dynamic Equation (1) in X-direction. In other words, when the initial conditions in X-direction, x 0 , .
x 0 , are the same and the intersection points of constraint lines and Z-axis, z s , are also the same, the three motions shown in Figure 3 are the same in X-direction.

LPM
During the DSP, both feet of the robot are in contact with the ground. It can be considered that there are two kick forces applied to the CoM. Assume that two kick forces are , , and their resultant force is f, as shown in Figure 4.

LPM
During the DSP, both feet of the robot are in contact with the ground. It can be considered that there are two kick forces applied to the CoM. Assume that two kick forces are f 1 , f 2 , and their resultant force is f, as shown in Figure 4.

LPM
During the DSP, both feet of the robot are in contact with the ground. It can be considered that there are two kick forces applied to the CoM. Assume that two kick forces are , , and their resultant force is f, as shown in Figure 4.   It is worth to note that is always negative. Similar to LIPM, the dynamic equation of the CoM is as follows: It seems that Equation (9) is the same as Equation (1). However, we should note that > 0, whereas < 0. To highlight < 0, Equation (9) can also be written as: The solution of the differential Equation (9) or (10) is:

Dynamic Equation with a Constant Height of CoM
The local coordinate system is established with the virtual suspension point as the origin, and (x, z d ) is the position of the CoM with respect to the virtual suspending point O. It is worth to note that z d is always negative. Similar to LIPM, the dynamic equation of the CoM is as follows: ..
It seems that Equation (9) is the same as Equation (1). However, we should note that z s > 0, whereas z d < 0. To highlight z d < 0, Equation (9) can also be written as: ..
The solution of the differential Equation (9) or (10) is: .

Orbital Energy
Similar to LIPM, LPM also maintains the conservation of orbital energy, multiplying both sides of (10) by .
x and then integrating it, we can obtain: The quantity E d is called the orbital energy of the LPM. Equation (13) shows the conservation of the LPM's orbital energy of the LPM.

Transfer Time
Sometimes, the time duration in which the CoM of the LPM moves from one point to another is required. Due to the periodic motion of the LPM, there are infinite solutions to the transfer time. We only want the time in which CoM reaches the final state for the first time. By doing some algebraic operations on (11) and (12), the transfer time can be obtained by:

The CoM Moving along a Constrained Line
If the CoM of the pendulum moves along an oblique line, as shown in Figure 6, the equation of the oblique line is: where k d is the slope of the line, z d is the intersection of the line and the Z-axis. This line is called the constraint line.  (11) and (12), the transfer time can be obtained by:

The CoM Moving along a Constrained Line
If the CoM of the pendulum moves along an oblique line, as shown in Figure 6, the equation of the oblique line is: where is the slope of the line, is the intersection of the line and the Z-axis. This line is called the constraint line. The dynamics of the CoM for LPM in X-direction can be obtained by: It is found that the equation of motion (16) is also independent of the slope of the constraint line, and when the CoM moves along the constraint line under the pull force , the dynamic Equation (16) is the same as the dynamic Equation (10) in the X-direction. In other words, when the initial conditions in X-direction, ( , ), are the same and the intersection points of constraint lines and Z-axis are also the same, the three motions shown in Figure 7 are the same in X-direction. The dynamics of the CoM for LPM in X-direction can be obtained by: ..
It is found that the equation of motion (16) is also independent of the slope of the constraint line, and when the CoM moves along the constraint line under the pull force f , the dynamic Equation (16) is the same as the dynamic Equation (10) in the X-direction. In other words, when the initial conditions in X-direction, x 0 , .
x 0 , are the same and the

Walking Stability
This section discusses the walking stability of LIPM and LPM. The ZMP stability criterion is chosen in this paper. When the ZMP/CoP remains in the support polygon, it can ensure that the robot's feet will not turn over.
During the SSP, the CoM is pushed by prismatic joint. Since it is assumed that there is no torque at the pivot, CoP is the support point of LIPM, as shown in Figure 8.

Walking Stability
This section discusses the walking stability of LIPM and LPM. The ZMP stability criterion is chosen in this paper. When the ZMP/CoP remains in the support polygon, it can ensure that the robot's feet will not turn over.
During the SSP, the CoM is pushed by prismatic joint. Since it is assumed that there is no torque at the pivot, CoP is the support point O 1 of LIPM, as shown in Figure 8.

Walking Stability
This section discusses the walking stability of LIPM and LPM. The ZMP stability criterion is chosen in this paper. When the ZMP/CoP remains in the support polygon, it can ensure that the robot's feet will not turn over.
During the SSP, the CoM is pushed by prismatic joint. Since it is assumed that there is no torque at the pivot, CoP is the support point of LIPM, as shown in Figure 8. During the DSP, two forces and are applied on CoM by two legs, and their resultant force is the pull force of the prismatic joint in LPM. Suppose that the position of CoM is P at a certain moment in the DSP, as shown in Figure 8, then the CoP is the intersection point between the line and the ground, where is the virtual suspending point of LPM. When the CoM gradually moves from to , the CoP gradually moves from to , where and are the support points of LIPM for the previous During the DSP, two forces f 1 and f 2 are applied on CoM by two legs, and their resultant force is the pull force f of the prismatic joint in LPM. Suppose that the position of CoM is P at a certain moment in the DSP, as shown in Figure 8, then the CoP is the intersection point between the line S 1 P and the ground, where S 1 is the virtual suspending point of LPM. When the CoM gradually moves from P 1 to P 2 , the CoP gradually moves from O 1 to O 2 , where O 1 and O 2 are the support points of LIPM for the previous step and next step, respectively. At the beginning of the DSP, f 2 = 0 and f 1 = f , this confirms an impactless landing of the swing foot; at the end of the DSP, f 1 = 0 and f 2 = f , the rear foot of the robot can be naturally lifted off to switch the next SSP.

Trajectory-Planning Method
According to the LIPM and LPM introduced in Section 3, the trajectory of the CoM in the SSP and DSP can be planned respectively. In this section, several trajectory-planning methods are presented for periodic walking, disturbed state recovery, variable speed walking, and walking terrain-blind.

Continuity Condition of CoM Acceleration at the Transition between LIPM and LPM
If a certain geometric relationship is satisfied at the moment of the switch between the SSP and DSP, the acceleration of the CoM of LIPM and LPM at transition can be guaranteed to be equal.
Suppose the robot walks from left to right using LIPM in the SSP and LPM in the DSP, as shown in Figure 9. The support point of LIPM in the SSP of the k-th step is O k , S k is the virtual suspending point of LPM in the DSP of the k-th step, and O k+1 is the support point of LIPM in the SSP of the (k+1)-th step. When the CoM moves to B k , it changes from the SSP to DSP; when the CoM moves to A k+1 , it changes from the DSP to the SSP of the next step. If three points, O k , B k, and S k are in a straight line, according to the geometric relationship and Equations (1) and (10), the acceleration of CoM will be continuous at B k . Similarly, if three points, O k+1 , A k+1 , and S k are in a straight line, the acceleration of CoM will be continuous at A k+1 . In the following trajectory-planning methods, the continuity condition of CoM acceleration at transition is satisfied by default. step and next step, respectively. At the beginning of the DSP, = 0 and = , this confirms an impactless landing of the swing foot; at the end of the DSP, = 0 and = , the rear foot of the robot can be naturally lifted off to switch the next SSP.

Trajectory-Planning Method
According to the LIPM and LPM introduced in Section 3, the trajectory of the CoM in the SSP and DSP can be planned respectively. In this section, several trajectory-planning methods are presented for periodic walking, disturbed state recovery, variable speed walking, and walking terrain-blind.

Continuity Condition of CoM Acceleration at the Transition between LIPM and LPM
If a certain geometric relationship is satisfied at the moment of the switch between the SSP and DSP, the acceleration of the CoM of LIPM and LPM at transition can be guaranteed to be equal.
Suppose the robot walks from left to right using LIPM in the SSP and LPM in the DSP, as shown in Figure 9. The support point of LIPM in the SSP of the k-th step is , is the virtual suspending point of LPM in the DSP of the k-th step, and is the support point of LIPM in the SSP of the (k+1)-th step. When the CoM moves to , it changes from the SSP to DSP; when the CoM moves to , it changes from the DSP to the SSP of the next step. If three points, , , and are in a straight line, according to the geometric relationship and Equations (1) and (10), the acceleration of CoM will be continuous at . Similarly, if three points, , , and are in a straight line, the acceleration of CoM will be continuous at . In the following trajectory-planning methods, the continuity condition of CoM acceleration at transition is satisfied by default.

Periodic Walking
Proposition 1. When the robot walks using LIPM and LPM in the SSP and DSP, respectively, the virtual suspension point of LPM and the support point of LIPM are symmetrical, i.e., |O k S k | = |O k+1 S k | or |O k D| = |DO k+1 | as shown in Figure 9.
Proof of Theorem 1. Using the method of proof by contradiction. Suppose |O k D| = |DO k+1 |, without loss of generality, let |O k D| > |DO k+1 |. According to geometric relationship, we have |A k P k | < |P k B k | and |B k C| = |CA k+1 |. According to conservation of orbital energy of LIPM for the SSP, .
x A k < .
x B k can be obtained. According to conservation of orbital energy of LPM for the DSP, .
x B k < .
x A k+1 can be obtained. Therefore, we obtain .
x A k < .
x A k+1 . This contradicts the periodic condition .
x A k = .
x A k+1 . Once the step length is determined, the distances that the CoM travels in the SSP and DSP are related to the height parameters h s and h d : The durations of the SSP and DSP are related to the walking speeds. Figure 10a shows the CoM velocity of period walking in X-direction. Height parameter of LIPM h s = 0.9 m, height parameter of LPM h d = 0.3 m, step length L = 0.4 m, CoM velocity of LIPM is 0.3 m/s 2 when it passes over the supporting point. In Figure 10a, the blue line is the one for the SSP, the red line is the one for the DSP. The time duration of the SSP is 0.77 s, the time duration of the DSP is 0.16 s. The average speed of period walking is 0.43 m/s. energy of LIPM for the SSP, < can be obtained. According to conservation of orbital energy of LPM for the DSP, < can be obtained. Therefore, we obtain < . This contradicts the periodic condition = .□ Once the step length is determined, the distances that the CoM travels in the SSP and DSP are related to the height parameters ℎ and ℎ : The durations of the SSP and DSP are related to the walking speeds. Figure 10a shows

Adjustment of Walking Speed with Fixed Step-Length
The walking speed of the biped can be adjusted by changing the touchdown time of the swing foot. Now we quantitatively analyze the touchdown timing of the swing foot when we adjust the robot's walking speed by defining the orbital energies of LIPM of two adjacent steps as E k and E k+1 respectively. Suppose that when the CoM passes through the point P k , which is above support point O k , its velocity is v k ; when the CoM reaches the point B k , the DSP begins and the velocity is v B ; when the CoM reaches the point A k+1 , the next SSP begins, and the velocity is v A ; In the next SSP, when the CoM passes through the point P k+1 , which is above the support point O k+1 , its velocity is v k+1 , as shown in Figure 9. We have: And: According to the conservation of orbital energy of LPM during the DSP, the following equation holds: Substituting (19) and (20) into (21) yields: According to the geometric relationship, we can get the following results: Substituting (23) into (22) yields: Let: Equation (25) can be simplified as: We can obtain: Then: Once the step length L, and height parameters h s and h d are known, the position of the virtual suspension point S k can be determined. Suppose O k (0, 0), then the coordinates of S k (x Sk , z Sk ), B k (x Bk , z Bk ) and A k+1 x A k+1 , z A k+1 can be calculated as follows: Figure 10b shows the adjustment of walking speed with fixed step-length. Height parameter of LIPM h s = 0.9 m, height parameter of LPM h d = 0.3 m, step length L = 0.4 m. The current step begins at the beginning of the SSP for period walking with the apex velocity 0.3 m/s, i.e., the velocity at which the CoM passes over the supporting point of LIPM. Through the adjustment of the DSP, apex velocity of the SSP in the next step becomes 0.35 m/s. In Figure 10b the blue line is the CoM's velocity during the SSP, and the red line is the CoM's velocity during the DSP.

Trajectory Planning for Unexpected Landing Time or Foot Placement
Sometimes, the swing foot lands aground at an unexpected time or in an unexpected foot placement. In these situations, the CoM enters the DSP with an unexpected state. Trajectory of the CoM during the DSP should be adjusted in real time to retain the orbital energy in the next SSP. Figure 11 shows one of these situations, i.e., the swing foot lands on the ground early and with a smaller step length. Under the previous plan, the swing foot lands on the point O 1 when the CoM reaches the point A 1 ; the CoM moves from  adjust the CoM state in the DSP to make it return to the periodic gait in the next SSP. In Figure 10c, the blue line (for the SSP) and the pink line (for the DSP) is the CoM velocity of period walking, the red line is the adjustment of the DSP for early landing, and the green line and yellow line are the CoM velocity of the next step after adjustment. It is observed that after adjustment, the CoM state returns to the period gait. Figure 11. Trajectory planning for unexpected landing time or foot placement.
Suppose (0,0) and CoM reaches ( , ℎ ) when the swing foot lands on the ground at ( , 0), then we can calculate the coordinate of ( , ) and ( , ℎ ) as follows:

Realtime Trajectory Planning for Terrain-Blind Walking
When a biped robot walks from laboratory to the outdoor environment, the walking surface is not flat any more. If the biped is not equipped with a visual sensing system, it cannot obtain the height of a new landing point predictively. However, the height of the   Figure 11. For some reason, when CoM moves to point A 2 , the swinging foot lands the ground on point O 2 . Suppose the origin of the local coordinate frame is at the support point O, x A2 = 0.13 m, x O2 = 0.36 m (see Figure 11). Due to the early landing of the swinging foot, it is necessary to adjust the CoM state in the DSP to make it return to the periodic gait in the next SSP. In Figure 10c, the blue line (for the SSP) and the pink line (for the DSP) is the CoM velocity of period walking, the red line is the adjustment of the DSP for early landing, and the green line and yellow line are the CoM velocity of the next step after adjustment. It is observed that after adjustment, the CoM state returns to the period gait.
Suppose O(0, 0) and CoM reaches A 2 (x A2 , h s ) when the swing foot lands on the ground at O 2 (x O2 , 0), then we can calculate the coordinate of S 2 (x S2 , z S2 ) and B 2 (x B2 , h s ) as follows:

Realtime Trajectory Planning for Terrain-Blind Walking
When a biped robot walks from laboratory to the outdoor environment, the walking surface is not flat any more. If the biped is not equipped with a visual sensing system, it cannot obtain the height of a new landing point predictively. However, the height of the new landing point can be obtained by the joint sensor at the moment of the robot's foot landing on the ground. In this situation, we also need to replan the trajectory of the CoM during the DSP. As shown in Figure 12, the CoM reaches the point P 2 at the moment of the swing foot landing on point O 2 . To make the figure clearer, we enlarged the height difference between two supporting points, O 1 and O 2 . Next, we introduce how to determine the position of the suspending point.  If the CoM of LPM moves from point D to point B, it is obvious that velocities at point D and B are the same in X-direction. As discussed in Section 2.2, CoM motion from to is the same as that from point D to point B in X-direction. Therefore = , | | = | | and | | = | | ensure that the motion of the next LIPM is the same as that of the previous LIPM. As a result the position of in X-direction is at the center of points and in X-direction. Suppose (0,0) and CoM reaches ( , ) when the swing foot lands on the ground at ( , ) , then we can calculate the coordinate of ( , ) and ( , ) as follows: Firstly, determine a point, P 4 , over the support point O 2 , to make |O 2 P 4 | = h s . Secondly, determine a point, P 3 , which satisfies the height of P 3 with respect to O 2 being h s and k O 2 P 3 = −k O 1 P 2 . Here k * is the slope of the specified line. Thirdly, determine point A, which is under P 3 and keeps a height h s with respect to O 1 . Point B is the midpoint of point P 3 and point A. Point D is over P 2 and has the same height as point B. Finally, determine the virtual suspending point S 1 , which satisfies k S 1 D = k O 1 P 2 and k S 1 B = k O 2 P 3 . If the CoM of LPM moves from point D to point B, it is obvious that velocities at point D and B are the same in X-direction. As discussed in Section 2.2, CoM motion from P 2 to P 3 is the same as that from point D to point B in X-direction. Therefore .
x P 3 , |P 1 P 2 | = |P 3 P 4 | and |O 1 P 1 | = |O 2 P 4 | ensure that the motion of the next LIPM is the same as that of the previous LIPM. As a result the position of S 1 in X-direction is at the center of points O 1 and O 2 in X-direction.
Suppose O 1 (0, 0) and CoM reaches P 2 (x P2 , z P2 ) when the swing foot lands on the ground at O 2 (x O2 , z O2 ), then we can calculate the coordinate of S 1 (x S1 , z S1 ) and P 3 (x P3 , z P3 ) as follows: x S1 = x O2 When the biped walks on uneven ground, its motion in X-direction is the same as the one walking on flat ground. The difference is the motion in Z-direction. For example, if we suppose the biped robot walks on uneven ground as illustrated in Figure 12, the height difference between points O 2 and O 1 is 0.04 m, step length is 0.4 m, other parameters are h s = 0.9 m, |P 1 P 2 | = 0.15 m, and the apex velocity of LIPM is 0.3 m/s. We chose these parameters to be the same as period walking deliberately. Figure 10d shows the CoM velocity in X-direction from point P 1 to point P 4 through points P 2 and P 3 . This is the same as the period walking of Figure 10a.

Flow Chart of Proposed Methods
In this section, we propose several trajectory-planning methods for the CoM, i.e., periodic walking, adjusting walking speed, recovery of unexpected landing, and walking on uneven ground. The flowchart of overall idea is shown in Figure 13.

Flow Chart of Proposed Methods
In this section, we propose several trajectory-planning methods for the CoM, i.e., periodic walking, adjusting walking speed, recovery of unexpected landing, and walking on uneven ground. The flowchart of overall idea is shown in Figure 13.  Figure 13. Flow chart of proposed methods including periodic walking, adjusting walking speed, recovery of unexpected landing, and walking on uneven ground. In the blocks of periodic walking and adjusting walking speed, Ak, Bk, Sk and Ak+1 are shown in Figure 9. In the block of unexpected landing time or placement, A2, B2 and S2 are shown in Figure 11. In the block of uneven ground, P2, P3 and S1 are shown in Figure 12.

Simulation
In this section, simulation results of the proposed trajectory-planning method using LIPM and LPM are illustrated. To demonstrate the implementation of the research, some simulations are produced in a physical scenario.

Modeling Reality
We develop the simulation platform in Matlab Simscape. The simulation model is built according to the prototype developed by our laboratory. The robot's height is 1.34 m, and its total mass is about 40 kg. It has 23 degrees of freedom. The structure of the biped is shown in Figure 14 and the detail physical parameters is listed in Table 1. The humanoid robot is treated as a tree-structure topology with a floating base. Here, we set Figure 13. Flow chart of proposed methods including periodic walking, adjusting walking speed, recovery of unexpected landing, and walking on uneven ground. In the blocks of periodic walking and adjusting walking speed, A k , B k , S k and A k+1 are shown in Figure 9. In the block of unexpected landing time or placement, A 2 , B 2 and S 2 are shown in Figure 11. In the block of uneven ground, P 2 , P 3 and S 1 are shown in Figure 12.

Simulation
In this section, simulation results of the proposed trajectory-planning method using LIPM and LPM are illustrated. To demonstrate the implementation of the research, some simulations are produced in a physical scenario.

Modeling Reality
We develop the simulation platform in Matlab Simscape. The simulation model is built according to the prototype developed by our laboratory. The robot's height is 1.34 m, and its total mass is about 40 kg. It has 23 degrees of freedom. The structure of the biped is shown in Figure 14 and the detail physical parameters is listed in Table 1. The humanoid robot is treated as a tree-structure topology with a floating base. Here, we set the pelvis as the floating base and number it link 0, as shown in Figure 14. Herein, we will fix some degrees of freedom as needed.    In this paper, the contact model in the normal direction between the feet and ground is assumed to be a distributed and nonlinear spring-damper model [34]. A number of contact points are located under the sole like the structure of the foot shoe, as shown in Figure 15. Each contact point is composed of two nonlinear springs and two nonlinear dampers connected in parallel pattern. The contact force can be expressed as follows: , 0 and z 0 , 0 and z 0 0, 0 Here, ( , ) is the current position of the contact point in the world coordinate system; ( , ) is the initial position of the contact point; ∆ = − denotes the tangential offset of the deformation; and are the stiffness coefficient and damper coefficient in the tangential direction, respectively; ∆ = − denotes the penetration depth of the ground; and are the stiffness coefficient and damper coefficient in the normal direction, respectively; and ζ describe the influences of the penetration depth on the stiffness and damper coefficients-these two parameters have distinct physical meanings: the greater the penetration depth, the closer the contact of body surfaces, and the more the stiffness and the damper weight. Experimentally, the parameters of the contact model are shown in Table 2.  In the simulation platform, a proportion-integration-differentiation (PID) controller is used for each joint. PID parameters of each joint are listed in Table 3.   Here, (x, z) is the current position of the contact point in the world coordinate system; x grd , z grd is the initial position of the contact point; ∆x = x − x grd denotes the tangential offset of the deformation; k x and c x are the stiffness coefficient and damper coefficient in the tangential direction, respectively; ∆z = z − z grd denotes the penetration depth of the ground; k z and c z are the stiffness coefficient and damper coefficient in the normal direction, respectively; ρ and ζ describe the influences of the penetration depth on the stiffness and damper coefficients-these two parameters have distinct physical meanings: the greater the penetration depth, the closer the contact of body surfaces, and the more the stiffness and the damper weight. Experimentally, the parameters of the contact model are shown in Table 2.  In the simulation platform, a proportion-integration-differentiation (PID) controller is used for each joint. PID parameters of each joint are listed in Table 3.

Desired Trajectories in Simulations
In Section 4, we only planned the trajectories of CoM during the SSP and DSP using LIPM and LPM, respectively. In actual robot walking, we need to control the movement of each joint, so we need to know the desired angle trajectory of each joint. In this paper, we only control the lower limbs of the robot and do not consider movements of the upper limbs, such as two arms and head. Because the weight of the robot legs is relatively small, we set the total CoM of the robot at the CoM of the trunk. Hence the position of the hip joint can be calculated according to the CoM position. The robot keeps its trunk upright and its feet parallel to the ground during walking. If we know the trajectories of the hip and two ankles, we can obtain the angle trajectories of hip joint, knee joint, and ankle joint using inverse kinematics.
In order to reduce the adverse impact of collision, the velocity of the swing ankle is designed to be zero when landing. In simulations, the desired trajectory of the swing ankle is as follows: Here, S 1 is the horizontal distance between stance ankle and swing ankle when swing foot takes off from the ground, and S 2 is the horizontal distance between stance ankle and swing ankle when swing foot lands on the ground; τ = t/T ss , t ∈ [0, T ss ] and T ss is the duration of the SSP; R is the height parameter of the swing ankle.

Walking on Level Ground
In the first simulation, the biped starts to walk from standing still on a level ground. The parameters for gait generation are: L = 0.25 m, h s = 0.85 m, h d = 0.3 m. The biped stands still at the beginning. At the moment t = 1.0 s, it starts to walk and increase the walking speed gradually in the next three steps. From the fourth step, it starts cyclic walking. Then at the seventh step, its foot lands on the ground early with a step length of 0.2 m. Then it returns to the cyclic gait in the next step and goes ahead until the eleventh step. Figure 16 shows the snapshots of the simulation.

Desired Trajectories in Simulations
In Section 4, we only planned the trajectories of CoM during the SSP and DSP using LIPM and LPM, respectively. In actual robot walking, we need to control the movement of each joint, so we need to know the desired angle trajectory of each joint. In this paper, we only control the lower limbs of the robot and do not consider movements of the upper limbs, such as two arms and head. Because the weight of the robot legs is relatively small, we set the total CoM of the robot at the CoM of the trunk. Hence the position of the hip joint can be calculated according to the CoM position. The robot keeps its trunk upright and its feet parallel to the ground during walking. If we know the trajectories of the hip and two ankles, we can obtain the angle trajectories of hip joint, knee joint, and ankle joint using inverse kinematics.
In order to reduce the adverse impact of collision, the velocity of the swing ankle is designed to be zero when landing. In simulations, the desired trajectory of the swing ankle is as follows: Here, is the horizontal distance between stance ankle and swing ankle when swing foot takes off from the ground, and is the horizontal distance between stance ankle and swing ankle when swing foot lands on the ground; = ⁄ , ∈ [0, ] and is the duration of the SSP; is the height parameter of the swing ankle.

Walking on Level Ground
In the first simulation, the biped starts to walk from standing still on a level ground. The parameters for gait generation are: L = 0.25 m, ℎ = 0.85 m, ℎ = 0.3 m. The biped stands still at the beginning. At the moment t = 1.0 s, it starts to walk and increase the walking speed gradually in the next three steps. From the fourth step, it starts cyclic walking. Then at the seventh step, its foot lands on the ground early with a step length of 0.2 m. Then it returns to the cyclic gait in the next step and goes ahead until the eleventh step. Figure 16 shows the snapshots of the simulation.  Figure 17 shows the desired trajectories of CoM and two ankles in X-and Z-directions, respectively. We can observe from Figure 17a that the actual CoM position is behind the desired one. This is because of the control error. From Figure 17b, we can observe that the actual CoM height is slightly lower than the desired one. This is mainly caused by the  Figure 17 shows the desired trajectories of CoM and two ankles in X-and Z-directions, respectively. We can observe from Figure 17a that the actual CoM position is behind the desired one. This is because of the control error. From Figure 17b, we can observe that the actual CoM height is slightly lower than the desired one. This is mainly caused by the control error, for example, the stance of the knee is more bent than desired, and there is foot subsidence with the nonlinear spring-damper ground. In addition, there is a decrease in the actual CoM height at the beginning of the simulation. This is because the biped starts the simulation with two legs parallel and straight, it needs to bend its knees to adjust the CoM. Ankle positions are similar to CoM positions. Due to the gravity and control error, the height of the swing ankle is lower than desired-swing foot lands on the ground earlier than the expected time. foot subsidence with the nonlinear spring-damper ground. In addition, there is a decrease in the actual CoM height at the beginning of the simulation. This is because the biped starts the simulation with two legs parallel and straight, it needs to bend its knees to adjust the CoM. Ankle positions are similar to CoM positions. Due to the gravity and control error, the height of the swing ankle is lower than desired-swing foot lands on the ground earlier than the expected time.  Figure 18 shows the ZMP trajectory in the first simulation. It can be observed that the ZMP stays in the ZMP region. In the simulation, the biped robot walks stably without falling down. Figure 19 shows the ground reaction force (GRF) during the whole simulation. It is noted that the GRF has very large magnitude at the beginning of the simulation. This is because our model assumes that the robot is suspended with zero height before the simulation, and there is a sudden gravity effect between the foot and the ground when the simulation starts. It is observed that the tangential component of the GRF is much smaller than the normal component of the GRF. It ensures no sliding between the foot and the ground.  Figure 18 shows the ZMP trajectory in the first simulation. It can be observed that the ZMP stays in the ZMP region. In the simulation, the biped robot walks stably without falling down. Figure 19 shows the ground reaction force (GRF) during the whole simulation. It is noted that the GRF has very large magnitude at the beginning of the simulation. This is because our model assumes that the robot is suspended with zero height before the simulation, and there is a sudden gravity effect between the foot and the ground when the simulation starts. It is observed that the tangential component of the GRF is much smaller than the normal component of the GRF. It ensures no sliding between the foot and the ground.

Walking on Uneven Ground
The second simulation is to verify the terrain-blind walking ability using the proposed method. Similar to the first simulation, the biped robot stands still at the beginning, then speeds up to a cyclic walking gait. From the fifth to the eighth step, the ground is no longer flat. The height differences between adjacent two feet locations are 0.03 m, 0.02 m, 0.01 m, and 0.01 m, respectively. Figure 20 shows the snapshots of the second simulation.
Similar to the first simulation, Figure 21 shows the ZMP trajectory in the second simulation and Figure 22 shows the GRF in the second simulation. From Figure 21, the ZMP measured by the contact model almost exits within the desired support polygon. Hence,

Walking on Uneven Ground
The second simulation is to verify the terrain-blind walking ability using the proposed method. Similar to the first simulation, the biped robot stands still at the beginning, then speeds up to a cyclic walking gait. From the fifth to the eighth step, the ground is no longer flat. The height differences between adjacent two feet locations are 0.03 m, 0.02 m, 0.01 m, and 0.01 m, respectively. Figure 20 shows the snapshots of the second simulation.
Similar to the first simulation, Figure 21 shows the ZMP trajectory in the second simulation and Figure 22 shows the GRF in the second simulation. From Figure 21, the ZMP measured by the contact model almost exits within the desired support polygon. Hence, Figure 19. The ground reaction forces of the left and right feet in the first simulation.

Walking on Uneven Ground
The second simulation is to verify the terrain-blind walking ability using the proposed method. Similar to the first simulation, the biped robot stands still at the beginning, then speeds up to a cyclic walking gait. From the fifth to the eighth step, the ground is no longer flat. The height differences between adjacent two feet locations are 0.03 m, 0.02 m, 0.01 m, and 0.01 m, respectively. Figure 20 shows the snapshots of the second simulation.
ankle position is lower than the desired one, caused by control error, which makes the swing foot touch the ground a little earlier than the desired moment. The unexpected touchdown of the swing foot generates an impact between the foot and the ground. Therefore, the actual ZMP is out of the ZMP margin temporarily at the moment of touchdown. Nevertheless, the biped robot will not lose its stability if the ZMP exceeds the stable region in a very short time. This situation can also be seen in other references, such as [35][36][37]. From Figure 22, it is observed that the biped can walk stably with no sliding.   Similar to the first simulation, Figure 21 shows the ZMP trajectory in the second simulation and Figure 22 shows the GRF in the second simulation. From Figure 21, the ZMP measured by the contact model almost exits within the desired support polygon. Hence, stable walking is guaranteed. At some moments, the ZMP escapes from the stable region, which may be caused by numerical and control errors. The ZMP margin is the analytic value computed in the stage of trajectory planning. In simulation, however, the swing ankle position is lower than the desired one, caused by control error, which makes the swing foot touch the ground a little earlier than the desired moment. The unexpected touchdown of the swing foot generates an impact between the foot and the ground. Therefore, the actual ZMP is out of the ZMP margin temporarily at the moment of touchdown. Nevertheless, the biped robot will not lose its stability if the ZMP exceeds the stable region in a very short time. This situation can also be seen in other references, such as [35][36][37]. From Figure 22, it is observed that the biped can walk stably with no sliding. stable walking is guaranteed. At some moments, the ZMP escapes from the stable region, which may be caused by numerical and control errors. The ZMP margin is the analytic value computed in the stage of trajectory planning. In simulation, however, the swing ankle position is lower than the desired one, caused by control error, which makes the swing foot touch the ground a little earlier than the desired moment. The unexpected touchdown of the swing foot generates an impact between the foot and the ground. Therefore, the actual ZMP is out of the ZMP margin temporarily at the moment of touchdown. Nevertheless, the biped robot will not lose its stability if the ZMP exceeds the stable region in a very short time. This situation can also be seen in other references, such as [35][36][37]. From Figure 22, it is observed that the biped can walk stably with no sliding.

Conclusions
In this paper, we proposed a trajectory-planning method that used LIPM for the SSP and LPM for the DSP. The dynamic equations of these two models have analytical solutions, so it is very easy to plan the CoM trajectory. The walking stability of LIPM and LPM is also discussed. Through dynamics analysis, walking stability of a planned trajectory can be ensured.
Using the proposed method, the trajectory of a biped robot can be generated in real time. Some trajectory-planning methods are presented. Periodic trajectory can be generated if some walking parameters are specified, such as step length, height parameters for LIPM and LPM, and apex velocity. We can also change the walking speed of the robot with a fixed step-length by adjusting the virtual suspending point of LPM. Disturbed by interference, the swing leg of the biped maybe lands at an unexpected time or on an unexpected foot placement. Through the adjustment of the DSP, the state of the biped can return to its original planned trajectory. We also present a trajectory-planning method for walking on uneven ground with unknown height. Its horizontal motion is similar to the situation of flat ground with a different motion in the vertical direction. At last, simulations are carried out. In both simulations, the biped robot can walk stably using our proposed method. Simulations verify the effectiveness of the proposed method.

Conclusions
In this paper, we proposed a trajectory-planning method that used LIPM for the SSP and LPM for the DSP. The dynamic equations of these two models have analytical solutions, so it is very easy to plan the CoM trajectory. The walking stability of LIPM and LPM is also discussed. Through dynamics analysis, walking stability of a planned trajectory can be ensured.
Using the proposed method, the trajectory of a biped robot can be generated in real time. Some trajectory-planning methods are presented. Periodic trajectory can be generated if some walking parameters are specified, such as step length, height parameters for LIPM and LPM, and apex velocity. We can also change the walking speed of the robot with a fixed step-length by adjusting the virtual suspending point of LPM. Disturbed by interference, the swing leg of the biped maybe lands at an unexpected time or on an unexpected foot placement. Through the adjustment of the DSP, the state of the biped can return to its original planned trajectory. We also present a trajectory-planning method for walking on uneven ground with unknown height. Its horizontal motion is similar to the situation of flat ground with a different motion in the vertical direction. At last, simulations are carried out. In both simulations, the biped robot can walk stably using our proposed method. Simulations verify the effectiveness of the proposed method.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.