3.1. Formulation and variables

Numerical trajectory optimization is a powerful and versatile tool for robotics. Optimizing a tailored cost function allows to generate a control trajectory for the robot so that it performs specific behaviors [Reference Mordatch, Todorov and Popovic47–Reference Todorov50]. The main advantage of this approach is the intuitiveness of setting the cost and constraints, which are strictly related to high-level goals that must be fulfilled. For the trajectory optimization problem on a discretized horizon with nodes $[0..N]$ , we use direct collocation with an augmented set of variables: $\mathbf{Z}=\left [\mathbf{X}, \mathbf{U}, \mathbf{A}, \mathbf{F}, \boldsymbol{\Gamma }\right ]$ where:

• • $\mathbf{X}$ is the decision vector collecting the evaluations of states of the robot $\mathbf{x}$ , each state includes the configuration and velocity of all its degrees of freedom. where $\mathbf{q}_{\text{b}}$ is the underactuated base position and $\mathbf{q}_{\text{a}}$ is the vector of actuated joint positions.

  (1) \begin{align} \mathbf{x} = [\underbrace{x, z, \theta, q_{1..n_u}}_{\mathbf{q}=[\mathbf{q}_{\text{b}}, \mathbf{q}_{\text{a}}] \in \mathbb{R}^{n_q} } \underbrace{ \dot{x}, \dot{z}, \dot{\theta }, \dot{q}_{1..n_u}}_{\mathbf{v} = [\dot{\mathbf{q}}_{\text{b}}, \mathbf{v}_{\text{a}}] \in \mathbb{R}^{n_v}} ] \in \mathbb{R}^{(n_x=n_q+n_v)} \end{align}

• • $\mathbf{U}$ contains the actuated joint torques $\mathbf{u} \in \mathbb{R}^{n_u}$ .

• • $\mathbf{A}$ is the vector of the joint accelerations $\mathbf{a} =\dot{\mathbf{v}} \in \mathbb{R}^{n_v}$ .

Just for the contact phase nodes $C=[N_{c,0}..N_{c,T}]\subseteq [0..N]$ , the foot position $\mathbf{p}_{c, f}$ for the feet $f$ in contact is fixed. For any foot $f$ in contact, we define additionally:

• • $\mathbf{F}$ : contact force vector, which stacks the contact forces on the contact $f$ for all the nodes in which a given contact phase $C$ is defined, $\mathbf{f} = \bigcup _{C}\bigcup _f (\mathbf{f}_{c, f} \in \mathbb{R}^{n_c})$ , where $n_c$ is the contact point dimension, which depends on the contact model. For instance, in planar models $n_c=2$ , while for three-dimensional and contact wrench models it equals $n_c=3$ and $n_c=6$ , respectively.

• • $\boldsymbol{\Gamma }$ : slack variable vector, which collects the contact slack variables acting on the contact $f$ for all the nodes in which a given contact phase $C$ is defined, $\boldsymbol{\gamma } = \bigcup _{C}\bigcup _{f} (\boldsymbol{\gamma }_f \in \mathbb{R}^{n_c})$ , following the formulation in ref. [Reference Posa, Kuindersma and Tedrake48], to impose constraints on the feet position and velocity. These variables are introduced only in the contact phases to avoid contact drift.

The choice of considering such an augmented set of variables, at torque and acceleration level, is motivated by the intention of imposing constraints which reflect the physical limitations of the actuator. This is not directly possible in the case of simplified models such as linear inverted pendulum [Reference Kajita, Kanehiro, Kaneko, Fujiwara, Harada, Yokoi and Hirukawa51, Reference Kajita, Nagasaki, Kaneko and Hirukawa52], spring-loaded inverted pendulum [Reference Kajita, Morisawa, Miura, Nakaoka, Harada, Kaneko, Kanehiro and Yokoi53], or centroidal model [Reference Dai, Valenzuela and Tedrake54, Reference Orin, Goswami and Lee55], which do not include the joint torques. Finally, highly dynamic behaviors are difficult to discover as they are often far from the simplified model assumptions.

3.2. Optimal control problem constraints

With the formulation outlined in Section 3.1, constraints can be imposed directly on the primal variables both in the form of equality and inequality constraints. This is an aspect of the utmost importance for co-design, as the feasibility of the motion needs to be guaranteed from the optimization stage.

Robot dynamics.

The robot state $\mathbf{x}$ evolves under the influence of the joint torques and contact forces as described by the constrained rigid body dynamics [Reference Featherstone56]:

(2) \begin{align} \begin{bmatrix} \mathbf{M} & \mathbf{J}_c^\top \\ \mathbf{J}_c & 0 \end{bmatrix} \begin{bmatrix} \mathbf{a}\\ -\mathbf{f} \end{bmatrix} = \begin{bmatrix} \mathbf{u} - \mathbf{b}\\ -\dot{\mathbf{J}}_c \mathbf{v} \end{bmatrix}, \end{align}

