3-D Biped Robot Walking along Slope with Dual Length Linear Inverted Pendulum Method (DLLIPM)

Abstract

A new design method to obtain walking parameters for a three-dimensional (3D) biped walking along a slope is proposed in this paper. Most research is focused on the walking directions when climbing up or down a slope only. This paper investigates a strategy to realize biped walking along a slope. In conventional methods, the centre of mass (CoM) is moved up or down during walking in this situation. This is because the height of the pendulum is kept at the same length on the left and right legs. Thus, extra effort is required in order to bring the CoM up to higher ground. In the proposed method, a different height of pendulum is applied on the left and right legs, which is called a dual length linear inverted pendulum method (DLLIPM). When a different height of pendulum is applied, it is quite difficult to obtain symmetrical and smooth pendulum motions. Furthermore, synchronization between sagittal and lateral planes is not confirmed. Therefore, DLLIPM with a Newton Raphson algorithm is proposed to solve these problems. The walking pattern for both planes is designed systematically and synchronization between them is ensured. As a result, the maximum force fluctuation is reduced with the proposed method.

1. Introduction

Robots are able to move around on the ground with different approaches such as by wheels [1], hybrid locomotion that combines wheels and legs [2] or merely by using legs. Furthermore, there are different groups of legged robots such as hexapod [3], quadruped [4] and biped [5–17]. The less the number of legs, the more difficult it is to ensure the dynamic balance of the robots. A robot with more than two legs may have a better dynamic balance but it may need more space for motions in comparison to a biped robot. A biped robot is more preferable as a platform since it can be used in small or wide areas. Therefore, a biped robot is chosen as the platform in this paper.

The most suitable robot to explore human environments is a biped robot because it has the same structure as a human being. The basic motion to explore this environment is by walking. However, the control and dynamic balance of a walking biped is not an easy task to carry out. There are many researchers who have investigated the dynamic balance motion criterion and biped robot walking patterns such as zero moment point (ZMP) [5–7], linear inverted pendulum (LIPM) [8] and virtual supporting point (VSP) [9]. Vukabrotovic et al. have introduced the ZMP concept in order to ensure the dynamic balance of a walking biped robot [5]. When the ZMP is located inside the support polygon, the robot dynamics equilibrium is achieved and therefore the robot will not fall down or tip over.

There is also research that has investigated biped robot motions on a sloped surface [10,11] and stairs [12–14]. However, all the research relates to biped walking across the slopes or stairs. Walking across the slope means climbing up or down it. A bipedal robot should be more versatile and be able to walk in various directions on the slope such as across, diagonally and as defined in Figure 1.

Figure 1.

Definitions of walking on a slope and on a step

Walking along a slope has not been investigated yet. For that reason, the main task of this paper is to investigate biped walking along a slope. In the conventional methods, the centre of mass (CoM) is suggested to be brought up or down in order to maintain the same pendulum length at the right foot (RF) and left foot (LF). In this paper, a method known as the dual length linear inverted pendulum method (DLLIPM) is proposed where the CoM always moves in a horizontal way, as shown in Figure 2.

Figure 2.

CoM trajectory in the lateral plane for walking along a slope with LIPM and DLLIPM

The conventional method requires extra effort to move the CoM up to the higher ground. In order to avoid this situation, the CoM always moves in a horizontal motion in the proposed method. It is expected that higher ground reaction force (GRF) will occur at the landing foot in the conventional methods, especially when the CoM is moving down from higher to lower ground. On the other hand, when different heights of the CoM are assigned to the RF and LF in DLLIPM, there is a possibility of unsmooth and non-symmetrical pendulum motions. Due to these non-symmetrical pendulum motions, the GRF fluctuation will increase. Therefore, the Newton Raphson algorithm is applied in order to find the optimal parameters for the walking pattern. This will be explained further in Section 4. In this paper, three-dimensional (3D) pattern generation for walking along a slope is investigated. The sagittal and lateral effects are analysed and their synchrony during walking is ensured.

There are several advantages expected from the proposed method. One of the benefits is that, since the CoM in the DLLIPM is maintained horizontally, the GRF is reduced. Less fluctuation of the forces will ensure a more stable walking motion. Besides that, smooth motions, straight walking direction and synchronization between the lateral and sagittal plane are confirmed with the proposed method.

The rest of this paper is organized as follows. In Section 2, the robot platforms are presented. In Section 3, the concept of DLLIPM is explained. The main contribution of this paper is written in Section 4. Simulation results are presented in Section 5. Finally, conclusion and future work are discussed in Section 6.

2. Robot platforms

The proposed method introduced in this paper is applicable to any general bipedal robot platforms. However, in this paper calculations and simulations are based on a bipedal robot called MARI-3 [15]. In order to verify the proposed method, a 3D dynamic simulator known as ROCOS (RObot COntrol Simulator) is used [16]. Figure 3 shows the actual robot and its ROCOS simulation model.

