{'fraction_non_alphanumeric': 0.03363120976265975, 'fraction_numerical': 0.008263668684539252, 'mean_word_length': 4.5654188948306595, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Structured contact force optimization for kino-dynamic motion generation

Optimal control approaches in combination with trajectory optimization have recently proven to be a promising control strategy for legged robots. Computationally efficient and robust algorithms were derived using simplified models of the contact interaction between robot and environment such as the linear inverted pendulum model (LIPM). However, as humanoid robots enter more complex environments, less restrictive models become increasingly important. As we leave the regime of linear models, we need to build dedicated solvers that can compute interaction forces together with consistent kinematic plans for the whole-body. In this paper, we address the problem of planning robot motion and interaction forces for legged robots given predefined contact surfaces. The motion generation process is decomposed into two alternating parts computing force and motion plans in coherence. We focus on the properties of the momentum computation leading to sparse optimal control formulations to be exploited by a dedicated solver. In our experiments, we demonstrate that our motion generation algorithm computes consistent contact forces and joint trajectories for our humanoid robot. We also demonstrate the favorable time complexity due to our formulation and composition of the momentum equations.


I. INTRODUCTION
Trajectory optimization and optimal control approaches have recently been very successful to control the locomotion of legged robots. Simplified linear models of the dynamics, usually variations of the linear inverted pendulum model [1], [2], [3], [4] have been widely used to compute the center of mass motion of legged robots, in particular biped robots. The linearity can be exploited to formulate optimal control problems that can be solved with quadratic programs, therefore allowing for fast computations of solutions in a model predictive control manner.
While these simplified models are useful, generalization to more complex situations where multiple non co-planar contact points might be desirable is limited. Moreover, they do not take into account how the generation of angular momentum affects the behavior of the robot, nor how individual contact points should be controlled to create a desired motion.
On the other hand, several contributions have shown that using the full dynamics model of the robot could be beneficial to generate even more complex behaviors [5], [6]. They have the advantage of explicitly taking into account the contact interactions with the environment and the dynamics of the robot. The problem with such approaches is that they usually require to solve non convex optimization problems which are computationally demanding. Several approaches have addressed the problem using variations of differential dynamic programming to compute optimal trajectories. The drawback of such algorithms is that they do not allow to easily add constraints on the state and controls. Other approaches formulate the problem as non-convex optimization problems, however in these cases the solvers used are generally off the shelf nonlinear solvers that do not exploit the structure of problem. This paper addresses the problem of planning robot motion and interaction forces of a legged robot given a set of predefined contact points. In [7], the authors showed that such an approach could be taken to compute with a mild computational complexity trajectories that could be then directly included in whole-body control approaches and executed on a simulated humanoid. With a similar approach [8], interesting behaviors were demonstrated on a robot, confirming the interest in resolving such mathematical problems more efficiently.
In this paper, we analyze the structure of the problem induced by the physical model inherent of legged robots modeled through rigid body dynamics. In particular, we decompose the planning problem into two alternating optimization phases: first the contact forces and overall robot momentum are computed to satisfy dynamic constraints and in a second step the kinematic plan is resolved to satisfy the dynamic plan. Alternating these two phases lead to a locally optimal solution for the full kino-dynamic plan.
We concentrate our efforts on the momentum optimal control problem. First, we show that by choosing the right representation of contact forces and momentum we can rewrite the problem as a Quadratically Constrained Quadratic Program for which a convex approximation of the constraints can be explicitly (and trivially) constructed. We then present two formulations, using either a simultaneous or sequential optimal control formulation. We show how the sequential formulation of the problem can reduce the size of the problem while preserving its sparse structure for more efficient computation. Finally, preliminary numerical experiments demonstrate that exploiting this structure in a dedicated solver can significantly improve computational efficiency and quality of the solutions when compared to the same problem solved with a state of the art general purpose nonlinear solver.
The remainder of this paper is structured as follows. In Sec. II we propose a motion generation algorithm that is based on optimization of the momentum equations. We write out different representations of contact forces in Sec. III and discuss their algebraic properties. These properties are then exploited in Sec. IV to construct different variations of optimal control problems on the momentum equations of a floating-base robot. Then, in Sec. V we show experiments with the proposed motion generation algorithm and conclude the paper.

II. KINO-DYNAMIC MOTION GENERATION
We are interested in the problem of motion generation for a robot in the presence of contact forces. In this section, we are going to write out an algorithm based on an optimal control formulation that is suited for motion generation in floating-base robots. Our underlying dynamics model is where q is the joint configuration, M, N the inertia matrix and nonlinear terms, τ q are the actuation torques, S T the selection matrix, J e the contact jacobian and λ = [. . . f T e τ T e . . . ] T is the vector of forces f e and torques τ e acting at contact e. The equations of motion of a floatingbase robot consist of two parts, the manipulator dynamics describing the torque at each joint, and the 6 rows of Newton-Euler equations that describe the change of the overall momentum in the system. The dynamics can be separated into the 6 rows of Newton-Euler equations and the actuated part as follows with the contact location x e (q) and center of mass (CoM) x CoM (q) computed through forward kinematics. H(q) is the centroidal momentum matrix that mapsq onto the overall linear and angular momentum of the system [9]. M,N,J T e represent quantities from Eq. (1) corresponding to the actuated joints. From Eqs. (2-3), we observe that the planning process can be decomposed in two steps. First, a vector of joint trajectories q(t) is found together with contact force profiles λ(t) that satisfy Equation (2). Then, suitable torques can be computed readily from Equation (3) assuming that sufficient joint torques can be provided by the actuators. In the following, we will discuss an optimal control approach that finds joint trajectories q(t) together with force profiles λ(t) that satisfy Eqs. (2)(3). Our motion generation algorithm will be phrased as an optimal control problem of the form min.
where we express objective functions on kinematic quantities J t (q) as well as on momentum and contact locations J t (h t , λ t , c). For instance could J t (q) enforce an endeffector motion from one contact location to another. Although contact locations c e,t and momentum h t can be written as functions of q explicitly, we introduce redundant variables in order to simplify the optimization process as will become clear shortly. Eqs. (5-6) state explicitly how contact locations and momentum relate through forward kinematics. The momentum equations of our robot model (2) are then expressed as constraint in Eq. (7) with discretization parameter ∆. S is the set of feasible contact forces as will be discussed in more details later. The use of redundant variables allows us to decompose the overall problem into two better structured mathematical programs that allow the application of better informed solvers. The two sub-problems will be solved iteratively until convergence resulting in a solution for the original optimization problem. The first subproblem is defined as This problem is optimized over the momentum and contact forces only. The resulting force and momentum profiles are dynamically consistent and achieve reference profiles h,λ,c as close as possible. However, the joint state q is ignored in Problem (9). Objectives on kinematic quantities are considered in the kinematic sub-problem min.
where additionally to optimizing J t (q t ), we put constraints on momentum and contact location (cf. Eqs. (5-6)) from the original problem into the cost of our kinematic sub-problem in form of soft constraints. We will solve the original Problem (4) by iteratively solving problems (9), (10) and use the reference trajectories to enforce consistency between the two independent sub-problems.
Each of the sub-problems minimizes parts of the original objective function and satisfies a subset of the original constraints. Additionally, the resulting momentum profiles from one sub-problem are optimized for consistency inside of the solve Problem (9) h,c e := h, c e until Solution of subproblems (10), (9) do not change other up until convergence. With the proposed decomposition we have to solve two easier problems, where Problem (10) is an unconstrained optimization that has no notion of contact forces. Similar formulations, usually without costs on momentum, are often addressed in trajectory optimization for manipulators [10]. Problem (9) is a non-convex, but a well structured optimal control problem as we will discuss in the remainder of this paper.

III. REACTION FORCE REPRESENTATIONS
It is well known that the (linear and angular) momentum of a system can only be changed by external forces. In this section, we will discuss forces acting at contact points on a robot structure and we will show how the choice of representation of contact forces results in different mathematical programs with beneficial properties to be exploited by optimization algorithms.

A. Contact Forces
A force f e and torque τ e acting at contact point p e can be represented equivalently as a force and torque acting at the CoM r where the torque acting at r transforms to and f e remains the same. The change of momentum can then be expressed equivalently using either of the two representations aṡ where , k are the linear and angular momentum. From Eq. (12) we can see that depending on how we chose the point of action the dynamics are either linear or not. One might be tempted to say that the linear dynamics lead to easier mathematical programs, which however is not the case when contact constraints have to be considered. Unilateral contacts, e.g. a hand pushing against a wall require that f e remains in a friction cone. On the other hand, point contacts cannot generate torques, i.e. τ e = 0. In both cases, the center of pressure p e at which f e is acting (effectively) needs to reside inside of a support polygon or at the point of contact c. We will now derive constraints on contact forces that generalize across various types of contacts, such as contact between surfaces (for instance a flat foot on the ground) or point contacts (for instance an elbow pushing on the table). For the following derivation, we define a contact surface for each contact e as shown in Fig. 1. Each contact is located on a surface described by a rotation matrix R e = R e,x,y R e,z and translation t e . We can then represent forces, torques and locations in local coordinates as follows where we drop the index e for better readability where the hat identifies local representations. We can now write contact constraints as where µ is a friction parameter of the surface and p max is a vector of support limits and τ max is a bound on torques that are normal to the contact surface. If we want to describe a point foot, we set τ max = 0, p max = 0. A contact between surfaces, e.g. floor and foot sole, allow for a wider range ofτ ,p. In Eq. (19) we restrict the force to remain in a pyramid which is used as an approximation for a friction cone. Note that in Eq. (18) we consider the support polygon to be rectangular for a clear notation, however, our derivation holds for more general polyhedra as well.
As expressed in Eqs. (17-19) contact constraints are linear if we express them as functions of f e , τ e acting at p e . On the other, if we chose to represent forces and torques as f e , κ e acting at the CoM r, we will end up with non-linear contact constraints. In both cases, either the constraints (18) become nonlinear or the momentum equations (12) do. In the following, we will investigate the type of non-linearity that we introduce in Eqs. (12), (18) and how we can exploit it for development of optimization algorithms.

B. Algebraic categorization
Depending on the choice of representation of the momentum dynamics (force acting atp or at r) the resulting optimization problem will be constructed from different sets of linear or non-linear functions. We are interested in tracking the algebraic properties of these functions as we construct mathematical programs using the dynamics and constraints from Eqs. (12), (18). Resulting objective functions and constraints will be classified as elements from the following sets of functions.
Definition 1: We define the set of affine functions We define the set of positive semi definite quadratic functions We define the set of differences of positive semi definite functions Throughout this paper we consider functions to be elements of Q ± , only if the quadratic parameter matrices Q i , P i can be accessed separately (instead of only having access to the indefinite matrix Q i − P i ).
It is worth noting that Q ± is closed under addition, scalar multiplication and composition with functions from A, i.e.
Specifically, if we have quadratic parameter matrices Q i , P i constructed separately then the operations in Eq. (20) preserve this separation. As can be noticed from Eqs. (12), (18) the nonlinearity introduced in the equations stems from cross products which can be classified as follows Theorem 1: The function ×(a, b) : R 6 → R 3 , ×(a, b) = a × b is an element of Q ± and we can write out the corresponding matrices Q i , P i explicitly.
Proof: It is straight forward to write a general cross product in a quadratic form with [a × b] i = a T b T H i a T b T T , i = 1, 2, 3. Using an eigenvalue decomposition (for symmetric real matrices), we can then Note that constructing the quadratic parameter matrices Q i , P i for a general cross product can be carried out offline as it does not require any problem data. Through out the construction in the remainder of this paper, we will not have to perform any eigenvalue decomposition as this is typically costly. Instead, we exploit Eq. (20) in order to construct functions from Q ± . It is straight forward to show that Eqs. (17-19) can be expressed asŝ(f ,p,τ ) ≤ 0, withŝ ∈ A, i.e. if we chose to represent contact forces atp then the contact constraints are affine. On the other hand, the angular momentum (12) will result in a non-linear function. In fact, we will show thatk(f ,p,τ ) ∈ Q ± and we can construct the parameters Q i , P i separately inexpensively. On the other hand, if we chose to represent the force torque at the CoM, we have the reverse effect, i.e. the angular momentum becomes an affine function of the force-torques and as we will show the contact constraints will be expressed as s(f , κ) ≤ 0, s ∈ Q ± .

C. Algebraic categorization of contact forces
In this section, we will categorize the nonlinear functions in Eqs. (12), (17-19) into one of the categories defined in Sec. III-B. Categorizing the CoM torque κ(f ,τ ,p) (cf. Eq. (11)) requires consequent usage of Eq. (20) and Theorem (1). As a result one can show that κ(f ,τ ,p) ∈ Q ± and thusk(f ,τ ,p) ∈ Q ± . Further, we will show that choosing to express external forces at the CoM will result in contact constraints that are expressed as inequalities on functions from Q ± . Assuming a unilateral contact, we can write the torque acting at the center of pressure aŝ From the definition of the CoP we know thatτ vanishes when projected on the contact surface, i.e.
In order to resolve forp, we first invert the premultiplied projector We note that from the friction constraint (cf. Eq. (19)) we have R T z f > 0. Next, we resolve Eq. (22) forp and substitute it into the CoP constraints (cf. Eq. (18)). For better readability we will only show the result for the right-hand side of Eq. (18) where we used Eq. (24). Applying Eq. (20) and Theorem (1), one can transform Eq. (27) into a function inequality s(f , κ) ≤ 0, s ∈ Q ± and again the quadratic parameter matrices can be constructed inexpensively.

IV. OPTIMAL CONTROL FORMULATIONS
Through out the following discussion, we will restrict ourselves to a predefined contact activation pattern (cf. α e in Eq. (12)). That means we predefine which endeffector is in contact with the environment at which time. This could for instance be computed with an acyclic contact planner [11]. The location of the contact and the force-torque are expressed as decision variables and will be computed by the optimizer. The goal of this section is to define an optimal control problem that finds sequences of states and controls under the dynamics constraint in Eq. (12) and contact constraints Eq. (17-19) thus solving problem (9).

A. Simultaneous Formulation
We start out with a simultaneous formulation, i.e. we optimize over both states and controls and express the dynamics constraint explicitly in the mathematical program. Here we represent contact forces at the center of pressure. We summarize the optimization variables into x 0:T = [. . .f T t,1:Ep T t,1:Eτ T t,1:E h T t+1 . . . ] T , t = 0 . . . T and c 1:E and write our optimal control problem as min.
This formulation is a Quadratically constrained Quadratic Program (QCQP), with a positive definite objective, linear inequality constraints and nonlinear equality constraints on functions from Q ± . As one can note from (28-31), summands of the objective as well as constraints all depend only on time-local variables x t:t+1 and on c 1:E . This property has the advantage that optimization algorithms can be applied that exploit the sparsity inside of objective Hessian and constraint jacobian. The constraint jacobian has the form  where we have a dense part on the right-hand side, which however does not depend on the chosen time discretization T , but only on the number of contacts E. The top left-hand side is a band-diagonal structure. The Lagrangian L [12] of the optimization program in Eq. (28) has a Hessian that can be decomposed into 4 blocks as follows We note that again the dense blocks on the left and bottom do not depend on the time discretization T . The top-left block on the other hand is block-diagonal. This sparsity is favorable in optimization algorithms, because typically an essential step in these solvers is to solve a linear system that inherits the sparsity properties from the constraint Jacobian and Lagrange Hessian. Especially, interior-point methods are known to preserve the band-diagonal sparsity structure [13].

B. Sparse Contact Force Representation
In the following we will show that the sequential optimal control formulations for our problem can be transformed to an equivalent problem with a sparse structure similar to the simultaneous formulation. We will assume that a contactphase e is defined as the time interval t ∈ [σ e , ε e ) at which an endeffector is touching the environment. Thus the contact activation from Eq. (12) can be refined to where contact-phase e starts at t = σ e and ends at t = ε e . For instance, a contact-phase could describe a foot step e where the heel hits the ground at t = σ e and the last time step at which the toe is touching the ground is t = ε e − 1. Two heel strikes with the same foot would correspond to two different contact phases e, e . Further, we define forces and CoM torques as linear functions of a new set of variables ϕ t , ψ t as f e,t = ϕ e,t − 2ϕ e,t−1 + ϕ e,t−2 (35) κ e,t = ψ e,t − 2ψ e,t−1 + ψ e,t−2 (36) ∀t < σ e : ϕ t = ψ t = 0 The variables ϕ t , ψ t correspond to the twice integrated contact forces and torques. General optimal control problems that are expressed in sequential form, usually lose their sparsity structure because the state variables need to be integrated out. In the following we will show that this is not the case for our momentum optimization problem. We will derive expressions for the state variables r, h,ḣ as functions of ϕ t , ψ t . The goal here is to show that the state variables at time step t do not depend on all previous time steps t = 0 . . . t, but only on a few ones. This has the advantage that a sequential formulation expressed with ϕ t , ψ t will maintain its sparse structure. We compose the linear momentum rate aṡ t (ϕ e,t−2:t ) = M g + ˙ e,t (ϕ e,t−2:t ) In order to construct expressions for the linear momentum and center of mass, we first define the sums of forces during a contact-phase over time , else With the derived expressions we can now construct the linear state terms as followṡ t = M g + e˙ e,t (ϕ e,t−2:t ) e e,t (ϕ e,t−2:t , ϕ e,εe−2:εe−1 ) M r e,t (ϕ e,t−2 , ϕ e,εe−2:εe−1 ) It is interesting to see, that the state variables at time t only depend on control variables at t ∈ {t − 2, t − 1, ε e − 1, ε e − 2|e = 1 . . . E}. This has the advantage that derivatives of the objective function and constraints remain sparse and as such result in more efficient optimization algorithms similar to the simultaneous formulation. In the previous derivation we obtained expressions for the linear momentum, but expressions for the angular momentum can be readily derived with the same sparsity properties.
x 0:T ,y 1:E T t J t (x t−2:t , y 1:E ) ∈ Q + (38) We can see that similar to the simultaneous formulation, summands of the objective and constraint functions only depend on x t−2:t , y 1:E rather than on all time steps. We will now write out the sparse constraint jacobian We will now write out the Hessian of the lagrangian

D. Comparison of sparsity patterns
As we analyzed in the previous sections, both the simultaneous as well as sequential formulation of the mathematical program (9) have sparse constraint Jacobians and Lagrangian Hessians. A common approach to solve optimal control formulations of that kind is to use interior point methods [14] as computation time benefits from these sparsity properties. Usually, the computation time is dominated by a decomposition of a augmented system matrix ∇ 2L +G T ΣG, where G is the constraint jacobian and ∇ 2L is a (preferably convex) approximation to the Lagrangian Hessian. Comparing the sparsity patterns of the two formulations from Sec. IV, we see that both problems result in augmented systems that have a block structure with a band-diagonal matrix on the top left. The computational complexity for decomposing the augmented system is in both formulations linear w.r.t. the discretization parameter T . Finding an approximation ∇ 2L that both makes progress towards a solution and can be robustly decomposed is usually not trivial. In our problem however, a convex approximation can be constructed exploiting the fact that ∇ 2 L can be decomposed into a difference of p.s.d matrices (cf. Eq. (47)). We can obtain ∇ 2L by simply dropping the negative part of ∇ 2 L from Eq. (47). Note that we do not require expensive operations such as eigenvalue decompositions to be carried out online, but obtain the convex part of ∇ 2 L through consequently applying Eq. (20) and Theorem (1) as we construct the problem. This can be done for both, the simultaneous as well as the sequential formulation. The number of variables is very different for the sequential and simultaneous formulations and has an impact on the computational complexity of resulting optimization algorithms. In the simultaneous case, we have to consider 6 force variables for each endeffector and 9 state variables. For instance a biped that remains in double support at all time requires 21T variables. On the other hand, the sequential method does not require variables for states and thus contains only 12T variables for the aforementioned example. As a consequence computational complexity is reduced with a sequential formulation as we solve for longer and finer grained time horizons T . It is also worth mentioning that since the momentum dynamics generalize the commonly used LIPM, the sparse structure in our sequential formulation will be similar using LIPM dynamics. This nicely complements the discussion in [13].

V. EXPERIMENTS
We demonstrate the proposed motion optimization Algorithm II on a stepping task for a humanoid robot. Further, we convince ourselves that the sparse structure in our momentum problem results in linear computational complexity w.r.t. time discretization T . All experiments were performed on a notebook with a 2.7 GHz intel i7 processor with 16gb RAM.

A. Optimization Algorithm
In this section we would like to discuss our choice of optimization algorithms for solving the two subproblems in (9), (10). In both algorithms, we exploit the sparse structure resulting from the optimal control formulation. Further, we condense the number of variables in the momentum optimal control problem (9) as we discussed in Sec. IV-C. We solve the optimization problem (10) with a Gauss-Newton method, where the Hessian is a block tridiagonal positive definite matrix that can be factorized with a dedicated cholesky decomposition [15]. The bottleneck in our implementation is a numerical differentiation of the momentum jacobian ∂h(q,q) ∂q , which should be constructed directly from the dynamics model. Nevertheless, with our naive implementation we were   kinematic pass optimizes for an unfeasible CoM lateral motion, which is corrected after the first momentum optimization. In the second iteration of the algorithm, the kinematic plan is adapted to realize the CoM sway. Horizontal angular momentum (bottom two plots) is also found to be consistent across kinematic as well as momentum optimization.
able to generate whole-body motion and force plans in 30s. For the momentum optimization (9) we implement a primaldual interior point method [12,Algo. 19.2], which allows us to preserve the band-diagonal structure in our problem. The algorithm is modified as discussed in Sec IV-D to exploit the convex approximation of the Lagrangian Hessian. We run the momentum optimization for increasing values of T and plot convergence of the Lagrangian gradient in Fig. 4. It becomes evident from the plot that computation time increases linearly as we increase T . Note that sequential methods in general scale cubicly with T , however, our sparse formulation allows to preserve linear complexity potentially leading to more efficient algorithms. In a side experiment we investigated the effect of our approximation ∇ 2L on convergence to a solution. We implemented a very basic sequential quadratic program [12, Algo. 18.1] that was using dense operations, meaning it was not exploiting the sparsity structure. Instead we focused on the benefits of the approximation ∇ 2L . Although much slower than the sparse interior-point method, our algorithm still converged by an order of magnitude faster than an offthe-shelf SQP method, SNOPT, [16] motivating the need for dedicated solvers in robotics.

B. Motion Generation
We construct a stepping scenario 1 for a simulation of our Sarcos humanoid robot (cf. Fig. 2). The humanoid is to step on stepping stones that increase in height and are angled as illustrated in the figure. Note that here the LIPM assumption does not hold, but instead non-coplanar contacts have to be considered together with angular momentum that is required for leg swing motions. In our experiments, we fix the footstep locations c e to predefined values for a simpler implementation and compute coherent joint trajectories and force profiles using Algorithm II. Every 0.7s contact is broken or created switching from single to double support or vice-versa. Our initial guess for the CoM and momentumh is to simply have zero momentum and keep the CoM centered above the feet. We introduce weights into the cost function in Eq. (10) giving more importance to footstep locations. In our example, Algorithm II converges already after 30s for T = 100, despite the inefficient implementation and requires only 2 passes to acquire consistent force and joint trajectories. The first kinematic plan realizes our naive CoM motion moderately well, which however cannot be accomplished given the contact constraints. Thus, the momentum optimization introduce a CoM sway in lateral direction, which is then realized by the kinematic optimizer in the consecutive iteration as can be seen in Fig. 3. Additionally, the kinematic plan introduces an arm sway motion after the first pass to compensate for angular momentum introduced by the swing legs. After the motion generator converges it yields a kinematic plan that is coherent with bounded force and CoP profiles (cf. Fig. 5). In this example, we keep the ankle joints fixed and do not account for object collision. However, collision avoidance can be added inside the kinematic optimization using trajectory optimization techniques [10]. Our experiments demonstrate that constrained force profiles and consistent whole-body motion plans can be generated efficiently with the proposed algorithm. The structure identified in the momentum optimization led to more efficient optimization algorithms.

CONCLUSION
We presented a contact-centric motion generation algorithm for floating-base robots. It consists of two subproblems optimizing for kinematic trajectories and momentum profiles iteratively. An analysis of the momentum optimization is presented and the relevance for optimal control formulations is discussed. In our experiments, we demonstrate the momentum-centric motion generator on a stepping task for a humanoid robot simulation. Indeed, consistent momentum and whole-body motion plans can be acquired efficiently. Further, our experiments show favorable scaling properties of our sparse optimal control formulation. 1 The result is visualized in a video: http://www-amd.is. tuebingen.mpg.de/˜herzog/structured_momentum_opt. mp4 Computation time of the momentum optimization (9). We execute the optimization for increasing values T and plot the resulting norm of the Lagrangian gradient (top) and the overall computation time (bottom). Although our implementation still requires improvements, we can already see a linear trend in the computation complexity confirming our analysis in Sec. IV-C