An experimental energy consumption comparison between trajectories generated by using the cart-table model and an optimization approach for the Bioloid robot

The aim of this article is to present a statistical comparison of the electric energy expenditure between two techniques for the generation of gait patterns in biped robots. The first one is to minimize the sum of the squared torques applied to the joints of the robot, and the second one is based on the cart-table model. For the experiments, we measured the energy delivered by the battery of the robot to the servomotors. We applied the two aforementioned methods for three velocities (0.5, 1.0, and 1.3 m/min). Additionally, each combination of method and velocity was performed by the robot 10 times. The energy expenditure for each method was compared by applying the Wilcoxon test. In all comparisons, the value of p was lower than 0.004, indicating that the differences were statistically significant. The optimization approach leads to a reduction in energy expenditure that ranged from 9.16 % to 13.35 %. The conclusion is that all the effort required to implement an approach that requires a complete dynamic model of the robot allows a significant reduction in energy consumption.


Introduction
Robots and mechatronic systems are required to be energetically efficient, especially when they are conceived for high-speed continuous operations. 1 In the case of robotic manipulators, this topic has received considerable attention of many researchers, and as a consequence, different optimization solutions have been proposed. Woolfrey et al. 2 proposed a torque minimization approach to resolve redundancy in manipulators that handle large external forces. In the study of Trigatti et al., 3 a trajectory planning strategy for spray painting robots is presented.
The problem of efficiency has also been addressed in biped robots because their autonomous operation requires independent energy sources such as batteries. These solutions can be oriented to modify the mechanical structure of the robot, for example, introducing torsion springs between the joints and the actuators. 1,4,5 However, as stated by Scalera et al., 1 the main drawback of these solutions is that the stiffness of the springs cannot be adapted to different motions. Other approaches are oriented to modify the motion patterns of the robot. Currently, there are two main approaches to the generation of biped robot trajectories. One of them is known as the cart-table approach, which involves only the kinematic model of the biped and its total mass. [6][7][8][9] The mass of the car symbolizes the total mass of the robot, and its location corresponds to the center of mass (CoM) of the biped. 10,11 This method is used by bipeds such as ASIMO, 12,13 HRP, [14][15][16] Nao, 17,18 HUBO, 19 and ATLAS. 20,21 The other approach is based on numerical optimization. It considers all the components of a Lagrangian model, such as the inertia matrix, and the vectors of centrifugal, gravitational, Coriolis, and friction forces. An example of this approach is the technique proposed by Chevallereau and Aoustin 22 based on the minimization of the average energy expended by the robot. In this approach, two cost functions are considered. One of them is the timeintegral of the squared Euclidean norm of the torques applied to robot joints. This variable is a surrogate of Joule effect losses 23,24 when the relationship between the current delivered to electric motors and the joint torques is linear. The second cost function is the time-integral of the absolute value of the dot product between the torques and the angular velocities. 25,26 The main limitation of the aforementioned optimization schemes is that the convergence to a global minimum cannot be guaranteed. 27 Furthermore, even in the case that this convergence could be guaranteed, the differences between the dynamic of the real robot and that of the model used by the optimization algorithm convert the resulting trajectories into suboptimal. One of the sources for these differences are model actuators. 22,28,29 For example, when magnetic saturation exists in an electrical motor, the nonlinear equation that relates flux and armature current is difficult to obtain and, as a consequence, is not considered in the dynamic model. 30,31 Another factor that hinders the adoption of the numerical optimization approaches is their high computational cost. In the case of the current article, for example, the minimization space is composed of 60 dimensions.
Despite all the difficulties presented in the previous paragraph, optimization approaches based on Lagrangian models have been applied to biped robots. The largest part of these works present simulation results, 27,32-34 whereas the number of experimental studies is still scarce. 25,35 This imbalance is due to the fact that most simulation studies use simplified models, which are very valuable as proof of concept, but which cannot be considered valid mathematical representations of real robots. For example, in the study of Roussel et al., 32 the Coriolis matrix is not considered, and the studies by Huang and Hase 33 and Tlalolini et al. 34 consider bidimensional prototypes. Unlike the previous ones, the study of Tlalolini et al. 27 presents a threedimensional model with 14 joints. Another factor reducing the number of experimental evaluations of energy consumption is the need to modify the biped to add a current and voltage monitoring system.
In this article, two techniques for the generation of gait patterns in biped robots are compared. The first one is to minimize the sum of the squared torques applied to the joints of the robot, and the second one is based on the cart-table model. The aim of this work is to contribute with experimental evidence to answer the following question: Does the generation of trajectories based on numerical optimization lead to significant energy savings compared to the cart-table method? In addition, a penalty-based optimization approach that allows to generate joint trajectories of minimal energy expenditure in biped robots is proposed. Another contribution of the present work is a sound statistical framework to compare the energy expenditure of gait trajectories in biped robots.

