Quadrupedal Robots’ Gaits Identiﬁcation via Contact Forces Optimization

: The purpose of the present paper is the identiﬁcation of optimal trajectories of quadruped robots through genetic algorithms. The method is based on the identiﬁcation of the optimal time history of forces and torques exchanged between the ground and the body, without any constraints on leg kinematics. The solutions show how it is possible to obtain similar trajectories to those of a horse’s walk but obtaining better performance in terms of energy cost. Finally, a map of the optimal gaits found according to the different speeds is presented, identifying the transition threshold between the walk and the trot as a function of the total energy spent.


Introduction
Nature has always been a source of inspiration for engineers and scientists who tried to replicate or mimick natural bio-mechanisms. In fact, this can be a smart strategy, since any natural living being is a highly optimized system over millions of years of evolution [1][2][3]. Magnificent examples of engineering bio-inspired devices date from very ancient Magna Greece times, as the Archita's dove automaton, and a large number of automata within the Erone's tradition. The masterful Leonardo's machines, such as the lion or the knight automaton, or the flying machines, some inspired by the flapping wings of birds, are also astonishing examples.
To date, the development of automatic systems, in particular quadrupedal robots, has found great motivation in the field of human assistance. To assist and cooperate with humans, robots must be able to move easily in workplaces originally designed according to human needs [4]. For example, climbing stairs, avoiding obstacles, or opening doors [5,6]. Among the many technical difficulties to be faced and solved, one of the most important aspects is related to portability and its duration. This scenario justifies the need to optimize the power and energy of locomotion to make these systems more efficient, compact, characterized by limited energy consumption, and using small size and weight actuators.
Interesting questions arise: is it possible to conceive more efficient quadrupeds than those proposed by nature? Is it possible to disclose a different kind of gaits for a quadrupedal mechanism as the product of a strict optimization process? For example, depending on speed, energy consumption, or jerk minimization. Many engineers and robot designers assume the gait of animals is optimal since they would have been able to survive the competition and natural selection [7][8][9][10]. For this reason, the quadruped robot tends to be designed inspired by the examples of nature [11,12]. However, biological kinematics of locomotion is not reproducible directly by a legged robot, since biological muscles are extraordinarily efficient components, and animals can produce energy for their actuation by grabbing nutritive elements from the environment, transformed into actuation energy by a digestive process. This permits them to have a relatively light energy storage system. Moreover, artificial sensors are still inferior, in terms of number, miniaturization, and sensitivity against the ones owned by animals [13].
To date, designers tend to create robots that resemble nature as closely as possible, such as imposing specific gaits such as walking, trotting, or galloping and developing mathematical optimization processes that only concern sub-systems, to overcome the technological limitations [14,15]. There are many different approaches for performing the optimization. In [16], for example, the authors use harmonic functions to describe foot trajectory and tune some parameters such as amplitude, phase, and frequency. In [17,18] the authors generate leg trajectories for overcoming small terrain irregularities by optimizing consumed energy. In [19] the authors optimize contact forces to walk on a high-slope terrain, but the gait is imposed a priori.
Evolutionary computation, such as genetic algorithms (GAs), is often a natural choice for the gait optimization of legged robots since it uses optimization methods based on evolutionary theory [20]. Also, the controls of the mentioned cases, tend to mimic or approach the natural movement of the quadrupeds found in nature. Typical controls are central pattern generators [21,22], hybrid zero dynamics [23], or the usual PD [24].
The present work breaks from all the problems relating to technological limitations and lays the foundations of how it is possible to identify, a priori, optimal gaits through parametric optimization algorithms. The authors propose an algorithm that does not depend either on the kinematic chains of the legs or power actuators. The method directly handles the contact forces that the hypothetical legs exchange with the ground, identifying the sequence that guarantees minimum energy consumption and smooth motion during a single stride cycle. The genetic algorithm (GA) is used to minimize the kinetic and potential energy as well as restraining the jerk, identifying the optimal profile of normal and friction forces. The equations of motion are solved through some periodicity constraints of the gait, thus ensuring a stable and continuous motion. The optimal gaits identified are very similar to those found in nature. The authors have optimized the gait of a heavy and bulky horse-like robot, obtaining a walking trot as the optimal gait, and then comparing it with experimental data from real horses. Additionally, a gait transition analysis is performed as a function of locomotion speed and gait cycle length, showing why it is convenient to walk or trot.
The paper is organized as follows: Section 2 describes the dynamic model and ground reaction force shapes; Section 3 defines the optimization parameters and objective functions while introducing stability constraints to assure a periodic motion; Section 4 provides the resulting optimal gait, along with the comparison with horse experiments and identification of gait transition; Section 5 discusses the limitations and future challenges of the method and finally the conclusions.