Figure 3.

Biped robot MARI-3

One of the advantages of DLLIPM is that it is a general method applicable to any type of bipedal robot platform. However, it is compulsory to calculate the CoM of the robot and it must be able to be tracked during walking. This is essential because all the derivations in this paper are developed based on the dynamics of a linear inverted pendulum, which is assigned to the motion of the CoM. The calculation of the CoM can be done from the kinematics of the robot. A picture of a general biped robot with its centre of mass position is shown in Figure 4. In this paper, the CoM position is controlled.

Figure 4.

MARI-3 robot structure

There are two different coordinate systems in this paper. For the trajectory planning, the centre of the landing foot is used as the origin of the coordinate system, where both the RF and LF are repeatedly used as the reference point, one after another. Whereas the origin of the coordinate system during simulations is always located at the centre of the RF at the initial position before walking begins. This coordinate system is known as the world coordinate system.

3. Dual length linear inverted pendulum method (DLLIPM)

DLLIPM is designed based on the LIPM concept [8]. DLLIPM can be used to solve biped walking along a slope or stairs. In both situations, the problems are solved by maintaining the CoM height so that it always moves horizontally. In order to realize this, the DLLIPM is used where the lengths of the pendulum at the RF and the LF are different, unlike in the conventional methods. In conventional methods, the CoM is always in parallel with the slope surface in order to let the height of the pendulum be the same for both feet. The concepts of both conventional and proposed methods are shown in Figure 5 and 6. In these figures, Zlf means the height of the CoM measured from the middle of LF when the CoM is just above the middle of LF vertically. Zrf means the height of the CoM measured from the middle of the RF when the CoM is just above the middle of the RF vertically. In the case of Figure 5b) and 6b), Zlf is also known as Zsp and Zrf is known as Zlp. ‘sp’ and ‘lp’ stand for short pendulum and long pendulum, respectively.

Figure 5.

CoM motion of a biped robot on a slope in the lateral plane

Figure 6.

CoM motion of a biped robot on a slope in the sagittal plane

3.1 Linear inverted pendulum mode (LIPM)

In this sub-section, the equations of the position and velocity of a single linear inverted pendulum are derived [8]. The single linear inverted pendulum models for sagittal and lateral motions are shown in Figures 7 and 8, respectively. In these figures, ‘i’ and ‘f’ stand for initial and final. In Figure 7, the initial position, initial velocity, initial acceleration, final position, final velocity and final acceleration of the CoM in the sagittal plane are shown. Whereas, in Figure 8, the initial position, initial velocity, initial acceleration, final position, final velocity and final acceleration of the CoM in the lateral plane are shown. x(t) and y(t) denote the instantaneous value of the CoM trajectory in the sagittal and lateral planes, correspondingly. Furthermore, S denotes the CoM stride motion between the LF and RF. Whereas Zc denotes the height of the pendulum, from the foot to the CoM position vertically.

Figure 7.

LIPM model for the sagittal plane

Figure 8.

LIPM model for the lateral plane

The motion equation of LIPM in the sagittal plane is shown as in (1). The same equation is also applicable for the lateral plane.

$\stackrel{..}{x}\left(t\right)=\frac{g}{Z_c}x\left(t\right)$

(1)

3.1.1 Position

In this sub-section, the derivation of the equation, which is related to position in single linear inverted pendulum model, is shown. From the analytical solution of (1), the equation of the position is obtained as shown in (2).

$x\left(t\right)=\left(\frac{x_i-T_c{\stackrel{.}{x}}_i}{2}\right)e^{-\frac{t}{T_c}}+\left(\frac{x_i+T_c{\stackrel{.}{x}}_i}{2}\right)e^{\frac{t}{T_c}}$

(2)

where the value of Tc in (2) is given as in (3).

$T_c=\sqrt{\frac{Z_c}{g}}$

(3)

The position equation shown in (2) is applicable for the sagittal and lateral planes.

3.1.2 Velocity

In this sub-section, the derivation of the equation, which is related to velocity in single linear inverted pendulum model, is shown. The time derivative of (2), gives an equation of the velocity, as shown in (4).

$\stackrel{.}{x}\left(t\right)=\left(\frac{{\stackrel{.}{x}}_i-\frac{x_i}{T_c}}{2}\right)e^{-\frac{t}{T_c}}+\left(\frac{{\stackrel{.}{x}}_i+\frac{x_i}{T_c}}{2}\right)e^{\frac{t}{T_c}}$

(4)

The velocity equation shown in (4) is applicable for the sagittal and lateral planes.

3.2 Trajectory planning method of DLLIPM

DLLIPM trajectory planning consists of 11 steps, as shown in this sub-section. This trajectory planning is based on an LIPM proposed by Kajita et al. with some changes [8].

