Over the past years, there has been a significant shift in research attention from industrial robotics towards the development of robots that are capable of a more dexterous and dynamic interaction with their environment. Such an interaction generally involves non-smooth dynamics that can arise from contact constraints, impact events, or frictional forces. Locomotion and dynamic manipulation lie at the forefront of such applications, where analyzing and exploiting the dynamics induced by contacts is fundamental for achieving the desired task at hand. The associated hybrid system dynamics make the design of motion plans and control policies an intricate process.
Some of the most prominent control techniques developed to deal with such problems rely on extensions of standard trajectory optimization (TO) methods (direct collocation, single shooting, multiple shooting…) [2, 8, 6]. These extensions incorporate concepts behind the modelling and simulation of non-smooth contact phenomena and impact dynamics. The primary integration schemes used for the forward simulation of hybrid systems can be identified as the event-driven or the time-stepping (event-capturing) techniques. The former solves for the different dynamic modes while requiring, at the same time, an event-detection plan that would indicate a transition from the current mode to the next one. On the other hand, time-stepping schemes directly discretize the overall non-smooth dynamics without the need for any detection-step. Adapting these concepts from simulation to control purposes yields two main classes of TO control methods: the multi-phase approach and the contact-implicit approach, which are analogous to the aforementioned event-driven and time-stepping schemes, respectively.
The multi-phase approach entails an a priori knowledge of the continuous phases or modes, as well as their sequence of occurrence over the full time horizon. The switching times for the discrete transitions can be a result of the optimization program. Researchers such as Patterson and Rao  have developed a MATLAB program (GPOPS-II) that utilizes the multi-phase strategy as a framework for solving generic hybrid optimal control problems. Other successful applications are repeatedly encountered in the legged-robotics community, such as for bipedal robot locomotion , and for quadrupedal locomotion . However, complications with the use of a predefined event-schedule begin to arise as the number of discrete transitions increases; and that is simply due to the exponential growth of the possible multi-phase sequences with respect to these transitions.
A remedy to this issue can be found in the mode-invariant approach, also referred to as contact-invariant, contact-implicit, or through-contact approach. The fundamental idea behind it is that, instead of having a predefined mode sequence along with inter-phase constraints, discrete transitions are not given any special treatment; hence, the contact constraints are enforced at every single grid-point. The method is also capable of incorporating discontinuous jumps in the states – such as in velocities during impact events – by treating the whole problem on an impulse level, thereby not differentiating between finite forces and impulsive forces.
Both multi-phase and contact-implicit schemes can be implemented along with a variety of underlying design choices. To elaborate, one could flexibly choose any of the direct TO methods, with either a soft-contact model or a hard-contact model assumption; and in case the latter was chosen, one would have to deal with the resulting set-valued force laws by either resorting to an augmented-Lagrangian approach  or by solving a linear complementarity problem (LCP) 
. Furthermore, one could decide to add the contact forces into the vector of decision variables, or to resolve them separately from the optimization iterations.
Some remarkable work has been done in the computer graphics and the robotics communities on optimal motion planning and control, under the framework of contact-implicit optimization (CIO): For instance, Mordatch et al. [13, 14, 12] use a contact-invariant approach for producing physics-based animations of anthropomorphic, contact-rich tasks like walking, climbing, crawling, moving objects, and multi-fingered manipulation. Although some variations of their original CIO approach are implemented, there is a common fundamental treatment of the contact constraints, which are imposed through penalty terms added to the objective function (soft-constraints); and this relaxation, while it certainly aids the gradient-based solver in efficiently converging to a local minimum, it also tolerates non-physical behavior caused by allowing contact forces to act at a distance. Another soft-constrained CIO formulation is achieved in  by a differential-dynamic-programming (DDP) optimal control method, where the hard unilateral contact constraints and frictional forces are dealt with by solving a stochastic linear complementarity problem rather than an LCP. With this, they get rid of the discontinuous dynamics and are left with a continuously differentiable system, while retrieving an optimal controller that is also robust against uncertainties. Neunert et al.  have developed a nonlinear model predictive control (NMPC) algorithm that solves, online and at high rates, the nonlinear optimal control problem using a generalization of the iLQR method. The NMPC scheme is tested on a quadruped, while adopting a contact-implicit formulation which assumes a soft-contact model with smoothing elements. On the other hand, Carius et al.  rely on the augmented-Lagrangian approach, specifically the Gauss-Seidel method combined with a proximal point projection (SORprox), in order to solve the set-valued force laws arising from the hard frictional contacts. This is performed as part of a Moreau time-stepping scheme which is used to simulate the non-smooth dynamics within an iterative-LQR (iLQR) optimal control plan.
In this paper, we propose a variation on the CIO formulation introduced by Posa et al. . Their method is originally motivated by the Stewart and Trinkle time-stepping scheme , where multi-contact dynamics are formulated as a linear complementarity problem (LCP). They make use of the simultaneous nature of direct collocation by incorporating the contact forces as optimization variables, while additionally imposing the unilateral hard-contact constraints in the form of complementarity conditions. In fact, the transcription of hybrid optimal control problems into finite-dimensional optimization programs including complementarity conditions has been initially proposed by Yunt Glocker ; but unlike , the optimal control inputs are retrieved separately from the contact forces and states through a bilevel optimization. We follow the same fundamental steps, while reformulating the nonlinear program (NLP) within the multiple-shooting framework instead. This modification is not immediately straightforward as it sacrifices the method’s physical correctness, unless some extra steps associated with the introduction of an impact-law are considered. The new NLP-structure enables us to exploit the computational power of the FORCES Pro  commercial tool, which generates highly customized and efficient optimization solvers based on the primal-dual interior point algorithm. Its effectiveness is mainly a result of an efficient linear-system solver that exploits the multi-staged problem-structure, in order to solve the resulting Karush-Kuhn-Tucker (KKT) system . we demonstrate the use of this method for solving dynamic object manipulation problems similar to that tackled in . Unlike , wherein Posa’s CIO scheme is also tested, we achieve optimal state-trajectories – at a significantly lower computational time – that are verified to be dynamically feasible through simulation as well as hardware experimentation. In fact, relatively simple examples in  are solved within a few minutes, compared to average convergence times ranging between 50 to 800 milliseconds for our setups. The resulting motion plans are shown to adhere to the imposed contact conditions, and result in realistic contact forces. These indeed form the crucial ingredients for successfully attaining the desired manipulation-task goal, reliably and accurately, without the need for any sort of post-optimization modifications.
Ii Problem Description
Ii-a Dynamic Object Manipulation
Dynamic manipulation belongs to the family of non-prehensile (graspless) manipulation, which generally includes tasks such as pushing, tumbling, throwing, and catching [10, 11]. Essentially, what makes the task “dynamic”, is that, unlike the quasi-static case, inertial effects have to be taken into account due to the relatively fast motions; and secondly, it is possible for the manipulator to lose contact with the object during the assigned manipulation period. This, in turn, makes it mandatory for one to reason about the dynamics induced by the environment on the uncontrollable object, while planning and controlling the robot motion that would satisfy the objective.
We consider a variety of dynamic pushing problems in this work. Ultimately, our goal is to formulate an NLP that would allow a six degrees of freedom (DOF) robot manipulator (section II-B) to dynamically push a block from a given initial state, to a final desired position that could lie outside the robot’s workspace (see fig. 1). The problem formulation is general enough to work on diverse manipulation problems as well, such as pushing a door open, or pushing a ball into a desired direction and up to an assigned goal (see fig. 2). We elaborate on the problem specifics and results of the block-task, as it comprises all relevant effects from a modelling and computational perspective, due to the action of the environment on the block in the form of dry-friction.
We test the proposed framework with ANYpulator [3, 5], a 6-DOF robot arm depicted in fig. 3. It comprises a set of robust and mechanically compliant series-elastic actuators (ANYdrive ) that feature precise position and torque control as well as high impact robustness. The six joints are responsible for: Shoulder Rotation (), Shoulder Flexion/Extension (), Elbow Flexion/Extension (), Wrist Flexion/Extension (), Wrist Deviation (), and Wrist Pronation (). For the sake of our dynamic pushing applications, a custom-built end-effector is used for all the object-manipulation problems. The specialized tool basically constitutes of a cylindrical rod with a spherical head, as shown in fig. 3. Consequently, the rotation of the last joint would have no contribution to the resulting solutions of the applications at hand, thus allowing us to extract a 5-DOF dynamic model from the manipulator equations of motion. Such a consideration helps decrease the optimization program-size significantly, which in turn reduces the associated computational effort.
Tackling the aforementioned manipulation problems with the contact-implicit approach requires an accurate dynamic model of the robot and object. We consider a frictionless end-effector to object contact in the following derivations, since the tangential components of the contact forces/impulses only play a minor role in our dynamic pushing tasks. This assumption helps reduce computational time at the expense of a slight loss in physical accuracy.
Iii-a Model Formulation
Noting that in our applications the generalized velocities coincide with the time-derivatives of the generalized positions , the fully-actuated dynamics of the robot are concatenated with the non-actuated dynamics of the object to form the following underactuated rigid body dynamic equations:
as well as the following impulse-momentum equations:
where the subscripts r and o refer to quantities corresponding to the robot and the object, respectively. Letting , is the generalized mass matrix, while contains Coriolis and centrifugal terms, as well as gravitational and viscous damping terms. The selection matrix maps joint input torques to generalized forces, while represents the generalized forces applied by the environment on the manipulated object, such as dry friction. We define the gap function , which is the signed distance between potentially colliding bodies. In our case, is a scalar-valued function indicating the minimum distance between the robot’s end-effector and the object. Consequently, the normal Jacobian and the magnitude of the normal contact force are defined accordingly. Moreover, the post- and pre-impact generalized velocities are given by and in (2), while and are the normal contact impulse-magnitude and external impulse, respectively.
Iii-B Nonlinear Program Setup
The generic finite-dimensional optimization program is defined using the following multi-staged framework (compatible with that of FORCES Pro’s high-level interface):
where is the vector of optimization variables at stage ; and are the inter-stage cost and the terminal cost functions, respectively. , , are selection matrices and is a state-transition function which maps the states from stage to stage ; it could be given as a result of any desired integration scheme, but to be more general, this function could also consist of other nonlinear equality constraints.
With the adopted CIO approach, we formulate and solve a discrete-time optimal control problem using the multiple-shooting method, along with a set of complementarity constraints. These are added to ensure that the unilateral contact conditions hold at every discrete grid-point. Such a setup results in a so-called mathematical program with complementarity constraints (MPCC).
We can now indicate and discuss the various components required to build our nonlinear program. Most of the upcoming considerations are generally applied to solve all of our manipulation tasks; nonetheless, there are still some problem-specifics associated with the object’s geometry and its configuration-space topology. As previously mentioned, the emphasis will be made on the ANYpulator-Block manipulation problem.
Transcription parameters: When choosing the discretization parameters, there is an inherent trade-off between computational speed and integration accuracy. We attain a fair compromise by using a time-horizon of and setting the number of stages ; hence, leaving us with a time-step size of .
State-transition function: The forward simulation of the dynamics is carried out on the basis of a semi-implicit Euler scheme. Similar to the classical explicit Euler method, it is a first-order integrator; however, it additionally holds inherent energy-conservation properties that make it more reliable. Therefore, the nonlinear equality constraints responsible for satisfying the dynamics are given as follows:
It is important to note the difference between , appearing in (4) and , in (1). The latter quantities are infinite during an impact, unlike and which are defined in relation to the normal percussion and frictional percussion – which arise as a result of discretizing the equality of measures  – as follows:
where is defined in terms of and , while in terms of and :
is such that at the time of a discontinuity and is zero elsewhere.
From here on out, the words contact force and frictional force refer to the quantities introduced in (4) and not (1), unless otherwise stated.
Optimization variables: The vector of optimization variables is given by the robot’s joint input torques, the system’s states, the normal contact force , and the dry friction force acting on the block by the environment. We define where is a unit vector indicating the assigned direction of motion for the block, expressed in the robot base frame, which is also used as our inertial frame of reference.
Therefore, where , , and , .
Upper/Lower bounds: The maximum joint velocity is set to , while the maximum allowable torque as . The joint position limits are also imposed in the program.
Boundary conditions: Initial boundary conditions are used to specify the starting configuration of both the manipulator and the object, while the satisfaction of the assigned manipulation task is ensured through final boundary conditions. It is worth noting that in the case of the block, the configuration variables are , and the associated manipulation task involves specifying values for and to indicate the desired displacement and direction of travel, respectively. This is done while also requiring that the block does not rotate at all throughout the full time-horizon.
Objective function: A common objective function that aims at minimizing energy and also puts a cost on the wrist joints velocities, is used. The latter term is introduced because it could be noticed that an unnecessary movement of the two wrist joints occurs if we only minimize energy, as their motion does not contribute much to the total energy expenditure. Furthermore, the torques and joint velocities are normalized as they do not have the same units. Therefore, the normalized cost function is given as follows:
where and are weighting matrices.
Complementarity conditions: The following constraints ensure no penetration between colliding bodies, no forces acting at a distance, and only allow forces that can push and not pull:
A subtle yet fundamental difference between the formulation in this paper and  can be seen in the relative indexing of and in (4) compared to (8). To elaborate, in  the state-transition function was obtained by applying a zero-order-hold (z.o.h) on instead of . Consequently, “the contact force over the whole time-step can be non-zero if and only if the distance function is zero at the end of the interval” . Whereas here, the condition holds if and only if the gap is closed at the start of the interval. This is indeed problematic since there is no more guarantee that the contact event is held over the entirety of the time-step, even though the contact force keeps acting unjustifiably at a non-zero distance and would only vanish at the end of it. Moreover, the complementarity conditions, which are purely imposed on a position-level, only ensure that the resulting normal contact force is large enough to prevent penetration in the next step. Nevertheless, there is certainly no upper limit on how big this force can be, thus indicating that it is quite arbitrary and merely directed towards satisfying the position-based conditions. To resolve this issue, we resort to the theory underlying the simulation of dynamics involving hard contacts and impact events, using the time-stepping technique. What we extract from that is the need for an impact law, such as Newton’s law of restitution:
where the constant is Newton’s coefficient of restitution, while is the relative separation velocity. Noting that (9) includes the case of a non-participating (superfluous) contact , which does not occur in any of our applications, we can reduce this expression to:
where the states at the end of the interval and can be retrieved from (4). A restitution coefficient of is adopted to ensure a purely inelastic collision event, thus effectively imposing a no-rebound condition.
We assume a static Coulomb model to describe the frictional force arising between the plane and the block
where is a positive scalar quantity (unidirectional displacement due to pushing and no pulling) that signifies the velocity along the desired direction of motion and is the static friction force or the breakaway force .
For the sake of avoiding the conditional if-statements appearing in (11), we introduce the following constraints instead:
Generalized distance function:
A generalized gap function is introduced for the ANYpulator-Object problem, while further adding geometry-dependent constraints to the NLP. The minimum signed distance between the robot’s end-effector and any convex-shaped object is defined as the projection of the position vector onto the object-frame’s -axis (refer to fig. 4). Therefore,
where is a body-attached frame (rotates and translates with the object), and its -axis is directed along the assigned direction of motion, while represents the orientation of frame with respect to frame . In order to ensure that the robot’s end-effector actually comes in contact with the block-face – and not only with the separating virtual plane – the following nonlinear inequality conditions are added
where and are the block-face’s width and height, respectively. It is assumed here that the block’s bottom surface and the robot base are at the same level.
Other considerations: An extra path constraint is introduced in order to ensure that the resulting optimal motion plan in joint space does not lead to collisions between the end-effector and the ground plane. Furthermore, the gradients of both the objective function and the nonlinear inequality constraints are provided analytically to reduce the computational load of the program.
Iii-C Control Policy
The output obtained from solving the optimization problem formulated in section III-B can be visualized to verify the satisfaction of the complementarity conditions, and to validate the attainment of the desired manipulation task. Before testing the motion plan on the real setup, we perform simulations within the Gazebo simulation environment , with Open Dynamics Engine (ODE) as its underlying physics engine (it consists of a rigid body dynamics module in addition to a collision detection engine).
The shift from pure visualization to either simulation or experimentation requires the introduction of a closed-loop control policy that is capable of tracking the optimal joint states. A simple yet robust approach is adopted in this work, where the optimal control problem is solved only once; then the resulting open-loop torque sequence is added as a feedforward term to a feedback term acting on the deviations from the optimal position and velocity trajectories. Hence, the control law is derived as follows,
where the index refers to an element from the set of actuated degrees of freedom; while and are the scalar feedback PD-gains. is obtained by applying a zero-order-hold on the sequence , while and
by linearly interpolating the grid-point-valuesand respectively.
A C++ class structure is developed – along with the use of libraries and tools provided by the Robot Operating System (ROS) middleware – to send commands to the ANYdrives with an update frequency of . These commands include the optimal feedforward torques, reference joint positions, and reference joint velocities. The PD gains are manually tuned with the following values during simulation: , ; and during experimentation: , .
Iv Experimental Setup
To have a somewhat accurate placement of the initial block pose, its starting position and orientation are provided as target references for the robot’s tool frame (frame ), which are then tracked with an inverse kinematics controller, to place the block accordingly. Moreover, for the sake of comparing the contact forces generated while performing our experiments to the forces predicted by the optimization, a force/torque sensor is mounted on the robot end-effector.
Several test runs are performed on variations of the dynamic block-pushing problem (i.e., different robot initial conditions, different block starting configurations, and different final desired positions for the block), both in simulation and hardware experimentation. We mainly emphasize three experiments that are tested five times each, and then analyze the results and draw conclusions from them. These experiments involve the same initial block pose and robot configuration; however, the desired displacement of the block varies among them, with experiments 1, 2, and 3 corresponding to displacements of , , and respectively.
The formulated CIO program needed to solve the ANYpulator-Block task consists of 896 optimization variables, 819 equality constraints, and 1955 inequality constraints. On average, it takes the solver around to converge to a local minimum, but this value can vary depending on the utilized initial guess, and on the manipulation-problem specifications. Proper tracking of the optimal reference joint-positions as well as the corresponding end-effector references, is demonstrated for the first experiment in fig. 5 and fig. 6, respectively. It could be noted from the cartesian coordinates of the tool that the robot indeed reaches the block’s initial position and tends to sustain the closed contact (the push occurs along the negative -axis of the inertial frame) as much as specified by the optimal plan, while also maintaining a tool height greater than throughout the full time-horizon. Moreover, the effect of the impact-event can be directly perceived in these plots as well, especially in the joint, where it is evident that a disturbance is introduced after around 0.5 seconds. Indeed, fig. 8 shows this to adhere with the optimal force plan, which involves a contact event occuring as a single impulsive push with a force of about acting at for . More interestingly, in the second experiment, we perceive a contact event that lasts for approximately (see fig. 8), while in the last one we encounter two separate contacts – one with a small magnitude and the other with a large one, held over a single time step each – as shown in fig. 9. The generated normal force is shown to decrease with the assigned displacement along the three experiments.
The resulting optimal trajectories associated with the contact force, indicate the richness of the CIO approach in general, as it yields contact schedules that could not be easily considered or planned for, within a multi-phase framework. In addition to that, the dynamic feasibility of the applied method is demonstrated in three different settings, by observing the properly-tracked motion plan, and by noting that the predicted force sequence is indeed compatible with the force sensor readings (see fig. 9), which is significant given that our control law does not include any closed-loop force tracking term. These two factors are fundamental for achieving the assigned manipulation task at hand, as presented for the third experiment in fig. 7.
One should keep in mind that eventhough the optimal and measured force peak magnitudes in the contact force plot are not equivalent, the resulting normal percussions turn out to be reasonably similar. It is also important to recall that each experiment is repeated 5 times in order to demonstrate the method’s repeatability and reliability. Indeed this is shown to be true, since the same results come up repeatedly, for various runs of the same experiment. To add an accuracy-based evaluation of the method, we refer to the plot shown in fig. 10, which indicates the final error (in absolute value) between the real block position and the desired one. First of all, it could be noticed that in all of our experiments, the block always travels a distance that is less than the expected one. This is likely related to the inaccuracy in our dynamic model in general and specifically the assumed static friction model, which seems to have underestimated the true frictional effects. Secondly, we observe a trend indicating that the higher the desired block displacement, the more erroneous is the outcome. This could be partly related to the fact that the higher the demanded displacement, the larger the generated contact force; and since there could be a slight deviation from the supposed contact location (the central line of the box), then this would lead to the block’s rotation – which is proportional to the contact force’s magnitude – along with its translation. As a result, it loses a part of its total kinetic energy, thus having less energy remaining for the translational motion. Moreover, the neglected frictional forces acting in the robot-to-object tangent contact plane may as well lead to the block’s undesired rotation. Even if the peak in the optimal force plan does not turn out to be larger (i.e., for a different local minimum) when assigning a greater displacement, contact is now established earlier, and for a longer time. Otherwise, a bigger displacement could not be attained by the end of the same time horizon. Consequently, the inaccuracy increases due to the greater accumulated integration error associated with the block dynamics.
The aim of this work has been to make use of a contact-implicit optimization approach to tackle a variety of dynamic pushing problems, performed by the ANYpulator robot-arm, with particular emphasis on the block-task. This was done by formulating the manipulation task as a multi-stage program, thus making it possible to integrate it within the FORCES Pro framework. The NLP is solved efficiently, and yields both optimal state and input trajectories, which, when used as components of a simple closed-loop control scheme, result in reasonable tracking of the reference motion and contact force. The resulting optimal trajectories indicate the richness of this approach in terms of its ability to discover different contact-schedules, while the experimental findings suggest the effectiveness of the control technique in terms of its simplicity, the achieved dynamic feasibility, and repeatability. As for the task-satisfaction accuracy of the method, it has been shown to be adequate, but decreases with the assigned block displacement.
There are several directions for potential future work and considerations, most notably is the further improvement of the method’s computational efficiency, in order to incorporate it within an NMPC framework. Possible ways of achieving that could be by either introducing effective warm-starting techniques for the currently adopted interior-point algorithm, or by exploring different algorithms in the first place, that could potentially enhance efficiency when solving the CIO program in an online-fashion. It would then be compelling to test the NMPC version as an integrated motion planner and tracking controller within other contact-rich domains, such as locomotion, where multiple contacting bodies are involved, and contacts are generally sustained for longer periods of time.
-  ANYdrive. Note: http://www.rsl.ethz.ch/robots-media/actuators/anydrive.html Cited by: §II-B.
Practical Methods for Optimal Control and Estimation Using Nonlinear Programming. 2 edition, Society for Industrial and Applied Mathematics. External Links: Cited by: §I.
-  ANYpulator: design and control of a safe robotic arm. In IEEE/RSJ International Conference on Intelligent Robots and Systems, IROS 2016, pp. 1119–1125. External Links: Cited by: §II-B.
-  (2018-10) Trajectory Optimization With Implicit Hard Contacts. IEEE Robotics and Automation Letters 3 (4), pp. 3316–3323. External Links: Cited by: §I.
-  (2018) Deployment of an autonomous mobile manipulator at MBZIRC. J. Field Robotics 35 (8), pp. 1342–1357. External Links: Cited by: Fig. 3, §II-B.
-  Cited by: §I.
-  (2014-07) FORCES Professional. Note: embotech GmbH (http://embotech.com/FORCES-Pro) Cited by: §I.
-  (2015) Transcription Methods for Trajectory Optimization. Tutorial, Cornell University, Feb. Cited by: §I.
Design and use paradigms for gazebo, an open-source multi-robot simulator. In IEEE/RSJ International Conference on Intelligent Robots and Systems, 2004, pp. 2149–2154. External Links: Cited by: §III-C.
-  Dynamic underactuated nonprehensile manipulation. In Proceedings of IEEE/RSJ International Conference on Intelligent Robots and Systems. IROS 1996, pp. 889–896. External Links: Cited by: §II-A.
-  Motion Planning of Robot Fingertips for Graspless Manipulation. In IEEE International Conference on Robotics and Automation. Proceedings. ICRA ’04. 2004, Vol. 3, pp. 2951–2956 Vol.3. External Links: Cited by: §II-A.
-  Contact-Invariant Optimization for Hand Manipulation. In Proceedings of the 2012 Eurographics/ACM SIGGRAPH Symposium on Computer Animation, pp. 137–144. External Links: Cited by: §I.
-  (2012) Discovery of Complex Behaviors Through Contact-invariant Optimization. ACM Trans. Graph. 31 (4), pp. 43:1–43:8. External Links: Cited by: §I.
-  (2013) Animating Human Lower Limbs Using Contact-invariant Optimization. ACM Trans. Graph. 32 (6), pp. 203:1–203:8. External Links: Cited by: §I.
-  (2017) Whole-Body Nonlinear Model Predictive Control Through Contacts for Quadrupeds. CoRR abs/1712.02889. Cited by: §I.
-  (2018) A Comparative Analysis of Contact Models in Trajectory Optimization for Manipulation. CoRR abs/1806.01425. External Links: Cited by: §I.
-  Hybrid direct collocation and control in the constraint-consistent subspace for dynamic legged robot locomotion. In Robotics: Science and Systems XIII, Massachusetts Institute of Technology, 2017, External Links: Cited by: §I.
-  (2014) GPOPS-II: A MATLAB Software for Solving Multiple-Phase Optimal Control Problems Using hp-Adaptive Gaussian Quadrature Collocation Methods and Sparse Nonlinear Programming. ACM Trans. Math. Softw. 41 (1), pp. 1:1–1:37. External Links: Cited by: §I.
-  (1996) Glocker, C.: Multibody Dynamics with Unilateral Contacts. Wiley, New York. External Links: Cited by: §III-B.
-  (2014) A Direct Method for Trajectory Optimization of Rigid Bodies through Contact. The International Journal of Robotics Research 33 (1), pp. 69–81. External Links: Cited by: §I, §III-B.
-  An Implicit Time-Stepping Scheme for Rigid Body Dynamics with Coulomb Friction. In Proceedings of the 2000 IEEE International Conference on Robotics and Automation, pp. 162–169. External Links: Cited by: §I, §I.
-  (2007) Solving Normal Cone Inclusion Problems in Contact Mechanics by Iterative Methods. Journal of System Design and Dynamics 1 (3), pp. 458–467. External Links: Cited by: §I, §III-B.
-  (2010) Stochastic Complementarity for Local Control of Discontinuous Dynamics. In Robotics: Science and Systems, Cited by: §I.
-  (2007) Feedback Control of Dynamic Bipedal Robot Locomotion. Cited by: §I.
-  (2006) Trajectory optimization of mechanical hybrid systems using SUMT. Vol. 2006, pp. 665 – 671. External Links: Cited by: §I.
-  (2017) FORCES NLP: An Efficient Implementation of Interior-Point Methods for Multistage Nonlinear Nonconvex Programs. International Journal of Control 0 (0), pp. 1–17. External Links: Cited by: §I.