Mathematical Model for the Gait Optimization
The considered quadrupedal model is sketched in Figure 1, where the body, suspended on four legs, transmits forces and moments thanks to the interaction with the ground.
The optimization model proposed in this paper consists of an optimal gait identification, capable of moving the quadrupedal body through the succession of four alternating thrusts generated by legs with some desirable characteristics. This approach is innovative. In fact, the kinematic linkages and actuators of the leg mechanism are not directly involved but left to a second engineering phase of design, not investigated in the present paper.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 3 of 2 ℎ ( ), while the one during which the leg is flying in the air movin towards the next contact point, is the ℎ ( ). In general, these times can be se differently for each leg.  Tait-Bryan angle  (yaw pitch and roll), is the angular velocity, and is the Jacobian matrix. Assuming small pitch and roll angles and vanishing yaw, the equations can b reduced into the same fixed reference frame: , are the total force and moment transmitted by the four legs , , and , respectively. When the quadruped moves on the ground, the leg-ground contact produces reaction force. Its vertical component supports the quadrupedal weight, the horizonta component drives the longitudinal motion of the body, and affects the lateral bod oscillations.
The vertical force for each leg acts only in the ; that implies ( ) = 0 and ( + ) = 0, that is, its time history, during the , starts and ends with vanish ing values, associated with the incipient contact and incipient detachment of the leg wit and from the ground, respectively. The force-time history along the is identified The legs are labelled as FL (front left), FR (front right), HL (hind left), HR (hind right). The period during which a single leg maintains its contact with the ground is the conctact phase (CP), while the one during which the leg is flying in the air moving towards the next contact point, is the swing phase (SP). In general, these times can be set differently for each leg.
The legs sequencing, sketched in Figure 2, originates a periodic motion of each leg, of period T. The time at which each leg touches the ground is t i ∈ [0, T], where i = FR,FL, HR, HL. The time duration of the contact phase is T i , and the time duration of the swing phase is T − T i .
Newton-Euler equations of rigid body dynamics are: where the first equation holds in the fixed reference frame, the second in the body reference frame (principal inertia axes), r = [x, y, z] identifies the gravity center position of the quadruped body of total mass m, f e is the total external force, I = diag I x , I y , I z the matrix of inertia, m e is the total external moment, θ = [ψ, θ, φ], the Tait-Bryan angles (yaw pitch and roll), ω is the angular velocity, and E is the Jacobian matrix.
Assuming small pitch and roll angles and vanishing yaw, the equations can be reduced into the same fixed reference frame: where F x , F y , F z , M x , M y are the total force and moment transmitted by the four legs FR, FL, HR, and HL, respectively.
When the quadruped moves on the ground, the leg-ground contact produces a reaction force. Its vertical component F z supports the quadrupedal weight, the horizontal component F x drives the longitudinal motion of the body, and F y affects the lateral body oscillations.
The vertical force F z i for each leg acts only in the CP; that implies F z i (t i ) = 0 and F z i (t i + T i ) = 0, that is, its time history, during the CP, starts and ends with vanishing values, associated with the incipient contact and incipient detachment of the leg with and from the ground, respectively. The force-time history along the CP is identified through three additional points P 1 i ; P 2 i ; P 3 i (besides the two previously identified) and is fitted by a suitable spline (Figure 3a).  The longitudinal and lateral ground reactions for each leg are assumed to belong t the static friction cone. With this assumption, any contact slip condition is absent, and th contact model is purely conservative. Therefore: where is the static friction coefficient. For the planar forces and , one can introduce suitable constitutive relationship ( Figure 3). Several experiments available in the technical literature [25,26] suggest the fol lowing gait dynamics for leg locomotion. When the leg touches the ground, in the firs phase, a backward longitudinal reaction brakes the body run, followed by a second phase in which the leg accelerates the body through a forward longitudinal reaction force ( Figure 3d). Moreover, the amplitude of the longitudinal reaction force is roughly propor tional to the normal vertical force. A similar force trend holds for the lateral force that i responsible for the unavoidable lateral staggering of the quadruped's body. For both lon gitudinal and lateral forces, one can introduce the suitable factorization form as constitu tive relationships:  The longitudinal and lateral ground reactions for each leg are assumed to belong to the static friction cone. With this assumption, any contact slip condition is absent, and the contact model is purely conservative. Therefore: where µ s is the static friction coefficient. For the planar forces F x i and F y i , one can introduce suitable constitutive relationships ( Figure 3). Several experiments available in the technical literature [25,26] suggest the following gait dynamics for leg locomotion. When the leg touches the ground, in the first phase, a backward longitudinal reaction brakes the body run, followed by a second phase, in which the leg accelerates the body through a forward longitudinal reaction force ( Figure 3d). Moreover, the amplitude of the longitudinal reaction force is roughly proportional to the normal vertical force. A similar force trend holds for the lateral force that is responsible for the unavoidable lateral staggering of the quadruped's body. For both longitudinal and lateral forces, one can introduce the suitable factorization form as constitutive relationships: that use the proportionality of the planar forces to the vertical one. Moreover, the role of the factors µ x i (t), µ y i (t) is that of reproducing the time history of the planar force, F x i (t) or F y i (t), as reported in Figure 3c,d, modulating the amplitude of the vertical reaction F z i (t).
In particular, µ x i (t), µ y i (t) are vanishing outside of the CP. Within the CP, in the first part of the contact, µ x i (t) takes negative values to reproduce the braking effect, while, in the second part, it becomes positive, to reproduce the acceleration effect. Appl. Sci. 2021, 11, x FOR PEER REVIEW 5 of 22 Combining Equations (4) and (3), the two modulating factors must satisfy the condition: Looking at Figure  ; selected by the designer within the contact cone (see Figure 3b). The expressions for the moments and follow directly from the leg's reactions. From Figure 4, let's consider, for example, the foot at the incipient contact phase and foot position ( , , ): Combining Equations (4) and (3), the two modulating factors must satisfy the condition: Looking at Figure 3b, this implies that at any time the point Q j i µ x j i (t), µ y j i (t) remains confined in the circle of radius µ s over the plane µ x i (t), µ y i (t). Moreover, we assume a correlation between the modulating factors during the gait. This implies, concerning Figure 3b, as time is spent, and Q completes the loop in a period T. More precisely, the time histories of µ x i (t), µ y i (t) are determined by a spline passing through the points Q 1 i ; Q 2 i ; Q 3 i selected by the designer within the contact cone (see Figure 3b).
The expressions for the moments M x and M y follow directly from the leg's reactions. From Figure 4, let's consider, for example, the FR foot at the incipient contact phase and foot position (x FR , y FR , z FR ): where ∆x FR (t) = x FR − x, ∆z FR (t) = z FR − z and ∆y FR = y FR − y.

