Structural Design and Kinematics Simulation of Hydraulic Biped Robot

: This paper presents a novel framework of hydraulic-driven biped robot. The target biped robot is a humanoid robot NWPUBR-1 with 12 degrees of freedom (DOF), and its dimensions are close to that of an average male. The joint axis adopts the modular structure design of sensor with angle measurement, and the force sensor is deployed on the side of hydraulic actuator to facilitate force and position control of the robot. Meanwhile, the ﬁnite element analysis of critical components is conducted to meet the requirements of mechanical strength in motion. Based on the screw theory, the forward kinematics and Jacobian matrix models are established, and the inverse kinematics of robot is solved by using the analytical geometric method. To achieve real-time control of the robot gait in 3D space, a three-dimensional linear inverted pendulum model (3D-LIPM) is built, and a 3D gait model is generated by combining the zero-moment point (ZMP) theory. In the virtual environment of MATLAB software, the results of programming simulation show that the biped robot can walk stably in the virtual environment, which proves the correctness of 3D-LIPM.


Introduction
In complex environments such as fields, mountains, hills, and battlefields, biped robots can actively adjust their body height to adapt to complex changing terrains and different ground conditions, and cross obstacles and ravines. Compared with the wheeled, crawler, and multi-legged robots, the bipedal robots have stronger environmental adaptability.
In recent years, humanoid robots have achieved great progress, and researchers have successively developed different types of biped robots and designed gait planning and balance control based on solid models. The WABOT-1 robot was developed by researches from Waseda University in 1973, which was the world's first biped robot that can visually recognize the position, direction, and distance of objects, and can receive commands through hearing. In 1986, Honda started to develop humanoid robots with intelligence and mobility. Since then, a series of humanoid walking mechanisms had been successively proposed. Honda developed seven biped robots which were designated E0 (Experimental Model 0) through E6 [1], and ASIMO [2,3], respectively. As robots have been developed from the "automatic machines" to the "autonomous machines" with decision-making capabilities, they can provide increasingly more practical applications in an environment where people can coexist with them. The HRP biped robot was developed by Kawada Industries from Japan, and the motor drive method MATLAB software. The conclusion and summarization are drawn in Section 6, and we also discuss possible future works.

Mechanical Design
For the biped robots with fast humanoid walking capabilities, the specific requirements in terms of kinematic structure, mass distribution, actuator, drive mechanism, sensor layout, and other aspects should be comprehensively considered. The shape of the humanoid robot is similar to that of the human body. To allow the biped robot to provide the motion range of the lower limbs of human body, the structural dimensions are based on the standard of "Human dimensions of Chinese adults" (GB/T 10000-1988), when design the hydraulically driven biped robot.
The hydraulically driven biped robot has 12 DOF, as shown in Figure 1, including 3DOF of hip joint, 1DOF of knee joint and 2DOF of ankle joint. To facilitate the calculation of the zero-moment point (ZMP), four miniature load cells (JHBM-M2 (300 kg)) are arranged on the sole. The baffle-force feedback electro-hydraulic flow control servo valve (FF102) is used as the three-position four-way reversing valve adopts, which is arranged in a distributed way to reduce the number of pipelines and the height of the robot's center of gravity, thus effectively improving the robot's walking stability.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 23 MATLAB software. The conclusion and summarization are drawn in Section 6, and we also discuss possible future works.

Mechanical Design
For the biped robots with fast humanoid walking capabilities, the specific requirements in terms of kinematic structure, mass distribution, actuator, drive mechanism, sensor layout, and other aspects should be comprehensively considered. The shape of the humanoid robot is similar to that of the human body. To allow the biped robot to provide the motion range of the lower limbs of human body, the structural dimensions are based on the standard of " Human dimensions of Chinese adults" (GB/T 10000-1988), when design the hydraulically driven biped robot.
The hydraulically driven biped robot has 12 DOF, as shown in Figure 1, including 3DOF of hip joint, 1DOF of knee joint and 2DOF of ankle joint. To facilitate the calculation of the zero-moment point (ZMP), four miniature load cells (JHBM-M2 (300 kg)) are arranged on the sole. The baffle-force feedback electro-hydraulic flow control servo valve (FF102) is used as the three-position four-way reversing valve adopts, which is arranged in a distributed way to reduce the number of pipelines and the height of the robot's center of gravity, thus effectively improving the robot's walking stability.