where $\mathbf{M}$ is the joint space inertia matrix, $\mathbf{b}$ is the vector containing the state-dependent nonlinear effects of gravity, centrifugal and Coriolis forces, and $\mathbf{J}_c$ is the contact Jacobian stacking the Jacobians of all the contact points. Based on this dynamics, the robot configuration $\mathbf{q}$ and its velocity $\mathbf{v}$ evolve under the control action $\mathbf{u}$ of the motors. Equation (2) can be solved using the forward dynamics [Reference Featherstone56], leading to a constraint on $\mathbf{a}$ . Joint accelerations must then be integrated numerically to obtain joint velocities and positions. To this aim, we introduce an integration function $\boldsymbol{\Phi }$ , which we used to formulate the following constraints

(3) \begin{align} \mathbf{x}^{+}= \boldsymbol{\Phi }( \mathbf{x}, \mathbf{a}, \mathbf{u}, \underbrace{\mathbf{f}, \boldsymbol{\gamma }}_{\text{if contact}}, \Delta t ) \end{align}

To improve numerical conditioning, the contact point velocity is corrected with a slack variable $\boldsymbol{\gamma }$ as proposed in ref. [Reference Posa, Kuindersma and Tedrake48]. This allows redundant constraints on the contact location and its velocity, avoiding drift. This modification is propagated in the integrator law (3), as shown in Eq. (13) of ref. [Reference Posa, Kuindersma and Tedrake48]. In the cost function, these slack variables are penalized for converging to physically accurate solution. Our integrator hence depends also on $\gamma$ because, instead of the state velocity $\mathbf{v}$ , the value $\tilde{\mathbf{v}}$ is used in the integration step, where $\tilde{\mathbf{v}}$ can be interpreted as the velocity projected in the kernel space of the contact velocity $\mathbf{J}_c \mathbf{v}$ :

(4) \begin{align} \tilde{\mathbf{v}} = \mathbf{v} + \mathbf{J}_c^\top \boldsymbol{\gamma } \end{align}

Contact constraints.

The rigid contact model leads to several constraints described in the following.

Forces

The non-sliding and unilaterality conditions impose the following constraints on any contact force (for flat ground) $\mathbf{f}_c = [f_{c,x}, f_{c,y}, f_{c,z}]^{\top }$ , given the friction coefficient $\mu$ :

(5) \begin{align} \begin{cases} \mu ^2{f}_{c,z}^2 \geq{f}_{c, x}^2 +{f}_{c,y}^2\\{f}_{c, z} \geq 0 \end{cases} \end{align}

Non-sliding contact points

During any contact phase of horizon $C=[N_{c,0}..N_{c,T}] \subseteq [0..N]$ , the position $\mathbf{p}_{c, f}$ of any foot $f$ in contact is constant for the whole phase. In particular, we set it equal to the value at the beginning of the phase:

(6) \begin{align} \mathbf{p}_{c, f}(\mathbf{q}_i) = \mathbf{p}_{c,f}(\mathbf{q}_{N_{c,0}}), \quad \forall \, i \in C, i\neq N_{c,0} \end{align}

Non-penetration

The $z$ coordinate of the contact point must be at ground level (flat ground assumption):

(7) \begin{align} \mathbf{p}_{c, f}(\mathbf{q}_{N_{c,0}})|_z = 0 \end{align}

Because of (6), this condition can be imposed just on the initial contact node $N_{c,0}$ .

Contact velocity

The velocity of the feet in contact must be zero:

(8) \begin{align} \mathbf{v}_{c, f}(\mathbf{q}_i) = 0, \quad \forall \, i \in C \end{align}

Collision avoidance with the ground.

To produce a feasible motion, constraints on the vertical position of some key-frames (e.g., shoulder and knee joints, indicated with the subscript $k$ ) need to be imposed in order to not penetrate the ground. This is enforced along the whole optimization horizon through inequalities of the type:

(9) \begin{align} \mathbf{p}_{k}(\mathbf{q}_i)|_{z} \geq 0, \quad \forall \, i \in [0..N] \end{align}

Cyclicity.

As noted in ref. [Reference Hubicki, Jones, Daley and Hurst57], limit cycles naturally arise from trajectory optimization in the simple cases of optimal locomotion. Cyclic motion patterns are then our chosen optimization target to obtain energy-efficient hardware. This choice also makes it possible to keep the optimization horizon per cycle short and therefore reduces computational complexity without sacrificing numerical precision. Once the motion primitive is obtained, a locomotion pattern that is representative of the robot’s main operation can be achieved by replicating the cycle multiple times. The periodicity of the solution is introduced in the OCP with non-Markovian constraints between the optimization variables at the initial and final nodes of the problem. Depending on the problem requirements, these constraints can involve the full set of decision variables $\mathbf{z}$ or just a subset of it.