Step 1) Deciding the length of the short pendulum and long pendulum: Suppose that the biped robot walks on a flat floor at a certain CoM height. This CoM height on the flat floor denotes the length of the long pendulum, Zlp. The height of the higher to lower slope surfaces is a known value, h. The length of the short pendulum, Zsp is given by (5). One leg will be on the lower surface and one leg will be on the higher surface. The pendulum on the lower surface is denoted by the longer pendulum and the pendulum on the higher surface is denoted by the shorter pendulum. The length of Zlp in the sagittal plane is the same as in the lateral plane. Also the length of Zsp in sagittal and lateral planes is equal.

$Z_{SP}=Z_{LP}-h$

(5)

Step 2) Deciding on an odd and even number of landing points, n^th: Which foot lands on the higher ground or on the lower ground needs to be determined based on the forward moving direction of the biped robot. Also, with which foot the robot starts to walk must be determined. For example, the robot is walking along a slope as shown in Figure 5b) and 6b). The robot moves in such a way that the left foot always lands on the higher ground, the right foot always lands on the lower ground and the motion starts from the right foot. Therefore, the even and odd numbers of landing points in this case are given as in (6) and (7).

$n=0,2,4\dots \left(RF,LP\right)$

(6)

$n=1,3,5\dots \left(LF,SP\right)$

(7)

Step 3) Deciding the walking parameters such as walking cycle, strides, initial CoM positions and initial landing positions: in the conventional LIPM method, these parameters should be decided before the robot starts to walk. In DLLIPM, these parameters should also be decided in the same manner. The walking cycles for the short and long pendulum are denoted by Δtsp and Δtlp, correspondingly.

Step 4) Setting the landing points' modifications in the sagittal and lateral planes, P_x* and P_y*, initial time T and initial number n, of landing points: T:=0 and n:=0. Initially P_x*=0 and P_y*=0.

Step 5) Using LIPM equations to obtain the CoM trajectory: For the long pendulum, which means n as in (6) in this example, the CoM is moved from T to T+Δtlp by using (8) for the sagittal plane and (9) for the lateral plane. Here the value of Zc is equal to Zlp. For the short pendulum, which means n as in (7), the CoM is moved from T to T+Δtsp also by using (8) for the sagittal plane and (9) for the lateral plane. In this case, the value of Zc is equal to Zsp.

$\begin{array}{l}x\left(t\right)=\left(x_i^{\left(n\right)}-P_x^{*}\right)\cosh \left(t/T_c\right)\\ \hspace{1em} +T_c{\stackrel{.}{x}}_i^{\left(n\right)}\sinh \left(t/T_c\right)+P_x^{*}\end{array}$

(8)

$\begin{array}{l}y\left(t\right)=\left(y_i^{\left(n\right)}-P_y^{*}\right)\cosh \left(t/T_c\right)\\ \hspace{1em} +T_c{\stackrel{.}{y}}_i^{\left(n\right)}\sinh \left(t/T_c\right)+P_y^{*}\end{array}$

(9)

where $T_c=\sqrt{\frac{Z_c}{g}}$

Step 6) Increment of number of landing points and time:

$n:=n+1$

(10)

If n is an even number, T:= T+Δtlp.

If n is an odd number, T:= T+Δtsp.

Step 7) Calculation of the next landing points in the sagittal and lateral planes, P_x⁽ⁿ⁾ and P_y⁽ⁿ⁾: The next landing point for the sagittal and lateral planes is calculated by using (11) and (12), respectively.

$P_x^{(n)}=P_x^{(n-1)}+s_x^{(n)}$

(11)

$P_y^{(n)}=P_y^{(n-1)}-{\left(-1\right)}^n{s}_y^{(n)}$

(12)

where s_x⁽ⁿ⁾ and s_y⁽ⁿ⁾ mean the foot strides of the n^th step in the sagittal and lateral directions.

Step 8) Calculation of the final position and final velocity of the CoM for each single support phase: If n is an even number, the final CoM position and the final velocity for the sagittal plane are calculated as in (13) and (14), respectively. On the other hand, the final CoM position and the final velocity for the lateral plane are calculated as in (15) and (16), respectively.

${\bar{x}}_{LP}^{\left(n\right)}=\left(\frac{{\bar{x}}_{LP}}{{\bar{x}}_{LP}+{\bar{x}}_{SP}}\right)s_x^{\left(n+1\right)}$

(13)

${\bar{v}}_{x-LP}^{\left(n\right)}=\left(\frac{C+1}{T_{c}S}\right){\bar{x}}_{LP}^{\left(n\right)}$

(14)

${\bar{y}}_{LP}^{\left(n\right)}={\left(-1\right)}^n\left(\frac{{\bar{y}}_{LP}}{{\bar{y}}_{LP}+{\bar{y}}_{SP}}\right)s_y^{\left(n+1\right)}$

(15)

