Abstract
Simulation is a virtual representation of a dynamic system. For the case of mechanical systems, the simulator is used to calculate the reaction forces between its base and the ground and other constraints. The intermittent nature of these forces and the mathematical inequalities that they must satisfy lead to models described by hybrid algebraic differential equations. In this paper, a simulator was developed for a seven degrees of freedom planar biped robot, which was modeled using the Euler-Lagrange formulation. This model allowed the design and implementation of a control strategy for balance management, and the monitoring of articular reference paths are tested in the simulator before proceeding to implementation on the actual prototype.
1. Introduction
Thanks to the evolution of computers during the last decades, simulation has become a tool indispensable for synthesizing control laws in mechanical systems. This tool allows to verify the behavior of a control law on a virtual model of the mechanical system known as a simulation model. It is important to highlight that this model allows to validate the hypotheses made during the controller synthesis phase. For this reason, the behavior described by the simulation model should be as close as possible to the actual behavior of the mechanical system. Bipedal robots are mechanical systems that are difficult to model and simulate due to impacts and intermittent contact between the foot and the ground. From the modeling point of view, the dynamics of a bipedal robot cannot be described in the same way as a classic mechanical system like a manipulator robot. Biped modeling requires the use of hybrid algebraic differential equations [1–3].
Bipedal robots have been studied from different approaches, like chaotic nonlinear systems [4–7]. A biped robot is modeled considering the effects of the impulses generated at the end of each step as a result of the impact of the free leg with the ground and the exchange of legs to start a new cycle; in these works, the closed-loop stability is tested via Poincare maps [8, 9]. With regard to the generation of gait trajectories, studies have been developed such as the trajectories obtained from human gait [10, 11].
The simulation model of a bipedal robot includes the calculation of the vector of generalized accelerations and the calculation of the vector of reaction forces between the foot and the ground. These forces, unlike generalized accelerations, must satisfy unilateral constraints described by mathematical inequalities: (i) normal reaction forces must be nonnegative. That means they prevent the foot from going through the ground but cannot prevent it from losing contact; (ii) the magnitude of the tangential components of the reaction forces must be less than or equal to the product between the coefficient of friction and the normal force. The simulation model of a biped is then constituted by a set of differential equations where a part of the unknowns must satisfy certain mathematical inequalities.
The main difficulty in simulating a bipedal robot comes from calculating the reaction forces between the foot and the ground. To calculate these forces, two types of methods have been proposed: the constraint methods [12] and the penalty methods [13]. For the first case, the reaction forces are determined by solving a special type of optimization problem known as a complementarity problem [14]. In the second case, models of springs and dampers are used in order to minimize the interpenetration of the bodies in contact. Both cases are based on the hypothesis of a perfectly rigid floor and nondeformable feet.
In the case of the constraint methods, there are two approaches. One called the acceleration-force approach and the other called the velocity-impulse approach. In the first case, the complementarity problem is formulated by expressing the accelerations of the contact points of the foot as a function of the reaction forces generated by the ground. For the second case, the direct dynamic model of the system is discretized using the Euler method. This makes it possible to transform the system of differential equations that governs the dynamics of the robot into a system of difference equations. Here, the complementarity problem is formulated by writing the velocities of the contact points in the instant as a function of the product between the reaction forces in and the discretization period .
The velocity-impulse approach, also known as time-stepping, was proposed in [15]. However, the most widespread formulation was presented in [16]. In the latter, the mathematical treatment is simpler; instead of differential inclusions, the problem is described using difference equations. In [16], the discretization of the dynamics of the system, it is carried out using Euler’s semi-implicit method. In [17], discretization is done using Euler’s implicit method. This modification allows to increase computational efficiency when simulating mechanical systems involving springs and dampers. As in the acceleration-force approach, in the time-stepping approach, the friction cone is approached by a polyhedron in order to get a linear problem of complementarity. In the time-stepping approach, the same complementarity formulation is employed to describe contacts and impacts regardless of whether the friction is static or dynamic. In this approach, the discretization of the dynamic model means that the simulations must be carried out with first-order fixed-step methods such as explicit and semi-implicit Euler. Nonadaptive, fixed-step simulation ensures that the accumulation of impact phenomenon cannot take place. First-order fixed-step methods, however, may require prohibitively small step sizes in order to ensure the numerical stability of the solution obtained. This, of course, may imply long simulation times. On the other hand, to compensate this deficiency, the linear problem can be solved with a low computational algorithm known as Lemke’s algorithm [18]. Current technologies for robotics simulation had their origin in the video games industry [19]. In this type of application, the computational cost plays a very important role; the avatars do not need to calculate the zero moment point (ZMP), to walk without falling. However, for bipedal robots, the simulation tools must calculate the reaction forces between the ground and the robot’s feet, in order to avoid the interpenetration of the feet with the ground and to evaluate the sliding conditions. These conditions make it necessary to include physical motors that respect the laws of physics, in each simulation step at a low computational cost. Under this premise, different simulation engines have been developed, such as ODE (Open Dynamics Engine), PhysX, Bullet, and MuJoCo. The most popular robotics simulation tools, [19] are Gazebo, V-Rep, Webots, and OpenHRP. Table 1 groups the different characteristics of the simulation tools.
Features such as open-source, multiplatform, customization, and documentation are desirable in all simulators. However, these are not found in all simulation tools; for these reasons, it was decided to program a simulator that allows evaluating the dynamics of gait in a biped robot. The rest of the paper is structured as follows: Section 2 shows the dynamic model of the robot and the constraints; Section 3 shows the discrete model by means of Euler’s method; Section 4 calculates the control law and the reference trajectories, followed by its numerical results in Section 5 and concluding the paper in Section 6.
2. Mathematical Model of the Robot
The robot’s dynamic model is given by
where is the vector of generalized coordinates. This vector contains three types of coordinates: (i) denoting the Cartesian position of the reference frame in the reference frame , (ii) denoting the orientation of the reference frame with respect to the reference frame , and (iii) denoting the six joint positions depicted in Figure 1. This model was reported in [10].
Vectors and are, respectively, the first and second derivatives of with respect to time, see the following:
is the robot’s inertia matrix. is the vector of centrifugal, gravitational, and Coriolis forces. Matrix contains only ones and zeros. The first three rows of are zero, indicating that has no direct influence on the acceleration of the first three generalized coordinates. The following 6 rows of form an identity matrix, indicating that directly affects is the vector of torques applied to the joints of the robot. This leads to fulfilling the following property:
is the Jacobian matrix relating the linear and angular velocities of the robot’s right foot in the reference frame with the vector
In Equation (4), are the linear velocities of the reference frame with respect to , whereas that is the angular velocity of the reference frame with respect to . The matrix is used to express the linear and angular velocities of the robot’s left foot to the ground (Figure 2).
In Equation (5), are the linear velocities of the reference frame with respect to , whereas that is the angular velocity of the reference frame with respect to .
The unknowns of the previous model are , , and . That is, there are 9 equations and 15 variables. The required additional equations can be obtained by using complementarity conditions related to normal and tangential foot reaction forces [20].
2.1. Complementarity Conditions
(1)Geometric Constraints. In order to prevent the foot from penetrating the ground during the simulation, the height of the vertices of the feet must be nonnegative (see Figure 2). To obtain these heights, the direct geometric model is used. This model calculates the Cartesian position of the vertices of the feet with respect to the absolute referent frame used on generalized positions defined by Equation (2).
The vector (Equation (6)) represents the Cartesian position of the vertex along the axes and from the inertial reference frame. Hereinafter, the component along the axes and of will be denoted as and , respectively. The vector with the heights of the vertices to will be denoted Equation (7). In order to ensure that the distances between the vertices of the foot and the ground are nonnegative, must satisfy the following inequality:
The equation above implies that all the 4 components of the vector must be nonnegative. (2)Inequality Constraints(i)Normal reaction forces are positive when there is contact between a foot and the ground. Then, the reaction forces prevent interpenetration between the foot and the ground but cannot stick the foot to the ground(ii)Tangential reaction forces must be maintained in the friction cone to avoid slipping ( is the Coulomb friction coefficient)
The scalars are the normal and tangential components of the reaction force between the contact point and the ground. The values are the components of , as follows: (3)Equality Constraints(i)The two corners of the stance foot must remain in contact with the ground and do not move with respect to the global reference frame. These two constraints correspond to the equations:
Vector contains the height of the two corners of the foot in contact with the ground, and were the three-element vector with the linear and the angular velocity of the foot defined in Equation (5).
Vector contains the Cartesian position of the four corners of the swing foot, and is the desired path (ii)The gait must be symmetrical with respect to the right and left legs. Hence, the final configuration of the right leg must be equal to the initial configuration of the left leg. In the case of the robot of Figure 1, this constraint implies(4)Complementarity Conditions. The height and normal reaction force of a given vertex of the foot are complementary quantities. This means that if the height is greater than zero, the normal reaction force is zero. In the same way, if the reaction force is positive, so the height is zero. Hence, the dot product between the vectors and must be zero. This condition and nonnegativity conditions on and can be written compactly using the symbol that indicates perpendicularity between vectors:
is derived twice with respect to time to make its dependency on explicit. However, the condition is equivalent to only for the contact points that satisfies and . To obtain the equivalent of Equation (14) but involving rather than , we define as the set of vertices of the foot satisfying the two aforementioned conditions, and as the vector with the normal acceleration of the elements of . Similarly, is the vector with the normal reaction forces of . The previous considerations lead to Equation (15).
Next, the complementarity equations that arise from Coulomb’s friction law are deduced.
Consider , the set of vertices of the foot in contact with the ground and with nonzero relative speed with respect to it. For the contact point , with , the complementarity condition between the tangential velocity and reaction force is [21]
Equation (16) ensures that (i) the tangential reaction force will be inside the friction cone ; (ii) if , then ; and (iii) scalars and will have opposite signs.
The solution of the complementarity problem , and is equivalent to solve Equation (17) [22, 23]:
where and are the elements of the vectors g(x) and x, respectively. The Fischer-Burmeister functional Equation (17) expresses the complementarity conditions Equations (15) and (16) as follows:
The unknowns in Equation (18) are and . In summary, the robot’s dynamic behavior is described by the differential equation (1) subject to the algebraic constraint Equation (18) like
2.2. Impact Model
Newton’s law of restitution will be used to model the impacts. The application of this law leads to discontinuities in the normal speeds of the contact points and consequently in the vector . To explain these discontinuities, the reaction forces must be assumed as Dirac distributions (the impulse forces generate impulse accelerations that when integrated produce discontinuous velocities). Around the instant of impact , the reaction force on the contact point is assumed to equal to , being the magnitude of the impulse and the Dirac distribution. The normal speed of the contact just after impact denoted and the magnitude of the impulse are subject to the following complementarity constraint (Equation (20)) [24]:
where is the set of the vertices of the foot that impact the ground in the instant . For the formulation of the impact model, it was assumed that the robot runs on completely inelastic soil. Consequently, the coefficient of restitution is assumed to equal to zero.
The tangential components of the reaction forces are also considered impulsive. The magnitude of these impulses, denoted pti, is related to tangential velocities , through
depends in turn on y de . To obtain this relationship, the direct dynamic model is integrated during an infinitesimal duration time interval around the moment of impact [25]. When making such integration is obtained like
Being in this case, is the unknown of the problem. If in the instant at least one of the vertices of the foot satisfies the conditions and , then a collision with the ground will have occurred. In such a case, the new generalized velocities qv+ are given by Equation (22).
The planar biped robot is described by Equations (19) and (22). These equations comprise different complementary formulations: one for impacts and the other two are applied depending on whether or not there is relative movement between the bodies in contact. The main drawback of this approach is the possible appearance of shock accumulation phenomena [26]. That is the occurrence of a succession of collisions whose separation interval progressively tends to zero. In the case of simulations carried out with numerical methods of variable size, this leads to the blocking of the simulation.
3. Discrete-Time Model Based on Euler’s Symplectic Method
The time-stepping approach is based on discretized models of a dynamic system. This means that the differential equation that describes the direct dynamic model must be transformed into a difference equation. For this, it is necessary to approximate the derivatives of the vectors and using a first-order difference:
y represents the values of vectors and in the instant The constant is defined as . In Euler’s symplectic method, is calculated using instead of . In the time-stepping approach, the Signorini condition is replaced by the following complementary condition:
denotes the set of vertices of the foot that satisfy the condition . The condition Equation (24) is interpreted as follows: if is equal to zero, the normal reaction force calculated at the instant must produce a positive or zero speed at the instant in order to avoid penetration of the foot into the ground. Equation (24) is valid for both, the calculation of contact and impact forces. Thus, the fundamental difference between the acceleration-force approach and the time-stepping approach is that in the latter, there is no explicit application of an impact model. Another difference between the acceleration-force approach and the time-stepping approach is that the same complementary formulation is used regardless of whether or not there is relative movement between the bodies in contact. The condition of complementarity between the tangential force in the instant and the tangential velocity in the instant is
The vectors comprises the reaction forces of the points contacting with the ground. The contact forces of the vertices that do not belong to the set are equal to zero. If the calculation of the reaction forces is formulated as an algebraic equation and combined with the discretized direct dynamic model, we obtain
The model described by Equation (26) shows the enormous simplification introduced by the discrete-time approach. The same model allows to represent robot behavior regardless of the presence or not impact and the type of contact between the vertices of the foot and the ground.
4. Control of the Biped Robot
In this section, a CTC for joint trajectory tracking is proposed. The reference articular trajectories were obtained by optimizing a performance index.
4.1. Nonlinear Control for the Single Support Phase
When the robot is resting on the left foot, the model Equation (1) is
The model Equation (27) has unknowns that are . That is, there are 9 equations and 12 variables. For the purpose of the control system, the model will be considered as subject to immobility constraints on the left foot
By deriving Equation (28) once:
By deriving Equation (29):
Combining Equations (27) and (30), a new model relating and as a function of , , and is obtained in
Rewriting Equation (30), we obtained
In Equation (33), is the vector containing the unactuated coordinates and is the actuated coordinate vector.
Matrix is composed by the first three columns of and by the last six columns. Solving of Equation (32):
Equation (34) indicates that the contact between the left foot and the ground leads to the acceleration of the coordinates not being directly driven by the acceleration of the coordinates actuated. We defined the following variables:
Substituting in the model Equation (27)
Equation (36) is multiplied by ; then
If is replaced by the desired acceleration , and is calculated by Equation (38), then Equation (37) becomes a classical CTC [27]
being the joint reference motion. By replacing in Equation (37), we obtain the following equations which allow us to calculate the joint torques required to the joint coordinates follow their respective reference values in
4.2. Optimal Reference Trajectories
The optimization approaches for the generation of gait patterns in bipedal robot permits find the articular trajectories that minimize a performance index, such as the torques exerted by the actuators or the mechanical energy expenditure [28]. The criterion used in this article was proposed in [29, 30]. It consists of minimizing Equation (40), which is related to the energy losses in the actuators of the robot
where represents the torque vector exerted by the actuators as a function of generalized positions , velocities , and accelerations . The numerical value of is calculated by using the inverse dynamic model of the robot. To transform problem defined by Equation (40) into a parametric optimization problem, the angular positions must be defined in terms of a finite set of parameters. In this way, minimizing no longer requires finding a function of time , but a vector of parameters, to find numerically. In [31], the author shows that spline interpolation as particularly simple, computationally efficient. In this approach, third-order polynomials will be used to parameterize like
Being , the vectors that contain the coefficients of the 6 polynomials. The coefficient vectors can be defined in terms of the initial and final positions and angular velocities
By solving the previous system, we obtained
where is an identity matrix of and 0 a zero matrix of . Of this way, instead of directly searching for 24 coefficients, the vectors and are sought.
4.3. Zero Moment Point
The zero moment point (ZMP) [32] is used for the control of humanoid robots. The ZMP specifies the point where the reaction forces generated by the foot contact with the ground produce no moment in the plane of contact; it is assumed that the surface of the foot is flat and has a coefficient of friction sufficient to prevent slippage. The zero moment point (ZMP) can be calculated by
where is the zero moment point, the angular momentum on axis , and the reaction force of the ground exerted on foot. To maintain the balance of robot, the condition described by Equation (45) must be satisfied.
If remains inside of the polygon support, then the robot remains in balance; otherwise, it falls. When the robot is in single support, this condition is , where is the foot length.
5. Numeric Results
Figure 3 shows the tracking of trajectories obtained by Equation (40) and the control law Equation (39).
The control law, Equation (39), the parameters like velocity, and joint acceleration measurements are necessary. To estimate the velocity and acceleration, first, the interpolation between each joint position is performed by using third-order polynomials. Then, the resulting sequence of polynomials is derived for the first time and second time to estimate the speed and acceleration, respectively. The joint velocities are shown in Figure 4. Hence, values of the desired joint coordinates and are calculated off-line.
The temporal evolution of the ZMP is shown in Figure 5. The ZMP fulfill the condition , and the control law ensures robot balance.
The tangential force is within the limits of the cone of friction for friction coefficient , like Figure 6. Then, the robot does not slip.
6. Conclusions
The simulation and control of mechanical systems subject to constraints are presented. The robot modeled is a planar biped with only six actuators and underactuated during the single support phases. The dynamic model is discretized directly from the system using Euler’s method. This allows transforming the differential equations governing the dynamics of robot in a system of difference equations. In this approach, the conditions of complementarity are obtained by writing the contact point speeds in the instant depending on the product between reaction forces at time and period discretization .
In the simulator developed the phenomenon of shock accumulation described in the introduction, it is presented when the two vertices of the same foot do not impact the ground simultaneously. In the case of the time-stepping approach, the discretization of the dynamic model implies that simulations must be performed at a fixed step. This fact ensures that the shock accumulation phenomenon may not occur. However, time-stepping methods may require step sizes prohibitively small in order to ensure the numerical stability of the solution obtained.
The trajectories were initially validated by simulation. This is crucial for the study of movement, rapid prototyping, controller design, and validation in a virtual environment before execution on a real robot. In the case of dynamic simulation, it offers visual information in the virtual world and the behavior of the robot in the real world, being the most efficient way to validate the trajectories of a gait cycle for a bipedal robot. The main difficulty in testing the dynamic simulator is that a complete mathematical model of the bipedal robot is needed, which limits its application to robots with many degrees of freedom.
In the future, I plan to prove the dynamic simulator with other biped robots with the mathematical model that was known and compare the results with commercial dynamic simulators.
Data Availability
Data of the mathematical model of the robot and the trajectories are available; please write an email to the corresponding author, Diego Bravo.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors would like to recognize and express their sincere gratitude to Universidad del Cauca (Colombia) for the financial support granted during this project (grant number 5411/2021).