*3.1. Design of a Simulation Model of a Bipedal Robot*

The starting point for the design of the robot was the kinematic model of the robot, which can be obtained in different ways depending on the principles used [19,20]. In this case, the modeling process was dependent on the analyzed walking model LIPM and ZMP, which was adapted to the number of degrees of freedom of the robot. Since it is a mobile robot, no part of the robot is firmly attached to the ground, so we cannot talk about a mechanism, but about a kinematic chain. The kinematic chain was divided into three main parts.


Based on the movement analysis, it follows that there are two degrees of freedom in the ankle and up to three in the hip joint. The implementation of several degrees of freedom in one place—such as in a spherical kinetic pair—was realized by dividing the joint into three rotary kinematic pairs with the smallest possible arms. These rotation pairs rotate around different axes of the coordinate system. In this way, the three degrees of freedom of the spherical pair were preserved with a simpler implementation in practice. One manipulator has six degrees of freedom, just like the analyzed template. Five degrees of freedom are controlled and one degree of freedom at both ends of the chain is uncontrolled. The robot has two arms on a frame with six degrees of freedom. The additional six degrees of freedom were granted because it is a mobile robot which can freely move and rotate in space in all directions. The number of freedom degrees of the kinetic chain

$$n\_{\text{robot}} = 6 + \sum\_{i=1}^{j} n\_i = 6 + 12 = 18 \ge 6,\tag{1}$$

where *ni* represents the number of degrees of freedom of the *i*-th kinematic pair and *j* is the total number of kinematic pairs, is greater than the minimum required number for the LIPM model.

The graph of the robot in Figure 4 shows the relationships between individual kinematic pairs. In the vertices of the graph, the number of degrees of freedom in each pair is indicated. The mobility of the robot is represented by the absence of a link between the ground and the kinematic chain. Furthermore, it can be observed that the graph of the robot is acyclic; therefore, we speak of an open kinematic chain.

**Figure 4.** Robot graph based on robot analysis (left), robot graph adapted for robot construction (right).

The kinematic model of the robot was derived from the robot graph (Figure 5). Kinematic pairs were defined in space and therefore in this model we are already talking about rotational joints with one degree of freedom. Relationships between joints are now represented by arms. This is a model with a combined type of walking (passive and controlled).

**Figure 5.** Kinematic model of the robot (left), visualization of LIPM and ZMP principles on the kinematic model (right).

Denavit–Hartenberg (DH) parameters were determined in the kinematic model (Figure 6). Conventional DH parameterization uses a uniform method for marking the axes of coordinate systems in individual joints. It uses two rules:

$$\text{DH1}: \mathbf{x}\_{i} \bot z\_{i-1}; \text{ DH2}: |\mathbf{x}\_{i} \cap z\_{i-1} \lor 1. \tag{2}$$

Therefore, in DH, the representation of the homogeneous transformation contains two parameters less than the usual homogeneous transformation, which facilitates the analysis of the n-arm manipulator [21].

In the process of identifying the DH parameters, the axes of rotation were first marked depending on the types of joints. Convention states that the axis of rotation is the zaxis. Based on the positive direction of rotation, the direction of the axis was determined. For some opposing joints, the direction of the z-axis was chosen differently to preserve the natural direction of rotation. Subsequently, two different origins of zero coordinate systems were chosen, specifically for the right and left leg. In both cases, the origin was chosen at the joint that connects the leg to the frame. Unlike the rest of the centers of the coordinate systems, these were not numbered. The center of the coordinate system in the joint connecting the left leg and the frame was marked O0, and the center in the joint connecting the right leg and the frame was marked O0'. A different zero-system naming convention for walking robots led to this decision. After finding the remaining origins of the coordinate systems, the *x* and *y* axes were determined by convention according to relations (3) and (4).

$$\mathbf{x} = z\_{i-1} \times z\_0 \tag{3}$$

$$y = z \times y \tag{4}$$

In effector coordinate systems, the z-axis is generally oriented in the same way as the axis of the last arm. Thus, the two non-controlled joints at the ends of the kinematic chain will rotate along the y-axis, not along the z-axis like the controlled joints. In both cases, the origins of the effectors coordinate systems were chosen at the point created by the intersection of the axis perpendicular to the ground, passing through the center of gravity, and the straight line created when the given effector touches the ground when the robot is in the starting position. Under ideal conditions, this point is identical to the ZMP. The analysis of the DH parameters of the robot's legs is shown in Figure 6. The specific values of the DH parameters of the kinematic model of the robot are listed in Table 3.

**Figure 6.** DH parameterization of the right (R) and the left (L) legs.


**Table 3.** DH parameters of the robot kinematic model.