(10) \begin{align} \mathbf{g}(\mathbf{z}_0, \mathbf{z}_N) = 0 \end{align}

For instance, some offsets or inequalities can be introduced just on specific parts of the state to enforce a given behavior (e.g., in a forward jump, we want that the base position translates at least of a given amount, but all the other variables match the values at the beginning of the trajectory).

(11) \begin{align} \mathbf{h}(\mathbf{z}_0, \mathbf{z}_N) \leq 0 \end{align}

These constraints enforce a dependency between the initial and final nodes. The major drawback is that the requirements to define a Markov chain are not respected anymore. This breaks the NLP sparsity pattern and makes faster iterative algorithms as DDP [Reference Mayne41] not viable.

Actuator model and limits.

The choice of the motor is modeled so to change cost function (power consumption) and limits in a physical consistent way. All the main actuator limits are taken into account as a set of constraints:

• • Position: joint position bounds are considered

  (12) \begin{align} \underline{\mathbf{q}} \leq \mathbf{q}_a \leq \overline{\mathbf{q}} \end{align}
  where $\mathbf{q}_a$ are the actuated joints angular positions and $\underline{\mathbf{q}}, \overline{\mathbf{q}}$ are their minimum and maximum allowable values (e.g., the knee joint angle is delimited by the presence of stoppers).

• • Velocity: each actuator speed limit is considered by imposing bounds on the joint angular velocity.

  (13) \begin{align} \underline{\mathbf{v}} \leq \mathbf{v}_a \leq \overline{\mathbf{v}} \end{align}
  where $\mathbf{v}_a$ are the actuated joints angular velocities and $\underline{\mathbf{v}}, \overline{\mathbf{v}}$ are their and maximum allowable values. For highly dynamic trajectories, this aspect is essential as these thresholds may easily be reached.

• • Torque: Generally, torque limits are modeled as fixed bounds on $\mathbf{u}$ . This is a necessary but not sufficient condition because the actuator cannot instantly provide arbitrary torque values: the intrinsic limitation due to the bandwidth of the actuation needs to be addressed. Approaches to treating it were proposed in refs. [Reference Grandia, Farshidian, Ranftl and Hutter58, Reference Gupta59] working in the frequency domain respectively on the cost function and to obtain feedback gains that can be applied to the real system. Our approach is to impose physically driven bounds on the torque values themselves. The rationale is that, considering the joint transmission, the elastic elements (particularly the transmission belt) can store energy through small deformations. This acts as a low-pass filter from the motor to the connected joint, which can be approximated by a first-order filter whose cutoff frequency depends on the actuator technology. In time domain, the filter presents a straightforward implementation. For each node $k \in [1..N]$ it results in the following constraints:

  (14) \begin{align} \begin{cases} \mathbf{u}_k \leq (1-\alpha ) \mathbf{u}_{k-1} + \alpha \overline{\mathbf{u}}\\ \mathbf{u}_k \geq (1-\alpha ) \mathbf{u}_{k-1} + \alpha \underline{\mathbf{u}}\\ \underline{\mathbf{u}} \leq \mathbf{u}_{0} \leq \overline{\mathbf{u}} \end{cases}, \end{align}
  where $\alpha \in [0,1]$ depends on $f_c$ and the discretization step $\Delta t$ as follows:
  (15) \begin{align} \alpha = \frac{2 \pi \Delta t f_c}{2 \pi \Delta t f_c + 1}. \end{align}
  $\underline{\mathbf{u}}, \overline{\mathbf{u}}$ are respectively the minimum and maximum torque that can be achieved by the actuator. By construction, (14) respects peak limits, as $\underline{\mathbf{u}} \leq \mathbf{u}_k \leq \overline{\mathbf{u}}$ $\forall k \in [1..N]$ .

3.3. Power cost function

For what concerns the cost function that is minimized, in the Lagrange term, the total electrical energy consumption is included as the time integral $\int _0^{T} P_{el}(t) dt$ of electrical power $P_{el}$ , as in refs. [Reference Fadini, Flayols, Prete, Mansard and Souères33, Reference Fadini, Flayols, Prete and Souères34]. $P_{el}$ is computed with the non-ideal dissipations of the actuators.

Joint friction

The power dissipation due to friction is computed from the identified values of static friction $\boldsymbol{\tau }_{\mu }$ and viscous friction $b$ .

(16) \begin{align} P_{f} = \mathbf{v}_a^\top \underbrace{(\tau _{\mu } \mbox{sign}(\mathbf{v}_a) + b \mathbf{v}_a)}_{\boldsymbol{\tau }_f}, \end{align}

where $\mathbf{v}_a$ is the velocity of the actuated joints. We denote as