Bioloid robot
The Bioloid humanoid (type A) from Robotis is an 18degree-of-freedom (DOF) robot, with 6 DOFs in each leg (see Figure 1). Each actuator is an AX-12A servo with position measurement of 0.29 resolution. These servos are controlled by the CM-530 board. When the stabilizing functionality is activated, a gyro sensor provides the angular rates around the anterior-posterior and medial-lateral axis to the CM-530, which modifies the trajectories of the servos. In Bioloid robots, the torques applied to the joints are not defined by the user but by proportional derivative (PD) controllers implemented as closed source firmware.

CAD model
The inertia moments and the CoMs of the robot were obtained from the computer-aided design (CAD) models provided by Robotis. Currently, it is unknown if these parameters include the contributions from electronics boards, wiring, and hollow spaces in plastic body. However, even if the CAD model is perfect, there are real-world phenomena that our model do not take into account or cannot be identified, as mechanical backlash in the gears, as well as the dead band of the PD controllers. The dynamical model of the robot was found using the Newton-Euler algorithm. 36,37 Dynamic model of the robot The Lagrangian representation of the robot is as follows where q 2 IR 18 is the vector of generalized coordinates. This vector contains (i) ð 0 x g ; 0 y g ; 0 z g Þ that are the Cartesian position of the trunk frame < x 0 ; y 0 ; z 0 > with respect to the ground frame < x g ; y g ; z g > (see Figure 2), (ii) s ¼ f; q; g are the Euler angles describing the orientation of the trunk frame < x 0 ; y 0 ; z 0 > with respect to the ground frame < x g ; y g ; z g >, and (iii) ðq 7 ; . . . ; q 18 Þ are the 12 joint positions of the legs q v 2 IR 18 is the vector of generalized velocities and q a 2 IR 18 is the vector of generalized accelerations. A 2 IR 18Â18 is the inertia matrix. H 2 IR 18 is the vector of centrifugal, gravitational, and Coriolis forces. B 2 IR 18Â12 is the actuation matrix that contains only ones and zeros. The first six rows of B are zero, indicating that G has no direct influence on components 1-6 of q a . The following 12 rows of B form an identity matrix, indicating that G i ði ¼ 7; Á Á Á ; 18Þ directly affects the component i of the vector q a . The vector G 2 IR 12 contains the torques applied to the joints of the robot. This leads B to satisfy the following property is the Jacobian matrix relating the linear and angular velocities of the robot's right foot with the vector q v . The super index g indicates that a variable is expressed in the ground frame In equation (4), g v 17 and g ! 17 are the linear and the angular velocities of the frame < x 17 ; y 17 ; z 17 > with respect to < x g ; y g ; z g >. J l 4 ¼ g J 18 2 IR 6Â18 is the Jacobian matrix relating the linear and angular velocities of the robot's left foot with the vector q v In equation (5), g v 18 and g ! 18 are the linear and angular velocities of the frame < x 18 ; y 18 ; z 18 > with respect to < x g ; y g ; z g >. The F r 2 IR 6 and F l 2 IR 6 contain the resulting of the contact forces between feet and ground. The six elements of the vector F r are the normal force f n r 2 IR, the tangential reaction forces f t r 2 IR 2 , and the ground reaction torques G c r 2 IR 3 . The components of F l are defined similarly to those of F r . When the left foot is in swing phase, then F l ¼ 0, and the unknowns of the model 1 are q a 2 IR 18 and F r 2 IR 6 . In such a case, there are 12 equations and 18 variables. The additional equations are obtained by taking into consideration that right foot does not move with respect to the ground reference frame. This constraint implies that the time-derivative of the right hand side of equation (4) is zero. Thus, when F l ¼ 0, the dynamical model is as follows The term _ J r ðqÞ q v is calculated by using the Newton-Euler algorithm. 36,37 The above described procedure can be applied to the other phases of gait. The kinematic and the dynamic parameters of the model equation (1) are presented in Tables 1 and 2, respectively.