${\bar{v}}_{y-LP}^{\left(n\right)}=\left(\frac{C-1}{T_{c}S}\right){\bar{y}}_{LP}^{\left(n\right)}$

(16)

If n is an odd number, the final CoM position and the final velocity for the sagittal and lateral planes can be solved in the same way as (13)–(16) but with the respect to the short pendulum.

The values of C and S in (14) and (16) are given in (17) and (18) (suppose that the instantaneous time at the end of each single support phase is denoted by t_end).

$C=\cosh \left(\frac{t_{end}}{T_c}\right)$

(17)

$S=\sinh \left(\frac{t_{end}}{T_c}\right)$

(18)

where the values of T_c in (14), (16), (17) and (18) are the same values used in Step 5.

Step 9) Calculation of the desired final position and desired final velocity of the CoM in the sagittal plane, x^d and ẋ^d, and in the lateral plane, y^d and ẏ^d, for each single support phase: If n is an even number, the calculation of the desired final position of the CoM in the sagittal plane is calculated from (11) and (13), as shown in (19). The calculation of the desired final position of the CoM in the lateral plane is calculated from (12) and (15), as shown in (20). Furthermore, the desired final velocity in the sagittal and lateral planes is calculated using (14) and (16), as shown in (21) and (22), correspondingly.

$x^d=P_x^{\left(n\right)}+{\bar{x}}_{LP}^{\left(n\right)}$

(19)

$y^d=P_y^{\left(n\right)}+{\bar{y}}_{LP}^{\left(n\right)}$

(20)

${\stackrel{.}{x}}^d={\bar{v}}_{x-LP}^{\left(n\right)}$

(21)

${\stackrel{.}{y}}^d={\bar{v}}_{y-LP}^{\left(n\right)}$

(22)

If n is an odd number, the desired final position and desired final velocity of the CoM for the sagittal and lateral planes can be solved in the same way as (19)–(22) but with the respect to the short pendulum.

Step 10) Calculation of foot landing modification: The foot landing modification in the sagittal plane is calculated by using (23) [8].

$\begin{array}{l}{P}_x^{*}=-\frac{a\left(C-1\right)}{D}\left(x^d-C{x}_i^{\left(n\right)}-T_{c}S{\stackrel{.}{x}}_i^{\left(n\right)}\right)\\ \hspace{1em} -\frac{bS}{T_{c}D}\left({\stackrel{.}{x}}^d-\frac{S}{T_c}{x}_i^{\left(n\right)}-C{\stackrel{.}{x}}_i^{\left(n\right)}\right)\end{array}$

(23)

where D = a(C-1)² + b(S/T_c)².

a and b are the weight parameters used in order to minimize the margin of error for position and velocity. The values of T_c, C and S in (23) are the same values as in Step 8. The foot landing modification in the lateral plane can be solved in the same manner as in (23) but with all the values from the lateral plane.

Step 11) Go back to Step 5: This step is repeated until all of the landing points, n, are completed.

4. DLLIPM with Newton Raphson (proposed method)

It is very important for biped robots to have a symmetrical and smooth motion during walking to reduce the GRF. In order to achieve this, the position, velocity and acceleration between single support phases must be connected smoothly. Furthermore, in 3D walking pattern generation, the sagittal and lateral planes must be synchronized properly. Thus, the use of the Newton Raphson method is proposed and will be presented in this section.

4.1 Biped modelling by DLLIPM for 3-D walking along a slope

For modelling, the long pendulum and short pendulum phases are represented by Phase 1 and Phase 2 as shown in Figure 9 and 10. Phase 1 is defined as the walking cycle of the biped robot when the robot is supported by the foot on the lower slope surface, which means the height of the pendulum in this cycle is referred to as the long pendulum (LP). On the other hand, Phase 2 is defined as the walking cycle of the biped robot when the robot is supported by the foot on the higher slope surface, which means the height of pendulum in this cycle is referred to as the short pendulum (SP).

Figure 9.

DLLIPM model for walking along slope in the sagittal plane

Figure 10.

DLLIPM model for walking along slope in the lateral plane

4.1.1 Definitions of parameters

All the parameters involved in Figure 9 and 10 are explained and defined in Table 1. Where k in Table 1 is either one or two, which is referred to as Phase 1 or 2, respectively.

Table 1.

Definitions of parameters, where k is 1 or 2 which is referred as phase 1 or 2, respectively

Parameters	Definitions
ZLP	CoM height of long pendulum
ZSP	CoM height of short pendulum
Sxk, Syk	Possible maximum stride of phase k
xk(t),yk(t)	Instantaneous position of phase k
xki, yki	Initial position of phase k
ẋki, ẏki	Initial velocity of phase k
ẍki, ÿki	Initial acceleration of phase k
x̄k = xkf, ȳk = ykf	Final position of phase k
ẋkf, ẏkf	Final velocity of phase k
ẍkf, ÿkf	Final acceleration of phase k
Tsk	Instantaneous time at the end of phase k