Modular Joint Axis Design
The degree of freedom of the robot is determined by the number of joints. By adopting the modular structure of modular joint axis, it can reduce the cost of design and production, the modules are interchangeable, and the assembly is also quick and easy, which is beneficial to the maintenance and repair of the robot. Various joint shafts of the hydraulic biped robot adopt the same structure, as shown in Figure 2. Joint link 2 is connected to the articulated shaft with set screws, which can perform a relative rotational movement with Joint link 1. The magnitude of the rotational angular displacement is measured by an incremental hollow sleeve type encoder (EB38K). The encoders are installed on both ends of the joint axis.

Modular Joint Axis Design
The degree of freedom of the robot is determined by the number of joints. By adopting the modular structure of modular joint axis, it can reduce the cost of design and production, the modules are interchangeable, and the assembly is also quick and easy, which is beneficial to the maintenance and repair of the robot. Various joint shafts of the hydraulic biped robot adopt the same structure, as shown in Figure 2. Joint link 2 is connected to the articulated shaft with set screws, which can perform a relative rotational movement with Joint link 1. The magnitude of the rotational angular displacement is measured by an incremental hollow sleeve type encoder (EB38K). The encoders are installed on both ends of the joint axis.
For the robot knee joint, the hydraulic cylinder is used as the driving device, as shown in Figure  3. To facilitate the measurement of the output force of the hydraulic cylinder, a bellows-type pull pressure sensor (JLBM-2 (1 T tension and compression)) is installed at the end of the piston rod. The two ends of the connecting rod are connected to the support lugs by joint bearings. Through the telescopic movement of the hydraulic cylinder, the joint link 1 and the joint link 2 can be rotated around the joint axis.

The Relationship between the Displacement of Hydraulic Cylinder and the Rotation Angle of Joint Axis
The drive input of the robot joint is the linear displacement of the hydraulic cylinder, and the position is detected by the rotation angle of the joint. In the robot controlling process, the relationship between the displacement of the hydraulic cylinder and the rotation angle of the joint axis needs to be calculated. The mechanism can be equivalent to a four-link mechanism, and taking the knee joint as an example, the mutual relationship is as shown in Figure 4. The relationship between the displacement of the hydraulic cylinder and the rotation angle of the joint axis is calculated as follows.

Design of Drive/Test Parts
For the robot knee joint, the hydraulic cylinder is used as the driving device, as shown in Figure 3. To facilitate the measurement of the output force of the hydraulic cylinder, a bellows-type pull pressure sensor (JLBM-2 (1 T tension and compression)) is installed at the end of the piston rod. The two ends of the connecting rod are connected to the support lugs by joint bearings. Through the telescopic movement of the hydraulic cylinder, the joint link 1 and the joint link 2 can be rotated around the joint axis.
3. To facilitate the measurement of the output force of the hydraulic cylinder, a bellows-type pull pressure sensor (JLBM-2 (1 T tension and compression)) is installed at the end of the piston rod. The two ends of the connecting rod are connected to the support lugs by joint bearings. Through the telescopic movement of the hydraulic cylinder, the joint link 1 and the joint link 2 can be rotated around the joint axis.

The Relationship between the Displacement of Hydraulic Cylinder and the Rotation Angle of Joint Axis
The drive input of the robot joint is the linear displacement of the hydraulic cylinder, and the position is detected by the rotation angle of the joint. In the robot controlling process, the relationship between the displacement of the hydraulic cylinder and the rotation angle of the joint axis needs to be calculated. The mechanism can be equivalent to a four-link mechanism, and taking the knee joint as an example, the mutual relationship is as shown in Figure 4. The relationship between the displacement of the hydraulic cylinder and the rotation angle of the joint axis is calculated as follows.

The Relationship between the Displacement of Hydraulic Cylinder and the Rotation Angle of Joint Axis
The drive input of the robot joint is the linear displacement of the hydraulic cylinder, and the position is detected by the rotation angle of the joint. In the robot controlling process, the relationship between the displacement of the hydraulic cylinder and the rotation angle of the joint axis needs to be calculated. The mechanism can be equivalent to a four-link mechanism, and taking the knee joint as an example, the mutual relationship is as shown in Figure 4. The relationship between the displacement of the hydraulic cylinder and the rotation angle of the joint axis is calculated as follows. The length of the joint link AD is determined by the expansion allowance of the hydraulic cylinder, and we can get, The relationship between the displacement of the hydraulic cylinder and the rotation angle of the joint axis is shown in Figure 5.

