2. State of the art

In most cases, the humanoid robot design and walking pattern generation are investigated separately. Some researchers stated that walking simulations were performed at the beginning of robot structure design. For example, Xia [Reference Xia, Liu, Jing, Qiang and Chen15] analysed the effects of the distance between the two hips and the height of the centroid on walking stability, and they attempted to optimize the leg morphology based on the simulation result. However, to the current knowledge of the authors, there are only a few coordinated morphology-locomotion optimization trials on humanoids.

The first such result can be found in Endo [Reference Endo, Maeno and Kitano16], where an EA is used to coevolve limb length and a neural-network-based controller for walking to avoid trial-and-error in mechanical design and motion generation. A robot that is capable of walking for the longest distance is finally found, and the result is used for designing the real bipedal walking robot OpenPINO. The robot with the best energy efficiency is also coevolved by EA/GA [Reference Endo, Yamasaki, Maeno and Kitano17]. These two pioneering studies brought the concept into being; however, the relationships among morphology, locomotion and the working environment have not been clearly investigated for living systems. The method appears to be limited to the particular robot and the neural networks and oscillators, and the GA/EA adopted is a computationally intensive method for optimizing such problems.

Eaton [Reference Eaton18] also uses a GA to generate the walking motion of bipeds. This research also varies some of the joints of the leg, and environmental factors are also considered to compare the different walking styles. Jouandeau [Reference Jouandeau and Hugel19] only optimized the separation between the two hips based on walking performance.

The research on passive walking [Reference Collins and Ruina20] is most likely to be a type of coordinated implementation of mechanical structure and motion generation. The mechanical structure is well optimized to realize a walking style with minimal power input. However, the thought behind the passive walking style appears to be a better motion generation primarily by optimizing the mechanical structure. Here, the active motion generator is less focused or has few chances to impact the walking quality. Therefore, the scheme of passive walking cannot be classified as a coordinated evolution of mechanical structure and motion generation.

Ref. [Reference Nolfi, Bongard, Husbands and Floreano21] provides a review of evolutionary robotics. The core is to develop various robot systems and control policies by GA. EA/GA is used to develop both the control system and body configuration in close interaction with the environment without human intervention. Behaviours for obstacle avoidance in complex circumstances, bipedal walking by neural networks, the coevolution of robot bodies and brains and vision and navigation are all developed in the framework of EAs.

Within all these studies, GA/EA plays a vital role. It generates the new candidate points in a stochastic manner, and it makes the optimization more likely to be similar to the evolution of biological species [Reference Endo, Maeno and Kitano16]. There are also several other stochastic-based approaches, such as particle swarm, pattern search and simulated annealing. All these approaches share several common features, as follows. Every execution of the optimization may follow a different path, and many parallel function evaluations are needed in every step towards convergence. Points previously evaluated do not participate in the later candidate generation. Because dynamics simulation for multiple objects with a relatively large amount of experiments can be computationally expensive, those evolutions would generate a large number of individuals for simulations to obtain a better machinery-motion couple.

In addition to these stochastic methods, there are also some global optimization methods that can reach the same goal. For example, design and analysis of computer experiments (DACE), established by Sacks [Reference Sacks, Welch, Mitchell and Wynn22, Reference Jones, Schonlau and Welch23], is widely used as the standard approach for deterministic black-box optimization problems. It belongs to a special type of method called a response-surface-based approach. Rather than the stochastic-based genetic and EAs, the response-surface-based approach does not perform optimization directly on the objective model. It first builds a surrogate model, which is actually the interpolation of existing experimental points. The interpolated model is actually an available derivative approximation of the target function, and the optimized point can be obtained using classical non-linear optimization methods. The surrogate model approach is preferable because of its far fewer demands in terms of experimental evaluation. A stochastic process model called the Kriging model is used to approximate the behaviour of the objective function with high non-linearity. Another kind of global optimization methods are walking pattern generation using reinforcement learning. Ref. [Reference Lee and Oh24] used reinforcement learning to find the proper initial position and velocity and the final position and velocity. The stable biped walking pattern is verified through a simulation. Silva [Reference Silva, Perico, Homem, Vilo and Bianchi25], [Reference Silva, Perico, Costa and Bianchi26] used reinforcement learning to improve the walking stability of humanoid robot on sloped terrain and optimize the parameter values for the gait pattern generation with temporal generalization. The methods based on reinforcement learning can optimize the gait parameters and are possible to identify the relationship between the parameters, which have great development potentials. However, the learning time is still very long at current stage.