(17) \begin{align} \boldsymbol{\tau } = \mathbf{u} + \boldsymbol{\tau }_f \end{align}

the total joint torque including the friction component $\boldsymbol{\tau }_f$ .

Joule effect

The Joule power losses are included on the motor side with the values of the motor constant $K$ coming from its specifications:

(18) \begin{align} P_{t} = \boldsymbol{\tau }^\top \mathbf{K} \boldsymbol{\tau }, \end{align}

where $\mathbf{K}$ is a diagonal matrix containing, for each joint, the reciprocal of the motor constant divided by its squared gear ratio.

Mechanical power

The mechanical power flow to/from the system can be obtained by considering the instantaneous power at the level of the actuator:

(19) \begin{align} P_m = \mathbf{u}^\top \mathbf{v}_a \end{align}

Total electrical power

The total electrical power is equal to the losses of the system plus the mechanical power:

(20) \begin{align} P_{e} = P_{f} + P_{t} + P_{m}, \end{align}

Mechanical energy invariance

For any periodic trajectory $\Omega$ , parametrized by the variable $\xi$ , the electrical energy equals the integral of the losses ( $P_f + P_t$ ). Therefore, it is not necessary to minimize explicitly the mechanical power $P_m$ because its circuitation is a conserved quantity and equals the difference in mechanical energy between the final and initial state (which is state-dependent and hence zero by definition of periodicity):

(21) \begin{align} \oint _{\Omega } \mathbf{u}^\top (\xi ) \mathbf{v}_a(\xi ) d\xi = (\mbox{E}_{mec} = \mbox{E}_{kin} + \mbox{E}_{pot})\rvert _{\mathbf{x}_0}^{\mathbf{x}_T = \mathbf{x}_0} \triangleq 0 \end{align}

This result can also be extended for semi-periodic trajectories. In particular, we consider the case in which joint velocities are the same and only the $x$ position of the robot base changes. Any translation of the base along $x$ is tolerated as it results in no net change in potential energy because:

• • the base lands at the same $z$ position it started from

• • the actuated joint position trajectories are cyclic

This is a sufficient condition: the final height of each link CoM is equal to the initial one, so no difference in potential energy is induced, and kinetic energy is conserved as there is no difference in state velocity (and the joint space inertia is invariant to base translations).

3.4. Problem specification example

Before presenting the optimization cases, for clarity, the notions introduced in Section 3 are visually represented and discussed in Fig. 2 with a toy problem involving a planar model. In such an example, the goal is to obtain a jump of at least 0.5 m. The specification of the task starts from the selection of several contact phases and timing. To perform this task, 4 distinct phases are chosen: (i) at first, the robot is in contact with both feet (rear and front, double support), then (ii) initiates the jump with both feet (flying phase), (iii) lands on the front legs (single support), and finally (iv) reaches double contact with the ground again (double support). In Fig. 2, the different phases are shown with different coloring, and for each of them, the node variables are reported together with the constraints and cost functions.

Figure 2. Example of a problem construction for a cyclic task of forward jumping of at least 0.5 m.

It is worth mentioning that some variables, costs, and constraints are phase-dependent. For instance, in the flying phase, the slack variables $\boldsymbol{\gamma }$ and the contact forces $\textbf{f}_c$ are not necessary, so they are not introduced. The same goes also for friction constraints, contact constraints, and slack variable penalization (as $\boldsymbol{\gamma }^\top \boldsymbol{\gamma }$ ). Similarly, also the penalization of the slack variables is only present when the slack variables are defined, that is only in the contact phases. The method consistently provides results with $|\boldsymbol{\gamma }| \lt 10^{-6}$ , when the penalization weight is $10^{5}$ . The only other cost component we introduce is related to energy consumption, which has a relatively low weight of $10^{-3}$ . This cost is present in all phases and is the main target of the optimization. Concerning the constraints, the dynamics constraint, the actuator limits, and the collision avoidance are imposed on all nodes, as can be seen by the row “Constraints.” The jumping task is enforced through an inequality in the $x$ base position among the initial and final nodes of the trajectory, in this example $x_N - x_0 \geq 0.5$ , while all the other optimization variables are constrained to be equal between initial and final nodes. The last constraint is non-Markovian and overarches the whole optimization problem. Starting from this simple example, similar contact and cyclic constraints allow to define also the problems which we present and solve in the results part. The strength of this problem formulation is that it can generate very different movements with just the selection of the final constraints and gait pattern as will be shown in the following results with bounding and backflip. Those will still be made up of the same OCP formulation that we have here defined. The different possible motions which are achievable can be seen also in the companion video.

4.1. The new quadruped prototype at DFKI