Finite Element Analysis of Critical Parts
When the robot is walking, deformation and fracture of the mechanism tend to occur at the weak link under the force of gravity and ground impact. The motion joint is the weakest link of the robot mechanism, so it is necessary to conduct a finite element analysis of the critical components in order to meet the application requirements. The material is 45# steel, the elastic modulus is 2.05 × 10 11 Pa, the Poisson's ratio is 0.29, the shear modulus is 800 MPa, the mass density is 7.85 × 10 3 kg/m 3 , and the  B is the knee joint, BC is the thigh, and BO is the calf, ∠CBO = θ, l AO = 65 mm, l BO = 72 mm, l BC = 385 mm, l DC = 45 mm. According to the actual structural relationship, we can obtain: ∠AOB = ∠BCD = 90 • , and further obtain, The length l AD of the joint link AD is determined by the expansion allowance of the hydraulic cylinder, and we can get, The relationship between the displacement of the hydraulic cylinder and the rotation angle of the joint axis is shown in Figure 5.  The length of the joint link AD is determined by the expansion allowance of the hydraulic cylinder, and we can get, The relationship between the displacement of the hydraulic cylinder and the rotation angle of the joint axis is shown in Figure 5.

Finite Element Analysis of Critical Parts
When the robot is walking, deformation and fracture of the mechanism tend to occur at the weak link under the force of gravity and ground impact. The motion joint is the weakest link of the robot mechanism, so it is necessary to conduct a finite element analysis of the critical components in order to meet the application requirements. The material is 45# steel, the elastic modulus is 2.05 × 10 11 Pa, the Poisson's ratio is 0.29, the shear modulus is 800 MPa, the mass density is 7.85 × 10 3 kg/m 3 , and the

Finite Element Analysis of Critical Parts
When the robot is walking, deformation and fracture of the mechanism tend to occur at the weak link under the force of gravity and ground impact. The motion joint is the weakest link of the robot mechanism, so it is necessary to conduct a finite element analysis of the critical components in order to meet the application requirements. The material is 45# steel, the elastic modulus is 2.05 × 10 11 Pa, the Poisson's ratio is 0.29, the shear modulus is 800 MPa, the mass density is 7.85 × 10 3 kg/m 3 , and the yield strength is 530 MPa. The simulation module of SolidWorks software is used for finite element analysis, the solid mesh type is selected, and the element size is 4.7 mm.
The maximum ground impact force during human walking and running is about three times the weight of the human body, and there is sufficient safety margin. The loading force is 5 kN. The analysis results are listed in Figure 6. The finite element analysis results show that the maximum stress is smaller than the yield stress, and the strength of critical components meets the requirements.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 23 yield strength is 530 MPa. The simulation module of SolidWorks software is used for finite element analysis, the solid mesh type is selected, and the element size is 4.7 mm. The maximum ground impact force during human walking and running is about three times the weight of the human body, and there is sufficient safety margin. The loading force is 5 kN. The analysis results are listed in Figure 6. The finite element analysis results show that the maximum stress is smaller than the yield stress, and the strength of critical components meets the requirements.

Kinematics Analysis
The Denavit-Hartenberg (D-H) modeling method needs to establish a coordinate system for each joint, and each coordinate system has different posture. When the robot configuration changes, the model needs to be reestablished. Therefore, the D-H method involves a complicated process, low calculation efficiency, and unclear geometric meaning [25]. To a certain extent, it's easier to use the screw theory to solve the kinematics of the humanoid robot is than to use the traditional D-H parameter method, because it does not need to establish one coordinate system for each link, and the entire system only has two coordinate systems (inertial coordinate system S and tool coordinate system T). Therefore, when analyzing complex spatial mechanisms, the screw theory can be employed to simplify and unify the description and solution of problems, which has more clear geometric meaning and higher numerical stability [26]. Based on the screw theory, the forward kinematics equation is obtained by using the product of exponential (POE) formula, and the velocity Jacobian matrix is derived.

Forward Kinematics
The coordinate system of the right leg is established in the form of the POE, as shown in Figure  7.

Kinematics Analysis
The Denavit-Hartenberg (D-H) modeling method needs to establish a coordinate system for each joint, and each coordinate system has different posture. When the robot configuration changes, the model needs to be reestablished. Therefore, the D-H method involves a complicated process, low calculation efficiency, and unclear geometric meaning [25]. To a certain extent, it's easier to use the screw theory to solve the kinematics of the humanoid robot is than to use the traditional D-H parameter method, because it does not need to establish one coordinate system for each link, and the entire system only has two coordinate systems (inertial coordinate system S and tool coordinate system T). Therefore, when analyzing complex spatial mechanisms, the screw theory can be employed to simplify and unify the description and solution of problems, which has more clear geometric meaning and higher numerical stability [26]. Based on the screw theory, the forward kinematics equation is obtained by using the product of exponential (POE) formula, and the velocity Jacobian matrix is derived.