There are also a few applications of finding the bipedal walking parameter for some of the approaches. Hemker [Reference Hemker, Stelzer and Stryk27] uses the Kriging model, a widely used surrogate model, to enhance the bipedal walking speed through sequential experiments on a physical robot. Our previous work seeks to improve the omnidirectional walking performance of this type of walking [Reference Zhang, Xia, Liu and Chen28].

In general, there is still no research that puts the mechanical structure and walking motion of humanoid robots into a surrogate model for collaborative optimization. The goal of this paper is to explore the relationship between structure and walking motion and to provide a systematic design guidance for humanoid robots.

3. Framework of humanoid robot evolution

In Fig. 1, body morphology, motion pattern and the designated tasks and constraints are the three decisive elements for designing a customized humanoid robot. The origin of this mapping relationship is detailed in the Appendix. However, their relationships are complex and hidden. It is extremely difficult to obtain mathematical expressions. Therefore, an evolutionary process that occurs in the biology and ecological worlds is adopted.

Figure 1. The mapping relationship of body morphology and walking pattern to the designated tasks (living environment).

In the evolution and natural selection processes, the tasks under certain circumstances and constraints serve as the driving force. To create such evolving environments, standard experiments are designed to represent the designated job and existing constraints. An evaluation method is also proposed to score every experiment. In addition to those external factors as the given conditions, many other internal factors, such as the walking control scheme and the optional joint driver for the robot, are also involved. These factors are also considered surroundings because of the limitations of available technologies. For example, a small-sized toy-like humanoid robot tends to adopt an open-loop walking scheme to avoid expensive sensor feedback, and the available joint servo motors are also limited to several known types. This walking scheme is vulnerable to uneven terrain and therefore exerts constraints on walking stability. The available driving motors have limits on their output torque, power and speed. These characteristics are also constraints and are treated as its surroundings and environments.

Body structure is the second key factor for coordinated morphology-locomotion optimization. Many parameters can distinguish one robot from another. Simply, a tiny modification of the bearing for limb joints can be treated as a variation in body structure. However, the concept of body structure here only indicates the morphology configuration rather than the structural details. It mainly includes limb length and the relative ratios of them, such as the ratio of the thigh and shank, the length of legs to its body and the relative pelvis width. The evolution ultimately creates a humanoid robot with the optimal size for these limbs.

The last included factor is the bipedal locomotion and walking pattern. Currently, there are many choices for a biped robot to walk. Some of them are even human-like: inverted pendulum mechanisms and stiff walking are used to achieve energy efficiency. Others are based on the 3DLIPM model [Reference Kajita, Hirukawa, Harada and Yokoi2], and the walking pattern is similar to the chimpanzee’s bending knee, bending hip (BKBH) walking style. The theory behind each walking style is different.

The parameters of walking styles are indispensable for a successful implementation, which is strongly related to the specific walking pattern. Consequently, the parameters selected for evolution also greatly vary. The evolution in this paper primarily focuses on the adjustment of parameters for a pre-defined walking scheme, such as the time duration of a single step and the ratio of the single leg support time. The walking scheme must be described by the same series of parameters among the experiments. In contrast, walking patterns with completely different basic concepts are unable to be compared and selected.

The process for finding the best configuration for locomotion and morphology is essentially the solution of a global optimization problem based on a type of mapping from the robot body configuration and walking pattern parameters to the walking experiments that cover the primary designated tasks. The problem is normally expressed in (1).