Optimization Model
Usually, in the legged locomotion, the stability of a gait is guaranteed by criterion such as zero-moment point [27,28]. Other studies face the stability probl the Poincare map [29,30] or ground reference points [31]. In this paper, unlike the approaches, the gait stability is guaranteed by satisfying periodic limit cycle cond The optimization operates over a single period using a genetic algorith which provides the stride length, frequency, and velocity, and additionally the t tory of the contact forces. Therefore, a vector of gait parameters is optimized the GA algorithm. Let us explicitly set . The relative phases of the leg motio the normalized time duration of the CP is = , the time history of the norm gitudinal and lateral contact forces are determined by the interpolation of the p and with = 1,2,3, respectively. Therefore: The structure of the optimization process involves the optimal vector an several constraints derived from the equations of motion, it looks for solutions th mize a multi-objective cost function related to the dissipated energy and loo smooth forces trends. Let us examine in detail these constraints.
Note, as a preliminary consideration, that we can impose the continuity k conditions over the stride period: We start with the vertical equation of motion (an identical procedure appli other equations of motion) that produces:

Optimization Model
Usually, in the legged locomotion, the stability of a gait is guaranteed by using a criterion such as zero-moment point [27,28]. Other studies face the stability problem with the Poincare map [29,30] or ground reference points [31]. In this paper, unlike the classical approaches, the gait stability is guaranteed by satisfying periodic limit cycle conditions.
The optimization operates over a single period T using a genetic algorithm (GA) which provides the stride length, frequency, and velocity, and additionally the time history of the contact forces. Therefore, a vector of gait parameters p GA is optimized by using the GA algorithm. Let us explicitly set p GA . The relative phases of the leg motion are t i , the normalized time duration T i of the CP is β i = T i T , the time history of the normal, longitudinal and lateral contact forces are determined by the interpolation of the points P j i and Q j i with j = 1, 2, 3, respectively. Therefore: The structure of the optimization process involves the optimal vector p GA and, using several constraints derived from the equations of motion, it looks for solutions that minimize a multi-objective cost function related to the dissipated energy and looking for smooth forces trends. Let us examine in detail these constraints.
Note, as a preliminary consideration, that we can impose the continuity kinematic conditions over the stride period: x(T); .
Moreover, the average speeds over the period T are where V is the average longitudinal speed of the quadruped body.
We start with the vertical equation of motion (an identical procedure applies to the other equations of motion) that produces: where g is the gravity acceleration and C is an integration constant. Therefore: .
and .
Using the condition .
z(T), it follows: The force-time history is set by randomly selecting the point P j i by the GA and using a spline interpolation. Further integration of Equation (10) produces: . and: Using the periodicity condition z(0) = z(T), it follows: .
Equation (15) automatically satisfies the condition of zero average speed . z = 0. In fact: .
The longitudinal dynamics takes into account two inertia effects: one related to the longitudinal motion of the robot body mass m B , and the second due to the reciprocating motion of the legs of mass m p i . The total robot mass m is consequently: Let us make the point about the leg motion. Looking at Figure 5, it appears that, during the contact phase, the horizontal speed of the leg is negative, and becomes positive during the swing phase. Therefore, along the gait period, the leg mass behaves as an oscillator that, attached to the primary body mass, generates a back and forth horizontal inertia force.
Let us make the point about the leg motion. Looking at Figure 5, it appears that, during the contact phase, the horizontal speed of the leg is negative, and becomes positive during the swing phase. Therefore, along the gait period, the leg mass behaves as an oscillator that, attached to the primary body mass, generates a back and forth horizontal inertia force.  The mathematical counterpart of this physical behaviour appears as follows: .. ..
x sw i (18) x sw i describes the leg mass longitudinal motion during the swing phase. The term ..
x sw i takes into account the inertia force during the swing phase, while, complementarily, m p i w i (t) ..
x accounts for the added mass effect during the contact phase, within which the leg motion depends essentially on the body motion. The switch function w i (t) is defined as: that activates alternatively the inertia effects of the contact and swing phase, respectively. The swing motion x sw i (t) is designed through β i and T by using a spline with initial and final slope equal to the average body velocity V (as an example see Figure 6).
Appl. Sci. 2021, 11, x FOR PEER REVIEW 8 of 22 The mathematical counterpart of this physical behaviour appears as follows: describes the leg mass longitudinal motion during the swing phase. The term 1 − ( ) takes into account the inertia force during the swing phase, while complementarily, ( ) accounts for the added mass effect during the contact phase within which the leg motion depends essentially on the body motion. The switch function ( ) is defined as: that activates alternatively the inertia effects of the contact and swing phase, respectively The swing motion ( ) is designed through and by using a spline with initia and final slope equal to the average body velocity (as an example see Figure 6). The analogous Equation (12) is determined for the longitudinal motion using the condition (0) = ( ): Then the average speed ̅ = over the period T is: that, together with the periodicity condition, produces: The analogous Equation (12) is determined for the longitudinal motion using the condition .
Then the average speed .
x = V over the period T is: .. that, together with the periodicity condition, produces: . ..
If the same procedure is applied to the lateral motion, it follows: derived under the condition .
We can apply the same procedure to the remaining cardinal equations. Using the conditions .
φ(T), the following constraints are obtained: These can be solved in closed form in terms of the initial conditions x 0 e y 0 . The arms of the moments are: Assuming a flat ground, z i = 0, and Finally, using again the same technique with the conditions φ(0) = φ(T), θ(0) = θ(T), the final equations yield: .
The previous constraints are based on the linear approximation (see Equation (2)) of the equations of motion, and they provide the initial guess of the kinematic variables at the initial and final times, 0 and T. Moreover, the force and the moments, through their spline approximations, satisfy the set of conditions from (7) to (27). With this starting guess, we can solve the following iterative nonlinear optimization problem: The unknown vector to be determined is p NLP = θ 0 , φ 0 , ψ 0 , .
. ω z 0 , x 0 , y 0 , z 0 . Figure 7 illustrates the flow chart of the GA solution's scheme. The setting of the GA which led to good balancing between learning and overfitting is a size population of 50 individuals, 110 generations, a crossover of 80% and mutation of 1%. Specifically, the genetic algorithm with multi-objective optimization is used for which simultaneous optimization of multiple, often competing, objectives take place [32,33]. The related pseudocode is shown below: Algorithms 1 Pseudo code for the gait optimization algorithm.