Energy expenditure measurement system
The energy consumption measurement system developed by authors (see Figure 3) is composed of four electronic modules and a software component. The electronic subsystems are: (1) The DC power supply module: It utilizes a 12 V and 2200 mAH lithium-polymer (LIPO) battery to feed the motors and all electronic components, including those added to measure energy consumption. The power supply module includes external regulators to ensure that the reference for the analog to digital converters (ADC) of the Arduino Nano V3 board remains constant despite the voltage fluctuations across the terminals of the LIPO battery. This aspect is of capital importance to guarantee the reliability of the measurements delivered by the ADCs.
(2) The current measurement module: It measures the current delivered by the LIPO battery. This is done by using a shunt resistance of 0:28 O (1 % of tolerance, 4 W) and a differential amplifier, following the scheme proposed in the study of Regan et al. 38 The output of the differential amplifier is connected to one of the ADCs of the Arduino Nano V3 board. (3) The voltage measurement module: It measures the voltage between the terminals of the LIPO battery.
As the battery voltage is 12 V and the input to the ADCs of the Arduino Nano must range from 0 to 5 V, a voltage divider of gain 2=5 followed by a unit gain amplifier were used. (4) The wireless transmission module: It uses a Bluetooth HC-05 board to send every 50 ms the voltage and the current measured by the Arduino Nano board to a computer running Matlab R2018b (The MathWorks, Inc., Natick, Massachusetts, USA).  The weight of the board containing the energy measurement system is 100 g, and its dimensions are 7 Â 10 cm 2 . The software component stores the voltage and current samples sent by the wireless transmission module and calculates power and energy consumption. This last variable is calculated by using the trapezoidal integration.

Joint trajectory convention
For trajectory generation, we considered the 12 joints belonging to the legs (see Figure 2). The six joints for the right leg are q 7 ; q 9 ; q 11 ; q 13 ; q 15 , and q 17 , and the joints for the left leg are q 8 ; q 10 ; q 12 ; q 14 ; q 16 , and q 18 . Each hip has three DOFs. The joints q 7 and q 8 allow hip rotation in  the transverse plane, q 9 and q 10 in the frontal plane, and q 11 and q 12 in the sagittal plane. The knee has one DOF in the sagittal plane, q 13 for the right leg and q 14 for the left leg. The ankle has two DOFs. The joints q 15 and q 16 allow ankle rotation in the sagittal plane and q 17 and q 18 in the frontal plane.