(1) \begin{equation}f\;:\;\left\{ {{{\boldsymbol{{P}}}_m},{{\boldsymbol{{P}}}_o}} \right\} \to e \,, \\\end{equation}

where ${\boldsymbol{P}}_m$ is the set of machinery parameters, ${{\boldsymbol{P}}_o}$ is the set of motion parameters and e is the evaluation after the survival test. The constraints are expressed implicitly in the experiments as stable walking without falling over or self-collision.

There is no derivative information in the target function. GA/EA and surrogate-model-based approaches are candidates for this type of problem. However, the surrogate model approach or more precisely the DACE method is preferable because of its far fewer demands in terms of experimental evaluation. These experiments and evaluations are also assumed to be performed by computer simulations to obtain definitive results with the same parameters. Refer to Section 2 for more information and comparisons of the GA/EA and surrogate-based methods. DACE runs iteratively as illustrated in Fig. 2.

Figure 2. Framework of the black-box optimization for morphology-locomotion coevolution.

The process consists of two main steps: the experimental step and the parameter-finding step. In the first step, a series of experiments are performed for survival tests, through which an evaluation value is generated. One or a number of species with different parameters is suggested by means of a certain heuristics function based on this value together with the last trials. Then, new simulations or experiments are performed using the new candidate constants. The entire process proceeds recursively or concurrently for the ultimate determination of the optimal values.

4. An Example of morphology-locomotion coevolution

The previous sections develop a framework and the necessary components for the morphology-locomotion coevolution, while this part presents an example to demonstrate how to obtain a customized design of a humanoid robot with fast and stable walking.

Suppose that a humanoid robot is supposed to be capable of walking rapidly for a race or for football games, such as Robocup or FIRA scenarios. The walking scheme is the model predictive control (MPC) of 3DLIPM [Reference Kajita, Hirukawa, Harada and Yokoi2]. It is chosen because the method has been developed for decades, and it is probably the only applicable, robust omnidirectional walking scheme. Other walking methods, such as the passive walking pattern, have not adopted because of simultaneous limitations on walking forward, sideways and turn. The 3DLIPM walking is also the chimpanzee’s BHBK style, which is much more stable due to its less peak reaction forces. This will bring less vibration, although less energy efficiency is considered. This is not the decisive consideration for a racing robot.

The given conditions are listed for later design.

• The robot height is supposed to be in a specified range, say $(H_l, H_u)$ ;

• The robot joint actuator is already defined, with a known weight, size, maximum speed and power capacity;

• According to the selected driver system, the onboard computing and control platform is already chosen, for example, the servo motor driver, the onboard battery and embedded computers.

• The upper body of a robot is then well designed. Its upper body is already designed with those control and driving systems encapsulated within its trunk. Consequently, its weight, centre of gravity and inertial matrix of the upper body are already well defined.

In the above list, the first line is to bound the body height of a robot, which is the key parameter for the upcoming robot. The height of the robot not only defines its application purpose but also decides the potential joint actuators and the volume of each body part. The joint actuators for a robot with heights lower than 60 cm will be greatly different than those with body heights larger than 170 cm.

Second, the joint actuators supply driving forces for walking locomotion. The joint actuation system, such as servo motors or hydraulic driving systems, is the key factor for building a physical robot. The second assumption is due to the special requirements of the torque/mass ratio for the current bipedal walking technologies. The optional actuators are strictly limited and are presumed to be already definite.

The trunk contains the servo drivers, controllers and onboard computers, as well as the power source and batteries. Sensors such as inertial measurement unit are also mounted within the trunk. The trunk can be well designed before determining the leg morphology. The arm extension can also affect walking ability; however, in this specific problem, we ignore its effect but focus on the vital role of legs. The limb length of arms can be added to form another new design problem, and it makes no essential difference for this sample design problem.