1:
t ← 0 2: Initialization parameters T, V, m, I x , I y , I z , L, h, m p 3: GA random population initialization of p GA vector 4: While {GA's stop criterion is not met} 5: Vertical ground force definition with continuity kinematic constraints (periodicity condition) 6: Longitudinal and lateral forces with continuity kinematic and static friction cone constraints 7: Torques with continuity kinematic constraints in the linear case (small angles) 8: Initialization guess parameters for Nonlinear Programming 9: Evaluate Objective function for Nonlinear Programming 10: While {Nonlinear Programming's stop criterion is not met and the constraints are satisfied} 11: Sequential Quadratic Algorithm 12: Evaluate Objective function 13: Endwhile 14: Evaluate fitness function for each individuals 15: Select of the individuals 16: Recombine of the individuals 17: Mutate of the individuals 18: Generations of new population 19: t ← t + 1 20: Endwhile Eventually, we introduce the multi-objective function to be minimized for the most efficient gait. The GA's fitness function is defined through the total specific energy , Eventually, we introduce the multi-objective function to be minimized for the most efficient gait. The GA's fitness function is defined through the total specific energy E d , that is, the energy per traveled distance x(T) − x 0 over the period, and the ground reaction jerk J GRF .
The specific energy E d includes the mechanical energy of the body and the swinging leg energy: ..
One can expect power fluctuations along the period with positive and negative values, perfectly balanced since these fluctuations have zero average (the system is conservative). Since E d integrates the absolute value of the power, it measures the intensity of the power fluctuations. Requiring a small value for E d means asking for small fluctuations of the power, which implies a benefit in terms of the actuator dimensions and power and thus of energy consumption. On the other hand, the request of a small jerk J GRF assures traditionally smooth contact forces: where F = F x ; F y ; F z . Multi-objective optimization is used to determine the Pareto frontier, the set of solutions that provides the optimal trade-off between the two criteria.