Trajectory generation using the cart-table approach
The cart-table approach, described in the literature, 10,11 is based on the following assumptions: (i) Robot walks on flat surface perpendicular to the gravity axis. (ii) The spatial trajectories of both feet are the same. (iii) There is a double support phase. (iv) The height of the CoM with respect to the ground is assumed constant. (v) In the contact stage, the foot reaches parallel to the ground. Under the above assumptions, the relationship between the Cartesian position of the zero moment point 39 (ZMP) and the projection on the ground of the CoM is described by equation (7) P where g is the Earth's gravity. The constant z c is the desired height of the CoM. The pair ðP ZMP x ; P ZMP y Þ is the Cartesian position of the ZMP and ðx; yÞ is the Cartesian position of the CoM. The axis x is parallel to the gait direction and y is perpendicular to it. The eigenvalues of the linear differential equation (7) are þ ffiffiffiffiffiffiffiffiffiffi g=Z c p and À ffiffiffiffiffiffiffiffiffiffi g=Z c p , and consequently the dynamics is unstable. That means that given the input P ZMP x ðtÞ and P ZMP y ðtÞ, the resulting solutions xðtÞ and yðtÞ will tend to infinity. An approximate but boundedoutput solution of equation (7) can be obtained by applying the preview control method described by Kajita et al. 10 The following is the procedure to obtain the joint angles: (1) Select the trajectory of each foot. This requires defining both the parameters presented in Table 3 and the dimensions of the each foot (10 Â 6 cm 2 ). The trajectories of the feet are the inverted parabolas presented in Figure 4.
(2) Define the ZMP trajectory ðP ZMP x ; P ZMP y Þ. This is done by fixing a set of intermediate points and then interpolating. Most works use linear interpolation, and therefore the ZMP trajectory has a triangular shape. 10,41,42 In the current case, we use cubic splines to interpolate the midpoints presented in Table 4. (3) Calculate the trajectory of the COM by applying the preview control and using the ZMP obtained in the previous step as the reference signal. (4) Obtain the joint angles by using the trajectories of the COM and the feet as inputs to the inverse kinematic model.