The listed items provide the internal condition discussed in Section 3 and avoid designing from scratch. Note that the problem is only a simplified one to illustrate the application of the coevolution framework. The method can be used in a much wider range in addition to this small-scale problem, but correspondingly at the cost of introducing more parameters and computing power.

As in the proposed sample design problem, three key components are indispensable for formalizing an evolution problem: the mechanical structure parameterization, the walking generation parameterization and the evaluation method in experiments. Therefore, these three factors are discussed. The surrogate-model-based optimization, which acts as the special evolution process, is also briefly presented.

4.1. Body morphology parameterization

Several variable morphology parameters are depicted in Fig. 3:

• Leg length L: how long will the robot legs be design constrained by the body height within $(H_l, H_u)$ ?

• Leg ratio r: what will be the ratio of length between the upper and lower leg?

• Pelvis width p: what is the offset between the two legs?

For the first one in the design task list, the robot tends to be less stable when walking when the robot height reaches $H_u$ , but it can walk with longer strides with longer legs. Although the same walking speed can be reached by changing the step frequency, other elements, such as the limit for power and torque exertion and rolling speed of the joints, may substantially vary from these leg length-step frequency companions. There appears to be a balance point for the design task.

Figure 3. Morphology parameters for the legs of a racing humanoid robot.

The ratio of thigh and shank length varies to some extent in human beings from 0.8 to 1.0, and these morphology variants affect the motion ability in human beings. The condition is similar for pelvis width p. If the pelvis width p is too large, then the step frequency must be lowered for a longer time duration for ZMP swapping between the feet. On the other hand, the robot tends to walk with a less stable margin for a lower pelvis width p.

Note that there are also some other parameters that affect body morphology, such as the height of the ankle, indicated as the dotted line in Fig. 3. They are ignored in this sample problem to have a limited set of involved parameters for morphology-locomotion coevolution, but they may be added to the coordinated optimization for other applications.

4.3. Experimental and scoring method for walking evaluation

The experiences of human beings have already told us that the environment and living pressure drive evolution, and the standard experiments for scoring different parameter trials must contain complete elements for the robot’s survival test. For a football player, the best experiment is in field play of the game, or at least some typical motion for playing football. However, in this sample problem, it is simplified as the ability to walk fast and stably.

A standard experiment is designed for this customized walking robot. The gradual accelerating walk procedure is designed to obtain the experimental scores, as in Fig. 4. For each experiment, the velocity-time curve is the same. Before $t_a$ , the robot is supposed to accelerate uniformly to maximum velocity, and then keep uniform motion. The experimental score is defined as the maximum time the robot can achieve under stable walking. On the basis of constant frequency, the robot can increase the moving speed by increasing the stride. The experiment serves as fair evaluation for different step frequencies and strides. A series of factors are considered in this standard experiment.

Figure 4. Standard experiments as the evaluation method for morphology-locomotion optimization.

Because bipedal walking is a type of discrete motion, it is first simulated as a continuously running vehicle and moving with a constant acceleration. The acceleration process is discretized as independent steps by calculating the translational distance during the time extension of one step. The time duration of one step can be obtained by the step frequency f discussed in Section 4.2. This method makes the comparison possible between the higher step frequency with small strides and lower step frequency with large strides.

Second, such acceleration is thresholded by the largest stride that the robot can reach, which is limited not only by the walking stability but also by the valid inverse kinematics solution of the leg.

Third, the key to obtain the maximum time under stable walking is to judge the last stable walking step. However, at times, it is difficult to distinguish between the stable steps and unstable steps. A robot is starting to fall over in step n, for example, with unexpected collision with the ground, but it can still continue walking for several steps until step $n+4$ before observable body inclination. One solution is to adopt contact forces of the feet to indicate these unstable risks and to predict falling over before observable inclination. Another is the direct observation of geometric status, whereas some delay exists for the latter method compared with the first one. A combination of these two methods is applied to obtain the maximum walking speed. The direct observation of physical inclination of the robot body indicates falling over. The maximum walking step is traced back to the point of contacting forceabnormality.