4.1.2 Assumptions

All derivations shown later are based on the following five assumptions:

Assumption 1: The initial position and initial velocity of Phase 1 are known.

Assumption 2: The height of the higher slope surface, h, which the robot will walk on, is a known value.

Assumption 3: The CoM height of the long pendulum and the short pendulum are known.

Assumption 4: The trajectories of the long and short pendulums are symmetrical. Reference of the CoM position for each phase is the middle of each supporting foot.

Assumption 5: There are six free parameters that need to be solved which are S_x1, S_x2, S_y1, S_y2, T_s1 and T_s2. S_x1, S_x2, S_y1, S_y2 are chosen because symmetrical pendulum motion can be determined when these parameters are made to be half as shown in the constraints section later. Furthermore, T_s1 and T_s2 are chosen because these walking cycle times are also important in order to determine symmetrical pendulum motion when the initial position and velocity have been decided.

4.1.3 Constraints

In order to ensure smooth motions and symmetrical pendulums, there are several constraints that must be enforced. In total, there are six constraints, which result in six equations, as shown by (24)-(29). The constraints are listed as follows:

• 1)

  By using the position equations, the pendulums for Phase 1 and 2 are made symmetrical by the following equations:

  $x_{1f}=+\frac{S_{x1}}{2}$

  (24)
  $x_{2f}=+\frac{S_{x2}}{2}$

  (25)
  $y_{1f}=+\frac{S_{y1}}{2}$

  (26)
  $y_{2f}=-\frac{S_{y2}}{2}$

  (27)

• 2)

  In order to ensure a smooth trajectory of the velocity and synchronization between the sagittal and lateral planes, a constraint as in (28) is enforced.

  ${\left({\stackrel{.}{x}}_{1f}\right)}^2+{\left({\stackrel{.}{y}}_{1f}\right)}^2={\left({\stackrel{.}{x}}_{2f}\right)}^2+{\left(-{\stackrel{.}{y}}_{2f}\right)}^2$

  (28)

• 3)

  In order to ensure a smooth trajectory of the acceleration and synchronization between the sagittal and lateral planes, a constraint as in (29) is enforced.

  ${\left({\stackrel{..}{x}}_{1f}\right)}^2+{\left({\stackrel{..}{y}}_{1f}\right)}^2={\left({\stackrel{..}{x}}_{2f}\right)}^2+{\left(-{\stackrel{..}{y}}_{2f}\right)}^2$

  (29)

4.2 Solution of S_x1, S_x2, S_y1, S_y2, T_s1 and T_s2 using the Newton-Raphson method

In order to find a solution for appropriate S_x1, S_x2, S_y1, S_y2, T_s1 and T_s2 parameters, a numerical method known as Newton-Raphson is chosen. There are several steps that need to be taken, which will be discussed further in this section.

4.2.1 Simplification of equations

Before the Newton-Raphson method can be implemented, long and complex equations must be simplified. Therefore, (30)-(39) are defined.

$\alpha =e^{\frac{T_{s1}}{T_{c1}}}$

(30)

$\beta =e^{\frac{T_{s2}}{T_{c2}}}$

(31)

$A=\left(\frac{x_{1i}+T_{c1}{\stackrel{.}{x}}_{1i}}{2}\right)$

(32)

$B=\left(\frac{x_{1i}-T_{c1}{\stackrel{.}{x}}_{1i}}{2}\right)$

(33)

$C=\left(\frac{x_{2i}+T_{c2}{\stackrel{.}{x}}_{2i}}{2}\right)$

(34)

$D=\left(\frac{x_{2i}-T_{c2}{\stackrel{.}{x}}_{2i}}{2}\right)$

(35)

$E=\left(\frac{y_{1i}+T_{c1}{\stackrel{.}{y}}_{1i}}{2}\right)$

(36)

$F=\left(\frac{y_{1i}-T_{c1}{\stackrel{.}{y}}_{1i}}{2}\right)$

(37)

$G=\left(\frac{y_{2i}+T_{c2}{\stackrel{.}{y}}_{2i}}{2}\right)$

(38)

$H=\left(\frac{y_{2i}-T_{c2}{\stackrel{.}{y}}_{2i}}{2}\right)$

(39)

The previous constraint in (24) is realized using (2), as shown in (40).

$\begin{array}{l}{x}_{1f}=\left(\frac{x_{1i}-T_{c1}{\stackrel{.}{x}}_{1i}}{2}{e}^{-\frac{T_{s1}}{T_{c1}}}\right)+\\ \hspace{1em}\hspace{1em}\left(\frac{x_{1i}+T_{c1}{\stackrel{.}{x}}_{1i}}{2}{e}^{\frac{T_{s1}}{T_{c1}}}\right)=+\frac{S_{x1}}{2}\end{array}$