The DFKI Robotics Innovation Center recently developed a new robot quadruped (see Fig. 3). The physical dimensions for the initial design are similar to the mjbots-quad [Reference Pieper60]. The validation and optimization results are based on the preliminary design, which is presented in this section. The goal of such research platform is to be capable of producing dynamic motions while keeping a certain degree of autonomy. The robot consists of a central body on which four legs of identical design are mounted. Each leg has 3 degrees of freedom (DoF), which are actuated by off-the-shelf quasi-direct drive actuators based on open-sourced MIT mini Cheetah actuators. The hip joint has one pitch, and one roll DoF and the knee joint can be rotated around a pitch axis. To keep the leg’s inertia low, the knee joint’s motor was shifted to the pitch axis of the hip joint and coupled to the knee joint via a toothed belt transmission with a ratio of $\frac{1}{2}$ . All structural elements of the robot were designed in such a way that they can be manufactured by waterjet cutting (including the gearbox pulleys). The components of the body consist of carbon fiber-reinforced plastic (FRP) plates with a thickness of 1 mm. This allows easier manufacturability and assembly, without sacrificing rigidity and lightness. The connections between the hip drives were made from 3 mm thick carbon fiber plates connected by aluminum parts. Likewise, the leg structures are made of carbon fiber plates connected by spacer bolts in the case of the upper leg and by a custom-designed plastic spacer in the case of the lower leg. Lastly, the feet of the robot are 3D-printed and exchangeable. This allows different material hardnesses to be tested for different substrates.

Figure 3. Quadruped prototype bounding tests at DFKI-RICFootnote 1.

Figure 4. The actuator model allows a close match between the ideal trajectories with friction compensation and the ideal torque applied to the system from measurement data.

4.2. Actuator model and power consumption validation

The trajectory optimization formulation introduced in Section 3 is used to produce an energy-optimal bounding motion (for more details on the task, see Section 5.2). By tracking the optimal reference trajectory with the prototype, the gap between the model and reality is assessed, and the models are validated. Figure 4 shows that the actuator model with the identified parameters, closely predicts the total joint torque $\boldsymbol{\tau }$ (including joint friction $\boldsymbol{\tau }_f$ ) as in (17). A jumping trajectory cycle lasts $0.8$ s, so it is repeated multiple times, with a phase in which the system resets to the initial position. To stabilize the trajectory, a PD joint position controller is used, with additional feedforward torques from the OCP. The value of the joint torques predicted by the model closely follows the measures, with the main difference in the flying phase [0.3 s, 0.5 s], which can be attributed to the unmodeled controller dynamics. In Fig. 5, the power prediction from joint data measurements (torque and velocity) is shown together with the measured data. The estimation of the total electrical power is given by $P_{e} = P_{m} + P_{t} + P_{f}$ , with the notation introduced in Section 3.3. To compute the values, $\mathbf{\tau }$ is inferred by our joint friction model. Figure 5 shows that the prediction, which solely uses joint measurements (velocities and commanded torques), follows the measurements of the electrical power provided by the power source, which is measured as the time average product of voltage $V$ and current $i$ , $\hat{P}_{e} = i V$ (at a lower sampling rate). Nonetheless, the integrated values of the total electrical consumption are accurate, despite the controller dynamics and the sim-to-real gap. These findings are reported for the energy-optimal trajectory in Table III. To give some context to these values, a heuristic baseline is used as a reference. Such baseline is based on a simple motion generation technique that is inspired by the energy-shaping approach [Reference Tedrake61], as implemented in Section II B of ref. [Reference Soni, Harnack, Isermann, Fushimi, Kumar and Kirchner62]. The trajectory was tuned to produce a similar jump on the prototype (in terms of time horizon and base displacement). Simple parabolic trajectories based on accelerations are computed for each leg and the configuration in stance and flight phases are pre-defined. However, it is indeed a hand-tuned technique that is error-prone, time-consuming, and, most of all, sub-optimal (for energy efficiency). Trajectory optimization overcomes these weaknesses, and we compare the margin of improvement over the baseline method available. On this handcrafted trajectory, the knowledge of the actuator was applied to assess the electrical power consumption. It was found that the power consumption of the motion obtained with the heuristic was higher (30% increase) than the energy-optimal trajectories. Even if this result is not conclusive to claim the optimality (we just test the optimization output against a simple baseline), it is still an important step to confirm that our method can accurately estimate and minimize energy expenditure before including it in the co-design optimization.

Figure 5. The electrical power estimation $P_{e}$ (blue) closely follows measured one $\hat{P}_e$ (orange).

5.1. Problem requirements and assumptions

High-level requirements for the platform are (i) to produce stable locomotion in the forward direction $x$ , (ii) to be capable of dynamic motions along the $z$ axis as shown in Fig. 6. In order to consider representative legged robot movements, we focus on the generation of (iii) stable and periodic motion patterns. Such movements need to be performed while being (iv) energy-efficient. Taking this into account, periodic bounding and backflip were selected as benchmark tasks to achieve.