4.4. Evolution of the robot via surrogate-model-based optimization

A surrogate-model-based approach for optimization is used to build a model that describes the behaviour of the related parameters to avoid the mathematical derivation of these complex functions. Many of these relationships are implicit and hidden or even interdisciplinary. The surrogate-model-based approach can be generalized after the parameterization of related factors to the problem. For example, in ref. [Reference Hemker, Stelzer and Stryk27], a set of walking pattern parameters is selected without a detailed dynamics analysis to optimize for a rapid walking capability, and in ref. [Reference Zhang, Xia, Liu and Chen28], many simplified model parameters are selected to obtain the best omnidirectional reactive walking. The details of the algorithm and its derivation are not discussed here; refer to refs. [Reference Jones, Schonlau and Welch23], [Reference Roustant, Ginsbourger and Deville29] for a complete description. The algorithm for morphology-locomotion evolution is presented in List 1.

In Algorithm 1, $\boldsymbol{x}$ is the parameter set $(L,p,r,f,t_r,h,s)^{\text{T}}$ by accumulating the body structure part in Section 4.1 and locomotion part in Section 4.2. Line 1 presents an original experimental design in which the Latin hypercube approach [Reference Bates, Sienz and Toropov30] is adopted for space filling. An initial Kriging model is created for the first time and is iteratively optimized in Line 4. The fundamental function of the Kriging model is as in Eq. (2).

(2) \begin{equation}\tau({\boldsymbol{x}}) = {{\boldsymbol{f}}}({\boldsymbol{x}})^{\text{T}}{\boldsymbol{\beta }} + z({\boldsymbol{x}})\,, \\\end{equation}

where $\;{{\boldsymbol{f}}}({\boldsymbol{x}})$ is the selected base function, which can be set to constant 1; $z({\boldsymbol{x}})$ is a variable corresponding to a stochastic process, which follows a normal distribution with a mean of $\textbf{0}$ and variance $\sigma$ ; and stochastic variables $z({\boldsymbol{x}})$ in different positions $\boldsymbol{x}$ share the same distribution. The covariance and the correlationship between any two points $\boldsymbol{x}^{(i)}$ and $\boldsymbol{x}^{(j)}$ are denoted in Eq. (3).

(3) \begin{equation}\begin{split}&\operatorname{cov} \left[{z\left({{\boldsymbol{{x}}}^{(i)}}\right),z\left({{{\boldsymbol{{x}}}^{(j)}}}\right)}\right] = {\sigma ^2}R\left( {{\boldsymbol{\theta }},{\boldsymbol{{p}}},{{\boldsymbol{{x}}}^{(i)}},{{\boldsymbol{{x}}}^{(j)}}} \right) \\[4pt] &{r_i}\left( {{{\boldsymbol{{x}}}^{(j)}}} \right) = \operatorname{corr}\left[ {z\left( {{{\boldsymbol{{x}}}^{(i)}}} \right),z\left( {{{\boldsymbol{{x}}}^{(j)}}} \right)} \right] = R\left( {{\boldsymbol{\theta }},{\boldsymbol{{p}}},{{\boldsymbol{{x}}}^{(i)}},{{\boldsymbol{{x}}}^{(j)}}} \right)\,,\\\end{split}\end{equation}

where $R({\boldsymbol \theta}, {\boldsymbol p}, {{\boldsymbol{x}}_1}, {{\boldsymbol{x}}_2})$ are the correlation type and are defined as Eq. (4).

(4) \begin{equation}R({\boldsymbol \theta}, {\boldsymbol p}, {{\boldsymbol{x}}^{(1)}}, {{\boldsymbol{x}}^{(2)}})=\exp \left( { - \sum\limits_{h = 1}^k {{\theta _h}{{\left| {x_h^{(1)} - x_h^{(2)}} \right|}^{{p_h}}}} } \right)\end{equation}

