# Design of a Dynamic Simulator for a Biped Robot

Authors: Diego A. Bravo M; Carlos F. Rengifo Rodas
Journal: Modelling and Simulation in Engineering
 (2021)

---

## 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.


---

## Body

## 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 instanttk+1 as a function of the product between the reaction forces in tk and the discretization period h≜tk+1−tk.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.Table 1
Features of the most popular robotics simulators.
SimulatorProgramming languagesSimulation enginesGraphic engineLicenseOperating systemsOpen sourceDocumentationGazeboPythonODEOGREGPL freeMac OS XYesGoodBulletLinuxSimbody DARTWindowsV-RepPythonODEOpenGLV-Rep Pro EduMac OS XYesGoodC/C++BulletFree licenseLinuxJava MatlabVortexWindowsWebotsPythonODECommercialMac OS XNoGoodC/C++LinuxJava MatlabWindowsOpenHRPPythonOwn engineEclipse Public LicenseWindowsYesGoodC/C++LinuxJava MatlabFeatures 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: Section2 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(1)Aqq¨+Hq,q̇=BΓ+JrTqFr+JlTqFl,whereq∈ℝ9 is the vector of generalized coordinates. This vector contains three types of coordinates: (i) x0,y0 denoting the Cartesian position of the reference frame x0,y0 in the reference frame xg,yg, (ii) q0 denoting the orientation of the reference frame x0,y0 with respect to the reference frame xg,yg, and (iii) q1,⋯,q6 denoting the six joint positions depicted in Figure 1. This model was reported in [10].Figure 1
7 DoF planar biped robot.Vectorsq̇∈ℝ9 and q¨∈ℝ9 are, respectively, the first and second derivatives of q with respect to time, see the following:
(2)q=gx0,gy0,q0,q1,⋯q6T,q̇=gẋ0,gẏ0,q̇0,q̇1,⋯q̇6T.A∈ℝ9×9 is the robot’s inertia matrix. H∈ℝ9 is the vector of centrifugal, gravitational, and Coriolis forces. Matrix B∈ℝ9×6 contains only ones and zeros. The first three rows of B are zero, indicating that Γ has no direct influence on the acceleration of the first three generalized coordinates. The following 6 rows of B form an identity matrix, indicating that Γii=1,⋯,6 directly affects q¨i,Γ∈ℝ6 is the vector of torques applied to the joints of the robot. This leads to B fulfilling the following property:
(3)BΓ=03×1Γ.Jr∈ℝ3×9 is the Jacobian matrix relating the linear and angular velocities of the robot’s right foot in the reference frame xg,yg with the vector q̇(4)gṗrxgṗryθ̇r=Jrqq̇.In Equation (4), gp˙rx,gp˙ry are the linear velocities of the reference frame x3,y3 with respect to xg,yg, whereas that θ̇r is the angular velocity of the reference frame x3,y3 with respect to xg,yg. The matrix Jl∈ℝ3×9 is used to express the linear and angular velocities of the robot’s left foot to the ground (Figure 2).
(5)gṗlxgṗlyθ̇l=Jlqq̇.Figure 2
Robot’s left foot.In Equation (5), gp˙lx,gp˙ly are the linear velocities of the reference frame x6,y6 with respect to xg,yg, whereas that θ̇l is the angular velocity of the reference frame x6,y6 with respect to xg,yg.The unknowns of the previous model areq¨∈ℝ9, Fr∈ℝ3, and Fl∈ℝ3. 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).(6)Φ=∅1⋮∅4=GqThe vectorφi∈ℝ2i=1⋯4 (Equation (6)) represents the Cartesian position of the vertex i along the axes x and y from the inertial reference frame. Hereinafter, the component along the axes x and y of φi will be denoted as φni and φti, respectively. The vector with the heights of the vertices v1 to v4 will be denoted Φn Equation (7). In order to ensure that the distances between the vertices of the foot and the ground are nonnegative, Φn must satisfy the following inequality:
(7)Φnq≥0.The equation above implies that all the 4 components of the vectorΦnq 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(8)fni≥0,i=1,⋯,4(ii)
Tangential reaction forces must be maintained in the friction cone to avoid slipping (μ is the Coulomb friction coefficient)(9)fti≤μfni,i=1,⋯,4The scalarsfni∈ℝyfti∈ℝ are the normal and tangential components of the reaction force between the contact point i and the ground. The values fniyfti are the components of FnyFt, as follows:
(10)Fn=fn1⋮fn4∈ℝ4,Ft=ft1⋮ft4∈ℝ4(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:(11)Hstn=0,Jlq̇=0.VectorHstan∈ℝ2 contains the height of the two corners of the foot in contact with the ground, and Jlq̇ were the three-element vector with the linear and the angular velocity of the foot defined in Equation (5).
(12)Pswn=Pswnd.VectorPswn∈ℝ12 contains the Cartesian position of the four corners of the swing foot, and Pswingd 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 Figure1, this constraint implies(13)q1T=−q40,q2T=−q50,q3T=−q60,q4T=−q10,q5T=−q20,q6T=−q30(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 Fn and Φn must be zero. This condition and nonnegativity conditions on Fn and Φn can be written compactly using the ⊥ symbol that indicates perpendicularity between vectors:(14)Φn≥0⊥Fn≥0Φn is derived twice with respect to time to make its dependency on Fn explicit. However, the condition Φn≥0 is equivalent to Φ¨n≥0 only for the contact points that satisfies φni=0 and φ̇ni=0. To obtain the equivalent of Equation (14) but involving Φ¨n rather than Φn, we define Nc as the set of vertices of the foot satisfying the two aforementioned conditions, and Φ¨nc as the vector with the normal acceleration of the elements of Nc. Similarly, Fnc is the vector with the normal reaction forces of Nc. The previous considerations lead to Equation (15).
(15)Φ¨nc≥0⊥Fnc≥0.Next, the complementarity equations that arise fromCoulomb’s friction law are deduced.ConsiderNs⊆Nc, 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 i, with i∈Ns, the complementarity condition between the tangential velocity and reaction force is [21]
(16)μfni−fti≥0⊥φ̇ti≥0,φ̇tifti–μfniφ̇ti=0,i∈Ns.Equation (16) ensures that (i) the tangential reaction force will be inside the friction cone fti≤μfni; (ii) if φ̇ti>0, then fti=μfni; and (iii) scalars fti and φ̇ti will have opposite signs.The solution of the complementarity problemgx≥0,x≥0, and gxTx=0 is equivalent to solve Equation (17) [22, 23]:
(17)gix2+xi2−gix+xi=0,wheregix and xi 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:
(18)Gq,q̇,Fl,Fr=0.The unknowns in Equation (18) are Fr and Fl. In summary, the robot’s dynamic behavior is described by the differential equation (1) subject to the algebraic constraint Equation (18) like
(19)Aqq¨+Hq,q̇=BΓ+JrTqFr+JlTqFl,Gq,q̇,Fl,Fr=0.


### 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 vectorq̇. 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 t=tI, the reaction force on the contact point i is assumed to equal to pniδt−tI, being pni the magnitude of the impulse and δt−tI the Dirac distribution. The normal speed of the contact i just after impact denoted φ̇ni+ and the magnitude of the impulse pni are subject to the following complementarity constraint (Equation (20)) [24]:
(20)φ̇ni++ϵφ̇ni−≥0⊥pni≥0,i∈Nc∪NI,whereNI is the set of the vertices of the foot that impact the ground in the instant tI. 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, denotedpti, is related to tangential velocities φ̇ti+, through
(21)μpni−pti≥0⊥φ̇ti+≥0,φ̇ti+pti−μpniφ̇ti+=0,i∈Nc∪NI.qv+ depends in turn on pnc y de ptc. 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
(22)qv+=qv−+A−1JrTFr+JlTFl.Being in this case,q̇v+ is the unknown of the problem. If in the instant t=tI at least one of the vertices of the foot satisfies the conditions φi=0 and φ̇i<0, 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.


## 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).(6)Φ=∅1⋮∅4=GqThe vectorφi∈ℝ2i=1⋯4 (Equation (6)) represents the Cartesian position of the vertex i along the axes x and y from the inertial reference frame. Hereinafter, the component along the axes x and y of φi will be denoted as φni and φti, respectively. The vector with the heights of the vertices v1 to v4 will be denoted Φn Equation (7). In order to ensure that the distances between the vertices of the foot and the ground are nonnegative, Φn must satisfy the following inequality:
(7)Φnq≥0.The equation above implies that all the 4 components of the vectorΦnq 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(8)fni≥0,i=1,⋯,4(ii)
Tangential reaction forces must be maintained in the friction cone to avoid slipping (μ is the Coulomb friction coefficient)(9)fti≤μfni,i=1,⋯,4The scalarsfni∈ℝyfti∈ℝ are the normal and tangential components of the reaction force between the contact point i and the ground. The values fniyfti are the components of FnyFt, as follows:
(10)Fn=fn1⋮fn4∈ℝ4,Ft=ft1⋮ft4∈ℝ4(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:(11)Hstn=0,Jlq̇=0.VectorHstan∈ℝ2 contains the height of the two corners of the foot in contact with the ground, and Jlq̇ were the three-element vector with the linear and the angular velocity of the foot defined in Equation (5).
(12)Pswn=Pswnd.VectorPswn∈ℝ12 contains the Cartesian position of the four corners of the swing foot, and Pswingd 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 Figure1, this constraint implies(13)q1T=−q40,q2T=−q50,q3T=−q60,q4T=−q10,q5T=−q20,q6T=−q30(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 Fn and Φn must be zero. This condition and nonnegativity conditions on Fn and Φn can be written compactly using the ⊥ symbol that indicates perpendicularity between vectors:(14)Φn≥0⊥Fn≥0Φn is derived twice with respect to time to make its dependency on Fn explicit. However, the condition Φn≥0 is equivalent to Φ¨n≥0 only for the contact points that satisfies φni=0 and φ̇ni=0. To obtain the equivalent of Equation (14) but involving Φ¨n rather than Φn, we define Nc as the set of vertices of the foot satisfying the two aforementioned conditions, and Φ¨nc as the vector with the normal acceleration of the elements of Nc. Similarly, Fnc is the vector with the normal reaction forces of Nc. The previous considerations lead to Equation (15).
(15)Φ¨nc≥0⊥Fnc≥0.Next, the complementarity equations that arise fromCoulomb’s friction law are deduced.ConsiderNs⊆Nc, 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 i, with i∈Ns, the complementarity condition between the tangential velocity and reaction force is [21]
(16)μfni−fti≥0⊥φ̇ti≥0,φ̇tifti–μfniφ̇ti=0,i∈Ns.Equation (16) ensures that (i) the tangential reaction force will be inside the friction cone fti≤μfni; (ii) if φ̇ti>0, then fti=μfni; and (iii) scalars fti and φ̇ti will have opposite signs.The solution of the complementarity problemgx≥0,x≥0, and gxTx=0 is equivalent to solve Equation (17) [22, 23]:
(17)gix2+xi2−gix+xi=0,wheregix and xi 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:
(18)Gq,q̇,Fl,Fr=0.The unknowns in Equation (18) are Fr and Fl. In summary, the robot’s dynamic behavior is described by the differential equation (1) subject to the algebraic constraint Equation (18) like
(19)Aqq¨+Hq,q̇=BΓ+JrTqFr+JlTqFl,Gq,q̇,Fl,Fr=0.


## 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 vectorq̇. 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 t=tI, the reaction force on the contact point i is assumed to equal to pniδt−tI, being pni the magnitude of the impulse and δt−tI the Dirac distribution. The normal speed of the contact i just after impact denoted φ̇ni+ and the magnitude of the impulse pni are subject to the following complementarity constraint (Equation (20)) [24]:
(20)φ̇ni++ϵφ̇ni−≥0⊥pni≥0,i∈Nc∪NI,whereNI is the set of the vertices of the foot that impact the ground in the instant tI. 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, denotedpti, is related to tangential velocities φ̇ti+, through
(21)μpni−pti≥0⊥φ̇ti+≥0,φ̇ti+pti−μpniφ̇ti+=0,i∈Nc∪NI.qv+ depends in turn on pnc y de ptc. 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
(22)qv+=qv−+A−1JrTFr+JlTFl.Being in this case,q̇v+ is the unknown of the problem. If in the instant t=tI at least one of the vertices of the foot satisfies the conditions φi=0 and φ̇i<0, 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 vectorsq and q̇ using a first-order difference:
(23)qk+1=qk+hq̇k,q̇k+1=q̇k+hA−1BΓ−H+JrTFrk+JlTFlk.qk y q̇k+1 represents the values of vectors q and q̇ in the instant tk. The constant h is defined as h≜tk+1−tk. In Euler’s symplectic method, q̇k+1 is calculated using qk+1 instead of qk. In the time-stepping approach, the Signorini condition is replaced by the following complementary condition:
(24)φ̇nik+1≥0⊥fnik,i∈Nυ.Nv denotes the set of vertices of the foot that satisfy the condition φni=0. The condition Equation (24) is interpreted as follows: if φni is equal to zero, the normal reaction force calculated at the instant tk must produce a positive or zero speed at the instant tk+1 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 tk and the tangential velocity in the instant tk+1 is
(25)μfnitk−ftitk≥0⊥φ̇titk+1≥0,φ̇titk+1ftitk–μfnitkφ̇titk+1=0,i∈Nν.The vectorsFnkyFtk comprises the reaction forces of the points contacting with the ground. The contact forces of the vertices that do not belong to the set Nv 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
(26)qk+1=qk+hq̇k,q̇k+1=q̇k+hA−1BΓ−H+JrTFrk+JlTFlk,0=Gqk,q̇k,Fnk,Ftk.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
(27)Aqq¨+Hq,qv=BΓ+JlTqFl.The model Equation (27) has unknowns that are q¨∈ℝ9andFl∈ℝ3. 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
(28)gplxq≡0,gplyq≡0,θlq≡0.By deriving Equation (28) once:
(29)gṗlxgṗlyθ̇l=Jlqqv=0.By deriving Equation (29):
(30)gp¨lxgp¨lyθ¨l=Jlqq¨v+Jlqqv=0.Combining Equations (27) and (30), a new model relating q¨ and Fl as a function of q, q̇, and Γ is obtained in
(31)Aq−JlTqJlq0q¨Fl=BΓ−Hq,qvJ̇lqqv.Rewriting Equation (30), we obtained
(32)Jluqq¨u+Jlaqq¨a+J̇lqqv=0.In Equation (33), qu is the vector containing the unactuated coordinates and qa is the actuated coordinate vector.
(33)q¨u≜gx¨0gy¨0q¨0∈ℝ3,q¨a≜gq¨1⋮q¨6∈ℝ6.MatrixJluq∈ℝ3×3 is composed by the first three columns of Jlq and Jlaq∈ℝ3×6 by the last six columns. Solving q¨u of Equation (32):
(34)q¨u=q¨uq¨a=−Jlu−1qJlaqIq¨a+−Jlu−1qJ̇lqqv0.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:
(35)Mq≜−Jlu−1qJlaqI,Nq,qv≜−Jlu−1qJ̇lqqv0.Substitutingq¨=Mq¨a+N in the model Equation (27)
(36)AMq¨a+AN+H=BΓ+JlTFl.Equation (36) is multiplied by MT; then
(37)MTAMq¨a+MTAN+MTH=Γ.Ifq¨a is replaced by the desired acceleration q¨ad, and q¨ad is calculated by Equation (38), then Equation (37) becomes a classical CTC [27]
(38)q¨ad=q¨r+kνqvr−q̇a+kpqr−qa.qar∈ℝ6,q̇ar∈ℝ6,andq¨ar∈ℝ6 being the joint reference motion. By replacing q¨ad 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
(39)Γ=MTAMq¨ad+MTAN+MTH,q¨ad=q¨r+kνqvr−q̇a+kpqr−qa.


### 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
(40)J=∫0TΓTq,qν,q¨Γq,qν,q¨dt,whereΓq,qv,q¨ represents the torque vector exerted by the actuators as a function of generalized positions q, velocities qv, and accelerations q¨. The numerical value of Γq,qv,q¨ 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 J no longer requires finding a function of time qt, 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 qat like
(41)qat≜a0+a1t+a2t2+a3t3,0≤t≤,0,t>T.Beingai∈ℝ6i=0,⋯,3, the vectors that contain the coefficients of the 6 polynomials. The coefficient vectors ai can be defined in terms of the initial and final positions and angular velocities
(42)qa0=a0,qaT=a0+a1T+a2T2+a3T3,q̇a0=a1,q̇aT=a1+2a2T+3a3T2.By solving the previous system, we obtained(43)a0a1a2a3=I,0,0,00,0,I,0−3T−2I,3T−2I,−2T−1I,−T−1I2T−3I,−2T−3I,T−2I,T−2Iqa0qaTq̇a0q̇aT,whereI is an identity matrix of 6×6 and 0 a zero matrix of 6×6. Of this way, instead of directly searching for 24 coefficients, the vectors qa0∈ℝ6,qaT∈ℝ6,q̇b0∈ℝ6, and q̇aT∈ℝ6 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
(44)px=mzfy,wherepx is the zero moment point, mz the angular momentum on axis z, and fy the reaction force of the ground exerted on foot. To maintain the balance of robot, the condition described by Equation (45) must be satisfied.
(45)mz−fypx=0,0−px1⏟Cpxfxfymz⏟Fl=0.Ifpx 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 px<0.5lf, where lf=0.2m is the foot length.


## 4.1. Nonlinear Control for the Single Support Phase

When the robot is resting on the left foot, the model Equation (1) is
(27)Aqq¨+Hq,qv=BΓ+JlTqFl.The model Equation (27) has unknowns that are q¨∈ℝ9andFl∈ℝ3. 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
(28)gplxq≡0,gplyq≡0,θlq≡0.By deriving Equation (28) once:
(29)gṗlxgṗlyθ̇l=Jlqqv=0.By deriving Equation (29):
(30)gp¨lxgp¨lyθ¨l=Jlqq¨v+Jlqqv=0.Combining Equations (27) and (30), a new model relating q¨ and Fl as a function of q, q̇, and Γ is obtained in
(31)Aq−JlTqJlq0q¨Fl=BΓ−Hq,qvJ̇lqqv.Rewriting Equation (30), we obtained
(32)Jluqq¨u+Jlaqq¨a+J̇lqqv=0.In Equation (33), qu is the vector containing the unactuated coordinates and qa is the actuated coordinate vector.
(33)q¨u≜gx¨0gy¨0q¨0∈ℝ3,q¨a≜gq¨1⋮q¨6∈ℝ6.MatrixJluq∈ℝ3×3 is composed by the first three columns of Jlq and Jlaq∈ℝ3×6 by the last six columns. Solving q¨u of Equation (32):
(34)q¨u=q¨uq¨a=−Jlu−1qJlaqIq¨a+−Jlu−1qJ̇lqqv0.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:
(35)Mq≜−Jlu−1qJlaqI,Nq,qv≜−Jlu−1qJ̇lqqv0.Substitutingq¨=Mq¨a+N in the model Equation (27)
(36)AMq¨a+AN+H=BΓ+JlTFl.Equation (36) is multiplied by MT; then
(37)MTAMq¨a+MTAN+MTH=Γ.Ifq¨a is replaced by the desired acceleration q¨ad, and q¨ad is calculated by Equation (38), then Equation (37) becomes a classical CTC [27]
(38)q¨ad=q¨r+kνqvr−q̇a+kpqr−qa.qar∈ℝ6,q̇ar∈ℝ6,andq¨ar∈ℝ6 being the joint reference motion. By replacing q¨ad 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
(39)Γ=MTAMq¨ad+MTAN+MTH,q¨ad=q¨r+kνqvr−q̇a+kpqr−qa.


## 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
(40)J=∫0TΓTq,qν,q¨Γq,qν,q¨dt,whereΓq,qv,q¨ represents the torque vector exerted by the actuators as a function of generalized positions q, velocities qv, and accelerations q¨. The numerical value of Γq,qv,q¨ 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 J no longer requires finding a function of time qt, 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 qat like
(41)qat≜a0+a1t+a2t2+a3t3,0≤t≤,0,t>T.Beingai∈ℝ6i=0,⋯,3, the vectors that contain the coefficients of the 6 polynomials. The coefficient vectors ai can be defined in terms of the initial and final positions and angular velocities
(42)qa0=a0,qaT=a0+a1T+a2T2+a3T3,q̇a0=a1,q̇aT=a1+2a2T+3a3T2.By solving the previous system, we obtained(43)a0a1a2a3=I,0,0,00,0,I,0−3T−2I,3T−2I,−2T−1I,−T−1I2T−3I,−2T−3I,T−2I,T−2Iqa0qaTq̇a0q̇aT,whereI is an identity matrix of 6×6 and 0 a zero matrix of 6×6. Of this way, instead of directly searching for 24 coefficients, the vectors qa0∈ℝ6,qaT∈ℝ6,q̇b0∈ℝ6, and q̇aT∈ℝ6 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
(44)px=mzfy,wherepx is the zero moment point, mz the angular momentum on axis z, and fy the reaction force of the ground exerted on foot. To maintain the balance of robot, the condition described by Equation (45) must be satisfied.
(45)mz−fypx=0,0−px1⏟Cpxfxfymz⏟Fl=0.Ifpx 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 px<0.5lf, where lf=0.2m is the foot length.


## 5. Numeric Results

Figure3 shows the tracking of trajectories obtained by Equation (40) and the control law Equation (39).Figure 3
Tracking of the reference motion.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 qrt,q̇rt and q¨rt are calculated off-line.Figure 4
Tracking of the velocity reference motion.The temporal evolution of the ZMP is shown in Figure5. The ZMP fulfill the condition px<0.5lf, and the control law ensures robot balance.Figure 5
Position of the ZMP for the CTC control. The boundaries of the robot’s support polygon are−0.1,0.1 represented by the dashed horizontal lines.The tangential force is within the limits of the cone of friction for friction coefficientμ=0.75, like Figure 6. Then, the robot does not slip.Figure 6
Tangential force.Ftt≤μFnt.


## 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 instanttk+1 depending on the product between reaction forces at time tk and period discretization h.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.



---

*Source: 5539123-2021-05-04.xml*