Robot model

Figure 6 shows a sketch of the joint placements on a complete robot. The general design choice is to place the motors as close as possible to the base to limit the reflected inertia of the leg links. Another preliminary design choice is to drive both the abduction joint and the hip joint directly, while using a belt transmission with a reduction factor of $2$ for the knee joint. Since the motion of leg abduction in the lateral plane $(y,z)$ is not needed for bounding or jumping, a planar model was used instead of the complete one. The motors are located on the robot’s trunk. The abduction of the leg (rotation around $x$ of the first leg joint) is blocked (the corresponding motor is located inside the base), while the motors actuating the hip and knee joints are shown in gray in Fig. 6. The robot model for this task exploits the symmetry of the motion with respect to the $(x,z)$ plane. The dynamical equivalence between the 3D model and the planar one is ensured by lumping each link inertia and control effort on a unique joint for each symmetrical hip and knee. The command torque on the joints, and their limits, are then doubled and need to be equally divided into two legs when controlling the real system.

Table III. Energy consumption values for the jump.

Figure 6. Complete robot model (left), its planar simplification (center), and scaling of the base, upper leg, and lower leg links.

Design variables.

For both co-design tasks, we optimize over the same set of variables, which is here reported.

Continuous design variables

Starting from the nominal design, the following continuous design parameters are as follows:

• • Lower leg link scaling $\quad \lambda _l \in [0.5, 1.5]$

• • Upper leg link scaling $\quad \lambda _u \in [0.5, 1.5]$

• • Base scaling $\quad \lambda _b \in [0.5, 1.5]$

Discrete design variables

For the implementation, two promising motor candidates were selected from the off-shelf Antigravity AK series as reported in Table IV. In particular, among AK80-6 and AK80-9, these two motors differ mainly from the reduction of the integrated rotary gear, which is respectively $6$ and $9$ . Negligible differences are found for the other parameters, especially concerning the motor constant and the winding resistance. In the co-optimization problem, the same leg design, and consequently actuator choice, is used for all four legs. The motor combinations for the hip and knee motors (respectively $m_{hip}, m_{knee}$ ) are then four. CMA-ES, which by default works on a continuous space, is modified to work on this discrete domain. This is not restrictive, as the discrete variables could also be handled by another GA selection, with minimal modifications to the framework. In our case, as the combinatorics of the problem are rather small, the discrete variables are internally considered continuous, and a remapping strategy is employed to pass from the continuous to the discrete counterpart before solving the problem. Thanks to this remapping, all of the motor characteristics are found with the catalog value without the need for an explicit parametrization as in refs. [Reference Yesilevskiy, Gan and Remy22, Reference Yesilevskiy, Gan and Remy23, Reference Fadini, Flayols, Prete, Mansard and Souères33, Reference Fadini, Flayols, Prete and Souères34]. This mechanism enables optimization even in cases in which a continuous parametrization is not viable (e.g., motor technologies are rather different from each other).

Table IV. Properties of the motor selection integer variables.

Actuator choice.

The actuator properties are taken into account by modifying the robot dynamics, the constraints of the OCP, and the cost function. The main effects of the actuator are as follows:

Joint limits

The choice of a motor and its transmission directly impacts the effort limit (torque) and speed limit which is achievable by the controlled joint.

Inertia

The added rotor inertia is considered in the model via the technique explained in Ch. 9.6 of ref. [Reference Featherstone56]: the rotor inertia is multiplied by the value of the squared reduction and added to the corresponding diagonal element of the joint space mass matrix.

Transmission friction

Given a motor and its transmission, the overall viscous and Coulomb friction are considered in the cost function that minimizes the overall energy, following the same approach as in ref. [Reference Fadini, Flayols, Prete, Mansard and Souères33].

Motor placement

The motor masses are also taken into account by the structural base scaling. Each motor is modeled as a localized mass, and its contribution adds to the parent link mass and rotational inertia.

Actuator bandwidth

From the current technology of the actuator and testing, the cutoff frequency can be estimated as $f_c \approx 20Hz$ . This is valid for both motors, and such choice constrains the torque trajectory conservatively.

Figure 7. Bounding task: (a) shows the different motion phases. Trajectories for the optimal and the nominal designs are respectively shown in (b) and (c). In both, from left to right, the plots show base, joint positions, and joint torque trajectories. Contact phases are highlighted with gray background.

5.2. Co-design for bounding

The first optimized task is a bounding motion, where the robot must perform a jump of at least $0.30$ m. The cyclic constraints enforce the robot state to be equal at the beginning and the end of the trajectory, except for the base $x$ position. Finally, a constraint is added to obtain zero joint velocities at the start and end of the trajectory. In this way, the system consumes just the energy required to perform the jump and decelerate to a full stop in the final part of the trajectory. The phases of such movement are as follows (Fig. 7(c)):