where $\theta_h$ denotes the correlation strength of each $x_h$ component and $p_h$ indicates the smoothness of the correlation function. Equations (2), (3) and (4) present the well-known Kriging model. This model describes the behaviour of the function mainly by the second term: the stochastic variable $z({\boldsymbol{x}})$ .

The parameters of the Kriging model $\left\{ {{\beta}, {\sigma^2}, {\boldsymbol \theta}, {\boldsymbol p} }\right\}$ are found by maximum likelihood estimation with a group of $\tau \sim \boldsymbol{x}$ evaluation tests. Prediction of function and the mean square error (MSE) of the estimation on any new point ${{\boldsymbol{x}}^ * }$ is written as:

(5) \begin{equation}\begin{split}\hat \tau({{\boldsymbol{{x}}}^ * }) &= \hat \beta + {{\boldsymbol{{r}}}^{\text{T}}}{{\boldsymbol{{R}}}^{ - 1}}({\boldsymbol{{y}}} - {\textbf{1}} \cdot \hat \beta )\,. \\{s^2}({{\boldsymbol{{x}}}^ * }) &= {\sigma ^2}\left[ {1 - {{\boldsymbol{{r}}}^{\text{T}}}{{\boldsymbol{{R}}}^{ - 1}}{\boldsymbol{{r}}} + \frac{{{{\left( {1 - {{\textbf{1}}^{\text{T}}}{{\boldsymbol{{R}}}^{ - 1}}{\boldsymbol{{r}}}} \right)}^2}}}{{{{\textbf{1}}^{\text{T}}}{{\boldsymbol{{R}}}^{ - 1}}{\textbf{1}}}}} \right]\,, \\\end{split}\end{equation}

where $\hat \beta$ is the coefficient of the estimated base constants, $\boldsymbol{R}$ is the correlation function matrix of tested points with ${r_i}\left( {{{\boldsymbol{x}}^{(j)}}} \right)$ in Eq. (3) as the $i^{th}$ row and $j^{th}$ column element and $\boldsymbol{\tau}$ is the resulting vector of experimented points. $\boldsymbol{r}$ is the correlation function vector of $\boldsymbol{x}^*$ and the tested points, where the element i of $\boldsymbol r$ is $r_i({\boldsymbol x}^*) = \operatorname{corr}\left( {\boldsymbol x}^*, {{\boldsymbol x}^{(i)}} \right)$ .

To balance global exploration and local exploitation, the expected improvement (EI) criterion is introduced in Line 5 in List 1 as follows:

(6) \begin{equation}\operatorname{E} \left( {I\left( {\mathbf{x}} \right)} \right) = \left( {{f_{\min }} - \hat \tau} \right)\Phi \left( {\frac{{{f_{\min }} - \hat \tau}}{s}} \right) + s\phi \left( {\frac{{{f_{\min }} - \hat \tau}}{s}} \right)\end{equation}

where $\phi(x)$ and $\Phi(x)$ are the probability density function and the probability distribution function of the Gaussian distribution, respectively. s is the MSE of the estimation in Eq. (5). The optimum of Eq. (6) is a new morphology-locomotion parameter couple for the survival test in Line 6. Line 9 presents the best result up to now, and after the designated number of iterations is performed, the final body structure and the corresponding motion parameters are obtained simultaneously as $\hat {\boldsymbol{x}}^{\star\star}$ and output as the current result.

5. Experiment and discussion

An experimental platform is built as illustrated in Fig. 5 to justify the strategy to design a customized humanoid robot. All packages and simulation experiments were performed on a workstation with an Intel Xeon E5-2630 CPU clocked at 2.4 GHz and 16.0 GB memory with Ubuntu 12.04.

Figure 5. Framework for experiment with open source platforms.