Trajectory generation using the optimization approach
The purpose of the optimization approaches for the generation of gait patterns in bipedal robots is to find the articular trajectories that minimize a performance index, such as the torques exerted by the actuators or the mechanical energy expenditure. The criterion used in this work was proposed by Chevallereau and Aoustin 22 and Saidouni. 23 It consists Step length 4 Â 10 À2 m W p Step   of minimizing the functional described by equation (8), which is related to the energy losses in the actuators of the robot The vector G q; q v ; q a ð Þrepresents the torques exerted by the actuators as a function of generalized positions q, velocities q v , and accelerations q a . The numerical value of G q; q v ; q a ð Þ is calculated by using the inverse dynamic model of the robot. To transform the variational calculus problem defined by equation (8) 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 qðtÞ, but a vector of parameters, which is easier to find numerically. The only difference between our approach and the one proposed by Chevallereau and Aoustin 22 is the parametrization of qðtÞ. In that study, 22 the authors use fourth-order polynomials, and in our case a truncated Fourier expansion, which we consider more appropriate for periodic signals. In the present work, the angular trajectory of the joint i is parametrized in terms of the 2N þ 1 coefficients defined by equation (9) q i ðtÞ ¼ a i;0 2 þ X N n¼1 a i;n cosðn! 0 tÞ þ b i;n sinðn! 0 tÞ ð9Þ ! 0 ¼ 2p=T is the fundamental frequency of q i ðtÞ. The variable T represents the duration of the step (see Table 3). As a part of this study, the value of N was changed from 1 to 5 to see its effect in the performance index J (see equation (8)). A set of coefficients for equation (9) is considered feasible if the resulting trajectories satisfy the following equality and inequality conditions for 0 t T . 22 Inequality constraints: I 1 : When there is contact between one foot and the ground, the resulting normal reaction forces must be positive. That means that reaction forces can prevent interpenetration between the foot and the ground but cannot stick the foot to the ground f n r ! 0 I 2 : When there is contact between one foot and the ground, the resulting tangential reaction forces must be in the friction cone to avoid slipping ( is the Coulomb friction) jjf t r jj f n r I 3 : The absolute value of the joint torques must be less than 15kg Á cm, which is the maximum torque that the Dynamixel AX12 servos can exert jG i j G max ði ¼ 7; . . . ; 18Þ Equality constraints: (1) E 1 : The four 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 following equations H stance 2 IR 4 contains the height of the four corners of the foot in contact with ground, and J r q v is the six elements vector with the linear and the angular velocity of the foot defined by equation (4).
(2) E 2 : The four corners of the swing foot must follow the paths presented in Figure 4 P swing ¼ P d swing P swing 2 IR 12 contains the Cartesian position of the four corners of the swing foot and P d swing is the desired path. There are two main approaches to solving constrained optimization problems such as equation (8). One of them is using gradient-based methods, such as sequential quadratic programming, interior-point and trust-regionreflective. In these methods, when the gradient of the function to optimize is not available, it is approximated by finite differences. Nevertheless, in the case of the problems with a large number of variables, 60 in our case (12 trajectories with 5 parameters each), this approximation is computationally expensive and introduces numerical errors. The other approach is based on penalty-based methods 43 that turn constrained into unconstrained problems by introducing an artificial penalty for violating constraints. For example, the minimization of the function 100=x subject to the constraint x 5 can become an unconstrained optimization problem by adding a penalty function that returns a huge value when x > 5 min x 100 x þ a max 0; The term max x À 5; 0 ð Þis equal to zero for x 5 and equal to a ðx À 5Þ when x > 5. In this way, if a is large enough, then no value of x that violates the constraint could be a minimum of the function 10.
The minimization of equation (8), subject to the previously defined constraints I 1 , I 2 , I 3 , E 1 , E 2 , and E 3 , was transformed in the following unconstrained problem min q J q ð Þ þ P qðqÞ; q v ðqÞ; q a ðqÞ ð Þ ½ ð 11Þ The vector q is composed by the coefficients of the truncated Fourier expansion defined by equation (9). The function P qðqÞ; q v ðqÞ; q a ðqÞ ð Þ is defined as follows P qðqÞ; q v ðqÞ; q a ðqÞ ð Þ ¼ P I qðqÞ; q v ðqÞ; q a ðqÞ ð Þ þP E qðqÞ; q v ðqÞ; q a ðqÞ ð Þ where P I and P E are the penalty functions for the inequality and the equality constraints, respectively. To simplify notation, the arguments qðqÞ, q v ðqÞ, and q a ðqÞ will be omitted. The expression for P I is given by equation (12) If, for example, the inequality constraint k is a b, then I k is defined as a À b. When I k 0, then the term max I k ; 0 ð Þ is zero, and as a consequence e max I k ;0 ð Þ À 1 is also zero. When the inequality I k is not satisfied, then its contribution to the penalty function P I is a I e I k À 1 ð Þ. The expression for P E is given by equation (13) If, for example, the equality constraint k is a ¼ b, then E k is defined as a À b. The constant Tol is the tolerance for the equality conditions. When jE k j Tol, then the term e max E k j jÀTol;0 ð Þ À 1 is also zero. When jE k j > Tol, then the contribution to the penalty function P E is a E e E k j jÀTol À 1 À Á . Thus, the penalty-based optimization process requires a trial and error procedure to select a I , a E , and Tol. The values used to calculate the optimal trajectories were a I ¼ 10 9 , a E ¼ 10 6 , and Tol ¼ 10 À3 .
As the purpose of using a penalty approach was to avoid gradient estimates, equation (11) was solved using Matlab's fminsearch function, which is based on the Nelder-Mead simplex method. 44 The main criticism to this method is that no proof of convergence has been proposed. Nevertheless, we obtained better results using the penalty approach and fminsearch than using fmincon, which is the function for constrained optimization problems.

Experimental protocol
The comparison of the energy consumption between the cart-table and the optimization approaches is carried out at three different walking speeds: 0.5, 1.0, and 1.3 m/min. When the walking speed is greater than 1.3 m/min, the robot falls. This is because in our experimental setup, there is no closed-loop control that guarantees the balance of the robot. Energy consumption was recorded for 60 robot walks, 20 for each of the three walking speeds considered. For each speed, 10 walks were implemented using the carttable approach, and the other 10 using the optimization approach. In each trial, the robot walked in a flat and obstacle-free surface along a straight path 0.75 m long.