(40)

By substituting (30), (32) and (33) into (40), f1(q) is defined as follows.

$f_1\left(q\right)=A{\alpha }^2-\frac{S_{x1}}{2}\alpha +B$

(41)

The previous constraint in (25) is also realized using (2), as shown in (42).

$\begin{array}{l}{x}_{2f}=\left(\frac{x_{2i}-T_{c2}{\stackrel{.}{x}}_{2i}}{2}{e}^{-\frac{T_{s2}}{T_{c2}}}\right)+\\ \hspace{1em}\hspace{1em}\left(\frac{x_{2i}+T_{c2}{\stackrel{.}{x}}_{2i}}{2}{e}^{\frac{T_{s2}}{T_{c2}}}\right)=+\frac{S_{x2}}{2}\end{array}$

(42)

By substituting (31), (34) and (35) into (42), f2(q) is obtained as follows.,

$f_2\left(q\right)=C{\beta }^2-\frac{S_{x2}}{2}\beta +D$

(41)

The previous constraint in (26) is realized using (2), as shown in (44).

$\begin{array}{l}{y}_{1f}=\left(\frac{y_{1i}-T_{c1}{\stackrel{.}{y}}_{1i}}{2}{e}^{-\frac{T_{s1}}{T_{c1}}}\right)+\\ \hspace{1em}\hspace{1em}\left(\frac{y_{1i}+T_{c1}{\stackrel{.}{y}}_{1i}}{2}{e}^{\frac{T_{s1}}{T_{c1}}}\right)=+\frac{S_{y1}}{2}\end{array}$

(44)

By substituting (30), (36) and (37) into (44), f3(q) is defined as follows.

$f_3\left(q\right)=E{\alpha }^2-\frac{S_{y1}}{2}\alpha +F$

(45)

The previous constraint in (27) is also realized using (2), as shown in (46).

$\begin{array}{l}{y}_{2f}=\left(\frac{y_{2i}-T_{c2}{\stackrel{.}{y}}_{2i}}{2}{e}^{-\frac{T_{s2}}{T_{c2}}}\right)+\\ \hspace{1em}\hspace{1em}\left(\frac{y_{2i}+T_{c2}{\stackrel{.}{y}}_{2i}}{2}{e}^{\frac{T_{s2}}{T_{c2}}}\right)=-\frac{S_{y2}}{2}\end{array}$

(46)

By substituting (31), (38) and (39) into (46), f4(q) is defined as follows.

$f_4\left(q\right)=G{\beta }^2+\frac{S_{y2}}{2}\beta +H$

(47)

The previous constraint in (28) is realized using (4), (30)–(39), as shown in (48). (48) is assigned as f5(q).

$\begin{array}{l}{\left(\frac{1}{T_{c1}}\left(2A\alpha -\frac{S_{x1}}{2}\right)\right)}^2+\\ {\left(\frac{1}{T_{c1}}\left(2E\alpha -\frac{S_{y1}}{2}\right)\right)}^2=\\ {\left(\frac{1}{T_{c2}}\left(2C\beta -\frac{S_{x2}}{2}\right)\right)}^2+\\ {\left(\frac{-1}{T_{c2}}\left(2G\beta +\frac{S_{y2}}{2}\right)\right)}^2\end{array}$

(48)

Then, f5(q) is expanded, as shown in (49).

$\begin{array}{l}{f}_5\left(q\right)=\\ \frac{4A^2{\alpha }^2}{T_{c1}{}^2}-\frac{2A\alpha S_{x1}}{T_{c1}{}^2}+\frac{S_{x1}{}^2}{4T_{c1}{}^2}+\\ \frac{4E^2{\alpha }^2}{T_{c1}{}^2}-\frac{2E\alpha S_{y1}}{T_{c1}{}^2}+\frac{S_{y1}{}^2}{4T_{c1}{}^2}-\\ \frac{4C^2{\beta }^2}{T_{c2}{}^2}+\frac{2C\beta S_{x2}}{T_{c2}{}^2}-\frac{S_{x2}{}^2}{4T_{c2}{}^2}-\\ \frac{4G^2{\beta }^2}{T_{c2}{}^2}-\frac{2G\beta S_{y2}}{T_{c2}{}^2}-\frac{S_{y2}{}^2}{4T_{c2}{}^2}\end{array}$

(49)

The previous constraint in (29) is realized using a time derivative of (4), (30)–(39), as simplified in (50). (50) is assigned as f6(q).

$\begin{array}{l}{\left(\frac{S_{x1}}{2T_{c1}{}^2}\right)}^2+{\left(\frac{S_{y1}}{2T_{c1}{}^2}\right)}^2=\\ {\left(\frac{S_{x2}}{2T_{c2}{}^2}\right)}^2+{\left(\frac{S_{y2}}{2T_{c2}{}^2}\right)}^2\end{array}$