• • Double support, with all feet in contact with the ground.

• • Flying phase, with no contact with the ground.

• • Double support, with all feet again in contact with the ground.

This task is symmetrical, meaning that the time left for each contact phase is the same. For the overall problem, the time window for each cycle of the jump is $0.7$ s, and the total number of nodes for the optimization horizon is 100.

Figure 8. Convergence of the algorithm along the evolution of the populations.

Outer loop hyper-parameters

For this optimization, the CMA-ES algorithm is initialized to evolve for 10 times a population of 1000 different individuals (different combinations of the design parameters). Depending on its complexity, each OCP problem’s computation time varies between $\approx 10$ s and $10$ min. Thanks to parallelization, the whole optimization routine is carried out in $\approx 10$ h. This time is comparable to the clock time of our previous contributions [Reference Fadini, Flayols, Prete, Mansard and Souères33, Reference Fadini, Flayols, Prete and Souères34], but for a much more complex system and task (more DoFs) and more accurate model of the system and dynamics (which is being used in the OCP). Figure 8 shows that this is sufficient to reach stationary values in the cost. It is clear from the trends that there is a diminishing return in exploring further combinations of parameters. In particular, in the same figure, we see different bands, which correspond to the various optimal designs for the 4 combinations of the motors. For simplicity, in the following the combinations will be called AK80 9-9, AK80 9-6, AK80 6-9, and AK80 6-6, where the first index is for the hip joint and the second for the knee joint.

Discussion

For this task, the optimal design is obtained for the values reported in Table V. We see that, with respect to the nominal leg design, the best solution is found for a smaller robot.

According to Table V, the optimization selects as best suited a smaller robot, with a different scaling of the thigh $\lambda _u$ and shank $\lambda _l$ , in particular $\lambda _u/ \lambda _l = 1.46$ . For jumping forward, it seems then that robots with longer thighs are performing better. The optimal solution is chosen to be very close to the lower bound of the variables $\lambda _b, \lambda _l$ . An additional effect of the choice of the base can be observed in Fig. 7(b): when the base scaling is reduced (Table V), the trajectories of the knee and hip joints are showing a higher degree of similarity. In the nominal case, the joint position is reaching the position limits of the actuator, which is no longer the case with the optimized hardware. Basically, we can explain this result as follows. The optimal quadruped for bounding tends to be shaped as a planar biped: since there is no advantage in carrying additional mass from an energetic point of view, the base length is chosen as short as possible. From the joint positions of the nominal design (Fig. 7(c)), the knee stopper can partially limit the robot’s motion. So, finding a solution that does not impose a limitation would be advisable. For both designs, the optimal joint trajectories are smooth and not hitting the velocity bounds. So, concerning the actuator selection for this task, a motor capable of quick motions is not really necessary. Conversely, the choice of a higher gear ratio allows to exert larger torques and to greatly decrease the Joule consumption. To produce the same output torque (as the motor constant is the same) the ratio of the Joule dissipation of the motor types AK80-9 and AK80-6 are equal to the quotient of the square of their gear ratio, so $2.25$ , even if a higher reduction also impacts the reflected inertia, reducing the actuator transparency and joint maximum speed. For bounding the higher gear ratio is selected.

Table V. Results of the optimization for the bounding task.

Figure 9. Backflip task: (a) shows the different motion phases. Trajectories for the optimal and the nominal designs are respectively shown in (b) and (b). In both, from left to right, the plots show base, joint positions, and joint torque trajectories. Contact phases are highlighted with a gray background.

5.3. Co-design for backflip

As a second task, we present the result of a backflip optimization, as shown in Fig. 9. This motion was selected as a complex and dynamic task, exploiting the dynamics of the whole system. The robot starts with zero velocity and has to perform a full rotation of the base before landing. In the landing phase, the excess velocity needs to be damped to reach a full stop at the end of the trajectory. For this task, all joint positions except the base $x$ component need to be equal at the beginning and the end of the trajectory. For this motion, the total time to perform the task is 1 s. As represented in Fig. 9(a), the different phases are as follows.

• • Double support, with all the feet in contact with the ground. In this phase, the motors need to accelerate the base to produce enough vertical velocity to break the contact with the ground. Moreover, the applied forces need to generate enough momentum for the upcoming rotation of the base.

• • Single support, with the front legs taking off. This phase is added to allow the robot to start the rotation of the base and still push the ground with the back legs.

• • Flying phase, in which there is no foot in contact with the ground and the base is following a ballistic movement. The motion of the legs is not contributing to the jump but is useful to get the feet in the right position before landing (preparing for the impact phase).

• • Single support, with the front legs touching the ground first.

• • Double support, with the rear legs reestablishing contact with the ground. The contact needs to be stable, so the forces are inside the friction cone and the motors bring the robot to a full stop at the end of the trajectory.