Optimal Gait Solution
For certain optimization variable sets p GA , the solver may not converge, or the solution may be unrealistic or impractical. A maximum foot stride s max has been imposed to avoid solutions that required unrealistic leg dimensions. The condition was applied to the upper bound of the duty factor: β i max = s max VT , where V is the mean target speed of the simulation: The optimization finds an efficient gait (in the sense specified in the previous section) that moves the body at an average speed of 1.35 m/s, with a gait period T of 1 s. The mass properties of the body are selected considering the characteristic parameters of quadrupeds in nature, in particular horses (Table 1, see reference [34]). In addition, some simplification but reasonable assumptions are made, namely: The obtained solutions are plotted in Figure 8, where a correlation between the energy cost and contact-forces jerk is shown (Pareto frontier). The resulting gaits present a specific leg sequence classified as walk and trot gaits. In nature, walking gaits involve overlapping contact phase such that β > 0.5. In particular, in a trot gait, the diagonal forelimb and hindlimb move in phase, while in a walking gait the limbs have a lag of 0.25 between the two legs. The best solution according to the Pareto frontier is presented in Figure 9. It can be seen how the body keeps a steady longitudinal speed ̇, about the target speed of = 1.35 / , maintaining a stable and restrained attitude and assuring the periodicity constraints. The best solution according to the Pareto frontier is presented in Figure 9. It can be seen how the body keeps a steady longitudinal speed .
x, about the target speed of V = 1.35 m/s, maintaining a stable and restrained attitude and assuring the periodicity constraints. The best solution according to the Pareto frontier is presented in Figure 9. It can be seen how the body keeps a steady longitudinal speed ̇, about the target speed of = 1.35 / , maintaining a stable and restrained attitude and assuring the periodicity constraints. This case is classified as a walking trot since it presents a duty cycle β FR = β FL = β HR = β HL = 0.55, equal for the hind and forelimbs, and t i = [0.7; 0.2; 0.2; 0.7] (Figure 10). The contact forces are instead shown in Figure 11, and the friction coefficients in Figure 12. It can be observed how the normal forces F z i have a flat double-hump time history, typical of quadrupeds' locomotion [2,35]. All other significant parameters are shown in Table 1.        In Figure 13 we can see the trajectory performed by the foot FR during one complete cycle.
In Figure 13 we can see the trajectory performed by the foot FR during one complete cycle. To better assess the variability of the results, a statistical analysis on 10 runs of the optimization algorithm was performed. In all cases, the best gait found was the walkingtrot and the ratio between the value of the objective function and its average differs by a maximum of 2.5%, as observed in Figure 14. Eventually, a comparison between the cumulative mechanical energy of the experimental and numerical optimization data is shown in Figure 15. For the experimental campaign [34], a group of horses with a weight ranging between 417 and 673 , was used. For each animal, ten sequences were considered within the speed range from 0.7 to 1.9 / . Five markers, positioned along the horses' backbones, identify the vertical and longitudinal motion of their center of mass ( ). The considered mechanical energy is defined by the expression (a special case of expression (29)): To better assess the variability of the results, a statistical analysis on 10 runs of the optimization algorithm was performed. In all cases, the best gait found was the walkingtrot and the ratio between the value of the objective function and its average differs by a maximum of 2.5%, as observed in Figure 14. To better assess the variability of the results, a statistical analysis on 10 runs of the optimization algorithm was performed. In all cases, the best gait found was the walkingtrot and the ratio between the value of the objective function and its average differs by a maximum of 2.5%, as observed in Figure 14. Eventually, a comparison between the cumulative mechanical energy of the experimental and numerical optimization data is shown in Figure 15. For the experimental campaign [34], a group of horses with a weight ranging between 417 and 673 , was used. For each animal, ten sequences were considered within the speed range from 0.7 to 1.9 / . Five markers, positioned along the horses' backbones, identify the vertical and longitudinal motion of their center of mass ( ). The considered mechanical energy is defined by the expression (a special case of expression (29)): Eventually, a comparison between the cumulative mechanical energy of the experimental and numerical optimization data is shown in Figure 15. For the experimental campaign [34], a group of horses with a weight ranging between 417 and 673 kg, was used. For each animal, ten sequences were considered within the speed range from 0.7 to 1.9 m/s. Five markers, positioned along the horses' backbones, identify the vertical and longitudinal motion of their center of mass (CoM rigid ). The considered mechanical energy is defined by the expression (a special case of expression (29)): It is clear how the resulting energies are comparable.