The framework is built based on several open source platforms, including R for statistics, OpenHRP for dynamics simulation and Openrobots for motion generation. Two R packages, DiceKriging and DiceOptim [Reference Roustant, Ginsbourger and Deville29], are adopted for Kriging modelling and efficient global optimization (EGO). DiceKriging performs estimation, simulation, prediction and validation of a Kriging model based on a series of experimental points, while DiceOptim performs sequential Kriging-based optimization based on EI criteria as in List 1. Every experiment corresponds to a walking simulation with a robot prototype and motion couple. Each point ever evaluated is obtained by the simulation process and recorded for later analysis along the optimization history. Latin hypercube initialization of the Kriging model is also realized in the R packages.

OpenHRP is used as the dynamics simulator and acts as the evolutionary experimental environment. It calculates the dynamics during the accelerated walking experiment, and it exports the external acting forces and posture of the robot. The sample robot available for the default installation of OpenHRP is used for the prototype (as in Fig. 6). Several modifications were made for this robot.

Figure 6. Sample robot in OpenHRP default installation used in morphology-locomotion evolution.

First, restrictions of power output on each joint are implemented besides the original Proportional Derivative (PD) controller used on the joint actuators because the unlimited power output on the sample robot is irrational on a physical robot. The restriction function is in Eq. (7).

(7) \begin{equation}T \leq \min ( {P_{\max}/\omega , \alpha \cdot i_{\max}} )\end{equation}

whereT is the torque to be exerted on the joint, $P_{\max}$ is the maximum power of the actuator, $\omega$ is the current turning speed and $\alpha \cdot i_{\max}$ denotes the maximum torque by current limitation in the motor actuator.

We also added a passive joint to each foot of the robot to extract the reaction force and momentum from the ground. This passive joint performs equivalently to a 6-D force sensing system on a physical robot. In this way, a state extractor post-processes the raw data provided by the virtual sensors and outputs meaningful results, such as force and posture. These results are evaluated to determine whether the walking motion is stable based on a previous falling-over determinator, and they are converted to the scores described in Section 4.3.

Key parameters of mechanical structure, such as the lengths of the upper leg and the lower leg, are changed by modifying the VRML file for the designated robot model. However, the mass and the distribution of the limb do not vary with its geometric size for simplifying the execution. Because the mass does not vary considerably when the length is different, the variation plays a much smaller role in gait planning when it is integrated into the whole body mass.

The motion planner is realized based on an open source library based on OpenRobots [Reference Mallet31]. It is able to calculate the body joint trajectories based on 3DLIPM. The new robot dynamics parameters are extracted with the same VRML file for OpenHRP simulation, and this guarantees the correspondence of motion generation and dynamics simulation. Additionally, the time cost for a single step and the ratio between the single support time and double support time of the flying foot height extreme are fed into the motion planner as the hyperparameters. The generated sequence of footsteps via Fig. 4 for the specific experiment is submitted to the joint trajectory generator to obtain the COM trajectory via 3DLIPM-based preview control. The joint trajectories of all joints, after inverse kinematics, are generated offline and are followed by those actuators of the robot. A complete dynamic simulation experiment is shown in Fig. 7 through snapshots.

Figure 7. Snapshots of a dynamic simulation experiment.

Thevariables of body configuration and walking pattern together with their limits are then as listed in Table I.

Table I. Body morphology-locomotion parameters involved in the coordinated evolution to build a humanoid robot.

5.1. Results and discussion

For the experiment, 70 initial points, that is, 10 times the vector dimension in Table I, are generated via the Latin hypercube method. These 70 initial points build up the first surrogate model. The leave-one-out experiments are depicted in Fig. 8.

Figure 8. Leave-one-out cross-validation for the previous Kriging metamodel.

Figure 9. Optimization process of 70 Latin hypercube experiments and 300 EGO experiments.

The sequential optimization process displaying the results of the initial 70 experiments and the subsequent 300 iterative EGO experiments is shown in Fig. 9. The lowest evaluation value appears after the 22nd experiment. This indicates that the convergence is realized much faster than expected. The optimized configuration with the lowest evaluation value is listed in Table II for comparison.