Results
The number of harmonics of equation (9) was selected by numerically solving the optimization problem (8) for different values of N. For N ¼ 1, no feasible solution was found. The values of the performance index (8) for N ¼ 2-7 are presented in Table 5. For N ! 5, the addition of harmonics to equation (9) does not represent a significant decreasing of equation (8), and by consequence we selected N ¼ 5.
From Table 6, it can be seen that for all joints the range of motion for the optimization method is smaller than that obtained with the cart-table approach.
The joint trajectories for the cart-table and the optimization methods are presented in Figures 5 and 6.
The energy consumption for the 60 trials is presented in Table 7. The last row of this table contains the cost of transportation (COT), which is calculated considering a mass of 1.7 kg and a distance of 0.75 m. From Figure 7, it can be seen that for the three considered speeds, the median of the energy consumption was lower for the optimization approach. At a speed of 0.5 m/min, the average energy for the optimization approach was 88.22% of the average energy for the cart-table method. For the other two speeds, these values were 91.61% (to 1.0 m/min) and 90.98% (to 1.3 m/min). Statistical analysis of experimental data can be performed using parametric or nonparametric methods. The former applies only when the data are normally distributed, and the latter does not require this assumption. The parametric method for comparing the means of two data sets is the Student's test, and its nonparametric equivalent is the Wilcoxon test. However, the Wilcoxon test does not compare means but medians. The Jarque-Bera test for normality was applied to the six data sets collected during the experiments (two trajectory generation methods and three gait speeds), and the results indicated that some of them are normally distributed but others are not. As a consequence, a nonparametric approach was used. The Wilcoxon test was applied to verify if the differences in energy consumption were significant. At each speed, the research hypothesis was that the energy consumption for the trajectory generated with the optimization approach was lower than that obtained by applying the cart-table method. By using the ranksum function of Matlab, the resulting p values were p < 0:001 for 0.5 m/min, p ¼ 0:0034 for 1.0 m/min, and p ¼ 0:0019 for 1.3 m/min. Figure 4 shows the joint angles as a function of time, for the two applied joint trajectory generation approaches, and considering a walking speed of 1 m/min.

The cost of transportation
During the experiments, we do not use the stabilizing mechanism provided by the gyroscope of the robot to avoid the corrections it introduces in the joint trajectories. In this manner, we compare the electrical energy expenditure of the original gaits generated by each method (cart-table and optimization) without introducing any bias in the experiment. Nevertheless, the main drawback of this open-loop motion is that the maximal gait velocity attained by the robot is only 1.3 m/min. Above this value, the robot falls at each trail. The consequence of the low gait speeds is the high COT we obtain for the Bioloid, which ranged from 28.88 to 72.39, depending on the gait speed and the trajectory generation approach. These values are by far greater than those reported for similar robots. For NAO, 18 the COT ranges between 2.4 and 5.8, and for the DARwIn-OP, 35 from 1.5 to 9.0. The high COT we obtain for the Bioloid is caused by the low gait speeds used during the experiments. The high energy consumption at low velocities has been reported for both humans and robots. In their study, Roberts et al. 35 present experimental results where the energy consumption at 0.025 m/s was nine times higher than that at 0.15 m/s.  Trial  Cart-table  Optimization  Cart-table  Optimization  Cart-table  Optimization   1  778  673  382  318  525  439  2  769  743  370  353  516  432  3  941  781  392  333  494  430  4  922  701  382  330  480  424  5  843  787  364  348  488  430  6  843  782  365  325  503  435  7  914  785  364  299  495  458  8  866  751  369  364  474  497  9  791  745  351  388  500  501  10  785  707  343  315  473  457  Mean  845  745  368  337  494   Another reason explaining the high COT obtained is the compliance margin and slope of the servomotors, which were configured to their maximal values to guarantee a precise tracking of the reference trajectory. Nevertheless, the higher the stiffness, the greater the current consumption. For example, in the case of the robot Aldebaran-NAO, 18 the experimental results show that the current supplied by the battery is almost double when the stiffness was changed from its minimum value to the maximum.