Gaits Analysis for Different Speed and Cycle Duration T
The optimization solutions are tested for different speeds and cycle period . Figure 16 depicts the trend of the walking gait representing , following Equation (30), as dependent on the cycle duration, for different speeds between 0.7 and 2 / . This map permits finding natural associations between convenient gait period and body velocity, to obtain gaits with a reasonably limited energy cost and smooth forces.
We notice that for < 0.5 the total energy is very large, because the swing phase has a very limited duration, needing very high energy to relocate the foot at its initial position. On the other hand, at high values for , again the energy cost increases. Therefore, the map suggests that for any desired speed, there is an intermediate region for the gait period that allows a low energy consumption. Minimum energy expenditures are identified at = 1.25 for speed 1.35 / , at = 1 for speed 2 / , and = 1.66 for speed 0.7 / (blue square, green triangle, and red diamond, respectively). In Figure 17, the trend of the trot gait is depicted. We notice the same behavior of low energy cost for the walking gait in the middle region, with an optimal association between It is clear how the resulting energies are comparable.

Gaits Analysis for Different Speed and Cycle Duration T
The optimization solutions are tested for different speeds and cycle period T. Figure 16 depicts the trend of the walking gait representing E d , following Equation (30), as dependent on the cycle duration, for different speeds between 0.7 and 2 m/s. This map permits finding natural associations between convenient gait period and body velocity, to obtain gaits with a reasonably limited energy cost and smooth forces.
We notice that for T < 0.5 s the total energy is very large, because the swing phase has a very limited duration, needing very high energy to relocate the foot at its initial position. On the other hand, at high values for T, again the energy cost increases. Therefore, the map suggests that for any desired speed, there is an intermediate region for the gait period that allows a low energy consumption. Minimum energy expenditures are identified at T = 1.25 s for speed 1.35 m/s, at T = 1 s for speed 2 m/s, and T = 1.66 s for speed 0.7 m/s (blue square, green triangle, and red diamond, respectively). It is clear how the resulting energies are comparable.

