Trajectory optimization based on the full dynamics of a robotic system provides a flexible tool to generate complex motion plans for robots. It enables the system to exploit the full dynamic capabilities of the robot to achieve a task. State-of-the-art approaches are able to rapidly find solutions while incorporating increasingly complex model descriptions, which allows using trajectory optimization in a model predictive control (MPC) fashion. However, relying on the specific structure of the model makes implementation of the synthesized motion plans prone to modeling errors. Executing motion plans on hardware has therefore proven to be nontrivial and often requires manual, task-dependent tuning of cost functions and constraints to achieve feasible motions.
A major source of modeling error is the treatment of actuators as perfect torque sources. Any real system is subject to bandwidth limits and as such is not an ideal torque source. Current approaches for trajectory optimization neglect this effect and find motion plans that change torque instantaneously, making it impossible to track the plans on the hardware. A similar modeling error occurs when assuming a rigid contact with the ground, which generally does not hold during locomotion in outdoor environments or on compliant surfaces as shown in Fig. 1. The rigid contact essentially provides the optimizer with infinite bandwidth on the contact forces. As a result, motion plans generated assuming idealized contact and actuator dynamics cannot be tracked by the hardware, leading to failure of the locomotion controller.
In this paper, we extend model predictive control (MPC) methods for legged robots to situations where the assumptions of rigid ground and perfect actuators are invalid. We address the issues of inherent bandwidth limits in real robots by adapting the cost function to be frequency-dependent, making it possible to penalize high frequencies in the motion plans. As a result, instantaneous changes to torque or ground reaction forces (GRF) are avoided. The solver, therefore, does not have to reason about the exact details of terrain and actuator dynamics but will produce solutions that are achievable under the bandwidth limits. We show that motion plans generated with our frequency-aware trajectory optimization can be followed by the hardware more accurately than those generated with a standard baseline and enables locomotion on compliant terrain.
I-a Related work
To increase the feasibility of optimized motions several aspects of the problem can be modified. In general, one can classify such efforts in adaptations of the cost function, model, or constraints. Independent of the strategy, one typically needs to strike a balance between ensuring compatibility with hardware and restricting the solution space.
In case of a series elastic actuator, the dynamics can be approximated and added to the model. The optimization algorithm is in those cases able to exploit the properties of the specific actuator and adding spring-damper elements to the joints is even known to result in motions that resemble those found in nature . For series elastic actuators, methods have been proposed to incorporate bandwidth, torque, and joint limits in a computationally efficient way . However, very often we don’t have exact details of actuators and modeling them would lead to high engineering effort. Moreover, since parts of the underlying actuator dynamics have very different time constants, smaller timesteps leading to slower update rates prevent such models from use in MPC with complex systems. Additionally, stability problems can arise when a limb with stiffly modeled actuators makes contact with the environment . We avoid such issues by not explicitly modelling the actuators, but by incorporating well-known bandwidth limitations up to which the perfect tracking assumption is valid.
While model parameters for actuators can be obtained from first principles or through repeated experiments, contact dynamics are considerably harder to model or predict. A combination of learning a terrain model and trajectory optimization has been proposed in . However, such methods have not reached real-time capabilities, and the question of how to optimize without prior knowledge of an unknown terrain remains unanswered.
In the context of contact invariant optimization, soft contact models  or contact smoothing  are used inside trajectory optimization. These models are selected and tuned for their numerical properties rather than physical accuracy. The models need to be smooth since the highly coupled interaction between stiff contact dynamics and slow dynamics of the robot lead to poor convergence of the algorithms. In the worst case, this numeric model tuning can lead to highly undesired effects when the optimizer starts exploiting dynamic properties of the terrain which are entirely wrong.
Instead of adapting the model, a smooth parameterization of the trajectories can be imposed to ensure feasibility. This common technique is found in spline based optimization [10, 11], collocation methods [12, 13], or when using dynamic motion primitives . This, however, limits the motions that can be expressed and can often require a problem specific, manual tuning procedure.
As an alternative to adapting the model or constraints, we propose to encode bandwidth limits through the cost function. A trivial way to do so is to put extra costs on input signals, but this results in slower response overall and goes against the desire to perform highly dynamic motions at the limits of the system. Instead, we intend to explore a slightly different approach and propose to use a frequency-dependent cost function . We penalize control actions only in the high frequency range and combine this idea with a modern optimal control framework [16, 17].
Perhaps the closest related idea in recent literature is found in , where a constraint is formed in the frequency domain. Since the locomotion problem already has numerous constraints, we prefer to use the cost function approach instead of adding more hard constraints that increase the risk of an infeasible solution. Additionally, the proposed window based constraints approach requires previous and future decision variables, which negatively affects the Riccatti-sparsity pattern exploited in optimal control methods.
We introduce frequency-dependent cost functions integrated into modern MPC strategies for legged locomotion. Through simulation experiments, we study the effect of such a cost function on the resulting solutions. The proposed method provides the user with an intuitive way to achieve robustness against unmodeled phenomena like actuator bandwidth limits and non-rigid contact dynamics. These findings were successfully validated in hardware experiments on different grounds. Using frequency-shaped cost functions, we could improve the robustness of ANYmal while locomoting under substantial external disturbances coming from external pushed or unmodeled soft ground.
First, we discuss uncertainty in the dynamics from a robustness point of view and motivate the particular choice of cost functions. Afterward, the integration with a sequential linear quadratic model predictive control (SLQ-MPC) method  is presented. This method, based on Differential Dynamic Programming (DDP), relies on a linear approximation of the dynamics and a quadratic approximation of the cost function around the latest trajectory.
For brevity of notation in the current section, and without loss of generality, we consider the following quadratic cost function without linear and mixed state-input costs,
where is the positive semi-definite state cost Hessian and is the positive definite input cost Hessian.
Ii-a High frequency robustness
Consider a linear plant with unstructured multiplicative uncertainty model ,
is the maximum singular value of the disturbance model andis a frequency-dependent upper bound. The closed loop stability condition for these models is ,
where is the minimum singular value, and is the transfer function of plant and controller together.
To be robust against large uncertainties at high frequencies, according to (4), the loop gain, , should be kept low. Intuitively, penalizing inputs at the high frequencies reduces the feedback gain at those frequencies, which allows for larger uncertainty magnitude, . We therefore propose to use the following frequency-dependent input weighting
where is the original input cost, and and are the zero and pole of the loopshaping transfer function. A visualization of such cost function is provided in Fig. 2. This cost function is of similar shape as the disturbance model proposed in , which highlights the close relation between disturbance modelling and cost function shaping.
Indeed, for Single-Input-Single-Output (SISO) systems Anderson et al.  established that the open loop gain at high frequency under the frequency-shaped cost function (5) is reduced, i.e. . According to (4), the resulting increase in permits a higher model uncertainty in the stopband. At the same time, some gain and phase margins in the passband are lost, with associated decrease in robustness in the passband.
Moreover, the condition on for the resulting controller to contain lag compensation is that the open-loop plant has zero, or an even number of, positive real poles greater than . The transfer function pole, , in Fig. 2 thus has to exceed the fastest unstable pole of the system.
Unfortunately, to the best of our knowledge, a robustness proof for MIMO systems is not available. Despite that, the intuition that penalizing high frequency input increases compatibility with actuator bandwidth limits remains. In this paper, we aim to empirically validate the effect of using such a cost function.
Ii-B frequency-shaped cost functions
MPC plans over a receding horizon. The cost function in (5) therefore needs to be expressed in the time domain as well. This can be achieved by a state augmentation as described in . We repeat the essential steps for completeness.
is the Hermitian transpose of the vector. Here, it becomes apparent that the standard costs over states and inputs are constant for all frequencies. To leave the possibility of having different loopshaping per input dimension, the frequency-dependent weight matrix in (5) is generalized to
where is the complex conjugate of . Every input directions can now have its own shaping function . In order to transfer this new cost function back into the time domain, a change of variables is required. The filtered inputs, , are defined as
In this work, we keep the original inputs as decision variables and thus require slight costs on
to obtain a positive definite cost matrix. We therefore interpolate between (6) and (11) and obtain
To compute the filtered inputs , the state space realization of the transfer functions is used,
The system dynamics should then be extended to
and cost matrices adapted as
Note that it is also possible to use the state space realization of , without the inversion. This, however, results in faster filter dynamics because would correspond to the fast pole instead of the slower pole of the inverse, . We prefer the slower filter dynamics to allow the continuous time forward integration in the SLQ-MPC to use larger time-steps.
Moreover, (14) can be substituted into and
directly instead of adding it as a constraint. However, keeping both decision variables gives the solver more degrees of freedom in intermediate steps and leaves the original system dynamics and constraints untouched. Since we want to interfere with the nominal problem as little as possible, we opt for the latter. Such a strategy, known aslifting , is a commonly used method to increase the region of convergence and contraction rates of Newton type algorithms.
A detailed study of the numerical properties of each equivalent reformulation is beyond the scope of this paper. We are rather interested in the effect that the cost function shaping has on the solutions of the optimal control problem.
The option to constrain the input rate can also be interpreted as a frequency-dependent weighting, where the weight is given by . Implementing this in an MPC framework requires one to use the original input as a state and have the input rate as the new inputs. This turns state-input constraints into pure state constraints, which are more computationally expensive to handle. Retaining the original input alongside the input rate analogous to (16) does not remedy this since for filters that are not of relative degree 0, which leads to an overconstrained problem.
Iii Experimental setup
Iii-a Problem Formulation
The proposed method is applied to the kinodynamic model of a quadruped robot, which describes the dynamics of a single free-floating body along with the kinematics for each leg . The Equations of Motion (EoM) are given by
where and are the rotation matrix of the base with respect the global frame and the transformation matrix from angular velocities in the base frame to the Euler angles derivatives in the global frame. is the gravitational acceleration in body frame, and
are the moment of inertia about the CoM and the total mass respectively. The inertia is assumed to be constant and taken at the default configuration of the robot.is the position of the foot with respect to CoM. is the orientation of the base in Euler angles, is the position of the CoM in world frame, is the angular rate, and is the linear velocity of the CoM. The inputs of the model are the joint velocity commands and end effector contact forces .
The constraints depending on the mode of a leg at that point in time are formulated as
where is the end effector velocity in world frame, which constrains a stance leg to remain on the ground and a swing leg to follow the predefined curve normal to the surface in order to avoid foot scoffing with zero contact force inputs .
The baseline cost is formulated as a quadratic function
where is a desired state consisting of commanded base pose and twist by the user and a default configuration for the joints. We use a simple diagonal cost on all state variables and contact force inputs. The joint velocities costs are a diagonal matrix pre and post multiplied by the end-effector Jacobians to weigh the costs in task space.
Iii-B System integration
In the model described in the previous section, the control inputs are end-effector forces and joint velocities. To translate the solution to torque commands, we extract a full position, velocity, and acceleration plan for CoM and end-effector trajectories, in addition to the planned contact forces. This plan is tracked by the hierarchical whole-body control (WBC) architecture described in . The tasks in decreasing priority are (1) satisfying the equation of motions and zero acceleration for contact feet (2) tracking CoM and swing leg trajectories, and (3) tracking the planned contact forces. The desired contact forces from the MPC are thus communicated in two ways: The CoM trajectory dictates the net acting forces, and force tracking task regulates the internal forces. Without the latter, contact forces would be redistributed among the contacts, which would override the planned smoothness of the trajectories. Additionally, on all priorities, torque limits and friction cone constraints are imposed as inequality constraints.
The optimal control problem in (18) is solved for a user-defined gait with the continuous time SLQ-MPC algorithm described in . We use a receding horizon length of s, which results in an MPC update rate of Hz for the baseline method and around Hz for the frequency-shaped method.
The MPC runs on a desktop PC with an Intel Core i7-8700K CPU@3.70GHz hexacore processor and continuously computes a motion plan from the latest known state through a real-time-iteration scheme. The WBC runs on the dedicated onboard PC and tracks the most recent plan. Here, the augmented filter state is propagated as well according to (13) with the currently available augmented input plan . Both nodes communicate over a local ROS network. An overview of the experimental setup is provided in Fig. 3.
We study the effects of adopting the cost function in (12) for the previously described setup under various locomotion tasks. To see the results at different levels of model errors, we conduct perfect model simulations, rigid-body simulations, and hardware experiments. When selecting different values for , is selected such that the frequencies in the stopband incur a cost of times the DC cost, i.e. .
Iv-a Perfect model
First, we investigate the effect of the loopshaping on the contact forces in a simulation that uses the same model as the MPC. This shows how the resulting trajectories are different already in the case of no modeling errors. The analyzed gait is a trot with a duty factor of 0.5, i.e. with no overlap in stance phases of the diagonal feet, while a forward velocity of is commanded. Figure 4 shows the ground reaction forces when selecting different values for . As seen from the plot, the baseline method instantaneously applies contact forces once a foot is in contact. As expected, lowering the frequency at which costs start to increase, i.e. lowering , results in increasingly smooth trajectories. The frequency-shaped method approaches the baseline as goes to infinity. The corresponding base height trajectories are plotted in Fig. 5. Smoother contact force trajectories require more vertical displacement of the base, while the baseline produces the exact amount of force to keep the base level.
A similar experiment is conducted to see how frequency-shaped solutions behave for different gait periods. The period of the trotting gait is varied for a fixed cost function with . The resulting ground reaction forces are shown in Fig. 6. For the fixed smoothing parameter, the contact profile transitions from a two hump solution to a single hump solution as the gait period decreases. Interestingly, a similar pattern is observed for humans and quadrupeds .
Iv-B Physics simulation
The combination of tracking controller and MPC is evaluated in the Open Dynamics Engine  rigid-body simulation, where we can vary ground properties in a controlled way. The model errors, in this case, come from the difference between rigid-body dynamics and the kinodynamic model, as well as the assumption of a rigid ground contract when the terrain is made compliant. ODE allows for simulation of soft contacts by modelling111ODE relaxes the rigid-contact solver such that it implicitly resembles a spring-damper interaction. the ground contact forces as a spring-damper system. Three different sets of spring-damper parameters and are selected to simulate hard, intermediate, and soft ground, respectively. For each terrain, three different cost functions are evaluated: the baseline without frequency shaping as well as frequency-dependent cost functions with and . These values were selected based on Fig. 4 to represent three levels of smoothness in the continuum of available cost functions. The resulting contact force profiles for all combinations are shown for a single stance phase in Fig. 7. Desired and measured contact forces are shown for a single leg during an in-place trotting motion with a stance duration of .
As the compliance of the terrain increases, the difference between desired forces generated from the WBC and resulting forces grows. The WBC uses rigid-body dynamics with a hard contact assumption to compute desired contact forces. The difference between desired and measured forces is therefore a measure of disturbances inserted by additional unmodeled dynamics, which in general includes the bandwidth limits of actuators and contact dynamics that we aim to avoid. Table I shows the force tracking performance averaged over six gait cycles in Mean Absolute Error (MAE) and Mean Squared Error (MSE) defined as
For all cost functions, tracking performance degrades as the model error increases. Qualitatively, the baseline controller suffers from a larger error at the beginning and end of the contact phase, due to its step-like change in the desired forces. Even on hard terrain, there is an apparent benefit of loopshaping. The smoother transition between contact phases mitigates the disturbance generated by contact timing mismatch. Differences become most apparent for the soft terrain, where the smoother trajectory has better performance in especially in MSE.
|Hard||6.1 (751.9)||4.1 (256.9)||4.2 (230.2)|
|Medium||7.0 (353.4)||4.7 (114.3)||5.3 (130.8)|
|Soft||17.7 (742.0)||17.4 (543.9)||11.9 (251.7)|
Furthermore, we examine the locomotion strategy under two extrema of cost functions. The commanded forward velocity during a trotting motion is gradually increased until failure occurs. The foot placement strategies are visualized in Fig. 8. The plots show footstep locations from a top-view with the robot starting at the origin. As the velocity increases towards the right side of the plot, the footstep locations start to differ. Interestingly, with the smoother cost function of , the foot placement strategy is significantly altered and becomes velocity depend. Where under the baseline costs the controller chooses a fixed width foot placement, the frequency-shaped solution places the feet increasingly inwards for higher velocities. Remarkably, this results in a significantly higher maximum velocity of versus .
On hardware, we aim to validate the simulation results for contact force tracking performance on different terrains. The floor of the lab, a foam block, and a mattress are selected to test a rigid, an intermediate and a very compliant terrain. A force-torque sensor is mounted on the right front leg to obtain direct measurements of the ground reaction forces. The resulting measured and desired forces for different cost function and terrain combinations are shown in Fig. 9. The plots show the difference between measured and desired contact forces of the right-front leg during the first three steps. The MAE and MSE averaged over those first three gait cycles are given in Table II. On hard terrain, all methods perform well, and slight differences in tracking performance occur at the beginning and end of the contact phase. In this area the controller provides the smoothest transition, resulting in the best MAE and MSE on all terrains. That controller, however, has in return less accurate tracking during the middle of the stance phase. Where in the simulation experiments we see that a medium amount of smoothing is best for hard and medium stiff terrain, we do not see this in the hardware experiments. The difference comes from the fact that ANYmal has series elastic actuators, causing model errors and bandwidth limits even on hard terrain. For a step input to the torque level reached during trot, a rise-time up to , equivalent to a bandwidth of around can be expected .
When further reducing the compliance by trotting on the mattress, the baseline and controllers suffer a substantial decrease in performance. The base height during the first part of the experiment is shown in Fig. 10. Due to the significant mismatch between the planned and resulting contact force with each footstep, the baseline controller loses base height in a few steps, causing it to fail. The cost function does achieve a trot, as shown in the video222https://youtu.be/Mk_AjNnp6t0, but strong oscillations are present between the terrain and the feet. The cost function, finally, results in both a stable trot and a smooth contact interaction.
|Hard||25.1 (2597.4)||21.2 (1696.8)||20.3 (1303.3)|
|Medium||22.1 (2044.6)||19.5 (1387.3)||18.5 (998.0)|
|Soft||37.9 (3470.6)||30.7 (2017.2)||22.1 (1022.4)|
In the accompanying video, we additionally show the behavior under disturbances. The robot trots in place and has costs on base deviation from the initial position. We disturb the robot in the horizontal plane. Qualitatively, the reactive stepping and push-back behavior differ. Under the baseline method, the robot reacts to a push by generating lateral forces and the user experiences instant resistance. The frequency-aware implementation with instead accepts deviation of the base trajectory and adapts future step location and force profile to smoothly return to the origin.
We have shown that a single parameter of a frequency-dependent cost function provides a handle on a rich variety of solutions. While the smoothness at the beginning of the stance phase is not surprising due to the filter, the anticipatory decrease in force before lifting the foot, as seen in Fig. 4, shows that the filter and planning are tightly working together. The additional filter states allow the Riccati-type algorithm to reason about future state-input constraints, in this case, zero contact forces during the swing, and adapts the strategy to approach them smoothly. This is remarkable because the backward pass is projected on those constraints only at the point in time that they are active. Interestingly we also find that the foot placement strategy changes significantly. These observations highlight the fact that embedding frequency awareness of the MPC is richer than simply filtering the obtained inputs. As a future work, we will explore ways to change to cost function online and in this way adapt the locomotion strategy to the terrain.
A further opportunity arises when the gait is part of the motion optimization. As shown in Fig. 6, there is an interplay between the gait and the obtained contact force profiles. frequency-dependent costs could therefore be used to drive gait selection through optimization.
We introduced frequency-aware MPC by combining frequency-dependent cost functions with modern MPC methods. With simulation experiments, we show that the resulting smoother force profiles improve tracking performance when the rigid terrain assumption is relaxed, without the need to explicitly model it. We validated these results on hardware and see a similar trend when comparing performance on various terrains. The method is shown to provide robustness against unmodeled dynamics of series elastic actuators and compliant terrains. We demonstrated that with this approach ANYmal is now able to execute dynamic motions even on highly compliant terrains.
-  M. Hutter, C. Gehring, D. Jud, A. Lauber, C. D. Bellicoso, V. Tsounis, J. Hwangbo, K. Bodie, P. Fankhauser, M. Bloesch, et al., “ANYmal - a highly mobile and dynamic quadrupedal robot,” in IROS, 2016, pp. 38–44.
-  G. Schultz and K. Mombaur, “Modeling and optimal control of human-like running,” IEEE/ASME Transactions on Mechatronics, vol. 15, no. 5, pp. 783–792, Oct 2010.
-  J. Nakanishi and S. Vijayakumar, “Exploiting passive dynamics with variable stiffness actuation in robot brachiation,” in Proceedings of Robotics: Science and Systems, Sydney, Australia, July 2012.
-  D. J. Braun, F. Petit, F. Huber, S. Haddadin, P. van der Smagt, A. Albu-Schäffer, and S. Vijayakumar, “Robots driven by compliant actuators: Optimal control under actuation constraints,” IEEE Transactions on Robotics, vol. 29, no. 5, pp. 1085–1101, Oct 2013.
-  R. Schlossman, G. C. Thomas, and L. Sentis, “Exploiting the natural dynamics of series elastic robots by actuator-centered sequential linear programming,” arXiv preprint arXiv:1802.10190, 2018.
-  J. Koenemann, A. D. Prete, Y. Tassa, E. Todorov, O. Stasse, M. Bennewitz, and N. Mansard, “Whole-body model-predictive control applied to the hrp-2 humanoid,” in IROS, 2015, pp. 3346–3351.
-  A. H. Chang, C. M. Hubicki, J. J. Aguilar, D. I. Goldman, A. D. Ames, and P. A. Vela, “Learning to jump in granular media: Unifying optimal control synthesis with gaussian process-based regression,” in ICRA, May 2017, pp. 2154–2160.
-  M. Neunert, F. Farshidian, A. W. Winkler, and J. Buchli, “Trajectory optimization through contacts and automatic gait discovery for quadrupeds,” Robotics and Automation Letters, vol. 2, no. 3, pp. 1502–1509, July 2017.
-  I. Mordatch, E. Todorov, and Z. Popović, “Discovery of complex behaviors through contact-invariant optimization,” ACM Transactions on Graphics (TOG), vol. 31, no. 4, p. 43, 2012.
-  M. Kalakrishnan, J. Buchli, P. Pastor, M. Mistry, and S. Schaal, “Fast, robust quadruped locomotion over challenging terrain,” in ICRA, May 2010, pp. 2665–2670.
-  C. D. Bellicoso, F. Jenelten, C. Gehring, and M. Hutter, “Dynamic locomotion through online nonlinear motion optimization for quadrupedal robots,” IEEE Robotics and Automation Letters, no. 99, 2018.
-  A. Hereid, E. A. Cousineau, C. M. Hubicki, and A. D. Ames, “3d dynamic walking with underactuated humanoid robots: A direct collocation framework for optimizing hybrid zero dynamics,” in ICRA, May 2016, pp. 1447–1454.
-  D. Pardo, M. Neunert, A. Winkler, R. Grandia, and J. Buchli, “Hybrid direct collocation and control in the constraint-consistent subspace for dynamic legged robot locomotion,” in RSS, 2017.
-  A. Werner, B. Henze, F. Loeffl, S. Leyendecker, and C. Ott, “Optimal and robust walking using intrinsic properties of a series-elastic robot,” in 2017 IEEE-RAS 17th International Conference on Humanoid Robotics (Humanoids), Nov 2017, pp. 143–150.
-  N. K. Gupta, “Frequency-shaped cost functionals-extension of linear-quadratic-gaussian design methods,” Journal of Guidance, Control, and dynamics, vol. 3, no. 6, pp. 529–535, 1980.
-  F. Farshidian, M. Neunert, A. W. Winkler, G. Rey, and J. Buchli, “An efficient optimal planning and control framework for quadrupedal locomotion,” in ICRA. IEEE, 2017, pp. 93–100.
-  F. Farshidian, E. Jelavic, A. Satapathy, M. Giftthaler, and J. Buchli, “Real-time motion planning of legged robots: A model predictive control approach,” in 2017 IEEE-RAS 17th International Conference on Humanoid Robotics (Humanoids), Nov 2017, pp. 577–584.
-  J.-H. Hours, M. N. Zeilinger, R. Gondhalekar, and C. N. Jones, “Constrained spectrum control,” IEEE Transactions on Automatic Control, vol. 60, no. 7, pp. 1969–1974, 2015.
-  J. Doyle and G. Stein, “Multivariable feedback design: Concepts for a classical/modern synthesis,” IEEE transactions on Automatic Control, vol. 26, no. 1, pp. 4–16, 1981.
-  F. Farshidian, E. Jelavić, A. W. Winkler, and J. Buchli, “Robust whole-body motion control of legged robots,” in Intelligent Robots and Systems (IROS), 2017 IEEE/RSJ International Conference on. IEEE, 2017, pp. 4589–4596.
-  B. Anderson and D. Mingori, “Use of frequency dependence in linear quadratic control problems to frequency-shape robustness,” Journal of Guidance, Control, and Dynamics, vol. 8, no. 3, pp. 397–401, 1985.
-  B. D. O. Anderson, J. A. Gibson, and H. R. Sirisena, “Phase lag and lead weighting in linear-quadratic control,” Optimal Control Applications and Methods, vol. 6, no. 3, pp. 249–263.
-  J. Albersmeyer and M. Diehl, “The lifted newton method and its application in optimization,” SIAM Journal on Optimization, vol. 20, no. 3, pp. 1655–1684, 2010.
-  C. D. Bellicoso, C. Gehring, J. Hwangbo, P. Fankhauser, and M. Hutter, “Perception-less terrain adaptation through whole body control and hierarchical optimization,” in Humanoid Robots, 2016 IEEE-RAS 16th International Conference on. IEEE, 2016, pp. 558–564.
-  Z. Gan, T. Wiestner, M. A. Weishaupt, N. M. Waldern, and C. D. Remy, “Passive dynamics explain quadrupedal walking, trotting, and tölting,” Journal of computational and nonlinear dynamics, vol. 11, no. 2, p. 021008, 2016.
-  R. Smith et al., “Open dynamics engine,” http://www.ode.org, 2005.