<sup>1</sup> perpendicular distance between axes zi−<sup>1</sup> and zi, measured along the axis xi. <sup>2</sup> the angle between zi−<sup>1</sup> and zi, measured from zi−<sup>1</sup> in a plane perpendicular to xi. <sup>3</sup> perpendicular distance between axes xi−<sup>1</sup> and xi, measured along the axis zi−1. <sup>4</sup> the angle between xi−<sup>1</sup> and xi, measured from xi−<sup>1</sup> in a plane perpendicular to zi−1.

In the DH parameterization, the centers of the coordinate systems and their relative positions were found. Despite the fact that the effectors seem to fulfill the role of joints in the kinematic chain, this phenomenon is not included in the DH parameters, since they are uncontrolled "joints" between the ground and the robot. The parameterization of the effector was developed by creating a three-dimensional model of the entire kinetic chain based on the analyzed kinematic model.

When modeling the effector, knowledge from the materials for the construction of the foot of a passive walking robot was used. In the starting position, when in contact with the ground, the effectors create support in the form of two parallel lines. Since the position of the effector itself is controlled, the length of the foot was designed to be as small as possible in order to maintain relative stability in the controlled direction of rotation. The foot width was determined to be as large as possible for the longest arch of its lower part. The resulting shape of the foot design is shown in Figure 7.

**Figure 7.** Proposal of a robot effector.

Due to the chosen ratio of the weights of the legs together with the effectors, with respect to the weight of the robot frame, the requirement for the use of the LIPM walking model for the designed three-dimensional model is fulfilled. The basic parameters of the robot are listed in Table 4.

**Table 4.** Parameters of the three-dimensional model.


The resulting dimensions from Table 4 were used in motion modeling using the inverse pendulum model. In the modeling process, the center of gravity of the robot represented the head of the pendulum in the form of a mass point *p* defined in space by default as *p* = (*x*, *y*, *z*) according to (5)–(7). This position is described in this model as a system of state variables *q* = (*θr*1, *θp*1, *z*) in which *r* represents the length of pendulum arm, or robot leg in this specific case [22].

$$\mathbf{x} = r\mathbf{S}\_p; \mathbf{S}\_p \equiv \sin \theta\_p \tag{5}$$

$$y = -rS\_{\mathbf{r}}; S\_{\mathbf{r}} \equiv \sin \theta\_{\mathbf{r}} \tag{6}$$

$$z = rD; \; D \equiv \sqrt{1 - Sr^2 - Sp^2} \tag{7}$$

*τr*, *τ<sup>p</sup>* and *f*, as vectors acting on the center of gravity, arm and head of the pendulum in Formula (8), represent the inputs to the system of equations describing the movement of the inverse pendulum in space, in the Cartesian right-handed system [22] as

$$m\begin{pmatrix} \ddot{\ddot{x}} \\ \ddot{y} \\ \ddot{z} \end{pmatrix} = \left(J^{T}\right)^{-1} \begin{pmatrix} \tau\_{r} \\ \tau\_{p} \\ f \end{pmatrix} + \begin{pmatrix} 0 \\ 0 \\ -mg \end{pmatrix},\tag{8}$$

where *m* is the mass of the center of gravity. Since the arm of the pendulum is considered as a link, the variable m belongs to the mass point *p*. The Jacobian *J* is an antisymmetric matrix dependent on the current configuration *q(t)* and for motion analysis, according to [21,22], has the form

$$J = \frac{\partial p}{\partial q} = \begin{pmatrix} 0 & r\mathbb{C}\_p & S\_p \\ -r\mathbb{C}\_r & 0 & -S\_r \\ \frac{-r\mathbb{C}\_r\mathbb{S}\_r}{D} & \frac{-r\mathbb{C}\_p\mathbb{S}\_p}{D} & D \end{pmatrix}; \mathbb{C}\_r \equiv \cos\theta\_r, \mathbb{C}\_p \equiv \cos\theta\_p \tag{9}$$

A similar Jacobian was later used in the recalculation of joint velocities of the robot in order to determine the inverse kinematics. In this case, it was used to create the dynamic equations of the pendulum, and therefore it is multiplied by the vector *(* .. *x*, .. *y*, .. *z)* from the left to create a mathematical model of the inverse pendulum in space [21,22] according to the relation

$$
\begin{pmatrix}
0 & r\mathbb{C}\_p & S\_p \\
\frac{-r\mathbb{C}\_r\mathbb{S}\_r}{D} & \frac{-r\mathbb{C}\_p\mathbb{S}\_p}{D} & D
\end{pmatrix}
\begin{pmatrix}
\ddot{\boldsymbol{x}} \\
\ddot{\boldsymbol{y}} \\
\ddot{\boldsymbol{z}}
\end{pmatrix} = \begin{pmatrix}
\tau\_r \\
\tau\_p \\
f
\end{pmatrix} - mg \begin{pmatrix}
\frac{-r\mathbb{C}\_p\mathbb{S}\_p}{D} \\
D
\end{pmatrix} \tag{10}
$$