The leg configuration in the first line, which is the optimized point that we obtained, indicated that the better choice of the leg design is not the equivalent lengths of the upper and lower legs. The upper leg is found to be 20% longer than the lower leg, and this leg configuration is drawn in Fig. 10. This result surprisingly corresponds to human beings, whose lower leg length is only 80% of that of the upper leg. We think it might be a coincidence, and this may be caused by the evaluation criteria we designed. If we change the optimization objective, the point with the lowest evaluation value may change. And this paper only provides one evaluation criterion described in Section 4.3. Users can redesign the evaluation criteria according to their own needs. A simple explanation is discovered for this configuration. Note that the walking style that we chose is based on 3DLIPM, which is characterized by its knee bending every moment during walking. The leg can never extend straight during this type of walking. The initial bending angle of the knee, which is also called the waist height for BHBK style before walking in Table I, defines the maximum step length, as illustrated in Fig. 11.

Table II. The optimized body morphology-locomotion parameters with the lowest evaluation value.

Figure 10. The robot with the optimized leg configuration.

The rationality of this configuration may be explained from another aspect. The force analysis of the robot leg is shown in Fig. 12. In the state on the left, the torque of the knee joint ( $M_1$ ) is larger, while the torque of the ankle joint ( $M_2$ ) is smaller. In the state on the middle, the ankle torque ( $M_2$ ) is large, while the knee torque ( $M_1$ ) is almost zero. To decrease the torque requirement for stable walking, the acting force line of the weight of the upper body together with the upper leg is supposed to be as close as possible to the ankle and knee joints. But the distance between these two joints can never be eliminated only if the robot is standing upright. However, this will conflict with the bending knee assumption. In the state on the right, in static balance, $M_1\approx G\cdot a$ and $M_2\approx G\cdot b$ . Therefore, a shorter lower leg can reduce the torques of the knee and ankle joints.

Figure 11. Waist height for BHBK style.

Figure 12. A simple force analysis of the robot leg.

With these considerations, when the lower leg is shortened, both the bending angle of the knees and the balance of the ankle joint and knee joint are achieved. However, the length of the lower leg cannot be reduced too much or the limbs tend to collide with ground and cause the robot to fall over.

The above explanation is only a heuristic proposition because the walking procedure itself is complex. It is impossible to predict that a short lower leg is preferred only for static analysis. However, the robot body successfully evolved to fit a better mobility through the approach of coordinated optimization of mechanical structure and motion generation.

Although some extent of evolution of the robot is achieved, the process runs only to a quite limited extent. The different leg ratio of the robot from human beings is mostly due to the special walking style. It is severely constrained by the 3DLIPM walking style, while human beings hardly walk in this way, and therefore, their legs are not necessarily of different lengths (the same length is better suited for the walking style of human beings). However, note that evolution or optimization is oriented to the toy problem as an example. More complex problems can be presented and solved in the framework of Kriging surrogate-model-based evolution. However, it may correspondingly introduce more concerned parameters.

Another factor worth noting is that DACE and EGO iterative optimization are performed for the computer simulations until now. Standard DACE assumes that the experiment is deterministic, which means that the same parameter will output the same experimental result, and that is the reason why the method can only be utilized in computer simulations rather than on a real robot, and this fact will undoubtedly introduce a gap between the simulation result and real robot execution. The method is used primarily in the design period, and such optimization on simulations may vary from a physical robot. This is because of the lack of determinacy in a real robot experiment and the fact that the performance of such a robot tends to be volatile even given the same model constants. This problem can be explained in two aspects: First, there are already several solutions available for noisy systems of this type, and these methods will be tested in future research. Second, a special variable leg configuration can be implemented on a physical robot to enable online optimization of mechanical structure and walking style. During the evolution process, the experimental condition is carefully controlled to obtain a relatively stable outcome for the same point.