Calculation of the Center of Mass Position of Each Link of Multibody Biped Robots

In this paper, a novel method to determine the center of mass position of each link of human-like multibody biped robots is proposed. A first formulation to determine the total center of mass position has been tested in other works on a biped platform with human-like dimensions. In this paper, the formulation is optimized and extended, and it is able to give as output the center of mass positions of each link of the platform. The calculation can be applied to different types of robots. The optimized formulation is validated using a simulated biped robot in MATLAB.


Introduction
One of the pioneers in the field of biped robots was the Waseda University of Tokyo.In 1973, research groups from Waseda University developed WABOT-I, and in 1984 WABOT-2 as to become a professional musician [1].In 1999, they developed a humanoid with a complete human configuration capable of biped locomotion, WABIAN (Waseda Bipedal humANoid), and in 2011 its Italian version SABIAN [2][3][4].In 2017, a more complex version of the WABIAN robot is presented in [5].In 2013, Google acquired eight advanced robot companies.Boston Dynamics, one of these, is known for its advanced robots including the world's fastest robot, Cheetah [6] (which can travel at 29 mph), Big Dog [7] (the all-rough and tough robot that walks, runs, climbs, and carries heavy loads), and the latest Atlas (a biped humanoid capable of walking in outdoor rough terrain with the upper limbs capable of performing other tasks while walking).Several versions of Atlas [8,9] have been prepared for the DARPA Robotics Challenge program.In 2017, Boston Dynamics presented a very innovative robot with higher locomotion capabilities including wheels in the feet [10].
The general robot design process includes many phases, like the design of every complex machine.A tentative method to define these phases is shown as follows: • PHASE 1-Determination of the technical specifications which define limits and characteristics that the robot should have.• PHASE 2-Conceptual design of the robot including analysis of the developed robots in the world; design of novel systems; definition of the whole system including mechanics, electronics, low and high level control.• PHASE 3-Functional design of the robot including interaction of the robot with the environment; theoretical formulation and optimization; software and hardware design of virtual models (virtual model prototyping using CAD tools, analytical simulations, finite elements analysis, multibody analysis).• PHASE 4-Development of the robot including rapid prototyping modeling and tests.
• PHASE 5-Realization of the final robot prototype and final tests.
All robots realized in Phase 5 should have, theoretically, the same architecture of the virtual models designed in Phase 3. In particular, the positions of the CoM (Center of Mass) of each link of the robot should have the same position defined in the virtual model.However, the total CoM of the real platform is not in the same calculated position of the total CoM of its respective virtual model.Why are these discrepancies are created?Are these errors (between virtual models and real robots) influences on the robot functionality?
In order to answer to these two questions we could show how the problem was evident in the SABIAN robot (height: 1500 mm; 64 kg).During our tests on the platform, we noted a weight difference of about 5 kg including differences of CoM position between real and virtual SABIAN.These discrepancies are created because during manufacturing, construction, and maintenance of the robot, the tolerances of the joints, cables, drivers, batteries, links, etc. could not be completely respected and errors are created.These errors influence the robot functionality and cannot be eliminated [11].Hence, the real center of mass (CoM) of the robot is not coincident with the CoM of its virtual model [12].In our experiments, the robot SABIAN [2,3] had about 5 kg of errors in an unknown position and during locomotion, the controller implemented on the virtual model architecture was not able to control the real platform with this unknown error position.Dynamic balance in locomotion is not simple to control in a biped platform and these errors may disturb biped stability.
Some researchers and specialists in the humanoid robotics field use a posture controller [11,12], in order to reduce the error between the robotic platform and its virtual model.Kwon et al. [13] (2007) proposed a method that uses a closed-loop observer based on a Kalman filter, adopted as estimation framework.Ayusawa et al. (2008) [14] proposed a method based on regression analysis models in order to estimate inertial parameters using a minimal set of sensors.In the work of Sujan and Dubowsky (2003) [15] the dynamic parameters of a mobile robot are calculated using an algorithm based on a mutual-information-based theoretic metric for the excitation of vehicle dynamics.Liu et al. (1998) [16], Khalil et al. (2007) [17], and Swevers et al. (1997) [18] show other methods oriented to improve the balancing performances of mobile biped robots when the center of mass is not precisely known.
In [2,3], the authors proposed a novel approach to determine the correct position of the center of mass in humanoid robots.In order to compensate the errors between the biped platform and its virtual model, an additional mass has been implemented in the virtual model of the humanoid robot.The value of this mass error is the analytical difference between the weight of the robot and the weight of its virtual model.Its position in the space is not known a priori, but it will be approximately calculated with the procedure described in [3].In order to define its position, the authors of the paper proposed an analytic formula that gives the real position of the CoM of the platform and is based on the application of a procedure that requires only the values of the force-torque sensors, applied on the feet of the humanoid robot, and the values of the motors torque.This procedure standardizes the calibration procedure in order to minimize the errors and it can be applied to every biped platform.The formulation approach has been implemented on the SABIAN robot with dimensions comparable to humans (height: 1500 mm; weight: 64 kg) giving very good results.
The limits underlined in the papers [2,3] were based on the approximation used to put the error mass in the determined CoM with the proposed formula.The approximation has been justified because the real position of the error mass is not known and if an external mass is positioned in the real CoM, its negative influence is reduced.However, the problem remains because an approximation has been used on behalf of exact calculation.
In this paper, the limits underlined in [2,3] are bypassed with the optimization of the formula based on the determination of the CoM position of each link of the robot.With the proposed solution, the error mass is distributed on each link of the biped robot.
Another advantage of the formulation presented in this paper is that if the total CoM position of the platform is known a priori, the first formulation proposed in [2,3] can be bypassed and the CoM positions of each link of the platform can be calculated analytically without using force-torque sensors and the motor's torque.
The paper is structured as follows: Section 2 presents, in synthesis, the first validated theoretical formulation proposed in [2,3]; Section 3 shows results and discussion on the second theoretical formulation to determine the center of mass position of each link of the platform.Section 4 presents validation of the second theoretical formulation.The paper ends with a conclusion and future works.