In order to use this model to generate step trajectories, restrictions on the movement of the pendulum along the z-axis were defined. These constraints ensure that the height of the pendulum head does not change under the influence of the inputs *q*(*t*). However, this height can be changed manually. This restriction is represented in a form of help vector (*kx*, *ky*, −1) of controlled inputs of height of the robot and the intersection *zc* of the z-axis. The step width was determined as the length of the arch of the foot; i.e., the maximum possible. The parameters for the design of the pendulum are listed in Table 5 [22].

**Table 5.** Three-dimensional model parameters for the LIPM model.


Step generation took place in three stages. In the first stage, the movement of the center of gravity and the ZMP in space was identified (Figure 8). According to the limitations of the model, the height of the center of gravity was anchored. This phenomenon simulated the fact that the legs are not extended when walking as in a fully passive walking robot. A step with bent knees reduces the movement of the center of gravity along the z-axis. At this stage, the pendulum consisted of a head and one arm, which alternately played the role of the robot's left and right legs. The weight of the pendulum head is preset to 1 and the weight of the arm is assumed to be zero.

**Figure 8.** Generated trajectories of the center of mass and the position of the ZMP.

In the second stage, the foot trajectory was generated based on the mutual position of the ZMP and the movement of the center of gravity. Trajectories were generated for eight cycle runs. The analysis showed that during the first two runs of the cycle, stride length and movement were less similar than during the rest of the runs. This phenomenon can also be observed in the first two steps of the model, when the feet in the initial position are located at a spacing of 91 mm from the robot axis, given by the construction. The model reaches the predetermined step width of 144 mm and length of 20 mm only during the third run of the cycle. The trajectories were stored in the data-structure *footinfos* 1 × 8 (Figure 9), where the column number determined the number of cycle run. One data element contained a 1 × 200 vector (*timevec*) with time values sampled for a one-second period. Furthermore, this element also contained the joint positions in space stored in two 6 × 200 matrices for the left leg (*footleft*) and the right leg (*footright*). The joints positions are represented by the amount of each joint rotation with respect to the defined direction using DH parameterization.


**Figure 9.** Data-structure covering the changes of DH parametrization in kinematic, 3D and simulation models.

In the third stage, the trajectories were recalculated to the controlled variables of the DH parameters using the inverse kinematics task, and the result is the motion model of the robot made from a 3D material model, which was created with Autodesk Inventor software.

#### *3.2. Simulation of a Bipedal Robot Walking*

The simulation model of the biped robot was created by exporting the three-dimensional model to the Matlab-Simscape simulation environment. The three-dimensional model was divided into 27 parts, which created links corresponding to the three-dimensional model in the simulation environment. This process was automated using the Simscape Multibody Link toolbox. Six degrees of freedom were added to the model for kinematic chain simulation. Contact forces on the lower part of the effectors were also created. These forces ensure that the bottom of the foot behaves like an ideal rigid body when in contact with another object.

The axes of rotation in the joints were adjusted to the z-axis in order to follow the convention of the DH parameterization and the positive direction of rotation was chosen according to the kinematic model. The resulting block diagram of the simulation model is shown in Figure 10.

**Figure 10.** Simulation model of bipedal robot walking.

The output of the inverse kinematic task is the required quantity in the unmodified form W'. The quantity W' must be filtered before entering the controller, because the generated trajectory is in the form of a function that the simulation solver cannot process correctly. Sharp local extremes of the trajectories were evaluated as noise and removed with a low-pass filter. The filter was implemented using the Filter Design Toolbox as a part of Matlab R2022a software, which allows the user to design the required filter in a short time even without expert experience in the field of signal processing.

From the diagram of the simulation model shown in Figure 11, it follows that the control system is based on the principle of regulating the speeds and positions of the servomotors. To increase the level of robot autonomy, an ideal orientation sensor was included in the system. The sensor detects the acceleration and position of the frame in three axes. The value on the z-axis was monitored at the output of the system. If the deviation from the planned position is exceeded by 30 mm, the robot will be turned off in

an emergency. The reason is to prevent damage to the robot due to redundant movements in an unwanted position on the ground.

**Figure 11.** Closed control loop for joint number 5.

The result of the simulation at this stage is the fact that the designed model of the robot can perform all eight planned steps from the starting position. The length and width of the first three steps contained negligible deviations from the movement model, which may be caused by using passive joints. The remaining seven steps were identical to the motion model. Figure 12 shows the simulation model in motion in different design phases from different perspectives.

**Figure 12.** Simulation of walking of the bipedal robot (**a**) basic kinematic model; (**b**) static 3D model; (**c**) walking 3D model—front view; (**d**) walking 3D model—side view.