Forward Kinematics
The coordinate system of the right leg is established in the form of the POE, as shown in Figure 7. is the rotation axis of the rigid body, is any point on the axis, and Matrix is the twist produced by pure rotation, shown as: Here, the antisymmetric matrix The rotation matrix is obtained by Rodrigues formula, namely The vector is calculated as follows: The directions of the line vectors of various joint axes are: Take one point on various joint axes as ω is the rotation axis of the rigid body, r is any point on the axis, and Matrix [V] is the twist produced by pure rotation, shown as: The rotation matrix is obtained by Rodrigues formula, namely The vector is calculated as follows: The directions of the line vectors of various joint axes are: Appl. Sci. 2020, 10, 6377 8 of 23 Take one point on various joint axes as Under the shape and position when θ = 0, the transformation matrix of the tool coordinate system {T} relative to the basic coordinate system {S} is: The twists of various joints are: It can be obtained from the POE formula that: Appl. Sci. 2020, 10, 6377 9 of 23 In the above formula, the values are: l 1 = 125 mm, l 2 = 258.5 mm, l 3 = 164 mm, l 4 = 520 mm, l 5 = 422 mm, l 6 = 60 mm. The POE formula for the left leg can be established similarly as above.

Inverse Kinematics
For the kinematics model built based on the spin theory, its inverse solution is mainly solved based on the Paden-Kahan subproblem. However, the simple Paden-Kahan sub-problem can only be adopted to solve the inverse solution of humanoid robots with low degrees of freedom (no more than three degrees of freedom) [27]. In this paper, the method of analytical geometry is used to solve the inverse kinematics of the robot.

Sagittal Inverse Kinematics
In the sagittal plane, the knee joint angle of the right leg is θ4, and the position coordinates of the hip and ankle joints are known, as shown in Figure 8.

Speed Jacobian Matrix Based on POE Method
The speed Jacobian matrix (hereinafter referred to as Jacobian) refers to the mapping matrix from the speed of the driving joint to the speed of the end effector. As the core of the robot speed analysis, the establishment of the Jacobian matrix is the basis for kinematics analysis, acceleration analysis, accuracy analysis, static/dynamic analysis, and synthesis. The traditional way to describe a robot's Jacobian matrix is to solve its forward kinematics using different equations, which involves a more complicated solving process and results in general. The screw theory can be employed to clearly From the inverse tangent function atan2, θ 4 can be obtained where, The ankle angle θ 5 is composed of θ a and θ b , θ a = a tan 2(x 3 − x 5 , z 3 − z 5 ), and θ b = a tan 2(l 4 sin θ 4 , l 5 + l 4 cos θ 4 ) After calculating the knee joint angle θ 4 and the ankle joint angle θ 5 , the next step is to calculate the hip joint angle θ 3 . During the movement of the robot, the upper body needs to maintain a vertical posture all the time, so the sum of the ankle joint angle, knee joint angle, and hip joint angle is zero, and θ 3 can be calculated by The inverse kinematics of the left leg is obtained in a similar way as for the right leg, and we get: where, θ l4 is the knee joint angle of the left leg, θ l5 is the ankle joint angle, and θ l3 is the hip joint angle.

Coronary Inverse Kinematics
In the coronal plane, as shown in Figure 9, θ 6 is the ankle joint angle of the right leg, the height h of the hip joint relative to the ankle joint is known, and the distance between the two feet is l. We can get,

Speed Jacobian Matrix Based on POE Method
The speed Jacobian matrix (hereinafter referred to as Jacobian) refers to the mapping matrix from the speed of the driving joint to the speed of the end effector. As the core of the robot speed analysis, the establishment of the Jacobian matrix is the basis for kinematics analysis, acceleration analysis, accuracy analysis, static/dynamic analysis, and synthesis. The traditional way to describe a robot's Jacobian matrix is to solve its forward kinematics using different equations, which involves a more complicated solving process and results in general. The screw theory can be employed to clearly describe the Jacobian matrix of the serial robot, which can highlight the geometric characteristics of the robot, and avoid the deficiencies of using local parameters in the differential method [28].
The screw theory and POE method can be used to easily solve the robot's Jacobian matrix. The POE formula of the robot kinematics equation is as follows: The instantaneous space velocity of the foot can be obtained.
According to the characteristics of the POE formula of the kinematics equation, we get: To ensure that the limb above the hip joint remains vertical, the angle of the knee joint in the coronal plane should not change, so the sum of the ankle joint angle and the hip joint angle is π/2. Therefore, the hip joint angle θ 2 can be calculated by the following equation, Similarly, the joint angle θ l2 and ankle joint angle θ l6 of left leg can be calculated.