Gaits Analysis for Different Speed and Cycle Duration T
The optimization solutions are tested for different speeds and cycle period . Figure 16 depicts the trend of the walking gait representing , following Equation (30), as dependent on the cycle duration, for different speeds between 0.7 and 2 / . This map permits finding natural associations between convenient gait period and body velocity, to obtain gaits with a reasonably limited energy cost and smooth forces.
We notice that for < 0.5 the total energy is very large, because the swing phase has a very limited duration, needing very high energy to relocate the foot at its initial position. On the other hand, at high values for , again the energy cost increases. Therefore, the map suggests that for any desired speed, there is an intermediate region for the gait period that allows a low energy consumption. Minimum energy expenditures are identified at = 1.25 for speed 1.35 / , at = 1 for speed 2 / , and = 1.66 for speed 0.7 / (blue square, green triangle, and red diamond, respectively). In Figure 17, the trend of the trot gait is depicted. We notice the same behavior of low energy cost for the walking gait in the middle region, with an optimal association between In Figure 17, the trend of the trot gait is depicted. We notice the same behavior of low energy cost for the walking gait in the middle region, with an optimal association between periods and speed: T = 1.25 s for 1.35 m/s, T = 1 s for 2 m/s, and T = 1.66 s for 0.7 m/s.

/ .
Then, taking the minimum points for each gait and comparing them with the ener cost and jerk performances (Figure 18), we see how the transition between walking an trotting takes place at a speed between 0.7 / and 1.35 / , by using a Pareto criterio In fact, for the lowest speed, 0.7 / , the walk is optimal, in a Pareto sense, since ener cost is lower, and jerk substantially similar. Increasing the speed at 1.35 / , a trot pe forms much better, since both energy and jerk indexes are lower with respect to the wal ing gait. For higher speeds, the trend of the two curves shows the trot gait is better th walking.   Then, taking the minimum points for each gait and comparing them with the energy cost and jerk performances (Figure 18), we see how the transition between walking and trotting takes place at a speed between 0.7 m/s and 1.35 m/s, by using a Pareto criterion. In fact, for the lowest speed, 0.7 m/s, the walk is optimal, in a Pareto sense, since energy cost is lower, and jerk substantially similar. Increasing the speed at 1.35 m/s, a trot performs much better, since both energy and jerk indexes are lower with respect to the walking gait. For higher speeds, the trend of the two curves shows the trot gait is better than walking. Then, taking the minimum points for each gait and comparing them with the energy cost and jerk performances (Figure 18), we see how the transition between walking and trotting takes place at a speed between 0.7 / and 1.35 / , by using a Pareto criterion. In fact, for the lowest speed, 0.7 / , the walk is optimal, in a Pareto sense, since energy cost is lower, and jerk substantially similar. Increasing the speed at 1.35 / , a trot performs much better, since both energy and jerk indexes are lower with respect to the walking gait. For higher speeds, the trend of the two curves shows the trot gait is better than walking.   Figure 18. Gait transition from walk to trot based on Pareto's frontier analysis.