The sources of variability
The coefficient of variation (COV) of the COT takes values between 3.46% and 7.60% (see the row COV in Table 7). This is because during the experimentation we could not guarantee that the robot would walk on a perfect straight path, and as a consequence, the distance traveled by the robot was different from 0.75 m. Our COVs are consistent with those reported in similar experiments. In their study, Roberts et al. 35 reported COVs ranging from 5.9% to 12.4% for the energy expenditure of the DARwIn-OP humanoid robot. In the study by Ogura et al., 45 the electrical energy expenditure of the WABIAN robot is measured for three walking styles and two different step lengths (a total of six experimental configurations). The walking patterns were generated using a zero-moment-point-based approach, and the difference between gait styles lay in the trajectories of the knees. Only one trial was made for each gait, and as a consequence, it is difficult to assess whether the differences in energy consumption were statistically significant or not.

The benefits of the optimization
In the study of Bhounsule et al., 46 an optimization based on the dynamic model of the four-legged robot Ranger was used to generate its gait trajectory. This approach led to a reduction in the total COT from 0.28 to 0.19, when compared with trajectories intuitively tuned. This reduction of 32.14% is larger than that obtained in our experiments, which ranged from 9.16% to 13.35%. We hypothesized that these differences are due to the optimization approach. In their study, Bhounsule et al. 46 optimize the control signal applied to the robot, and in our case, the parameters used to describe the trajectory as a function of time.
The experimental results presented in this article show that for the three considered walking speeds, the average energy consumption obtained by applying the optimization approach is lower than that obtained with the cart-table method. From the results presented in Table 6, it can be seen that the main reason for these differences is the range of joint motion, which for the cart-table approach is always greater than that obtained by optimization. For the knee, a joint that can be an important source of energy consumption, 45,47 the optimization reduces the range of motion to 75.65% of the value obtained by the cart-table method.
The results of Wilcoxon test show that in all cases the differences are not due to chance but to the method used to generate the joint trajectories. To the best of our knowledge, an experimental comparison of electrical energy consumption between the cart-table and the optimization approaches has not been reported elsewhere.
Despite the encouraging results presented in this article, the optimization approach has some disadvantages when compared with the cart-table method. One of them is the high computational cost associated with solving the minimization problem described by equation (8). In the case of the trajectories presented here, 3-6 h of calculations were required on a laptop CORE-I7 with 32 GB of RAM to find a local minimum (the calculation time depended largely on the walking speed). Nevertheless, the main barrier to the adoption of the optimization approach is the need for a dynamic model of the robot, which is a procedure that requires vast experience and is prone to errors.
It is important to clarify that the limitation in the maximal velocity reduces the generality of our results. This is because the minimization of equation (8) implies the minimization of the electrical energy expenditure only when the relationship between the joint torques and the current supplied to the actuators is linear. When the joint velocities are increased, the equation that relates these two variables could become nonlinear, for example, due to magnetic saturation, 48 and consequently, the minimization of equation (8) does not guarantee energy expenditure minimization. In such a case, it could not be ensured that the optimal trajectories lead to a less energy expenditure than the carttable approach.

Future work
The electronic system presented in Figure 3 measures the energy expenditure of a given trajectory. Thanks to this system, we conclude that a trajectory minimizing the square torques leads to a lower energy expenditure than a trajectory generated by using the cart-table method. Nevertheless, for future works, we need to develop a current and voltage measurement system for each of the actuators to know which joints have the highest energy expenditure. In addition, we want to experimentally compare energy expenditure for innovative techniques for the generation of trajectories in biped robots such as the divergent component of motion 49 and the essential model. 50

Conclusion
The optimization approach requires a full mathematical model describing the dynamics of the biped. In practice, this model is just an approximation of the actual dynamic of the robot, neglecting the dynamics of the actuators and using an oversimplified representation of friction phenomena. Despite these limitations, the optimization approach leads to an important saving in energy consumption for low gait speeds.