Speed Jacobian Matrix Based on POE Method
The speed Jacobian matrix (hereinafter referred to as Jacobian) refers to the mapping matrix from the speed of the driving joint to the speed of the end effector. As the core of the robot speed analysis, the establishment of the Jacobian matrix is the basis for kinematics analysis, acceleration analysis, accuracy analysis, static/dynamic analysis, and synthesis. The traditional way to describe a robot's Jacobian matrix is to solve its forward kinematics using different equations, which involves a more complicated solving process and results in general. The screw theory can be employed to clearly describe the Jacobian matrix of the serial robot, which can highlight the geometric characteristics of the robot, and avoid the deficiencies of using local parameters in the differential method [28].
The screw theory and POE method can be used to easily solve the robot's Jacobian matrix.
The POE formula of the robot kinematics equation is as follows: The instantaneous space velocity of the foot can be obtained.
According to the characteristics of the POE formula of the kinematics equation, we get: The twist coordinate can be represented as: where, the matrix S T J s (θ) ∈ R 6×n is the spatial Jacobian matrix of the foot.
The operator [ · ] in square bracket indicates an antisymmetric matrix, and ∨ is its inverse operator. It maps the joint velocity vector to the corresponding twist coordinate of the foot. Let [V i ] ∈ se(3) be the unit twist, so any column of the spatial Jacobian matrix is as follows.
Convert it to twist coordinates: Then, the space Jacobian matrix of the robot can be obtained: By conducting adjoint transformation to the twist coordinate V i of the robot's i-th joint, the i-th column V i of the Jacobian matrix can be obtained.
where, r i is the position vector of a point on the lower axis of current configuration, and ω i is the unit vector of the joint axis direction of current configuration.
Here, r i (0) is the position vector of a point on the lower axis of the initial configuration.

Gait Model
Humanoid robot systems are complex three-dimensional systems. To better understand their dynamic behavior, many simplified models have been proposed. The linear inverted pendulum (LIP) model proposed by S. Kajita [29] is often used to study gait. The model assumes that the vertical acceleration of the centroid is zero, thereby defining the analytical expression of the centroid and decoupling the sagittal and coronal equations.
In three-dimensional space, we approximately simplify the biped robot into a virtual three-dimensional linear inverted pendulum model (3D-LIPM), which is composed of the robot's centroid m and the legs connected to the centroid whose mass is ignored, as shown in Figure 10. Assuming the inverted pendulum can rotate freely around the support point, the length of the leg can change according to the support force f.

Dynamic Modeling
The motion equation of the centroid is: where, r is the distance from the centroid to the support point, and x, y, and z are the projection lengths of r in the three directions of the coordinate axis. When controlling the centroid of the robot, a constraint plane needs to be added. The equation is: To make the centroid move along the plane, the acceleration of the centroid needs to be orthogonal to the normal vector of the constraint plane, so we get − Figure 10. 3D-LIPM.

Dynamic Modeling
The motion equation of the centroid is: where, r is the distance from the centroid to the support point, and x, y, and z are the projection lengths of r in the three directions of the coordinate axis. When controlling the centroid of the robot, a constraint plane needs to be added. The equation is: To make the centroid move along the plane, the acceleration of the centroid needs to be orthogonal to the normal vector of the constraint plane, so we get Solve the equation, substitute it into Equation (43), and we can get By ensuring that the supporting force f is proportional to the leg length r, the centroid can move along the constraint plane, as shown in Figure 11.
where, = is the reciprocal of the support time constant T of the gait cycle, which is related to the height of the centroid and the acceleration of gravity. The displacement and velocity in the y-direction can be obtained using the same method.