Discussion
Based on the previous models and results, we can conjecture the optimal gait observed in nature can be derived as an evolutionary effect that combines minimum energy requirements and minimum jerk contact forces, to obtain a gait that is energy-saving and with a moderated mechanical impact on the legs articulations. This point of view has importance also for a quadrupedal robot design that would be aimed at imitating or even enhancing nature's results.
Nevertheless, the previous analysis introduces several approximations with related inaccuracies in the representation of quadrupedal mechanics.
As a first remark, dissipation effects are not taken into considerations. Credible dissipation models need indeed a fine-tuning and require specific and non-trivial experimental tests. This approximation introduces larger errors when increasing the leg's speed since the non-modeled effects of dissipation increase with some power of the velocity. Furthermore, one of the advantages of the previous analysis is the absence of a detailed leg model, that is on the other hand a limitation in the accuracy of the model. Additionally, the absence of a head in the model and its coordinated motion with legs is another approximation when analyzing the quadruped gaits.
Finally, the presence of flexible and deformable elements in the animal mechanical structure is not considered in the model, and these elements absorb elastic energy that plays a role in the global energy balance.
Actuators and their related power are not explicitly included in the animal mechanical model. Besides their obvious importance in the energy balance in the robot design, they guarantee the optimal gait is suitably tracked by the leg motion, with no obvious implications for their performances in terms of supplied maximum force and power, weights, and volumes.
Finally, the contact model is an approximate empirical model, and it will be necessary to check its reliability on the experimental ground.
Future work challenges will be related to overcoming these limitations. More specifically related to quadruped robot design implications, the kinematics and dynamics of leg implementation will be further optimized to be able to pursue the trajectories determined by the optimal gait generator. Non-linear optimal control algorithms recently developed by the authors [36,37] will be used to track kinematic references and forces exchanged on the ground. Such algorithms can be combined with reference trajectory generators [38,39] but can also incorporate the dynamics model and the power actuators by variational feedback control methods [40,41]. Therefore, it is possible to verify and confirm the validity of the results and then extend the work to the optimization of speed changes, changes of gaits, changes of direction, and optimal handling in the presence of obstacles and irregularities of the terrain [42,43].
The present analysis lays the groundwork for the design of new highly efficient quadrupedal robots.

Conclusions
This paper presents an investigation related to the optimal gaits for a quadrupedal robot. The interest of this analysis clearly relies on the use of such legged robots in many engineering areas, where the replacement of wheels by legs is highly desirable in an environmental scenario that has an uneven ground with obstacles and in the absence of specific infrastructure that permits a wheeled robot to operate comfortably.
The analysis presents several elements of originality. The optimal gait problem is formulated independently of any specific architecture for the leg mechanisms. The focus is the leg-ground reaction forces that control the quadruped's body motion. A detailed analysis of the contact introduces constitutive relationships for the reaction forces that permit their use in the formulation of the dynamic equation of the body. With the help of the periodicity conditions, that characterize the state variables of the robot involved in its gait and using them an initial guess based on a linearization of the equation of motion, a GA optimization algorithm produces the optimal gaits of the quadruped. A multiobjective function introduces two basic principles of optimization. The first relies on the minimization of the maximum actuators' power and energy consumption. The second relies on the minimization of the contact jerks; that implies a requirement for preventing overload and stressful conditions for the robot mechanisms and its kinematic joints. Their combination, subjected to the Pareto frontier analysis, suggests optimal gaits.
The method is applied to a quadruped the size and mass of which are those of a horse, and the optimal gaits obtained by this optimization method are compared with results found in the technical literature for a real animal. The comparison shows that the different gaits, such as walk and trot, and the convenient transition from the first to the second, depends on the desired speed of the quadruped. The need for a transition comes out as a requirement of energy and jerk optimization; when the speed overcomes a given threshold (located between 0.7 m/s and 1.35 m/s) it is clearly convenient to change the gait. One may speculate this result could explain the origin of the quadruped gaits, as driven by nature, under a very long evolutionary process.

Author Contributions:
The optimal gaits of a quadruped robot are part of the Ph.D. Thesis of M.L., under the guidance of G.P. (tutor). The methodology and conceptualization are by G.P. and M.L.; software, and validation, are by M.L. and G.P.; initial suggestions and some conceptualizations about optimal gait problem are by A.C., supervisor of the thesis project, who edited the writing of the manuscript and by N.P.B., former professor at Sapienza, now at University Roma 3. All authors have read and agreed to the published version of the manuscript.