Dynamics of Multibody Biped Robots
Figure 1 shows global and local reference Cartesian system (respectively G-XYZ and P-XpYpZp and three points in the space (0, 1 and 2).The three points can be considered as belonging to a rigid body in the space; furthermore, the rigid body can be compared to a humanoid platform, or multibody robot, with its center of gravity in the Point 2 and its feet in the Points 0 and 1.A humanoid robot is indeed composed of a trunk and articulated kinematic chains such as legs and arms, connected to the trunk with joints and motors, and with force-torque sensors positioned on the feet and on the hands.The Points 0 and 1 represent the feet, where the force-torque sensors are positioned, and are shown with a light blue colour; the center of mass represented with the Point 2 is shown in red; the other black points indicate the center of mass of the links of the robot.

Dynamics of Multibody Biped Robots
Figure 1 shows global and local reference Cartesian system (respectively G-XYZ and P-XpYpZp and three points in the space (0, 1 and 2).The three points can be considered as belonging to a rigid body in the space; furthermore, the rigid body can be compared to a humanoid platform, or multibody robot, with its center of gravity in the Point 2 and its feet in the Points 0 and 1.A humanoid robot is indeed composed of a trunk and articulated kinematic chains such as legs and arms, connected to the trunk with joints and motors, and with force-torque sensors positioned on the feet and on the hands.The Points 0 and 1 represent the feet, where the force-torque sensors are positioned, and are shown with a light blue colour; the center of mass represented with the Point 2 is shown in red; the other black points indicate the center of mass of the links of the robot.
The dynamics of the system is described by the two equations (see Figure 1): where: ] respectively represent the torques ( , , ) and the forces ( , , ) acting at the Points 0, 1 and 2.
= [ , , ] is the resultant moment calculated with respect to the local Cartesian system (see Figure 1).The direction of the force and torque vectors, shown in Figure 1, is only indicative; the positive direction of the force and torque vectors has been considered with the same positive direction of the global reference Cartesian system (G-XYZ).The dynamics of the system is described by the two equations (see Figure 1): M 2 includes only the internal torques; M 0 and M 1 include the ground reaction torques.The Equation ( 3) is based on the assumption that the Jacobian Matrix is equal to the identity Matrix.This assumption is correct if the joints axis of the robot remain parallel, during motion, to the Y axis of the global reference Cartesian system (G-XYZ).It means that the motion of the robot for the determination of the formula is performed in a 2D plane.In this paper, the XZ plane is considered.∑ M j_roll , ∑ M j_pitch , ∑ M j_yaw are the torques of all roll, pitch, and yaw motors of the robot [11] and M j is obtained from the equation: The accuracy in the estimation of M j , calculated in (4) depends on the accuracy of the K value that is a constant parameter set for each motor and on the accuracy of the current I necessary for the motor function.In particular, the resolution of the used A/D converters is a fundamental parameter to define the accuracy of the current I.

Equilibrium
Considering the robot equilibrium (m 2 → a 2 = 0; → M P = 0) with respect to the Point P as shown in the Figure 1, Equations (1) and ( 2) can be modified.Moving the point of view from the vector shape to the scalar one, the values of the three components x, y, and z of the force and the torque can be obtained.The new system consists of six equations (five linearly independent) in six unknown values The forces and torques in Points 0 and 1 can be calculated by means of the load cells.The torques M X2 , M Y2 , M Z2 are determined using (3) and (4).
In order to simplify the system, the robot is positioned in two different configurations.The two configurations are chosen in order to have a simplified geometry using l = q, n = o = e = f = 0 (see Figure 1) obtaining In the first step the robot is placed on a walking surface and the platform should be kept in a first balance configuration (Scheme A, see Figure 2), allowing a measurement of the forces ( ) by means of the force-torque sensors on the feet, and of the armature currents (∑ M j_roll , ∑ M j_pitch , ∑ M j_yaw ).In a second step, the robot is placed in a second balance configuration (Scheme B, see Figure 2), and in the same way, forces, torques, and motor currents associated with this new balance configuration are measured.The two balance configurations can be performed as the reader prefers underlining that the robot should be in a balance position.In particular, the coordinates X 2A and Z 2A are relative to the position of the center of mass of the platform in the configuration A along the first straight common line, which is chosen (in this paper) orthogonal to the plane of standing, and then parallel to the Z axis.m u , r u , l u , are the mass and the position of the center of mass of the robot ankle link from the floor to the ankle joint; U is the length of the ankle.m w , r w , l w , are the mass and the position of the center of mass of the remaining links of the platform.In the second balance configuration B, the components of the body are aligned according to a second straight common line, inclined to the vertical line with an angle θ t .While m u , m w , r u , l u remain constant, r w and l w change their values.In this case, X 2B and Z 2B identify the coordinates of the center of mass of the body in the second balance configuration B. The coordinates of the two feet are the same because we chose this configuration as input.The implementations have been done positioning the robot in this initial position using a leveller and the encoders of the motors.

Proposed Coefficients
In following, four novel coefficients ( , , , ) are introduced Rewriting Equations ( 8) and ( 9) (for Using the parameters mu, mw, ru, lu, rw and lw, the Equations ( 16) and ( 17) can be obtained.Thus, the Equations ( 14), ( 16) and ( 17) can be seen as a system composed of six equations in six unknown variables (for i = A and i = B) X2A, Z2A, X2B, Z2B, rw, lw; the relation between mw and mu is given by the Equation (18).θt is fixed by the user (θt = 0 in i = A), in a way that does not allow to tilt the platform.

Proposed Coefficients
In following, four novel coefficients (α i , β i , γ i , δ i ) are introduced Rewriting Equations ( 8) and ( 9) (for I = A or I = B) Using the parameters m u , m w , r u , l u , r w and l w , the Equations ( 16) and ( 17) can be obtained.Thus, the Equations ( 14), ( 16) and ( 17) can be seen as a system composed of six equations in six unknown variables (for i = A and i = B) X 2A , Z 2A , X 2B , Z 2B , r w , l w ; the relation between m w and m u is given by the Equation ( 18).θ t is fixed by the user (θ t = 0 in i = A), in a way that does not allow to tilt the platform.16) Solving the equations system constituted by ( 14), ( 16) and ( 17), the positions of the center of mass are calculated in both the configurations A and B (for i = A and i = B).It must be underlined that only X as a function of Z has been considered, but the same result can be obtained considering Y as a function of Z.

Determination of the Partial Center of Mass Position
Placing θ t = 0 (then i = A) and substituting ( 16) and ( 17) into ( 14) and placing θ t = 0 (then i = B) and substituting ( 16) and ( 17) into (14), two different equations will be obtained.Finally, combining these two equations l w and r w are obtained.
Placing θ t = 0 (then i = B) and rewriting (17) with the latter values given by l w and r w , Z 2B is obtained.
Z 2B represents the general position of the height Z of the center of mass for any value of θ t .Placing i = B in ( 14) and ( 15) and substituting the found value of Z 2B , the general formula of the position of the center of mass in (19) is given.Considering the equilibrium configuration, A (i = A) and then θ t equal to zero, the system (20) is obtained.
In particular, the calculation of the center of mass position is strictly related to the proposed coefficients α A , β A , γ A , δ A , α B , β B , γ B , δ B , that are numerical values associated with the first and second measurements on the robot.In order to calculate these coefficients, it is necessary to consider the arbitrary θ t associated with the second balance configuration, in addition to other parameters such as the above mentioned position of the center of mass of the feet (r u and l u ) and its mass (m u ) calculated using the CAD model.These parameters have a lower weight with respect to other links of the platform and then a lower inertial influence [3].The two formulations (19) and (20) were tested and validated on the SABIAN platform giving very good results and presented in [2,3].

Procedure for n Degrees of Freedom
In order to find the CoM positions of each link of the platform, the following procedure and formulation must be used.In particular, if n are the degrees of freedom of the platform; i represents the used configurations to calculate the CoM of each link.j represents the initial configuration.The number (k) of the configurations necessary to calculate the CoM position for each link of the system is calculated in the following The total CoM position (X 2i , Y 2i , Z 2i ) of the complete system for each configuration i can be obtained with the following formulas where the coefficients α i , β i , γ i , δ i , α j , β j , γ j , δ j , are calculated using respectively ( 10)- (13).
The following formulas from ( 23) to (26) allow to determine the CoM positions of each link of the robot.m w , m 2 , m u , r u , q u , l u , U, are the input of the system as shown in the Section 2. r wi , q wi , l wi , are the components of the vector position → t i = r wi q wi l wi + U T respectively in X P , Y P , and Z P directions (see Figure 1).χ i and ε i are the angles of the vector position → t i = r wi q wi l wi + U T as shown in Figure 3.
Appl.Sci.2017, 7, 724 7 of 14 The following formulas from ( 23) to (26) allow to determine the CoM positions of each link of the robot.mw, m2, mu, ru, qu, lu, U, are the input of the system as shown in the Section 2. rwi, qwi, lwi, are the components of the vector position = [ + ] respectively in XP, YP, and ZP directions (see Figure 1).and are the angles of the vector position = [ + ] as shown in Figure 3.

Procedure for n = 2 Degrees of Freedom
Figure 4 shows a sketch of a robot in a plane XP, ZP with two degrees of freedom (n = 2) and three links (La, Lb, Lc).From (21) k = 3 is obtained.This result means that two configurations must be used to calculate the coefficients , , , , and one configuration must be used to calculate the coefficients , , , using respectively ( 10)-( 13) in order to find the CoM position of the links of the robot.In particular, the following iterative procedure should be used: The robot is placed on a walking surface and the platform should be kept in a first balance configuration j allowing a measurement of the forces ( , ) and the torques ( , ) by means Appl.Sci.2017, 7, 724 8 of 15 →

Procedure for n = 2 Degrees of Freedom
Figure 4 shows a sketch of a robot in a plane X P , Z P with two degrees of freedom (n = 2) and three links (L a , L b , L c ). From (21) k = 3 is obtained.This result means that two configurations must be used to calculate the coefficients α i , β i , γ i , δ i , and one configuration must be used to calculate the coefficients α j , β j , γ j , δ j using respectively ( 10)-( 13) in order to find the CoM position of the links of the robot.In particular, the following iterative procedure should be used:

•
The robot is placed on a walking surface and the platform should be kept in a first balance configuration j allowing a measurement of the forces ( → F 0 , → F 1 ) and the torques ( ) by means of the force-torque sensors on the feet, and of the armature currents (∑ M j_roll , ∑ M j_pitch , ∑ M j_yaw ).These values are used to calculate the coefficients α j , β j , γ j , δ j , of the (22) using respectively ( 10)-( 13); The robot is placed in a second and third balance configuration i, and in the same way, forces, torques, and motor currents associated with each balance configuration (second and third) are measured.These values are used to calculate the coefficients α i , β i , γ i , δ i using respectively (10)-( 13).

•
For each balance configuration i, the total CoM position is calculated using (22); • Each total CoM position allows to determine r wi , q wi , l wi , using (23) and the CoM position of each link of the robot using (27).For each configuration i, an equation using (26) is created.n linearly independent equations are used to find n vector positions.
Appl.Sci.2017, 7, 724 8 of 14 of the force-torque sensors on the feet, and of the armature currents ( ∑ _ , ∑ _ , ∑ _ ).These values are used to calculate the coefficients , , , , of the (22) using respectively (10)-( 13); The robot is placed in a second and third balance configuration i, and in the same way, forces, torques, and motor currents associated with each balance configuration (second and third) are measured.These values are used to calculate the coefficients , , , using respectively (10)-( 13).

•
For each balance configuration i, the total CoM position is calculated using (22); • Each total CoM position allows to determine rwi, qwi, lwi, using (23) and the CoM position of each link of the robot using (27).For each configuration i, an equation using (26) is created.n linearly independent equations are used to find n vector positions.

Implementation of the Analytical Formulation
The example shown in Figure 4 has four unknown parameters (rb, lb, rc, lc), but four linearly independent equations can be obtained by the formulas shown in Section 3.1.In this case, we suppose to have ma, mb, mc, ra, la, as input.Figure 5 shows in details the configurations b and c shown in Figure 4.In particular, the four unknown parameters (rb, lb, rc, lc) which should be determined using analytical formulation are found using polar coordinates.and are the position vectors of the mass mw

Implementation of the Analytical Formulation
The example shown in Figure 4 has four unknown parameters (r b , l b , r c , l c ), but four linearly independent equations can be obtained by the formulas shown in Section 3.1.In this case, we suppose to have m a , m b , m c , r a , l a , as input.Figure 5 shows in details the configurations b and c shown in Figure 4.In particular, the four unknown parameters (r b , l b , r c , l c ) which should be determined using analytical formulation are found using polar coordinates.
Using ( 23)-(25), rwb, lwb, and rwc, lwc, can be determined (see Figure 5) and the modules and angles of the two position vectors ( and ) are found Using ( 26) and ( 27), the modules of the position vectors and (respectively the position vectors of the masses mb and mc respect to the local reference system of each link) are determined.In the following, detailed equations are shown which can determine the four unknown parameters (rb, lb, rc, lc).Using (22), Z 2b , X 2b , and Z 2c , X 2c , can be calculated as shown in following Using ( 23)-(25), r wb , l wb , and r wc , l wc , can be determined (see Figure 5) and the modules and angles of the two position vectors ( Using ( 26) and ( 27), the modules of the position vectors    The proposed representation is general and can be implemented on robots with different types of joints (prismatic, revolute, spherical, helical, etc.).In case of robot conceived in an unconventional way, such as passive or flexible robots or robot with wheels [19][20][21], the procedure (21) and the formulas from ( 22) to (27) can be implemented if the two following points are satisfied: 1.
Only two force-torque sensors must be the contact elements between the robot and the ground; 2.
Joint sensors must give relative position of motion and current values to produce join motion.
M P = [M XP , M YP , M ZP ] T torque vector of the resultant moment calculated with respect to the local Cartesian system P-XpYpZp ∑ M j_roll , ∑ M j_pitch , ∑ M j_yaw torques of all roll, pitch, and yaw motors of the robot X 2A , Y 2A , Z 2A position of the centre of mass of the platform in the configuration A X 2B , Y 2B , Z 2B position of the centre of mass of the platform in the configuration B m u , r u , q u , l u mass and position of the centre of mass of the robot ankle link from the floor to the ankle joint U length of the ankle m w , r w , q w , l w mass and the position of the centre of mass of the remaining links of the platform θ t angle used in the configuration B K constant parameter set for each motor I current necessary for the motor function α i , β i , γ i , δ i four novel coefficients for the configuration i α j , β j , γ j , δ j four novel coefficients for the configuration j n degrees of freedom of the platform k number of the configurations to calculate positions of the CoM for each link of the system i, j used configurations to calculate the CoM of each link → t i = r wi q wi l wi + U T vector position of m w r wi , q wi , l wi components of the vector position → t i = r wi q wi l wi + U T respectively in X P , Y P and Z P directions χ i and ε i angles of the vector position → t i = r wi q wi l wi + U T with the the local reference system P-X P Y P Z P the acceleration of the CoM; m2 is the total mass of the robot without feet; = [− , − , ] and = [ , , ] are the position vectors shown in Figure 1.= [ , , ] is the CoM position that will be determined with the proposed formula in Section 3

Figure 1 .
Figure 1.Scheme for determining the total CoM (Center of Mass) position.The direction of the force and torque vectors is only indicative.Points 0 and 1 represent the feet, where the force/torque sensors are positioned; the center of mass is represented with the Point 2; the other black points indicate the mass of the links of the robot.

Figure 1 .
Figure 1.Scheme for determining the total CoM (Center of Mass) position.The direction of the force and torque vectors is only indicative.Points 0 and 1 represent the feet, where the force/torque sensors are positioned; the center of mass is represented with the Point 2; the other black points indicate the mass of the links of the robot.

Figure 2 .
Figure 2. Two examples of balance configurations; scheme A and scheme B. m2 is the total CoM.mu and mw are the CoMs of the two links.

Figure 2 .
Figure 2. Two examples of balance configurations; scheme A and scheme B. m 2 is the total CoM.m u and m w are the CoMs of the two links.

Figure 3 .
Figure 3. Scheme and vector position of the mass mw.

Figure 3 .
Figure 3. Scheme and vector position of the mass m w .

Figure 4 .
Figure 4. Sketch of a robot in a plane with 2 degrees of freedom.Three Configurations (a)-(c).m2 is the total mass of the system; ma, …, mc are the masses of the links and ra, …, rc and la, …, lc are respectively the x and z components of the center of mass of the links respect to the relative revolute joint; La, …, Lc are the lengths of the links.

Figure 4 .
Figure 4. Sketch of a robot in a plane with 2 degrees of freedom.Three Configurations (a)-(c).m 2 is the total mass of the system; m a , . . ., m c are the masses of the links and r a , . . ., r c and l a , . . ., l c are respectively the x and z components of the center of mass of the links respect to the relative revolute joint; L a , . . ., L c are the lengths of the links.

→
t b and → t c are the position vectors of the mass m w respect to the local reference system P-X P Y P Z P and respectively of the configuration b and c. → s b and → s c are respectively the position vectors of the masses, m b and m c , respect to the local reference system of each link.

Figure 5 .
Figure 5. Example for implementing the formula.

Figure 5 .
Figure 5. Example for implementing the formula.

→ s b and →
s c (respectively the position vectors of the masses m b and m c respect to the local reference system of each link) are determined.In the following, detailed equations are shown which can determine the four unknown parameters (r b , l b , r c , l c ). m

)
Solving the system using the Cramer's rule, we obtain the following equations where sin θ b = sθ b ; cos θ b = cθ b ; sin θ c = sθ c ; cos θ c = cθ c ; A = m b cθ b m b sθ b −m b sθ b m b cθ b m c cθ b m c sθ b −m c sθ b m c cθ b m b 0 0 m b −m c cθ c m c sθ c −m c sθ c m c sθ c (38)

a 1 m b sθ b a 2 m b cθ b a 1 −m b sθ b a 2 m c cθ b m c sθ b −m c sθ b m c cθ b m b a 3 0 a 4 −
m b cθ b m c cθ b m c sθ b −m c sθ b m c cθ b a 3 0 a 4 m b −m c cθ c m c sθ c −m c sθ c m c sθ c det(m c cθ c m c sθ c −m c sθ c m c sθ c det(A) (41) r c = m b cθ b m b sθ b −m b sθ b m b cθ b a 1 m c sθ b a 2 m c cθ b m b 0 0 m b a 3 m c sθ c a 4 m c sθ c det(A) (42) l c = m b cθ b m b sθ b −m b sθ b m b cθ b m c cθ b a 1 −m c sθ b a 2 m b 0 0 m b −m c cθ c a 3 −m c sθ c a 4 det(A) (43)

→
the mass m w respect to the local reference system P-X P Y P Z P and respectively of the configurations b and c the masses m b and m c respect to the local reference system of each link