Gait Support Time
In the walking process, the time required by the centroid from the moment when foot falls on the ground to the moment when foot leaves the ground needs to be found out. Assuming that the initial conditions ( , ) and termination conditions ( , ) are known, then we can get according to Equation (47): By utilizing the following hyperbolic function relation: + − Solving the first formula of Equation (46), we can obtain: where, ω = g z c is the reciprocal of the support time constant T of the gait cycle, which is related to the height of the centroid and the acceleration of gravity.
The displacement and velocity in the y-direction can be obtained using the same method.

Gait Support Time
In the walking process, the time required by the centroid from the moment when foot falls on the ground to the moment when foot leaves the ground needs to be found out. Assuming that the initial conditions (x 0 , .

Orbital Energy
If the initial speed is insufficient, the centroid of the inverted pendulum will return without reaching the position directly above the support point; if the initial speed is sufficient, the centroid of the inverted pendulum will continue to move forward after reaching the position directly above the support point. Therefore, the centroid reaches the support point at the minimum speed, that is, the robot centroid reaches directly above the ankle joint, and the relationship between the potential energy and the centroid motion can be calculated at this position.
Multiply both sides of the first formula of Equation (46) by .
x, and obtain the following equations by integration: we can get where, E is the sum of kinetic energy and potential energy, which is called the orbital energy. The first term on the left is kinetic energy, and the second term is potential energy. In this calculation, it is assumed that unit mass has energy, so the mass of the centroid is ignored. When E ≤ 0, the centroid will return before reaching the support point. When E > 0, we can get the velocity of the centroid passing directly above the support point. .

Exchange of Support Legs
Although the motion of a linear inverted pendulum is determined only by its initial condition, we can still control its speed, because the initial condition can be modified by the ground-contact time. For walking with fixed step, the support exchange of quickened ground-contact will decelerate the walking speed, while the support exchange of a delayed ground-contact will accelerate the walking speed.
The relationship between the exchange of support legs and the movement of inverted pendulum can be calculated. As shown in Figure 11, the step size is s, the position of the centroid relative to the previous support point is x f , and the speed of the centroid during the exchange is v f . The support is exchanged instantly, so v f is both the final speed of the previous support phase and the initial speed of the new support phase.
By defining the orbital energies before and after the exchange as E 1 and E 2 respectively, we can get: when the orbital energies E 1 and E 2 were given, we can calculate the necessary exchange condition by eliminating v f from Equations (57) and (58): The speed at the exchange moment can be calculated by (57), and we get:

Generation of 3D Gait Model
Assuming that the robot walks at a constant support exchange speed, T sup represents the support period of each step. A trajectory of the centroid during the support period is shown in Figure 12, which is symmetric around the y-axis, and the support time is [0 T sup ]. The end position is (x, y), the end speed is (v x , v y ), according to the characteristics of symmetry, the initial condition along the x-axis is (−x, v x ), and the end position is x. According to the analytical solution of the linear inverted pendulum, we can get: condition by eliminating vf from Equations (57) and (58): The speed at the exchange moment can be calculated by (57), and we get:

Generation of 3D Gait Model
Assuming that the robot walks at a constant support exchange speed, Tsup represents the support period of each step. A trajectory of the centroid during the support period is shown in Figure 12, which is symmetric around the y-axis, and the support time is [0 Tsup]. The end position is ( ̅ , ), the end speed is ( ̅ , ̅ ), according to the characteristics of symmetry, the initial condition along the xaxis is (− ̅ , ̅ ), and the end position is ̅ . According to the analytical solution of the linear inverted pendulum, we can get: Solve the above equation to obtain the end velocity ̅ : Similarly, the initial condition along the y-axis is ( , − ̅ ), and the end position is . We can get the end velocity:

Walking Parameters
In real situations, such as climbing stairs or avoiding obstacles, it is often necessary to directly specify the position of the feet.
As shown in Figure 13, the position of a simple walking foot is represented by length and width of the step. For a normal person, the step length is generally within the range 500~800 mm, and the step width is within 80 ± 35 mm. The used data are listed in Table 1.
Similarly, the initial condition along the y-axis is y, −v y , and the end position is y. We can get the end velocity:

Walking Parameters
In real situations, such as climbing stairs or avoiding obstacles, it is often necessary to directly specify the position of the feet.
As shown in Figure 13, the position of a simple walking foot is represented by length and width of the step. For a normal person, the step length is generally within the range 500~800 mm, and the step width is within 80 ± 35 mm. The used data are listed in Table 1.   In Table 1, sx represents the step length in the walking direction, sy represents the step width in the lateral direction, superscript (n) represents the n-th step, and the position of foot at the n-th step can be expressed as where, ( ) , ( ) is the position of the first supporting foot, and Figure 13 shows the right foot. If the first supporting foot is the left foot, then −(−1) n needs to be replaced with +(−1) n in the formula.
The foot position at n-th step is: It is worth noting that the foot position at the n-th step is determined by the step length and step width of the (n + 1)-th step, which is necessary for proper adjustment of the foot position and walking movement. The velocity formula at the end point is:

Modification of Foot Placements
For a robot walking with a fixed gait cycle, its speed can be controlled by adjusting the foot placements. Define the adjusted foot drop position as Px, which is also the zero-torque point (ZMP) in the x-direction, and the centroid position is xc, as shown in Figure 14.   In Table 1, s x represents the step length in the walking direction, s y represents the step width in the lateral direction, superscript (n) represents the n-th step, and the position of foot at the n-th step can be expressed as (p (n) is the position of the first supporting foot, and Figure 13 shows the right foot. If the first supporting foot is the left foot, then −(−1) n needs to be replaced with +(−1) n in the formula.
The foot position at n-th step is: It is worth noting that the foot position at the n-th step is determined by the step length and step width of the (n + 1)-th step, which is necessary for proper adjustment of the foot position and walking movement. The velocity formula at the end point is:

Modification of Foot Placements
For a robot walking with a fixed gait cycle, its speed can be controlled by adjusting the foot placements. Define the adjusted foot drop position as P x , which is also the zero-torque point (ZMP) in the x-direction, and the centroid position is x c , as shown in Figure 14.
From Equation (42), the dynamic equation along the x-direction is ..
Its analytical solution is, Its analytical solution is, Equation (68) is written in the form of a state equation, and the relationship between the position and the end state of the n-th step can be obtained Here, = is the reciprocal of the support time constant T of the gait cycle. x ( ) and x ( ) are the initial position and velocity of the n-th step, which are called the initial conditions. As the target state, it can be represented by the end state in the ground coordinate system.
The calculated final position of the foot is close to the target position, and the evaluation function defined below is adopted, The calculated final position of the foot is close to the target position, and the evaluation function defined below is adopted, where, a and b are positive. By substituting (69) into the evaluation function, let ∂N/∂p x = 0, and we can get the position of the foot when N is minimum.
where, it is assumed that C = cosh(ωt), and S = sinh(ωt). The method of walking pattern generation can be organized as an algorithm shown as follows.
Step 1. Set support period T sup and walk parameters s x , s y . Set initial position (x, y) of CoM and initial foot placement p x , p y = (p (0) y ).
Step 3. Integrate the linear inverted pendulum Equation (67) (and the equation for the y-axis) from T to T + T sup .
Step 7. Calculate the target state (x d , . x d ) using (70), and calculate the target state (y d , . y d ) using corresponding equation.
Step 8. Calculate the modified foot position p x , p y by (72) (as well as y component).

Simulation Analysis
Based on the above model and algorithm analysis, simulation was conducted in the MATLAB software through programming. The simulated virtual environment is flat ground without any barrier, and the simplified model has the weight of 27.13 kg, the height of 1.42 m and the centroid height of 0.76 m. In the simulation, the following parameters are set: the step size is 0.5 m, the stride width is 0.25 m, the foot swing height is 0.05 m, the static friction coefficient between foot bottom and ground is 1, and the dynamic friction coefficient is 0.8. The generated centroid and the feet trajectory are shown in Figure 15. During the walking process, the angle changes of various joints are presented in Figure 16. Through simulation, it was found that the biped robot can walk steadily, as shown in Figure 17, which verifies the correctness of kinematics and 3D-LIPM and has laid the foundation for further analysis of stability.
Step 8. Calculate the modified foot position , by (72) (as well as y component).

Simulation Analysis
Based on the above model and algorithm analysis, simulation was conducted in the MATLAB software through programming. The simulated virtual environment is flat ground without any barrier, and the simplified model has the weight of 27.13 kg, the height of 1.42 m and the centroid height of 0.76 m. In the simulation, the following parameters are set: the step size is 0.5 m, the stride width is 0.25 m, the foot swing height is 0.05 m, the static friction coefficient between foot bottom and ground is 1, and the dynamic friction coefficient is 0.8. The generated centroid and the feet trajectory are shown in Figure 15. During the walking process, the angle changes of various joints are presented in Figure 16. Through simulation, it was found that the biped robot can walk steadily, as shown in Figure 17, which verifies the correctness of kinematics and 3D-LIPM and has laid the foundation for further analysis of stability.  Step 6. Set the ideal position parameter ( ̅ ( ) , ̅ ( ) ) for next step using (65) and (66).

Simulation Analysis
Based on the above model and algorithm analysis, simulation was conducted in the MATLAB software through programming. The simulated virtual environment is flat ground without any barrier, and the simplified model has the weight of 27.13 kg, the height of 1.42 m and the centroid height of 0.76 m. In the simulation, the following parameters are set: the step size is 0.5 m, the stride width is 0.25 m, the foot swing height is 0.05 m, the static friction coefficient between foot bottom and ground is 1, and the dynamic friction coefficient is 0.8. The generated centroid and the feet trajectory are shown in Figure 15. During the walking process, the angle changes of various joints are presented in Figure 16. Through simulation, it was found that the biped robot can walk steadily, as shown in Figure 17, which verifies the correctness of kinematics and 3D-LIPM and has laid the foundation for further analysis of stability.

Discussion
The hydraulic biped robot with innovative design has the characteristics of simple structure, fast motion speed and strong load capacity. The configuration of distributed hydraulic valve and valve block has effectively reduced the height of centroid and improved the walk stability of biped robot, and the position or force feedback of robot can be realized using this configuration. The shortage is that the overall weight of robot is higher. Restricted by the structural size, the motion range of mechanical joints has not completely reached the motion range of human joins at present.

Discussion
The hydraulic biped robot with innovative design has the characteristics of simple structure, fast motion speed and strong load capacity. The configuration of distributed hydraulic valve and valve block has effectively reduced the height of centroid and improved the walk stability of biped robot, and the position or force feedback of robot can be realized using this configuration. The shortage is that the overall weight of robot is higher. Restricted by the structural size, the motion range of mechanical joints has not completely reached the motion range of human joins at present.
For the forward kinematics solution of robot, if the D-H method is used to establish the coordinate system, it involves a complicated process, low efficiency, and no universality. Therefore, we employ the product of exponentials (PoE) method based on the screw theory to build the forward kinematics model of robot, which is easy to understand with clear geometric significance. When solving the inverse kinematics of robot, in order to prevent excessive dependency on the Paden-Kahan sub-problem, the analytical geometry method is adopted when the motion trail of hip joint and ankle joint is known, which can easily obtain the inverse kinematics solution. This method has provided a reference for the forward kinematics and inverse kinematics solutions of similar biped robots, which has great universality.
For the gait analysis of humanoid robot, the 3D-LIPM with broad application is used, and the ZMP theory is combined for gait planning. This method has simplified the model to a certain extent. Through simulation, the problems in the mathematical model and 3D model can be found in a timely manner, so that corresponding modification and improvement can be made. In order to demonstrate the simulation program in a physical model, and the model also needs to be further optimized. In addition, falling-prevention measures should also be taken to prevent damaging the robot.

Conclusions
We propose a new hydraulic-driven biped robot, and because it involves too many components, we only conduct finite element analysis of the critical components at modular joints under bigger force and with relatively weaker structure. The loading force is the maximum impact on the robot during walking, and certain safety factors are kept. According to the analysis results, the maximum stress on critical components is lower than their allowable stress, which satisfies the strength condition.
The forward kinematics and Jacobian matrix formulas are inferred based on the screw theory. The analytical geometry method is adopted to easily obtain the inverse kinematics solution, which has provided the mathematical basis for gait planning and simulation. This method can also be used in biped robot models with similar structures.
A three-dimensional linear inverted pendulum model (3D-LIPM) is built for gait planning in the 3D space. Based on the MATLAB software platform, programming and simulation are adopted to verify the accuracy of kinematics and 3D-LIPM. The simulation results show that this biped robot can walk steadily on flat ground, which has not only laid down the basis for further verification in a physical model, but also provided basic data to conduct gait training of biped robot based on deep reinforcement learning theory in the future.

Design of Mechanical Structure
The varus and valgus motions of ankle joint is underactuated, and in order to improve the controllability of mechanism, we will consider adding this actuation device in the next step. The angle sensor of modular joint is external, although it can satisfy the measurement requirement, it affects the appearance of structure, and the sensor might also be damaged if the robot falls, so we consider designing a built-in sensor in the next step.

Gait Model
The gait model was only simulated on flat ground, and further gait planning should also be conducted on slope, turning, stairs and uneven ground to improve the theoretical model. The 3D-LIPM