(50)

Then, f6(q) is defined by expanding (50) as follows.

$f_6\left(q\right)=\frac{S_{x1}{}^2}{T_{c1}{}^4}+\frac{S_{y1}{}^2}{T_{c1}{}^4}-\frac{S_{x2}{}^2}{T_{c2}{}^4}-\frac{S_{y2}{}^2}{T_{c2}{}^4}$

(51)

All the functions are arranged as in (52). In (52) the values of q are the variables or parameters that need to be solved, as shown in (53).

$\overrightarrow{f}\left(q\right)={\left[f_1\left(q\right) f_2\left(q\right) f_3\left(q\right) f_4\left(q\right) f_5\left(q\right) f_6\left(q\right)\right]}^T$

(52)

$\overrightarrow{q}={\left[\alpha \beta S_{x1} S_{x2} S_{y1} S_{y2}\right]}^T$

(53)

4.2.2 Newton Raphson algorithm

The Newton Raphson method involves a process of iterations. The objective is to obtain the roots from several complex functions. In this method, the roots are calculated using (54).

$q_{n+1}=q_n-\frac{f\left(q_n\right)}{f^{\text{'}}\left(q_n\right)}$

(54)

f'(q_n) in (54) is the derivative of the functions, which is obtained as in (55).

$f^{\text{'}}\left(q_n\right)={\left[\begin{array}{ccc}\frac{\partial f_1}{\partial q_1}& \cdots & \frac{\partial f_1}{\partial q_6}\\ ⋮& \ddots & ⋮\\ \frac{\partial f_6}{\partial q_1}& \cdots & \frac{\partial f_6}{\partial q_6}\end{array}\right]}^{6\times 6}$

(55)

The value of each element in the matrix above is then obtained.

5. Simulation

5.1 The known parameters and initial guess values for Newton Raphson calculations

In the Newton Raphson method there are known parameters and initial guess values for iterations that must be given. These parameters are shown in Table 2. The meaning of each parameter in Table 2 can be determined from Table 1.

Table 2.

Parameters for Newton Raphson calculations

Parameters	Values	Type
ZLP	0.5376 m	Known value
ZSP	0.4376 m	Known value
Sx 10	0.3280 m	Initial guess
x 1i	− Sx10/2 = − 0.164 m	Known value
ẋ 1i	0.7887 m/s	Initial guess
Sx 20	0.2800 m	Initial guess
x 2i	− Sx20/2 = − 0.1027 m	Known value
ẋ 2i	0.7887 m/s	Initial guess
Sy 10	0.2053 m	Initial guess
y 1i	Sy10/2 = 0.1027 m	Known value
ẏ 1i	- 0.33 m/s	Initial guess
Sy 20	0.1616 m	Initial guess
y 2i	− Sy20/2 = − 0.0808 m	Known value
ẏ 2i	0.33 m/s	Initial guess
Ts1	0.9 s	Initial guess
Ts2	0.6 s	Initial guess
Tc 1	$\sqrt{(Z_{LP}/g)}$	Known value
Tc2	$\sqrt{(Z_{SP}/g)}$	Known value
α 0	e (Ts1/Tc1)	Initial guess
β 0	e (Ts2/Tc2)	Initial guess

5.2 Results of the iterations with the Newton Raphson method

From the values used in Sub-section 5.1, iterations are done by using (54). The results are shown in Figure 10 and 11. α is obtained as 10.5198 at the 5^th iteration and β is obtained as 10.2831 at the 4^th iteration. Furthermore, S_x1, S_x2, S_y1 and S_y2 are obtained as 0.1842m, 0.2438m, 0.2841m and 0.1287m, correspondingly. During the iterations process of the Newton Raphson method, different initial settings and different numbers of iterations may give different solutions. However, the solutions for all the parameters, α, β, S_x1, S_x2, S_y1 and S_y2 must always be positive. It should not be negative because α and β represent the time, which cannot be negative. S_x1, S_x2, S_y1 and S_y2 are the distance of the centre of mass maximum trajectory, which should also always be positive.

Figure 10.

Results of α and β after several iterations by using the Newton Raphson algorithm

Figure 11.

Results of S_x1, S_x2, S_y1 and S_y2 after several iterations by using the Newton Raphson algorithm

In order to verify all the functions from the Newton Raphson solutions, Figure 12 is plotted. In this figure, all functions are obtained as ‘0’ values after all the parameters are converged. This validates that all the constraints have been fulfilled.

Figure 12.

Results of all functions from f1(q) to f6(q) by using the Newton Rahpson algorithm

5.3 Simulation of 3D walking pattern with DLLIPM (without Newton Raphson)

In this sub-section, a simulation of a biped robot walking along a slope with the height of 10cm is presented. All simulations in this paper are done for a slope, which means the floor is not horizontally flat. Therefore, the foot should be made parallel to the slope surface and realized with an orientation command [15]. The strides are decided as in Table 3.

Table 3.

Strides table

n	1	2	3	4	5	6	7	8	9	10
Sy(n)[m]	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2	0.2
Sx(n)[m]	0.3	0.3	0.3	0.3	0.3	0.3	0.3	0.3	0.3	0.3

Zlp, Δt-lp and Δt-sp are 0.5376m, 0.9s and 0.6s correspondingly. Here, the parameter settings are chosen intuitively without the Newton Raphson method. With these settings, the trajectory planning method in sub-Section 3.2 is used. The result for the sagittal plane is obtained as in Figure 13. The result for lateral plane is shown in Figure 14.

Figure 13.

Walking pattern of ten steps in the sagittal plane (without proposed method)

Figure 14.

Walking pattern of ten steps in the lateral plane (without proposed method)

As for the sagittal plane, it is noticed from Figure 13 that the CoM position trajectory shown by the blue line is accelerated and decelerated repeatedly. Thus, the impact force may increase. It is also noticed that there are multiple levels of maximum and minimum velocity for the long pendulum and short pendulum cycles. These multiple levels of velocity will create unsmooth motion during walking.

As for the lateral plane, it is observed from Figure 14 that the CoM position trajectory becomes deviated from its initial straight position. The deviated problem is more critical during the first few steps. It is also observed in this figure that the shape of the velocity trajectory, shown by the red line is not symmetrical.

5.4 Simulation of 3D walking pattern with improved DLLIPM (with Newton Raphson)

By using DLLIPM algorithm shown earlier in Sub-section 3.2, the walking pattern for the environment or situation in Sub-section 5.3 is simulated again with the new walking parameters obtained from Sub-section 5.2. The result for the sagittal and lateral planes is obtained as in Figure 15 and 16, respectively.

Figure 15.

Walking pattern of ten steps in the sagittal plane (with proposed method)

Figure 16.

Walking pattern of ten steps in the lateral plane (with proposed method)

As for the sagittal plane, it is noticed in Figure 15 that the trajectory of the CoM does not accelerate and decelerate repeatedly as much as the previous situation, shown in the Figure 13. Thus, the impact force may be reduced. Besides that, it is noticed that the minimum velocity for the long pendulum and short pendulum cycles is now at the same value, which results in a smoother walking motion.

As for the lateral plane, it is observed that the motion of the CoM does not deviate from its initial straight position as shown in Figure 16. Furthermore, it is also noticed from this figure that the problem of the unsymmetrical shape of velocity trajectory has been solved. Smooth trajectories will help a bipedal robot to achieve safer and more stable walking. From all the results illustrated, it is verified that with the proposed method, the biped walking pattern for walking along a slope has been improved.

5.5 Ground reaction force (GRF)

Simulations of biped robot walking with consideration of the lateral and sagittal planes are conducted for walking along a slope. The simulations are conducted for walking with LIPM and DLLIPM approaches, as explained earlier in Figure 5 and 6. In other words, these simulations are done for comparison between the conventional method and the proposed method. In these figures, the RF always lands on the lower slope surface. Whereas the LF always lands on the higher slope surface during walking. The GRF during these walking simulations is measured with the force sensors located at the RF and LF of the biped robot. The results are shown in Figure 17 and 18.

Figure 17.

GRF data measured from force sensors during walking simulation (without proposed method)

Figure 18.

GRF data measured from force sensors during walking simulation (proposed method)

The solid lines in Figure 17 and 18 represent forces measured at the RF of the biped robot. Furthermore, the dashed lines in these figures are the forces measured at the LF of the biped robot. It is noticed from the red circles shown in Figure 17 that the maximum force fluctuation is 415 [N] at 5.4[s] and 4.38[N] at 8.2[s]. The maximum force fluctuation is reduced to 380[N] at 5.4[s] and 410[N] at 8.3[s] with the proposed method, as shown in Figure 18. These facts demonstrate the effectiveness of the proposed method in order to reduce the GRF.

6. Conclusions and future work

In this paper, a generalized method known as DLLIPM with the Newton Raphson algorithm is proposed in order to determine parameters of 3D biped walking along a slope. By using DLLIPM, different heights of pendulum are applied at the left and right legs.

However, when different heights of pendulum are applied, difficulty occurred in obtaining symmetrical and smooth pendulum motions. Furthermore, synchronization between the sagittal and lateral planes is not confirmed. Therefore, DLLIPM with the Newton Raphson algorithm is proposed in order to solve these problems. As a result, smoother walking trajectories and more symmetrical pendulum motions are achieved. Furthermore, it is verified that the maximum force fluctuation is reduced with the proposed method. In the future, the concept of DLLIPM will be expanded for other possible motions such as diagonal walking, turning on a slope and stairs.