Outer loop hyper-parameters

CMA-ES is initialized so that each generation is made up of $10^3$ individuals, and the number of overall evolutions of the population is fixed to $10$ . In this problem, as the task is more challenging, some designs could not physically satisfy the constraints and perform the motion within the problem constraints. IpOpt provides debug information on the infeasibility of the problem. If an individual is unfeasible, an arbitrary high-cost value, higher than the other feasible designs, is assigned to it. The outer loop algorithm is elitist, meaning that when generating a new population it will automatically discard the outlier designs.

Discussion

The results of this optimization are reported in Table VI. Running the optimization routine, we notice that the leg size is reduced while the base dimensions are slightly increased.

Table VI. Results of the optimization for the backflip.

The optimization selects a smaller robot, but interestingly a different optimal scaling of thigh $\lambda _u$ and shank $\lambda _u$ is found (with ratio $\lambda _u/ \lambda _l = 0.73$ ) with respect to the bounding task. For backflips, it seems then that robots with longer shanks are performing better. The optimal solution is very close to the lower bound of the variable $\lambda _u$ . Figure 9(b) and (c) report the optimal and nominal design trajectories. The optimal base scaling is obtained with a bigger base without reaching the upper bound. This can be explained as there is a tradeoff between the base inertia and the capability to apply momentum to perform a full base rotation around $\theta$ of $- 2 \pi$ rad. For the same applied contact forces, the longer the base, the easier the backflip can be performed. However, there is still a tradeoff as a bigger base increases also the inertia of the rigid body. For the backflip, it seems that the most critical constraint is the maximum joint velocity. The nominal design, featuring a higher reduction, can exert more torque but reaches the joint velocity limit. A motor that can produce faster motion is needed for task completion. The optimization leads to a smaller reduction to achieve a higher joint velocity, differently from what was obtained for the bounding task. In this case, as the motion is quicker, the effect of the rotor inertia is higher and a smaller reduction helps in accelerating the joint. Compared to the torque required for bounding, in this case, saturation is reached. In Fig. 9(b) and (c), the low-pass filter effect can be noticed from the smooth torque trajectories that do not exceed the upper and lower torque bounds of the actuator (shown with dotted lines). For this task, the knee joint is reaching the position limits of the actuator, in both designs (nominal and optimized). Such limit needs to be taken carefully considered for the final robot design.

5.4. Landscape analysis for multiple objectives

As rather different designs were produced for the two tasks by the co-design optimization (Tables V and VI, Fig. 10), an additional grid search was performed to better understand the impact of the design selection. In this case, the base was kept to the nominal length $\lambda _b = 1$ , and leg design alone was studied for the two tasks presented before. The scaling of the upper and lower leg link is then modified together with the motor selection. For the scaling, a uniform grid of 50 $\times$ 50 was studied within the range $[0.5, 1.5]$ . The results are depicted in Fig. 11(a) and (b) for the bounding and the backflip tasks, respectively. The plots show the value of the cost against the scaling of the upper link $\lambda _u$ and lower link $\lambda _l$ once a specific motor combination is chosen. From the trends of the optimal value $\mathcal{L}$ , some orthogonality emerges between the two tasks in the design space.

Figure 10. Center: nominal, left: optimized backflip, right: optimized bounding.

Figure 11. Cost landscape for the different motor combinations and link scaling. In (b), white regions are associated with unfeasible problems, which occur consistently when the robots have short shanks and long thighs.

Figure 12. Pareto front approximation for the two tasks’ cost. The designs are visually superimposed in (a). The one highlighted in orange is the design which requires the least modifications to the nominal prototype, shown in red. In the side Table (b) the different points are reported.

With the values obtained from the grid search, the Pareto front is approximated (Fig. 12(a)) for the two different task costs. The result is reported in Fig. 12(b): it constitutes a reduced set of candidates that can perform both bounding and backflip reasonably well. As a second-order criterion, designs that involve fewer modifications to the nominal prototype are preferred. Practical considerations drive this: modifying the shank link is easier than the thigh, as the modification of the latter involves a re-design of the transmission, which is more challenging. Moreover, an optimized robot for bounding is preferred if this implies a sacrifice of performance for backflips (locomotion on the $x$ direction is the main movement mode). So the closest options are the 2 $^{\text{nd}}$ and 3 $^{\text{rd}}$ rows of Table (b) in Fig. 12, which use the same motors and a lower link very close to the nominal one. Among the two, the one with the scaling parameter of the thigh at $\lambda _u=0.66$ was selected. The chosen design decreases the cost for the motion, as shown in Fig. 12(a). The relative improvement of this design with respect to the nominal one is 52% for backflip and 67% for bounding. Some performance was sacrificed for the sake of versatility as, for a single task, we found improvements of 87% and 77% respectively.